CN106447704A - 基于显著区域特征和边缘度的可见光‑红外图像配准方法 - Google Patents

Abstract

本发明公开了一种基于显著区域特征和边缘度的可见光‑红外图像配准方法，用于解决现有可见光‑红外图像配准方法提取的图像间特征相似度低导致配准精度差的技术问题。技术方案是采用基于自适应局部对比度的显著性区域检测方法从待配准图像对中提取显著区域，利用Zernike旋转不变矩描述显著区域形状特征，保证了所利用特征的平移、灰度、尺度和旋转不变性，通过RANSAC求取同名点对估计得到初始变换参数；最后利用边缘度构造的三维联合直方图对初始变换参数进行优化得到最终结果。该方法通过异源图像间相似度较高的显著区域和边缘度特征及Zernike矩的不变性，提高了图像配准的精度。

Description

基于显著区域特征和边缘度的可见光-红外图像配准方法

技术领域

本发明涉及一种可见光-红外图像配准方法，特别是涉及一种基于显著区域特征和边缘度的可见光-红外图像配准方法。

背景技术

多源图像配准技术是指将不同类型的传感器针对同一场景，从不同时间、从不同角度所获得的两幅或多幅图像数据进行预处理后，采用一定的算法、模型将它们的互补性信息或优势信息结合起来，变换到同一坐标系下的过程。

将多源图像中的数据进行可选择性的有机结合，能够更加全面真实的反映出事物本质信息和特征的多源融合技术的需求日益增多。多源图像配准为多源图像提供了统一的空间基准，是多源图像融合的基础，同时也可用于图像拼接、目标检测、运动估计、特征匹配等应用，具有非常重要的研究和实用价值。

可见光和红外图像间的边缘方向基本保持一致，且克服了梯度带来的差异。

文献1“J.H.Lee,Y.S.Kim,D.Lee,etal.Robust CCD and IR image registrationusing gradient-based statistical information,IEEE Signal Process.Lett.17(4)(2010)347–350”公开了一种基于边缘度构造三维联合直方图的可见光-红外图像配准方法。该方法首先计算图像边缘度和梯度信息，然后构建三维联合直方图，根据信息熵作为相似性测度进行优化得到图像变换矩阵。文献所述方法在平均目标配准误差(TargetRegistration Error)和配准成功率上明显优于其他方法，但由于构建三维联合直方图计算量较大，优化搜索空间大，受初始值影响容易陷入局部最优解，导致配准结果适应度不高。

人类的视觉容易被场景或图像中的较显著区域所吸引，对其他区域“视而不见”，称该区域具有显著性。

文献2“Cheng M M,Zhang G X,Mitra N J,et al.Global contrast basedsalient region detection[C].Computer Vision and Pattern Recognition(CVPR),2011IEEE Conference on.IEEE,2011:409-416”中分别基于全局直方图统计对比度和局部区域对比度进行了显著性区域的检测。其方法在正确率和召回率上明显优于其他方法，简单高效，适合于对图像进行显著性区域检测与分割。在对多源图像的显著性区域进行检测分割后，得到的相对应显著性区域的形状大体上相似。基于此，可利用Zernike不变矩对显著性区域的形状进行分析比较。1979年，根据正交多项式理论，M.R.Teague提出了Zernike矩，Zernike矩能够很容易地构造图像的任意高阶矩，并能够使用较少的矩来重建图像，且具有平移、灰度、尺度和旋转不变性。

发明内容

为了克服现有可见光-红外图像配准方法提取的图像间特征相似度低导致配准精度差的不足，本发明提供一种基于显著区域特征和边缘度的可见光-红外图像配准方法。该方法采用基于自适应局部对比度的显著性区域检测方法从待配准图像对中提取显著区域，利用Zernike旋转不变矩描述显著区域形状特征，保证了所利用特征的平移、灰度、尺度和旋转不变性，通过RANSAC求取同名点对估计得到初始变换参数；最后利用边缘度构造的三维联合直方图对初始变换参数进行优化得到最终结果。该方法通过异源图像间相似度较高的显著区域和边缘度特征及Zernike矩的不变性，提高了图像配准的精度。

本发明解决其技术问题所采用的技术方案是：一种基于显著区域特征和边缘度的可见光-红外图像配准方法，其特点是包括以下步骤：

步骤一、基于改进的局部对比度的显著性区域检测。

首先，使用基于图的图像分割算法将输入图像分割成若干区域，计算各区域的颜色直方图，则区域r_k显著性值S(r_k)的计算方法为式(1)。

式中，ω(r_i)为区域r_i的权值，D_r(r_k,r_i)为区域r_k和r_i之间的颜色距离度量，其计算方法见式(2)。

式中，f(c_k,i)为第i个颜色c_k,i在第k个区域r_k的所有颜色中出现的概率，该概率密度函数即为归一化的颜色直方图。D(c_1,i,c_2,j)为颜色c_1,i和颜色c_2,j的Lab空间颜色距离。

