opencv(4.5.3)-python(二十七)--傅里叶变换

用户9875047

发布于 2023-02-26 15:14:49

7750

发布于 2023-02-26 15:14:49

文章被收录于专栏：机器视觉全栈er

翻译及二次校对：cvtutorials.com

目标

在本节中，我们将学习:

• 使用OpenCV找到图像的傅里叶变换
• 利用Numpy中的FFT函数
• 傅立叶变换的一些应用
• 我们将看到以下函数：cv.dft(), cv.idft()等。

理论

傅里叶变换被用来分析各种过滤器的频率特性。对于图像，二维离散傅里叶变换（DFT）被用来寻找频域。一种叫做快速傅里叶变换（FFT）的快速算法被用来计算DFT。关于这些的细节可以在任何图像处理或信号处理教科书中找到。请看其他资源部分。

对于一个正弦信号，x(t)=Asin(2πft)，我们可以说f是信号的频率，如果取其频域，我们可以在f处看到一个尖峰。如果信号被采样形成离散信号，我们得到相同的频域，但在[-π,π]或[0,2π]范围内是周期性的（或者对于N点DFT来说是[0,N]）。你可以把图像看作是一个在两个方向上被采样的信号。因此，在X和Y两个方向上进行傅里叶变换，就可以得到图像的频率表示。

更直观地说，对于正弦信号，如果振幅在短时间内变化得很快，你可以说它是一个高频信号。如果它变化缓慢，它就是一个低频信号。你可以把同样的想法延伸到图像上。在图像中，哪里的振幅变化剧烈？在边缘点，或噪音。所以我们可以说，边缘和噪音是图像中的高频内容。如果振幅没有太大的变化，那就是低频成分。(一些链接被添加到附加资源中，它用例子直观地解释了频率变换)。

现在我们来看看如何找到傅里叶变换。

Numpy中的傅里叶变换

首先我们将看到如何使用Numpy找到傅立叶变换。np.fft.ft2()为我们提供了频率变换，它将是一个复数。它的第一个参数是输入图像，它是灰度的。第二个参数是可选的，决定输出数组的大小。如果它大于输入图像的大小，在计算FFT之前，输入图像将被填充零。如果它小于输入图像，输入图像将被裁剪。如果没有传递参数，输出数组的大小将与输入相同。

现在一旦你得到结果，零频率分量（直流分量）将在左上角。如果你想把它带到中心，你需要把结果在两个方向上移位N/2。这可以通过函数np.fft.fftshift()简单地完成。(它更容易分析）。一旦你找到了频率变换，你就可以找到幅值谱了。

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.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()

结果如下:

看，你可以看到中心的白色区域更多，显示出低频内容更多。

所以你找到了频率变换现在你可以在频域做一些操作，比如高通滤波和重建图像，即找到反DFT。为此，你只需用一个大小为60x60的矩形窗口进行遮蔽，以去除低频。然后用np.fft.ifftshift()进行反移位，使直流成分再次出现在左上角。然后使用np.ifft2()函数找到反FFT。其结果也是一个复数。你可以取其绝对值。

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)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])
plt.show()

结果如下：

结果显示高通滤波是一种边缘检测操作。这就是我们在图像梯度一章中看到的情况。这也表明大部分的图像数据存在于频谱的低频区域。总之我们已经看到了如何在Numpy中找到DFT、IDFT等。现在让我们看看如何在OpenCV中实现。

如果你仔细观察结果，特别是最后一张JET颜色的图像，你可以看到一些伪影（其中一个例子我已经用红色箭头标出）。它显示了一些类似波纹的结构，这被称为振铃效应。这是由我们用于遮蔽的矩形窗口造成的。这个掩膜被转换为sinc形状，导致了这个问题。所以矩形窗口不能用于滤波。更好的选择是高斯窗口。

OpenCV中的傅立叶变换

OpenCV为此提供了cv.dft()和cv.idft()函数。它返回的结果与前面的相同，但有两个通道。第一个通道将有结果的实部，第二个通道将有结果的虚部。输入的图像应该首先被转换为np.float32。我们将看到如何做到这一点。

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
img = cv.imread('messi5.jpg',0)
dft = cv.dft(np.float32(img),flags = cv.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)
magnitude_spectrum = 20*np.log(cv.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()

注意：你也可以使用cv.cartToPolar()，它可以一次性返回幅值和相位。

所以，现在我们要做反DFT。在上一节课中，我们创建了一个HPF，这次我们将看到如何去除图像中的高频内容，即我们对图像应用LPF。它实际上模糊了图像。为此，我们首先创建一个掩膜，在低频处用高值（1），即我们通过低频内容，而在高频区域用0。

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 = cv.idft(f_ishift)
img_back = cv.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函数cv.dft()和cv.idft()比Numpy的对应函数要快。但Numpy函数更方便用户使用。关于性能问题的更多细节，请看下面的章节。

DFT的性能优化

DFT计算的性能对于某些数组大小来说是比较好的。当数组大小为2的幂时，它是最快的。对于那些大小为2、3、5的乘积的数组，处理起来也相当有效。因此，如果你担心你的代码的性能，你可以在寻找DFT之前将数组的大小修改为任何最佳大小（通过填充零）。对于OpenCV，你必须手动填充零。但是对于Numpy来说，你指定FFT计算的新大小，它就会自动为你填充零。

那么我们如何找到这个最佳尺寸呢？OpenCV为此提供了一个函数，cv.getOptimalDFTSize()。它同时适用于cv.dft()和np.fft.fft2()。让我们用IPython的神奇命令timeit来检查它们的性能。

In [16]: img = cv.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print("{} {}".format(rows,cols))
342 548
In [19]: nrows = cv.getOptimalDFTSize(rows)
In [20]: ncols = cv.getOptimalDFTSize(cols)
In [21]: print("{} {}".format(nrows,ncols))
360 576

看，大小（342,548）被修改为（360, 576）。现在让我们用零来填充它（对于OpenCV来说），并找到它们的DFT计算性能。你可以通过创建一个新的零数组并将数据复制到其中，或者使用cv.copyMakeBorder()来完成。

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

或者：

right = ncols - cols
bottom = nrows - rows
bordertype = cv.BORDER_CONSTANT #just to avoid line breakup in PDF file
nimg = cv.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)

现在我们计算一下Numpy函数的DFT性能比较。

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= cv.dft(np.float32(img),flags=cv.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv.dft(np.float32(nimg),flags=cv.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop

它还显示了4倍的速度。你还可以看到OpenCV函数比Numpy函数快3倍左右。这也可以用于反FFT的测试，而这也是留给你的一个练习。

为什么Laplacian是一个高通滤波器？

在一个论坛上也有一个类似的问题。问题是，为什么Laplacian是高通滤波器？为什么Sobel是高通滤波器？而给出的第一个答案是傅里叶变换。就拿拉普拉斯的傅里叶变换来说吧，它的FFT大小较高。对它进行分析。

import cv2 as cv
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 = cv.getGaussianKernel(5,10)
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()

结果：