Image in Frequency Domain

By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Source

• An Introduction to Digital Image Processing with MATLAB
• by Alasdair McAndrew at Victoria University of Technology

Table of Contents

0. Filter on Images

0.1. Linear System: Convolution

• Computational photography at Gatech

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import cv2
from skimage.util import random_noise

%matplotlib inline

# noise image

img = cv2.imread('.\image_files\lena_sigma25.png', 0)

print(img.shape)

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

plt.imshow(img, cmap = 'gray')
plt.axis('off')
plt.show()

(512, 512)

M = np.ones([3,3])/9
print(M)

[[0.11111111 0.11111111 0.11111111]
 [0.11111111 0.11111111 0.11111111]
 [0.11111111 0.11111111 0.11111111]]

# smoothing or noise reduction

img_conv = signal.convolve2d(img, M, 'valid')

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

plt.subplot(1,2,1)
plt.imshow(img, cmap = 'gray')
plt.title('Noisy Image')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(img_conv, cmap = 'gray')
plt.title('Smoothed Image')
plt.axis('off')

plt.show()

# original image

img = cv2.imread('.\image_files\lena.png', 0)

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

plt.imshow(img, cmap = 'gray')
plt.axis('off')
plt.show()

# guess what kind of image will be produced after convolution

Mv  = [[-1, 0, 1],
       [-2, 0, 2],
       [-1, 0, 1]]
Mv2 = [[1, 0, -1],
       [2, 0, -2],
       [1, 0, -1]]

img_conv = signal.convolve2d(img, Mv, 'valid')
img_conv2 = signal.convolve2d(img, Mv2, 'valid')

plt.figure(figsize = (12,6))

plt.subplot(1,2,1)
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.title('Mv')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(img_conv2, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv2))
plt.title('Mv2')
plt.axis('off')

plt.show()       

# guess what kind of image will be produced after convolution

Mh = [[-1, -2, -1],
      [0, 0, 0],
      [1, 2, 1]]

img_conv = signal.convolve2d(img, Mh, 'valid')

