一次学会Fourier变换图像处理

Fourier Transform Demonstrations

Some demonstrations of how Fourier Transforms relate to image processing.

Contents

• Setup
• Display of 2-D Fourier components
• FT of an image
• Take the FT
• Get the frequency values for the indices of the FT array
• Partially reconstruct the image
• Making a mask to pick out particular frequencies
• Smoothing the image via the Fourier transform
• Differentiating the image via the Fourier transform
• Experimenting yourself

Setup

You can execute the code in this file a section at a time to see the effects. (If you run the whole script, the results will overwrite each other.) You can execute a few lines by copying them into the Matlab Command Window. Alternatively, you can view the script in the editor, and activate Cell Mode from the Cell menu. Then you can click in a cell (a section between two headings) to highlight it and click the "Evaluate cell" icon, or press CTRL-ENTER, to execute it.

First you must set your Matlab path to find the Sussex Matlab computer vision teaching libraries. At the time of writing, the command, which you give at the Matlab command prompt, is

 addpath H:\teach\ComputerVision\matlab

If you are running this as a script on Informatics machines, you can uncomment the line above to execute it. You will also need to execute it at the start of each future session. You can do this by copying it to a startup.m file in the directory (folder) in which you start up Matlab, or you can just type it in or cut and paste it each time.

Display of 2-D Fourier components

A 1-D FT decomposes a signal into the sum of a set of sine/cosine waves. The output of the FT is the set of amplitudes and phases of these waves.

For image processing, we need the 2-D equivalent. The image is broken down into patterns which have a sinusoidal variation along each axis. We can display these patterns by taking the inverse FT of an array that is zero everywhere, except for a single value of 1. Thus we pick out a single Fourier component.

Actually it isn't quite this simple. Fourier transforms of real images (i.e. those whose grey levels are represented by real, not complex, numbers - the usual case) are symmetrical about zero. We have to set the "negative frequencies" as well. In Matlab, the FT has zero frequency at row 1, column 1, and the frequencies wrap round the matrix, resulting in the code below.

In the figure, the component with M cycles on the y-axis and N cycles on the x-axis is at row M+1, column N+1.

s = 128;   % Size of each matrix to display
n = 4;     % number of rows and cols in the figure
for M = 0:n-1           % no. of cycles on y-axis
    for N = 0:n-1       % no. of cycles on x-axis
        F = zeros(s);   % F is the Fourier transform array ...

        % ... which is all zeros apart from a single element ...
        F(M + 1,          N + 1         ) = 1;  % at this point
        % ... and the symmetrical elements for the negative frequencies
        F(mod(-M, s) + 1, N + 1         ) = 1;
        F(M + 1,          mod(-N, s) + 1) = 1;
        F(mod(-M, s) + 1, mod(-N, s) + 1) = 1;

        % Now it's set up, do the inverse FT
        imcomp = ifft2(F);       % f is the Fourier component in the spatial domain

        subplot(n, n, n*M + N + 1);
        imshow(imcomp, []);
    end
end

FT of an image

We read in an image, reduce the size, and display it. We make each dimension a number that has many small factors, as the Fast Fourier Transform algorithm can exploit this to speed up the calculation.

im = teachimage('chess1.bmp');
im = im(1:2:384, 1:2:512);
figure; imshow(im);

Take the FT

We now compute the Fourier Transform of the image and display the amplitudes of the components. (That's what the call to abs gives.) We take the logarithm before display as otherwise the zero-frequency component dominates and we can't see any detail.

The low-frequency components are at the corners of the plot and the high-frequencies are in the middle. There are two mirror symmetries, as mentioned earlier. The bright streak is a result of the high spatial frequencies in the carpet in the background.

F = fft2(im);
figure;
imshow(log(abs(F)), []);

Get the frequency values for the indices of the FT array

To do more, we need a mapping between indices into the array of FT components and the actual spatial frequencies. We generate this here, taking account of the symmetries.

First we find the frequencies in units of "cycles across the image". This is what cxrange and cyrange are. Note that the frequencies become negative as we cross the centre of the image.

Then we convert to "radians per pixel", and store these as fxrange andfyrange. These are the correct units for mathematical manipulation of the FT.

We use meshgrid to produce 2-D arrays containing the frequency values for every Fourier component.

nx = size(F, 2);
ny = size(F, 1);
cxrange = [0:nx/2, -nx/2+1:-1];        % cycles across the image
cyrange = [0:ny/2, -ny/2+1:-1];
[cx, cy] = meshgrid(cxrange, cyrange);

fxrange = cxrange * 2*pi/nx;           % radians per pixel
fyrange = cyrange * 2*pi/ny;
[fx, fy] = meshgrid(fxrange, fyrange);

Partially reconstruct the image

We new reconstruct the image from some low-frequency components, in order to see how they build up an increasingly accurate representation of the image as higher frequencies are added.

Each of the components displayed in the first section is given the correct amplitude and phase for the chessboard image. Each image is made by adding together all the components with vertical frequency less than or equal to M, and horizontal frequency less than or equal to N, measured in cycles across the image. M and N depend on position in the display: they go 0, 1, 2, 4, 8 down the rows or across the columns.

With 8 cycles in each direction, we have enough detail to see the squares on the chessboard.

figure;
cycles = [0 1 2 4 8];
n = length(cycles);
for p = 1:n
    for q = 1:n
       keep = abs(cx) <= cycles(p) & abs(cy) <= cycles(q);
       partF = F .* keep;
       partim = ifft2(partF);
       subplot(n, n, n*(p-1) + q);
       imshow(partim);
    end
end

Making a mask to pick out particular frequencies

We can suppress the high spatial frequencies by multiplying by an array that has large values only for the low spatial frequences. Here, we do this by using a Gaussian function of the frequency.

We display the logarithm of the mask, so that its structure is more visible.

sigma = 0.3;        % Width of the Gaussian, in radians/pixel
ms = exp(-(fx.^2 + fy.^2)/(2*sigma^2));
figure; imshow(log(ms), []);

Smoothing the image via the Fourier transform

We now simply mutliply the FT point by point by this mask, and the take the inverse FT to get an image.

The result is a smoothed image, because the higher spatial frequencies have been suppressed. In fact, it can be shown that this is equivalent to convolving the image with a Gaussian mask.

smoothF = F .* ms;          % multiply FT by mask
smooth = ifft2(smoothF);    % do inverse FT
imshow(smooth, []);

Differentiating the image via the Fourier transform

We can perform other convolutions by multiplying the FT by a suitable mask. It turns out that mutliplying by the x-frequency and also the imaginary unit i is equivalent to differentiating along the x-axis. (This doesn't work for the highest frequency, which has to be suppressed.)

This can also be understood as suppressing the low frequencies (on the x-axis) and emphasising the high frequencies, whilst preserving the symmetries of the FT.

ftd = ft .* fx .* i;    % multiply FT by mask
ftd(:, nx/2+1) = 0;     % suppress top frequency
d = ifft2(ftd);         % do inverse FT
imshow(d, []);

Experimenting yourself

You can experiment with this demonstration yourself, by downloading this html document and using Matlab's grabcode function to extract the original M-file. You can then edit it to change the parameters or try other kinds of image filtering.

This document may be accessible from outside Sussex University, but to use functions from the local library you need to be a student or member of staff of the University.

Copyright University of Sussex, 2006

转载：http://www.cogs.susx.ac.uk/courses/acv/matlab_demos/fft_image_demo.html