Matlab中使用FFT进行模板匹配

Question

我在Matlab中使用傅里叶域中的模板匹配遇到了困难。这里是我的图片（艺术家是DeviantArt上的RamalamaCreatures）：

我的目标是在负鼠的耳朵周围放置一个边界框，就像这个例子一样（我使用normxcorr2执行了模板匹配）：

这是我正在使用的Matlab代码：

clear all; close all;

template = rgb2gray(imread('possum_ear.jpg'));
background = rgb2gray(imread('possum.jpg'));

%% calculate padding
bx = size(background, 2); 
by = size(background, 1);
tx = size(template, 2); % used for bbox placement
ty = size(template, 1);

%% fft
c = real(ifft2(fft2(background) .* fft2(template, by, bx)));

%% find peak correlation
[max_c, imax]   = max(abs(c(:)));
[ypeak, xpeak] = find(c == max(c(:)));
figure; surf(c), shading flat; % plot correlation 

%% display best match
hFig = figure;
hAx  = axes;
position = [xpeak(1)-tx, ypeak(1)-ty, tx, ty];
imshow(background, 'Parent', hAx);
imrect(hAx, position);

代码未按预期运行 - 它未能识别正确的区域。这是失败的结果 - 错误的区域被框选：

这是匹配失败情况下相关性的表面绘图：

希望你能提供帮助！谢谢。

- Chris Parry

---

如果您查看图像处理工具箱中normxcorr2的实现，其中一个代码路径是在频域中计算的。您应该将其作为参考进行查看。 - Alex Taylor

Answer 1

你在代码中所做的实际上根本不是相关性。你正在使用模板并对输入图像进行卷积。如果您回想傅里叶变换，两个信号频谱的乘积相当于时间/空间域中两个信号的卷积。基本上，您所做的是将模板用作内核，并使用其来过滤图像。然后，您会找到此输出的最大响应，这就是被认为是模板所在位置的地方。响应被框定的区域之所以有意义，是因为该区域完全是白色的，并且使用模板作为内核和一个完全是白色的区域将为您提供非常大的响应，这就是为什么它很可能确定该区域为最大响应的原因。具体而言，该区域将具有许多高值（约255左右），自然使用模板补丁和该区域进行卷积将为您提供非常大的输出，因为该操作是加权和。因此，如果您在图像的暗区域中使用模板，则输出将很小-这是错误的，因为模板也由黑色像素组成。

---

然而，您可以使用傅里叶变换来定位模板的位置，但我建议您使用相位相关。基本上，您不是计算两个频谱的乘积，而是计算交叉功率谱。在频域中，两个信号之间的交叉功率谱R的定义如下：

^{来源: 维基百科}

Ga和Gb分别是原始图像和频率域中的模板，*表示共轭。 o是所谓的Hadamard乘积或逐元素乘积。我还想指出，分子和分母的除法也是逐元素的。使用交叉功率谱，如果您在此处找到产生绝对最大响应的(x,y)位置，则模板应位于背景图像中的该位置。

因此，您只需要更改计算“相关性”的代码行，以便它计算交叉功率谱。但是，我想指出非常重要的一点。当您执行normxcorr2时，相关性从图像的左上角开始。模板匹配从此位置开始，并与大小为模板的窗口进行比较，其中左上角是原点。在找到模板匹配的位置时，该位置是相对于匹配窗口的左上角的。一旦计算了normxcorr2，您通常会将最大响应的半行和半列加起来，以找到中心位置。

因为我们在模板匹配方面（滑动窗口、相关等）基本上都是使用FFT/频域执行相同操作，所以当您完成在此相关数组中找到峰值时，您必须也考虑到这一点。但是，您调用imrect来绘制一个矩形框，围绕着模板匹配的位置，无论如何都会输入一个边界框的左上角，所以这里没有必要进行偏移。因此，我们将稍微修改该代码，但在使用此代码查找匹配中心位置时，请记住偏移逻辑。

---

我也修改了你的代码，直接从StackOverflow读取图像，以便可以复现：

clear all; close all;

template = rgb2gray(imread('http://i.stack.imgur.com/6bTzT.jpg'));
background = rgb2gray(imread('http://i.stack.imgur.com/FXEy7.jpg'));

%% calculate padding
bx = size(background, 2); 
by = size(background, 1);
tx = size(template, 2); % used for bbox placement
ty = size(template, 1);

