OpenCV-Python中的圖像處理-傅里葉變換
- 傅里葉變換
- Numpy中的傅里葉變換
- Numpy中的傅里葉逆變換
- OpenCV中的傅里葉變換
- OpenCV中的傅里葉逆變換
- DFT的性能優化
- 不同濾波算子傅里葉變換對比
傅里葉變換
- 傅里葉變換經常被用來分析不同濾波器的頻率特性。我們可以使用 2D 離散傅里葉變換 (DFT) 分析圖像的頻域特性。實現 DFT 的一個快速算法被稱為快速傅里葉變換（ FFT）。
- 對于一個正弦信號：x (t) = A sin (2πft), 它的頻率為 f，如果把這個信號轉到它的頻域表示，我們會在頻率 f 中看到一個峰值。如果我們的信號是由采樣產生的離散信號組成，我們會得到類似的頻譜圖，只不過前面是連續的，現在是離散。你可以把圖像想象成沿著兩個方向采集的信號。所以對圖像同時進行 X 方向和 Y 方向的傅里葉變換，我們就會得到這幅圖像的頻域表示（頻譜圖）。
- 對于一個正弦信號，如果它的幅度變化非常快，我們可以說他是高頻信號，如果變化非常慢，我們稱之為低頻信號。你可以把這種想法應用到圖像中，圖像那里的幅度變化非常大呢？邊界點或者噪聲。所以我們說邊界和噪聲是圖像中的高頻分量（注意這里的高頻是指變化非常快，而非出現的次數多）。如果沒有如此大的幅度變化我們稱之為低頻分量。
Numpy中的傅里葉變換
Numpy 中的 FFT 包可以幫助我們實現快速傅里葉變換。函數 np.fft.fft2() 可以對信號進行頻率轉換，輸出結果是一個復雜的數組。本函數的第一個參數是輸入圖像，要求是灰度格式。第二個參數是可選的, 決定輸出數組的大小。輸出數組的大小和輸入圖像大小一樣。如果輸出結果比輸入圖像大，輸入圖像就需要在進行 FFT 前補0。如果輸出結果比輸入圖像小的話，輸入圖像就會被切割。
f = np.fft.fft2(img)
現在我們得到了結果，頻率為 0 的部分（直流分量）在輸出圖像的左上角。如果想讓它（直流分量）在輸出圖像的中心，我們還需要將結果沿兩個方向平移 N/2 。函數 np.fft.fftshift() 可以幫助我們實現這一步。（這樣更容易分析）。進行完頻率變換之后，我們就可以構建振幅譜了。
fshift = np.fft.fftshift(f)
import numpy as np
import cv2
from matplotlib import pyplot as pltimg = cv2.imread('./resource/opencv/image/messi5.jpg', cv2.IMREAD_GRAYSCALE)
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.subplot(122), plt.imshow(magnitude_spectrum, cmap='gray'), plt.title('Magnitude Spectrum')
plt.show()
我們可以看到輸出結果的中心部分更白（亮），這說明低頻分量更多。
Numpy中的傅里葉逆變換
- 對圖像進行FFT變換之后得到頻域圖像數據，然后再做IFFT變換又可以得到原始圖像。相關函數：np.fft.ifftshift(),np.fft.ifft2()
fishift = np.fft.ifftshift(fshift)
img_ifft = np.fft.ifft2(fishift) - 我們可以對頻域圖像數據進行操作以實現一些圖像處理效果，如在頻域內將低頻分量的值設為0，可以實現對圖像的高通濾波處理：
rows, cols = img.shape
crow, ccol = int(rows/2) , int(cols/2)
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
import numpy as np
import cv2
from matplotlib import pyplot as pltimg = cv2.imread('./resource/opencv/image/messi5.jpg', cv2.IMREAD_GRAYSCALE)# 1.在Numpy內對圖像進行傅里葉變換，得到其頻域圖像
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# 這里是構建振幅圖，顯示圖像頻譜
magnitude_spectrum = 20*np.log(np.abs(fshift))# 2.IFFT 將頻域圖像還原成原始圖像，這里只是驗證FFT的逆運算
fishift = np.fft.ifftshift(fshift)
img_ifft = np.fft.ifft2(fishift)
img_ifft = np.abs(img_ifft) # 取絕對值，否則不能用imshow()來顯示圖像# 3.在頻域內將低頻分量的值設為0，實現高通濾波。
rows, cols = img.shape
crow, ccol = int(rows/2) , int(cols/2)
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0 # 4.對高通濾波后的圖像頻域數據進行逆傅里葉變換，得到高通濾波后圖像。
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back) # 取絕對值，否則不能用imshow()來顯示圖像
# 構建高通濾波后的振幅圖，顯示圖像頻譜
after_sepctrum = 20*np.log(np.abs(fshift))plt.subplot(231), plt.imshow(img, cmap='gray'), plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(232), plt.imshow(magnitude_spectrum, cmap='gray'), plt.title('Input Image Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(233), plt.imshow(img_ifft, cmap='gray'), plt.title('Input IFFT'), plt.xticks([]), plt.yticks([])
plt.subplot(234), plt.imshow(after_sepctrum, cmap='gray'), plt.title('After HPF Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(235), plt.imshow(img_back, cmap='gray'), plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(236), plt.imshow(img_back), plt.title('Result in JET'), plt.xticks([]), plt.yticks([])
plt.show()
OpenCV中的傅里葉變換
OpenCV 中相應的函數是 cv2.dft() 和 cv2.idft()。和前面輸出的結果一樣，但是是雙通道的。第一個通道是結果的實數部分，第二個通道是結果的虛數部分。輸入圖像要首先轉換成 np.float32 格式。
import numpy as np
import cv2
from matplotlib import pyplot as pltimg = cv2.imread('./resource/opencv/image/messi5.jpg', cv2.IMREAD_GRAYSCALE)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()
OpenCV中的傅里葉逆變換
前面的部分我們實現了一個 HPF（高通濾波）高通濾波其實是一種邊界檢測操作。現在我們來做 LPF（低通濾波）將高頻部分去除。其實就是對圖像進行模糊操作。首先我們需要構建一個掩模，與低頻區域對應的地方設置為 1, 與高頻區域對應的地方設置為 0。
import numpy as np
import cv2
from matplotlib import pyplot as pltimg = cv2.imread('./resource/opencv/image/messi5.jpg', cv2.IMREAD_GRAYSCALE)
# 1.OpenCV中做DFT
dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)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('Output Image'), plt.xticks([]), plt.yticks([])
plt.show()
DFT的性能優化
- 當數組的大小為某些值時 DFT 的性能會更好。當數組的大小是 2 的指數時 DFT 效率最高。當數組的大小是 2， 3， 5 的倍數時效率也會很高。所以如果你想提高代碼的運行效率時，你可以修改輸入圖像的大小（補 0）。對于OpenCV 你必須自己手動補 0。但是 Numpy，你只需要指定 FFT 運算的大小，它會自動補 0。
- OpenCV 提供了一個函數:cv2.getOptimalDFTSize()來確定最佳大小。它可以同時被 cv2.dft() 和 np.fft.fft2() 使用。
import numpy as np
import cv2img = cv2.imread('./resource/opencv/image/messi5.jpg', cv2.IMREAD_GRAYSCALE)
rows,cols = img.shape
print('原始圖像大小:',rows, cols)
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
print('DFT最佳大小:',nrows, ncols)
原始圖像大小: 342 548
DFT最佳大小: 360 576
import numpy as np
import cv2
from matplotlib import pyplot as pltimg = cv2.imread('./resource/opencv/image/messi5.jpg', cv2.IMREAD_GRAYSCALE)
rows,cols = img.shape
print('原始圖像大小:',rows, cols)
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
print('DFT最佳大小:',nrows, ncols)# Numpy數組操作，原圖擴大到最佳DFT size
nimg = np.zeros((nrows, ncols))
nimg [:rows, :cols] = img#
right = ncols - cols
bottom = nrows - rows
# just to avoid line breakup in PDF file
mimg = cv2.copyMakeBorder(img, 0, bottom, 0, right, cv2.BORDER_CONSTANT, value=0)plt.subplot(231), plt.imshow(img, cmap='gray')
plt.subplot(232), plt.imshow(nimg, cmap='gray')
plt.subplot(233), plt.imshow(mimg, cmap='gray')
plt.show()
不同濾波算子傅里葉變換對比
為什么拉普拉斯算子是高通濾波器？為什么 Sobel 是 HPF？等等。對于第一個問題的答案我們以傅里葉變換的形式給出。我們一起來對不同的算子進行傅里葉變換并分析它們：
import numpy as np
import cv2
from matplotlib import pyplot as plt# simple averaging filter whitout 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()
