School of Computer Science
& Software Engineering
Computer Vision CITS4240

Laboratory 4

Aim

The aim of this laboratory session is to introduce you 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.

Building up an Image from its Fourier Components

Taking the Fourier Transform of an image splits the image up into different components, one for each frequency value. If we add the components back together again we reconstruct the image.

It is a bit like splitting a colour image up into different colour components, one each for the red, green and blue components. If we add the colour components back together again we reconstruct the image.

When we split an image up into components (colour or frequency) what we are interested in is how 'big' each component is (eg. how 'red' is a point? - or what is the magnitude and/or phase of a frequency component?)

What makes things awkward is that when we split an image into its different frequency components there are lots and lots of them. Displaying the properties of each component is usually done in the form of an image - the spectrum, and interpreting this image/spectrum takes a bit of getting used to. The position of a value in a spectrum does not correspond to any position in the original image. Instead the position indicates the sine wave frequency that the value is providing data on.

Each element in the Fourier transform represents the magnitude and phase of a basis function. We can reconstruct an image from its transform if we take each Fourier component, multiply it by its corresponding basis function, and then sum them all.

Fourier components come in pairs. For a real valued signal, the Fourier component at a frequency ω is F(ω), which is the complex conjugate of the Fourier component, F(-ω), at -ω. (Note: Taking the complex conjugate simply involves negating the imaginary component.)

In the case of our 2D images where we have a 2D transform F(ω,ν) = F(-ω,-ν)^*, where the superscript ^* indicates the complex conjugate.

If you take the Fourier transform of your image and zero out all the components except for one component pair, F(ω,ν) and F(-ω,-ν), and invert this you will obtain a partial reconstruction of your image that corresponds to this component pair. For the 1D case, it can be easily shown that the partial reconstruction is a linear combination of the basis sine and cosine functions (see also the Matlab demo).

A program freqcomp has been written to (hopefully) give you a feeling for what goes on. From within MATLAB load your image of yourself, convert it to greyscale, and halve its size using imresize. (freqcomp is happier running with small images)

 >> im = imresize(im,0.5);

Then invoke freqcomp

  >> freqcomp(im, 50);

this tells freqcomp we want to look at the first 50 Fourier components of this image.

The program will display:

The image.
The Fourier transform (spectrum) of the image with two red dots marked on it. These are a pair of Fourier components (Fourier components come in pairs placed symmetrically about the centre of the spectrum)
The linear combination of the basis functions that corresponds to the Fourier transform pair in the image above.
The reconstruction of the image (which will look nothing like it at the moment).

If you hit any key the next pair of Fourier components will be shown, along with its sine wave pattern. This sine wave is added to the previous one and displayed as the reconstruction.

If you hold a key down you will start to see the image sequentially built up from its component sine waves. The Fourier components close to the centre of the spectrum correspond to the low frequencies. As you move out the frequencies get higher and the detail in the image starts to emerge. Try running the following

  >> freqcomp(im, 1000, .01);

This runs an animation of the process, the last parameter is the time delay between successive reconstructions.

Filtering images in the frequency domain corresponds to blocking out or letting through specific ranges of frequencies in the Fourier transform of the image. What you see in the animation is the effect of starting with only the lowest frequencies in the image and successively letting through the higher frequencies.

Spatial Convolution

We start off by introducing spatial convolution. A test image with a simple intensity pattern, called 'test2.png', has been placed in the Images directory. Download this image, we will use it to see the effects of different kinds of filtering.

Look at the profile through the image. Assuming that your image is displayed in 'Figure 1' type the following

>> figure(1), improfile       % The figure(1) command makes it the 'active' figure

With the left mouse button click on a point on, say, the left hand side of the image. Then click on the right-side of the image with the right mouse button to digitise a line across the image. A new window will pop up displaying a cross section through the image corresponding to the two end points that you selected in the image.

