OpenCV中的图像处理 —— 傅里叶变换+模板匹配

现在也在逐渐深入啦，希望跟大家一起进步越来越强

1. 傅里叶变换

关于傅里叶变换最重要的两个概念：时域与频域。以时间作为参照来观察动态世界的方法我们称其为时域分析，而频域是什么呢，它是描述信号在频率方面特性时用到的一种坐标系，频域图显示了在一个频率范围内每个给定频带内的信号量。贯穿时域与频域的方法之一就是大名鼎鼎的傅里叶分析，它可以分为傅里叶级数和傅里叶变换，傅里叶变换也就是我们这一部分要说的东西

傅里叶变换是分析线性系统的一个有力工具。它告诉我们任何周期函数，都可以看作是不同振幅，不同相位正弦波的叠加。从数学意义上说，傅里叶变换将一个任意的周期函数分解成为无穷个正弦函数的和的形式；从物理效果上看，傅里叶变换实现了将信号从空间域到频率域的转换

在计算机视觉中傅立叶变换用于分析各种滤波器的频率特性，对于图像，使用2D离散傅里叶变换(DFT)查找频域（还有一种称为快速傅立叶变换(FFT)的快速算法）这一段文字是不是不太好理解，因为里面涉及太多比较深奥的东西了，傅里叶变换本身是比较难的一个点，在这里我就不细说了，我们只说说在计算机视觉领域我们是怎么用它的，想要深入了解的同学来看看这篇文章：深入浅出的讲解傅里叶变换（真正的通俗易懂）

对于正弦信号，如果幅度在短时间内变化比较快，则可以说它是高频信号，如果变化缓慢，则为低频信号，我们可以将相同的想法扩展到图像，图像中的振幅在哪里急剧变化？当然是在边缘点或噪声，因此，可以说边缘和噪声是图像中的高频内容

1.1 Numpy实现傅里叶变换

Numpy提供了FFT软件包来查找傅里叶变换，np.fft.fft2()为我们提供了频率转换，它将是一个复杂的数组，它的第一个参数是输入图像（灰度图像），第二个参数是可选的，它决定输出数组的大小。如果它大于输入图像的大小，则在计算FFT之前用零填充输入图像。如果小于输入图像，将裁切输入图像。如果未传递任何参数，则输出数组的大小将与输入的大小相同，但是现在获得的结果它的零频率分量（DC分量）将位于左上角，为了便于分析我们要把它居中，居中处理关系到np.fft.fftshift()函数

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread(r'E:\image\test16.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，结果同样是一个复数

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread(r'E:\image\test15.png', 0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20 * np.log(np.abs(fshift))
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.abs(img_back)
plt.subplot(131), plt.imshow(img, cmap='gray'),
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])
plt.show()

1.2 OpenCV实现傅里叶变换

