图像匹配是计算机视觉和模式识别中一项重要而广泛的应用，利用相关的匹配算法在两幅或者多福图像之间识别同一特征点。随着图像匹配技术的逐步发展，国内外的学者都对其进行了深入的研究，并提出了许多图像匹配算法^[1-2]。经典的匹配算法有：归一化互相关(normalized cross-correlation，NCC)匹配算法、差的绝对值和(sum of absolute difference，SAD)相关算法、差的平方和(sum of square differences，SSD)相关算法等。其中NCC算法在现有相关匹配算法中应用最为广泛，具有较好的鲁棒性，对光照强度的线性变化不太敏感，抗白噪声干扰能力强，但其缺点是计算量大、匹配速度慢、错误匹配较多，并不符合系统实时性的要求^[3-5]。

因此如何降低计算的复杂度，减少运算过程中的计算量，提高匹配精度，已成为现在研究的难点和热点问题。Lewis提出用两个加和表的方法来来简化NCC分母的计算量，但没有改变分子的运算复杂度^[6]；Tsai等通过改进了NCC的定义，构造三个和表法来同时简化分子和分母的计算量，但匹配精度并不高^[7]；汪华琴针对归一化互相关算法精度不高，采用分层匹配来提高搜索匹配的效率，但计算仍很复杂^[8]。孙祖鑫等采用递推与多模板思想构建归一化互相关快速算法，但匹配率并不是很理想^[9]。针对以上问题，本文首先对传统的归一化互相关算法进行改进，减少计算量，同时在搜索策略上利用小波金字塔分层匹配的思想，分辨率由低到高，匹配由粗到细，提出一种基于小波金字塔搜索策略的快速NCC图像匹配算法。

1 快速归一化相关算法 1.1 归一化互相关匹配算法

归一化互相关算法(NCC)是一种常用的图像特征点匹配算法，其原理是根据两幅图像中特征点邻域像素灰度值的相似性来匹配的，对于左图像中的一点，计算其与右图像中所有特征点的归一化互相关系数，当得到其中最大值的点就为最佳匹配位置。

设W₁和W₂分别是图像I₁和图像I₂的两个大小相同的匹配窗口，图像大小为M×N，匹配窗口大小为m×n，u₁和u₂是匹配窗口内像素灰度的均值。归一化相关系数的定义形式通常有去均值和不去均值两种^[10-11]：

不去均值：

$r\left( x,y \right)=\frac{\sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{{{W}_{1}}}\left( x+i,y+j \right){{W}_{2}}\left( x+i,y+j \right)}}{\sqrt{\sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{{{W}_{1}}}{{\left( x+i,y+j \right)}^{2}}}}\sqrt{\sum\limits_{i=0}^{m-1}{\sum\limits_{j=0}^{n-1}{{{W}_{2}}}{{\left( x+i,y+j \right)}^{2}}}}}$

(1)

$r\left( x,y \right)=\frac{\sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{\left( {{W}_{1}}\left( x+i,y+j \right)-{{u}_{1}} \right)}\left( {{W}_{2}}\left( x+i,y+j \right)-{{u}_{2}} \right)}}{\sqrt{\sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{{{\left( {{W}_{1}}\left( x+i,y+j \right)-{{u}_{1}} \right)}^{2}}}}}\sqrt{\sum\limits_{i=0}^{m-1}{\sum\limits_{j=0}^{n-1}{{{\left( {{W}_{2}}\left( x+i,y+j \right)-{{u}_{2}} \right)}^{2}}}}}}$

(2)

式中：(x，y)∈M×N；r(x，y)的取值范围为－1~1，值越大表示相关程度越高。

虽然去均值的归一化互相关图像匹配算法的计算量比较大，但是它对于灰度和较小的几何畸变存在不变性，同时具有很强的抗噪声干扰的能力，因此本文选用去均值的归一化互相关匹配算法。在本文中，对两幅同样大小的图像，使用Harris算法进行特征点检测，传统的归一化互相关方法行进特征点匹配，但匹配结果精度不高，且所用时间长。本文在现有算法的基础上，针对NCC算法，通过构建3个加和表的方法，减少分子分母计算量，并引入图像小波金字塔分层结构，提高搜索匹配的效率。