You will see that the section has a centre step section surrounded by a ramp, which in turn is surrounded by a spike. If you look at the image you may notice a dark band and light band at the bottom and top points of the ramp. These perceptual artifacts are known as Mach bands - if you look at the image profile you will see that they are not actually there.

Blurring

You will find this operation on most image processing programs. On xv you can find this under the 'Algorithms' menu on the main control panel. Typically what is meant by blur is spatial averaging - this leads to the specification of a window size (usually an odd number so that there is a 'middle pixel').

Under MATLAB we can create an averaging filter of size 21x21 as follows:

>> avfilter = fspecial('average', [21 21]);   % Do a help on fspecial

The filter avfilter is simply a 21 x 21 uniform array of values, these values are scaled so that their sum is 1.

You can look at this filter as an image

>> imagesc(avfilter);      % imagesc rescales all values to the range 0-255

Or look at it as a surface

>> surf(avfilter);

This particular filter is not very exciting to look at because it is flat.

The filter is applied to an image 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 image is replaced with the average of all the pixel values within the surrounding 21 x 21 pixel grid.

This process of moving a filter 'over' an image is convolution (Actually strictly speaking the process described above is correlation but the two are equivalent if the filter is symmetric).

Under MATLAB We can then apply this filter by spatial convolution to the image as follows:

>> newim = filter2(avfilter,im);   % do a help on filter2
>> imagesc(newim);                 % use 'imagesc' or 'show' on all floating point data
                                   % (newim will be floating point, imshow will misbehave)

Display the new image and plot a slice through it. Notice how the features have changed. The step changes to a ramp (producing Mach bands - an artifact that really should not be there). The cross section of the line (a spike) changes to match the shape of the filter, a rectangle. The ramp becomes 'rounded' in section

Gaussian smoothing

We can now repeat the process above using a Gaussian filter. Create a Gaussian filter of size 25x25 with a standard deviation of 4 using the command

>> gfilter = fspecial('gaussian', [25 25], 4);

(Note that the standard deviation controls how broad the gaussian shape is. You need to make the size of the filter about 6 times the size of the standard deviation to get its full shape).

Visualising this filter using 'show' or 'surf' will be more rewarding. Try using

>> surfl(gfilter), shading interp, colormap(gray)

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. Apply this filter to the original test image and look at its cross section. You will notice that the smoothing is smoother! The cross section of the line (a spike) changes to match the shape of the filter, now a Gaussian. The Mach bands that were created using the averaging filter on the step should not be present on this image.

MATLAB's Fourier transform routines

The MATLAB function call y = fft(x) calculates the discrete Fourier transform of a vector x of data, using the fast Fourier transform algorithm.

The function ifft defines the inverse Fourier transform. For any x, the value ifft(fft(x)) equals x to within roundoff error.

The function fft2 calculates the 2D fast Fourier transform.

The function call z = fftshift(y) shifts the data in y around so that the 0 frequency point is moved from the 'ends' of the array to the middle. The array y can be 1 or 2D.

Convolution in the Frequency Domain

Make a Gaussian filter, this time make its size match that of the test image (256 x 256 pixels), though still with a standard deviation of 4. (The standard deviation is the 'active ingredient' that controls how broad the Gaussian shape is and specifies what effect the filter will have - the overall size is secondary)

 >>  gfilt = fspecial('gaussian', [256 256], 4);

If we apply this filter to an image using spatial convolution we will smooth and blur the image. Fine detail, corresponding to high frequency information, in the image will be removed. We can see why this is the case if we look at the Fourier Transform of the filter.

We generate the 2D Fourier Transform of the filter in MATLAB using the command

 >> gfiltfft = fft2(gfilt);   % Use fft2 for 2D data, fft for 1D data

To look at the Fourier Transform we typically only display its amplitude. The MATLAB function abs, when applied to complex values, returns their magnitude.

 >> imagesc(abs(gfiltfft)), colormap(gray)

Again make sure you use the imagesc or show function which rescales an array of floating point values to the range 0-255 and displays them. Note that data corresponding to the low frequencies are positioned around the edge of the displayed image. Use the fftshift function to move the 0 frequency point to the middle.

