OpenCV-Python官方文檔關于圖像傅里葉變換和反變換的教程網址：
https://docs.opencv.org/4.1.0/de/dbc/tutorial_py_fourier_transform.html?
目標
我們將要學習：
? 使用 OpenCV 對圖像進行傅里葉變換（DFT）：cv2.dft()，cv2.idft()
? 使用 Numpy 中 FFT（快速傅里葉變換）函數：
? 傅里葉變換的一些用處
? 我們將要學習的函數有： cv2.dft()，cv2.idft() 等
原理
把圖像想象成沿著兩個方向采集的信號。所以對圖像同時進行 X 方向和 Y 方向的傅里葉變換，我們就會得到這幅圖像的頻域表示（頻譜圖）。更直觀一點，對于一個正弦信號，如果它的幅值變化非常快，我們可以說它是高頻信號，如果變化非常慢，我們稱之為低頻信號。把這種想法應用到圖像中，圖像哪里的幅度變化非常大呢？邊界點或者噪聲。所以我們說 邊界和噪聲是圖像中的高頻分量 。
1.1 Numpy 中的傅里葉變換?
Numpy 中的 FFT 包 可以幫助我們實現快速傅里葉變換。 函數 np.fft.fft2() 可以對信號進行頻率轉換，輸出結果是一個復雜的組。本函數的第一個參數是輸入圖像，要求是灰度格式。第二個參數是可選的, 決定輸出數組的大小。輸出數組的大小和輸入圖像大小一樣。如果輸出結果比輸入圖像大，輸入圖像就需要在進行 FFT 前補0。如果輸出結果比輸入圖像小的話，輸入圖像就會被切割。?
現在我們得到了結果，頻率為 0 的部分（直流分量）在輸出圖像的左上角。如果想讓它（直流分量）在輸出圖像的中心，我們還需要將結果沿兩個方向平移 N/2 。 函數 np.fft.fftshift() 可以幫助我們實現這一步。（這樣更容易分析）。進行完頻率變換之后，我們就可以構建振幅譜了。?
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)
# 這里構建振幅圖的公式沒學過
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()
結果：
然后就可以在頻域對圖像進行操作了，比如進行高通濾波，再傅里葉反變換回空間域。比如我們可以使用一個60x60 的矩形窗口對圖像進行掩模操作從而去除低頻分量。然后再使用 函數np.fft.ifftshift() 進行逆平移操作，所以現在直流分量又回到左上角了，然后使用 函數 np.ifft2() 進行 FFT 逆變換。同樣又得到一堆復雜的數字，我們可以對他們取絕對值：
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)
# 取絕對值
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()
結果如下：
上圖的結果顯示高通濾波其實是一種邊界檢測操作。這就是我們在前面圖像梯度那一章看到的。同時我們還發現圖像中的大部分數據集中在頻譜圖的低頻區域。
如果你觀察仔細的話，尤其是最后一章 JET 顏色的圖像，你會看到一些不自然的東西（如我用紅色箭頭標出的區域）。看上圖那里有些條帶裝的結構，這被成為 振鈴效應 。這是由于我們使用矩形窗口做掩模造成的。這個掩模被轉換成正弦形狀時就會出現這個問題。所以 一般我們不使用矩形窗口濾波。最好的選擇是高斯窗口。
1.2 OpenCV 中的傅里葉變換
OpenCV 中相應的函數是 cv2.dft() 和 cv2.idft()。和前面輸出的結果一樣，但是是雙通道的。第一個通道是結果的實數部分，第二個通道是結果的虛數部分。輸入圖像要首先轉換成 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()
?注意： 你可以使用 函數 cv2.cartToPolar() ，它會同時返回幅度和相位。
現在我們來做逆 DFT。在前面的部分我們實現了一個 HPF（高通濾波），現在我們來做 LPF（低通濾波）將高頻部分去除。其實就是對圖像進行模糊操作。首先我們需要構建一個掩模，與低頻區域對應的地方設置為 1, 與高頻區域對應的地方設置為 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()
?結果：
注意： OpenCV 中的函數 cv2.dft() 和 cv2.idft() 要比 Numpy 快 。但是Numpy 函數更加用戶友好。關于性能的描述，請看下面的小節。
1.3 DFT 的性能優化
當數組的大小為某些值時 DFT 的性能會更好。當數組的大小是 2 的指數時 DFT 效率最高。當數組的大小是 2， 3， 5 的倍數時效率也會很高。所以如果你想提高代碼的運行效率時，你可以修改輸入圖像的大小（補 0）。對于OpenCV 你必須自己手動補0。但是 Numpy，你只需要指定 FFT 運算的大小，它會自動補 0。
那我們怎樣確定最佳大小呢？OpenCV 提供了一個 函數 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 548
In [19]: nrows = cv2.getOptimalDFTSize(rows)
In [20]: ncols = cv2.getOptimalDFTSize(cols)
In [21]: print nrows, ncols
360 576
看到了吧，數組的大小從（342， 548）變成了（360， 576）。現在我們為它補 0，然后看看性能有沒有提升。你可以創建一個大的 0 數組，然后把我們的數據拷貝過去，或者使用 函數 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)
?現在我們看看 Numpy 的表現：
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 倍。我們再看看 OpenCV 的表現：
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 的速度是 Numpy 的 3 倍 。你也可以測試一下逆 FFT 的表現。
?
1.4 為什么拉普拉斯算子是高通濾波器？
在論壇中遇到了一個類似的問題。為什么拉普拉斯算子是高通濾波器？為什么 Sobel 是 HPF？等等。對于第一個問題的答案我們以傅里葉變換的形式給出。我們一起來對不同的算子進行傅里葉變換并分析它們：
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 為矩陣轉置
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()
?結果：
從圖像中我們就可以看出每一個算子允許通過那些信號。從這些信息中我們就可以知道那些是 HPF 那是 LPF。
?
更多資源：
- An Intuitive Explanation of Fourier Theory?by Steven Lehar
- Fourier Transform?at HIPR
- What does frequency domain denote in case of images?
另：為了更清楚的觀察振鈴效應，換用OpenCV的logo圖像做測試，如下：?
"""
使用 numpy 進行 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 進行 FFT 和 IFFT
"""
# cv2.dft( )函數返回值是雙通道的，第一個通道是實數部分，第二個是虛數部分
# cv2.dft( )函數輸入圖像要先轉為 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]))
# 創建低通濾波器掩模、使用掩模濾波、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])
# 顯示結果圖像
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()
結果如下：
從上圖 img_back2 可以很明顯看出 使用矩形窗口濾波產生的振鈴效應 。?
?
參考資料：
1.?OpenCV-Python官方文檔
2.《OpenCV-Python 中文教程》（段力輝 譯）
更多文章、技術交流、商務合作、聯系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號聯系： 360901061
您的支持是博主寫作最大的動力，如果您喜歡我的文章，感覺我的文章對您有幫助，請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點擊下面給點支持吧，站長非常感激您！手機微信長按不能支付解決辦法：請將微信支付二維碼保存到相冊，切換到微信，然后點擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對您有幫助就好】元
