OpenCV-Python官方文檔關(guān)于圖像傅里葉變換和反變換的教程網(wǎng)址：
https://docs.opencv.org/4.1.0/de/dbc/tutorial_py_fourier_transform.html?
目標(biāo)
我們將要學(xué)習(xí)：
? 使用 OpenCV 對(duì)圖像進(jìn)行傅里葉變換（DFT）：cv2.dft()，cv2.idft()
? 使用 Numpy 中 FFT（快速傅里葉變換）函數(shù)：
? 傅里葉變換的一些用處
? 我們將要學(xué)習(xí)的函數(shù)有： cv2.dft()，cv2.idft() 等
原理
把圖像想象成沿著兩個(gè)方向采集的信號(hào)。所以對(duì)圖像同時(shí)進(jìn)行 X 方向和 Y 方向的傅里葉變換，我們就會(huì)得到這幅圖像的頻域表示（頻譜圖）。更直觀一點(diǎn)，對(duì)于一個(gè)正弦信號(hào)，如果它的幅值變化非常快，我們可以說(shuō)它是高頻信號(hào)，如果變化非常慢，我們稱(chēng)之為低頻信號(hào)。把這種想法應(yīng)用到圖像中，圖像哪里的幅度變化非常大呢？邊界點(diǎn)或者噪聲。所以我們說(shuō) 邊界和噪聲是圖像中的高頻分量 。
1.1 Numpy 中的傅里葉變換?
Numpy 中的 FFT 包 可以幫助我們實(shí)現(xiàn)快速傅里葉變換。 函數(shù) np.fft.fft2() 可以對(duì)信號(hào)進(jìn)行頻率轉(zhuǎn)換，輸出結(jié)果是一個(gè)復(fù)雜的組。本函數(shù)的第一個(gè)參數(shù)是輸入圖像，要求是灰度格式。第二個(gè)參數(shù)是可選的, 決定輸出數(shù)組的大小。輸出數(shù)組的大小和輸入圖像大小一樣。如果輸出結(jié)果比輸入圖像大，輸入圖像就需要在進(jìn)行 FFT 前補(bǔ)0。如果輸出結(jié)果比輸入圖像小的話，輸入圖像就會(huì)被切割。?
現(xiàn)在我們得到了結(jié)果，頻率為 0 的部分（直流分量）在輸出圖像的左上角。如果想讓它（直流分量）在輸出圖像的中心，我們還需要將結(jié)果沿兩個(gè)方向平移 N/2 。 函數(shù) np.fft.fftshift() 可以幫助我們實(shí)現(xiàn)這一步。（這樣更容易分析）。進(jìn)行完頻率變換之后，我們就可以構(gòu)建振幅譜了。?
import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv2.imread('messi5.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# 這里構(gòu)建振幅圖的公式?jīng)]學(xué)過(guò)
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()
結(jié)果：
然后就可以在頻域?qū)D像進(jìn)行操作了，比如進(jìn)行高通濾波，再傅里葉反變換回空間域。比如我們可以使用一個(gè)60x60 的矩形窗口對(duì)圖像進(jìn)行掩模操作從而去除低頻分量。然后再使用 函數(shù)np.fft.ifftshift() 進(jìn)行逆平移操作，所以現(xiàn)在直流分量又回到左上角了，然后使用 函數(shù) np.ifft2() 進(jìn)行 FFT 逆變換。同樣又得到一堆復(fù)雜的數(shù)字，我們可以對(duì)他們?nèi)〗^對(duì)值：
rows, cols = img.shape
crow,ccol = int(rows/2) , int(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)
# 取絕對(duì)值
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é)果顯示高通濾波其實(shí)是一種邊界檢測(cè)操作。這就是我們?cè)谇懊鎴D像梯度那一章看到的。同時(shí)我們還發(fā)現(xiàn)圖像中的大部分?jǐn)?shù)據(jù)集中在頻譜圖的低頻區(qū)域。
如果你觀察仔細(xì)的話，尤其是最后一章 JET 顏色的圖像，你會(huì)看到一些不自然的東西（如我用紅色箭頭標(biāo)出的區(qū)域）。看上圖那里有些條帶裝的結(jié)構(gòu)，這被成為 振鈴效應(yīng) 。這是由于我們使用矩形窗口做掩模造成的。這個(gè)掩模被轉(zhuǎn)換成正弦形狀時(shí)就會(huì)出現(xiàn)這個(gè)問(wèn)題。所以 一般我們不使用矩形窗口濾波。最好的選擇是高斯窗口。
1.2 OpenCV 中的傅里葉變換
OpenCV 中相應(yīng)的函數(shù)是 cv2.dft() 和 cv2.idft()。和前面輸出的結(jié)果一樣，但是是雙通道的。第一個(gè)通道是結(jié)果的實(shí)數(shù)部分，第二個(gè)通道是結(jié)果的虛數(shù)部分。輸入圖像要首先轉(zhuǎn)換成 np.float32 格式。
import numpy as np
import cv2
from matplotlib import pyplot as plt
img = 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()
?注意： 你可以使用 函數(shù) cv2.cartToPolar() ，它會(huì)同時(shí)返回幅度和相位。
現(xiàn)在我們來(lái)做逆 DFT。在前面的部分我們實(shí)現(xiàn)了一個(gè) HPF（高通濾波），現(xiàn)在我們來(lái)做 LPF（低通濾波）將高頻部分去除。其實(shí)就是對(duì)圖像進(jìn)行模糊操作。首先我們需要構(gòu)建一個(gè)掩模，與低頻區(qū)域?qū)?yīng)的地方設(shè)置為 1, 與高頻區(qū)域?qū)?yīng)的地方設(shè)置為 0。
rows, cols = img.shape
crow,ccol = int(rows/2) , int(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()
?結(jié)果：
注意： OpenCV 中的函數(shù) cv2.dft() 和 cv2.idft() 要比 Numpy 快 。但是Numpy 函數(shù)更加用戶(hù)友好。關(guān)于性能的描述，請(qǐng)看下面的小節(jié)。
1.3 DFT 的性能優(yōu)化
當(dāng)數(shù)組的大小為某些值時(shí) DFT 的性能會(huì)更好。當(dāng)數(shù)組的大小是 2 的指數(shù)時(shí) DFT 效率最高。當(dāng)數(shù)組的大小是 2， 3， 5 的倍數(shù)時(shí)效率也會(huì)很高。所以如果你想提高代碼的運(yùn)行效率時(shí)，你可以修改輸入圖像的大小（補(bǔ) 0）。對(duì)于OpenCV 你必須自己手動(dòng)補(bǔ)0。但是 Numpy，你只需要指定 FFT 運(yùn)算的大小，它會(huì)自動(dòng)補(bǔ) 0。
那我們?cè)鯓哟_定最佳大小呢？OpenCV 提供了一個(gè) 函數(shù) cv2.getOptimalDFTSize() 。它可以同時(shí)被 cv2.dft() 和 np.fft.fft2() 使用。讓我們一起使用 IPython的魔法命令%timeit 來(lái)測(cè)試一下吧。
In [16]: img = cv2.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print rows,cols
342 548
In [19]: nrows = cv2.getOptimalDFTSize(rows)
In [20]: ncols = cv2.getOptimalDFTSize(cols)
In [21]: print nrows, ncols
360 576
看到了吧，數(shù)組的大小從（342， 548）變成了（360， 576）。現(xiàn)在我們?yōu)樗a(bǔ) 0，然后看看性能有沒(méi)有提升。你可以創(chuàng)建一個(gè)大的 0 數(shù)組，然后把我們的數(shù)據(jù)拷貝過(guò)去，或者使用 函數(shù) cv2.copyMakeBoder() 。
nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img
或者：
right = ncols - cols
bottom = nrows - rows
#just to avoid line breakup in PDF file
bordertype = cv2.BORDER_CONSTANT
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)
?現(xiàn)在我們看看 Numpy 的表現(xiàn)：
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 倍。我們?cè)倏纯?OpenCV 的表現(xiàn)：
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 倍，同時(shí)我們也會(huì)發(fā)現(xiàn) OpenCV 的速度是 Numpy 的 3 倍 。你也可以測(cè)試一下逆 FFT 的表現(xiàn)。
?
1.4 為什么拉普拉斯算子是高通濾波器？
在論壇中遇到了一個(gè)類(lèi)似的問(wèn)題。為什么拉普拉斯算子是高通濾波器？為什么 Sobel 是 HPF？等等。對(duì)于第一個(gè)問(wèn)題的答案我們以傅里葉變換的形式給出。我們一起來(lái)對(duì)不同的算子進(jìn)行傅里葉變換并分析它們：
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)
# x.T 為矩陣轉(zhuǎn)置
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 range(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()
?結(jié)果：
從圖像中我們就可以看出每一個(gè)算子允許通過(guò)那些信號(hào)。從這些信息中我們就可以知道那些是 HPF 那是 LPF。
?
更多資源：
- An Intuitive Explanation of Fourier Theory?by Steven Lehar
- Fourier Transform?at HIPR
- What does frequency domain denote in case of images?
另：為了更清楚的觀察振鈴效應(yīng)，換用OpenCV的logo圖像做測(cè)試，如下：?
"""
使用 numpy 進(jìn)行 FFT 和 IFFT
"""
import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv2.imread('opencv_logo.jpg', 0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20 * np.log(np.abs(fshift))
rows, cols = img.shape
crow, ccol = int(rows / 2), int(cols / 2)
fshift[crow - 20: crow + 20, ccol - 20: ccol + 20] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
"""
使用 OpenCV 進(jìn)行 FFT 和 IFFT
"""
# cv2.dft( )函數(shù)返回值是雙通道的，第一個(gè)通道是實(shí)數(shù)部分，第二個(gè)是虛數(shù)部分
# cv2.dft( )函數(shù)輸入圖像要先轉(zhuǎn)為 np.float32 格式
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
magnitude_spectrum2 = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
# 創(chuàng)建低通濾波器掩模、使用掩模濾波、IDFT
mask = np.zeros((rows, cols, 2), np.uint8)
mask[crow - 30:crow + 30, ccol - 30:ccol + 30] = 1
f_dft_shift = dft_shift * mask
f_idft_shift = np.fft.ifftshift(f_dft_shift)
img_back2 = cv2.idft(f_idft_shift)
img_back2 = cv2.magnitude(img_back2[:, :, 0], img_back2[:, :, 1])
# 顯示結(jié)果圖像
imgList = [img, magnitude_spectrum, img_back, img, magnitude_spectrum2, img_back2]
imgName = ['img', 'magnitude_spectrum', 'img_back', 'img', 'magnitude_spectrum2', 'img_back2']
for i in range(6):
plt.subplot(2, 3, i + 1), plt.imshow(imgList[i], cmap='gray'), plt.title(imgName[i])
plt.xticks([]), plt.yticks([])
plt.show()
結(jié)果如下：
從上圖 img_back2 可以很明顯看出 使用矩形窗口濾波產(chǎn)生的振鈴效應(yīng) 。?
?
參考資料：
1.?OpenCV-Python官方文檔
2.《OpenCV-Python 中文教程》（段力輝 譯）
更多文章、技術(shù)交流、商務(wù)合作、聯(lián)系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號(hào)聯(lián)系： 360901061
您的支持是博主寫(xiě)作最大的動(dòng)力，如果您喜歡我的文章，感覺(jué)我的文章對(duì)您有幫助，請(qǐng)用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點(diǎn)擊下面給點(diǎn)支持吧，站長(zhǎng)非常感激您！手機(jī)微信長(zhǎng)按不能支付解決辦法：請(qǐng)將微信支付二維碼保存到相冊(cè)，切換到微信，然后點(diǎn)擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對(duì)您有幫助就好】元
