OpenCV 24: 傅立葉變換
開啟掘金成長之旅！這是我參與「掘金日新計劃 · 2 月更文挑戰」的第 26 天，點選檢視活動詳情
目標
在本節中，將學習
* 使用OpenCV查詢影象的傅立葉變換
* 利用Numpy中可用的FFT函式
* 傅立葉變換的某些應用程式
* 函式：cv2.dft()，cv2.idft()等
理論
傅立葉變換用於分析各種濾波器的頻率特性。對於影象，使用2D離散傅立葉變換(DFT)查詢頻域。一種稱為快速傅立葉變換(FFT)的快速演算法用於DFT的計算。關於這些的詳細資訊可以在任何影象處理或訊號處理教科書中找到。
對於正弦訊號$x(t)=A\sin(2\pi ft)$, 可以說f是訊號的頻率，如果採用其頻域，則可以看到f的尖峰。如果對訊號進行取樣以形成離散訊號，將獲得相同的頻域，但是在[ − π ， π ]或[ 0 , 2 π ]範圍內（對於N點DFT為[ 0 ， N ]）是週期性的。可以將影象視為在兩個方向上取樣的訊號。因此，在X和Y方向都進行傅立葉變換，可以得到影象的頻率表示。
更直觀地說，對於正弦訊號，如果幅度在短時間內變化非常快，則可以說它是高頻訊號。如果變化緩慢，則為低頻訊號。可以將相同的想法擴充套件到影象。影象中的振幅在哪裡急劇變化？在邊緣點或噪聲。因此，可以說邊緣和噪聲是影象中的高頻內容。如果幅度沒有太大變化，則它是低頻分量。
現在，將學習如何進行傅立葉變換。
Numpy中的傅立葉變換
首先，將看到如何使用Numpy進行傅立葉變換。Numpy具有FFT軟體包來執行此操作。np.fft.fft2()提供了頻率轉換，它將是一個複雜的陣列。
np.fft.fft2() * 第一個引數是輸入影象，即灰度影象。 * 第二個引數是可選的，決定輸出陣列的大小。 如果它大於輸入影象的大小，則在計算FFT之前用零填充輸入影象。如果小於輸入影象，將裁切輸入影象。如果未傳遞任何引數，則輸出陣列的大小將與輸入的大小相同。
現在，一旦獲得結果，零頻率分量（DC分量）將位於左上角。如果要使其居中，則需要在兩個方向上將結果都移動$\frac{N}{2}$ 。也可以通過函式np.fft.fftshift()完成。找到頻率變換後，就可以找到幅度譜。
```python
import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv2.imread('messi.png', 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() ```
結果看起來像下面這樣:
可以看到，在中心看到更多白色區域，這表明低頻內容更多。
因此，發現了頻率變換，可以在頻域中進行一些操作，例如高通濾波和重建影象，即找到逆DFT。為此，只需用尺寸為60x60的矩形視窗遮罩即可消除低頻。然後，使用np.fft.ifftshift()應用反向移位，以使DC分量再次出現在左上角。然後使用np.ifft2()函式找到逆FFT。同樣，結果將是一個複數。可以採用其絕對值。
```python
ifft
import cv2 import numpy as np from matplotlib import pyplot as plt
img = cv2.imread('messi.png', 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+31, ccol-30:ccol+31] = 0 f_ishift = np.fft.ifftshift(fshift) img_back = np.fft.ifft2(f_ishift) img_back = np.real(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, cmap='gray') plt.title('Result in JET') plt.xticks([]) plt.yticks([])
plt.show() ```
結果看起來像下面這樣：
結果表明高通濾波是邊緣檢測操作。這就是在“影象漸變”一章中看到的。這也表明大多數影象資料都存在於頻譜的低頻區域。 如果仔細觀察結果，尤其是最後一張JET顏色的影象，會看到一些偽像（用紅色箭頭標記的一個例項）。它在那裡顯示出一些波紋狀結構，稱為振鈴效應(ringings effects)。這是由用於遮罩的矩形視窗引起的。此掩碼轉換為正弦形狀，從而導致此問題。因此，矩形視窗不用於過濾。更好的選擇是高斯視窗。
OpenCV中的傅立葉變換
OpenCV為此提供了cv2.dft()和cv2.idft()函式。它返回與前一個相同的結果，但是有兩個通道。
* 第一個通道是結果的實部
* 第二個通道是結果的虛部。
輸入影象首先應轉換為np.float32
```python
opencv
import cv2 import numpy as np from matplotlib import pyplot as plt
img = cv2.imread('messi.png', 0) img_32 = np.float32(img) dft = cv2.dft(img_32, 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 = rows//2, cols//2 print(crow, ccol)
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的冪時，速度最快。對於大小為2、3和5的乘積的陣列，也可以非常有效地進行處理。因此，如果擔心程式碼的效能，可以在找到DFT之前將陣列的大小修改為任何最佳大小(通過填充零)。對於OpenCV而言，必須手動填充零。但是對於Numpy，指定FFT計算的新大小，它將自動填充零。
那麼如何找到最優的大小呢?OpenCV為此提供了一個函式，cv2.getOptimalDFTSize()。它同時適用於cv2.dft()和np.fft.fft2()。使用IPython魔術命令timeit來檢查它們的效能。
```python
performance
img = cv2.imread('messi.png', 0) rows, cols = img.shape print("{}, {}".format(rows, cols))
259, 419
nrows = cv2.getOptimalDFTSize(rows) ncols = cv2.getOptimalDFTSize(cols) print("{}, {}".format(nrows, ncols))
270, 432
可以看到，將大小(259, 419)修改為(270，432)。現在用零填充（對於OpenCV），並找到其DFT計算效能。可以通過建立一個新的零陣列並將資料複製到其中來完成此操作，或者使用`cv2.copyMakeBorder()`。python
nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img
```
或者:
python
right = ncols - cols
bottom = nrows - rows
bordertype = cv2.BORDER_CONSTANT ＃只是為了避免PDF檔案中的行中斷
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)
現在，計算Numpy函式的DFT效能比較： ```python %timeit fft1 = np.fft.fft2(img)
15.4 ms ± 1.22 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit fft2 = np.fft.fft2(img,[nrows,ncols])
7.55 ms ± 742 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```
通過用0填充改變尺寸，效能有了2倍的加速。
現在將嘗試使用OpenCV函式。 ```python %timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
3.21 ms ± 217 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
989 µs ± 91.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
```
可以看到，效能有了4倍的加速。還可以看到OpenCV函式比Numpy函式快5倍左右。
為什麼拉普拉斯運算元是高通濾波器？
為什麼拉普拉斯變換是高通濾波器? 為什麼Sobel是HPF?。第一個答案是關於傅立葉變換的。只需採取Laplacian的傅立葉變換，以獲得更高尺寸的FFT： ```python
laplacian is high pass filter
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 gaussian filter
x = cv2.getGaussianKernel(5, 10) gaussian = x*x.T
different edge detecting filters
scharr in x-direction
scharr_x = 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 directio
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_x] 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 = [20*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