使用归一化互相关匹配对象外形

Question

我正在使用 Matlab 中的 normxcorr2 函数进行模板匹配。但是，我想做的与normxcorr2 做的不同。内置的normxcorr2 计算互相关，考虑到矩形模板中的所有像素。但我只希望某些像素参与归一化互相关过程。

例如，我希望在关联时仅将下图中的环状白色区域用作模板。 （里面的黑色区域不应该用于计算）

因此，如果我将上面的模板与下图相关联，我可能会得到大约 1.0 的归一化值（考虑到两个图像中的外圆直径相同）

我已经检查了这个解决方案：- matlab template matching only for 0 (or 1) in matrix 但它不是通用的。

任何人都可以帮助我在 matlab 中提供更通用的解决方案吗？这可用于通过外部形状跟踪对象

编辑：- 这是给想看的人的完整图像。我只想通过外圆边缘检测球，而不是内部细节。

Answer 1

这是上一篇文章的衍生品，我在这里提供了答案：matlab template matching only for 0 (or 1) in matrix

然而，这个解决方案使用了for 循环，效率很低。因此，我们可以使用im2col、bsxfun 和col2im 来帮助我们更快地执行此操作。 im2col 获取图像中的重叠区域并将它们分别放置到单独的列中。本质上，这会像使用任何类型的空间图像过滤一样获取图像的滑动窗口，并收集滑动窗口内的所有像素并将每个窗口放置为单独的列。

假设您的模板大小为M x N，您要搜索的图片大小为R x C，假设您的图片模板名为imTemplate，而您要搜索的图片为imSearch，我们可以进行以下设置。我们还假设两个图像都是二进制的。

[M, N] = size(imTemplate);
[R, C] = size(imSearch);

%// Cast to double for precision
imTemplate = im2double(imTemplate);
imSearch = im2double(imSearch);

neigh = im2col(imSearch, [M, N]);
templateCol = imTemplate(:); %// Ensures we place template into single column

现在，您希望排除在圆形边界内的所有像素。因此，我们可以做的是反转图像，使黑色像素变为白色，然后移除边框周围的所有像素。这应该会给我们圆的内部。

imInvert = ~imTemplate;
imInvertNoBorder = imclearborder(imInvert, 8); %// Search 8-pixel neighbourhood

我们将使用它来确定要从搜索中删除哪些像素。这可以通过以下方式完成：

rowsToRemove = imInvertNoBorder(:) == 1;

现在，我们可以做的是最终移除那些在圆圈内部的像素，这些像素不会在我们的相关方案中被搜索到。

neigh(rowsToRemove,:) = [];

我们现在可以做的是计算所有这些列的 NCC。回想一下，两个信号之间的 NCC 如下：

_{（来源：www.jot.fm）}

因此，我们需要从每个邻域中减去平均值，我们还需要从每一列中减去平均值。然后我们计算如上所示的公式。我们可以像这样在 MATLAB 中轻松实现向量化：

neighMeanSubtract = bsxfun(@minus, neigh, mean(neigh));
templateMeanSubtract = templateCol - mean(templateCol);

我们可以如下计算每个邻域的 NCC 分子（在求和之前）：

numerator = bsxfun(@times, neighMeanSubtract, templateMeanSubtract);

现在，我们要做的就是对所有列求和，这将给我们最终的分子：

sumNumerator = sum(numerator);

分母可以这样计算：

denominator1 = sqrt(sum(neighMeanSubtract.*neighMeanSubtract));
denominator2 = sqrt(sum(templateMeanSubtract.*templateMeanSubtract));
sumDenominator = denominator1 .* denominator2;

最后，我们的 NCC 可以这样计算：

NCC = sumNumerator ./ sumDenominator;

您会注意到这是单行值。每行对应于邻域定义的输出。因此，我们还需要将其重新整形为矩阵，因此您可以使用col2im：

finalOutput = col2im(NCC, [M, N], [R, C]);

上述语句将采用 NCC 中定义的 M x N 的重叠邻域，并对其进行重塑，使其成为 R x C 矩阵。有时，您会遇到除以零错误，尤其是当邻域搜索窗口全一致时。因此，您将获得NaN 号码。假设没有变化的区域在图像处理中没有相关性，所以让我们将这些位置归零：

finalOutput(isnan(finalOutput)) = 0;

如果您想找到相关性最高的位置，只需执行以下操作：

[rowNCC, colNCC] = find(finalOutput == max(finalOutput(:)));

