A5: Fourier Transform Model of Image Formation

The Fourier Transform (FT) has many uses in imaging. We're to demonstrate its basic uses in this activity.

Discrete FFT

Let us first look at the equation for the FT F(f_x,f_y) of an image f(x,y):

This image, initially in the x,y domain, is basically expressed in the inverse or 1/x,1/y domain. Hence, FT reports the frequency at which a certain spatial pattern appears.

As in many integration algorithms done on the machine, discretization is needed. On Scilab, the efficient Fast Fourier Transform (FFT) algorithm of Cooley and Turkey is used in the fft2 function (for two-dimensional images). However, this function outputs an FFT with the diagonals in reverse order. The fftshift function corrects this. fft2 also outputs a matrix of complex numbers. To get only the intensity for plotting, we take only the modulus with the use of the abs command. Hence, to get the FFT of an image on Scilab:

img = imread('ImageDirectory\ImageFilename.jpg'); //import image
grayimg = im2gray(img); //convert to grayscale
FFTimg = fft2(grayimg); //get FFT
imshow(abs(FFTimg),[]); //view FFT

We see in Figure 1 and 2 that the resulting FFT has peaks at the corners because of the wrongly positioned diagonals on the matrix. After correction using fftshift, we find the expected peaks in the middle. For the circle we find a circle with rings that decrease in intensity as you go farther from the peak, which is expected. For the letter "A", we find what appears to be a six-pointed star.

The inverse FT simply returns an inverted version of the image. We arrive at the inverse FFT by applying fft2 on the FFT of the image (FFTing raw image twice). From Figure 1, we do not see any change after FFT and inverse FFT because it is symmetric. But from Figure 2, we see that the process inverts the image.

Figure 1. Left to right: input, FFT, shifted FFT, and inversed FFT images for a circle.


Figure 2. Left to right: input, FFT, shifted FFT, and inverse FFT images for the letter "A".

Convolution and the lens

In convolving two two-dimensional arrays, we use the following equation:

In imaging systems, f = object being imaged, g = imaging system transfer function, and h = image produced by imaging system.

We note that if we let H, F, and G be the FT of h, f, and g, respectively, the convolution theorem says that H = FG. Hence, if we just get the FFT of the object and the imaging system, we can get the resulting image's FFT by simple multiplication, which may then be inverse FFTed to get the resulting image function. We demonstrate this on Scilab using an image containing the text "VIP" as the object and the a circle as the imaging system or lens. We implement this on Scilab using the following code:

img = imread('ObjectDirectory\ObjectFilename.jpg');
lens = imread('LensDirectory\LensFilename.jpg');
imggray = im2gray(img); //convert to grayscale
lensgray = im2gray(lens); //convert to grayscale
FFTimg = fft2(imggray); //get FFT of object
FFTlens = fftshift(lensgray); //lens is Fourier plane and need not be FFT'ed
convimglens = FFTimg.*(FFTlens); //convolution (element per element multiplication)
invconvimglens = fft2(convimglens); //inverse FFT
modimg = abs(invconvimglens);
imshow(modimg, []);

From Figure 4 below, we see that resolution increases with increasing aperture size. This is expected since a larger aperture allows more light from the object through.

Figure 3. "VIP" text (left) as object and circle (right) as lens aperture.

Figure 4. Resulting images for aperture diameters: 20px, 50px, and 126px (left to right).

Correlation and template matching

In correlating two two-dimensional arrays, we use the following equation:

And as in convolution, if we let H, F, and G be the FT of h, f, and g, respectively, the correlation theorem says that P = F^*G, where the asterisk means that we take the complex conjugate of F. Hence, if we just get the conjugate of the FFT of the object and the FFT of the imaging system, we can get the resulting image's FFT by simple multiplication, which may then be inverse FFTed to get the resulting image function. This multiplication effectively gives information on locations where the pattern in g may be found in f. We demonstrate this on Scilab using two images that we know share a pattern. We implement this on Scilab using the following code:

//note that we use bmp files for sharper text edges
text = imread('fDirectory\fFilename.bmp'); //import f
A = imread('gDirectory\gFilename.bmp'); //import g/pattern
textgray = im2gray(text); //convert to grayscale
Agray = im2gray(A); //convert to grayscale
FFTtext = fft2(textgray); //get FFT
FFTA = fft2(Agray); //get FFT
corrtextA = FFTA.*(conj(FFTtext)); //FFT of correlation
invcorrtextA = fft2(corrtextA); //inverse FFT to get correlation
modimg = abs(invcorrtextA); //get modulus since complex
imshow(fftshift(modimg), []); //visualize

From Figure 5 (rightmost) we find that, as expected, the correlation peaks at parts where the letter/pattern "A" is found in the long text.

Figure 5. Left to right: text that contains A's with font style Arial and font size 16px, pattern/letter "A" with font style Arial and font size 16px, and the resulting correlation image.

Convolution and edge detection

Convolution may be done easily using Scilab's imcorrcoef function. Convolution may be used in edge detection:

vip = imread('VIPDirectory\VIPFilename.bmp'); //import image
vipgray = im2gray(vip); //convert to grayscale
//set patterns
//pat = [-1 -1 -1; 2 2 2; -1 -1 -1]; //horizontal
//pat = [-1 2 -1; -1 2 -1; -1 2 -1]; //vertical
//pat = [-1 -1 2; -1 2 -1; 2 -1 -1]; //increasing diagonal
//pat = [2 -1 -1; -1 2 -1; -1 -1 2]; //decreasing diagonal
pat = [-1 -1 -1; -1 2 -1; -1 -1 -1]; //spot
conv = imcorrcoef(vipgray,pat);
imshow(conv);

It is clear in Figure 6 that the edges detected depended on the pattern used. Vertical edges were missing on the result for the horizontal edge pattern, horizontal edges were missing on the result for the vertical edge pattern, and so on. The spot pattern was most effective in tracing the edge, as expected because all edges are made up of points.

Figure 6. Edges detected for (clockwise): horizontal, vertical, increasing diagonal, decreasing diagonal, and spot patterns.

I give myself a grade of 9 because all required images for the demonstration of the uses of FFT were produced.

I would like to thank Mr. Luis Buño III for his help in troubleshooting some errors in getting the FFTs.