Fast Normalized Cross-Correlation

本文提出了一种从变换域卷积中获得归一化互相关的新算法,该算法在某些情况下提供了一个数量级的加速。归一化互相关(NCC)在特征匹配中是一种常用方法,但在图像尺度、旋转和透视畸变方面不是不变的。本文讨论了NCC的优缺点,并与其他特征跟踪方法进行了比较。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

Abstract

虽然互相关可以在变换域中有效实现,但特征匹配应用中首选的互相关归一化形式并没有一个简单的频域表达式。为此,在空间域中计算了归一化互相关。通过对搜索窗口中图像和image2的预计算积分,证明了非归一化交叉相关可以有效地归一化。

1 Introduction

两个信号之间的相关性(互相关)是特征检测的标准方法[6,7],也是更复杂的技术(如[3])的组成部分。相关的教科书描述了卷积定理,以及使用快速傅立叶变换在频域中高效计算相关的可能性。不幸的是,模板匹配中首选的相关(相关系数)的归一化形式并没有相应的简单有效的频域表达式。因此,在空间域(如[7],p. 585)计算归一化互相关。由于空间域卷积的计算成本较高,因此也开发了几种不精确但快速的空间域匹配方法。本文介绍了一种从变换域卷积中获得归一化互相关的新算法[10]。在某些情况下,新算法在归一化互相关的空间域计算上提供了一个数量级的加速(第5节)。

由于我们展示的是一种熟悉且广泛使用的算法版本,因此我们不打算调查关于特征选择、美白、快速卷积技术、扩展、替代技术或应用程序的文献。关于这些主题的文献可以通过介绍性文本和手册[16,7,13]以及最近的论文,如[1,19]。然而,由于所提倡的各种特征跟踪方案,可能有必要确定,即使不是所有的应用程序,标准化的相互关联仍然是一种可行的选择。这是在第3节中完成的。

为了使本文具有自主性,第2节介绍了归一化互相关,第4节简要回顾了变换域和其他快速卷积方法以及相位相关技术。大多数读者可以跳过这些部分。第5节描述如何从相关的转换域计算中获得规范化的相互关联。第6节给出性能结果。

2 Template Matching by Cross-Correlation

模板匹配中使用的互相关是由距离度量(平方欧几里得距离)驱动的

\small d_{f,t}^2(u,v)=\sum_{x,y}[f(x,y)-t(x-u,y-v)]^2

(其中f是图像,并且总和在包含位于u,v的特征t的窗口下的x,y之上)。 在扩展d2

\small d_{f,t}^2(u,v)=\sum_{x,y}[f^2(x,y)-2f(x,y)t(x-u,y-v)+t^2(x-u,y-v)]

\small \sum t^2(x-u,y-v)是常数。如果\small \sum f^2(x,y)近似为常数,则剩余的互相关项

\small c(u,v)=\sum_{x,y}f(x,y)t(x-u,y-v)(1)

是图像与特征之间相似性的度量。

使用(1)进行模板匹配有几个缺点:

.如果图像能量\small \sum f^2(x,y)随位置变化,使用(1)匹配可能失败。例如,特征与图像中精确匹配区域的相关性可能小于特征与亮点的相关性。

.c(u,v)的范围取决于特征的大小。

.Eq.(1)对于图像幅值的变化不是不变的,例如在整个图像序列中由于光照条件的变化而引起的图像幅值的变化。

相关系数克服了这些困难,将图像和特征向量归一化为单位长度,得到一个余弦相关系数

\small \gamma (u,v)=\frac{\sum_{x,y}[f(x,y)-\bar{f}_{u,v}][t(x-u,y-v)-\bar{t}]} {\{\sum_{x,y}[f(x,y)-\bar{f}_{u,v}]^2 \sum_{x,y}[t(x-u,y-v)-t]^2\}^{0.5}}(2)

其中\small \bar{t}是特征的平均值,\small \bar{f}_{u,v}是特征下区域中f(x,y)的平均值。 我们将(2)称为归一化互相关。

