CN101763633B - 基于显著性区域的可见光图像配准方法 - Google Patents

Abstract

本发明涉及图像配准领域，特别是一种基于显著性区域的可见光图像配准方法，包括：(1)加载图像。(2)提取图像的显著性区域。(3)对提取出的显著性区域，计算区域特征描述子，根据区域特征描述子的相似性进行显著性区域的匹配。(4)对步骤(3)中初步匹配上的显著性区域进行局部刚性配准。(5)采用局部刚性配准后的显著性区域中心作为控制点，进行全局二次多项式变换配准。本发明方法是一种快速、精确、鲁棒的自动图像配准方法，在图像配准方面有重大的应用价值。

Description

基于显著性区域的可见光图像配准方法

技术领域

本发明涉及图像处理，模式识别技术，特别涉及一种基于显著性区域的自动图像配准技术。

背景技术

目前主流的自动图像配准方法，主要有基于特征点的配准，基于图像灰度值的配准，基于互信息的配准方法等。这些方法都还存在一些不足之处，基于特征点的配准，对一些质量比较差的可见光图像，特征点难以准确提取，基于图像灰度值的配准要求两幅图像的灰度值必须要一致，对受光照等环境影响的图像配准精度不高，基于互信息的配准，配准需要的时间比较长，而且可能会陷入局部极值，无法得到精确配准结果。因此，对于低质量的图像很多还是采用手动配准，手动配准的成功率和精度都比较高，但是它加大了操作者的负担，配准的速度比较慢。

发明内容

针对现有技术的缺陷，本发明的目的是提供一种快速、精确、鲁棒、基于显著性区域的可见光图像配准方法。

为达成所述目的，本发明提供一种基于显著性区域的自动图像配准方法，该方法的步骤如下：

步骤1：在计算机上加载两幅待配准图像，选择一幅作为参考图像，另一幅作为浮动图像；

步骤2：将参考图像和浮动图像分成M×N个矩形区域，计算每一个区域R的局部显著性函数Ls(R)；对局部显著性函数Ls(R)进行高斯拟合，计算拟合后的局部显著性函数值Fls(R)；选择Fls(R)的局部极值区域的中心作为显著性区域中心；对每一个显著性区域中心，根据邻域内的Fls(R)分布计算区域的半径，提取参考图像和浮动图像的显著性区域R；

步骤3：对提取出的显著性区域R构建一个72维的尺度无关特征描述子Lfd(R)；定义一个距离度量函数为Dist(Lfd(R₁)，Lfd(R₂))，衡量两个特征描述子Lfd(R₁)，Lfd(R₂)间的相似性；对两幅图像的任一显著性区域对C(i，j)，计算C(i，j)的区域匹配相似度，采用由粗到精的匹配策略，进行显著性区域的匹配；

步骤4：对初步匹配上的显著性区域对Cmp(i，j)，采用基于归一化相关系数的相似性度量进行局部刚体配准；

步骤5：对局部刚体配准后的区域，采用聚类分析的方法，提取精确匹配成功的区域中心点作为控制点进行全局二次多项式变换配准，实现两幅图像的精确配准。

其中，所述局部显著性函数Ls(R)是：

Ls(R)＝Av(R)·Lge(R)，式中：Av(R)是区域R的归一化区域差异函数表示为Av(R)＝σ/μ，σ是区域R的标准差，μ是区域R的均值；Lge(R)是区域R的梯度场熵表示如下：

$\mathrm{Lge}\left(R\right)=-\underset{i=1}{\overset{36}{\mathrm{\Σ}}}{p}_i\left(R\right){\log }_2{p}_i\left(R\right),$

p_i(R)是梯度方向位于第i个扇形区域的像素点集在区域R中占的梯度幅值比例表示如下：

$p_i\left(R\right)=\frac{{\∫}_{R_i}\left|g\left(X_i\right)\right|d{X}_i}{{\∫}_R\left|g\left(X\right)\right|\mathrm{dX}},$ R_i是由所有梯度方向位于第i个扇形区域内的像素点构成的点集， $\left|g\left(X\right)\right|=\sqrt{G_x{\left(X\right)}^2+G_y{\left(X\right)}^2}$ 是像素点X梯度幅值，X_i是点集R_i中的像素点。

其中，所述拟合后的局部显著性函数值Fls(R)的计算为：

