1. Introduction
This document introduces the Fourier transform of an image, then the discrete Fourier transform (DFT) of a sampled image. The calculation of the DFT of an image with Python is explained. We will see how to represent the spectrum of the image and how to perform filtering in the frequency space, by multiplying the DFT by a filtering function.
2. Fourier transform and discrete Fourier transform
We consider a monochrome image (gray levels) represented by a function of two real variables, with complex values, denoted u (x, y).
The Fourier transform of this image is the function with two real variables and with complex values defined by: S (fx, fy) = ∫-∞∞∫-∞∞u (x, y) exp (-i2π (fxx + fyy )) dxdy
We now consider a sampling of the image on a rectangular domain of dimensions (Lx, Ly) centered in (x, y) = (0,0) and comprising (Nx, Ny) points. The matrix U is defined by: xk = -Lx2 + kLxNx0≤k≤Nx-1ym = -Ly2 + mLyNy0≤m≤Ny-1Um, k = u (xk, ym)
The matrix U therefore has Nx columns and Ny rows.
Consider an approximation of the Fourier transform obtained by the rectangle method on the domain considered: S (fx, fy) ≃LxLyNxNyexpiπ (fxLx + fyLy) ∑k = 0Nx-1∑m = 0Ny-1Um, kexp-i2πfxkLxNxexp-i2
In practice, this approximation is calculated for the following spatial frequencies: fx, n = nLx0≤n≤Nx-1fy, l = lLy0≤l≤Ny-1
The discrete Fourier transform (DFT) associates with the matrix U a matrix V of the same dimensions, defined by: Vn, l = 1NxNy∑k = 0Nx-1∑m = 0Ny-1Um, kexp-i2πknNxexp-i2πmlNy
Note: the DFTs calculated by the software FFT functions sometimes omit the inverse of NxNy as a factor.
The approximation of the integral is: S (fx, n, fy, l) ≃LxLyexpiπ (n + l) Vn, l
The sampling frequency is constituted by the pair: fe = NxLx, NyLy
We easily verify the following property of the DFT: Vn + Nx, l = Vn, lVn, l + Ny = Vn, l
This property corresponds to the periodicity of the spectrum of the sampled image. The period of the spectrum is equal to the sampling frequency. In general, we try to calculate an approximation of the Fourier transform that does not show this effect of sampling. To do this, limit the highest frequency to half the sampling frequency. For simplicity, we assume that Nx and Ny are even. The DFT gives an approximation of the Fourier transform for the following frequencies: fx, n = nLx0≤n≤Nx2fy, l = lLy0≤l≤Ny2
In addition, the frequency of the DFT allows access to the opposites of these frequencies. Indeed: V-n, l = V-n + Nx, lVn, -l = Vn, -l + Ny
The Nx indices n correspond to the following frequencies: 0.1Lx, 2Lx,…, Nx2Lx, -Nx2-1Lx,…, -2Lx, -1Lx
The same relationship is valid for the frequencies of the y axis.
For certain calculations (filtering, diffraction), it is useful to place the zero frequency in the middle of the matrix, so as to obtain increasing frequencies as a function of the indices. Denote by VC this centered matrix. In this matrix, the point of zero frequency must be placed at the following indices: (Ny2-1, Nx2-1)
For the x axis, there are therefore Nx2 + 1 indices corresponding to positive (or zero) frequencies and Nx2-1 indices corresponding to strictly negative frequencies.
To build the VC matrix, the following index transformation must be used: nc = Nx2-1 + nsi0≤n≤Nx2nc = -Nx2-1 + nsiNx2 + 1≤n≤Nx-1
3. Fourier transform of an image with Python
3.a. Principe
The imageio module allows you to read image files. The following image has 200×100 pixels.
The image is obtained as a three-dimensional numpy array with the imageio.imread function. We remove the red layer and we convert the array to floats.
import imageio
import numpy
imA = imageio.imread("imageA.png")
U=numpy.array(imA[:,:,0],dtype=numpy.float64)
import numpy
import math
from pylab import *
The following function converts a matrix into an RGB image matrix, performing normalization and gamma correction. Gamma correction makes it possible to simulate the response of the eye for diffraction figures.
def matriceImage(matrice,gamma,rgb):
s = matrice.shape
a=1.0/gamma;
norm=matrice.max()
m = numpy.power(matrice/norm,a)
im = numpy.zeros((s[0],s[1],3),dtype=float64)
im[:,:,0] = rgb[0]*m
im[:,:,1] = rgb[1]*m
im[:,:,2] = rgb[2]*m
return im
The TFD modulus decreases very rapidly with frequency. To better visualize the spectrum, we use a logarithmic scale. The following function generates an RGB image from a matrix, by applying a logarithm function:
def matriceImageLog(matrice,rgb):
s = matrice.shape
m = numpy.log10(1+matrice)
min = m.min()
max = m.max()
m = (m-min)/(max-min)
im = numpy.zeros((s[0],s[1],3),dtype=float64)
im[:,:,0] = rgb[0]*m
im[:,:,1] = rgb[1]*m
im[:,:,2] = rgb[2]*m
return im
The next function displays the spectrum. The axes are graduated with the spatial frequencies. The dimensions of the image are given as an argument:
def plotSpectre(image,Lx,Ly):
(Ny,Nx,p) = image.shape
fxm = Nx*1.0/(2*Lx)
fym = Ny*1.0/(2*Ly)
imshow(image,extent=[-fxm,fxm,-fym,fym])
xlabel("fx")
ylabel("fy")
from TfdImage import *
from numpy.fft import fft2, ifft2, fftshift, ifftshift
The discrete Fourier transform is calculated with the function numpy.fft.fft2. The zero frequency centering is obtained with the function numpy.fft.fftshift. We calculate the power, that is to say the squared modulus of the TFD and we plot the corresponding image in red.
V=fft2(U)
VC = fftshift(V)
P = numpy.power(numpy.absolute(VC),2)
img = matriceImage(P,2.0,[1.0,0.0,0.0])
plotSpectre(img,200.0,100.0)
We chose the dimensions of the image equal to the number of pixels, which gives a maximum frequency of 0.5 on the spectrum. Frequency 1 is the sample rate, which corresponds to 1 per pixel. We thus have a square spectrum (the pixels are square), although the DFT is a rectangular matrix (of the same dimensions as the starting image). In the present case, the image has a length Lx twice as great as its height Ly, so the resolution of the spectrum is twice as great on the X axis.
Voici l’affichage en échelle logarithmique :
img = matriceImageLog(P,[1.0,0.0,0.0])
plotSpectre(img,200.0,100.0)
To perform filtering on the Fourier transform of the image, we must define a transfer function. Here is for example a low pass filtering:
def transfert(fx,fy,p):
fc=p[0]
return 1.0/numpy.power(1+(fx*fx+fy*fy)/(fc*fc),2)
Frequencies are relative to the sampling frequency. For example, a frequency fx = 1/8 means a period of 8 pixels. Filtering in the frequency domain consists in multiplying the DFT by a matrix H which is obtained (in first approach) by sampling of the transfer function.
The following function performs the multiplication of a matrix by the sampling of a transfer function. The zero frequency of the matrix is centered.
def matriceFiltre(matrice,transfert,p):
s = matrice.shape
Nx = s[1]
Ny = s[0]
nx = Nx/2
ny = Ny/2
Mat = zeros((Ny,Nx),dtype=numpy.complex)
for n in range(Nx):
for l in range(Ny):
fx = float(n-nx-1)/Nx
fy = float(l-ny-1)/Ny
Mat[l,n] = matrice[l,n]*transfert(fx,fy,p)
return Mat
Let’s see the use of this function for the transfer function defined above. Initially, we must create a matrix H of the same size as the image and filled with 1.
H = numpy.ones(VC.shape,dtype=numpy.complex)
H = matriceFiltre(H,transfert,[0.1])
VCF = VC*H
PF = numpy.power(numpy.absolute(VCF),2)
imPF = matriceImageLog(PF,[1.0,0,0])
plotSpectre(imPF,200,100)
To obtain the filtered image, it is first necessary to restore the initial order of the frequencies with the function ifftshift then apply the inverse discrete Fourier transform with the function ifft2:
VF = ifftshift(VCF)
UF=ifft2(VF)
imageF = matriceImage(UF,1.0,[1.0,1.0,1.0])
figure()
imshow(imageF)
As a second example of a filter, we define a high-pass filter:
def transfert(fx,fy,p):
fc=p[0]
return 1.0-1.0/numpy.power(1+(fx*fx+fy*fy)/(fc*fc),2)
H = numpy.ones(VC.shape,dtype=numpy.complex)
H = matriceFiltre(H,transfert,[0.4])
VCF = VC*H
VF = ifftshift(VCF)
UF=ifft2(VF)
imageF = matriceImage(UF,1.0,[1.0,1.0,1.0])
figure()
imshow(imageF)
3.b. Application: filtering a photograph
imA = imageio.imread("imageC.png")
U=numpy.array(imA[:,:,0],dtype=numpy.float64)
s = U.shape
figure(figsize=(12,6))
imshow(U,cmap='gray')
V=fft2(U)
VC = fftshift(V)
P = numpy.power(numpy.absolute(VC),2)
img = matriceImageLog(P,[1.0,1.0,1.0])
figure(figsize=(6,6))
plotSpectre(img,s[1],s[0])
This photograph has a periodicity in the X direction because of the repetition of the windows. There is also the repetition of the bars of the balustrades and the junctions on the zinc roof. The periodicity of the windows is 140 pixels, ie a frequency of 1/140 = 0.00714, which can be seen on the spectrum in the form of very close vertical lines. The bars of the balustrades are spaced 6 pixels apart, i.e. a frequency of 1/6 = 0.16. The pattern of the frieze located under the first balcony has a period of 11 pixels, i.e. a frequency of 1/11 = 0.091.
As an example of a filter, consider a Gaussian notch filter that acts on horizontal frequencies:
def transfert(fx,fy,p):
f0 = p[0]
sigma = p[1]
ux = numpy.absolute(fx)-f0
return 1.0-numpy.exp(-(ux*ux)/(sigma*sigma))
This filter is used to remove the lines located at a horizontal frequency multiple of 0.091, which correspond to a periodicity of 11 pixels.
H = numpy.ones(VC.shape,dtype=numpy.complex)
H = matriceFiltre(H,transfert,[0.091,0.01])
H = matriceFiltre(H,transfert,[0.091*2,0.01])
H = matriceFiltre(H,transfert,[0.091*3,0.01])
VCF = VC*H
PF = numpy.power(numpy.absolute(VCF),2)
imPF = matriceImageLog(PF,[1.0,1.0,1.0])
figure(figsize=(6,6))
plotSpectre(imPF,s[1],s[0])
VF = ifftshift(VCF)
UF=ifft2(VF)
imageF = matriceImage(UF,1.0,[1.0,1.0,1.0])
figure(figsize=(12,6))
imshow(imageF)
The filtering made the frieze located under the first balcony disappear. There is also an overall loss of sharpness.
It is interesting to calculate the impulse response of this filter. To do this, you have to create a black image with a single white pixel.
impuls = numpy.zeros(s)
y0 = s[0]/2
x0 = s[1]/2
impuls[y0,x0] = 1.0
V=fft2(impuls)
VC = fftshift(V)
VCF = VC*H
P = numpy.power(numpy.absolute(VCF),2)
img = matriceImage(P,1.0,[1.0,1.0,1.0])
figure(figsize=(6,6))
plotSpectre(img,s[1],s[0])
This filter only acts on frequencies fx, so its impulse response has only one line. We can therefore plot it with plot:
reponse = ifft2(ifftshift(VCF))
figure()
h = reponse[y0,x0-100:x0+100]
plot(h,'.')
From the impulse response, we can get the frequency response in the X direction, as we do for a univariate signal:
import scipy.signal
w,g=scipy.signal.freqz(h)
figure()
plot(w/(2*numpy.pi),numpy.absolute(g))
xlabel("f/fe")
ylabel("G")
grid()
When we have the impulse response, we can do the filtering by convolution:
from scipy.signal import convolve2d
h1 = numpy.zeros((1,h.size))
h1[0,:] = h
imgF = convolve2d(U,h1,mode='same')
figure(figsize=(12,6))
imshow(imgF,cmap='gray')
The result is similar to that obtained by filtering in frequency space. However, there is a more significant edge effect with convolution filtering, in the form of vertical bands to the left and to the right of the image.
Spectral analysis of the image is necessary to design very selective filters like the previous one. For performing the filtering, there is a choice between filtering in the frequency domain (also called filtering by Fourier transform) or filtering in the spatial domain (filtering by convolution). To compare the effectiveness of these two methods, consider a square image with n rows and n columns. The fast Fourier transform for a line takes a time an ln (n). For n lines, that makes an2 ln (n). After application of the TFD on each row, the TFD must be applied on each column. The total time for the FFT of the image is therefore 2an2 ln (n). Filtering in the frequency domain consists in multiplying each element of the DFT, which takes a time bn2. By taking into account the inverse FFT, we finally obtain for the filtering by Fourier transform a time: τ1 = n2 (4aln (n) + b)
Let Q be the size of the impulse response, ie the convolution mask. For a filter having a certain response with respect to the sampling frequency, the mask size is fixed. The number of operations in a convolutional filtering is equal to the number of pixels in the image multiplied by the size of the mask: τ2 = cQn2
The ratio of the two times is therefore: τ1τ2 = 4aln (n) + bcQ (1)
Within the limit of large images (n tending to infinity), convolutional filtering is much more efficient. However, for a fixed image size, the Fourier transform filtering can be faster when Q is large, that is to say when the impulse response has a size of the same order of magnitude as the image.
In addition to its greater efficiency, convolution filtering also has another advantage: the convolution mask can vary from one area of the image to another. In other words, we can define a filter with a local impulse response. Conversely, filtering in the frequency domain can only act on the image taken as a whole.
To continue with the previous example, we can apply convolutional filtering only on part of the image:
partie = U[550:580,0:s[1]]
partieF = convolve2d(partie,h1,mode='same')
img = U.copy()
img[550:580,0:s[1]] = partieF
figure(figsize=(12,6))
imshow(img,cmap='gray')