The cross-section of the FFT of a Gaussian filter will also be a Gaussian (though, in general, with a different width).

What does this mean?
The Fourier Transform of a filter tells you which frequencies in the original image will be 'let through' by the filter.

The FFT of this Gaussian filter is 'high' or 'bright' in the middle indicating that it will 'let through' the low frequency parts of the image. Away from the middle the FFT of the filter is zero, indicating that these (high) frequency components in the image will be eliminated.

Now apply the filter to the test image. We read in the test image and take its Fourier Transform.

 >> im = imread('test2.png');    
 >> imfft = fft2(im);

If we look at imfft

 >> imagesc(fftshift(abs(imfft)));

all you will see is a small dot in the middle of the image. This is because in most images the magnitudes of the low frequency components are very large relative to the high frequency components, and totally dominate what is displayed. To make things easier to see we typically take the logarithm of the Fourier amplitudes before displaying them. This is like making a gamma adjustment to a very dark image in order to bring out the dark portions.

It is worthwhile adding a small constant to the magnitude of the Fourier transform before taking the log to avoid log of 0 problems, log of 0 tends to minus infinity and this can muck up the display of your Fourier transform!

 >> imagesc(fftshift(log(abs(imfft)+eps)));

(eps is a predefined value in MATLAB giving you the effective machine precision. It is defined as the smallest value you can add to 1 that will give you a value different from 1.)

Use improfile to pick out the cross section of the Fourier transform.

Now can apply the filter to the image by multiplying the FFT of the image and the FFT of the filter together with the following command.

 >> newimfft = imfft.*gfiltfft;    % Note you must use '.*' to multiply the
                                   % two together.

This last operation takes each pixel in the FFT of the image and multiplies it by the corresponding pixel in the FFT of the filter. Note that MATLAB is doing point-wise complex multiplication for you here. This produces a new Fourier Transform, 'newimfft' - the Fourier Transform of the filtered image. Look at newimfft (the log of its amplitude). The line that went across the width of the FFT will change:

Where the FFT of the filter is zero, newimfft will also be zero.

Where the FFT of the filter is large, newimfft will roughly correspond to the FFT of the original image.

In doing this we are using the convolution theorem to perform spatial convolution via a simple multiplication in the frequency domain.

To see the convolution result we have to apply the inverse Fourier Transform to 'newimfft' to move ourselves from the frequency domain back to the (normal) spatial domain to obtain 'newim'.

 >> newim = fftshift(real(ifft2(newimfft)));

Here the steps are:

Take the inverse 2D Fourier transform of newimfft.
Extract the real part of this. The inverse fft operation returns a complex valued array of values. In theory our image should be real-valued, in practice one finds that very small imaginary components do exist due to roundoff errors.
Perform a quadrant shifting operation to put all the bits of the image where they should be. Why do we need to do this? Remember that as far as the Fourier transform is concerned the 'origin' of the data is at the edge of the array of data and not in the centre. Thus ideally when we constructed the Gaussian filter we should have quadrant swapped it before taking its Fourier Transform. This will have moved the origin of the Gaussian (the centre) to the edge of the 2D array. This final quadrant swapping done above compensates for us not applying the quadrant swapping to the filter in the first place.

The result should look just like the blurred image you generated via spatial filtering earlier on. (Actually the results will be a bit different around the edges because the way the boundaries of the image are treated is different between the two techniques of filtering)

If you do not believe this you can apply 'gfilt' via convolution in the spatial domain.

 >> newim2 = filter2(gfilt, im);

Actually if you try this it will take a long time to compute because gfilt is very large (256 x 256 pixels). We had to make it the same size as the image (test.png) to allow us to multiply the Fourier Transforms of gfilt and the image together (they must be identical in size).

Most of this filter's values are 0 so we can construct a smaller filter (but with the same standard deviation of 4) of size 25 x 25 that achieves the same effect - but in much less time.

