Matlab 模板匹配使用 FFT

Posted 2023-02-25

【中文标题】Matlab 模板匹配使用 FFT

【英文标题】：Matlab Template Matching Using FFT

【发布时间】：2015-09-19 05:10:40

【问题描述】：

我在 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);

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

这是失败匹配的相关性的表面图：

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

【问题讨论】：

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

【参考方案1】：

您在代码中所做的实际上根本不是相关性。您正在使用模板并对输入图像执行卷积。如果你回忆一下傅里叶变换，两个信号的频谱相乘相当于两个信号在时/空域中的卷积。

基本上，您所做的是将模板用作内核并使用它来过滤图像。然后，您将找到此输出的最大响应，这就是模板所在的位置。响应被装箱的地方是有道理的，因为该区域完全是白色的，并且使用模板作为具有完全白色区域的内核会给您一个非常大的响应，这就是为什么它很可能将该区域确定为最大值回复。具体来说，该区域将有很多高值（~255 左右），并且自然地使用模板补丁执行卷积，并且由于该操作是加权和，因此该区域将为您提供非常大的输出。因此，如果您在图像的暗区使用模板，则输出会很小 - 这是错误的，因为模板也由暗像素组成。

---

不过，您当然可以使用傅立叶变换来定位模板所在的位置，但我建议您改用Phase Correlation。基本上，不是计算两个光谱的乘积，而是计算交叉功率谱。频域中两个信号之间的交叉功率谱R定义为：

^{来源：Wikipedia}

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

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

因为我们在模板匹配（滑动窗口，相关等）与FFT /频域或多或少做相同的操作，当你在这个相关数组中找到峰值时，你还必须采取考虑到这一点。但是，您调用imrect 在模板匹配的位置周围绘制一个矩形无论如何都会占用边界框的左上角，因此无需在此处进行偏移。因此，我们将稍微修改该代码，但在以后使用此代码时要记住偏移逻辑，以便找到匹配的中心位置。

---

我还修改了您的代码以直接从 *** 读取图像，以便重现：

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);

这样，我得到以下图像：

在显示互功率谱的曲面图时，我还得到以下信息：

有一个明确定义的峰值，其中其余输出的响应非常小。这实际上是相位相关的一个属性，所以很明显，最大值的位置是明确定义的，这就是模板所在的位置。

---

希望这会有所帮助！

【讨论】：

@ChrisParry - 我的荣幸！顺便说一句，我知道为什么匹配给了你左下角。我将此添加到我的答案中，但是由于您将模板用作内核，因此使用此内核执行卷积以及非常白的区域自然会给您很大的响应，因此找到最大响应会更自然在那个区域。这就是相位相关更适合这里的原因，因为它自然会发现两个图像发生了多少变化。我很高兴能帮上忙。祝你好运！

@ChrisParry - 我还添加了交叉功率谱曲面图的样子。最好看到它，这样你就可以说服自己我们得到的位置肯定是模板所在的位置。

完美，谢谢！我相信这个答案也会对其他人有所帮助。

@ChrisParry - 在我解释模板匹配位置的方式上有一个小错误。匹配的坐标是窗口的左上角，因此只需偏移一半的行和列即可获得中心。但是，您的imrect 电话无论如何都需要左上角，因此我修改了该逻辑。结果还是一样，我编辑了帖子。

@RajnikantSharma 该模板的判别力可能不足以提供良好的结果（即其傅立叶变换不够复杂）。结果，您可能还会得到除以零的错误，这可能就是您得到NaN 的原因。顺便说一句，如果您喜欢这个答案，请转发并投票。祝你好运。

【参考方案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()

此外，下面的动画显示了在从一组鸟的视频中提取的一组（目标）帧中定位鸟的模板图像时获得的结果。

需要注意的一点：输出在很大程度上取决于要与模板匹配的对象的大小和形状的相似性，如果与模板图像差异很大，则模板可能不是完全匹配。

【讨论】：

如果你能包含你的 Python 源代码将会很有用。问题当然是 MATLAB，但如果能缩小差距并看到 Python 实现也会很好。经常有人要求我将上述逻辑翻译成 Python，但我没有时间。