$\mathrm{Fls}\left(R_{\mathrm{ab}}\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{\mathrm{Ls}\left(R_{\mathrm{ij}}\right)}{2\mathrm{\π}{\mathrm{\σ}}^2}{\exp }^{-({(a-i)}^2+{(b-j)}^2)/\left(2{\mathrm{\σ}}^2\right)},$ 式中：用σ＝1.5的高斯核函数进行拟合，M是图像在X方向上划分的矩形区域数，N是图像在Y方向上划分的矩形区域数，(a，i)∈{1，2，...，M}，是区域R在M×N矩形区域阵列中X方向上的坐标，(b，j)∈{1，2，...，N}是区域R在M×N矩形区域阵列中Y方向上的坐标，矩形区域为R_ab。

其中，所述每一个显著性区域中心的区域半径的计算包括：以显著性区域中心所在矩形区域R_ab为中心，构建一个半径最大的正方形矩形区域集合Ω，Ω必须满足以下条件：

Fls(R_ij)≥λ·Fls(R_ab)， $\∀R_{\mathrm{ij}}\∈\mathrm{\Ω},$

式中：Fls(R_ij)是区域R_ij拟合后的局部显著性函数值，λ是矩形区域半径控制参数，经验取值为0.75；选择Ω的长度和宽度中的较小值作为显著性区域半径。

其中，所述构建一个72维的尺度无关特征描述子为：Lfd(R)＝(p₁(R)，...p₃₆(R)，da₁(R)，...da₃₆(R))，式中p_i(R)是梯度方向位于第i个扇形区域的像素点集在区域R中占的梯度幅值比例，da_i(R)∈[0，2π)是梯度方向位于第i个扇形区域的像素点集的几何中心到显著性区域中心的方向角；然后定义一个距离度量函数Dist(Lfd(R₁)，Lfd(R₂))，衡量两个特征描述子Lfd(R₁)，Lfd(R₂)间的相似性为：

$\mathrm{Dist}(\mathrm{Lfd}\left(R_1\right),\mathrm{Lfd}\left(R_2\right))=\underset{i=1}{\overset{36}{\mathrm{\Σ}}}\{\mathrm{Eud}{({\mathrm{da}}_i\left(R_1\right),{\mathrm{da}}_i\left(R_2\right))}^2\·\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right))\}$

$\·\underset{i=1}{\overset{36}{\mathrm{\Σ}}}\{\log \frac{\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right)}{\mathrm{Min}(p_i\left(R_1\right),p_i\left(R_2\right)}\·\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right))\},$

式中Eud(da_i(R₁)，da_i(R₂))是R₁和R₂中对应第i个方向角间的夹角。

其中，所述显著性区域的匹配包括：1)遍历两幅图像的每一个可能的显著性区域匹配C(i，j)，式中i表示参考图像的第i个显著性区域R_i，j表示浮动图像的第j个显著性区域R_j；满足以下条件的C(i，j)则认为是粗匹配上的显著性区域对Cmp(i，j)：

$\frac{\mathrm{Min}(\mathrm{Av}\left(R_i\right),\mathrm{Av}\left(R_j\right))}{\mathrm{Max}(\mathrm{Av}\left(R_i\right),\mathrm{Av}\left(R_j\right))}\·\frac{\mathrm{Min}(\mathrm{Lge}\left(R_i\right),\mathrm{Lge}\left(R_j\right))}{\mathrm{Max}(\mathrm{Lge}\left(R_i\right),\mathrm{Lge}\left(R_j\right))}>T,$

式中Min(·)函数是求最小值，Max(·)函数是求最大值，T是粗匹配控制参数，经验取值为0.6；

2)对每一个粗匹配区域对Cmp(i，j)，按如下方法计算R_i和R_j间的相似性S(i，j)以及旋转角度θ_ij： ${\mathrm{\θ}}_{\mathrm{ij}}=\frac{2\mathrm{k\π}}{36},$

S(i，j)＝Dist(Lfd(R_i)，Lfd(R_jk))，

k＝arg_k Min(Dist(Lfd(R_i)，Lfd(R_jk)))，k∈{0，1，...35}，

式中，R_jk是将R_j逆时针旋转10·k度得到的新区域；每一个粗匹配区域对Cmp(i，j)确定三个全局刚性变换参数为：二维平移(由R_i和R_j的区域中心决定)以及旋转参数θ_ij；

3)按S(i，j)升序排列所有的Cmp(i，j)，选取前2000(如果提取的粗匹配区域对数目不足则以实际的数目为准)个粗匹配区域对Cmp(i，j)作为输入样本集合，设置合适的类内距离阈值，对Cmp(i，j)进行全局刚性变换参数空间上的聚类，选取聚类内数目最多的类作为初步匹配上的显著性区域对F(i，j)；按S(i，j)的大小来剔除重复的区域，保证初步匹配上的显著性区域对F(i，j)中不包含重复的区域，减小后续的计算量。

其中，所述局部刚体配准是对初步匹配上的显著性区域对F(i，j)，以显著性区域R_i，R_j的中心位置和旋转角度θ_ij作为初始配准参数，进行局部刚体配准。

