Lab (3) - Image Processing
Frequency Domain Processing (Fourier Transform)

Objectives:
To learn about the Discrete Fourier Transform (DFT)
To learn about the Fast Fourier Transform (FFT)

Please DO NOT Print this lab as it will not appear correctly

Preparations:
As usual, create a MS Word document and submit your answers and the requested figures as one file for the in-lab activities and another one for the post-lab.

Fourier Transform
In this lab, we will see some of the applications of DFT and FFT.  In class we talked about the 1-D and 2-D discrete Fourier Transforms. The 1-D Discrete Fourier Transform (DFT) is defined F(u) as:

The 1-D Inverse Discrete Fourier Transform (IDFT) is defined as:

In the first equation i = 0, 1,..., M-1 represent the number of samples of the continuous function and in the second equation those corresponds to their corresponding Discrete Fourier Transform coefficients.

Any x, (xi) can be computed using:

where the second term (delta x) is the width of the grid on x direction. Any value of u, ui, can be computed using:

where the second term (delta u) is defined as:

The 2-D DFT is defined as:

The 2-D IDFT is defined as:

For the 2-D case, in these equation i = 0, 1,..., N-1 and j = 0, 1, ..., M-1 represent the number of samples of the continuous function in the x and y direction respectively and in the second equation these values correspond to the Discrete Fourier Transform of those points.

Similar to the 1-D case, any x or y can be computed using the initial value and the number of grids that that point is away from the origin. Thus,

,

from these two, we will find:

,

where,
MATLAB has some built-in functions for computing the 1-D and 2-D DFT.   If the grid size is a power of 2, the calculation is done in a faster scheme called Fast Fourier Transform.  By default, MATLAB attempts to solve the problem using FFT, however, in the case that the grid size is not a power of 2, then it uses DFT.

That being said, let's solve some examples.

Fourier Transform in MATLAB
We want to first conduct some experiments with Fourier Transform.  If you run the following three commands you will produce a function displayed on Figure (1).

f = zeros(512,512);
f(245:265,235:275)=1;
imshow(f)

Figure (1)

We can find the Fourier Transform of this function:
F = fft2(f);
fft2 produces real and imaginary values of a 2-dimensional function such as f.  The real values are as small as -5.3731e-007 and as large as 183.   As we mentioned in class one way to avoid the problem dealing with both the imaginary and real values is to compute the spectrum of the result.  The spectrum is defined as:

|F(u,v)| = (R(u,v)2 + I(u,v)2 )1/2

Where R represents the real part and I denotes the imaginary part. We can find the real and imaginary part of the result using the following two commands respectively:

r = real(F) and  i = imag(F).

As you may have noticed it is not possible to display negative values. Also displaying large values would not produce a good grayscale image.  Sometimes to display the result of a Fourier Transform we use a log transformation of the positive values to convert the large values to a more reasonable range.  Note that log will reduce the values between:
[10-100] to [1-2],
[100-1000] to [2-3],
and so forth.

>> F2 = log(abs(F));
>> imshow(F2)

This produces:

Figure (2)

You can see the periodicity behavior on this image.  The four bright corners come from periodicity property of Fourier Transform.  Usually, we wish to see the very bright spots (main phase) at the center.  One way to do this is to multiply every element of the mesh by (-1)x+y, i.e., by  (-1) to the power of  the sum of  indices on x and y axes.  This will shift the phase to the center and distribute other points around the center.  For example a pixel at row 2 and column  3 would be multiplied by (-1)2+3= -1.   Following is a codes segment that does this for us.

One way to do this is to use a for loop (you can create a temp.m file and copy and paste these lines there to run):

% temp.m
F = fft2(f);
F2 = log(abs(F));
imshow(F2)
for i = 1:512
for j = 1:512
f(i,j) = (-1)^(i+j)*f(i,j);
end
end
F = fft2(f);
F2 = log(abs(F));
figure,imshow(F2)

This produces:

Figure (3)

There is also a command in MATLAB that has the same effect:

Fc = fftshift(F)  % to shift the values after computing the Fourier transform

The result of the shift should be the same as that we obtained by multiplying pixels by (-1)x+y.

Lab Activity (1)
Assuming you have gone through above example and have obtained Figure (2) and  (3).  We now want to conduct an experiment.  By executing the following commands, you will move the rectangle on the original image by 40 pixels to the left.

>> f = zeros(512,512);
>> f(245:265,195:235)=1;
>> imshow(f)

Figure (4)

How do you think the centered Fourier transform of the new rectangle will look like?  You can do this by modifying the temp.m file by replacing the line that produces the rectangle with:
f(245:265,195:235)=1;
that moves the line to the left.

Run the script.
Question1: What difference moving to the left by 40 pixels had on your output?

In addition to moving to the left by 40 pixels, try 3 more cases by moving the rectangle
1) move the original line up by 40 pixels,
2) move the original line up and left by 40 pixels on each direction

Question2: Explain the affect of moving in each direction on the output Fourier transform. Include your answers and the resulting 4 images in your in-lab file.

Lab Activity (2)

In this activity, we will compute a 1-D signal and will compute the spectrum and phase for that signal.  We begin by creating 20 values between 0 to 720 degree:

t = [0: 10 : 720]*pi/180;     % 20 angles between 0 to 720 in Radian (2 full period)
s = sin(t);
plot(t,s);

Figure (5)

Fsin =  fft(s);
F2 = abs(Fsin);
figure, plot(t, F2);

Figure (6)

Fc = fftshift(F2);   % to center it
figure, plot(t, Fc);

