Fourier Transforms and image filtering

In the previous post we have seen how to compute the fft of a one dimensional signal. In this post we will see how to compute and plot the Fourier spectrum of an image (we can think at images as a bidimensional signals) and how to use the spectrum in order to make an high pass filter in the frequency domain. The following function implements an high pass filter in three steps:

• computes the frequency spectrum
• cuts the low frequencies from the spectrum
• and finally, the filtered time domain signal is calculated from the modified frequency spectrum

from pylab import *
from numpy import *

def freqHighPass(img, r):
 """
  Applies an High Pass filter to img in the frequency domain.
  In:  img, image to filter
       r, percentage of spectrum to keep (in [0 1])
  Out: ifiltered, filtered version of the image
       P, filtered spectrum
 """
 I = fftshift(fft2(img)) # entering to frequency domain
 # fftshift moves zero-frequency component 
 # to the center of the array
 P = zeros(I.shape,dtype=complex)

 r = int((r*min(img.shape))/2);
 c1 = I.shape[0]/2 # spectrum center
 c2 = I.shape[1]/2

 for i in range(c1-r,c1+r):  # frequency cutting
  for j in range(c2-r,c2+r): # around the center
   P[i,j] = I[i,j]

 ifiltered = real(ifft2(ifftshift(P)))
 return ifiltered,P # back to the spatial domain

This is the first application of the function freqHighPass():

i = imread('diane.jpg') # load an image
i = mean(i,2) # to get a 2-D array

ifiltered,P = freqHighPass(i,0.4) # 40% of spectrum preserved

subplot(2,2,1)
title('Original - Fourier Specturm')
imshow(log(abs(fftshift(fft2(i))))**2)
subplot(2,2,2)
title('Filtered - Fourier Specturm')
imshow(log(abs(P))**2)
subplot(2,2,3)
title('Original image')
imshow(flipud(i))
gray()
subplot(2,2,4)
title('Filtered image')
imshow(flipud(ifiltered))
show()

The result of the code above is the following plot. We can see the spectrum of the original image, the filtered spectrum (the matrix P returned by freqHighPass), the original image and the filtered image.

And this is another application that preserve more spectrum than the first:

i = imread('castle.png')

ifiltered,P = freqHighPass(i,0.7) # 70% of spectrum preserved
figure()
subplot(1,2,1)
title('Original image')
imshow(i)
gray()
subplot(1,2,2)
title('Filtered image')
imshow(ifiltered)
show()

In this case the original image is affected by noise. As we can see the filter is able to reduce noise.