3 Feature Tracking Approaches and Issues

很明显,归一化互相关(NCC)并不是一种理想的特征跟踪方法,因为它在图像尺度、旋转和透视畸变方面不是不变的。这些限制已经在各种方案中得到解决,包括一些将NCC作为一个组成部分的方案。本文不提倡选择NCC而不是其他方法。相反,下面的讨论将指出各种特性跟踪方法中涉及的一些问题,并得出结论,NCC是某些应用程序的合理选择。

SSDA。序列相似度检测算法(SSDA)[2]的基础是观察到只有在互相关联函数的最大值附近才需要全精度,而降低的精度可以在其他地方使用。[2]的作者描述了几种实现降低精度的方法。交叉关系的SSDA实现以随机顺序计算(1)的和,并使用部分计算作为蒙特卡洛估计,以确定特定匹配位置是否接近相关面的最大值。如果估计值表明某个位置与较差匹配,则在完成和之前终止特定位置的计算。

SSDA算法简单,在空间域的相互关联上提供了显著的加速。它的缺点是不能保证找到相关曲面的最大值。当相关曲面具有较浅的斜率和较大的最大值时,SSDA表现较好。虽然这一条件在许多应用中可能得到了满足,但显然,包含对象数组(卵石、砖块、其他纹理)的图像可以在相关表面生成多个窄极值,从而误导了SSDA方法。SSDA的第二个缺点是它有需要确定的参数(用于形成相关系数估计数的术语的数量,以及该估计数的早期终止阈值)。

Snake。蛇(活动轮廓模型)的缺点是它们不能跟踪没有可定义轮廓的物体。有些物体没有一个明确的边界(无论是由于内在的模糊还是由于光照条件),但仍然有一个颜色的特征分布,可以通过相互关联来跟踪。主动轮廓模型解决了一个比简单模板匹配更普遍的问题,因为它们提供了随时间变化的变形轮廓的表示。相互关联可以跟踪随时间而变形的对象,但是具有这里不讨论的明显而重要的条件。交叉关联还可以很容易地跟踪一个特征,该特征在帧间移动的幅度是它自身大小的很大一部分,而这种移动量可以将一条蛇置于其收敛盆地之外。

Wavelets and other multi-resolution schemes。虽然小波卷积定理的存在性还有待讨论(例如,[11];在一些方案中,小波卷积实际上是使用傅立叶卷积定理实现的),但利用粗到细多分辨率搜索,小波和其他多分辨率表示可以实现有效的特征跟踪。然而,多分辨率技术要求图像包含足够的低频信息,以指导搜索的初始阶段。正如在第6节中所讨论的,理想的特征有时是不可用的,人们必须求助于定义不明确的特征,这些特征可能只有很少的低频信息,例如在一个均匀表面上的小点的配置。

上面讨论的每一种方法都是由不同的作者提倡的,但是不同方法之间的比较较少。参考[19]在梯度搜索框架中推导出一种最优的特征跟踪方案,但该框架的局限性并未得到解决。通过对五种模板匹配算法在各种图像畸变情况下的实证研究,[4]发现,NCC在所有图像类别中表现最佳,虽然其中一种较便宜的算法对某些类型的畸变表现几乎与之相当。在[1]中讨论了运动跟踪的一般层次结构框架。在考虑梯度法的基础上,选择了一种基于相关性的匹配方法。

尽管NCC算法的时代已经过去,而且出现了更多的新技术来解决它的各种缺点,但可以说,一个合适的替代方案还没有得到普遍认可。NCC对图像序列的要求很少,没有用户需要搜索的参数。NCC可用来提供简单的特性跟踪,也可以作为更复杂(可能是多分辨率)匹配方案的组件使用,该匹配方案可以处理缩放和旋转不变性、特性更新和其他问题。在绝对差和等可选匹配准则中选择相关系数也被证明为最大似然估计[18]。我们承认NCC在许多应用程序中是默认的选择,在这些应用程序中,特征跟踪本身并不是研究的主题,也不是视觉和模式识别研究(例如[3])中的偶尔构建块。因此,一种快速算法很有趣。