其中，所述实现两幅图像的精确配准包括：对局部刚体配准后的区域，设置更精细的类内距离阈值，进行全局刚性变换参数空间上的聚类分析，提取精确匹配上的区域对；并以精确匹配上的区域中心点作为控制点进行全局二次多项式变换配准，实现两幅图像的精确配准。

本发明的有益效果：本发明利用提取图像中的显著性区域，计算区域特征描述子，然后采用由粗到精的显著性区域匹配策略，将初步匹配上的显著性区域对进行局部刚体配准，最后提取局部刚体配准成功的区域中心点作为控制点进行全局二次多项式变换配准，实现两幅图像的精确配准。由于我们很好的定义了区域显著性函数和特征描述子，采用由粗到精的显著性区域匹配策略以及图像配准策略，整个算法的计算量大大减少，同时对低质量的图像也能很好的完成精确配准，算法具有非常高的鲁棒性。实验结果表明，本方法可以在4s左右完成一般图像的配准，低质量图像配准时间会长一些，在10s左右也可完成精确配准，对1548×1260的图像配准精度可以达到2个像素。因此，具有重大的应用价值。

附图说明

图1是本发明方法执行流程示意图；

图2(a)浮动图像；

图2(b)参考图像；

图2(c)是提取的浮动图像显著性区域；

图2(d)是提取的参考图像显著性区域：

图2(e)是初步匹配上的显著性区域示意图；

图3是梯度方向划分示意图；

图4(a)是显著性区域R的梯度分布示意图；

图4(b)是显著性区域特征描述子前36维特征向量示意图；

图4(c)是显著性区域特征描述子后36维特征向量示意图；

图5是配准后的浮动图像示意图；

图6是融合后的图像示意图；

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

如图1示出本发明方法执行流程示意图；在计算机101上加载两幅待配准图像，选择一幅作为参考图像，另一幅作为浮动图像；利用计算机101实现以下四个顺序处理单元对可见光图像配准：显著性区域提取单元102，用来提取可见光图像中的显著性区域；区域特征描述子计算及显著性区域匹配单元103，用来计算显著性区域的特征描述子，以及基于特征描述子间的相似性进行显著性区域匹配；局部刚体配准单元104，对初步匹配上的显著性区域进行局部刚体配准；全局二次多项式变换配准单元105，提取局部刚体配准后的区域中心点作为控制点进行全局二次多项式变换，实现图像的精确配准。

所述显著性区域提取单元102，具体实现是应用计算机101，并可以利用编程语言C++编写程序，根据我们定义的拟合后的局部显著性函数值Fls(R)，提取两幅图像中的显著性区域。

区域特征描述子计算及显著性区域匹配单元103，具体实现是应用计算机101，并可以利用编程语言C++编写程序，实现如下功能：对提取出的显著性区域，计算它的72维的尺度无关特征描述子Lfd(R)，定义距离度量函数Dist(Lfd(R₁)，Lfd(R₂))，衡量两个特征描述子Lfd(R₁)，Lfd(R₂)间的相似性，对两幅图像的任一显著性区域对，计算它们的区域匹配相似度，采用由粗到精的匹配策略，进行显著性区域的匹配。

所述局部刚体配准单元104，具体实现是应用计算机101，并可以利用编程语言C++编写程序，实现如下功能：对初步匹配上的显著性区域对，采用基于归一化相关系数的相似性度量进行局部刚体配准。

所述全局二次多项式变换配准单元105，具体实现是应用计算机101，并可以利用编程语言C++编写程序，实现如下功能：对局部刚体配准后的区域，设置更精细的类内距离阈值，进行全局刚性变换参数空间上的聚类分析，提取精确匹配成功的区域中心点作为控制点进行全局二次多项式变换配准，实现两幅图像的精确配准。

本发明的配准方法主要包括以下步骤：

步骤1：加载两幅待配准图像，应用计算机读取图像，将图像转化为二维数组，存储在计算机中，以方便后续单元进行处理，选择一幅作为参考图像，另一幅作为浮动图像，如图2(a)，图2(b)所示。

步骤2：运行显著性区域提取单元102，提取两幅图像的显著性区域。

显著性区域的提取主要由以下几步完成：

1、将参考图像和浮动图像分成M×N个矩形区域，M，N的取值跟图像大小相关，在我们的方法中，将1548×1260的图像分成100×100个矩形区域，计算每一个矩形区域R的局部显著性函数Ls(R)：

Ls(R)＝Av(R)·Lge(R)，

其中：Av(R)是区域R的归一化区域差异函数，表示为：Av(R)＝σ/μ，式中σ是区域R的标准差，μ是区域R的均值；

Lge(R)是区域R的梯度场熵，表示如下：