如果您想解释负相关，这完全取决于您的应用程序。如果您想确保您的模板匹配算法是旋转不变，那么您实际上应该检查绝对值的最大值。负相关仅仅意味着模板和邻域之间的匹配被简单地旋转。因此，找到最佳社区的更好方法是：

maxCoeff = max(abs(finalOutput(:)));
[rowNCC, colNCC] = find(abs(finalOutput) == maxCoeff);

为了您的复制和粘贴乐趣，这里是完整的代码：

function [rowNCC, colNCC] = testCorr(imTemplate, imSearch)
    [M, N] = size(imTemplate);
    [R, C] = size(imSearch);

    %// Cast to double for precision
    imTemplate = im2double(imTemplate);
    imSearch = im2double(imSearch);

    neigh = im2col(imSearch, [M, N]);
    templateCol = imTemplate(:); %// Ensures we place template into single column

    imInvert = ~imTemplate;
    imInvertNoBorder = imclearborder(imInvert, 8); %// Search 8-pixel neighbourhood
    rowsToRemove = imInvertNoBorder(:) == 1;
    neigh(rowsToRemove,:) = [];

    neighMeanSubtract = bsxfun(@minus, neigh, mean(neigh));
    templateMeanSubtract = templateCol - mean(templateCol);

    numerator = bsxfun(@times, neighMeanSubtract, templateMeanSubtract);
    sumNumerator = sum(numerator);

    denominator1 = sqrt(sum(neighMeanSubtract.*neighMeanSubtract));
    denominator2 = sqrt(sum(templateMeanSubtract.*templateMeanSubtract));
    sumDenominator = denominator1 .* denominator2;

    NCC = sumNumerator ./ sumDenominator;    

    finalOutput = col2im(NCC, [M, N], [R, C]);
    finalOutput(isnan(finalOutput)) = 0;

    maxCoeff = max(abs(finalOutput(:)));
    [rowNCC, colNCC] = find(abs(finalOutput) == maxCoeff);
end

祝你好运！

Answer 2

好的，让我们试一试...这个解决方案尝试使用现有的normxcorr2 实现并对其进行修改以解决yoru问题。

归一化互相关的公式为：

在这种情况下，您希望更改每个窗口的集成边界。这反过来会影响标准偏差和相关性本身。让我们分几个步骤来解决它：

第 1 步：正确建立相关性

我们可以通过修改模板图片来做到这一点：

template_fix = template;
mean_template_mask = mean(template(mask == 1));
template_fix(mask == 0) = mean_template_mask;
result = normxcorr2(template_fix, query)

请注意，通过进行此更改，我们使模板的平均值等于掩码侧模板的平均值。这样，掩码外的所有模板像素都不会对积分做出贡献，因为它们等于平均值​​。

步骤 #2：修复模板标准

size_mask = sum(mask(:));
size_template = numel(template);
std_template = std2(template);
std_template_masked = sqrt(sum((template_fix(:) - mean_template_mask).^2)/size_mask);
result = result * (std_template/std_template_masked);

第 3 步：修复查询标准

sum_filt = ones(size(template));
std_query = filter2(query.^2, sum_filt) - filter2(query, sum_filt).^2/size_template;
std_query = sqrt(std_query/size_template);

std_query_mask = filter2(query.^2, mask) - filter2(query, mask).^2/size_mask;    
std_query_mask = sqrt(std_query_mask/size_mask);

result = result .* std_query ./ std_query_mask;

我的 Matlab 没有响应，所以我没有机会在实践中对其进行测试。除非我错过了一些错误，否则它应该在数学上是等价的。

此解决方案会进行一些额外的卷积，但不会多次处理重叠像素。

如果您多次使用相同的模板图像，那么您可以重构第 1 步和第 2 步，使其只运行一次以进行预处理。虽然两者的计算成本都不应该很高。

不同的方法：直截了当

这是一种不同的、直接的方法，它不使用原始的normxcorr2 函数。这段代码可以很容易地针对内存使用进行优化，但会牺牲可读性。

% q for query, t for template and mask for mask
% shape = 'full' or 'same' or 'valid'

t_mask = t .* mask;
n      = numel(mask);
tq_sum = filter2(t_mask,q, shape);

q_sum  = filter2(mask, q, shape);    
q_mean = q_sum/n;
t_sum  = sum(t_mask(:));
t_mean = t_sum/n;

res1 = tq_sum - t_mean*q_sum - t_sum*q_mean + t_mean*q_mean*n;

t_std = sqrt((sum(t_mask(:).^2) - sum(t_mask(:)).^2/n)/n);
q_std = sqrt((filter2(mask, q.^2, shape) - q_sum.^2/n)/n);

res = res1 ./ (n * t_std * q_std)