School of Computer Science
& Software Engineering
Computer Vision CITS4240

Laboratory 5

Aim

The aim of this laboratory session is to reinforce the concepts associated with the Fourier Transform. This will involve constructing filters directly in the frequency domain in order to implement some image enhancement techniques.

Constructing Filters in the Frequency Domain

One way to process an image is to devise a spatial filter, look at its Fourier spectrum to see if it has the desired characteristics and, assuming it is appropriate, apply the filter by convolution.

The alternative approach is to devise a filter directly in the frequency domain, apply it to the Fourier transform of the image by point-wise multiplication, and then apply the inverse Fourier transform to see the final result.

The Butterworth Filter

The Butterworth filter is simply a convenient equation that allows you to specify a 'bump', of height 1, centred on the origin.

                  1.0
  F  =   ______________________   ,
                           2*n
         1.0 + (r / cutoff)

where r is the distance from the origin, the radius of the bump is specified by cutoff and the sharpness of its boundary is controlled by the parameter n (a positive integer).

This 'bump' centred on the origin acts as a low-pass filter - letting through low frequencies in the image and blocking out the high frequency components. The cut off radius specifies the cut off frequency - the point at which frequencies start being 'cut off'.

Notes:

A high-pass filter can be constructed from a low-pass filter by simply subtracting a low-pass filter from the constant 1.
A high-boost filter is simply a high-pass filter with a constant offset added to it. Unlike a high-pass filter, which blocks out all the low frequency components, the offset value in a high-boost filter means that all the frequency components in the image get transferred to the output to some degree, and in addition the high frequency components get amplified.

Download the functions lowpassfilter.m, highpassfilter.m, and highboostfilter.m. Construct a few low-pass, high-pass and high-boost Butterworth filters and experiment with different values of radius and order. To view the filters that you generate, use the imagesc or show functions and use improfile to check their cross sections. For 'nice' 3D plots, use the MATLAB function surfl and use the shading interp option. Note that these functions return filters with the 0 frequency point at the corners, use fftshift to view them. Note also that the highest frequency that you can specify for the cutoff frequency is 0.5 -- the Nyquist frequency, corresponding to a waveform with a wavelength of 2 pixels.

For instance, the Matlab statements below call the lowpassfilter function to generate a low-pass filter and then display it using the surfl function.

  >> F = lowpassfilter(size(im), cutoff, n);
  >> surfl(fftshift(F)), shading interp;

Perform low-pass, high-pass and high-boost filtering on your image by taking the Fourier transform of your image, point-wise multiplying it with the filter, and then back transforming to get the result.

You should also calculate the inverse Fourier transforms of the filters you used so that you can see the shape of the filter in the spatial domain. Note that when you take the inverse Fourier transform you should be extracting the real component of the result (if the imaginary component is not nearly zero then you have a problem). You will also need to do a 'fftshift' on the result to move the centre of the filter to the middle of the image. I suggest that you view the filter as a surface to see its shape properly.

'Black-spot' filtering

Sometimes one wants to 'black out' small spots within a Fourier spectrum

Load up image 'finger3.png'. (thanks to Bruce Comber of the Australian Federal Police for these images). This is a fingerprint on a half-toned magazine picture. Look at its Fourier Transform - you will see the repeated pattern of the underlying half-tone pattern appear as a series of bright dots in the FFT spectrum. These bright dots represent all the harmonics of the repetitive half-tone pattern. Each bright dot tells you that that frequency component is very strong. If we could black out these bright spots in the Fourier spectrum and then invert back into the spatial domain, we would get rid of this half-tone pattern.

Write a function that meets the following specification


% BLACKSPOT
%
% Function returns a high-pass Butterworth filter centred at position
% (xc,yc) in an image
%
% usage: f = blackspot(im, cutoff, n, xc, yc)
% 
% where: im      is the image for which the filter is to be used on.
%        cutoff  is the cut off frequency of the filter.  The range
%                of values allowed is [0,0.5].
%        n       is the order of the filter, the higher n is the sharper
%                the transition is.
%        xc,yc   is the desired centre of the filter.

function f = blackspot(im, cutoff, n, xc, yc)