$\mathrm{Lge}\left(R\right)=-\underset{i=1}{\overset{36}{\mathrm{\Σ}}}{p}_i\left(R\right){\log }_2{p}_i\left(R\right)$ 式中：p_i(R)是梯度方向(像素点X的梯度方向由该点的梯度向量g(X)＝[G_x(X)，G_y(X)]确定，如图3所示，我们将整个G_xG_y平面分成36等分)位于第i个扇形区域的像素点集在区域R中占的梯度幅值比例表示如下：

$p_i\left(R\right)=\frac{{\∫}_{R_i}\left|g\left(X_i\right)\right|d{X}_i}{{\∫}_R\left|g\left(X\right)\right|\mathrm{dX}},$

式中R_i是由所有梯度方向位于第i个扇形区域内的像素点构成的点集， $\left|g\left(X\right)\right|=\sqrt{G_x{\left(X\right)}^2+G_y{\left(X\right)}^2}$ 是像素点X梯度幅值，X_i是点集R_i中的像素点。如果区域的归一化区域差异函数Av(R)的值比较小，那么区域的像素值分布比较一致，区域同质性比较强，区域显著性不明显，如果Av(R)的值比较大，区域像素值分布比较复杂，区域异质性比较强，区域显著性比较明显。另一方面，如果区域是同质的，那么局部梯度场分布应该是规则的，梯度场熵Lge(R)的值比较小，区域显著性低，如果区域是异质的，那么局部梯度场的分布就比较复杂，梯度场熵Lge(R)的值比较大，区域显著性增强。结合Av(R)和Lge(R)的区域显著性度量，Ls(R)能更好的度量区域的显著性水平。

2、在完成局部显著性函数Ls(R)的计算之后，我们采用σ＝1.5的高斯核函数进行拟合，计算拟合后的局部显著性函数值Fls(R)：

$\mathrm{Fls}\left(R_{\mathrm{ab}}\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{\mathrm{Ls}\left(R_{\mathrm{ij}}\right)}{2\mathrm{\π}{\mathrm{\σ}}^2}{\exp }^{-({(a-i)}^2+{(b-j)}^2)/\left(2{\mathrm{\σ}}^2\right)},$

M是图像在X方向上划分的矩形区域数，N是图像在Y方向上划分的矩形区域数，(a，i)∈{1，2，...，M}，是区域R在M×N矩形区域阵列中X方向上的坐标，(b，j)∈{1，2，...，N}是区域R在M×N矩形区域阵列中Y方向上的坐标，通过对Ls(R)进行高斯拟合，Fls(R)可以反映区域R的邻域的显著性水平，从而能更准确的度量R的区域显著性，所以我们选择Fls(R)的局部极值区域的中心点作为显著性区域中心。

3、在完成显著性区域中心提取后，根据邻域内的Fls(R)分布计算每一个显著性区域中心的区域半径，包括：

以显著性区域中心所在矩形区域R_ab为中心，构建一个半径最大的正方形矩形区域集合Ω，Ω必须满足以下条件：

Fls(R_ij)≥λ·Fls(R_ab)， $\∀R_{\mathrm{ij}}\∈\mathrm{\Ω},$

其中：Fls(R_ij)是区域R_ij拟合后的局部显著性函数值，λ是区域半径控制参数，经验取值为0.75。选择Ω的长度(X方向包含的像素数)和宽度(Y方向包含的像素数)中的较小值作为显著性区域半径，完成图像的显著性区域提取，如图2(c)，图2(d)所示。

步骤3：运行区域特征描述子计算及显著性区域匹配单元103，完成两幅图像的显著性区域匹配。

显著性区域匹配主要由以下几步完成：

1、对每一个显著性区域R，我们遍历区域R的每一个像素，计算其梯度幅值和梯度方向。然后统计区域R的梯度向量分布和每个梯度方向上包含的像素点集合，构建一个72维的尺度无关特征描述子Lfd(R)：

Lfd(R)＝(p₁(R)，...p₃₆(R)，da₁(R)，...da₃₆(R))，

其中p_i(R)是梯度方向位于第i个扇形区域的像素点集在区域R中占的梯度幅值比例，da_i(R)∈[0，2π)是梯度方向位于第i个扇形区域的像素点集的几何中心到显著性区域中心C的方向角。Lfd(R)的前36维特征描述了区域R内36个梯度方向上的像素点梯度幅值分布(示意图见图4(b)，限于图像大小并未画齐所有36个方向上的特征分布)，后36维特征描述了36个梯度方向的像素点几何中心相对于区域中心C的方位特征(示意图见图4(c)，限于图像大小并未画齐所有36个方向上的特征分布)。一个示意性的尺度无关特征描述子Lfd(R)的构造如图4(a)，图4(b)，图4(c)所示，其中，图4(a)中m示意性表示第三个梯度方向上的像素点集的(梯度方向以虚线表示)几何中心，C是区域中心，da₃(R)表示从C到m的方向角。

2、定义一个距离度量函数Dist(Lfd(R₁)，Lfd(R₂))，衡量两个特征描述子Lfd(R₁)，Lfd(R₂)间的相似性为：

$\mathrm{Dist}(\mathrm{Lfd}\left(R_1\right),\mathrm{Lfd}\left(R_2\right))=\underset{i=1}{\overset{36}{\mathrm{\Σ}}}\{\mathrm{Eud}{({\mathrm{da}}_i\left(R_1\right),{\mathrm{da}}_i\left(R_2\right))}^2\·\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right))\}$

