【翻译：OpenCV-Python教程】傅里叶变换

最新推荐文章于 2021-09-08 14:39:58 发布

宇翔1024

最新推荐文章于 2021-09-08 14:39:58 发布

阅读量385

点赞数

分类专栏： OpenCV 翻译文章标签： OpenCV 计算机视觉人工智能中文文档

OpenCV 同时被 2 个专栏收录

37 篇文章 8 订阅

订阅专栏

翻译

37 篇文章 0 订阅

订阅专栏

⚠️由于自己的拖延症，3.4.3翻到一半，OpenCV发布了4.0.0了正式版，所以接下来是按照4.0.0翻译的。

⚠️除了版本之外，其他还是照旧，Fourier Transform，附原文。

目标

在这一节，我们将学习：

用OpenCV求图像的傅里叶变换
利用Numpy中可用的FFT函数(FFT下面有解释)
傅里叶变换的一些应用
我们会遇到以下这些函数：cv.dft()，cv.idft() 等等

理论

利用傅立叶变换是用来分析了各种滤波器的频率特性。对于图像来说，二维离散傅里叶变换(2D Discrete Fourier Transform (DFT))用来找出频域。一个较快的算法叫做快速傅里叶转换(Fast Fourier Transform (FFT))用来计算DFT。关于这些的细节，可以在任何图像处理或者信息处理的教科书里找到。请看额外资源的那一部分。

对于一个正弦信号，x(t)=Asin(2πft)，我们可以认为f就是信号频率。如果我们取到它的频域，那我们可以看到在f处有一个峰值。如果信号被采样形成一个离散的信号，我们会取到相同的频域，但却是周期性的出现在区间[−π,π]或者[0,2π] (或者[0,N]，对于多点离散傅里叶变换来说)。你可以把一张图像看作是在两个方向上采样的信号。在X和Y方向上进行傅里叶变换，就可以描述出图像的频率。

更直观的讲，对于正弦信号，如果振幅在短时间内变化非常大，你可以说它是一个高频信号。如果它变化很慢，那它就是个低频信号。你可以把同样的想法延伸到图像上，图像上什么地方的振幅会激烈变化呢？在边缘，或者噪音点上啊。所以我们可以说边缘和噪音是一张图像中的高频内容。如果没有那么大的振幅变化，那它就属于低频分量。(一些链接已经被添加在了额外资源的部分，它通过一个示例直观的解释了频率转换)。

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

Numpy里的傅里叶变换

首先我们看看如何用Numpy求傅里叶变换。Numpy有一个FFT包来做这件事。np.fft.fft2()提供给我们频率转换的算法，返回一个复杂的数组。它的第一个参数是输入图像，图像得是灰度图。第二个参数是可选的，它决定了输出数组的大小。如果它大于输入图像，输入图像在做FFT算法之前，就要用零来填充(到和这个输出大小一样)。如果它小于输入图像，输入图像就会被裁剪。如果你不传这个参数，输出数组的大小就会和输入图像一致。

现在一旦你得到这个结果，零频分量(直流分量(译者附：DC components))将会出现在左上角，如果你想要把它们移动到中心，你需要同时在两个方向把结果移动2N。这可以通过函数 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色图模式(译者注：JET从蓝到红，中间经过青绿、黄和橙色。它是HSV色图的一个变异。)来看最后一张的话。你可以看到一些副产物(我已经用红色箭头标记了其中一个实例)。它显示了一些波纹状的结构，这被称为振铃效应。它是由我们用来遮蔽的矩形窗口(60X60那个)引起的。这个遮蔽层被转换成辛格函数的形状(译者注：即sinc(x)，频域中的矩形对应时域中的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，High Pass Filtering)这次我们将看到如何删除图像中的高频内容，即我们将对图像应用低通滤波器(LPF)。它实际上会模糊图像。为此，我们首先创建一个在低频下具有高值(1)的遮罩层，就是说，(这个遮罩层的样子应该是)我们在低频区域传入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的s次幂，3的s次幂，以及5的s次幂时效率也很不错。所以如果您担心代码的性能，可以在找到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倍。这也可以拿逆FFTT来做测试，把这个测试作为你的练习吧。

为什么拉普拉斯算子是高通滤波器？

类似的问题在论坛上常被提出。问题大概是，为什么拉普拉斯算子是个高通滤波器，为什么索贝尔是个高通滤波器？等等。(译者附，你可以往前翻图像梯度的那一章。) 论坛上对此给出的最高赞回答，就是用傅里叶变换的术语来解释的。只要对拉普拉斯算子做傅里叶变换就可以得到更大范围的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 xrange(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()

看结果：

从图像中，你可以看到每个算子块的频率区域，以及它经过的区域。从这些信息中，我们可以知道为什么每个内核是高通滤波器或者低通滤波器。

额外资源

练习

上篇：【翻译：OpenCV-Python教程】OpenCV里的直方图

下篇：【翻译：OpenCV-Python教程】模板匹配

宇翔1024

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
1
评论
【翻译：OpenCV-Python教程】傅里叶变换

⚠️由于自己的拖延症，3.4.3翻到一半，OpenCV发布了4.0.0了正式版，所以接下来是按照4.0.0翻译的。⚠️除了版本之外，其他还是照旧，Fourier Transform，附原文。目标在这一节，我们将学习：用OpenCV求图像的傅里叶变换利用Numpy中可用的FFT函数(FFT下面有解释) 傅里叶变换的一些应用我们会遇到以下这些函数：cv.dft()，c...
复制链接

扫一扫