Opencv筆記（二十一）——傅里葉變換

參考

Numpy 中的傅里葉變換

    首先我們看看如何使用 Numpy 進行傅里葉變換。Numpy 中的 FFT 包可以幫助我們實現快速傅里葉變換。函數 np.fft.fft2() 可以對信號進行頻率轉換,輸出結果是一個復雜的數組。本函數的第一個參數是輸入圖像,要求是灰度格式。第二個參數是可選的, 決定輸出數組的大小。輸出數組的大小和輸入圖像大小一樣。如果輸出結果比輸入圖像大,輸入圖像就需要在進行 FFT 前補0。如果輸出結果比輸入圖像小的話,輸入圖像就會被切割。

   頻率為0 的部分(直流分量)在輸出圖像的左上角。(2D傅里葉變換F(x,y)的F(0,0)位置在圖像的左上角,F(0,0)表示的是圖像灰度的均值)如果想讓它(直流分量)在輸出圖像的中心,我們還需要將結果沿兩個方向平移 N/2 。函數 np.fft.fftshift() 可以幫助我們實現這一步。

# coding=utf-8
import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread("/home/wl/3.jpg", 0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# 這里構建振幅圖的公式沒學過
magnitude_spectrum = 20*np.log(np.abs(fshift))#先取絕對值,表示取模。取對數,將數據范圍變小 
print magnitude_spectrum
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()

我們可以看到輸出結果的中心部分更白(亮),這說明低頻分量更多。現在我們可以進行頻域變換了,我們就可以在頻域對圖像進行一些操作了,例如高通濾波和重建圖像(DFT 的逆變換)。比如我們可以使用一個60x60 的矩形窗口對圖像進行掩模操作從而去除低頻分量。然后再使用函數np.fft.ifftshift() 進行逆平移操作,所以現在直流分量又回到左上角了,左后使用函數 np.ifft2() 進行 FFT 逆變換。同樣又得到一堆復雜的數字,我們可以對他們取絕對值:

# coding=utf-8
import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread("/home/wl/3.jpg", 0)
f = np.fft.fft2(img)#得到結果為復數矩陣
fshift = np.fft.fftshift(f)#直接取中心
rows, cols = img.shape
crow,ccol = rows/2 , cols/2
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0#蒙板大小60×60
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)#使用FFT逆變換,此時結果仍然是復數
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()

上圖的結果顯示高通濾波其實是一種邊界檢測操作。這就是我們在前面圖像梯度那一章看到的。同時我們還發現圖像中的大部分數據集中在頻譜圖的低頻區域。

OpenCV 中的傅里葉變換

OpenCV 中相應的函數是 cv2.dft() 和 cv2.idft()。和前面輸出的結果一樣,但是是雙通道的。第一個通道是結果的實數部分,第二個通道是結果的虛數部分。輸入圖像要首先轉換成 np.float32 格式。

# coding=utf-8
import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread("/home/wl/3.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。

# coding=utf-8
import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread("/home/wl/3.jpg", 0)
dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
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 中的函數 cv2.dft() 和 cv2.idft() 要比 Numpy 快。但是Numpy 函數更加用戶友好。

DFT 的性能優化

    當數組的大小為某些值時 DFT 的性能會更好。當數組的大小是 2 的指數時 DFT 效率最高。當數組的大小是 2,3,5 的倍數時效率也會很高。所以如果你想提高代碼的運行效率時,你可以修改輸入圖像的大小(補 0)。對於OpenCV 你必須自己手動補 0。但是 Numpy,你只需要指定 FFT 運算的大小,它會自動補 0。那我們怎樣確定最佳大小呢?OpenCV 提供了一個函數:cv2.getOptimalDFTSize()。它可以同時被 cv2.dft() 和 np.fft.fft2() 使用。

# coding=utf-8
import cv2
import numpy as np

img = cv2.imread("/home/wl/3.jpg", 0)
dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
rows, cols = img.shape
print rows,cols
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
print nrows,ncols

1420 946
1440 960