Copy the lowpassfilter code and adapt it so that instead of offsetting the x and y matrices by fix(rows/2)-1 etc. you should use the values xc and yc to shift the origin to the specified position. Note also that you want to change the low-pass filter equation so that it forms a high-pass filter. Also, this function should not quadrant shift the filter before returning it.

Now write a function

% EDITSPECT -  Remove periodic interference patterns in an image.
%
% This function allows one to edit the amplitude spectrum of an image by
% 'blacking out' points in the spectrum with the mouse.  
%
% Usage:  newim = editspect(im, cutoff,n)
%
% Arguments:   im      - The image to be processed.
%              cutoff  - The radius of the blackout spot generated by a mouse
%                        click.  A typical value is in the range [0.02,0.1].
%              n       - order of the Butterworth function used to construct the spot.
%
% Returns:  newim      - The image reconstructed with the new amplitude
%                        spectrum.

function newim = editspect(im, cutoff, n)

This code should:

calculate the Fourier transform of the image and quadrant shift it,
display the log of the magnitude of the Fourier transform, and then
enter a while loop that allows one to repeatedly digitise points in the image of the Fourier spectrum blacking out points as we go.

You would find the MATLAB function ginput useful for getting the coordinates of the position where the mouse button was pressed. The function also allows you to identify which mouse button was pressed: left = 1, middle (left + right) = 2, and right = 3. For instance, your code could be so constructed that the clicking of the left mouse button would mean continuously collecting the digitised points in a while loop and the clicking of the right mouse button would mean terminating the loop.

You can halve the number of points you have to digitise by exploiting the symmetry that occurs in the Fourier transform. If you have a spot to be blacked out at coordinates (x,y) relative to the zero frequency point, there will be an identical point at coordinates (-x,-y).

Inverse Filtering

Load the image 'blur.png' from the Images directory. This is an out-of-focus image of some text. We can model this out-of-focus image as being the convolution of the ideal, sharp image with an averaging filter of some kind. If we can undo this convolution then we can reconstruct a sharp image.

The filter that represents this blurring process is the point spread function. It is literally the 'shape' that an ideal point in the sharp image would be spread into in the blurred image.

This convolution with the point spread function can be inverted using the convolution theorem. If F and G are the Fourier transforms of the sharp and blurred images respectively, and H is the Fourier transform of the point spread function, we can write the following:

   G = F.*H

If we want to calculate, F, an estimate of the Fourier transform of the sharp image, then you might think we can do the following

   F = G./H

We can think of 1/H as being the inverse filter.

Unfortunately there are three problems to be solved:

We do not know what the point spread function is, and hence we do not know H.
Even if we do know H, it is likely that some values in H will be very small or zero, making 1/H undefined.
It is also likely that the blurred image had some noise added to it in the imaging and discretisation process. Noise typically has a fairly uniform amplitude spectrum and can be the dominant component of the high frequency components of an image. Where H is very small it is likely that the noise component will be greatly amplified in the reconstructed image.

One way of approaching these problems is to use the Wiener filter. In its simplest form, a Wiener filter, W, can be constructed by multiplying the inverse filter 1/H by a modifying factor:

           2
        |H|
    ____________
         2
    ( |H|  + K )

giving

       
                    2
         1       |H|
    W = ___  ____________,
                  2
         H   ( |H|  + K )

where the term K is a constant, representing the ratio of the noise power spectrum to the signal power spectrum.

Note that when the magnitude of H is small or zero, the modifying term above tends to zero, making the overall filter value zero (rather than infinity). On the other hand, when the magnitude of H is large relative to K, the modifying term tends to one, returning the overall filter value to the ideal of 1/H. Thus, the modifying term in the Wiener filter zeros out the magnitude of the inverse filter at the locations where otherwise it would go off to plus or minus infinity. This means, however, that the original image spectrum in these regions can never be restored.