plt.figure(figsize = (10,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')

plt.show()

M = np.array(Mv + Mh)/2

img_conv = signal.convolve2d(img, M, 'valid')

plt.figure(figsize = (10,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')

plt.show()

M = [[0, -1, 0],
     [-1, 4, -1],
     [0, -1, 0]]

img_conv = signal.convolve2d(img, M, 'valid')

plt.figure(figsize = (10,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')

plt.show()

0.2. Non-linear System: Median Filter

$$y[n] = \text{median}\{x[n-k],\cdots,x[n+k]\}$$

There are nonlinear neighborhood operations that can be perfromed for the purpose of noise reduction that can do a better job of preserving edges than simple smoothing filters.

Median filters can do an excellent job of rejecting certain types of noise, in particular, "shot" or impulse noise (outlier in a time series) in which some individual pixels or signals have extreme values.

Example - removing shot noise (salt & pepper noise) in image

• Image and video processing at Duke (https://www.youtube.com/watch?v=xFaddafLbcg)

img = cv2.imread('.\image_files\lena.png', 0)
img_shot = random_noise(img, mode = 's&p')

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

plt.subplot(1,2,1) # original image
plt.imshow(img, cmap = 'gray', vmin = 0, vmax = np.amax(img))
plt.axis('off')

plt.subplot(1,2,2) # Degraded image with Salt & Pepper noise
plt.imshow(img_shot, cmap = 'gray', vmin = 0, vmax = np.amax(img_shot))
plt.axis('off')

plt.show()

# see the performance

M = np.ones([3,3])/9

img_avg = signal.convolve2d(M, img_shot, 'valid')

plt.figure(figsize = (10,6))
plt.imshow(img_avg, cmap = 'gray', vmin = 0, vmax = np.amax(img_avg))
plt.axis('off')

plt.show()

# Median Filter

img_median = signal.medfilt(img_shot,5)

plt.figure(figsize = (10,6))
plt.imshow(img_median, cmap = 'gray', vmin = 0, vmax = np.amax(img_median))
plt.axis('off')

plt.show()

1. Fourier Transforms of Images

1.1. Synthetic Images

# generate an image

N = 2**8
a = np.hstack((np.zeros((N,int(N/2))), np.ones((N,int(N/2)))))
a = np.asmatrix(a)

plt.figure(figsize = (10,6))
plt.imshow(a, cmap = 'gray')

# just for plotting a boundary

m = a.shape[0]

w = [0, 0]
x = [0, m-1]
y = [m-1, m-1]
z = [m-1, 0]

boundary = np.vstack([w, x, y, z])
boundary = np.asmatrix(boundary)

plt.plot(boundary[:,0], boundary[:,1], 'k-')

plt.show()

# Fourier Transformation with Image

Ak = np.fft.fft2(a)
As = np.fft.fftshift(Ak)

plt.figure(figsize = (10,6))
plt.imshow(np.log(1+np.absolute(As)), cmap = 'gray')
plt.axis('off')
plt.show()

L = np.absolute(As)

plt.figure(figsize = (10,6))
plt.plot(np.arange(L.shape[0]), L[int(N/2),:])
plt.xlim(int(N/2)-50, int(N/2)+50)
plt.ylim(bottom = 0)

plt.show()

N = 2**8
x, y = np.meshgrid(np.arange(N), np.arange(N))

k = 10
m = np.cos(2*np.pi/N*k*x)

M = np.fft.fft2(m)
Ms = np.fft.fftshift(M)

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

plt.subplot(1,2,1)
plt.imshow(m, cmap = 'gray')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(np.log(1+np.absolute(Ms)), cmap = 'gray')
plt.axis('off')

plt.show()

N = 2**8
x, y = np.meshgrid(np.arange(N), np.arange(N))

k = 5
m = np.cos(2*np.pi/N*k*y)

M = np.fft.fft2(m)
Ms = np.fft.fftshift(M)

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

plt.subplot(1,2,1)
plt.imshow(m, cmap = 'gray')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(np.log(1+np.absolute(Ms)), cmap = 'gray')
plt.axis('off')

plt.show()

N = 2**8
x, y = np.meshgrid(np.arange(N), np.arange(N))

k = 7
m = np.cos(2*np.pi/N*k*(x-y))

M = np.fft.fft2(m)
Ms = np.fft.fftshift(M)

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

plt.subplot(1,2,1)
plt.imshow(m, cmap = 'gray')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(np.log(1+np.absolute(Ms)), cmap = 'gray')
plt.axis('off')

plt.show()

# sqaure (box) image and its fft

N = 2**8
a = np.zeros([N,N])
a[int(N/2)-15 : int(N/2)+15, int(N/2)-15 : int(N/2)+15] = 1

As = np.fft.fftshift(np.fft.fft2(a))

# inverse fft (recovered image)

Ar = np.fft.ifftshift(As)
ar = np.real(np.fft.ifft2(Ar))

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original')
plt.imshow(a, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(As)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Recovered')
plt.imshow(ar, cmap = 'gray')
plt.axis('off')

plt.show()

# a box rotated 45 degree

x, y = np.meshgrid(np.arange(N), np.arange(N))
b = (x + y < 329) & (x + y > 182) & (x - y > -67) & (x - y < 73)
Bs = np.fft.fftshift(np.fft.fft2(b))

# inverse fft

Br = np.fft.ifftshift(Bs)
br = np.fft.ifft2(Br)

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original')
plt.imshow(b, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(Bs)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Recovered')
plt.imshow(np.real(br), cmap = 'gray')
plt.axis('off')

plt.show()

# circle

x, y = np.meshgrid(np.arange(-int(N/2), int(N/2), 1), np.arange(-int(N/2), int(N/2), 1))
z = np.sqrt(x**2 + y**2)

c = z < 15
Cs = np.fft.fftshift(np.fft.fft2(c))

# inverse fft

Cr = np.fft.ifftshift(Cs)
cr = np.fft.ifft2(Cr)

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original')
plt.imshow(c, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(Cs)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Recovered')
plt.imshow(np.real(cr), cmap = 'gray')
plt.axis('off')

plt.show()

Note "ringing" in the Fourier transform. This artifact is associated with the sharp cutoff of the circle. Think of a sinc function.

1.2. Circle with a Gentle Cutoff

Butterworh

$$ \frac{1}{1+\alpha(x^2+y^2)} $$

x, y = np.meshgrid(np.arange(-int(N/2), int(N/2), 1), np.arange(-int(N/2), int(N/2), 1))
z = np.sqrt(x**2 + y**2)

d = 1/(1 + (z/15)**2)
Ds = np.fft.fftshift(np.fft.fft2(d))

# inverse fft

Dr = np.fft.ifftshift(Ds)
dr = np.fft.ifft2(Dr)

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original')
plt.imshow(d, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(Ds)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Recovered')
plt.imshow(np.real(dr), cmap = 'gray')
plt.axis('off')

plt.show()

$$e^{\left(-\alpha_1 x^2 - \alpha_2 y^2\right)}$$

# Gaussian in time

x, y = np.meshgrid(np.arange(-int(N/2), int(N/2), 1), np.arange(-int(N/2), int(N/2), 1))
z = np.exp(- 0.01*x**2 - 0.05*y**2)

# Gaussain in freq.

Cs = np.fft.fftshift(np.fft.fft2(z))

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

plt.subplot(1,2,1)
plt.title('Original')
plt.imshow(z, cmap = 'gray')
plt.axis('off')

plt.subplot(1,2,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(Cs)), cmap = 'gray')
plt.axis('off')

plt.show()

1.3. Real Images

• Download the image

# fft with an image

cm = cv2.imread('.\image_files\cameraman.png', 0)

CMs = np.fft.fftshift(np.fft.fft2(cm))

CMr = np.fft.ifftshift(CMs)
cmr = np.fft.ifft2(CMr)

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original')
plt.imshow(cm, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(CMs)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Recovered')
plt.imshow(np.real(cmr), cmap = 'gray')
plt.axis('off')

plt.show()

2. Image Filtering in Frequency

2.1. Low-Pass Filtering in Frequency

• ideal filter

x, y = np.meshgrid(np.arange(-int(N/2), int(N/2), 1), np.arange(-int(N/2), int(N/2), 1))
z = np.sqrt(x**2 + y**2)

C = z < 15

CMFs = CMs*C

CMFr = np.fft.ifftshift(CMFs)
cmfr = np.fft.ifft2(CMFr)

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

plt.subplot(2,2,1)
plt.title('Original Image')
plt.imshow(cm, cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(CMs)), cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,3)
plt.title('Filtered image')
plt.imshow(np.real(cmfr), cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,4)
plt.title('Filter in Freq.')
plt.imshow(np.log(1+np.absolute(C)), cmap = 'gray')
plt.axis('off')

plt.show()

• Butterworth $ \frac{1}{1+\alpha(x^2+y^2)} $ type

N = 2**8
x, y = np.meshgrid(np.arange(N), np.arange(N))
z = np.sqrt((x - N/2)**2 + (y - N/2)**2)

D = 1/(1 + (z/10)**2) # filter in the frequency domain

CMFs = CMs*D

CMFr = np.fft.ifftshift(CMFs)
cmfr = np.fft.ifft2(CMFr)

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original Image')
plt.imshow(cm, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(CMs)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Filtered Image')
plt.imshow(np.real(cmfr), cmap = 'gray')
plt.axis('off')

plt.show()

• Guassian type

x, y = np.meshgrid(np.arange(-int(N/2), int(N/2), 1), np.arange(-int(N/2), int(N/2), 1))
Z = np.exp(- 0.0005*x**2 - 0.0005*y**2) # Gaussian

CMFs = CMs*Z

CMFr = np.fft.ifftshift(CMFs)
cmfr = np.fft.ifft2(CMFr)

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original Image')
plt.imshow(cm, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('Filtered FFT')
plt.imshow(np.log(1+np.absolute(CMFs)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Filtered Image')
plt.imshow(np.real(cmfr), cmap = 'gray')
plt.axis('off')

plt.show()

Finding Abraham Lincoln

• Gala Contemplating the Mediterranean Sea Which at Twenty Meters Becomes the Portrait of Abraham Lincoln
• http://thedali.org/exhibit/gala-contemplating-mediterranean-sea/

• Download the image

imbgr = cv2.imread('./image_files/Dali.jpg') # bgr (color) image
b, g, r = cv2.split(imbgr)
imrgb = cv2.merge([r, g, b]) # change bgr image to rgb image
imgray = cv2.imread('./image_files/Dali.jpg', 0) # gray image

plt.figure(figsize = (12,6))

plt.subplot(1,2,1)
plt.title('Colored Original Painting')
plt.imshow(imrgb)
plt.axis('off')

plt.subplot(1,2,2)
plt.title('Grayed Painting')
plt.imshow(imgray, cmap = 'gray')
plt.axis('off')

plt.show()

IM = np.fft.fft2(imgray)
IMs = np.fft.fftshift(IM)

# Gaussian

N = IMs.shape[0]
M = IMs.shape[1]
x, y = np.meshgrid(np.arange(-int(M/2), int(M/2), 1), np.arange(-int(N/2), int(N/2), 1))
Z = np.exp(- 0.005*x**2 - 0.005*y**2)

IMFs = IMs*Z

IMFr = np.fft.ifftshift(IMFs)
imfr = np.fft.ifft2(IMFr)

plt.figure(figsize = (12,6))

plt.subplot(1,2,1)
plt.title('Grayed Painting')
plt.imshow(imgray, cmap = 'gray')
plt.axis('off')

plt.subplot(1,2,2)
plt.title('Blurred Image')
plt.imshow(np.real(imfr), cmap = 'gray')
plt.axis('off')

plt.show()

2.2. High-Pass Filter in Freq.

• sharpen (or shows the edges of) an image

N = 2**8
x, y = np.meshgrid(np.arange(-int(N/2), int(N/2), 1), np.arange(-int(N/2), int(N/2), 1))
z = np.sqrt(x**2 + y**2)

C = z > 40
CMFs = CMs*C

CMFr = np.fft.ifftshift(CMFs)
cmfr = np.fft.ifft2(CMFr)

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

plt.subplot(2,2,1)
plt.title('Original Image')
plt.imshow(cm, cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(CMs)), cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,3)
plt.title('Filtered Image')
plt.imshow(np.real(cmfr), cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,4)
plt.title('Filter in Freq.')
plt.imshow(np.log(1+np.absolute(C)), cmap = 'gray')
plt.axis('off')

plt.show()

$$1-e^{\left(-\alpha_1 x^2 - \alpha_2 y^2\right)}$$

x, y = np.meshgrid(np.arange(-int(N/2), int(N/2), 1), np.arange(-int(N/2), int(N/2), 1))
Z = 1 - np.exp(- 0.0005*x**2 - 0.0005*y**2) # Gaussian

CMFs = CMs*Z

CMFr = np.fft.ifftshift(CMFs)
cmfr = np.fft.ifft2(CMFr)

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original Image')
plt.imshow(cm, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('Filtered FFT')
plt.imshow(np.log(1+np.absolute(CMFs)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Filtered Image')
plt.imshow(np.real(cmfr), cmap = 'gray')
plt.axis('off')

plt.show()

$$ 1- \frac{1}{1+\alpha(x^2+y^2)} $$

N = 2**8
x, y = np.meshgrid(np.arange(N), np.arange(N))
z = np.sqrt((x - N/2)**2 + (y - N/2)**2)

D = 1 - 1/(1 + (z/10)**2) # filter in the frequency domain

CMFs = CMs*D

CMFr = np.fft.ifftshift(CMFs)
cmfr = np.fft.ifft2(CMFr)

plt.figure(figsize = (20,6))

plt.subplot(1,3,1)
plt.title('Original Image')
plt.imshow(cm, cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,2)
plt.title('Filtered FFT')
plt.imshow(np.log(1+np.absolute(CMFs)), cmap = 'gray')
plt.axis('off')

plt.subplot(1,3,3)
plt.title('Filtered Image')
plt.imshow(np.real(cmfr), cmap = 'gray')
plt.axis('off')

plt.show()

2.3. Edge Detection

• https://www.youtube.com/watch?v=JQoOL2v4YSE

cm = cv2.imread('.\image_files\cameraman.png', 0)
cme = cv2.Canny(cm, 100, 200)

plt.figure(figsize = (10,10))
plt.imshow(cme, cmap = 'gray')
plt.axis('off')
plt.show()

cmd = cv2.dilate(cme, np.ones([2,2]), iterations = 1)

plt.figure(figsize = (10,10))
plt.imshow(cmd, cmap = 'gray')
plt.axis('off')
plt.show()

2.4. Notch Filter

• band-stop filter
• periodic noise reduction

• Download the image

# load a test image (with patterned noises)

im = cv2.imread('.\image_files\ClownOrig.jpg', 0)

plt.figure(figsize = (10,6))
plt.imshow(im, cmap = 'gray')
plt.axis('off')
plt.show()

# transform to Fourier Domain

IM = np.fft.fft2(im)
IMs = np.fft.fftshift(IM)

plt.figure(figsize = (10,6))
plt.imshow(np.log(1+np.absolute(IMs)), cmap = 'gray')
plt.axis('off')
plt.show()

# run it in Python for interactive demo

N = IMs.shape[0]
x, y = np.meshgrid(np.arange(N), np.arange(N))

# notch filter generation (need to understand)

a1 = 0.008
a2 = 0.008

NF1 = 1 - np.exp(-a1*(x-190)**2 - a2*(y-123)**2) # Gaussian
NF2 = 1 - np.exp(-a1*(x-104)**2 - a2*(y-172)**2) # Gaussian
NF3 = 1 - np.exp(-a1*(x-126)**2 - a2*(y-135)**2) # Gaussian
NF4 = 1 - np.exp(-a1*(x-168)**2 - a2*(y-161)**2) # Gaussian

Z = NF1*NF2*NF3*NF4
IMFs = IMs*Z

IMFr = np.fft.ifftshift(IMFs)
imfr = np.fft.ifft2(IMFr)

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

plt.subplot(2,2,1)
plt.title('Original Image')
plt.imshow(im, cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,2)
plt.title('Filter in Freq')
plt.imshow(np.log(1+np.absolute(Z)), cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,3)
plt.title('Filtered Image')
plt.imshow(np.real(imfr), cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,4)
plt.title('Filtered FFT')
plt.imshow(np.log(1+np.absolute(IMFs)), cmap = 'gray')
plt.axis('off')

plt.show()

• Download the image

# load a test image (with patterned noises)

im = cv2.imread('.\image_files\lena2.jpg', 0)

IM = np.fft.fft2(im)
IMs = np.fft.fftshift(IM)

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

plt.subplot(1,2,1)
plt.title('Original')
plt.imshow(im, cmap = 'gray')
plt.axis('off')

plt.subplot(1,2,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(IMs)), cmap = 'gray')
plt.axis('off')

plt.show()

• Download the image

# load a test image (with patterned noises)

im = cv2.imread('.\image_files\script.jpg', 0)

IM = np.fft.fft2(im)
IMs = np.fft.fftshift(IM)

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

plt.subplot(1,2,1)
plt.title('Original')
plt.imshow(im, cmap = 'gray')
plt.axis('off')

plt.subplot(1,2,2)
plt.title('FFT')
plt.imshow(np.log(1+np.absolute(IMs)), cmap = 'gray')
plt.axis('off')

plt.show()

• Download the image

im = cv2.imread('.\image_files\halftone.png', 0)

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

plt.title('Original')
plt.imshow(im, cmap = 'gray')
plt.axis('off')

plt.show()

IM = np.fft.fft2(im)
IMs = np.fft.fftshift(IM)

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

plt.title('FFT')
plt.imshow(np.log(1+np.absolute(IMs)), cmap = 'gray')
plt.axis('off')

plt.show()

N = IMs.shape[0]
x, y = np.meshgrid(np.arange(N), np.arange(N))

# notch filter generation (need to understand)

a1 = 0.005
a2 = 0.005

NF1  = 1 - np.exp(-a1*(x-400)**2 - a2*(y-134)**2) # Gaussian
NF2  = 1 - np.exp(-a1*(x-400)**2 - a2*(y-267)**2) # Gaussian
NF3  = 1 - np.exp(-a1*(x-400)**2 - a2*(y-400)**2) # Gaussian
NF4  = 1 - np.exp(-a1*(x-400)**2 - a2*(y-533)**2) # Gaussian
NF5  = 1 - np.exp(-a1*(x-266)**2 - a2*(y-134)**2) # Gaussian
NF6  = 1 - np.exp(-a1*(x-266)**2 - a2*(y-267)**2) # Gaussian
NF7  = 1 - np.exp(-a1*(x-266)**2 - a2*(y-400)**2) # Gaussian
NF8  = 1 - np.exp(-a1*(x-266)**2 - a2*(y-533)**2) # Gaussian

NF9  = 1 - np.exp(-a1*(x-533)**2 - a2*(y-134)**2) # Gaussian
NF10 = 1 - np.exp(-a1*(x-533)**2 - a2*(y-267)**2) # Gaussian
NF11 = 1 - np.exp(-a1*(x-533)**2 - a2*(y-400)**2) # Gaussian
NF12 = 1 - np.exp(-a1*(x-533)**2 - a2*(y-533)**2) # Gaussian
NF13 = 1 - np.exp(-a1*(x-133)**2 - a2*(y-134)**2) # Gaussian
NF14 = 1 - np.exp(-a1*(x-133)**2 - a2*(y-267)**2) # Gaussian
NF15 = 1 - np.exp(-a1*(x-133)**2 - a2*(y-400)**2) # Gaussian
NF16 = 1 - np.exp(-a1*(x-133)**2 - a2*(y-533)**2) # Gaussian

NF17 = 1 - np.exp(-a1*(x-66)**2 - a2*(y-65)**2)   # Gaussian
NF18 = 1 - np.exp(-a1*(x-66)**2 - a2*(y-200)**2)  # Gaussian
NF19 = 1 - np.exp(-a1*(x-66)**2 - a2*(y-334)**2)  # Gaussian
NF20 = 1 - np.exp(-a1*(x-66)**2 - a2*(y-468)**2)  # Gaussian
NF21 = 1 - np.exp(-a1*(x-66)**2 - a2*(y-602)**2)  # Gaussian

NF22 = 1 - np.exp(-a1*(x-200)**2 - a2*(y-65)**2)  # Gaussian
NF23 = 1 - np.exp(-a1*(x-200)**2 - a2*(y-200)**2) # Gaussian
NF24 = 1 - np.exp(-a1*(x-200)**2 - a2*(y-334)**2) # Gaussian
NF25 = 1 - np.exp(-a1*(x-200)**2 - a2*(y-468)**2) # Gaussian
NF26 = 1 - np.exp(-a1*(x-200)**2 - a2*(y-602)**2) # Gaussian

NF27 = 1 - np.exp(-a1*(x-334)**2 - a2*(y-65)**2)  # Gaussian
NF28 = 1 - np.exp(-a1*(x-334)**2 - a2*(y-200)**2) # Gaussian
NF29 = 1 - np.exp(-a1*(x-334)**2 - a2*(y-334)**2) # Gaussian
NF30 = 1 - np.exp(-a1*(x-334)**2 - a2*(y-468)**2) # Gaussian
NF31 = 1 - np.exp(-a1*(x-334)**2 - a2*(y-602)**2) # Gaussian

NF32 = 1 - np.exp(-a1*(x-466)**2 - a2*(y-65)**2)  # Gaussian
NF33 = 1 - np.exp(-a1*(x-466)**2 - a2*(y-200)**2) # Gaussian
NF34 = 1 - np.exp(-a1*(x-466)**2 - a2*(y-334)**2) # Gaussian
NF35 = 1 - np.exp(-a1*(x-466)**2 - a2*(y-468)**2) # Gaussian
NF36 = 1 - np.exp(-a1*(x-466)**2 - a2*(y-602)**2) # Gaussian

NF37 = 1 - np.exp(-a1*(x-600)**2 - a2*(y-65)**2)  # Gaussian
NF38 = 1 - np.exp(-a1*(x-600)**2 - a2*(y-200)**2) # Gaussian
NF39 = 1 - np.exp(-a1*(x-600)**2 - a2*(y-334)**2) # Gaussian
NF40 = 1 - np.exp(-a1*(x-600)**2 - a2*(y-468)**2) # Gaussian
NF41 = 1 - np.exp(-a1*(x-600)**2 - a2*(y-602)**2) # Gaussian

# Z = NF1*NF2*NF3*NF4*NF5*NF6*NF7*NF8
# Z = Z*NF9*NF10*NF11*NF12*NF13*NF14*NF15*NF16
# Z = Z*NF17*NF18*NF19*NF20*NF21*NF22*NF23*NF24*NF25*NF26
# Z = Z*NF27*NF28*NF29*NF30*NF31*NF32*NF33*NF34*NF35*NF36
# Z = Z*NF37*NF38*NF39*NF40*NF41

Z = NF2*NF3*NF6*NF7
IMFs = IMs*Z

IMFr = np.fft.ifftshift(IMFs)
imfr = np.fft.ifft2(IMFr)

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

plt.subplot(2,2,1)
plt.title('Original Image')
plt.imshow(im, cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,2)
plt.title('Filter in Freq')
plt.imshow(np.log(1+np.absolute(Z)), cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,3)
plt.title('Filtered Image')
plt.imshow(np.real(imfr), cmap = 'gray')
plt.axis('off')

plt.subplot(2,2,4)
plt.title('Filtered FFT')
plt.imshow(np.log(1+np.absolute(IMFs)), cmap = 'gray')
plt.axis('off')

plt.show()

plt.figure(figsize = (10,10))
plt.imshow(np.real(imfr), cmap = 'gray')
plt.axis('off')
plt.show()

%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')