其次，考虑空间影响，将空间信息作为权值代入式(1)中，则区域r_k的显著性值S(r_k)为：

式中，D_S(r_k,r_i)为区域r_k和r_i的质心的欧氏距离；σ_S是控制空间权值大小的阈值，σ_S越小，空间权值影响越大，较远区域的对比度对当前区域的显著性值影响较小。由于可见光图像和红外图像分辨率不一致，因此采用自适应的空间权值σ_S，计算公式如式(4)：

式中，D_ang是图像对角线的长度，num是区域数量。

步骤二、基于Zernike矩的显著性区域形状特征描述。

Zernike矩是基于Zernike多项式构造出来的，m阶n次的Zernike多项式的定义范围为x²+y²≤1。其定义为式(5)。

V_mn(x,y)＝V_mn(ρ,θ)＝R_mn(ρ)e^jnθ (5)

式中，m为非负整数，m-|n|为偶数且m≥|n|；ρ为原点到点(x,y)的矢量长度；θ为x轴与ρ矢量逆时针方向夹角；R_mn(ρ)为点(x,y)的径向多项式，计算公式如下所示：

由此给出Zernike矩的定义为式(9)。

将Zernike矩的实部和虚部分别记为C_mn和S_mn，则：

将Zernike矩用于处理二维图像数据，因此需要对式(10)和式(11)进行离散化。引入两个新的参数(r,σ)，分别定义为：

r＝max(|x|,|y|) (12)

对于N×N的图像f(x,y)，以图像中心为坐标原点，则图像中的坐标值都满足-N/2≤(x、y)≤N/2。由式(12)和式(13)得知，r和σ的取值范围分别为1≤r≤N/2，1≤σ≤8r。则式(6)和式(7)分别改写为式(14)和式(15)。

θ＝4πσr (15)

经过新的变换，式(10)和式(11)可分别离散化为式(16)和式(17)。

综上所述，Zernike矩计算方法如下：

①确定图像f(x,y)的大小，从而确定式(16)和式(17)中N的值以及r和σ的取值范围；

②利用式(8)计算出R_mn(ρ)的值，并结合式(16)和式(17)计算出Zernike矩的实部C_mn和虚部S_mn；

③对C_mn和S_mn求模，得到Zernike矩的幅值|Z_mn|。

由于尺度差异问题，还使用区域r_k最小外接矩形的对角线l_k对|Z_mn|尺度归一化，计算方法如式(18)：

将Zernike矩用于显著性图中各显著性区域的形状分析，由于Zernike矩能够构造任意高阶，但考虑到计算复杂度，采用前36阶Zernike矩对区域特征进行描述。

综上所述，对区域r_k的Zernike矩描述子为：其中m∈{0,1,2,3,4,5,6,7}，m-|n|为偶数且m≥|n|。

步骤三、基于形状特征相似性度量和RANSAC的同名点求取和初始变换参数估计。

(1)同名点求取。

通过对图像各显著性区域计算质心和Zernike矩描述后，产生两组未匹配的质心点。在匹配时，正确同名点所对应的Zernike矩最具相似性，同时，由于得到的匹配点对中存在误匹配对，因此，利用RANSAC算法对匹配点对进行筛选，剔除错误的匹配点对。

设表示图像I_p中的第k个区域，参考图像为I₁，待配准图像为I₂，对于区域和区域其Zernike矩描述子分别为和则定义两区域的相似距离为式(19):

相似距离越小，则区域和区域越相似，匹配度越高。将相似距离最小的两个区域认为是匹配区域，其对应的质心为匹配点对。区域r_k的质心定义为式(20):

式中，(x,y)为区域r_k中的像素坐标，f(x,y)为像素值。即为区域r_k的质心坐标。由此得到一组匹配点对集合，但由于其中仍然包含误匹配对，因此接下来采用RANSAC算法将其中的错误匹配点对剔除。

根据误差度量函数把所有数据分为内点和外点。内点的集合又称为一致集，经过多次随机采样，当一致集达到最大时，则对应的模型即为所求结果。RANSAC算法的步骤如下：

①从候选点对集合S中随机抽取一个含有n个点对的样本集P，并估计初始化模型；

②对集合S中的剩余点对计算变换模型，与初始化模型M的距离小于误差阈值τ的点对标记为内点，形成一致集S*；

③重复上述步骤，直到达到最大采样次数N，或者集合S*的元素个数大于某一个阈值T。

④根据最大一致集S*，计算变换模型的参数。

RANSAC算法包含了四个输入参数：随机抽样数n、误差阈值τ、最大采样次数N和最优一致集阈值T。

1)随机抽样数n需要满足变换模型估计的最小配置要求，n取4组不共线的数据点对。

2)误差阈值τ对RANSAC性能有很大的影响，根据期望的内点概率来进行设置。τ为经验值，经验值取7，

