#!/usr/bin/env python
# coding: utf-8
#
# ESF projekt Západočeské univerzity v Plzni reg. č. CZ.02.2.69/0.0/0.0/16 015/0002287
# # Fourier transform
# Spojitá fourierova tranformace
#
# $$ F(u,v)=\int \int_{-\infty}^{\infty} f(x,y) e^{-j2\pi(ux+vy)} dx dy $$
#
# Kde $u$ a $v$ jsou souřadnice ve frekvenčním prostoru a $x$ a $y$ jsou souřadnice v obrazovém prostoru. $F(u,v)$ je 2D spektrum obrazové funkce $F(x,y)$
#
# Diskrétní fourierova transformace
#
# $$ F(u,v)=\sum_{m=-\infty}^{\infty} \sum_{n=-\infty}^{\infty} f[m,n]e^{-j2\pi(umx_o+vny_o)} $$
# [Interactive 1D Fourier Transform](http://www.jezzamon.com/fourier/)
#
# [YouTube Video - 2D Fourier Transform](https://www.youtube.com/watch?v=qa1ZxK9Y1Tw)
#
# In[1]:
get_ipython().run_line_magic('pylab', 'inline')
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
import skimage.io
# Výstup fourierovy transformace je (relativně) souměrný podle středu a má běžně nízké frekvence na okrajích. Pro vizualizaci a lepší představu se kvadranty obrazu přehazují pomocí funkce fftshift()
#
# ![fft shift](http://i.stack.imgur.com/I6XDQ.gif)
# In[2]:
im = data.camera()
ft = np.fft.fft2(im)
ftshift = np.fft.fftshift(ft)
spek = 20*np.log(np.abs(ftshift))
# vizualizace
plt.subplot(121)
plt.imshow(im, cmap='gray')
plt.subplot(122)
plt.imshow(spek, cmap = 'gray')
plt.title('Spektrum')
# Example: Show fourier transform of the camera image.
# # Filtration
#
# ## Lo-pass
# In[15]:
x, y = np.indices(im.shape)
center1 = (256, 256)
radius1 = 20
mask = (x - center1[0])**2 + (y - center1[1])**2 < radius1**2
plt.imshow(mask)
# In[16]:
#ftshift * mask
ftshift_mask = ftshift*(mask)
ftishift_back = np.fft.ifftshift(ftshift_mask)
im_back = np.fft.ifft2(ftishift_back)
im_back = np.abs(im_back)
# Visualization
plt.figure(figsize=(9, 4))
plt.subplot(121)
plt.imshow(im_back, cmap='gray')
plt.subplot(122)
plt.imshow(np.log(np.abs(ftshift_mask)), cmap='gray')
# ## Hi-pass
# In[17]:
# Vytvoření filtru
x, y = np.indices(im.shape)
center2 = (256, 256)
radius2 = 80
mask2 = (x - center1[0])**2 + (y - center1[1])**2 < radius1**2
# Filtrace a zpětná FT
ftshift_mask = ftshift*(1-mask2)
ftishift_back = np.fft.ifftshift(ftshift_mask)
im_back2 = np.fft.ifft2(ftishift_back)
im_back2 = np.abs(im_back2)
# Visualization
plt.figure(figsize=(9, 4))
plt.subplot(121)
plt.imshow(im_back2, cmap='gray')
plt.subplot(122)
plt.imshow(np.log(np.abs(ftshift_mask)), cmap='gray')
# # Convolution theorem
#
# $$ f*g = FT^{-1} \left( FT(f) \cdot FT(G) \right) $$
# # Example
# Make restauration of the image
# http://www.imagemagick.org/Usage/fourier/twigs.jpg
#
# ![clown](http://www.imagemagick.org/Usage/fourier/twigs.jpg)
#
#
# In[33]:
clown = skimage.io.imread('http://www.imagemagick.org/Usage/fourier/clown.jpg')
twigs = skimage.io.imread('http://www.imagemagick.org/Usage/fourier/twigs.jpg')