Welcome to STARKFX.com page on Using the Discrete Fourier Transform

On this page, I will show my matlab code for taking advantage of the DFT (discrete fourier transform) to process images, allowing me to choose a given set of spatial frequencies to allow in reconstructing an image.

Image Basics

Looking at an Image as a 2D Array of Pixels, we can use a coordinate system to display the Brightness value for a given pixel

I(x,y)
Our eyes (and brain) will interpret this array of brightness values as an image. In order to simplify our treatment of the Discrete Fourier Transform in image processing, we will only use images in gray scale - no color.

So, let's read in an image file:

A=imread(frathouse.jpg);

If we display this image, we see it is in color:

figure; imshow(A);

So, we will change it to grayscale using a matlab provided function that uses the RGB values in the image array to determine a grayscale value to replace it. That way we have a single 2D array rather than 3 such 2D arrays representing the red (R), green (G) and blue (B) component values that make up most jpeg images.

A=rgb2gray(A);

We will then do a little conditioning to put the values for the pixels all in the range from 0 to 1:

A=double(A);
A=A-min(A(:));
A=A/max(A(:));

figure; imshow(A,[0 1]); colormap gray; title('Original image');

Now we can see it is a grayscale image:

Take the 2D discrete fourier transform:

fftA=fft2(A);

Split out magnitude and phase:

magfftA=abs(fftA);
phasefftA=angle(fftA);