3)最大采样次数N只要足够大，就能够保证匹配的正确性，N的取值计算公式如式(21)所示。

其中，n为随机抽样数；ω为数据局内点的概率；p表示迭代过程中从数据集内随机选取出的点均为局内点的概率。通常数据中局内点的概率ω是未知的，但能够通过自动更新ω来自适应决定所需要采样的次数。

4)较大的一致集S*使得变换模型参数更为精确，为了得到表征数据集S的正确模型，对T不做限制，最大一致集S*取局内点最多的一致集。

(2)初始变换参数估计。

获得筛选后的正确同名点对后，选择合适的变换模型并进行变换模型参数估计。对于透视投影变换模型，待配准图像I₂中匹配点(x',y')经其变换后得到配准后图像中点(x,y)的计算公式如式(22)所示，对应的参数矩阵H₀为式(23)：

式中，H₀为待配准图像I₂到参考图像I₁坐标系的变换矩阵，至少使用4对不共线同名点对计算出H₀中所有未知量的值。

步骤四、基于边缘特征的变换参数优化求解。

(1)边缘特征。

采用边缘度的概念。对于一个结构张量Q_ρ，其定义式如下式(24)：

式中，*是卷积算子，K_ρ表示一个标准差为ρ的二维高斯核，I_x和I_y分别表示x和y方向的梯度分量。通过在每个像素上对Q_ρ特征值分解，得到分别表示最大和最小的亮度变化的特征向量ω₁和ω₂，ω₁和ω₂分别对应两个特征值μ₁和μ₂，μ₁≥μ₂≥0。边缘度C的定义如下式(25)：

C的取值范围是[0,1]，ε为设定的阈值，用来防止同向区域像素取值过大。由此得到可见光-红外图像的边缘度图，它对边缘描述比常规的边缘检测方法差异更小，增加了可见光和红外图像边缘的相似性，提高了正确配准的概率。

边缘像素点(x,y)处的广义梯度向量流通过求取下式(26)的能量函数最小值得到：

式中，为像素点(x,y)处的梯度值，g(C)＝e^-C/K，C为边缘度，K为梯度场平滑程度的正则化因子，K取1。

(2)相似性度量。

当两幅图像配准时，边缘度的互信息越大，广义梯度向量流的方向误差不确定性越小，配准精度越高。因此，定义如下相似性度量函数(27)：

F(S,T)＝M(C_S,C_T)-H(D|C_S,C_T) (27)

式中，M(·)表示互信息，H(·)表示Shannon熵，C_S和C_T分别表示图像S和图像T的边缘度，D为C_S和C_T相应像素之间方向差异，其计算公式如下式(28)：

式中，和分别是图像S和图像T的广义梯度向量流。然后，以上一步估计得到的初始配准参数H₀为初值，相似性度量函数(27)为目标函数，采用单纯形法进行优化搜索，当F(S,T)达到最大值时，配准精度最高，以此时的H为最终变换矩阵。

最后，以H为变换矩阵对待配准图像进行图像变换，然后与参考图像加权融合，得到最终结果。

本发明的有益效果是：该方法采用基于自适应局部对比度的显著性区域检测方法从待配准图像对中提取显著区域，利用Zernike旋转不变矩描述显著区域形状特征，保证了所利用特征的平移、灰度、尺度和旋转不变性，通过RANSAC求取同名点对估计得到初始变换参数；最后利用边缘度构造的三维联合直方图对初始变换参数进行优化得到最终结果。该方法通过异源图像间相似度较高的显著区域和边缘度特征及Zernike矩的不变性，提高了图像配准的精度。由于利用基于自适应局部区域对比度的显著性区域检测方法，使得异源图像间的图像区域特征相似度更高。在对检测出的显著性区域进行形状分析时，利用了Zernike矩进行形状分析，保证了所利用特征的平移、灰度、尺度和旋转不变性。最后加入图像的边缘信息，使用边缘度对配准结果优化求解，使得配准结果更加精确。在真实的可见光-红外图像数据集上的试验表明，本发明所获得的可见光-红外图像配准精度达到了实际应用的需求。

下面结合具体实施方式对本发明作详细说明。

具体实施方式

本发明基于显著区域特征和边缘度的可见光-红外图像配准方法具体步骤如下：

1、基于改进的局部对比度的显著性区域检测。

首先，使用基于图的图像分割算法将输入图像分割成若干区域，计算各区域的颜色直方图，则区域r_k显著性值S(r_k)的计算方法为式(1)。

式中，ω(r_i)为区域r_i的权值，D_r(r_k,r_i)为区域r_k和r_i之间的颜色距离度量，其计算方法见式(2)。

式中，f(c_k,i)为第i个颜色c_k,i在第k个区域r_k的所有颜色中出现的概率，该概率密度函数即为归一化的颜色直方图。D(c_1,i,c_2,j)为颜色c_1,i和颜色c_2,j的Lab空间颜色距离。由于每个小区域中包含的颜色种类较少，故采用稀疏直方图来提高计算效率。