$\·\underset{i=1}{\overset{36}{\mathrm{\Σ}}}\{\log \frac{\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right)}{\mathrm{Min}(p_i\left(R_1\right),p_i\left(R_2\right)}\·\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right))\},$

其中Eud(da_i(R₁)，da_i(R₂))是R₁和R₂中对应第i个方向角间的夹角。Dist(Lfd(R₁)，Lfd(R₂))的第二项类似于K-L散度的定义，但是与K-L散度相比，我们的距离度量函数定义具有对称性的优点。

3、遍历两幅图像的每一个可能的显著性区域匹配C(i，j)，其中i表示参考图像的第i个显著性区域R_i，j表示浮动图像的第j个显著性区域R_j；满足以下条件的C(i，j)则认为是粗匹配上的显著性区域对Cmp(i，j)：

$\frac{\mathrm{Min}(\mathrm{Av}\left(R_i\right),\mathrm{Av}\left(R_j\right))}{\mathrm{Max}(\mathrm{Av}\left(R_i\right),\mathrm{Av}\left(R_j\right))}\·\frac{\mathrm{Min}(\mathrm{Lge}\left(R_i\right),\mathrm{Lge}\left(R_j\right))}{\mathrm{Max}(\mathrm{Lge}\left(R_i\right),\mathrm{Lge}\left(R_j\right))}>T,$

其中Min(·)函数是求最小值，Max(·)函数是求最大值，T是粗匹配控制参数，T值设定越高，则粗匹配条件要求越严格，粗匹配上的显著性区域对Cmp(i，j)也越少，经验取值为0.6。

4、对每一个粗匹配区域对Cmp(i，j)，按如下方法计算R_i和R_j间的相似性S(i，j)以及旋转角度θ_ij：

${\mathrm{\θ}}_{\mathrm{ij}}=\frac{2\mathrm{k\π}}{36},$

S(i，j)＝Dist(Lfd(R_i)，Lfd(R_jk))，

k＝arg_k Min(Dist(Lfd(R_i)，Lfd(R_jk)))，k∈{0，1，...35}，

其中，R_jk是将R_j逆时针旋转10·k角度得到的新区域。

每一个粗匹配区域对Cmp(i，j)都可以确定三个全局刚性变换参数为：浮动图像中心相对参考图像中心的平移t_x，t_y(t_x＝O_jx-O_ix，t_y＝O_jy-O_iy)以及绕浮动图像中心的旋转角度θ_ij。其中O_ix，O_iy是显著性区域R_i中心Q_i的X和Y方向上的坐标，O_jx，O_jy是显著性区域R_j中心O_j的X和Y方向上的坐标。

5、按S(i，j)升序排列所有的Cmp(i，j)，选取前2000(如果提取的粗匹配区域对数目不足则以实际的数目为准)个粗匹配区域对Cmp(i，j)作为输入样本集合，设置合适的类内距离阈值，对Cmp(i，j)进行全局刚性变换参数空间上的聚类，选取聚类内数目最多的类作为初步匹配上的显著性区域对F(i，j)。具体的聚类方法如下：固定浮动图像的中心Cm，通过每一个粗匹配对Cmp(i，j)确定的三个刚体变换变换参数，将Cm映射到参考图像上的Cf_ij点，具体表达式如下：

$\left[\begin{array}{c}{\mathrm{Cf}}_{\mathrm{ij}}x\\ {\mathrm{Cf}}_{\mathrm{ij}}y\end{array}\right]=\left[\begin{array}{cc}\cos {\mathrm{\θ}}_{\mathrm{ij}}& -\sin {\mathrm{\θ}}_{\mathrm{ij}}\\ \sin {\mathrm{\θ}}_{\mathrm{ij}}& \cos {\mathrm{\θ}}_{\mathrm{ij}}\end{array}\right]\×\left[\begin{array}{c}\mathrm{Cmx}-O_{j}x\\ \mathrm{Cmy}-O_{j}y\end{array}\right]-\left[\begin{array}{c}{t}_x\\ t_y\end{array}\right]+\left[\begin{array}{c}{O}_{j}x\\ O_{j}y\end{array}\right]$

