%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
import scipy.misc
face = scipy.misc.face()
image = face[:,:,1]
# image size, square side length, number of squares
ncols, nrows = 768, 1024
sq_size, nsq = 10, 20

plt.figure(),plt.imshow(image,cmap='gray')

# Take the 2-dimensional DFT and centre the frequencies
ftimage = np.fft.fft2(np.fft.fftshift(image))
ftimage = np.fft.fftshift(ftimage)
plt.figure(),plt.imshow(np.abs(ftimage),cmap='gray',vmin=0, vmax=10000)

# Build and apply a Gaussian filter.
sigmax, sigmay = 50, 50
cx,cy = nrows/2, ncols/2
x = np.linspace(0, nrows, nrows)
y = np.linspace(0, ncols, ncols)
X, Y = np.meshgrid(x, y)
gmask = np.exp(-(((X-cx)/sigmax)**2 + ((Y-cy)/sigmay)**2))

ftimagep = ftimage * gmask
plt.figure(),plt.imshow(gmask,cmap='gray')
plt.figure(),plt.imshow(np.abs(ftimagep),cmap='gray',vmin=0, vmax=10000)


# Finally, take the inverse transform and show the blurred image
imagep = np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(ftimagep)))
plt.figure(),plt.imshow(np.abs(imagep),cmap='gray')

ftimage_restore = ftimagep / gmask

plt.figure(),plt.imshow(gmask,cmap='gray')
plt.figure(),plt.imshow(np.abs(ftimage_restore),cmap='gray',vmin=0, vmax=10000)


# Finally, take the inverse transform and show the blurred image
imagep = np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(ftimage_restore)))
plt.figure(),plt.imshow(np.abs(imagep),cmap='gray')

!dir