OpenCV 24: 傅里叶变换

开启成长之旅！这是我参与「日新计划 2 月更文挑战」的第 26 天，点击查看活动详情

方针

在本节中，将学习

• 运用OpenCV查找图画的傅立叶改换
• 运用Numpy中可用的FFT函数
• 傅立叶改换的某些应用程序
• 函数：cv2.dft()，cv2.idft()等

理论

傅立叶改换用于剖析各种滤波器的频率特性。关于图画，运用2D离散傅里叶改换(DFT)查找频域。一种称为**快速傅立叶改换(FFT)**的快速算法用于DFT的核算。关于这些的详细信息能够在任何图画处理或信号处理教科书中找到。

关于正弦信号$x(t)=A\sin (2ft)$, 能够说f是信号的频率，假如选用其频域，则能够看到f的尖峰。假如对信号进行采样以构成离散信号，将取得相同的频域，可是在[ − ， ]或[ 0 , 2 ]范围内（关于N点DFT为[ 0 ， N ] ）是周期性的。能够将图画视为在两个方向上采样的信号。因而，在X和Y方向都进行傅立叶改换，能够得到图画的频率标明。

更直观地说，关于正弦信号，假如起伏在短时间内改动非常快，则能够说它是高频信号。假如改动缓慢，则为低频信号。能够将相同的想法扩展到图画。图画中的振幅在哪里急剧改动？在边际点或噪声。因而，能够说边际和噪声是图画中的高频内容。假如起伏没有太大改动，则它是低频重量。

现在，将学习怎么进行傅立叶改换。

Numpy中的傅里叶改换

首要，将看到怎么运用Numpy进行傅立叶改换。Numpy具有FFT软件包来执行此操作。np.fft.fft2()供给了频率转化，它将是一个杂乱的数组。

np.fft.fft2()

• 第一个参数是输入图画，即灰度图画。
• 第二个参数是可选的，决定输出数组的巨细。 假如它大于输入图画的巨细，则在核算FFT之前用零填充输入图画。假如小于输入图画，将裁切输入图画。假如未传递任何参数，则输出数组的巨细将与输入的巨细相同。

现在，一旦取得成果，零频率重量（DC重量）将坐落左上角。假如要使其居中，则需要在两个方向大将成果都移动N2\frac{N}{2}2N​ 。也能够通过函数np.fft.fftshift()完结。找到频率改换后，就能够找到起伏谱。

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。为此，只需用尺度为60×60的矩形窗口遮罩即可消除低频。然后，运用np.fft.ifftshift()应用反向移位，以使DC重量再次出现在左上角。然后运用np.ifft2()函数找到逆FFT。相同，成果将是一个复数。能够选用其绝对值。

# 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

# 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来查看它们的功能。

# 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()。

nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img

或许:

right = ncols - cols
bottom = nrows - rows
bordertype = cv2.BORDER_CONSTANT ＃只是为了避免PDF文件中的行中止
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)

现在，核算Numpy函数的DFT功能比较：

%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函数。

%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：

# 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

附加资源

• docs.opencv.org/4.1.2/de/db…
• 傅里叶改换的直观解说
• 傅里叶改换
• 图画中的频率域指什么？
• 频域低通滤波、频域高通滤波