1.2 快速归一化互相关算法

传统的归一化互相关匹配算法进行特征点匹配时，需要对特征点进行遍历搜索，共需要3mnMN加运算和2mnMN乘运算，计算量非常庞大。因此在本文中通过构建加和表^[12]来减少匹配过程中计算量，加快运算速度。

快速NCC匹配算法在计算归一化互相关系数时，将公式变换为

$r\left( x,y \right)=\frac{\sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{{{W}_{1}}\left( x+i,y+j \right){{W}_{2}}\left( x+i,y+j \right)-{{u}_{2}}}-mn{{u}_{1}}{{u}_{2}}}}{\sqrt{\sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{{{W}_{1}}{{\left( x+i,y+j \right)}^{2}}-mnu_{1}^{2}\times }}}\sqrt{\sum\limits_{i=0}^{m-1}{\sum\limits_{j=0}^{n-1}{{{W}_{2}}{{\left( x+i,y+j \right)}^{2}}-mnu_{2}^{2}}}}}$

(3)

其中，u₁=$\frac{1}{mn}\sum\limits_{i=0}^{m\text{-}1}{\sum\limits_{j=0}^{n\text{-}1}{{{W}_{1}}}}$(x+i，y+j)；u₂=$\frac{1}{mn}\sum\limits_{i=0}^{m\text{-}1}{\sum\limits_{j=0}^{n\text{-}1}{{{W}_{2}}\left( x+i,y+j \right)}}$。

从式(3) 可以看出来，方程的归一化互相关操作涉及源图像W(x，y)和模板图像T(x，y)的平均值，平方以及W和T之间的互相关的计算。

在本文中，构建的和表法计算图像均值S₁(x，y)，图像方差S₂(x，y)和图像间的相关S₃(x，y)为^[13]

$\begin{align} &{{S}_{1}}\left( x,y \right)=W\left( x,y \right)+{{S}_{1}}\left( x-1,y \right)+ \\ &\quad \quad \quad \quad {{S}_{1}}\left( x,y-1 \right)-{{S}_{1}}\left( x-1,y-1 \right) \\ \end{align}$

(4)

$\begin{align} &{{S}_{2}}\left( x,y \right)=W{{\left( x,y \right)}^{2}}+{{S}_{2}}\left( x-1,y \right)+ \\ &\quad \quad \quad \quad {{S}_{2}}\left( x,y-1 \right)-{{S}_{2}}\left( x-1,y-1 \right) \\ \end{align}$

(5)

$\begin{align} &{{S}_{3}}\left( x,y \right)=W\left( x,y \right)\times T\left( x,y \right)+{{S}_{3}}\left( x-1,y \right)+ \\ &\quad \quad \quad \quad {{S}_{3}}\left( x,y-1 \right)-{{S}_{3}}\left( x-1,y-1 \right) \\ \end{align}$

(6)

当x, y＜0时，S₁(x, y)=S₂(x, y)=0。则W在x, y位置上与模板同样大小的区域的累积和与平方和为

$\begin{align} & \sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{W\left( x+i,y+j \right)=}}{{S}_{1}}\left( x+m-1,y+n-1 \right)- \\ & \quad \quad \quad \quad {{S}_{1}}\left( x-1,y+n-1 \right)-\rm{ }{{\mathit{S}}_{1}}\left( \mathit{x}+\mathit{m}-1,\mathit{y}-1 \right)+ \\ & \quad \quad \quad \quad {{S}_{1}}\left( x-1,y-1 \right) \\ \end{align}$

(7)

$\begin{align} & \sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{W{{\left( x+i,y+j \right)}^{2}}=}}{{S}_{2}}\left( x+m-1,y+n-1 \right)- \\ & \quad \quad \quad \quad {{S}_{2}}\left( x-1,y+n-1 \right)-\rm{ }{{\mathit{S}}_{2}}\left( \mathit{x}+\mathit{m}-1,\mathit{y}-1 \right)+ \\ & \quad \quad \quad \quad {{S}_{2}}\left( x-1,y-1 \right) \\ \end{align}$