其中Cmx，Cmy，Cf_ijx，Cf_ijy分别是Cm和Cf_ijX和Y方向上的坐标。

我们对映射后的Cf_ij点在二维欧式空间上进行聚类，通过设置合适的类内阈值(跟图像的尺寸大小相关，在我们的实验中选取的类内阈值t为50)，选取聚类内数目最多的类作为初步匹配上的显著性区域对F(i，j)。二维欧式空间上的聚类步骤如下：

1.初始化分类数目N为0。

2.按顺序遍历映射后的Cf_ij点，如果N＝0，将当前Cf_ij作为第一类的中心，N＝N+1；反之，依次计算Cf_ij与第k(k＝1，...，N)类中心的距离d_k，选取最小值d_kmin，如果d_kmin＜t，则将Cf_ij归入第kmin类并更新第kmin类的中心以及类内元素数，如果d_kmin≥t，则将Cf_ij作为第N+1类的中心，N＝N+1。

3.选取类内数目最多的类对应的粗匹配Cmp(i，j)对作为初步匹配上的显著性区域对F(i，j)，然后，我们按S(i，j)的大小来剔除重复的区域，如果F(3，5)和F(3，7)都是F(i，j)中的元素，我们比较S(3，5)和S(3，7)的大小，剔除值比较大的匹配。保证初步匹配上的显著性区域对F(i，j)中不包含重复的区域，减小后续的计算量。初步匹配上的显著性区域对如图2(e)所示。

步骤4：运行局部刚体配准单元104。对步骤3中初步匹配上的显著性区域对F(i，j)，以区域R_i，R_j的中心确定初始平移参数t_x，t_y，以旋转角度θ_ij作为初始旋转参数，以归一化相关系数作为区域相似性度量，进行局部刚体配准。归一化相关系数的定义如下：

$\mathrm{NCC}(F,M)=-\frac{{\mathrm{\Σ}}_{i=1}^N(F_i\left(x\right)-\stackrel{\‾}{F\left(x\right)})\·(M_i\left(x\right)-\stackrel{\‾}{M\left(x\right)})}{\sqrt{{\mathrm{\Σ}}_{i=1}^N{(F_i\left(x\right)-\stackrel{\‾}{F\left(x\right)})}^2}\·\sqrt{{\mathrm{\Σ}}_{i=1}^N{(M_i\left(x\right)-\stackrel{\‾}{M\left(x\right)})}^2}},$

其中N是参考图像和浮动图像重叠区域内包含的像素数，F_i(x)和M_i(x)分别是参考图像和浮动图像第i个像素值，F(x)和M(x)分别是参考图像和浮动图像区域内像素的均值。

步骤5：运行全局二次多项式变换配准单元105。在步骤4中完成局部刚体配准后，可能存在错误的区域配准结果，在进行全局二次多项式变换前需要剔除错误的区域配准结果，在这里我们对局部刚体配准后的区域及对应的全局刚性配准参数，采用步骤3中第5步所描述的聚类分析方法，设置更精细的类内距离阈值(和图像大小相关，在我们的实验中取为20)，进行全局刚性变换参数空间上的聚类分析，选择聚类内数目最多的类作为精确匹配上的区域对。并以精确匹配上的区域中心点作为控制点进行全局二次多项式变换配准，实现两幅图像的精确配准。全局二次多项式变换配准的数学模型如下：

X_C＝A·B，

$A=\left[\begin{array}{cccccc}{a}_{00}& a_{10}& a_{01}& a_{11}& a_{20}& a_{02}\\ b_{00}& b_{10}& b_{01}& b_{11}& b_{20}& b_{02}\end{array}\right],$

$B={\left[\begin{array}{cccccc}1& x_D& y_D& x_D{y}_D& x_D^2& y_D^2\end{array}\right]}^T,$

X_D＝[x_D，y_D]^T是浮动图像中的点坐标，X_C＝[x_C，y_C]^T是参考图像中的点坐标，A是二次多项式变换矩阵，A可以通过以下方法进行求解，其中K是控制点个数：

$R=\left[\begin{array}{cccc}{x}_{C1}& x_{C2}& \·\·\·& x_{\mathrm{CK}}\\ y_{C1}& y_{C2}& \·\·\·& y_{\mathrm{CK}}\end{array}\right],$