其次，考虑空间影响，将空间信息作为权值代入式(1)中，则区域r_k的显著性值S(r_k)为：

式中，D_S(r_k,r_i)为区域r_k和r_i的质心的欧氏距离；σ_S是控制空间权值大小的阈值，σ_S越小，空间权值影响越大，较远区域的对比度对当前区域的显著性值影响较小。由于可见光图像和红外图像分辨率不一致，因此采用自适应的空间权值σ_S，计算公式如下(4)：

式中，D_ang是图像对角线的长度，num是区域数量。

2、基于Zernike矩的显著性区域形状特征描述。

Zernike矩是基于Zernike多项式构造出来的，m阶n次的Zernike多项式的定义范围为x²+y²≤1。其定义为式(5)。

V_mn(x,y)＝V_mn(ρ,θ)＝R_mn(ρ)e^jnθ (5)

式中，m为非负整数，m-|n|为偶数且m≥|n|；ρ为原点到点(x,y)的矢量长度；θ为x轴与ρ矢量逆时针方向夹角；R_mn(ρ)为点(x,y)的径向多项式，它们的计算公式如下所示：

由此给出Zernike矩的定义为式(9)。

将Zernike矩的实部和虚部分别记为C_mn和S_mn，则：

将Zernike矩用于处理二维图像数据，因此需要对式(10)和式(11)进行离散化。引入两个新的参数(r,σ)，分别定义为：

r＝max(|x|,|y|) (12)

对于N×N的图像f(x,y)，以图像中心为坐标原点，则图像中的坐标值都满足-N/2≤(x、y)≤N/2。由式(12)和式(13)可知，r和σ的取值范围分别为1≤r≤N/2，1≤σ≤8r。则式(6)和式(7)可分别改写为式(14)和式(15)。

θ＝4πσr (15)

经过新的变换，式(10)和式(11)可分别离散化为式(16)和式(17)。

综上所述，Zernike矩计算方法如下：

①确定图像f(x,y)的大小，从而确定式(16)和式(17)中N的值以及r和σ的取值范围；

②利用式(8)计算出R_mn(ρ)的值，并结合式(16)和式(17)计算出Zernike矩的实部C_mn和虚部S_mn；

③对C_mn和S_mn求模，得到Zernike矩的幅值|Z_mn|。

由于尺度差异问题，还使用区域r_k最小外接矩形的对角线l_k对|Z_mn|尺度归一化，计算方法如式(18)：

本发明将Zernike矩用于显著性图中各显著性区域的形状分析，由于Zernike矩可以构造任意高阶，但考虑到计算复杂度，本发明采用前36阶Zernike矩对区域特征进行描述。

综上所述，本发明中对区域r_k的Zernike矩描述子为：其中m∈{0,1,2,3,4,5,6,7}，m-|n|为偶数且m≥|n|。

3、基于形状特征相似性度量和RANSAC的同名点求取和初始变换参数估计。

(1)同名点求取。

通过对图像各显著性区域计算质心和Zernike矩描述后，产生两组未匹配的质心点。在匹配时，正确同名点所对应的Zernike矩最具相似性，同时，由于得到的匹配点对中存在误匹配对，因此，利用RANSAC算法对匹配点对进行筛选，剔除错误的匹配点对。

设表示图像I_p中的第k个区域，参考图像为I₁，待配准图像为I₂，对于区域和区域其Zernike矩描述子分别为和则定义两区域的相似距离为式(19):

相似距离越小，则区域和区域越相似，匹配度越高。将相似距离最小的两个区域认为是匹配区域，其对应的质心为匹配点对。区域r_k的质心定义为式(20):

式中，(x,y)为区域r_k中的像素坐标，f(x,y)为像素值。即为区域r_k的质心坐标。由此得到一组匹配点对集合，但由于其中仍然包含误匹配对，因此接下来采用RANSAC算法将其中的错误匹配点对剔除。

RANSAC的基本思想是随机选择一组样本计算变换模型，并根据误差度量函数把所有数据分为内点(Inlier)和外点(Outlier)。内点的集合又称为一致集，经过多次随机采样，当一致集达到最大时，则对应的模型即为所求结果。RANSAC算法的基本步骤如下：

①从候选点对集合S中随机抽取一个含有n个点对的样本集P，并估计初始化模型；

②对集合S中的剩余点对计算变换模型，与初始化模型M的距离小于误差阈值τ(又称为误差容忍度)的点对标记为内点，形成一致集S*；

③重复上述步骤，直到达到最大采样次数N，或者集合S*的元素个数大于某一个阈值T。

④根据最大一致集S*，计算变换模型的参数。

根据上述的RANSAC算法基本步骤可以看出，该算法包含了四个输入参数：随机抽样数n、误差阈值τ、最大采样次数N和最优一致集阈值T。

