本文翻译自论文：/publication/2378357_Fast_Normalized_Cross-Correlation
摘要
虽然众所周知，互相关可以在变换域中有效实现，但特征匹配应用程序首选的归一化形式的互相关没有简单的频域表达式。为此，在空间域中计算了归一化互相关。这篇短文表明，使用搜索窗口上的图像和图像 2 ^2 2的预计算积分，可以有效地归一化非规范化互相关。
1 介绍
两个信号之间的相关性（互相关）是特征检测的标准方法，也是更复杂技术的组成部分。TextBook presentations of correlations描述卷积定理以及使用快速傅立叶变换在频域有效计算相关性的可能性。不幸的是，模板匹配中首选的归一化相关形式（相关系数）没有相应的简单有效的频域表达式。因此，在空间域中计算了归一化互相关。由于空间域卷积的计算量大，还开发了几种不精确但快速的空间域匹配方法。本文描述了最近引入的一种从变换域卷积中获得归一化互相关的算法。在某些情况下，新算法在归一化互相关的空域计算上提供了一个数量级的加速（第5节）。
由于我们正在介绍一种熟悉且广泛使用的算法，因此不会试图调查有关特征选择、白化、快速卷积技术、扩展、替代技术或应用的文献。有关这些主题的文献可以通过介绍性文本和手册来获取。然而，由于所提倡的特征跟踪方案的多样性，可能有必要确定归一化互相关对于某些应用（如果不是所有应用）仍然是可行的选择。这在第3节中完成。为了使本文内容完整，第2节介绍了归一化互相关，第4节简要回顾了变换域和其他快速卷积方法以及相位相关技术。大多数读者都可以跳过这些部分。第5节描述了如何从相关的变换域计算中获得归一化互相关。第6节介绍了性能结果。
2 互相关模板匹配
互相关用于模板匹配的动机是距离度量（平方欧几里德距离）：
d f , t 2 ( u , v ) = ∑ x , y [ f ( x , y ) − t ( x − u , y − v ) ] 2 d_{f,t}^{2}(u,v) = \sum_{x,y}[f(x,y)-t(x-u,y-v)]^2 df,t2(u,v)=x,y∑[f(x,y)−t(x−u,y−v)]2
f f f是图像，而位于 u u u, v v v处的特征 t t t的窗口与图像逐个像素进行作差取平方再求和得到 d 2 d^2 d2
将 d 2 d^2 d2进行展开得到，
d f , t 2 ( u , v ) = ∑ x , y [ f 2 ( x , y ) − 2 f ( x , y ) t ( x − u , y − v ) + t 2 ( x − u , y − v ) ] 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)] df,t2(u,v)=x,y∑[f2(x,y)−2f(x,y)t(x−u,y−v)+t2(x−u,y−v)]
其中， ∑ t 2 ( x − u , y − v ) \sum{t^2(x-u,y-v)} ∑t2(x−u,y−v)是常量。那么若 ∑ f 2 ( x , y ) \sum{f^2(x,y)} ∑f2(x,y)近似为常量，则剩余的互相关项为：
c ( u , v ) = ∑ x , y f ( x , y ) t ( x − u , y − v ) c(u,v) = \sum_{x,y}{f(x,y)t(x-u,y-v)} c(u,v)=x,y∑f(x,y)t(x−u,y−v)
于是，该公式(1)为图像和特征之间的相似性度量
但是，使用该公式来模板匹配有几个缺点：
如果图像能量 ∑ f 2 ( x , y ) \sum{f^2(x,y)} ∑f2(x,y)随位置变化，则使用进行匹配可能会失败。例如，特征与图像中精确匹配区域之间的相关性可能小于特征与亮点之间的相关性。 c ( u , v ) c(u,v) c(u,v)的范围取决于特征的大小对于图像振幅的变化（如图像序列中的照明条件变化引起的变化）不是不变的
相关系数通过将图像和特征向量归一化为单位长度来克服这些困难，从而产生类余弦相关系数，公式(2)：
γ ( u , v ) = ∑ x , y [ f ( x , y ) − f ˉ u , v ] [ t ( x − u , y − v ) − t ˉ ] { ∑ x , y [ f ( x , y ) − f ˉ u , v ] 2 [ ∑ x , y t ( x − u , y − v ) − t ˉ ] 2 } 0.5 \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)-\bar{t}]^2}\}^{0.5}} γ(u,v)={∑x,y[f(x,y)−fˉu,v]2[∑x,yt(x−u,y−v)−tˉ]2}0.5∑x,y[f(x,y)−fˉu,v][t(x−u,y−v)−tˉ]
其中， t ˉ \bar{t} tˉ是特征的平均值， f ˉ u , v \bar{f}_{u,v} fˉu,v是特征下区域 f ( x , y ) f(x,y) f(x,y)的平均值
3 特征跟踪方法和问题
很明显，归一化互相关（NCC）并不是特征跟踪的理想方法，因为它在成像比例、旋转和透视失真方面并不是不变的。这些限制已在各种方案中得到解决，包括一些将NCC作为一个组件合并的方案。本文不提倡选择NCC的替代方法。相反，下面的讨论将指出各种特征跟踪方法所涉及的一些问题，并得出结论，对于某些应用来说，NCC是一个合理的选择。
SSDA.序列相似性检测算法 (SSDA) 的基础是观察到，只有在互相关函数的最大值附近才需要全精度，而降低的精度可以在其他地方使用。作者描述了实现“降低精度”的几种方法。SSDA的互相关实现通过以随机顺序计算公式(1)中的和来进行，并使用部分计算作为特定匹配位置是否接近相关曲面最大值的蒙特卡罗估计。如果估计结果表明某个位置对应的匹配不好，则在完成求和之前终止特定位置的计算。
SSDA算法简单，在空域互相关上提供了显著的加速。它的缺点是不能保证找到相关曲面的最大值。当相关面具有较浅的斜率和较宽的最大值时，SSDA表现良好。虽然在许多应用中可能满足此条件，但很明显，包含阵列对象（卵石、砖块、其他纹理）的图像会在相关曲面中生成多个窄极值，从而误导SSDA方法。SSDA的第二个缺点是，它具有需要确定的参数（用于形成相关系数估计的术语数量，以及该估计的提前终止阈值）。
Gradient Descent Search.如果假设相邻帧之间的特征平移很小，则可以通过梯度下降获得平移（和仿射扭曲的参数）。成功的梯度下降搜索要求帧间平移小于匹配误差曲面最小值周围的盆地半径。这一条件在许多应用中都可以满足。然而，手持摄像机的图像序列违反了要求：摄像机的小旋转会导致较大的物体平移。小的或（与SSDA一样）纹理模板会导致使用窄极值匹配错误曲面，从而限制可以成功跟踪的帧间平移范围。梯度下降技术的另一个缺点是搜索本质上是串行的，而NCC允许并行实现。
Snakes.Snakes（活动轮廓模型）的缺点是无法跟踪没有可定义轮廓的对象。一些“对象”没有明确定义的边界（无论是由于固有的模糊性还是由于照明条件），但仍然具有可通过互相关跟踪的颜色特征分布。活动轮廓模型解决了比简单模板匹配更普遍的问题，因为它们提供了随时间变化的变形轮廓的表示。互相关可以跟踪随时间变形的对象，但具有明显和重要的限制条件，此处将不讨论这些限制条件。互相关还可以很容易地跟踪一个特征，该特征在帧之间移动了其自身大小的一个显著分形，而这种平移量可能会使蛇超出其收敛范围。
Wavelets and other multi-resolution schemes.虽然对于小波来说，一个有用的卷积定理的存在仍然是一个讨论的问题，在某些方案中，小波卷积实际上是使用傅里叶卷积定理实现的，有效的特征跟踪可以使用小波和其他多分辨率表示（使用从粗到细的多分辨率搜索）来实现。然而，多分辨率技术要求图像包含足够的低频信息，以指导搜索的初始阶段。如第6节所述，理想特征有时不可用，人们必须求助于定义不清的“特征”，这些特征可能几乎没有低频信息，例如均匀表面上的小点配置。
以上讨论的每种方法都是由不同的作者提倡的，但不同方法之间的比较较少。参考文献[19]在梯度搜索框架内导出了一个最佳特征跟踪方案，但没有解决该框架的局限性。在存在各种图像失真的情况下，对五种模板匹配算法的实证研究发现，NCC在所有图像类别中都提供了最好的性能，尽管其中一种更便宜的算法在某些类型的失真中表现几乎一样好。文献[1]讨论了一种用于运动跟踪的通用层次结构框架。虽然也考虑了梯度方法，但选择了基于相关的匹配方法。
尽管 NCC 算法的时代已经到来，并且存在解决其各种缺点的最新技术，但可以公平地说，合适的替代方案尚未得到普遍认可。NCC对图像序列的要求很少，没有用户需要搜索的参数。NCC可以“按原样”提供简单的特征跟踪，也可以作为更复杂（可能是多分辨率）匹配方案的一个组成部分，该方案可以解决缩放和旋转不变性、特征更新和其他问题。选择相关系数而不是替代匹配标准（例如绝对差之和）也被证明是最大似然估计。我们承认 NCC 是许多应用程序的默认选择，其中特征跟踪本身不是研究主题，以及视觉和模式识别研究中的偶尔构建块。因此，一种快速算法很有意义。
4 变换域计算
考虑公式(2)中的分子，假定我们有图像 f ′ ( x , y ) = f ( x , y ) − f ˉ u , v f'(x,y)=f(x,y)-\bar{f}_{u,v} f′(x,y)=f(x,y)−fˉu,v和特征 t ′ ( x , y ) = t ( x , y ) − t ˉ t'(x,y)=t(x,y)-\bar{t} t′(x,y)=t(x,y)−tˉ，其中均值可以看作被移除，公式(3)：
γ ( u , v ) = ∑ x , y f ′ ( x , y ) t ′ ( x − u , y − v ) \gamma(u,v) = \sum_{x,y}{f'(x,y)t'(x-u,y-v)} γ(u,v)=x,y∑f′(x,y)t′(x−u,y−v)
对于大小为 M 2 M^2 M2 的搜索窗口和大小为 N 2 N^2 N2 的特征大约需要 N 2 ( M − N + 1 ) 2 N^2(M − N + 1)^2 N2(M−N+1)2 次加法和 N 2 ( M − N + 1 ) 2 N^2(M − N + 1)^2 N2(M−N+1)2 次乘法。
公式(3)是图像与反转特征 t ′ ( − x , − y ) t'(-x,-y) t′(−x,−y)的卷积，可以由下式计算，公式(4)：
F − 1 { F ( f ′ ) F ∗ ( t ′ ) } F^{-1}\{F(f')F^{*}(t')\} F−1{F(f′)F∗(t′)}
其中 F F F是傅里叶变换。复共轭通过傅里叶变换特性实现特征的反转 F ∗ ( − x ) = F ∗ ( ω ) F^{*}(−x) =F^{*}(ω) F∗(−x)=F∗(ω)
FFT 算法的实现通常需要将 f ′ f' f′ 和 t ′ t' t′ 用零扩展为 2 的共同幂。变换计算 (3) 的复杂度是 12 M 2 l o g 2 M 12M^2log_2M 12M2log2M 实数乘法和 18 M 2 l o g 2 M 18M^2log_2M 18M2log2M 实数加法/减法。当 M M M远大于 N N N时，直接“空间”计算的复杂度约为 N 2 M 2 N^2M^2 N2M2乘法/加法，直接方法比变换方法快。当 N N N接近 M M M时或者两者都较大时，变换方法变得相对更有效。
4.1 快速卷积
有几种众所周知的“快速”卷积算法不使用变换域计算。这些方法分为两类：将乘法换成附加加法的算法，以及通过将一维卷积的部分嵌入到较小的多维卷积的单独维度中来找到（一维）卷积 O ( N 2 ) O(N^2) O(N2)特征的较低点的方法。虽然比直接卷积更快，但这些算法在中等大小下仍比变换域卷积慢，而且在任何情况下，它们都不涉及公式（2）分母的计算。
4.2 相位相关
由于公式（4）可以在变换域中有效地计算，因此开发了几种近似公式（2）中图像能量归一化的变换域方法。通过在互相关之前对图像进行高通滤波，可以减少模板下图像能量的变化。该滤波可以方便地添加到频域处理中，但是截止频率的选择是有问题的——低截止可能会留下显著的图像能量变化，而高截止可能会删除对匹配有用的信息。
一种更稳健的方法是相位相关法。在这种方法中，在频域中计算相关性之前，将变换系数归一化为单位幅值。因此，相关性仅基于相位信息，对图像强度的变化不敏感。虽然经验表明这种方法是成功的，但它有一个缺点，即所有转换组件的权重都是相等的，而人们可能会认为不重要的组件应该被赋予较少的权重。原则上应该选择频谱预滤波，以便在给定信号和信号噪声的预期二阶矩的情况下最大化预期的相关信噪比。这种方法在[16]中进行了讨论，类似于经典的匹配滤波随机信号处理技术。对于典型的 (ρ ≈ 0.95) 图像相关性，最好的预滤波是近似拉普拉斯滤波，而不是纯白化滤波。
5 归一化
再次检查公式（2）的分子，我们注意到特征的平均值可以预先计算，留下
γ ( u , v ) = ∑ x , y f ′ ( x , y ) t ′ ( x − u , y − v ) − f ˉ u , v ∑ x , y t ′ ( x − u , y − v ) \gamma(u,v) = \sum_{x,y}{f'(x,y)t'(x-u,y-v)}-\bar{f}_{u,v}\sum_{x,y}{t'(x-u,y-v)} γ(u,v)=x,y∑f′(x,y)t′(x−u,y−v)−fˉu,vx,y∑t′(x−u,y−v)
由于 t ′ t' t′的平均值为零，因此术语 f ˉ u , v ∑ x , y t ′ ( x − u , y − v ) \bar{f}_{u,v}\sum_{x,y}{t'(x-u,y-v)} fˉu,v∑x,yt′(x−u,y−v)也为零，因此可以使用公式（4）计算归一化互相关的分子。通过检查公式（2）的分母，可以在大约 3 N 2 3N^2 3N2操作中预计算特征向量的长度（与互相关的成本相比较小），事实上，可以将特征预规范化为长度1。
有问题的量是表达式中的 ∑ x , y [ f ( x , y ) − f ˉ u , v ] 2 \sum_{x,y}[f(x,y)-\bar{f}_{u,v}]^{2} ∑x,y[f(x,y)−fˉu,v]2。必须在每个u，v，即在 ( M − N + 1 ) 2 (M-N+1)^2 (M−N+1)2处计算图像均值和本地能量，导致近似 3 n 2 ( M − N + 1 ) 2 3n^2(M−N+1)^2 3n2(M−N+1)2个运算（加法、减法、乘法各算一个运算）。这种计算比（3）的直接计算所需的要多，当变换方法适用时，它可能大大超过（4）所示的计算。因此，需要一种更有效的方法来计算特征下的图像平均值和能量。
这些量可以从包含搜索区域上图像和图像平方的积分（运行和）的表格中有效地计算出来，例如：
s ( u , v ) = f ( u , v ) + s ( u − 1 , v ) + s ( u , v − 1 ) − s ( u − 1 , v − 1 ) s(u,v) = f(u,v)+s(u-1,v) +s(u,v-1) -s(u-1,v-1) s(u,v)=f(u,v)+s(u−1,v)+s(u,v−1)−s(u−1,v−1)
s 2 ( u , v ) = f 2 ( u , v ) + s 2 ( u − 1 , v ) + s 2 ( u , v − 1 ) − s 2 ( u − 1 , v − 1 ) s^2(u,v) = f^2(u,v)+s^2(u-1,v) +s^2(u,v-1) -s^2(u-1,v-1) s2(u,v)=f2(u,v)+s2(u−1,v)+s2(u,v−1)−s2(u−1,v−1)
当 u u u或 v v v小于0时， s ( u , v ) = s 2 ( u , v ) = 0 s(u,v)=s^2(u,v) =0 s(u,v)=s2(u,v)=0。然后，位于u，v的特征下的图像能量为
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 ) 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) ef(u,v)=s2(u+N−1,v+N−1)−s2(u−1,v+N−1)−s2(u+N−1,v−1)+s2(u−1,v−1)
对于特征下的图像和是一样的。
问题量 ∑ x , y [ f ( x , y ) − f ˉ u , v ] 2 \sum_{x,y}[f(x,y)-\bar{f}_{u,v}]^{2} ∑x,y[f(x,y)−fˉu,v]2现在只需很少的运算就可以计算，因为它扩展为一个表达式，只涉及图像和以及特征下的和平方。表的构造需要大约 3 M 2 3M^2 3M2的运算，这比计算分子的成本（4）要小得多，并且大大小于需要计算在每个 u , v u,v u,v计算 ∑ x , y [ f ( x , y ) − f ˉ u , v ] 2 \sum_{x,y}[f(x,y)-\bar{f}_{u,v}]^{2} ∑x,y[f(x,y)−fˉu,v]2导致的 3 n 2 ( M − N + 1 ) 2 3n^2(M−N+1)^2 3n2(M−N+1)2。
这种从预先计算的运行和计算确定和的技术已在许多领域中独立使用；[5]中开发了一个计算机图形学应用程序。如果以系统行扫描顺序搜索相关曲面的最大值，则可以通过状态变量组合表构造和引用，从而避免显式存储表。然而，当在通用计算机上实现时，表的大小不是主要考虑因素，搜索相关曲面的灵活性可能是有利的。请注意， s ( u , v ) s(u,v) s(u,v)和 s 2 ( u , v ) s^2(u,v) s2(u,v)表达式略微稳定，这意味着它们的 z z z变换 H ( z ） = 1 / ( 1 − Z − 1 ) H(z）=1/(1− Z^{−1}) H(z）=1/(1−Z−1)（此处为一维版本）在 z = 1 z=1 z=1处有一个极点，而稳定性要求极点严格位于单位圆内。因此，计算应该使用大整数而不是浮点运算。