$D=\left[\begin{array}{cccc}1& 1& \·\·\·& 1\\ x_{D1}& x_{D2}& \·\·\·& x_{\mathrm{DK}}\\ y_{D1}& y_{D2}& \·\·\·& y_{\mathrm{DK}}\\ x_{D1}{y}_{D1}& x_{D2}{y}_{D2}& \·\·\·& x_{\mathrm{DK}}{y}_{\mathrm{DK}}\\ x_{D1}^2& x_{D2}^2& \·\·\·& x_{\mathrm{DK}}^2\\ y_{D1}^2& y_{D2}^2& \·\·\·& y_{\mathrm{DK}}^2\end{array}\right],$

A＝RD^T(DD^T)^-1。

运行结果：

为了验证本发明方法，我们选取了20对图像作为实验样本，其中包含5对低质量的图像，实验结果表明，本算法可以在4s左右完成一般图像的配准，低质量的图像配准时间会长一些，在10s左右也可完成精确配准，对1548×1260的图像配准精度可以达到2个像素。具体的配准结果如图5，图6所示。实验表明，我们的方法是快速，精确，鲁棒的，具有巨大的应用价值。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (8)

1.基于显著性区域的可见光图像配准方法，其特征在于，包括以下步骤：

步骤1：在计算机上加载两幅待配准图像，选择一幅作为参考图像，另一幅作为浮动图像；

步骤2：将参考图像和浮动图像分成M×N个矩形区域，计算每一个区域R的局部显著性函数Ls(R)；对局部显著性函数Ls(R)进行高斯拟合，计算拟合后的局部显著性函数值Fls(R)；选择Fls(R)的局部极值区域的中心作为显著性区域中心；对每一个显著性区域中心，根据邻域内的Fls(R)分布计算区域的半径，提取参考图像和浮动图像的显著性区域R；

步骤3：对提取出的显著性区域R构建一个72维的尺度无关特征描述子Lfd(R)；定义一个距离度量函数为Dist(Lfd(R₁)，Lfd(R₂))，衡量两个特征描述子Lfd(R₁)，Lfd(R₂)间的相似性；对两幅图像的任一显著性区域对C(i，j)，计算C(i，j)的区域匹配相似度，采用由粗到精的匹配策略，进行显著性区域的匹配；

步骤4：对初步匹配上的显著性区域对Cmp(i，j)，采用基于归一化相关系数的相似性度量进行局部刚体配准；

步骤5：对局部刚体配准后的区域，采用聚类分析的方法，提取精确匹配成功的区域中心点作为控制点进行全局二次多项式变换配准，实现两幅图像的精确配准。

2.根据权利要求1所述的方法，其特征在于，所述局部显著性函数Ls(R)是：Ls(R)＝Av(R)·Lge(R)，其中：

Av(R)是区域R的归一化区域差异函数，表示为Av(R)＝σ/μ，式中σ是区域R的标准差，μ是区域R的均值；

Lge(R)是区域R的梯度场熵表示如下：

$\mathrm{Lge}\left(R\right)=-\underset{i=1}{\overset{36}{\mathrm{\Σ}}}{p}_i\left(R\right){\log }_2{p}_i\left(R\right),$

式中：p_i(R)是梯度方向位于第i个扇形区域的像素点集在区域R中占的梯度幅值比例表示如下：

$p_i\left(R\right)=\frac{{\∫}_{R_i}\left|g\left(X_i\right)\right|d{X}_i}{{\∫}_R\left|g\left(X\right)\right|\mathrm{dX}},$

式中R_i是由所有梯度方向位于第i个扇形区域内的像素点构成的点集，

是像素点X梯度幅值，X_i是点集R_i中的像素点。

3.根据权利要求1所述的方法，其特征在于，所述拟合后的局部显著性函数值Fls(R)的计算为：