1)随机抽样数n需要满足变换模型估计的最小配置要求，本发明中n取4组不共线的数据点对。

2)误差阈值τ对RANSAC性能有很大的影响，可以根据期望的内点概率来进行设置。τ为经验值，本发明中取7，

3)最大采样次数N只要足够大，就可以保证匹配的正确性，N的取值计算公式如式(21)所示。

其中，n为随机抽样数；ω为数据局内点的概率；p表示迭代过程中从数据集内随机选取出的点均为局内点的概率。通常数据中局内点的概率ω是未知的，但可以通过自动更新ω来自适应决定所需要采样的次数。

4)较大的一致集S*可以使得变换模型参数更为精确，为了得到表征数据集S的正确模型，本发明中对T不做限制，最大一致集S*取局内点最多的一致集。

(2)初始变换参数估计。

获得筛选后的正确同名点对后，需要选择合适的变换模型并进行变换模型参数估计。对于透视投影变换模型，待配准图像I₂中匹配点(x',y')经其变换后得到配准后图像中点(x,y)的计算公式如式(22)所示，对应的参数矩阵H₀为式(23)：

式中，H₀为待配准图像I₂到参考图像I₁坐标系的变换矩阵，至少使用4对不共线同名点对计算出H₀中所有未知量的值。

4、基于边缘特征的变换参数优化求解。

(1)边缘特征。

一般的边缘提取算法对于边缘对比度很敏感，为了解决这个问题，本发明采用边缘度的概念。对于一个结构张量Q_ρ，其定义式如下式(24)：

式中，“*”是卷积算子，K_ρ表示一个标准差为ρ的二维高斯核，I_x和I_y分别表示x和y方向的梯度分量。通过在每个像素上对Q_ρ特征值分解，得到分别表示最大和最小的亮度变化的特征向量ω₁和ω₂，它们分别对应了两个特征值μ₁和μ₂(μ₁≥μ₂≥0)。边缘度C的定义如下式(25)：

C的取值范围是[0,1]，ε为设定的阈值，用来防止同向区域像素取值过大。由此可得到可见光-红外图像的边缘度图，它对边缘描述比常规的边缘检测方法差异更小，增加了可见光和红外图像边缘的相似性，提高了正确配准的概率。

边缘像素点(x,y)处的广义梯度向量流可以通过求取下式(26)的能量函数最小值得到：

式中，为像素点(x,y)处的梯度值，g(C)＝e^-C/K，C为边缘度，K为梯度场平滑程度的正则化因子，本发明中K取1。

(2)相似性度量。

当两幅图像配准时，边缘度的互信息越大，广义梯度向量流的方向误差不确定性越小，配准精度越高。因此，定义如下相似性度量函数(27)：

F(S,T)＝M(C_S,C_T)-H(D|C_S,C_T) (27)

式中，M(·)表示互信息，H(·)表示Shannon熵，C_S和C_T分别表示图像S和图像T的边缘度，D为C_S和C_T相应像素之间方向差异，其计算公式如下式(28)：

式中，和分别是图像S和图像T的广义梯度向量流。然后，以上一步估计得到的初始配准参数H₀为初值，相似性度量函数(27)为目标函数，采用单纯形法进行优化搜索，当F(S,T)达到最大值时，配准精度最高，以此时的H为最终变换矩阵。

最后，以H为变换矩阵对待配准图像进行图像变换，然后与参考图像加权融合，得到最终结果。

Claims (1)

1.一种基于显著区域特征和边缘度的可见光-红外图像配准方法，其特征在于包括以下步骤：

步骤一、基于改进的局部对比度的显著性区域检测；

首先，使用基于图的图像分割算法将输入图像分割成若干区域，计算各区域的颜色直方图，则区域r_k显著性值S(r_k)的计算方法为式(1)；

$S\left(r_k\right)=\underset{r_k\≠r_i}{\Σ}\mathrm{\ω}\left(r_i\right)D_r(r_k,r_i)---\left(1\right)$

式中，ω(r_i)为区域r_i的权值，D_r(r_k,r_i)为区域r_k和r_i之间的颜色距离度量，其计算方法见式(2)；

$D_r(r_1,r_2)=\underset{i=1}{\overset{n_1}{\Σ}}\underset{j=1}{\overset{n_2}{\Σ}}f\left(c_{1,i}\right)f\left(c_{2,j}\right)D(c_{1,i},c_{2,j})---\left(2\right)$

式中，f(c_k,i)为第i个颜色c_k,i在第k个区域r_k的所有颜色中出现的概率，该概率密度函数即为归一化的颜色直方图；D(c_1,i,c_2,j)为颜色c_1,i和颜色c_2,j的Lab空间颜色距离；

其次，考虑空间影响，将空间信息作为权值代入式(1)中，则区域r_k的显著性值S(r_k)为：

