|| Lab2 || Lab3 || Lab4 || Lab5 || Lab6 || Lab7 || Lab8 ||
 Back to Main Menu

________________________________________________________________________________________________________________________

Aim of this lab:

The aim of this laboratory session is to introduce us to MATLAB's Fourier transform routines, learn how to interpret Fourier transform data, and then to see how convolution is implemented in the frequency domain using simple averaging filters as an example.
________________________________________________________________________________________________________________________

The first thing we did for this lab is to create an partial reconstruction of our own image from Fourier component. A function has been provided to us in order to do this operation, that is freqcomp.m. Since this function is happier running with small images, we than resize our image to half the size.

MATLAB Command:
im = imread('me.jpg');
g = rgb2gray(im);         %convert to grey scale
im = imresize(im,0.5);    %half the image size

Executions:
freqcomp(im,50);      %frequency = 50, get the first 50 Fourier component of this image
freqcomp(im,500,0.1)  %frequency = 500, delay = 0.1 , get the first 500 Fourier component of this image.
freqcomp(im,2000,0.1) %frequency = 2000, delay = 0.1, get the first 2000 Fourier component of this image.

I've run this function with various frequencies which include 50,500 and 2000. To make things more convenient, I have executed the function with an extra parameter, that is the delay between successive reconstructions which is set to 0.1. The results of the executions are as follows:

       
Figure1: The first 50 Fourier components of my image                  Figure2: The first 500 Fourier components of my image

Figure3: The first 2000 Fourier components of my image

Spatial Convolution

The next task we did in this lab to experiment the spatial convolution with 2 different filters on 'test2.png'. The filters are averaging filter and Gaussian filter.

MATLAB Command:
test = imread('test2.png');     %load 'test2.png'
test(1), improfile              %get the profile of 'test2.png'

       
Figure4: 'test2.png'                                                              Figure5: Image profile of 'test2.png'

In the image itself (Figure4), we can notice that there is a dark band and light band at the bottom and top of the ramp. These are known as the Mach bands. However, if we look at the profile of the image (Figure5), these features does not exist.

Blurring

The first spatial convolution we are going to experiment is blurring or also known as Spatial Averaging. We have created an averaging filter of size 21x21 (Figure6). This filter is simple a 21x21 uniform array of values, and these values are scaled so that their sum is 1. Below are the figures of the 21x21 averaging filter and also the figures of 'test2.png' after applying the averaging filter to it:

MATLAB Command:
avfilter = fspecial('average',[21,21]);     %create a averaging filter of size 21x21
surf(avfilter);                             %display the averaging filter as a surface (Figure6)
newtest1 = filter2(avfilter,test);           %apply the averaging filter on 'test2.png'
show(newtest1);                      %display the new 'test2.png' after applying the averaging filter (Figure7)
testim(1),improfile;                 %get the image profile of 'test2.png' after averaging

       
Figure6: Averaging filter of size 21x21                               Figure7: 'test2.png' after applying the averaging filter


Figure8: Image profile of 'test2.png' after applying the averaging filter

The filter is applied to 'test2.png'  by placing it 'over' each location in the image, multiplying the image pixel values under each grid point of the filter by the value of the filter at that point, adding up all the results, and assigning the answer to the pixel under the centre location of the filter. The filter is then indexed to the next point in the image and the process repeated. Thus each pixel in the images is replaced with the average of all the pixels values within the surrounding 21x21 pixel grid.

Gaussian smoothing

Next, we go through the next type spatial convolution that is the Gaussian smoothing. By using the same image as above ('test2.png'), we now apply a Gaussian filter on it. First, we create a Gaussian filter of size 25x25 with a standard deviation of 4. The standard deviation controls how broad the Gaussian shape is. Following are the results after applying the Gaussian filter to 'test2.png':

MATLAB Command:
gfilter = fspecial('gaussian',[25,25],4);       %create a Gaussian filter of size 25x25, standard deviation = 4
surfl(gfilter), shading interp, colormap(gray)  %display the Gaussian filter
newtest2 = filter(gfilter,test);                %apply the Gaussian filter on 'test2.png'
show(newtest2)                       %display the new 'test2.png' after applying the Gaussian filter (Figure10)
newtest2(1), improfile               %get the image profile of 'test2.png' after averaging

       
Figure9: Gaussian filter of size 25x25, std dev = 4             Figure10: 'test2.png' after applying Gaussian filter


Figure11: Image profile of 'test2.png' after applying the Gaussian filter

This filter will perform a weighted averaging of pixel values within the filter window, with the centre pixels having more weight in the averaging operation. After applying the filter to the original test image, we can notice that the smoothing is smoother compare to the average filter. The cross section of the line (a spike) changes to match the shape of the filter, now a Gaussian. Notice that the Mach Bands that were created using the averaging filter are disappeared now.

Fourier Transform

The next section of this lab is dealing with the Fourier Transform. Firstly, we will apply a Fast Fourier Transform (FFT) to our image and then multiply with two different types of filter that is the Gaussian Smoothing Filter and the Averaging filter. The product of the Fourier Transform of my image and the filter will create a Fourier Transform of a new image. Then we will take this result and apply the Inverse Fourier Transform to obtain the new image.

Applying Gaussian Filter:-

