13.4. The Fourier Transform#

The 2D discrete Fourier Transform decomposes a 2D array into a sum of complex exponentials:

\[ f(x,y) = \frac{1}{mn} \sum_{u=0}^{m-1} \sum_{v=0}^{n-1} F(u,v) \exp \left[2\pi i \left(\frac{xu}{m} + \frac{yv}{n}\right) \right] \]

where the Fourier coefficients are defined by

\[ F(u,v) = \sum_{x=0}^{m-1} \sum_{y=0}^{n-1} f(x,y) \exp \left[-2\pi i \left(\frac{xu}{m} + \frac{yv}{n}\right) \right] \]

They are useful in image processing for performing fast filtering, compression, detection of periodic features, etc.

First let us plot the Fourier coefficients of some simple images:

# Code from previous section
using PyPlot
using Statistics
A = imread("sample_photo.png")
B = mean(A, dims=3)[:,:,1]
function imshow_scale(A)
    # Like imshow(A) but scales the values to [0,1] and supports grayscale
    
    A .-= minimum(A)            # Scale and shift to [0,1]
    A ./= maximum(A)
    if ndims(A) < 3
        A = reshape(A, size(A,1), size(A,2), 1)
    end
    if size(A,3) == 1
        A = repeat(A, 1, 1, 3)  # Set R=G=B for grayscale
    end
    imshow(A)
    end;
# Uncomment below if the package is not already installed
#using Pkg; Pkg.add("FFTW") 
using FFTW

function imagefft_demo(A)
    AF = fftshift(fft(A))
    subplot(1,2,1); imshow_scale(A);
    subplot(1,2,2); imshow_scale(log.(1 .+ abs.(AF)));
    return
end
imagefft_demo (generic function with 1 method)
G = [zeros(256,128) ones(256,128)]
imagefft_demo(G)
../../_images/aeccf504daefceea58f8ecb6139c0568b881162ffe8029d7fd62969f0162cef3.png
G = zeros(256, 256)
G[78:178, 78:178] .= 1.0
imagefft_demo(G)
../../_images/52a1f121d006a50da47b0263832d04212f1c083debe758ec3dd69fb7e39ae1f7.png
G = Float32[ (i+j<329) && (i+j>182) && (i-j>-67) && (i-j<73) for i = 1:256, j = 1:256 ]
imagefft_demo(G)
../../_images/10647fbc9ac31a0c7b843838b4458f7316810c83a5bfc2673f86b727304d1813.png
G = Float32[ sin(2pi*(10i + 20j)/256) for i = 1:256, j = 1:256 ]
imagefft_demo(G)
../../_images/a2f7275115ba9bda9317172278ad1cd9bcd1b23f0b61da7677c405eaa689749a.png

Many of these patterns can be understood from the underlying Fourier expansion. However, for a general image the pattern of the Fourier spectrum is less clear. We would expect the coefficients from a relatively smooth image to decay away from the center though:

imagefft_demo(B)
../../_images/7c2da13ee20aad5da71ffad22903d3370b5a736afd672be159309d5cea9fafdc.png

13.4.1. Removing periodic noise#

One application of image processing using the Fourier transform is to remove periodic noise. Below we demonstrate this using a made-up example with a given frequency and direction of the noise, but it can be made more general.

Bpernoise = copy(B)
Bpernoise = B + 0.5*Float32[sin(2π*10j / size(B,2)) for i = 1:size(B,1), j = 1:size(B,2) ]
imagefft_demo(Bpernoise)
../../_images/b073c54859615ca620ac9eea56cf81d7799aae29ec004131ef0081c66b521af2.png

Now compute the Fourier transform, and set the (known) noise frequencies to zero:

# Filter
BF = fftshift(fft(Bpernoise))
mid = size(B)  2 .+ 1
BF[mid[1], mid[2] + 10] = 0
BF[mid[1], mid[2] - 10] = 0
Bfiltered = real.(ifft(ifftshift(BF)))

imagefft_demo(Bfiltered)
../../_images/6b11804b8a7857603a505c5e9f00403b5c5212796e1ea38ccb22d5ec9317c745.png