$S\left(r_k\right)=\underset{r_k\≠r_i}{\Σ}\exp (-\frac{D_S(r_k,r_i)}{{\mathrm{\σ}}_S^2})\mathrm{\ω}\left(r_i\right)D_r(r_k,r_i)---\left(3\right)$

式中，D_S(r_k,r_i)为区域r_k和r_i的质心的欧氏距离；σ_S是控制空间权值大小的阈值，σ_S越小，空间权值影响越大，较远区域的对比度对当前区域的显著性值影响较小；由于可见光图像和红外图像分辨率不一致，因此采用自适应的空间权值σ_S，计算公式如式(4)：

$\begin{array}{c}{\mathrm{\σ}}_s=\stackrel{\‾}{D}/D_{ang}\\ \stackrel{\‾}{D}=\underset{r_k\≠r_i}{\Σ}{D}_s(r_k,r_i)/(num-1)\end{array}---\left(4\right)$

式中，D_ang是图像对角线的长度，num是区域数量；

步骤二、基于Zernike矩的显著性区域形状特征描述；

Zernike矩是基于Zernike多项式构造出来的，m阶n次的Zernike多项式的定义范围为x²+y²≤1；其定义为式(5)；

V_mn(x,y)＝V_mn(ρ,θ)＝R_mn(ρ)e^jnθ (5)

式中，m为非负整数，m-|n|为偶数且m≥|n|；ρ为原点到点(x,y)的矢量长度；θ为x轴与ρ矢量逆时针方向夹角；R_mn(ρ)为点(x,y)的径向多项式，计算公式如下所示：

$\begin{array}{cc}\mathrm{\&rho;}=\sqrt{x^2+y^2}& (-1<x,y<1)\end{array}---\left(6\right)$

$\begin{array}{cc}\mathrm{\&theta;}=arctan\left(\frac{y}{x}\right)& (-1<x,y<1)\end{array}---\left(7\right)$

$R_{mn}\left(\mathrm{\&rho;}\right)=\underset{s=0}{\overset{(m-|n\left|\right)/2}{\&Sigma;}}\frac{{(-1)}^s(m-s)!{\mathrm{\&rho;}}^{m-2s}}{s!(\frac{m+\left|n\right|}{2}-s)!(\frac{m+\left|n\right|}{2}-s)!}---\left(8\right)$

由此给出Zernike矩的定义为式(9)；

$\begin{array}{c}{Z}_{mn}=\frac{m+1}{\mathrm{\&pi;}}\underset{0}{\overset{1}{\&Integral;}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}\&lsqb;V_{mn}(\mathrm{\&rho;},\mathrm{\&theta;})\&rsqb;f(\mathrm{\&rho;},\mathrm{\&theta;})\mathrm{\&rho;}\frac{dy}{dx}d\mathrm{\&rho;}d\mathrm{\&theta;}\\ =\frac{m+1}{\mathrm{\&pi;}}\&Integral;\&Integral;R_{mn}\left(\mathrm{\&rho;}\right)e^{jn\mathrm{\&theta;}}f(\mathrm{\&rho;},\mathrm{\&theta;})d\mathrm{\&rho;}d\mathrm{\&theta;}\end{array}---\left(9\right)$

将Zernike矩的实部和虚部分别记为C_mn和S_mn，则：

$C_{mn}=\frac{2m+2}{\mathrm{\&pi;}}\underset{0}{\overset{1}{\&Integral;}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}{R}_{mn}\left(\mathrm{\&rho;}\right)cos\left(n\mathrm{\&theta;}\right)f(\mathrm{\&rho;},\mathrm{\&theta;})\mathrm{\&rho;}d\mathrm{\&rho;}d\mathrm{\&theta;}---\left(10\right)$

$S_{mn}=\frac{2m+2}{\mathrm{\&pi;}}\underset{0}{\overset{1}{\&Integral;}}\underset{0}{\overset{2\mathrm{\&pi;}}{\&Integral;}}{R}_{mn}\left(\mathrm{\&rho;}\right)sin\left(n\mathrm{\&theta;}\right)f(\mathrm{\&rho;},\mathrm{\&theta;})\mathrm{\&rho;}d\mathrm{\&rho;}d\mathrm{\&theta;}---\left(11\right)$

将Zernike矩用于处理二维图像数据，因此需要对式(10)和式(11)进行离散化；引入两个新的参数(r,σ)，分别定义为：

r＝max(|x|,|y|) (12)