Similar to the low-pass, high-pass, and high-boost filters, the Wiener filter is defined directly in the frequency domain. One may therefore be interested to inspect the shapes of these filters in the spatial domain. To do so, you can apply the inverse Fourier transform (i.e., the Matlab's ifft function) to these filters. The cross sections (i.e., the one-dimensional profiles) of these filters can be visualised using the improfile function. Alternatively, one-dimensional low-pass, high-pass, high-boost, and Wiener filters can also be constructed.

In practice, K is not a constant but varies across the spectrum. Sophisticated Wiener filtering algorithms attempt to infer K automatically and make it vary across the spectrum.

Your task...

Restore the image 'blur.png' with a Wiener filter. You will need to guess and generate a suitable point spread function, and you will need to determine a suitable value of K. I suggest that you use the function psf2.m to generate point spread functions of varying sizes and sharpness of shape. Carefully read the help information for psf2.m. Note that the parameters to this function allow one to generate elongated and angled point spread functions. For the 'blur.png' image we would expect the point spread function to be round. Thus, the following parameter values to psf2.m are recommended: the 'sze' parameters should be matched to the size of the image; the 'lngth' and 'width' parameters of the point spread function should be the same; the 'sqrness' parameter should be 2; the 'ang' parameter can be set to 0.

Write a function to the following specification

% WIENFILT
%
% Usage:  newim = wienfilt(im, leng, width, ang, order, K)
%
% Arguments:   im      - The image to be processed.
%              leng    - The length of the point spread function to use.
%              width   - The width of the PSF
%              ang     - angle of the PSF
%              order   - order of the Butterworth function used to construct
%                        the point spread function. 
%              K       - The estimated noise/signal power ratio.
%
% Returns:  newim      - The reconstructed image

Your code should construct a point spread function matching the size of the image using psf2.m and the parameters supplied. This should be quadrant shifted prior to taking the Fourier transform to calculate H. The quadrant shifting is required because as far as the Fourier transform is concerned the origin is at the top left corner. So to 'centre' the point spread function we need to quadrant shift it. If you do not do this you will produce a quadrant shifted final image.

Once you have written this function you will need to experiment by trial and error to find the parameters that will restore the image as best as possible. As a guide, the radius of the point spread function should be somewhere between 2 and 6 pixels (this is affected by the order of filter that you choose). The value of K should be kept as small as possible. The smaller you can get it the sharper the reconstruction, but the greater the noise amplification. You usually have to put up with a reasonable amount of noise amplification to get a reconstruction you can read. Try a value of K around 0.01 to start with. Later try reducing K to 0 to see what happens.

You should display the real part of the point spread function's Fourier transform and the real part of the Wiener filter that was generated (the imaginary parts should be approximately 0). You should also display their point-wise product to see how well it is able to restore the spectrum of the image. Ideally the product should be one everywhere; however, where the modifying factor of the Wiener filter becomes significant the product will fall to zero -- these are regions of the image spectrum that can never be restored.

Now try to deblur the number plate on the motion_blurred.png image...
You will need to generate an elongated point spread function at the appropriate angle. I suggest you display the image and then, with 'pixval on', determine the start and end points of one of the corners of the blurred image of the number plate to determine the length and angle of the point spread function.

Portfolio Requirements

The following should be included in your portfolio:

Listings of your blackspot, editspect and wienfilt code and any Matlab scripts that you may have written for testing and generating figures. Ensure that you submit all your .m files as normal text files (do not reformat them into html or Word), as they will be tested and run under Matlab.
The low-pass, high-pass, and high-boost versions of your own image. Include the following and provide brief explanation/comments to each of them:
- the surface plots of the filters that you used in both the frequency domain and the spatial domain.
- the parameters that you used to produce each image.
- the results that you obtain when you use a high order and a low order of the Butterworth filter. By high order I mean a value of n that is over 50.
- the behaviour of high-boost filtering when the boost is very large.
Images showing your best efforts at enhancing the finger3 and pidgie images using 'black spot' filtering. Note the parameters that you used.
For each of the blur and motion_blurred images, include
- images of the point spread function and the real part of its Fourier transform (the imaginary part should be approximately zero).
- image of the real part of the Wiener filter that was generated (the imaginary part should be approximately zero).
- plots of the cross sections of the real part of the Fourier transform of the point spread function and the real part of the Wiener filter and their point-wise product. You should also include the cross section of the inverse Fourier transform of the Wiener filter.
- the resultant image showing your best effort at enhancing the image. Note the parameters that you used.

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 412

School of Computer Science & Software Engineering Computer Vision CITS4240