理論
傅里葉變換用來分析多種過濾器的頻率特征。對于圖片，2D離散傅里葉變換(DFT)用來找頻率范圍。一個(gè)快速算法叫快速傅里葉變換（FFT）用來計(jì)算DFT植旧。
對于正弦信號，x(t) = Asin(2πft).我們可以說f是信號的頻率离唐，如果頻率范圍給定病附，我們可以看到f的峰值。如果信號是從離散信號采樣亥鬓，我們還是一樣的頻率范圍完沪，但是是周期范圍[-π, π] 或[0, 2π]。你可以認(rèn)為一個(gè)圖像是從兩個(gè)方向采樣的信號嵌戈。所以在X和Y方向做傅里葉變換得到圖像的頻率表現(xiàn)覆积。
對于正弦信號，如果振幅變化很快熟呛，你可以說它是高頻信號宽档。如果變化很慢，就是低頻信號庵朝。這個(gè)也可以擴(kuò)展到圖像領(lǐng)域吗冤，圖像里的振幅變化又厉，在邊緣，或者噪點(diǎn)變化最大欣孤。所以我們可以說，邊緣和噪點(diǎn)是圖像里的高頻內(nèi)容昔逗。如果在振幅上沒特別大變化降传，就是低頻內(nèi)容
我們來看怎么做傅里葉變換
Numpy 里的傅里葉變換
首先我們來看怎么在Numpy里找傅里葉變換。Numpy有一個(gè)FFT包來做這個(gè)勾怒。np.fft.fft2()讓我們做頻率變換婆排。第一個(gè)參數(shù)是輸入圖像，是灰度的笔链，第二個(gè)參數(shù)是可選的段只，決定輸出數(shù)組的大小。如果比輸入圖像大的話鉴扫，輸入圖像會(huì)補(bǔ)全0然后再做轉(zhuǎn)換赞枕。如果比輸入圖像小，輸入圖像會(huì)被裁切坪创，如果不傳炕婶，輸出圖像大小和輸入一樣
當(dāng)你得到了結(jié)果，0頻率內(nèi)容會(huì)在左上角莱预，如果你想拿到中間來柠掂，你需要N/2在兩個(gè)方向上來移動(dòng)結(jié)果。這個(gè)是用函數(shù)np.fft.fftshift().當(dāng)你找到了頻率轉(zhuǎn)換依沮，你可以找到振幅譜線涯贞。
import cv2
import numpy as np
from matplotlib import pyplot as pltimg = cv2.imread('messi5.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
你可以看到白色區(qū)域在中心，表示低頻內(nèi)容更多危喉。
現(xiàn)在你達(dá)到了頻率轉(zhuǎn)換宋渔，就可以在頻率范圍內(nèi)做更多運(yùn)算了，比如高通過濾和重建圖像辜限，比如反向離散傅里葉變換傻谁。你只需要通過一個(gè)60x60的矩形窗口去掉低頻內(nèi)容，然后用np.fft.ifftshift()做反向變換列粪，0頻率內(nèi)容又在左上角了审磁。然后做反向快速傅里葉變換 np.ifft2()。結(jié)果是復(fù)數(shù)岂座。
rows, cols = img.shape
crow,ccol = rows/2 , cols/2
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap ='gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])plt.show()
結(jié)果是
結(jié)果顯示高通濾波器是一個(gè)邊緣檢測運(yùn)算态蒂，這個(gè)我們在圖片地圖的章節(jié)看過，這還顯示了大多數(shù)凸顯該數(shù)據(jù)是在光譜的低頻區(qū)域费什。
如果你仔細(xì)看結(jié)果钾恢，特別是JET色下的圖像手素，你可以看到一些（紅色箭頭指向的位置），這些波紋狀的結(jié)構(gòu)叫做ringing effects(振鈴效應(yīng))瘩蚪。這是由于我們用來掩圖的矩形窗口導(dǎo)致的泉懦。所以矩形窗口不用來做過濾，更好的選擇是高斯窗口
OpenCV里的傅里葉變換
OpenCV提供了函數(shù)cv2.dft()和cv2.idft()疹瘦。它返回和前面同樣的結(jié)果但是是雙通道的崩哩。第一個(gè)通道包含了結(jié)果的實(shí)數(shù)部分，第二個(gè)通道包含虛數(shù)部分言沐，輸入圖像應(yīng)該先被轉(zhuǎn)換成np.float32邓嘹。我們來看看：
import numpy as np
import cv2
from matplotlib import pyplot as pltimg = cv2.imread('messi5.jpg',0)
dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
你也可以用cv2.cartToPolar()同時(shí)返回magnitude和phase
現(xiàn)在我們改做反向DFT了，在前面的例子里我們創(chuàng)建了HPF险胰，這次我們看如何用LPF來去除圖像里的高頻內(nèi)容汹押，它實(shí)際上會(huì)模糊圖片。我們用高值（1）在低頻創(chuàng)建一個(gè)掩圖起便，傳入低頻內(nèi)容棚贾，在高頻區(qū)域用0.
rows, cols = img.shape
crow,ccol = rows/2 , cols/2# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1# apply mask and inverse DFT
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()
注意：
和平常一樣，OpenCV函數(shù)cv2.dft()和cv2.idft()比Numpy要快榆综，但是Numpy函數(shù)更友好鸟悴。
DFT性能優(yōu)化
DFT計(jì)算的性能在數(shù)組大小在有些情況下會(huì)好一些，當(dāng)數(shù)組大小是2的乘方的時(shí)候是最快的奖年。大小是2,3,5的城際是也處理的比較有效细诸。所以如果你擔(dān)心你的代碼的性能，你可以修改數(shù)組的大小到某些優(yōu)化的大新亍（通過補(bǔ)0）震贵，之后再做DFT，對于OpenCV水评，你必須手動(dòng)補(bǔ)0猩系，對于Numpy，你指定FFT計(jì)算的大小中燥，它會(huì)自動(dòng)補(bǔ)0.
我們怎么找優(yōu)化大小呢寇甸？OpenCV提供了函數(shù)cv2.getOptimalDFTSize()。它可以給cv2.dft()和np.fft.fft2()用疗涉，我們可以用IPython魔法命令%timeit
In [16]: img = cv2.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print rows,cols
342 548In [19]: nrows = cv2.getOptimalDFTSize(rows)
In [20]: ncols = cv2.getOptimalDFTSize(cols)
In [21]: print nrows, ncols
360 576
(342,548)被修改成了(360,576)∧妹梗現(xiàn)在用0補(bǔ)全（對于OpenCV）然后看計(jì)算DFT的性能。你可以創(chuàng)建一個(gè)新的大的全0的數(shù)組咱扣，然后吧數(shù)據(jù)拷貝進(jìn)去绽淘，使用cv2.copyMakeBorder()。
nimg=np.zeros((nrows,ncols))
nimg[:rows,:cols]=img
或者
right = ncols - cols
bottom = nrows - rows
bordertype = cv2.BORDER_CONSTANT ? ?#just to avoid line breakup in PDF file
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype,value=0)
現(xiàn)在我們計(jì)算DFT性能和Numpy函數(shù)對比
In [22]: %timeit fft1 = np.fft.fft2(img)
10 loops, best of 3: 40.9 ms per loop
In [23]: %timeit fft2 = np.fft.fft2(img,[nrows,ncols])
100 loops, best of 3: 10.4 ms per loop
顯示有4倍的性能差別闹伪，現(xiàn)在我們用OpenCV的函數(shù)
In [24]: %timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop
也是4倍的差別沪铭。你還可以看到OpenCV函數(shù)比Numpy的函數(shù)快3倍壮池。在反向FFT也是一樣。
為什么拉普拉斯是高通濾波器杀怠？
有很多人提類似的問題椰憋，為什么拉普拉斯是一個(gè)高通濾波器？為什么Sobel是高通濾波器等赔退，第一個(gè)回答是通過傅里葉變換橙依。對于一些高的FFT取拉普拉斯的傅里葉變換，分析它
import cv2
import numpy as np
from matplotlib import pyplot as plt# simple averaging filter without scaling parameter
mean_filter = np.ones((3,3))# creating a guassian filter
x = cv2.getGaussianKernel(5,10)
gaussian = x*x.T# different edge detecting filters
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
[-10,0,10],
[-3, 0, 3]])# sobel in x direction
sobel_x= np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]])# sobel in y direction
sobel_y= np.array([[-1,-2,-1],
[0, 0, 0],
[1, 2, 1]])# laplacian
laplacian=np.array([[0, 1, 0],
[1,-4, 1],
[0, 1, 0]])filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
'sobel_y', 'scharr_x']fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]for i in xrange(6):
? ? plt.subplot(2,3,i+1),plt.imshow(mag_spectrum[i],cmap = 'gray')
? ? plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])plt.show()
從圖像里你能看到每個(gè)kernel的屏蔽的頻率區(qū)域离钝，什么區(qū)域被通過了票编。從這個(gè)信息我們可以知道哪個(gè)kernel是HPF哪個(gè)是LPF褪储。
