注释:本文翻译自OpenCV3.0.0 document->OpenCV-Python Tutorials,包括对原文档种错误代码的纠正
3.10 OpenCV中的图像变换
第一节:傅里叶变换(Fourier Transform)
1.目标
- 使用OpenCV查找图像的傅里叶变换
- 利用Numpy中可用的FFT函数
- 傅里叶变换的一些应用
- 我们将学习以下函数:cv2.dft()、cv2.idft()等
2.原理
傅里叶变换用于分析各种滤波器的频率特性。对于图像,使用2D离散傅里叶变换(DFT)来查找频域。称为快速傅里叶变换(FFT)的快速算法用于计算DFT。有关这些细节可以在任何图像处理或信号处理教科书中找到。请参考其它资源部份。
对于一个正弦信号, ,我们可以说f是信号的频率,如果取其频率域,我们可以在f处看到一个尖峰。如果信号被采样以形成离散信号,我们得到相同的频域。但是在范围 或 或 ,你可以将图像视为在两个方向上采样的信号。因此,在X方向和Y方向进行傅里叶变换将为你提供图像的频率表示。
更直观地说,对于正弦信号,如果幅度在短时间内变化得如此之快,则可以说它是高频信号。如果变化缓慢,则是低频信号。你可以将相同的想法扩展到图像。图像中振幅的变化幅度如何?在边缘点,或噪音?所以我们可以说,边缘和噪音是图像中的高频内容,如果振幅没有太大的变化,则它是低频分量。(一些链接被添加到附加资源中,通过示例直观地解释频率转换)。
现在我们将看到如何找到傅里叶变换。
3.Numpy中的傅里叶变换
首先,我们将看到如何使用Numpy查找傅里叶变换。Numpy有一个FFT软件包来做到这一点。np.fft.fft2()为我们提供了频率变换,这将是一个复杂的数组。它的第一个参数是灰度的输入图像,第二个参数是可选的,它决定输出数组的大小。如果它大于输入图像的大小,则在计算FFT之前,输入图像用0填充。如果它小于输入图像,输入图像将被裁剪。如果没有参数传递,则输出数组大小将与输入相同。
现在一旦得到结果,零频率分量(直流分量)将位于左上角。如果要将它置于中心,则需要在两个方向上将结果移动 。这只是通过np.fft.fft2()完成的(分析起来更容易)。一旦你找到频率变换,就可以找到幅度谱。
示例如下:
import cv2 import numpy as np from matplotlib import pyplot as plt img = cv2.imread('5.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, 'gray') plt.title('Input Image'), plt.xticks([]), plt.yticks([]) plt.subplot(122), plt.imshow(magnitude_spectrum, 'gray') plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([]) plt.show()
结果:
看,你可以在中心看到更多更白的区域,显示低频内容更多。
所以你找到了频率变换,现在你可以在频域中做一些操作,如高通滤波和重构图像,即找到逆DFT。为此,你只需通过使用尺寸为60*60的矩形窗口进行遮盖来移除低频。然后使用np.fft.ifftshit()应用反向偏移,以便DC分量再次出现在左上角。然后使用np.ifft2()函数来找到反FFT。结果也是一个复杂的数字。你可以拿它的绝对值。
import cv2 import numpy as np from matplotlib import pyplot as plt img = cv2.imread('5.jpg', 0) f = np.fft.fft2(img) # 快速傅里叶变换算法得到频率分布 fshift = np.fft.fftshift(f) # 默认结果中心点位置是在左上角,转移到中间位置 magnitude_spectrum = 20 * np.log(np.abs(fshift)) # 结果是复数,求绝对值才是振幅 # 找到频域变换后,可以在频域中做一些操作,如高通滤波和重构图像,即找到逆DFT rows, cols = img.shape crow, ccol = rows / 2, cols / 2 fshift[crow - 30:crow + 30, ccol - 30:ccol + 30] = 0 f_ishift = np.fft.ifftshift(fshift) img_back = np.fft.ifft2(f_ishift) img_back = np.abs(img_back) plt.subplot(131), plt.imshow(img, 'gray') plt.title("Input Image"), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(img_back, '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()
结果:
结果显示高通滤波是边缘检测操作。这就是我们在Image Gradients章节看到的。这也表明大部分图像数据存在于频谱的低频区域。无论如何,我们已经看到如何在Numpy中找到DFT,IDFT等。现在让我们看看如何在OpenCV中完成它。
如果仔细观察结果,尤其是JET颜色中的最后一幅图像,可以看到一些工作,它显示了一些波纹状结构,它被称为震荡效应。这是由于我们用于掩蔽的矩形窗口引起的。这个蒙版被转换为正弦形状,导致了问题。所以矩形窗口不用于过滤。更好的选择是高斯窗口。
4.OpenCV中的傅里叶变换
OpenCV为此提供了cv2.dft()和cv2.idft().它返回与以前相同的结果,但具有两个通道。第一个通道将具有结果的实际部分,第二个通道将具有结果的虚部。输入图像应首先转换为np.float32.我们将看到如何去做。
注意:你也可以使用cv2.cartToPolar(),它可以一次性返回幅值和相位。
所以,现在我们必须做逆DFT。在之前的会话中,我们创建了一个HPF,这次我们将看到如何去除图像中的高频内容,即我们将LPF应用于图像。它实际上模糊了图像,为此,我们首先在低频处创建一个高值(1)的掩膜,即我们通过LF内容,在HF区域传递0.
# -*- coding: utf-8 -*- ''' OpenCV中的傅里叶变换: 1.cv2.dft()和cv2.idft(),输入图像应先转换成np.float32 attention:也可以使用cv2.cartToPolar(),可以一次性返回幅值和相位 cv2.dft()和cv2.idft()比Numpy对应的函数更快,但是Numpy函数更友好。 ''' import cv2 import numpy as np from matplotlib import pyplot as plt img = cv2.imread('5.jpg', 0) 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, 'gray') plt.title('Input Image'), plt.xticks([]), plt.yticks([]) plt.subplot(122), plt.imshow(magnitude_spectrum, 'gray') plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([]) plt.show() ''' rows, cols = img.shape crow, ccol = int(rows / 2), int(cols / 2) # 这里还是必须转换为整数哦 # 首先创建一个掩膜,中心点为1,其余值为0 mask = np.zeros((rows, cols, 2), np.uint8) mask[crow - 30:crow + 30, ccol - 30: ccol + 30] = 1 # 应用掩膜和逆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, 'gray') plt.title('Input Image'), plt.xticks([]), plt.yticks([]) plt.subplot(122), plt.imshow(img_back, 'gray') plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([]) plt.show()
结果:
5.DFT的性能优化
对于某些数组大小,DFT计算的性能更好。当阵列大小是2的幂时,它是最快的。大小为2,3和5的乘积的数组也被相当有效地处理。因此,如果你担心代码的性能,可以在查找DFT之前将数组的大小修改为任何最佳大小(通过填充零)。对于OpenCV,你必须手动填充零。但对于Numpy,你指定了FFT计算的新大小,并且它会自动为你填充零。
那么我们如何寻找这个最佳尺寸?OpenCV为此专门提供了一个函数cv2.getOptimalDFTSize().它适用于cv2.dft()和np.fft.fft2()。让我们使用Python的命令好%timeit来检查它们的性能。
In [16]: img = cv2.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print rows,cols
342 548
In [19]: nrows = cv2.getOptimalDFTSize(rows)
In [20]: ncols = cv2.getOptimalDFTSize(cols)
In [21]: print nrows, ncols
360 576
看,(342,548)被修改为(360,576).现在让我们用零来填充(对于OpenCV)并且找到它们的DFT计算性能。你可以通过创建一个新的空数组并将数据复制给它,或者使用cv2.copyMakeBorder()来完成。
nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img
或者:
right = ncols - cols
bottom = nrows - rows
bordertype = cv2.BORDER_CONSTANT #just to avoid line breakup in PDF file
nimg = cv2.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
它显示了四倍的加速,现在我们尝试OpenCV函数。
In [24]: %timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop
它也显示了四倍加速。你还可以看到OpenCV函数比Numpy函数大约快倍。这也可以针对逆FFT进行测试。对你来说是一个练习。
6.为什么拉普拉斯算子是高通滤波器
在论坛上提出了类似的问题。问题是,为什么拉普拉斯是高通滤波器?为什么Sobel是HPF,等等。第一个答案就是傅里叶变换。只需对拉普拉斯算子进行傅里叶变换即可获得更高的FFT大小。分析它:
# -*- coding: utf-8 -*- ''' 集中不同滤波器的对比:均值滤波,高斯滤波、scharr滤波、sobel_x,sobel_y,拉普拉斯算子 ''' import cv2 import numpy as np from matplotlib import pyplot as plt #简单的均值滤波器没有参数 mean_filter=np.ones((3,3)) #创建一个高斯滤波器 x=cv2.getGaussianKernel(5,10) gaussian=x*x.T #不同的边缘检测滤波器 #x方向上的scharr scharr=np.array([[-3,0,3], [-10,0,10], [-3,0,3]]) #x方向上的sobel sobel_x=np.array([[-1,0,1], [-2,0,2], [-1,0,1]]) #y方向上的sobel sobel_y=np.array([[-1,-2,-1], [0,0,0], [1,2,1]]) #拉普拉斯算子 laplacian=np.array([[0,1,0], [1,-4,1], [0,1,0]]) filters=[mean_filter,gaussian,scharr,sobel_x,sobel_y,laplacian] filter_name=['mean_filter','gaussian','scharr','sobel_x','sobel_y','laplacian'] fft_filter=[np.fft.fft2(x) for x in filters] fft_shift=[np.fft.fftshift(y) for y in fft_filter] #地图谱 map_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(map_spectrum[i],'gray') plt.title(filter_name[i]),plt.xticks([]),plt.yticks([]) plt.show()
结果:
从图像中,可以看到每个内核阻塞的频率区域以及它传递的区域。从这些信息中,我们可以说出为什么每个内核都是HPF或LPF。
3.11 模板匹配
1.目标
- 使用模板匹配查找图像中的对象
- 学会这些函数:cv2.matchTemplate()、cv2.minMaxLoc()
2.原理
模板匹配是一种在较大图像中搜索和查找模板图像位置的方法。OpenCV为此提供了函数cv2.matchTemplate().它只是将模板图像滑过输入图像(如2D卷积),并比较模板图像下的输入图像的模板和补丁。OpenCV中实现了集中比较方法。(你可以查看文档了解更多详情)。它返回有一个灰度图像,其中每个像素表示该像素的领域和模板匹配的程度。
如果输入 图像的大小是(W*H)和模板图像的大小(w*h),则输出图像的大小为(W-w+1,H-h+1)。一旦得到结果,就可以使用cv2.minMaxLoc()函数来查找最大值/最小值。将它作为矩形的左上角,并将(w,h)作为矩形的宽度和高度。该矩形是你的模板区域。
注意:如果你使用cv2.TM_SQDIFF作为比较方法,最小值会给出最佳匹配。
3.OpenCV中的模板匹配
在这个栗子中将在图像中搜索下面同学的这张脸,所以创建模板如下:
我们将尝试所有比较方法,以便我们可以看到它们的结果如何:
# -*- coding: utf-8 -*- ''' OpenCV中的模板匹配 1.cv2.matchTemplate()、cv2.minMaxLoc() 2.OpenCV中实现了集中比较方法,返回一个灰度图像,其中每个像素表示该像素的领域和模板匹配的程度 3.输入图像大小若是(w*h),则输出图像大小为(W-w+1,H-h+1).一旦得到结果,可以使用cv2.minMaxLoc()函数来查找最大值和最小值 下面的栗子:在图像中搜索一张脸 可以看到cv2.TM_CCORR的结果并不如预期那样好 ''' import cv2 import numpy as np from matplotlib import pyplot as plt img = cv2.imread('1.jpg', 0) img2 = img.copy() template = cv2.imread('template.jpg', 0) w, h = template.shape[::-1] # 6中类型的比较 methods = ['cv2.TM_CCOEFF', 'cv2.TM_CCOEFF_NORMED', 'cv2.TM_CCORR', 'cv2.TM_CCORR_NORMED', 'cv2.TM_SQDIFF', 'cv2.TM_SQDIFF_NORMED'] for meth in methods: img = img2.copy() method = eval(meth) # 应用模板匹配 res = cv2.matchTemplate(img, template, method) min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(res) # 如果方法是cv2.TM_SQDIFF或cv2.TM_SQDIFF_NORMED if method in [cv2.TM_SQDIFF, cv2.TM_SQDIFF_NORMED]: top_left = min_loc else: top_left = max_loc bottom_right = (top_left[0] + w, top_left[1] + h) cv2.rectangle(img, top_left, bottom_right, 255, 2) plt.subplot(121), plt.imshow(res, 'gray') plt.title('Matching Result'), plt.xticks([]), plt.yticks([]) plt.subplot(122), plt.imshow(img, 'gray') plt.title('Detected Point'), plt.xticks([]), plt.yticks([]) plt.suptitle(meth) plt.show()
结果:
cv2.TM_CCOEFF
cv2.TM_CCOEFF_NORMED
cv2.TM_CCORR
cv2.TM_CCORR_NORMED
cv2.TM_SQDIFF
cv2.TM_SQDIFF_NORMED
你可以看到cv2.TM_CCORR的结果并不像我们预期的那样好。
4.多目标的模板匹配
在上一节中,我们搜索了梅西的脸部图像,该图像只在图像中出想过一次。假设你正在搜索一个多次出现的对象,cv2.minMaxLoc()不会给你所有的位置。在这种情况下,我们将使用阈值。所以在这个例子中,我们将使用游戏Mario的屏幕截图,并在其中找到硬币。
# -*- coding: utf-8 -*- ''' 多目标模板匹配 ''' import cv2 import numpy as np from matplotlib import pyplot as plt img_rgb = cv2.imread('2.jpg') img_gray = cv2.cvtColor(img_rgb, cv2.COLOR_BGR2GRAY) template = cv2.imread('temp.jpg', 0) w, h = template.shape[::-1] res = cv2.matchTemplate(img_gray, template, cv2.TM_CCOEFF_NORMED) threshold = 0.8 loc = np.where(res >= threshold) for pt in zip(*loc[::-1]): cv2.rectangle(img_rgb, pt, (pt[0] + w, pt[1] + h), (0, 0, 255), 2) cv2.imwrite('res.jpg', img_rgb) test = cv2.imread('2.jpg') cv2.imshow('original', test) cv2.imshow('res', img_rgb) cv2.waitKey(0) & 0xFF cv2.destroyAllWindows()
结果: