Analyzing Images in Fourier Space

Last activity, we played around with images in Fourier space. Given a periodic function \(f(t)\) with fundamental period \(T\) such that \(f(t+T)=f(t)\) that has a finite number of finite discontinuities and finite number of extreme values within the interval \([a,a+T]\) for all \(a\in\mathbb{R}\), then \(f(t)\) can be rewritten as:
\begin{equation}\label{fseries}f(t)=\dfrac{a_0}{2}+\sum_{n=1}^\infty a_n\cos(nt)+\sum_{n=1}^\infty b_n\sin(nt)\end{equation}
[1]. This means that any periodic function, or signal, can be decomposed as a sum of sines and cosines of different frequencies with different contributions. To get the contribution of each frequency of sines and cosines to get \(f(t)\), what we can do is to obtain its Fourier transform, \(\mathscr{F}\{f(t)\}=F(\omega)\) given by:
\begin{equation}\label{1d_ft}F(\omega)=\int_{-\infty}^\infty f(t)\mathrm{e}^{-\mathrm{i}2\pi\omega t}\,\mathrm{d}t\end{equation}
[2]. This can even be extended to work for functions dependent on two functions \(f(t_x,t_y)\), where \(\mathscr{F}\{f(t_x,t_y)\}=F(\omega_x,\omega_y)\) is given by:
\begin{equation}\label{2d_ft}F(\omega_x,\omega_y)=\int_{-\infty}^\infty f(t_x,t_y)\mathrm{e}^{-\mathrm{i}2\pi(\omega_xt_x+\omega_yt_y}\,\mathrm{d}t_x\,\mathrm{d}t_y\end{equation}
[2]. By taking the Fourier transform of \(f(t)\) or \(f(t_x,t_y)\), we are bringing ourselves from real space, to Fourier space, or sometimes called as frequency space, since we look at the frequency distribution of our functions.

However, \eqref{1d_ft} and \eqref{2d_ft} assume \(f(t)\) and \(f(t_x,t_y)\) to be continuous functions of \(t\) and \((t_x,t_y)\). If we want to work with real signals, we have to consider discrete functions \(f[t]\)  and \(f[t_x,t_y]\). Fortunately, they too have a Fourier transform - the Discrete Fourier transform. Assuming \(f[t]\) is sampled over \(M\) samples, and \(f[t_x,t_y]\) is sampled over \((M,N)\) samples, then \eqref{1d_ft} and \eqref{2d_ft} become:
\begin{align}\label{1D_FT} F[\omega_k]&=\dfrac{1}{M}\sum_{m=0}^{M-1}f[t_m]\mathrm{e}^{-\mathrm{i2\pi km/M}} \\ \label{2D_FT} F[\omega_k,\omega_l]&=\dfrac{1}{MN}\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}f[t_m,t_n]\mathrm{e}^{-\mathrm{i}2\pi(km/M+ln/N)}\end{align}
[2].

In this activity, \eqref{1D_FT} and \eqref{2D_FT} was implemented in code, analyzed the Fourier transform of different functions and images, and compared its performance with that of the Fast Fourier Transform (FFT) implemented in computers.

We started by taking the Fourier transform of waves: square, sawtooth, and modulated sine wave. To get their Fourier transform, we have to implement \eqref{1D_FT} to code. Doing so, we have:

def get_1D_fourier_coefficient(y):
    N = len(y)
    n = np.arange(N)
    return np.array([np.sum(y*np.exp(-1j*2*np.pi*k*n/N))/N for k in range(len(y))])

Having these, we can finally get the waves' Fourier transform. Applying the code to the given waves, we get Figure 1.

Figure 1: (a)-(c) The waves undergoing Fourier transform. (d)-(f) Their Fourier transform
as obtained using the given code, and (g)-(i) that obtained using FFT in computer.

In Figure 1, we see the waves that underwent Fourier transform, then that obtained using the code presented above, and that obtained using fft of scipy.fftpack. Notice that the result of both methods are the same, hence the code presented is reliable in obtaining the Fourier transform of any discrete function \(f[t]\). Notice too that since the Fourier transform involves exponential with \(\mathrm{i}=\sqrt{-1}\), then it is expected that \(F(\omega)\in\mathbb{C}\), thus to plot the resulting Fourier transform, we have to get the absolute value of the modulus, \(|F(\omega)|\).

Having done Fourier transform in 1D and since we want to look at Fourier space of 2D images, we go higher into 2D. To begin, we were tasked by taking the Fourier transform of a white circle embedded in a black background using fft2 of scipy.fftpack, then take fftshift, and take the fft2 of the fft2 of the circle. After doing so, we get Figure 2.

Figure 2: (a) The original image, (b) its fft2, (c) fftshift-ed, and (d) the fft2
of the fft2 of the image. (e) and (f) are the real and imaginary parts
of the fft2 of the image, respectively.

