Performing Fourier Analysis using Python for image reconstruction

Question

I'm trying to perform a Fourier analysis on some shapes I produced using Python. I managed to obtain a 2D Fourier transform on the images as well as applying a Gaussian filter, however the inverse of the image with the Gaussian filter is being shifted and I don't know how to resolve it. I tried to take the real values only in the hopes that it would solve it but nothing changed. I'm demonstrating my code below, does someone have any ideas on what I could do?

#importing necessary libraries 
from PIL import Image, ImageDraw
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter


#function to apply Gaussian filter to the Fourier transform
def apply_gaussian_filter(fourier_transform, sigma):
    return gaussian_filter(fourier_transform, sigma)

##code to produce an image of a box and perform Fourier analysis on it

#defining the size bounds of the box
w, h = 65, 65
box = [(50, 29), (w-50, h-30)]

#creating the image with a white background
im = Image.new("I", (w,h))
img = ImageDraw.Draw(im)

#producing the box shape
img.rectangle(box, fill = "white", outline = "white")

#convering image to numpy library
im_array = np.array(im)

#taking the 2D Fourier transform
fourier_box = np.fft.fft2(im_array)

#shifting the zero frequency component to the center
fourier_box_shifted = np.fft.fftshift(fourier_box)

#applying Gaussian filter to the Fourier transform of the box
sigma_box = 5
filtered_fourier_box = apply_gaussian_filter(np.abs(fourier_box_shifted), sigma_box)

##code to produce an image of a circle and perform Fourier analysis on it

#defining the size bounds of the shape
w_1, h_1 = 65, 65
circle = [(10, 10), (w_1-10, h_1-10)]

#creating the image with a white background
im_1 = Image.new("I", (w_1,h_1))
img_1 = ImageDraw.Draw(im_1)

#producing the circle shape
img_1.ellipse(circle, fill = "white", outline = "white")

#convering image to numpy library
im_1_array = np.array(im_1)

#taking the 2D Fourier transform
fourier_circle = np.fft.fft2(im_1_array)

#shifting the zero frequency component to the center
fourier_circle_shifted = np.fft.fftshift(fourier_circle)

#applying Gaussian filter to the Fourier transform of the circle
sigma_circle = 5
filtered_fourier_circle = apply_gaussian_filter(np.abs(fourier_circle_shifted), sigma_circle)

#applying inverse Fourier transform to the filtered Fourier data
inverse_filtered_box = np.fft.ifft2(np.fft.ifftshift(filtered_fourier_box)).real
inverse_filtered_circle = np.fft.ifft2(np.fft.ifftshift(filtered_fourier_circle)).real

#adjusting size of the graphs 
plt.figure(figsize=(13, 8))

#displaying the Fourier transform of the box 
plt.subplot(231)  
plt.imshow(np.log(np.abs(fourier_box_shifted) + 1))
plt.title('2D Fourier Transform Box')

#displaying the Fourier transform of the box with Gaussian filter
plt.subplot(232) 
plt.imshow(np.log(filtered_fourier_box + 1))
plt.title('Gaussian Filtered 2D Fourier Transform Box')

#displaying the inverse of the filtered Fourier box
plt.subplot(233)  
plt.imshow(np.abs(np.real(inverse_filtered_box)))
plt.title('Inverse Filtered Box')

#displaying the Fourier transform of the circle
plt.subplot(234)
plt.imshow(np.log(np.abs(fourier_circle_shifted) + 1))
plt.title('2D Fourier Transform Circle')

#displaying the Fourier transform of the circle with Gaussian filter
plt.subplot(235)
plt.imshow(np.log(filtered_fourier_circle + 1))
plt.title('Gaussian Filtered 2D Fourier Transform Circle')

#displaying the inverse of the filtered Fourier circle
plt.subplot(236)  
plt.imshow(np.abs(np.real(inverse_filtered_circle)))
plt.title('Inverse Filtered Circle')

The image I attached are the results I'm obtaining, the issue pertains to the inverse graphs.

Answer 1

After reading your code and trying to compare each step, I managed to identify a problem with your code: apply_gaussian_filter(np.abs(fourier_circle_shifted), sigma_circle). Applying the gaussian on the absolute of fourier_circle_shifted will make you lose phase information, hence reconstruction would not work. Here are two functions that apply the Gaussian on the whole fft or on the absolute of it:

def filterAndReconstruct_Abs(im, sigma = 1):
    fftOfIm = np.fft.fft2(im) # fft2 on image
    fftOfIm_Shifted = np.fft.fftshift(fftOfIm) # shift
    fftOfIm_Filtered = gaussian_filter(np.abs(fftOfIm_Shifted), sigma) # gaussian on absolute!!
    imReconstructed = np.fft.ifft2(np.fft.ifftshift(fftOfIm_Filtered)).real # reconstruct
    return imReconstructed
def filterAndReconstruct_WithPhase(im, sigma = 1):
    fftOfIm = np.fft.fft2(im) # fft2 on image
    fftOfIm_Shifted = np.fft.fftshift(fftOfIm) # shift
    fftOfIm_Filtered = gaussian_filter(fftOfIm_Shifted, sigma) # gaussian on fft with phase!!
    imReconstructed = np.fft.ifft2(np.fft.ifftshift(fftOfIm_Filtered)).real # reconstruct
    return imReconstructed

And now let's give it an example of a rectangle:

im = np.zeros((1000, 1000)) # create empty array
im[200:800, 200:800] = 255 # draw a rectangle
sigma = 0.2 # define sigma
fig, axs = plt.subplots(ncols = 3, nrows = 1, sharex = True, sharey = True)
axs[0].imshow(im)
axs[1].imshow(filterAndReconstruct_Abs(im,sigma)+1, norm=LogNorm())
axs[2].imshow(filterAndReconstruct_WithPhase(im,sigma)+1, norm=LogNorm())
for ax in axs:
    ax.axis("off")

The results will be like this:

As you can see, the reconstruction with the FFT is the correct variation. Playing with your example, if you use apply_gaussian_filter(fourier_circle_shifted, sigma_circle) and a lower sigma,you will get some nice results as well.

Hope this helps you further, also the imports of my example:

%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter
from matplotlib.colors import LogNorm