(8)

$\begin{align} & \sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{W\left( x+i,y+j \right)\times =}}T\left( x+i,y+j \right)= \\ & \quad \quad \quad \quad {{S}_{3}}\left( x+m-1,y+n-1 \right)-\rm{ }{{\mathit{S}}_{3}}\left( \mathit{x}-1,\mathit{y+n}-1 \right)- \\ & \quad \quad \quad \quad {{S}_{3}}\left( x+m-1,y-1 \right)+{{S}_{3}}\left( x-1,y-1 \right) \\ \end{align}$

(9)

将式(7)~(9) 代入式(3)，令δ_W2=$\sqrt{\sum\limits_{i=0}^{m\rm{-}1}{\sum\limits_{j=0}^{n\rm{-}1}{{{W}_{2}}{{\left( x+i,y+j \right)}^{2}}}}-m\cdot n\cdot u_{2}^{2}}$得，

$\begin{array}{l} r\left( {x,y} \right) = \\ \frac{{\left( \begin{array}{l} {S_3}\left( {x + m - 1,y + n - 1} \right) - {\rm{ }}{\mathit{S}_3}\left( {\mathit{x} - 1,\mathit{y + n} - 1} \right) - \\ {S_3}\left( {x + m - 1,y - 1} \right) + {S_3}\left( {x - 1,y - 1} \right) \end{array} \right) - \left( \begin{array}{l} {S_1}\left( {x + m - 1,y + n - 1} \right) - {\rm{ }}{\mathit{S}_1}\left( {\mathit{x} - 1,\mathit{y + n} - 1} \right) - \\ {S_1}\left( {x + m - 1,y - 1} \right) + {S_1}\left( {x - 1,y - 1} \right) \end{array} \right)\bar T}}{{\sqrt {\left( \begin{array}{l} {S_3}\left( {x + m - 1,y + n - 1} \right) - {\rm{ }}{\mathit{S}_3}\left( {\mathit{x} - 1,\mathit{y + n} - 1} \right) - \\ {S_3}\left( {x + m - 1,y - 1} \right) + {S_3}\left( {x - 1,y - 1} \right) \end{array} \right){\delta _{{W_2}}}} }} \end{array}$

(10)

对于大小为M×N的两个比较图像和大小为m×n的邻域窗口，提出的和表法计算复杂度从ο(m×n×M×N)显著降为ο(M×N)，只需要18MN加运算和2MN乘运算。

2 基于小波金字塔的图像匹配算法 2.1 小波金字塔搜索策略

金字塔结构能够减少图像匹配搜索的时间，而小波分解变换具有多分辨率的优点，能够有效地保留图像中的大部分信息，实现对图像的自适应分析。因此本文采用小波金字塔分层结构进行搜索。

其思想是^[14-15]：首先用L表示低通滤波器，H表示高通滤波器，分别对图像的行列进行卷积，然后进行2取1的亚抽样，则离散小波变换可将原图像分解成4个子带，即LL₁，LH₁，HL₁，HH₁。其中LL₁反应图像的低频成分，是由水平和垂直两个方向的低通滤波器获得的子带，LH₁则反应图像的水平边缘细节，是由水平方向低通滤波器和垂直方向的高通滤波器获得的子带，类似地，HL₁为水平方向高频和垂直方向低频获得的子带，HH₁则是由水平方向高频和垂直方向高频获得的子带。其次将分辨率设为原图的1/2，对低频分量LL₁进行再一步分解，又能够获得LL₂、LH₂、HL₂和HH₂共4个子带。按照这个过程反复，能够实现图像的多级分解。则小波变换三级分解如图 1所示。

考虑到对图像原有信息的保护以及算法的稳定性，分层的层数一般选择3~5层。在本文中，选用3层金字塔分层的结构进行匹配。其原理是^[16]：利用小波变换的理论对原始图像进行逐级分解，得到一个尺寸规模由小到大、分辨率从低到高的金字塔分层结构；匹配从小波分解的最低分辨率即金字塔顶层开始，利用匹配算法首先确定两幅图像上粗匹配的大致位置(x₁，y₁)和(x₂，y₂)，然后根据子带树型关系，在下一层映射位置中的中心点(2x₁，2y₁)和(2x₂，2y₂)邻域内找到精确匹配点，最终在金字塔底层图像上得到满足匹配精度要求的特征点。