This splitting process works because the result of the fft2 function is a complex array. The magnitude given by the absolute value of this array is the magnitude of the spatial frequency component represented by that location in the array. The angle gives the phase for that spatial frequency component. (Phase actually turns out to be the more important value for reconstructing an image since it tells the phases relative to each other for all of the spatial frequency components so that when they add together, we get the particular light and dark pattern that makes up our image.

Now, we are going to apply a little trick in order to visualize our FFT results a bit more intuitively. The spatial frequencies represented by the FFT in a 2D array span from f=[-N/2,-(N-1)/2, -(N-2)/2....0.....(N-2)/2, (N-1)/2]*Fs/N where Fs is the sampling frequency (1 sample per pixel) and N is the number of pixels in that dimension. So as in our example with 964 pixels in the x direction and 686 pixels in the y direction we get frequencies from negative 1/2 cycle per pixel to a positive 1/2 cycle per pixel. Negative and positive frequencies really represent the same frequencies so we don't need to keep both and you will see that the fft results are always symmetric with respect to the negative and positvie frequencies. Resolution in frequency space, or the step size in spatial frequency represented as we move from one point in the array to an adjacent point in the array is 1/N or, in our case, just looking at the x dimension, 1/964 cycle per pixel. These frequencies are nicely represented in an array with the zero frequency component (average level of the whole image) in the center of the array with increasing frequencies stepping at a rate of 1/964th cycle per pixel as we step in the x direction and 1/686th cycle per pixel as we step in the y direction away from the center of the image.

But, that is not exacty what we get from the output of the fft2 function. Instead, we get an array with the zero frequency components in the corners of the array and the highest frequency at the center (1/2 cycle per pixel in either direction). So, we use a function that simply, directly swaps our quadrants such that the frequency of one half cycle per pixel is at the edges of the array, while the zero component is at the center. This is really just for our benefit in displaying the magnitude and phase information. The function that accomplishes this switch is:

fftshift();

So, let's look at the phase information in an array form:

figure; imshow(fftshift(phasefftA),[-pi pi]); colormap gray; title('FFT Phase');

Notice that the values of the phase run from -pi to pi.

We are also going to use a little trick to enhance the contrast in displaying the magnitude at each spatial frequency component. We will take the log of the value so that the huge component magnitude at very low spatial frequencies does not completely drown out the high spatial frequency components in our display of the magnitudes at the various spatial frequencies.

sfftA=fftshift(fftA);
osa=log(abs(sfftA));
osa=osa-min(osa(:));
osa=osa/max(osa(:));
figure; imshow(osa,[0 1]); colormap gray; title('FFT magnitude');

Now, if we want to see what the image looks like if we use only the high spatial frequency components or the low spatial frequency components, or even some region between a high and low cutoff, we simply create a mask to remove the unwanted spatial frequencies. That mask can be an array of ones and zeros to multiply our quadrant shifted fft array, (sfft(x,y)):

[rows, cols]=size(sfftA); %get the number of rows and columns.
centrow=round(0.5+rows/2);
centcol=round(0.5+cols/2);
dist=sqrt((centrow-1)^2 + (centcol-1)^2);
%Create the masks for filters
hpfilt=zeros(rows, cols); % Hi pass filter
lpfilt=zeros(rows, cols); % Low pass filter
rpfilt=ones(rows, cols); % Region pass filter

Now, we set the values in the filter arrays and simply multiply the shifted fft array by the mask, undo the fft shift before we reconstruct the image.

In filling the masks, we can choose to use straight up ones or zeros, to turn on or off completely any given spatial frequency component. Or, we can use any number between 0 and 1 to choose how much of a given spatial frequency component to use. In the case of using a chosen amount of a frequency, I like to use a gaussian filter with a maximum of one and approaching a minimum of zero.

My Gaussian mask making function:

function [ SGmask ] = SGaussMask( SzRows, SzCols, Width, Invert)
%
% for inverting, either nothing or N is false, anything else is true;
if nargin < 4
Invert=false;
else
if strcmpi(Invert,'N')
Invert = false;
else
Invert = true;
end
end
SzRows=round(abs(SzRows));
SzCols=round(abs(SzCols));
if nargin < 3
Width=0.3*min(SzRows, SzCols);
end
Width=double(Width);

SGmask=double(zeros(SzRows, SzCols));
centrow=double(SzRows)/2;
centcol=double(SzCols)/2;
for i=1:SzRows
for j=1:SzCols
SGmask(i,j)=exp(-((i-centrow)^2+(j-centcol)^2)/Width^2);
end
end
if Invert
SGmask=double(ones(SzRows, SzCols))-SGmask;
end
end

I can use that function to create my gaussian type filters.

w=cop*min(rows,cols); %where cop is the cutoff percent.
hpfilt=SGaussMask(rows,cols,w,'Y');
lpfilt=SGaussMask(rows,cols,w,'N');
wl=lpcf*min(rows,cols);%where lpcf is the low cut-on
wh=hpcf*min(rows,cols);%where hpcf is the high cut-off
rpfilt=SGaussMask(rows,cols,wh,'Y').*SGaussMask(rows,cols,wl,'N');
rmx=max(rpfilt(:));
rpfilt=rpfilt/rmx;

This creates filters that look like the following:

High Pass Filter from Gaussian

Low Pass Filter from Gaussian

Region Pass Filter from Gaussian

Application of these filters is shown this way:

Applying the gaussian filters:

hpfftA=ifftshift(sfftA.*hpfilt);
lpfftA=ifftshift(sfftA.*lpfilt);
rpfftA=ifftshift(sfftA.*rpfilt);

With these filtered fourier transforms we can recontruct the image:

hpim=abs(ifft2(hpfftA));
lpim=abs(ifft2(lpfftA));
rpim=abs(ifft2(rpfftA));

And, now, show them off.

figure; imshow(hpim,[0 1]); colormap gray; title('hp filtered image Gaussian Used');

figure; imshow(lpim,[0 1]); colormap gray; title('lp filtered image Gaussian Used');

figure; imshow(rpim,[0 1]); colormap gray; title('region filtered image Gaussian Used');

Results are a little different if use the simple on-or-off filters.

High Pass Filter on/off

Low Pass Filter on/off

Region Pass Filter on/off

Application of these filters is shown this way:

Applying the on/off filters:

hpfftA=ifftshift(sfftA.*hpfilt);
lpfftA=ifftshift(sfftA.*lpfilt);
rpfftA=ifftshift(sfftA.*rpfilt);

With these filtered fourier transforms we can recontruct the image:

hpim=abs(ifft2(hpfftA));
lpim=abs(ifft2(lpfftA));
rpim=abs(ifft2(rpfftA));

And, now, show them off.

figure; imshow(hpim,[0 1]); colormap gray; title('hp filtered image');

figure; imshow(lpim,[0 1]); colormap gray; title('lp filtered image');

figure; imshow(rpim,[0 1]); colormap gray; title('region filtered image');

That was my quick show of how to use the DFT for image processing, using matlab in particular. I'll end this article by showing the results of a gaussian filter applied to the magnitude of an fft before recontructing the image using an image from my windows OS background options.

Now, I want to do something totally nonsensical. Instead of doing any filtering, I want to use the phase from my xmas girls picture:

Changed to gray scale:

with the magnitude information from my garden path picture:

Changed to gray scale:

The result looks like a strange filtered version of the one we used the phase information for!

A=imread('xmasgirls.jpg');
B=imread('gardenpath.jpg');
A=rgb2gray(A);
B=rgb2gray(B);
A=double(A);
B=double(B);
A=A/255;
B=B/255;
fftA=fft2(A);
fftB=fft2(B);
magfftA=abs(fftA);
phasefftA=angle(fftA);
magfftB=abs(fftB);
phasefftB=angle(fftB);
imafterfftAB=abs(ifft2(magfftB.*exp(i*phasefftA)));
figure: imshow(imafterfftAB)