Homework 7
1. Implement Radon Transform & then Filtered Back
Projection to reconstruct back the imag
a) Take any image, f, first implement Radon Transform in the discrete
domain. To do this accurately,
First compute 2D Fourier
transform of f, estimate approximate bandwidth in x and y
directions, choose zeta_0 = maximum of x and y bandwidths
d = 1/(2 zeta_0)
Choose Delta = pi/180
Compute Discrete Radon
Transform g(md, n Delta), m goes from -M/2 to M/2-1,
M = 2 D d where D is diameter of the object (discretize the
integral given in equation 10.6)
b) Reconstruct back the image by implementing Filtered Back Projection
(as I explained in class, also read the book) with g(md, n Delta) as
input.
Be very careful: you will have to
implement an approximation of the high pass filter H(f) = |f|, some
examples are given in the book
c) Compute the reconstruction errror. Increase zeta_0, and see the
trend. Vary Delta and see the trend
2. AK Jain Chapter 10, Problem 10.6
3. Extra credit: AK
Jain Chapter 10, Problem 10.16
Homework 6 (Gaussian Mixture Model:
Alternating Maximization and EM Algorithm)
You are given N data points Y1, Y2,... YN which are samples from
a Gaussian mixture with K components. All Yi's are assumed independent.
Each Yi is a P-dimensional vector. For e.g. each Yi can be
the R,G,B color intensity of an image pixel, in that case, P=3 and N =
number of pixels in the image and K is the number of different regions
in the image (if 2 objects and background, then K=3). Also, for
simplicity, assume Sigma_c for c= 1,2,.. K are diagonal covariance
matrices.
The above problem can be reformulated as follows (in order to use Alt
Max or EM): There is hidden data Xi which denotes the class membership
of Yi. Conditioned on Xi, Yi is Gaussian distributed with mean mu_{Xi}
and covariance Sigma_{Xi}. Also the prior probability of Xi taking the
value c is p_c, i.e.
p(Yi|Xi) = N ( Yi; mu_c,
Sigma_c) and p(Xi = c) = p_c for c = 1,2,... K
Denote X = [X1,
X2,... XN] and denote Y
= [Y1, Y2,... YN] and denote the parameters as theta = { mu_c,
Sigma_c, p_c, c = 1,2,...K}
A. Use alternating maximization to compute max_{theta, X} p_{theta}( Y, X ), i.e. derive the
expressions for both steps of the alternating maximization.
I gave the main idea for both steps in class. When theta is held fixed,
you just get a Maximum Likelihood classifier. When X (class membership)
is held fixed, you need to take first derivative and set it to
zero to get the max.
B. Use EM algorithm to compute max_{theta} sum_{ X } p_{theta}( Y, X ), i.e. derive the
expressions for the E step and the M step.
The E step was done in class, for M step also I gave the main idea in
class.
C. Extra: Pick a colored image of any object (K=2),
and apply both algorithms. Manually changed the object intensity so it
is not too different from background intensity. Then apply both
methods. You will see that EM will work better in that case (when
clusters close to each other).
Note both
algorithms only give local maximas
Homework 5 (Extract contours, estimate
geometric distortion & correct for it)
First use the circle image
(convert to greyscale).
Apply a geometric distortion to the image. Try rotation by 40 degrees,
try uniform scaling, and try any type of affine deformation (x' = ax +
by+c, y'=dx + ey + f). For applying any of these, first shift the
origin to the center of the image, e.g. a 21x21 image will have
coordinates x,y belonging to [-10, 10] x [-10, 10]. Also, remember you
need to interpolate grey scale values for the distorted image
(interp2 command can be used or write your own code).
Note rotation will have almost no effect on circle (except due to
interpolation errors).
Extract contour of the object from the original image and extract
contour of the distorted image. If you're using the circle image, just
threshold it to get a binary image, followed by using either
bwtraceboundary or use the contour command in Matlab (try both). Both
perform boundary tracing, but contour gives the boundary in sub-pixel
coordinates, while bwtraceboundary gives the innermost boundary pixel
locations on the object.
In general, the two contours will have different number of points and
will also be of different lengths (remember arclength is different from
number of points for the contour, I will explain this in next class).
Use interp1 to resample both contours to M points where M = round( (L1+
L2)/2) where L1 and L2 are the lengths of the two contours. Once
you have M points on either contour, estimate the geometric distortion
as taught in class, and then correct for it as taught in class. Again
remember to do grey scale interpolation (interp2 command can be used or
write your own code).
Now, try the same thing
for any real image. Here you may need to use edge detection followed by
edge linking to get the contours or other segmentation methods. Or use
feature matching (KLT tracker) to get corresponding points.
Extra: Understand/experiment
with using the Matlab command imtransform.
Homework 4 (Image Enhancement & Wiener
filtering) (you can use all
MATLAB
functions in this homework)
Do at least 4 of the following 5 parts. You're welcome to do all!
By the way, there is a command called dicomread in Matlab - to read
.dcm images
A. Demonstrate
using MATLAB that the covariance matrix of a zero mean periodic
stationary random signal is circulant.
B. Histogram Specification & Equalization.
Take any image. First do a histogram specification to get the pixels to
follow a discrete version of the truncated Gaussian PDF N(20,3)
truncated between [10, 40]. For e.g. the probability that a pixel takes
value 12 is the truncated Gaussian integrated between 11.5, 12.5.
Matlab has standard functions to compute Gaussian CDFs which can be
used to evaluate this.
Then do a histogram equalization on it, so that the final histogram
approximates a uniform between 0,255.
Also try the same histogram equalization method on the following
image: brain MRI
C. Do edge enhancement (using either high pass filtering or unsharp
masking or bandpass filtering) on the following image: brain MRI
Try the same operations in DFT domain.
D. Compare Low Pass Filtering (say averaging or Gaussian filter) with
Median filtering & Directional smoothing for the following
images:
(i) brain MRI , (ii) same image with salt
and pepper noise added, (iii) same image with Gaussian noise added.
See the benefit of Directional Smoothing in this image which is already
pretty blurred, so if you do smoothing in all directions, edges will
get even more blurred.
E. Perform FIR Wiener smoothing (assume no blurring) for
all the three images used in part D. Compare with algorithms of part D.
I
think you should write your own code for this one. Start with designing
a 5x5 Wiener smoother for every 16x16 block in the image.
Increase/decrease Wiener smoother size or block size as needed.
Try changing window size to keep output
SNR constant. Compare performance
using MSE or visual inspection.
Homework 3 (postponed, due Feb 16)
Find an application of DFT in your research area (or in any
image processing application) and implement it. Of course it doesn't
necessarily have to be a "new" one. Some examples are:
- use of DFT to detect regions and times where there is fast
motion:
apply DFT to image differences and look for high frequency regions
- use DFT (or DCT) to replace eigenvectors of coviance matrix (e.g.
PCA) in any application. Can you define something like LDA but using
DFT or DCT
- study Fourier descriptors literature (for feature extraction,
recognition, retrieval...). FD is the Fourier series of some feature of
the boundary contour of an object as a function of contour arclength,
e.g. tangent angle, x,y location (complex FD), ... F series (defined
only for periodice signals) is the same as DFT of one period
- study the use of face asymmetry as a biometric for face
recognition
(a recent PAMI paper uses the energy in the imaginary part of the DFT
coefficients of each row as a measure of asymmetry).
- directional Fourier transforms... connection with Radon transform
Extra (this may be part of HW 4): Demonstrate
using MATLAB that the covariance matrix of a zero mean periodic
stationary random signal is circulant?
Extra: Find an application of
circulant matrices outside of signal/image processing and implement it.
One domain is graph theory.
Homework 2 (Due Friday Feb 2)
- Problem 4.9 of AK Jain.
- Problem 4.11, part (a). Also verify using MATLAB. A continuous
function can be simulated by choosing x = [0:1/(100 zeta_s):
5/zeta_0] (using a very small sampling interval).
- Demonstrate Moire effect for a DC source (constant grey scale
image) using MATLAB. Note what I said in class wasn't exactly correct.
You will need to decimate the image first and then interpolate using
different size averaging filters.
- For the chessboard image of homework 1 (or any one row of it),
write out equations for every image you got after decimation and after
interpolation, i.e. do the entire thing on paper. Some of you haven't
done it in the submission. One row of the chessboard can be written out
as a difference of two unit step functions: u(x) - 2u(x-20)
Homework 1 (Decimation &
Interpolation of Images) Try to finish by Jan 19 or 22.
Use as I(m,n) a chessboard image
or a fence image.
A. Decimation. First
decimate without applying the low pass filter (Decimate in two steps:
first do zeroing in the original image and then downsample). See the
effect of aliasing in the decimated image. Start with decimating by a
factor of 2 and keep going until you "see" the effect of aliasing. Let
the chosen decimation factor is D. Now apply a low pass filter (say a
3x3 or a 5x5 averaging filter) to I(m,n) and then decimate by D. Keep
increasing size of the averaging filter until the effect of aliasing is
"negligible".
B. DFT. Take a 16x16 piece of
the image and take its DFT (or you may even just take a row of the
image and do 1dim DFT). Do this for all the images you get in the above
steps. See the aliasing in the frequency domain and see the aliasing
getting reduced if you low pass filter first.
Note DFT
should help you decide what factor to decimate by to "see" aliasing,
and then also decide how big an averaging filter to apply in order that
the effect of aliasing is "negligible"
Details for decimation by 2
(similarly do for other factors D)
- Take the image I(m,n), Set every alternate row and column to
zero. Call this I_Z(m,n). Now downsample to get I_D(m,n) which is half
the size.
- Compute DFT of I, I_Z and I_D . Either do 1dim DFT or do 2dim DFT
of a 16x16 or 32x32 region where you can "see" aliasing most easily
- Low pass filter the original image and repeat the above steps
- Do the same thing for different decimation rates D and for
different
C. Interpolation. Take the
decimated image I_D(m,n), add D-1 zeros to along rows and columns to
get an original size image I_U(m,n). Interpolate the upsampled image
(try zero-order hold, linear or cubic interpolation) to get the
reconstructed image I_hat(m,n). Do this both for the decimation of the
original image and decimation of the LPF'ed image. Notice that the
effect of aliasing will not go away by interpolation. To get a
quantitative measure, you may want to look at the mean squared error
between the reconstructed and original image in both cases. Note for the LPF'ed case, the original
image will be the LPF'ed one.
D. DFT. Compute DFT
(again you may choose to just do a 1dim DFT of a relevant row since its
easier to plot and see magnitudes) of all images (I_U and I_hat in both
cases) and connect with what you learn in class.
Submission: Send me via email
(or put in your webspace and send me the link, need not create an
explicit webpage)
- Send MATLAB code.
- Send a short
report "explaining" what you did and what you observed and why based on
what you've learnt about decimation and interpolation in class.
- Can be a handwritten report if you like. You won't lose any points
for that. But I'll suggest learning LaTex!