First of all, we take the FFT of our image from the grey scale image then we must get the size of the image because the filter has to be the same size as the image. The size of the image is 276 x 368. Therefore, we create a Gaussian Filter of the size 276 x 368 with a standard deviation of 4. All of this process are executed with the following command:

MATLAB Command:
im = imread('me.jpg');
g = rgb2gray(im);
imfft = fft2(g);
imagesc(fftshift(log(abs(imfft)+eps))), colormap(gray);   %display the FFT of my image
size(imfft)                                               %get the size of my image

%create the Gaussian Filter
gfilter = fspecial('gaussian',[276,368],4);
surfl(gfilter),shading interp, colormap(gray);         %display the Gaussian Filter

%get the FFT of the filter
gfft = fft2(gfilter);                                  %get the 2D Fourier Transform of the Gaussian Filter
imagesc(fftshift(log(abs(gfft)+eps))), colormap(gray); %display the FFT of the filter (Figure15)

%multiply the FFT of the image and the FFT of the filter
newfft = gfft.*imfft;     
imagesc(fftshift(log(abs(newfft)+eps))),colormap(gray); %display the product newfft

%apply inverse Fourier Transform to get from frequency domain back to spatial domain (normal)
newimfft = fftshift(real(ifft2(newft)));
show(newimfft);       %display new image (Figure 17)

        
Figure12: My grey scale image                          Figure13: FFT of my image                            Figure14: The Gaussian Filter with std dev = 4

           
Figure15: FFT of the Gaussian Filter                Figure16: Product of the FFT of my image     Figure17: My new image after applying inverse FFT
                                                                                         and Gaussian Filter

After that, another test with the same approach as above has been carried out. But this time we will apply a Gaussian Filter with a standard deviation of 10 instead of 4. The result of this test is as follows:


Figure18: New image after using Gaussian Filter
               with std dev = 10

Applying Averaging Filter:-

Here is the result of my image after applying an averaging filter. In order to create the averaging filter, we first get the size of the image and apply the following command:

MATLAB Command:
size(imfft)    %get the size of the image

%create the Averaging Filter
avfilter = zeros(276,368)                  %An array of zeros
avfilter(128:148,174:194) = ones(21)       %Fill a centre region with ones
avfilter = avfilt.(sum(sum(avfilter)));    %Normalise so that sum of values is one

%get the FFT of the filter
avfft = fft2(avfilter);                    
imagesc(fftshift(log(abs(avfft)+eps))),colormap(gray); %display FFT of the filter (Figure21)

%multiply the FFT of the image and the FFT of the filter
newfft = avfft.*imfft;                               
imagesc(fftshift(log(abs(newfft)+eps))),colormap(gray); %display the product newfft (Figure22)

%apply inverse Fourier Transform to get from frequency domain back to spatial domain (normal)
newim = fftshift(real(ifft2(newfft)));
show(newim);    %display the new image (Figure23)

       
Figure19: My grey scale image                          Figure20: FFT of my image

           
Figure21: FFT of the Averaging Filter              Figure22: Product of the FFT of my image     Figure23: My new image after applying inverse FFT
                                                                                         and Averaging Filter

Phase and Amplitude Swapping

In this section, we will take 2 different images and swap their phase and amplitude. I've took 2 different pictures of myself which has already being resize since both images must be the same size. Below are the results of the swapping and commands used to do this task:

Image A and Image B are both the original images and after the swapping, the new image A (Figure26) has the Image A's amplitude and Image B's phase information. Whereas, the new Image B (Figure27) has the Image B's amplitude and Image A's phase information.

MATLAB Command:
im1 = imread('resize1.jpg');
im2 = imread('resize2.jpg');
g1 = rgb2gray(im1);
g2 = rgb2gray(im2);
im1fft = fft2(g1);
im2fft = fft2(g2);
im1fftamp = abs(im1fft);
im2fftamp = abs(im2fft);
im1fftphase = im1fft./im1fftamp;
im2fftphase = im2fft./im2fftamp;
newim1fft = im1fftamp.*im2fftphase;
newim2fft = im2fftamp.*im1fftphase;
newim1 = real(ifft2(newim1fft));
newim2 = real(ifft2(newim2fft));
show(newim1);
show(newim2);

       
Figure24: Image A                                                     Figure25: Image B

       
Figure26: New Image A                                    Figure27: New Image B

Remember each complex valued Fourier component is given by
F(w) = magnitude * (cos(phase) + i*sin(phase))

rather than working with the sine and cosines of the phase angles we can note that
(cos(phase) + i*sin(phase)) = F(w)/magnitude

Thus we can extract the phase part as an array of complex values simply by dividing (point-wise) the Fourier transform by its magnitude. This phase part will be an array of complex values all having a magnitude of one and the phase angle being represented by relatives sizes of the real and imaginary components.

Phase-only Images

If we take the inverse FFT of this array of unit magnitude, complex valued, phase data we will construct the phase-only image. The results are as follows:

MATLAB Command:
phaseonly1 = real(ifft2(im1fftphase));
phaseonly2 = real(ifft2(im2fftphase));
show(phaseonly1);    %Figure28
show(phaseonly2);

   
Figure28: Phase-only for Image A                      Figure29: Phase-only for Image B