Figure (7)
Just to highlight a few things, Let's repeat this for sin(2x).  Sin(2x) is the same wave function as sin(x) but the frequency is twice larger.
s = sin(2*t);
figure, plot(t,s);

Fsin =  fft(s);
F2 = abs(Fsin);
figure, plot(t, F2);

What kind of difference do you observe, if any?

I now will create a signal: sin(t) + sin(2t) + sin(4t) and will see what the above analysis will do on this one.

s = sin(t) + sin(2*t) + sin(4*t);
figure, plot(t,s);

Fsin =  fft(s);
F2 = abs(Fsin);
figure, plot(t, F2);

What do you learn from this spectrum?

Question3:
Consider Figure (5), what is the frequency (number of cycles per unit time (sec) ) of the signal?  What is the amplitude (the largest absolute value) of this signal?

Question4:
What kind of relationship do you observe between the
sin signal and its Fourier Transform?

Question5:
On Figure (6), What will the distance between the two maximums suggest?  What is the relation between these peaks and the behavior of the signal on Figure (5).

Lab Activity  (3) - If you didn't finish this part in the lab, do it for post-lab
One of the applications of Fourier Transform is to locate features in an image by performing convolution.  For instance occurrences of a letter in a text.  The MATLAB Image Processing Toolbox has an example in which the letter a has been extracted from a sample text.  The procedure for extracting that a letter is explained below.  Before you try this please copy three image files: text.jpg , e.jpg and n.jpg to your directory.   Please do not chnage the size of the images.  Save and use them as they are.

These two images are the images of letters "n" and "e" with the same size as they appeared in the text.jpg image file. The idea (if it works) is to find where the location of "e" or "n" is in the image shown below.

text.jpg Image

The correlation of the image e.jpg and text.jpg can be computed by:
1) Read the text.jpg image and store it in an array called text. Read the n.jpg and e.jpg images and store them in corresponding n and e arrays.

2) use imshow to display the text image.

3) use figure and imshow to display one of the letter images (let's say n).

4) Rotate the letter image by 180:
rn = rot90(n, 2); % this rotates n by 90 twice (180)

5) Compute the Fourier Transform of the rotated image and bring the size to 256x256 by padding. Note that the text image is 256x256.
rn = double(rn);  % convert to double
fn = fft2(rn, 256, 256);

6) Compute the Fourier Transform of the text image:
ftext = fft2(double(text));

7) Multiply element-by-element the resulting two transfroms.
invboth = ifft2( ftext .* fn);

8) Take the real part of the result.
realboth = real(invboth);

9) Find the MAX value in invboth. You need this value for Thresholding. Suppose you found that the MAX was 2700000. Then set the threshold value a little below this. Let's say 2500000.
maxt = max(max(realboth))

Suppose you got 2700000.  This is just an example. Set a threshold value based on the maximum value, for example:
thres = 2500000;

Note: I usually take the MAX and subtract 10% from it to get the threshold.  Thus, thres = maxt - 0.1*maxt

10) Display the result with pixel values above the threshold ONLY.
figure, imshow( realboth > thres);

You should see a blank page with bright dots on it. You may need to play a little bit with the threshold value to fine tune the image. The bright dots locate where that letter has appeared in the text.

Try the second letter. HAVE FUN.  Submit the list of commands and the resulting images with your file.

Post-lab - Due Wed Nov. 1

Problem1. In MATLAB we have a function imrotate that rotates an image by a given angle.  For this problem, you will try different rotation angles between 0 to 180 to find out the effect of rotation on Fourier Transform.  We will use the rectangle image (Figure 1) that we used in the lab to complete this lab.

f = zeros(512,512);
f(245:265,235:275)=1;
imshow(f)

Perhaps the best way to complete this post-lab is modify the temp.m script.

Produce the centered Fourier Transform of the rectangle after
1.1) 45 degree rotation,
1.2) 90 degree rotation
1.3) 135 degree rotation, and
1.4) 180 degree rotation.

The syntax for imrotation is:

B = IMROTATE(A, ANGLE, METHOD)

Where A is the given image, ANGLE is the counter-clockwise rotation angle in degree, and METHOD is the interpolation method.  As you may have noticed when you rotate an image by an angle, let's say 45, the edge of the shape may become the diagonal of some pixels.  A pixels may only have one gray value representing it.  Thus, you will settle for a value that best describes the pixel.  There are several METHODS to determine such a value for a questionable pixel. Three are used by MATLAB:
1)  'nearest'  (default) nearest neighbor interpolation
2)  'bilinear' bilinear interpolation
3)  'bicubic'  bicubic interpolation

For now you can use the second method bilinear interpolation.

Important: Also, when you rotate the 512x512 image the size may no longer be 512x512 depending on the rotation angle.  For example for 45 degree you will not produce a 512x512 image.  In such a case you need to cut the center 512x512 part of the image to use in your transformation.  The best way to make sure the image is correct before you apply the Fourier transform is to display it first.  For example, display the 45 degree rotated image and then cut the 512x512 center part.

Example (if you rotate the 512x512 image by 45):
g = imrotate(f, 45, 'bilinear');
size(g)     % This gives the size of g
725  725

So the size of this image f after a 45 degree rotation is 725.  We need to make square of 512x512 from this before we apply the Fourier transform.

If you run the imshow on the rotated image, you will get a warning:

>> imshow(g)

Warning: Image is too big to fit on screen; displaying at 56% scale.
In truesize>Resize1 at 308
In truesize at 44
In imshow at 161
A warning is displayed about the dimension.

After I cut the center 512x512 pixels, and used imshow to display the result, I got this image. Note: you want to cut the center 512x512 image.