4 Transform Domain Computation

考虑(2)中的分子并假设我们有图像\small f^{'}(x,y)\equiv f(x,y)-\bar f(u,v)\small t^{'}(x,y)\equiv t(x,y)-\bar{t}其中已经删除了平均值:

\small \overset{num}{\gamma}(u,v)=\sum_{x,y}f^{'}(x,y)t^{'}(x-u,y-v)(3)

对于大小为\small M^2的搜索窗口,大小为\small N^2(3)的特征需要大约\small N^2(M-N+1)^2个加法和\small N^2(M-N+1)^2个乘法。

式(3)是具有反转特征\small t^{'}(-x,-y)的图像的卷积,并且可以通过计算

\small F^{-1}\{F(f^{'}) F^{\ast }(t^{'})\}(4)

其中F是傅里叶变换。复共轭通过傅里叶变换性质\small Ff^{\ast }(-x)=F^{\ast}(w)完成特征的反转。

FFT算法的实现通常需要将\small f^'\small t^'扩展为0到2的公共幂。变换计算(3)的复杂度为\small 12M^{2}log_{2}M实乘法和\small 18M^{2}log_{2}M实加减。当M远大于N时,直接空间计算(3)的复杂度约为\small N^{2}M^{2}的乘法/加法,且直接法比变换法快。当N接近M时,随着M、N的增大,变换方法变得更加有效。

4.1 Fast Convolution

有几种著名的快速卷积算法不使用变换域计算[13]。这些方法可以分为两类:用乘法交换额外加法的算法,以及通过将一维卷积的各部分嵌入到较小的多维卷积的各个独立维度中,从而在(一维)卷积的O(n2)特征上找到一个较低点的方法。虽然这些算法比直接卷积快,但是在中等大小的[13]下,它们比变换域卷积慢,而且在任何情况下,它们都不能解决(2)的分母的计算问题。

4.2 Phase Correlation(相位相关)

由于(4)可以在变换域中进行有效的计算,因此提出了几种近似(2)中图像能量归一化的变换域方法。在互相关之前,通过高通滤波可以减小模板下图像能量的变化。这种滤波可以方便地添加到频域处理中,但是选择截止频率是有问题的,低截止频率可能会留下显著的图像能量变化,而高截止频率可能会删除对匹配有用的信息。

一种更可靠的方法是相位相关[9]。该方法首先将变换系数归一化为单位量值,然后在频域内计算相关系数。因此,这种相关性仅基于相位信息,对图像强度的变化不敏感。尽管经验表明这种方法是成功的,但它有一个缺点,即所有转换组件的权重都是相等的,而人们可能会认为不重要的组件应该得到更少的权重。在给定信号二阶矩和信号噪声的前提下,原则上应选择频谱预滤波,使期望的相关信噪比达到最大。该方法在[16]中进行了讨论,与经典的匹配滤波随机信号处理技术相似。与典型(ρ0:95)图像相关性最好的预滤器大约是拉普拉斯算子,而不是一个纯粹的美白。

5 Normalizing

再次检查(2)的分子,我们注意到特征的均值可以被预先计算,留下

\small \overset{num}{\gamma}(u,v)=\sum f(x,y)t^{'}(x-u,y-v)-\bar{f}_{u,v}\sum t^{'}(x-u,y-v)

由于\small t^{'}具有零均值并且因此零和和项\small \bar{f}_{u,v}\sum t^{'}(x-u,y-v)也是零,因此可以使用(4)来计算归一化互相关的分子。

通过检验(2)的分母,特征向量的长度可以在大约\small 3N^2次操作中预先计算出来(相对于互相关关系的代价来说,这是很小的),事实上,特征可以被预先归一化为长度为1。

有问题的量是表达式\small \sum_{x,y}[f(x,y)-\bar{f}_{u,v}]^2。图像平均值和局部能量(RMS)必须在每个u,v,即。在\small (M-N+1)^2个位置,产生了几乎\small 3N^{2}(M-N+1)^2个操作(将加、减、乘作为一个操作计算)。该计算超过了(3)的直接计算所需的计算,并且当变换方法适用时,它可能显着地超过由(4)指示的计算。 期望一种更有效的计算特征下的图像平均值和能量的方法。

这些量可以从包含图像和图像平方在搜索区域上的积分(运行和)的表中有效地计算出来,即

\small s(u,v)=f(u,v)+s(u-1,v)+s(u,v-1)-s(u-1,v-1)

\small s^{2}(u,v)=f^{2}(u,v)+s^{2}(u-1,v)+s^{2}(u,v-1)-s^{2}(u-1,v-1)

\small s(u,v)=s^{2}(u,v)=0当任意一个u,v < 0。则定位于u, v处的特征下图像的能量为

\small e_{f}(u,v)=s^{2}(u+N-1,v+N-1)-s^{2}(u-1,v+N-1)-s^{2}(u+N-1,v-1)+s^{2}(u-1,v-1)

并且类似于特征下的图像总和。

现在可以用很少的操作来计算有问题的量\small \sum_{x,y}[f(x,y)-\bar {f}_{u,v}]^2,因为它扩展成仅涉及特征下的图像和和平方的表达式。构建这些表需要大约\small 3M^2次运算,这比计算分子(4)的成本要低,也比计算每个u,v处的\small \sum_{x,y}[f(x,y)-\bar{f}_{u,v}]^2所需要的\small 3N^{2}(M-N+1)^2少得多。

这种从预先计算的运行和计算定和的技术已经在许多领域单独使用;在[5]中开发了一个计算机图形应用程序。如果以系统的行扫描顺序搜索关联曲面的最大值,则可以通过状态变量将表构造和引用结合起来,从而避免显式存储表。然而,在通用计算机上实现时,表的大小不是主要考虑的因素,搜索相关曲面的灵活性可能是有利的。注意s(u,v)和\small s^{2}(u,v)表达式是略微稳定的,这意味着它们的z变换\small H(z)=1/(1-z^{-1})(这里是一维的版本)在z = 1处有一个极点,而稳定性要求极点严格在单位圆[14]内。因此,计算应该使用大整数而不是浮点运算。

6 Performance

该算法的性能将在特殊效果图像处理的背景下进行讨论。将合成图像和经过处理的图像集成到特效序列中,往往需要对序列运动和特征进行准确的跟踪。在特技效果中使用自动特征跟踪是在电影中首创的,比如《扣人心弦》(Cliffhanger)、《阿甘正传》(Forest Gump)和《速度》(Speed)。最近,基于互相关的特征跟踪器被引入到商业图像合成系统中,如Flame/Flint[20]、Matador、Advance[21]和After Effects[22]。

本文所描述的算法是为电影《阿甘正传》(1994)开发的,并在后续的几个项目中得到了应用。这部电影的特效序列包括替换各种移动元素,在历史电影和视频序列中加入当代演员。从序列的一帧中手动选择的特征会自动跟踪到其余帧;这些信息被用作进一步处理的基础。

算法的相对性能是搜索窗口大小和特征大小与搜索窗口大小之比的函数。相对性能沿窗口尺寸轴增加(图1);更高分辨率的图将显示一个额外的波纹,反映搜索窗口大小和两个边界功率之间的关系。在更大的问题上,相对性能更大是可取的。表1说明了在特效特征跟踪应用程序中获得的性能。表2比较了我们的算法与高端商业图像合成包的性能。

注意,虽然小的特征尺寸(例如102)在理想的数字图像中已经足够,但在实践中,有时需要或更喜欢大得多的特征尺寸和搜索窗口:

.在电影和视频中使用的图像序列有时是从移动摄像机中获得的,由于相机抖动,在帧之间可能有相当大的平移。由于表示数字电影所需的高分辨率,即使是帧间的微小移动也可能对应许多像素的距离。

.所选择的特征当然受限于图像中可用的特征;不同的特征并不总是在首选的规模和位置可用。

.在典型的数字化图像中,由于相机或物体的运动,许多潜在的特征要么失焦,要么模糊(图2)。在存在模糊和噪声的情况下,大的特征更准确。

由于这些考虑,通常使用\small 20^2或更大的特征大小和\small 50^2或更大的搜索窗口。

在某些情况下,这种快速算法可以将高分辨率的特征跟踪从一夜之间减少到午餐时间。由于具有较低的(代理)分辨率和更快的机器,半自动化的特征跟踪在交互式系统中是可以接受的。其他领域中的某些应用程序也可以从本文描述的算法中获益。

/* * SPDX-License-Identifier: MIT * * Copyright (C) 2013-2024 OpenMV, LLC. * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal * in the Software without restriction, including without limitation the rights * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN * THE SOFTWARE. * * Template matching with NCC (Normalized Cross Correlation) using exhaustive and diamond search. * * References: * Briechle, Kai, and Uwe D. Hanebeck. "Template matching using fast normalized cross correlation." Aerospace * Lewis, J. P. "Fast normalized cross-correlation." * Zhu, Shan, and Kai-Kuang Ma. "A new diamond search algorithm for fast block-matching motion estimation." */ #include <stdio.h> #include <float.h> #include <limits.h> #include "imlib.h" static void set_dsp(int cx, int cy, point_t *pts, bool sdsp, int step) { if (sdsp) { // Small DSP // 4 // 3 0 1 // 2 pts[0].x = cx; pts[0].y = cy; pts[1].x = cx + step / 2; pts[1].y = cy; pts[2].x = cx; pts[2].y = cy + step / 2; pts[3].x = cx - step / 2; pts[3].y = cy; pts[4].x = cx; pts[4].y = cy - step / 2; } else { // Large DSP // 7 // 6 8 // 5 0 1 // 4 2 // 3 pts[0].x = cx; pts[0].y = cy; pts[1].x = cx + step; pts[1].y = cy; pts[2].x = cx + step / 2; pts[2].y = cy + step / 2; pts[3].x = cx; pts[3].y = cy + step; pts[4].x = cx - step / 2; pts[4].y = cy + step / 2; pts[5].x = cx - step; pts[5].y = cy; pts[6].x = cx - step / 2; pts[6].y = cy - step / 2; pts[7].x = cx; pts[7].y = cy - step; pts[8].x = cx + step / 2; pts[8].y = cy - step / 2; } } static float find_block_ncc(image_t *f, image_t *t, i_image_t *sum, int t_mean, uint32_t t_sumsq, int u, int v) { int w = t->w; int h = t->h; int num = 0; uint32_t f_sumsq = 0; if (u < 0) { u = 0; } if (v < 0) { v = 0; } if (u + w >= f->w) { w = f->w - u; } if (v + h >= f->h) { h = f->h - v; } // Find the mean of the current patch uint32_t f_sum = imlib_integral_lookup(sum, u, v, w, h); uint32_t f_mean = f_sum / (w * h); // Find the normalized sum of squares of the image for (int y = v; y < v + h; y++) { for (int x = u; x < u + w; x++) { int a = (int) f->data[y * f->w + x] - f_mean; int b = (int) t->data[(y - v) * t->w + (x - u)] - t_mean; num += a * b; f_sumsq += a * a; } } // Find the normalized cross-correlation return (num / (fast_sqrtf(f_sumsq) * fast_sqrtf(t_sumsq))); } float imlib_template_match_ds(image_t *f, image_t *t, rectangle_t *r) { point_t pts[9]; // Integral images i_image_t sum; imlib_integral_image_alloc(&sum, f->w, f->h); imlib_integral_image(f, &sum); // Normalized sum of squares of the template int t_mean = 0; uint32_t t_sumsq = 0; imlib_image_mean(t, &t_mean, &t_mean, &t_mean); for (int i = 0; i < (t->w * t->h); i++) { int c = (int) t->data[i] - t_mean; t_sumsq += c * c; } int px = 0; int py = 0; // Initial center point int cx = f->w / 2 - t->w / 2; int cy = f->h / 2 - t->h / 2; // Max cross-correlation float max_xc = -FLT_MAX; // Start with the Large Diamond Search Pattern (LDSP) 9 points. bool sdsp = false; // Step size == template width int step = t->w; while (step > 0) { // Set the Diamond Search Pattern (DSP). set_dsp(cx, cy, pts, sdsp, step); // Set the number of search blocks (5 or 9 for SDSP and LDSP respectively). int num_pts = (sdsp == true)? 5: 9; // Find the block with the highest NCC for (int i = 0; i < num_pts; i++) { if (pts[i].x >= f->w || pts[i].y >= f->h) { continue; } float blk_xc = find_block_ncc(f, t, &sum, t_mean, t_sumsq, pts[i].x, pts[i].y); if (blk_xc > max_xc) { px = pts[i].x; py = pts[i].y; max_xc = blk_xc; } } // If the highest correlation is found at the center block and search is using // LDSP then the highest correlation is found, if not then switch search to SDSP. if (px == cx && py == cy) { // Note instead of switching to the smaller pattern, the step size can be reduced // each time the highest correlation is found at the center, and break on step == 0. // This makes DS much more accurate, but slower. step--; } // Set the new search center to the block with highest correlation cx = px; cy = py; } r->x = cx; r->y = cy; r->w = t->w; r->h = t->h; if (cx < 0) { r->x = 0; } if (cy < 0) { r->y = 0; } if (cx + t->w > f->w) { r->w = f->w - cx; } if (cy + t->h > f->h) { r->h = f->h - cy; } imlib_integral_image_free(&sum); //printf("max xc: %f\n", (double) max_xc); return max_xc; } /* The NCC can be optimized using integral images and rectangular basis functions. * See Kai Briechle's paper "Template Matching using Fast Normalized Cross Correlation". * * NOTE: only the denominator is optimized. * */ float imlib_template_match_ex(image_t *f, image_t *t, rectangle_t *roi, int step, rectangle_t *r) { int den_b = 0; float corr = 0.0f; // Integral images i_image_t sum; i_image_t sumsq; imlib_integral_image_alloc(&sum, f->w, f->h); imlib_integral_image_alloc(&sumsq, f->w, f->h); imlib_integral_image(f, &sum); imlib_integral_image_sq(f, &sumsq); // Normalized sum of squares of the template int t_mean = 0; imlib_image_mean(t, &t_mean, &t_mean, &t_mean); for (int i = 0; i < (t->w * t->h); i++) { int c = (int) t->data[i] - t_mean; den_b += c * c; } for (int v = roi->y; v <= (roi->y + roi->h - t->h); v += step) { for (int u = roi->x; u <= (roi->x + roi->w - t->w); u += step) { int num = 0; // The mean of the current patch uint32_t f_sum = imlib_integral_lookup(&sum, u, v, t->w, t->h); uint32_t f_sumsq = imlib_integral_lookup(&sumsq, u, v, t->w, t->h); uint32_t f_mean = f_sum / (float) (t->w * t->h); // Normalized sum of squares of the image for (int y = v; y < (v + t->h); y++) { for (int x = u; x < (u + t->w); x++) { int a = (int) f->data[y * f->w + x] - f_mean; int b = (int) t->data[(y - v) * t->w + (x - u)] - t_mean; num += a * b; } } uint32_t den_a = f_sumsq - f_sum * (f_sum / (float) (t->w * t->h)); // Find normalized cross-correlation float c = num / (fast_sqrtf(den_a) * fast_sqrtf(den_b)); if (c > corr) { corr = c; r->x = u; r->y = v; r->w = t->w; r->h = t->h; } } } imlib_integral_image_free(&sum); imlib_integral_image_free(&sumsq); return corr; }帮我找到我想要的函数
最新发布
07-25
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值