Looking at Figure 2b, we see that most of the detail of the Fourier transform is seen at the corners of the image. That's why we fftshift-ed them, so that we bring them back to the center of the image. We have to fftshift the fft2 of the image so that we bring the scale of the image from real space \(x\), to Fourier space \(1/x\) [2]. Looking at the Fourier shifted image, Figure 2c, we see a familiar diffraction pattern, the Airy pattern [3], which is obtained by letting light through a small circular aperture [3]. This suggests that taking the Fourier transform of a white circle in a black background is similar to letting light through that white circle, the aperture.

Taking the Fourier transform of the Fourier transform of the circle, we get back the original circle, as in Figure 2d. This may be because of the symmetry of the Fourier transform of the circle.

Figures 2e and 2f show the real and imaginary parts of the Fourier transform of the circle, as expected since like 1D Fourier transform, the 2D Fourier transform involves an exponential with \(\mathrm{i}=\sqrt{-1}\).

We now took the Fourier transform of different images. The original images and their Fourier transforms are shown in Figure 3.

Figure 3: Other images whose Fourier transform was taken. (a) Letter 'C', (b) corrugated roof,
(c) double slit, (d) square slit, and (e) Gaussian slit. (f)-(j) show their corresponding Fourier transform

We took the Fourier transform of letter 'C', a corrugated roof, double slits, square slit, and Gaussian slit. Notice that Fourier transform of 'C' is just series of fringes. That of corrugated roof is series of dots with distinct spacing, much like the Fourier transform of a square wave but in 2D. The Fourier transform of double slits give us the familiar diffraction pattern for Young's double slit experiment. The same can be said for square slit, where there are diffraction patterns in \(1/x\) and \(1/y\) directions. Lastly, notice that the Fourier transform of a Gaussian slit is also Gaussian, which is expected because Gaussian is \(f(t_x,t_y)=\mathrm{e}^{-(t_x^2+t_y^2)^2/(2w^2)}\) and so \(f(t)\mathrm{e}^{-\mathrm{i}2\pi(\omega_xt_x+\omega_yt_y)}\) is just an exponential, thus can be factored in terms of \(t_x\) and \(t_y\).

Having background on what to expect when taking Fourier transform of images, we implemented \eqref{2D_FT} in code. Doing so we get

def get_2D_fourier_coefficient(y):
    M, N = y.shape[0], y.shape[1]
    m, n = np.arange(M), np.arange(N)
    M_m, N_m = np.meshgrid(m, n)
    return np.reshape(np.array([np.sum(y*np.exp(-1j*2*np.pi*(k*M_m/M+l*N_m/N))) for l in range(N) for k in range(M)]), y.shape)

To test this code, we took the Fourier transform of the previous images, and compared that to the Fourier transform of the same images obtained using fft2. Doing so, we get Figure 4.

Figure 4: (a)-(f) Fourier transform of images in Figures 2a and 3a through 3e
obtained using fft2, (g)-(l) That obtained using the code shown.

Notice that yet again, both methods yield the same result, thus showing the reliability of the code made.

Comparing the performance of each method though, we get Figure 5.

Figure 5: Time it takes to process Fourier transform through self-made code and
built-in code, i.e. FFT.

Notice that it takes significantly lesser time to take the Fourier transform of the image using fft2 rather than the code I made, which makes sense since it is known as the Fast Fourier Transform.

Lastly. we convolved images in Fourier space. In my case, I convolved a square with a circular aperture with variable pixel radius. Doing this we get Figure 6.

Figure 6: (a)-(c) are images in real space. (d)-(f) are images in Fourier space.
We convolved cross with circle in Fourier space, giving (f) in Fourier space
and (c) in real space.

Notice that in convolving images in Fourier space, we take the Fourier transform of the images we have, multiply their results, and take its Fourier transform. Doing so we get Figure 6c. Looking at it, we see that the cross is recognizable, although it has blurred edges, and weird remnants can be seen around it. This might be because of the effect of convolution.

Doing this convolution for circles with different pixel radii, we get Figure 7.

Figure 7: Resulting image when convolving cross with circular aperture with pixel radius \(w\).

From the figure, we see that as radius of circle is increased, the more the cross becomes well defined, and the remnants become much smaller but more frequent. This is reminiscent of imaging systems as it shows how well an imaging system is able to resolve the signals it is sensing/detecting. In this case the larger the aperture, the more we are able to resolve the detail of the image we want to detect.

This activity was rather tedious compared to the others as a lot of images had to be pumped out. At least this one's done!

References:
[1] G. B. Arfken, H. J. Weber, and F. E. Harris, Mathematical methods for physicists: A comprehensive guide (Academic Press, USA, 2013).
[2] R. C. Gonzales, and R. E. Woods, Digital image processing (Prentice Hall, USA, 2002).
[3] n. a., Circular aperture diffraction (Web, Accessed last 09/14/18 on http://hyperphysics.phy-astr.gsu.edu/hbase/phyopt/cirapp2.html)