在本文中，邻域的选取影响匹配的准确度，邻域选择过小时，包含信息少，邻域选择过大时，包含特征点多，都易造成误匹配，因此本文选择7×7的邻域进行匹配，小波金字塔搜索示意图如图 2所示^[17]。

2.2 本文算法具体步骤

采用小波金字塔搜索策略，每一层的结果都是以前一层搜索结果作为约束，与直接用原图像进行匹配相比，缩小了搜索的范围，因此减少运算中的计算量，并且提高了匹配的准确度。该算法的具体步骤如下^[18-19]：

1) 对左右两幅图像用小波变换进行2层分解，构成三层图像金字塔。在本文中，采用参数设置相同的Harris特征点检测算法对每层金字塔图像进行特征点提取，获得图像的特征点金字塔。

2) 因为低频分量图像中集中了原始图像的大部分能量，所以匹配从顶层金字塔图像的低频分量图像开始，利用改进后的快速NCC算法进行匹配，得到该层的最佳匹配点。

3) 将上层的匹配点作为下层图像匹配的中心点，在左右两幅图像的中心点的邻域内重新搜索进行归一化互相关匹配计算，得到本层的最佳匹配点。

4) 重复进行3)，随着分辨率的提高，互相关匹配的搜索范围被限定在一个比较小的范围，随着分辨率的提高，匹配点对的精度逐渐提高。经过匹配得到原始图像上的最佳匹配点。根据所得到的最终结果，将待匹配的准图像进行相应处理，完成匹配过程。

本文算法的框架图为图 3。

3 实验结果与分析

为了便于分析，采用本文算法对两组图像进行匹配，并与传统NCC算法和文献[9]匹配算法进行比较。对于A组图像进行详细的分析，实验采用的左右两幅原始图像大小都为320×400，经小波金字塔分层处理后，中层图像大小为160×200，顶层为80×100，在金字塔分层搜索过程中，归一化互相关邻域大小统一设置为7×7。采用传统NCC算法对两幅图像进行匹配，结果如图 4(a)所示；采用文献[9]匹配算法获得的结果图如图 4(b)所示, 采用本文算法得到的原始图像的匹配点连线如图 4(c)所示。对于B组图像，大小为224×344，采用传统NCC算法，文献[9]算法和本文算法进行匹配，结果图如图 5(a)、5(b)、5(c)所示。

由实验结果对比图可以看出，传统的NCC算法特征点过多，匹配不精确，而文献[9]匹配算法虽然能减少特征点的数目，但误匹配率仍然很高。本文算法采用由粗到精的匹配模式，随着分辨率的提高，匹配点对的精度逐渐提高，从而产生更多匹配对，误匹配连线较少，提高了匹配精度。

传统NCC算法、文献[9]匹配算法与本文算法采用时间结果如表 1所示。

从表 1可以看出，本文的算法虽然增加了图像特征点检测的时间，由于采用小波金字塔分层搜索策略，缩小搜索范围，以前一层匹配的结果作为下一层匹配的约束，在上层匹配点的邻域内搜索匹配，不用遍历整幅图像，大大地减少了匹配的时间。

4 结论

本文通过构造三个加和表的方法较少计算过程中乘法的计算量，降低运算的复杂度，同时采用小波三级分解获得的图像金字塔，根据由粗到精的匹配模式，使得特征点定位更加精确。对于A、B两组不同的实验图像，与采用传统的算法和文献[9]算法相比，该算法虽然在进行图像特征点检测阶段增加了1.28 s和1.43 s，但是在特征点匹配阶段也大大缩短了时间，因此总体时间上分别减少了3.5 s、0.47 s和3.52 s、0.3 s；同时该算法获得的误匹配较少，匹配连线效果更好，具有实际应用价值。

文中针对静止的图像进行的匹配，后续可以针对改进的算法对运动目标的匹配测试的准确度和实时性进行进一步研究。