$\left\{\begin{array}{cc}\mathrm{\&sigma;}=\frac{2(r-x)y}{\left|y\right|}+\frac{xy}{r},& r=\left|x\right|\\ \mathrm{\&sigma;}=2y-\frac{xy}{r},& r=\left|x\right|\end{array}\right.---\left(13\right)$

对于N×N的图像f(x,y)，以图像中心为坐标原点，则图像中的坐标值都满足-N/2≤(x、y)≤N/2；由式(12)和式(13)得知，r和σ的取值范围分别为1≤r≤N/2，1≤σ≤8r；则式(6)和式(7)分别改写为式(14)和式(15)；

$\mathrm{\&rho;}=\frac{2r}{N}---\left(14\right)$

θ＝4πσr (15)

经过新的变换，式(10)和式(11)可分别离散化为式(16)和式(17)；

$C_{mn}=\frac{2m+2}{N^2}\underset{r=1}{\overset{N/2}{\&Sigma;}}{R}_{mn}\left(\frac{2r}{N}\right)\underset{\mathrm{\&sigma;}=1}{\overset{8r}{\&Sigma;}}\cos \frac{\mathrm{\&pi;}n\mathrm{\&sigma;}}{4r}f(r,\mathrm{\&sigma;})---\left(16\right)$

$S_{mn}=-\frac{2m+2}{N^2}\underset{r=1}{\overset{N/2}{\&Sigma;}}{R}_{mn}\left(\frac{2r}{N}\right)\underset{\mathrm{\&sigma;}=1}{\overset{8r}{\&Sigma;}}sin\frac{\mathrm{\&pi;}n\mathrm{\&sigma;}}{4r}f(r,\mathrm{\&sigma;})---\left(17\right)$

综上所述，Zernike矩计算方法如下：

①确定图像f(x,y)的大小，从而确定式(16)和式(17)中N的值以及r和σ的取值范围；

②利用式(8)计算出R_mn(ρ)的值，并结合式(16)和式(17)计算出Zernike矩的实部C_mn和虚部S_mn；

③对C_mn和S_mn求模，得到Zernike矩的幅值|Z_mn|；

由于尺度差异问题，还使用区域r_k最小外接矩形的对角线l_k对|Z_mn|尺度归一化，计算方法如式(18)：

$\left|Z_{mn}^{\&prime;}\right|=\frac{\left|Z_{mn}\right|}{l_k}---\left(18\right)$

将Zernike矩用于显著性图中各显著性区域的形状分析，由于Zernike矩能够构造任意高阶，但考虑到计算复杂度，采用前36阶Zernike矩对区域特征进行描述；

综上所述，对区域r_k的Zernike矩描述子为：其中m∈{0,1,2,3,4,5,6,7}，m-|n|为偶数且m≥|n|；

步骤三、基于形状特征相似性度量和RANSAC的同名点求取和初始变换参数估计；

(1)同名点求取；

通过对图像各显著性区域计算质心和Zernike矩描述后，产生两组未匹配的质心点；在匹配时，正确同名点所对应的Zernike矩最具相似性，同时，由于得到的匹配点对中存在误匹配对，因此，利用RANSAC算法对匹配点对进行筛选，剔除错误的匹配点对；

设表示图像I_p中的第k个区域，参考图像为I₁，待配准图像为I₂，对于区域和区域其Zernike矩描述子分别为和则定义两区域的相似距离为式(19):

相似距离越小，则区域和区域越相似，匹配度越高；将相似距离最小的两个区域认为是匹配区域，其对应的质心为匹配点对；区域r_k的质心定义为式(20):

$m_{pq}=\underset{(x,y)\&Element;r_k}{\&Sigma;}{x}^p{y}^{q}f(x,y)$

$\left\{\begin{array}{c}\stackrel{\&OverBar;}{x}=m_{10}/m_{00}\\ \stackrel{\&OverBar;}{y}=m_{01}/m_{00}\end{array}\right.---\left(20\right)$

式中，(x,y)为区域r_k中的像素坐标，f(x,y)为像素值；即为区域r_k的质心坐标；由此得到一组匹配点对集合，但由于其中仍然包含误匹配对，因此接下来采用RANSAC算法将其中的错误匹配点对剔除；

根据误差度量函数把所有数据分为内点和外点；内点的集合又称为一致集，经过多次随机采样，当一致集达到最大时，则对应的模型即为所求结果；RANSAC算法的步骤如下：

①从候选点对集合S中随机抽取一个含有n个点对的样本集P，并估计初始化模型；

②对集合S中的剩余点对计算变换模型，与初始化模型M的距离小于误差阈值τ的点对标记为内点，形成一致集S*；

③重复上述步骤，直到达到最大采样次数N，或者集合S*的元素个数大于某一个阈值T；

④根据最大一致集S*，计算变换模型的参数；

RANSAC算法包含了四个输入参数：随机抽样数n、误差阈值τ、最大采样次数N和最优一致集阈值T；

1)随机抽样数n需要满足变换模型估计的最小配置要求，n取4组不共线的数据点对；

2)误差阈值τ对RANSAC性能有很大的影响，根据期望的内点概率来进行设置；τ为经验值，经验值取7，

3)最大采样次数N只要足够大，就能够保证匹配的正确性，N的取值计算公式如式(21)所示；

$N=\frac{log(1-p)}{log(1-{\mathrm{\&omega;}}^n)}---\left(21\right)$