The Fourier Transform of an averaging filter

Construct an averaging filter, but instead of using the 'fspecial' function use the following commands to construct a filter with the same size as the image (256 x 256 in this case), with a (say) 21 x 21 square averaging region in the middle

 >> avfilt = zeros(256);                   % An array of zeros.
 >> avfilt(118:138, 118:138) = 1;          % Fill a centre region with ones.
 >> avfilt = avfilt/(sum(sum(avfilt)));    % Normalise so that sum of values is one.

Look at the cross section of the amplitude of the Fourier transform of the averaging filter - it will look like the path of a bounding ball. Where the section 'bounces up' actually corresponds to regions where the real part of the Fourier Transform is negative.

The fact that this FFT has some negative values means that when we apply the FFT of this filter to the FFT of the test image by multiplying the two together we will end up negating/reversing the phase of some frequency components in the image. This corruption of some phase values in the image interferes with our perception of the image. This is responsible for producing the Mach bands on the step edge when before there were none.

In general, when we are seeking to preserve perceptually important features in the image, we should use filters that have an FFT where all the values are 0 or greater. This means that while some frequency components might be reduced to zero, none are actually reversed, or phase shifted by some other angle.

Note that if you do not centre the averaging filter properly the convolution result will produce an image wrapped around on itself. Try it if you like!

Phase-only Images and Phase Swapping

Select your image from the ClassOf2007 images directory and some other image (If you do not have your own image there pick any two images from the class of a previous year.). Ensure that the two images that you have chosen are of the same dimensions. Convert the images to grey valued ones

>> greyim = rgb2gray(im);

Calculate the 2D FFT for each image and extract the phase and amplitude data from each of the FFT images (you can get the phase information from an array of complex values using the MATLAB function angle). Now swap the phase information from the two images. Apply an inverse Fourier transform to get back to the spatial domain (you will want to take the real part of this inverse).

Remember each complex valued Fourier component is given by

F(ω) = magnitude * (cos(phase) + i*sin(phase))

Note that i (and j) is a predefined variable in MATLAB to represent sqrt(-1).

Rather than working with the sine and cosines of the phase angles you can note that

(cos(phase) + i*sin(phase)) = F(ω)/magnitude

Thus you 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 the relative sizes of the real and imaginary components.

If you take the inverse fft of this array of unit magnitude, complex valued, phase data, then you will construct what is know as a phase-only image. Here you will see a ghostly, but recognizable, version of the original image. It will look very much like a high-pass version of the image because the high frequency components will have much higher magnitudes than in a normal image. This type of image is called a phase-only image because all magnitude imformation has been removed.

Portfolio Requirements

The following should be included in your portfolio:

A partial reconstruction of your image using freqcomp. Note the number of frequency components that you used in the reconstruction.
An image of yourself, your FFT, the FFT of a Gaussian smoothing filter, the product of your FFT and the filter, and the final smoothed image. (Choose a Gaussian filter standard deviation other than 4 pixels)

Note: There are some small catches in saving images of floating point values. See the notes on saving images and plots at the end of the document How to construct your portfolio.

A sequence of images as above but using an averaging filter.
An image of your two face images and the two phase swapped images.
A phase-only version of your image.
A listing of the MATLAB commands you used to do the phase swapping (or the function/script you might have written to do this). These Matlab commands should be saved to appropriate script files and included in your portfolio so that I can test them on your images, which should be submitted also with your portfolio.
Do not forget to include all other Matlab functions (excluding the Matlab library functions but including Dr Peter Kovesi's) that are needed by your Matlab functions/scripts.

You might find it more tidy to save your MATLAB code and input/output test image files into subdirectories away from your HTML file. Ensure that all the links in your HTML file are correct. You should also provide some explanation in your HTML file to demonstrate your understanding about the problem and the work that you carried out.

Return to Computer Vision 4240

School of Computer Science & Software Engineering Computer Vision CITS4240