Contact

Frequency Domain of Images – Fourier Transform and Filtering

frequential_donike
In [1]:
from skimage import io
from skimage import color
from skimage.restoration import denoise_nl_means, estimate_sigma

import numpy as np
from numpy.fft import fft, fftfreq, ifft

from scipy import ndimage as nd
from scipy.fft import fft, ifft
from scipy import fftpack

from PIL import Image, ImageDraw
import matplotlib.pyplot as plt
import time
import cv2
from math import sqrt
from math import exp

1.1 Load an image

In [2]:
img = io.imread("images/baseballdiamond08.tif")
plt.imshow(img)
plt.show()

1.2 Verify the invertibility of the FFT

In [3]:
# https://scipy-lectures.org/intro/scipy/auto_examples/solutions/plot_fft_image_denoise.html

# fft
img_fft =  fftpack.fft2(img)


# inverse of signal
img_ifft = fftpack.ifft2(img_fft).real

plt.figure(figsize=(10,10))

plt.subplot(131)
plt.imshow(img)
plt.title("Original")

plt.subplot(132)
plt.title("Spectrum")
plt.imshow(img_fft.astype(np.uint8))
plt.subplot(133)

plt.title("Reconstructed")
plt.imshow(img_ifft.astype(np.uint8))
plt.show()

/anaconda/lib/python3.7/site-packages/ipykernel_launcher.py:18: ComplexWarning: Casting complex values to real discards the imaginary part

The reconstruction based on the spectrum returns the original image.

1.3 Compare the effect of the main low-pass filters

Closely modeled (euphemism for copied) after this awesome Tutorial that helped me understand the process. Trying to implement these steps manually proved to be very difficult, with mixed/non-sensical results.

In [4]:
def rgb2gray(rgb):
    return np.dot(rgb[...,:3], [0.2989, 0.5870, 0.1140])
img = rgb2gray(img)

plt.imshow(img,cmap="gray")
plt.title("Gray")
plt.show()
In [5]:
plt.figure(figsize=(8, 6), constrained_layout=False)


img_fft = np.fft.fft2(img)
img_fftshift = np.fft.fftshift(img_fft)
img_ifftshit = np.fft.ifftshift(img_fftshift)
img_ifft = np.fft.ifft2(img_ifftshit)

plt.subplot(231), plt.imshow(img, "gray"), plt.title("Original Image")
plt.subplot(232), plt.imshow(np.log(1+np.abs(img_fft)), "gray"), plt.title("Spectrum, FFT")
plt.subplot(233), plt.imshow(np.log(1+np.abs(img_fftshift)), "gray"), plt.title("Centered Spectrum")
plt.subplot(234), plt.imshow(np.log(1+np.abs(img_ifftshit)), "gray"), plt.title("Decentralized IFFT")
plt.subplot(235), plt.imshow(np.abs(img_ifft), "gray"), plt.title("Reversed Image")

plt.show()
In [6]:
def distance(point1,point2):
    return sqrt((point1[0]-point2[0])**2 + (point1[1]-point2[1])**2)

def gaussianLP(D0,imgShape):
    base = np.zeros(imgShape[:2])
    rows, cols = imgShape[:2]
    center = (rows/2,cols/2)
    for x in range(cols):
        for y in range(rows):
            base[y,x] = exp(((-distance((y,x),center)**2)/(2*(D0**2))))
    return base
In [7]:
def try_d0s_lp(d0):
    plt.figure(figsize=(25, 5), constrained_layout=False)


    plt.subplot(161), plt.imshow(img, "gray"), plt.title("Original Image")

    original = np.fft.fft2(img)
    plt.subplot(162), plt.imshow(np.log(1+np.abs(original)), "gray"), plt.title("Spectrum")

    center = np.fft.fftshift(original)
    plt.subplot(163), plt.imshow(np.log(1+np.abs(center)), "gray"), plt.title("Centered Spectrum")

    LowPassCenter = center * gaussianLP(d0,img.shape)
    plt.subplot(164), plt.imshow(np.log(1+np.abs(LowPassCenter)), "gray"), plt.title("Centered Spectrum multiply Low Pass Filter")

    LowPass = np.fft.ifftshift(LowPassCenter)
    plt.subplot(165), plt.imshow(np.log(1+np.abs(LowPass)), "gray"), plt.title("Decentralize")

    inverse_LowPass = np.fft.ifft2(LowPass)
    plt.subplot(166), plt.imshow(np.abs(inverse_LowPass), "gray"), plt.title("Processed Image")

    plt.suptitle("D0:"+str(d0),fontweight="bold")
    plt.subplots_adjust(top=1.1)
    plt.show()
    
    
for i in [100,50,30,20,10]:
    try_d0s_lp(i)