其中，n为随机抽样数；ω为数据局内点的概率；p表示迭代过程中从数据集内随机选取出的点均为局内点的概率；通常数据中局内点的概率ω是未知的，但能够通过自动更新ω来自适应决定所需要采样的次数；

4)较大的一致集S*使得变换模型参数更为精确，为了得到表征数据集S的正确模型，对T不做限制，最大一致集S*取局内点最多的一致集；

(2)初始变换参数估计；

获得筛选后的正确同名点对后，选择合适的变换模型并进行变换模型参数估计；对于透视投影变换模型，待配准图像I₂中匹配点(x',y')经其变换后得到配准后图像中点(x,y)的计算公式如式(22)所示，对应的参数矩阵H₀为式(23)：

$\left\{\begin{array}{c}x=\frac{a_{11}{x}^{\&prime;}+a_{12}{y}^{\&prime;}+a_{13}}{a_{31}{x}^{\&prime;}+a_{32}{y}^{\&prime;}+1}\\ y=\frac{a_{21}{x}^{\&prime;}+a_{22}{y}^{\&prime;}+a_{23}}{a_{31}{x}^{\&prime;}+a_{32}{y}^{\&prime;}+1}\end{array}\right.---\left(22\right)$

$H_0=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& 1\end{array}\right]---\left(23\right)$

式中，H₀为待配准图像I₂到参考图像I₁坐标系的变换矩阵，至少使用4对不共线同名点对计算出H₀中所有未知量的值；

步骤四、基于边缘特征的变换参数优化求解；

(1)边缘特征；

采用边缘度的概念；对于一个结构张量Q_ρ，其定义式如下式(24)：

$Q_{\mathrm{\&rho;}}=\left(\begin{array}{cc}{K}_{\mathrm{\&rho;}}*I_x^2& K_{\mathrm{\&rho;}}*I_x{I}_y\\ K_{\mathrm{\&rho;}}*I_x{I}_y& K_{\mathrm{\&rho;}}*I_y^2\end{array}\right)---\left(24\right)$

式中，*是卷积算子，K_ρ表示一个标准差为ρ的二维高斯核，I_x和I_y分别表示x和y方向的梯度分量；通过在每个像素上对Q_ρ特征值分解，得到分别表示最大和最小的亮度变化的特征向量ω₁和ω₂，ω₁和ω₂分别对应两个特征值μ₁和μ₂，μ₁≥μ₂≥0；边缘度C的定义如下式(25)：

$C=\left\{\begin{array}{c}({\mathrm{\&mu;}}_1-{\mathrm{\&mu;}}_2)/{\mathrm{\&mu;}}_1,{\mathrm{if\&mu;}}_1>\mathrm{\&epsiv;}\\ 0,otherwise\end{array}\right.---\left(25\right)$

C的取值范围是[0,1]，ε为设定的阈值，用来防止同向区域像素取值过大；由此得到可见光-红外图像的边缘度图，它对边缘描述比常规的边缘检测方法差异更小，增加了可见光和红外图像边缘的相似性，提高了正确配准的概率；

边缘像素点(x,y)处的广义梯度向量流通过求取下式(26)的能量函数最小值得到：

$E(u,v)=\&Integral;\&Integral;\{g\left(C\right)(u_x^2+u_y^2+v_x^2+v_y^2)+\&lsqb;1-g\left(C\right)\&rsqb;\&CenterDot;{\left|\stackrel{\&RightArrow;}{v}-\&dtri;f\right|}^2\}dxdy---\left(26\right)$

式中，为像素点(x,y)处的梯度值，g(C)＝e^-C/K，C为边缘度，K为梯度场平滑程度的正则化因子，K取1；

(2)相似性度量；

当两幅图像配准时，边缘度的互信息越大，广义梯度向量流的方向误差不确定性越小，配准精度越高；因此，定义如下相似性度量函数(27)：

F(S,T)＝M(C_S,C_T)-H(D|C_S,C_T) (27)

式中，M(·)表示互信息，H(·)表示Shannon熵，C_S和C_T分别表示图像S和图像T的边缘度，D为C_S和C_T相应像素之间方向差异，其计算公式如下式(28)：

$D=arccos\left(\frac{{\stackrel{\&RightArrow;}{v}}_S\&CenterDot;{\stackrel{\&RightArrow;}{v}}_T}{\left|{\stackrel{\&RightArrow;}{v}}_S\right|\&CenterDot;\left|{\stackrel{\&RightArrow;}{v}}_T\right|}\right)---\left(28\right)$

式中，和分别是图像S和图像T的广义梯度向量流；然后，以上一步估计得到的初始配准参数H₀为初值，相似性度量函数(27)为目标函数，采用单纯形法进行优化搜索，当F(S,T)达到最大值时，配准精度最高，以此时的H为最终变换矩阵；

最后，以H为变换矩阵对待配准图像进行图像变换，然后与参考图像加权融合，得到最终结果。