%% fft
%c = real(ifft2(fft2(background) .* fft2(template, by, bx)));

%// Change - Compute the cross power spectrum
Ga = fft2(background);
Gb = fft2(template, by, bx);
c = real(ifft2((Ga.*conj(Gb))./abs(Ga.*conj(Gb))));

%% find peak correlation
[max_c, imax]   = max(abs(c(:)));
[ypeak, xpeak] = find(c == max(c(:)));
figure; surf(c), shading flat; % plot correlation    

%% display best match
hFig = figure;
hAx  = axes;

%// New - no need to offset the coordinates anymore
%// xpeak and ypeak are already the top left corner of the matched window
position = [xpeak(1), ypeak(1), tx, ty];
imshow(background, 'Parent', hAx);
imrect(hAx, position);

有了那个，我得到了以下图片：

当显示交叉功率谱的曲面图时，我也会得到以下内容：

在相位相关中，存在一个明确定义的峰值，其余输出响应非常小。这实际上是相位相关的一种特性，因此最大值的位置是明确定义的，这就是模板所在的位置。

---

希望这可以帮到你！

- rayryeng

---

@ChrisParry - 不客气！顺便说一下，我知道为什么匹配会给你底部左角。我已经在我的答案中添加了这个，但是因为你正在使用模板作为内核，使用此内核和非常白色的区域进行卷积自然会给您一个大的响应，因此找到最大响应更自然地位于该区域。这就是为什么相位相关在这里更合适，因为它自然地发现了两个图像移动了多少。我很高兴能够帮助。祝你好运！ - rayryeng

@ChrisParry - 我还添加了交叉功率谱表面图的外观。最好看一下，这样你就可以确信我们得到的位置绝对是模板所在的位置。 - rayryeng

太好了，谢谢！我相信这个答案也会帮助很多其他人。 - Chris Parry

这对我来说很有效，谢谢。因为匹配速度非常快，在不知道模板方向的情况下，我可以循环遍历各种旋转。然而，匹配似乎也非常敏感于模板的比例尺度。除了循环遍历模板的各种变换之外，有没有好的方法来处理这个问题？ - Chris Parry

@RajnikantSharma 模板可能不具备足够的区分能力，无法产生良好的结果（即其傅里叶变换不够复杂）。您可能也会因此出现除以零的错误，这可能是您得到“NaN”的原因。顺便说一句，如果您喜欢这个答案，请向前支付并点赞。祝你好运。 - rayryeng

Answer 2

最终我使用Python实现了与@rayryeng类似的思路，使用scipy.fftpack.fftn() / ifftn()函数，在相同的目标和模板图像上得到了以下结果：

import numpy as np
import scipy.fftpack as fp
from skimage.io import imread
from skimage.color import rgb2gray, gray2rgb
import matplotlib.pylab as plt
from skimage.draw import rectangle_perimeter

im = 255*rgb2gray(imread('http://i.stack.imgur.com/FXEy7.jpg'))    # target
im_tm = 255*rgb2gray(imread('http://i.stack.imgur.com/6bTzT.jpg')) # template

# FFT 
F = fp.fftn(im)                   
F_tm = fp.fftn(im_tm, shape=im.shape)

# compute the best match location
F_cc = F * np.conj(F_tm)
c = (fp.ifftn(F_cc/np.abs(F_cc))).real
i, j = np.unravel_index(c.argmax(), c.shape)
print(i, j)
# 214 317

# draw rectangle around the best match location
im2 = (gray2rgb(im)).astype(np.uint8)
rr, cc = rectangle_perimeter((i,j), end=(i + im_tm.shape[0], j + im_tm.shape[1]), shape=im.shape)
for x in range(-2,2):
    for y in range(-2,2):
        im2[rr + x, cc + y] = (255,0,0)

# show the output image
plt.figure(figsize=(10,10))
plt.imshow(im2)
plt.axis('off')
plt.show()

同时，下面的动画展示了在一段鸟群视频中提取出的一组（目标）帧中定位鸟类模板图像所得到的结果。

需要注意的一点是：输出结果非常依赖于要匹配的对象与模板图像的大小和形状的相似度，如果它与模板图像差别很大，可能根本无法匹配。

- Sandipan Dey

---

如果你能提供你的Python代码，那将会很有用。当然这个问题是关于MATLAB，但是看到Python的实现也是不错的。我经常被要求把上述逻辑翻译成Python，但我没有时间。 - rayryeng