$\mathrm{Fls}\left(R_{\mathrm{ab}}\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{\mathrm{Ls}\left(R_{\mathrm{ij}}\right)}{2\mathrm{\π}{\mathrm{\σ}}^2}{\exp }^{-({(a-i)}^2+{(b-j)}^2/\left({2\mathrm{\σ}}^2\right))},$

式中：用σ＝1.5的高斯核函数进行拟合，M是图像在X方向上划分的矩形区域数，N是图像在Y方向上划分的矩形区域数，(a，i)∈{1，2，...，M}，是区域R在M×N矩形区域阵列中X方向上的坐标，(b，j)∈{1，2，...，N}是区域R在M×N矩形区域阵列中Y方向上的坐标。

4.根据权利要求1所述的方法，其特征在于，所述每一个显著性区域中心的区域半径的计算，包括：

以显著性区域中心所在矩形区域R_ab为中心，构建一个半径最大的正方形矩形区域集合Ω，Ω必须满足以下条件：

Fls(R_ij)≥λ·Fls(R_ab)， ${\∀R}_{\mathrm{ij}}\∈\mathrm{\Ω},$

式中：Fls(R_ij)是区域R_ij拟合后的局部显著性函数值，λ是矩形区域半径控制参数，经验取值为0.75；选择Ω的长度和宽度中的较小值作为显著性区域半径。

5.根据权利要求1所述的方法，其特征在于，所述构建一个72维的尺度无关特征描述子为：

Lfd(R)＝(p₁(R)，...p₃₆(R)，da₁(R)，...da₃₆(R))，式中p_i(R)是梯度方向位于第i个扇形区域的像素点集在区域R中占的梯度幅值比例，da_i(R)∈[0，2π)是梯度方向位于第i个扇形区域的像素点集的几何中心到显著性区域中心的方向角；然后定义一个距离度量函数Dist(Lfd(R₁)，Lfd(R₂))，衡量两个特征描述子Lfd(R₁)，Lfdd(R₂)间的相似性为：

$\mathrm{Dist}(\mathrm{Lfd}\left(R_1\right),\mathrm{Lfd}\left(R_2\right))=\underset{i=1}{\overset{36}{\mathrm{\Σ}}}\{\mathrm{Eud}{({\mathrm{da}}_i\left(R_1\right),{\mathrm{da}}_i\left(R_2\right))}^2\·\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right))\}$

$\·\underset{i=1}{\overset{36}{\mathrm{\Σ}}}\{\log \frac{\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right))}{\mathrm{Min}(p_i\left(R_1\right),p_i\left(R_2\right))}\·\mathrm{Max}(p_i\left(R_1\right),p_i\left(R_2\right)),$

式中Eud(da_i(R₁)，da_i(R₂))是R₁和R₂中对应第i个方向角间的夹角。

6.根据权利要求1所述的方法，其特征在于，所述显著性区域的匹配包括：

1)遍历两幅图像的每一个可能的显著性区域匹配C(i，j)，其中i表示参考图像的第i个显著性区域R_i，j表示浮动图像的第j个显著性区域R_j；满足以下条件的C(i，j)则认为是粗匹配上的显著性区域对Cmp(i，j)：

$\frac{\mathrm{Min}\left(\mathrm{Av}(R_i,\mathrm{Av}\left(R_j\right))\right)}{\mathrm{Max}(\mathrm{Av}\left(R_i\right),\mathrm{Av}\left(R_j\right))}\·\frac{\mathrm{Min}(\mathrm{Lge}\left(R_i\right),\mathrm{Lge}\left(R_j\right))}{\mathrm{Max}(\mathrm{Lge}\left(R_i\right),\mathrm{Lge}\left(R_j\right))}>T,$

式中Av(R)是区域R的归一化区域差异函数，Lge(R)是区域R的梯度场熵，Min(·)函数是求最小值，Max(·)函数是求最大值，T是粗匹配控制参数，经验取值为0.6；

2)对每一个粗匹配区域对Cmp(i，j)，按如下方法计算R_i和R_j间的相似性S(i，j)以及旋转角度θ_ij：

${\mathrm{\θ}}_{\mathrm{ij}}=\frac{2\mathrm{k\π}}{36},$

S(i，j)＝Dist(Lfd(R_i)，Lfd(R_jk))，

k＝arg_k Min(Dist(Lfd(R_i)，Lfd(R_jk)))，k∈{0，1，...35}，

式中，R_jk是将R_j逆时针旋转10·k度得到的新区域；

每一个粗匹配区域对Cmp(i，j)确定三个全局刚性变换参数为：二维平移以及旋转参数θ_ij，二维平移由R_i和R_j的区域中心决定；

3)按S(i，j)升序排列所有的Cmp(i，j)，选取前2000个粗匹配区域对Cmp(i，j)作为输入样本集合，如果提取的粗匹配区域对数目不足则以实际的数目为准，设置合适的类内距离阈值，对Cmp(i，j)进行全局刚性变换参数空间上的聚类，选取聚类内数目最多的类作为初步匹配上的显著性区域对F(i，j)；按S(i，j)的大小来剔除重复的区域，保证初步匹配上的显著性区域对F(i，j)中不包含重复的区域，减小后续的计算量。

7.根据权利要求1所述的方法，其特征在于，所述局部刚体配准是对初步匹配上的显著性区域对F(i，j)，以显著性区域R_i，R_j的中心位置和旋转角度θ_ij作为初始配准参数，进行局部刚体配准。

8.根据权利要求1所述的方法，其特征在于，所述实现两幅图像的精确配准包括：对局部刚体配准后的区域，设置更精细的类内距离阈值，进行全局刚性变换参数空间上的聚类分析，提取精确匹配上的区域对；并以精确匹配上的区域中心点作为控制点进行全局二次多项式变换配准，实现两幅图像的精确配准。

Images