OpenCV为此提供了cv.dft()和cv.idft()函数。它返回与前一个相同的结果，但是有两个通道。第一个通道是结果的实部，第二个通道是结果的虚部。输入图像首先应转换为np.float32

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread(r'E:\image\test17.png', 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.dft()函数的参数怎么传递？cv.magnitude()函数是怎么用的？

cv.dft()函数的作用是对一维或者二维浮点数数组进行正向或反向离散傅里叶变换，其中包括4个参数，第一个即源图像，第二个参数是OutputArray类型的dst，函数调用后返回的运算结果存在这里，它的尺寸和类型取决于第三个参数flags转换标识符，它的默认值为0（参考自：opencv：dft()函数详解）

cv.magnitude()函数用来计算二维矢量的幅值，其中包括3个参数，第一个是InputArray类型的x，表示矢量的浮点型X坐标值，也就是实部，第二个参数是InputArray类型的y，表示矢量的浮点型Y坐标值，也就是虚部，第三个参数是输出的幅值

接下来我们需要做OpenCV中DFT的逆变换，上一节用了高通滤波器HPF，这一部分我们会将低通滤波器LPF应用到图像中

注意：通常，OpenCV函数cv.dft()和cv.idft()比Numpy函数更快，大约快3倍，但是Numpy函数更容易使用

我们把这一部分的代码放在后面，与DFT的性能优化放在一起更容易理解

1.3 DFT的性能优化

对于某些数组尺寸，DFT的计算性能较好，例如当数组大小为2的幂时，速度最快，对于大小为2、3和5的乘积的数组，也可以非常有效地进行处理，关于代码的性能问题，我们可以在找到DFT之前将数组的大小修改为任何最佳大小(通过填充零)，对于OpenCV，我们必须手动填充零，但是对于Numpy，指定FFT计算的新大小，它将自动为您填充零

关于寻找最优大小，OpenCV为此提供了一个函数：cv.getOptimalDFTSize()

import cv2
import numpy as np
from matplotlib import pyplot as plt
img = cv2.imread(r'E:\image\test17.png', 0)
rows, cols = img.shape
print(rows, cols)
# 计算DFT效率最佳的尺寸
nrows = cv2.getOptimalDFTSize(rows)
ncols = cv2.getOptimalDFTSize(cols)
print(nrows, ncols)
nimg = np.zeros((nrows, ncols))
nimg[:rows, :cols] = img
img = nimg
# OpenCV计算快速傅里叶变换，输入图像应首先转换为np.float32，然后使用函数cv2.dft()和cv2.idft()。
# 返回结果与Numpy相同，但有两个通道。第一个通道为有结果的实部，第二个通道为有结果的虚部。
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, 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()
rows, cols = img.shape
crow, ccol = rows // 2, cols // 2
# 首先创建一个mask，中心正方形为1，其他均为0
# 如何删除图像中的高频内容，即我们将LPF应用于图像。它实际上模糊了图像。
# 为此首先创建一个在低频时具有高值的掩码，即传递LF内容，在HF区域为0。
mask = np.zeros((rows, cols, 2), np.uint8)
mask[crow - 30:crow + 30, ccol - 30:ccol + 30] = 1
# 应用掩码Mask和求逆DTF
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()

2. 模板匹配

模板匹配是一种用于在较大图像中搜索和查找模板图像位置的方法。为此，OpenCV带有一个函数cv.matchTemplate()，它只是将模板图像滑动到输入图像上（就像在2D卷积中一样），然后在模板图像下比较模板和输入图像的拼图，它返回一个灰度图像，其中每个像素表示该像素的邻域与模板匹配的程度。但是在这个方法的内涵中，到底是用什么样的方法去做做匹配的呢？这个就由函数的其中一个参数来决定了

如果输入图像的大小为 (W * H) ，而模板图像的大小为 (w * h) ，则输出图像的大小将为(W-w + 1，H-h + 1) ，在我们得到结果后，可以使用cv.minMaxLoc()函数查找最大/最小值在哪，将其作为矩形的左上角，并以 (w，h) 作为矩形的宽度和高度

匹配方式（比较方法）与cv.matchTemplate()的参数有关，我们先来看看都有哪些参数，并且应用不同的参数会有什么样不同的效果。该函数第一个参数是源图像，第二个参数是模板图像，第三个参数是匹配的结果图像，第四个参数是用于指定比较的方法

cv::TM_SQDIFF：该方法使用平方差进行匹配，因此最佳的匹配结果在结果为0处，值越大匹配结果越差
cv::TM_SQDIFF_NORMED：该方法使用归一化的平方差进行匹配，最佳匹配也在结果为0处
cv::TM_CCORR：相关性匹配方法，该方法使用源图像与模板图像的卷积结果进行匹配，因此，最佳匹配位置在值最大处，值越小匹配结果越差
cv::TM_CCORR_NORMED：归一化的相关性匹配方法，与相关性匹配方法类似，最佳匹配位置也是在值最大处
cv::TM_CCOEFF：相关性系数匹配方法，该方法使用源图像与其均值的差、模板与其均值的差二者之间的相关性进行匹配，最佳匹配结果在值等于1处，最差匹配结果在值等于-1处，值等于0直接表示二者不相关
cv::TM_CCOEFF_NORMED：归一化的相关性系数匹配方法，正值表示匹配的结果较好，负值则表示匹配的效果较差，也是值越大，匹配效果也好

资料摘自：【OpenCV3】模板匹配——cv::matchTemplate()详解

2.1 单对象的模板匹配

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread(r'E:\image\test15.png', 0)
img2 = img.copy()
template = cv.imread(r'E:\image\temple.png', 0)
w, h = template.shape[::-1]
# 列表中所有的6种比较方法
methods = ['cv.TM_CCOEFF', 'cv.TM_CCOEFF_NORMED', 'cv.TM_CCORR',
           'cv.TM_CCORR_NORMED', 'cv.TM_SQDIFF', 'cv.TM_SQDIFF_NORMED']
for meth in methods:
    img = img2.copy()
    method = eval(meth)
    # 应用模板匹配
    res = cv.matchTemplate(img, template, method)
    min_val, max_val, min_loc, max_loc = cv.minMaxLoc(res)
    # 如果方法是TM_SQDIFF或TM_SQDIFF_NORMED，则取最小值
    if method in [cv.TM_SQDIFF, cv.TM_SQDIFF_NORMED]:
        top_left = min_loc
    else:
        top_left = max_loc
    bottom_right = (top_left[0] + w, top_left[1] + h)
    cv.rectangle(img, top_left, bottom_right, [255, 0, 0], 2)
    plt.subplot(121), plt.imshow(res, cmap='gray')
    plt.title('Matching Result'), plt.xticks([]), plt.yticks([])
    plt.subplot(122), plt.imshow(img, cmap='gray')
    plt.title('Detected Point'), plt.xticks([]), plt.yticks([])
    plt.suptitle(meth)
    plt.show()

2.2 多对象的模板匹配

上一节我们匹配了梅大人的面部，但是如果图像中有很多满足匹配条件的模板呢？这个时候cv.minMaxLoc()不会为我们提供所有位置，我们会使用阈值化

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img_rgb = cv.imread(r'E:\image\num.png')
img_gray = cv.cvtColor(img_rgb, cv.COLOR_BGR2GRAY)
template = cv.imread(r'E:\image\temple2.png', 0)
w, h = template.shape[::-1]
res = cv.matchTemplate(img_gray, template, cv.TM_CCOEFF_NORMED)
threshold = 0.8
loc = np.where(res >= threshold)
for pt in zip(*loc[::-1]):
    cv.rectangle(img_rgb, pt, (pt[0] + w, pt[1] + h), (0, 0, 255), 2)
cv.imshow('res.png', img_rgb)
cv.waitKey(0)
cv.destroyWindow()

代码解析：第8行的shape()函数是numpy.core.fromnumeric中的函数，它的功能是查看矩阵或者数组的维数，第11行的locloc是满足“res >= threshold”的像素点的索引集合，第12行的函数zip()用可迭代的对象作为参数，将对象中对应的元素打包成一个个元组，然后返回由这些元组组成的列表

（注：文章内容参考OpenCV4.1中文官方文档）
如果文章对您有所帮助，记得一键三连支持一下哦

标签 NC, OpenCV, PE

OpenCV中的图像处理 —— 傅里叶变换+模板匹配