(三)OpenCV中的图像处理之图像变换及模板匹配

注释：本文翻译自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()

结果：