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.
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)