CN101315698A - 基于直线特征图像配准中的特征匹配方法 - Google Patents

Abstract

本发明公开了一种基于直线特征图像配准中的特征匹配方法，目的是解决基于直线特征图像配准中的特征匹配问题。技术方案是首先从参考图像和待配准图像中提取直线编组作为配准基元，定义匹配矩阵M；接着计算这两幅图像中直线编组之间的可靠性测度、相似性测度、空间关系一致性测度，得到可靠性测度矩阵R、相似性测度矩阵B、空间关系一致性测度矩阵S；然后建立以M为自变量的能量函数，将直线编组的匹配问题转变成能量函数最优化问题；最后采用两步迭代优化算法对能量函数最优化问题进行迭代优化。本发明在特征匹配过程中有机组合了直线编组的可靠性、空间关系一致性和相似性，克服了现有方法存在的计算复杂、容易造成误匹配以及同名像点获取不准确的缺陷。

Description

基于直线特征图像配准中的特征匹配方法

技术领域

本发明属于图像处理技术领域的图像配准方法，尤其是图像配准过程中基于两幅图像中提取的直线特征进行匹配的方法。

背景技术

图像配准是各种遥感图像处理软件、医学图像处理软件中的重要模块，是遥感图像融合、环境监视、变化检测、图像拼接、天气预测以及医学图像处理应用中的关键一步，目的是将不同时期、不同视点或不同传感器获得的同一地域或物体的两幅或多幅图像进行几何校正，使它们达到空间位置上的一致性。多幅图像之间的配准问题往往可以分解成多个两幅图像之间的配准问题，因此，图像配准问题主要都是研究两幅图像之间的几何校正问题，其中一幅图像称为参考图像，另一幅图像称为待配准图像。由于图像配准的重要性，近三十年来，人们对其进行了广泛深入的研究，取得了很多重要的成果，提出了大量的方法。

已有的图像配准方法主要分为两大类：基于区域的方法和基于特征的方法。

基于区域的方法是从待配准图像中取一小窗口的图像块与参考图像中同样大小的图像块进行统计地相关比较，把相似性测度(一般采用灰度差的平方和、相关系数等，也可能是变换域的一些属性，如傅立叶变换、小波变换等)最大的两个窗口的中心作为一对同名像点，然后利用获得的所有同名像点和最小二乘算法估计变换模型参数向量，完成配准。由于采用的相似性测度一般与图像辐射特性和变换模型有关，因而基于区域的方法不适合辐射特性和几何形变较大的两幅图像之间的配准。

基于特征的方法通常包括三个步骤：首先分别从参考图像和待配准图像中提取特征，这些特征一般是空域特征，例如点、边缘、轮廓、直线段、区域等；然后利用提取特征的属性信息进行特征匹配，由匹配的特征对获得同名像点；最后根据这些同名像点和最小二乘算法估计变换模型参数向量，完成配准。跟基于区域的方法相比，基于特征的方法对图像辐射失真不太敏感，能适应于不同传感器图像之间的配准，已成为图像配准技术研究的重点和热点。

基于特征的图像配准方法的一个关键和难点是如何寻找参考图像和待配准图像中的对应特征，即特征匹配问题。现有文献中已提出多种策略来解决这一问题，主要分为利用空间关系、利用特征相似性、利用序贯组合空间关系和特征相似性这三类方法。利用空间关系的方法是首先从两幅图像提取的每个特征中抽取出一个或多个控制点(如直线的中点、区域的重心、轮廓上的突出点等)，然后寻找两幅图像中各个控制点在空间关系达到全局一致性时的最佳变换模型参数向量。这类方法的最大优点是一般都能获得较好的配准结果，不过它的最大问题是计算复杂，同时要求正确的匹配特征数必须大于任何一类局外情况的特征数目，才能获得比较可靠的结果，而且随着正确匹配特征数量的减少，算法的复杂度将会显著增加，配准的成功率却大为下降。利用特征相似性的方法是首先采用一组参数描述两幅图像中的每个特征，然后利用这些参数计算不同特征之间的相似性测度，根据相似性测度最大的原则寻找匹配特征对。特征描述参数一般要满足不变性、唯一性、稳定性以及独立性等几个条件，常用的特征描述参数有相关系数、矩不变量、形状矩阵、链码等。这类方法的特点是计算简单，但当图像中存在多个形状相同或相近的特征时，容易造成特征之间的误匹配，使得配准失败。正因为单独利用空间关系或特征相似性进行特征匹配都有局限性，人们自然想到组合这两种匹配准则，这就是利用序贯组合空间关系和特征相似性的方法，这种方法首先利用特征相似性剔除一些不太可能匹配的特征组合，以减少基于空间关系方法的搜索空间，然后利用空间关系得到最终的特征匹配结果。这类方法结合了前两类方法的优点，在配准速度和精度的兼顾性上优于前两类方法，但是这种序贯组合方式实际上仍是一种基于空间关系的方法，不能克服基于空间关系的固有缺陷：即随着两幅图像中正确匹配特征数量的减少，算法的复杂度将会显著增加、配准的成功率却大为下降。

此外，现有的绝大多数基于特征的图像配准方法在特征匹配阶段都等同地看待每一个待匹配特征，没有考虑特征提取结果的可靠性。由于图像噪声和特征提取算法等方面的原因，特征提取结果或多或少存在一些不确定性(如位置提取不准确，特征提取不完整等)。这样，如果将一些不可靠的特征作为匹配特征来估计变换模型参数向量，就会导致获得的配准结果精度不高，甚至配准失败。

基于特征的图像配准方法的关键是从参考图像和待配准图像中提取的点、线、面特征中获得同名像点。点特征本身可用作同名像点，因而它成为基于特征的图像配准方法中最常用的一种特征，但当图像分辨率较低时，点特征往往定位不准确，不适合高精度的图像配准，且点特征含有的信息量有限，匹配也比较困难。相比而言，线特征被认为更适合用于图像配准，一方面它含有比较丰富的信息量，利于匹配；另一方面，人造目标的图像中往往含有大量的线特征。面特征尽管也含有丰富的信息量，但图像中并不是总含有面特征，而且面特征通常可以用一系列的线特征表示。因此，线特征是基于特征的图像配准方法中的首选特征。现有的大部分基于直线特征的配准方法往往都是直接将直线特征作为配准基元进行匹配，然后选用匹配直线特征对(参考图像和待配准图像中对应的两条直线组成)的中点作为同名像点，当直线特征提取不完整(即直线的端点位置提取不准确)时，这样获得的同名像点往往是不正确的。

发明内容

本发明要解决的技术问题是：提出一种新的直线特征匹配方法，克服已有方法计算复杂、容易造成误匹配以及同名像点获得不准确的缺陷，解决基于直线特征图像配准中的特征匹配问题。

本发明的技术方案是：

第一步，分别从参考图像和待配准图像中提取直线编组，作为配准基元。

一般从客观景物得到的图像都是二维的，可以用一个二维数组表示。参考图像和待配准图像分别记为I₁(x₁，y₁)和I₂(x₂，y₂)，它们分别表示两幅图像的灰度值，其中(x₁，y₁)、(x₂，y₂)分别表示两幅图像中的像素点坐标。采用稳健的直线提取方法从参考图像提取N₁条直线特征，从待配准图像中提取N₂条直线特征，分别记为 $L_1=\{l_1^1,l_1^2,\·\·\·,l_1^{N_1}\}$ 和 $L_2=\{l_2^1,l_2^2,\·\·\·,l_2^{N_2}\},$ 其中N₁、N₂均为正整数(具体直线提取方法参见：文贡坚，***。一种稳健的直线提取算法。软件学报，12(11)：1660-1667，2001.)。

如背景技术中所述，直接利用直线特征进行匹配，当直线特征提取不完整时获得的同名像点往往不准确。为了克服这一缺点，本发明利用直线编组(由L₁或L₂中的任意两条直线形成的一个直线对)作为配准基元进行匹配。如果直接将参考图像和待配准图像中提取的任意两条直线特征都组成一个直线编组，获得的直线编组数量太大，增加了后续匹配的计算量。为此，只选择那些满足一定约束条件的直线编组作为配准基元，约束条件是：直线编组的提取交点到组成它的两条直线之间的Hausdroff距离之和小于一个预设门限T，T取10□30个象素。

假设参考图像中提取了N′₁个直线编组，待配准图像中提取了N′₂个直线编组，分别记为 $Z_1=\{z_1^1,z_1^2,\·\·\·,z_1^{N_1^{\′}}\}$ 和 $Z_2=\{z_2^1,z_2^2,\·\·\·,z_2^{N_2^{\′}}\},$ 其中N′₁、N′₂均为正整数。定义一个大小为N′₁×N′₂的匹配矩阵M来表示参考图像和待配准图像中直线编组之间的匹配关系，即

由于实际中参考图像中的一个直线编组只能与待配准图像中的一个直线编组相匹配，因此匹配矩阵M应满足一一匹配准则，即

$\left\{\begin{array}{c}\underset{j=1}{\overset{N_2^{\′}}{\mathrm{\Σ}}}{M}_{\mathrm{ij}}\≤1(1\≤i\≤N_1^{\′})\\ \underset{i=1}{\overset{N_1^{\′}}{\mathrm{\Σ}}}{M}_{\mathrm{ij}}\≤1(1\≤j\≤N_2^{\′})\end{array}\right.---\left(2\right)$

因此解决基于直线特征图像配准中的特征匹配问题就变为如何得到匹配矩阵M。用参数 $V=[l_H^{k1},l_H^{k2},{\mathrm{length}}_H^{k1},{\mathrm{length}}_H^{k2},S_H^{k1},S_H^{k2},E_H^{k1},E_H^{k2},M_H^{k1},M_H^{k2},O_H^k,x_H^k,y_H^k]$ 来描述参考图像和待配准图像中的任意一个直线编组z_H ^k(H＝1，k＝1，2，…，N′₁；H＝2，k＝1，2，…，N′₂)的属性，其中l_H ^k1，l_H ^k2表示组成直线编组z_H ^k的两条直线特征，length_H ^k1，length_H ^k2，S_H ^k1，S_H ^k2，E_H ^k1，E_H ^k2，M_H ^k1，M_H ^k2分别表示它们的长度、四个端点和两个中点，O_H ^k是直线编组的提取交点，x_H ^k，y_H ^k是提取交点坐标。

第二步，计算参考图像和待配准图像中直线编组之间的可靠性测度，得到可靠性测度矩阵R。

直线编组的可靠性由组成它的两条直线特征的可靠性决定。正如背景技术中提到的，实际中由于图像噪声和特征提取算法本身的原因，直线特征的提取位置与它的实际位置可能不一致，这种不一致性主要体现在两个方面：端点和方向。在方向一定的情况下，两条直线特征的端点位置对直线编组交点的提取位置没有影响，因此，只须计算直线特征的方向偏移对直线编组交点提取位置的影响。

对于L₁或L₂中任意一条直线特征l_H ^m(H＝1，m＝1，2，…，N₁；H＝2，m＝1，2，…，N₂)，假设过l_H ^m的中点M_H ^m的与l_H ^m夹角均为Δθ_H ^m的两条虚线是l_H ^m可能偏离的最大方向，Δθ_H ^m称为方向偏移。这样，l_H ^m的真实位置只可能落在两条虚线之间。一般情况下，直线特征长度越长，其提取结果越稳健，即方向偏移Δθ_H ^m越小，当直线特征达到一定长度时，其提取结果不发生方向偏移，即 $\mathrm{\Δ}{\mathrm{\θ}}_H^m=0;$ 另一方面，直线特征的方向偏移Δθ_H ^m往往会小于一定的范围。因此，可将直线特征l_H ^m的方向偏移Δθ_H ^m定义为：

$\mathrm{\Δ}{\mathrm{\θ}}_H^m=\mathrm{\αexp}[-\frac{{\left({\mathrm{length}}_H^m\right)}^2}{{{2\mathrm{\σ}}_1}^2}]---\left(3\right)$

其中α是一个预设门限，用来限定方向偏移Δθ_H ^m的最大范围，一般取α＝0.175弧度，即直线特征的方向偏移一般不超过10度；length_H ^m是l_H ^m的长度；σ₁是高斯函数的方差，它决定多长的直线特征将不发生方向偏移，一般取20□40，具体取决于图像的分辨率，图像分辨率越高，σ₁越靠近40。

对于参考图像或待配准图像中任意一个直线编组z_H ^k，z_H ^k由直线特征l_H ^k1，l_H ^k2组成，过l_H ^k1的中点M_H ^k1的与l_H ^k1夹角均为Δθ_H ^k1的两条虚线是l_H ^k1可能发生的最大方向偏移，过l_H ^k2的中点M_H ^k2的与l_H ^k2夹角均为Δθ_H ^k2的两条虚线是l_H ^k2可能发生的最大方向偏移，四条虚线交于A，B，C和D四点。Δθ_H ^k1和Δθ_H ^k2是直线特征l_H ^k1和l_H ^k2的方向偏移，它们的大小由公式(3)确定。这样，直线特征l_H ^k1和l_H ^k2的真实位置只能落在过它们各自中点的两条虚线之间，因此直线编组z_H ^k真实交点Q_H ^k必将落在由A，B，C和D四点围成的四边形r_ABCD中，且可能落在r_ABCD内的任何位置上。因此，z_H ^k的真实交点Q_H ^k是一个随机变量，记它的概率密度函数为f(Q_H ^k)。用d_H ^k表示z_H ^k的真实交点Q_H ^k和提取交点O_H ^k之间的距离，则d_H ^k也是一随机变量，记d_H ^k的期望表示为d_H ^k，则

${\stackrel{\‾}{d}}_H^k=E\left(d_H^k\right)=E\left(\right||O_H^k-Q^H|\left|\right)={\∫}_{r_{\mathrm{ABCD}}}\left|\right|O_H^k-Q_H^k\left|\right|f\left(Q_H^k\right)d{Q}_H^k---\left(4\right)$

其中E(·)表示对·求期望，‖·‖表示欧式距离，H＝1时，k＝1，2，…，N′₁，H＝2时，k＝1，2，…，N′₂。z_H ^k的真实交点Q_H ^k可以是高斯分布、均匀分布等，本发明采用均匀分布，即

$f\left(Q_H^k\right)=\frac{1}{S_{r_{\mathrm{ABCD}}}}---\left(5\right)$

其中

表示r_ABCD的面积，其大小可根据过l_H ^k1的中点M_H ^k1和过l_H ^k2的中点M_H ^k2的四条虚线方程求出，四条虚线方程通过Δθ_H ^k1和Δθ_H ^k2以及提取的直线特征l_H ^k1、l_H ^k2获得。这样，(4)式可写成

${\stackrel{\‾}{d}}_H^k=\frac{1}{S_{r_{\mathrm{ABCD}}}}{\∫}_{r_{\mathrm{ABCD}}}\left|\right|O_H^k-Q_H^k\left|\right|d{Q}_H^k---\left(6\right)$

将直线编组z_H ^k的可靠性记为c_H ^k，它与d_H ^k有关，d_H ^k越小，z_H ^k越可靠，即c_H ^k越大。c_H ^k采用下式求得：

$c_H^k=\exp [-\frac{{\left({\stackrel{\‾}{d}}_H^k\right)}^2}{{{2\mathrm{\σ}}_2}^2}]---\left(7\right)$

其中σ₂是高斯函数的方差，实际中一般取2□4。

直线编组的可靠性测度矩阵R是一个大小为N′₁×N′₂的矩阵，它的每一项R_ij表示直线编组对LP_ij的可靠性测度，LP_ij为Z₁中的任意一个直线编组z₁ ⁱ与Z₂的任意一个直线编组z₂ ^j组成的直线编组对。一般情况下，如果组成直线编组对LP_ij的两个直线编组都是可靠的，即直线编组z₁ ⁱ的可靠性c₁ ⁱ和直线编组z₂ ^j的可靠性c₂ ^j都比较大，则LP_ij也是可靠的，即R_ij也比较大。因此，可以采用D-S证据理论方法，对c₁ ⁱ和c₂ ^j进行融合计算，得到R_ij，即

$R_{\mathrm{ij}}=\frac{c_1^i\×c_2^j}{c_1^i\×c_2^j+(1-c_1^i)(1-c_2^j)}---\left(8\right)$

可见，可靠性测度矩阵R可以事先计算得到，是一个常数矩阵。

第三步，计算参考图像和待配准图像中直线编组之间的相似性测度，得到相似性测度矩阵B。

相似性测度矩阵B是一个大小为N′₁×N′₂的矩阵，它的每一项B_ij表示参考图像和待配准图像中任意两个直线编组z₁ ⁱ和z₂ ^j的相似性测度。相似性测度表示两个直线编组之间的相似程度，一般特征编组越相似，它们就越可能匹配。相似性测度B_ij定义为

$B_{\mathrm{ij}}=\left|\right|\frac{{\mathrm{length}}_1^{i1}-{\mathrm{length}}_1^{i2}}{{\mathrm{length}}_1^{i1}+{\mathrm{length}}_1^{i2}}|-|\frac{{\mathrm{length}}_2^{j1}-{\mathrm{length}}_2^{j2}}{{\mathrm{length}}_2^{j1}+{\mathrm{length}}_2^{j2}}\left|\right|---\left(9\right)$

由(9)式可知，B_ij越小表示直线编组z₁ ⁱ和z₂ ^j越相似，即它们匹配的可能性也越大。为了与直线编组之间的可靠性测度保持一致，即直线编组之间的相似性测度越大，两个直线编组越可能匹配，并且要求它的值落在0-1之间，因此需要对(9)式进行变换，即

B_ij＝F₁(B_ij)    (10)

其中F₁(·)是一个单调递减函数，采用高斯函数，即

$F_1\left(x\right)=\exp (-\frac{x^2}{{2\mathrm{\σ}}_3^2})---\left(11\right)$

其中σ₃是高斯函数的方差，在实际中一般取σ₃＝1。将B_ij代入(11)式，得

$B_{\mathrm{ij}}=\exp (-\frac{{B_{\mathrm{ij}}}^2}{2{\mathrm{\σ}}_3^2})---\left(12\right)$

可见，相似性测度矩阵B也可以事先计算得到，也是一个常数矩阵。

第四步，计算参考图像和待配准图像中直线编组之间的空间关系一致性测度，得到空间关系一致性测度矩阵S。

直线编组之间的空间关系一致性就是指参考图像中的直线编组与从待配准图像中变换而来的直线编组之间的几何临近程度。本发明采用直线编组的提取交点在变换模型约束下的欧式距离来度量它们之间的空间关系一致性。

空间关系一致性测度矩阵S也是大小N′₁×N′₂的矩阵，它的每一项S_ij表示参考图像和待配准图像中任意两个直线编组z₁ ⁱ和z₂ ^j之间的空间关系一致性测度。S_ij通过下式计算

$S_{\mathrm{ij}}=\sqrt{{(x_1^i-T_x\left(x_2^j\right))}^2+{(y_1^i-T_y\left(y_2^j\right))}^2}---\left(13\right)$

其中T_x(·)、T_y(·)分别表示参考图像和待配准图像之间的变换模型在x方向和y方向上的变换函数。参考图像和待配准图像之间的变换模型通常可以用一组变换模型参数向量P来描述，P与采用的变换模型类型有关，例如，当采用相似变换模型时，P是4维向量；当采用仿射变换模型时，P是6维向量；当采用投影变换模型时，P是8维向量。

为了与直线编组之间的可靠性测度、相似性测度保持一致，即空间关系一致性测度越大，两个直线编组越可能匹配，需要对(13)式做一变换，本发明也采用(11)式的高斯函数，即

S_ij＝F₁(S_ij)    (14)

将S_ij代入(14)式，得

$S_{\mathrm{ij}}=\exp (-\frac{{S_{\mathrm{ij}}}^2}{2{\mathrm{\σ}}_3^2})---\left(15\right)$

空间关系一致性测度矩阵S与T_x(·)、T_y(·)有关，因此与变换模型参数向量P有关。

第五步，建立能量函数

参考图像和待配准图像中的任意两个直线编组如果是正确匹配的，按照上面的计算方法，则它们之间的可靠性测度、相似性测度和空间关系一致性测度都应该比较大。为此，定义如下能量函数：

$E(M,P)=\underset{i=1}{\overset{N_1^{\′}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_2^{\′}}{\mathrm{\Σ}}}{M}_{\mathrm{ij}}{R}_{\mathrm{ij}}{B}_{\mathrm{ij}}{S}_{\mathrm{ij}}---\left(16\right)$

如果匹配矩阵M中包含错误的匹配直线编组对，则利用M获得的同名像点和最小二乘算法计算的变换模型参数向量P将是不准确的甚至是错误的，这时大部分实际匹配的直线编组对之间的空间关系一致性测度将会明显的下降，导致通过该M计算的E(·)比通过不包含错误的匹配直线编组对的M计算的E(·)要小。另一方面，如果在M中加入正确的匹配直线编组对，则利用M获得的同名像点和最小二乘算法计算的变换模型参数向量P变化很小，因此，尽管一些正确匹配直线编组对的空间关系一致性测度可能会有所减小，但由于在E(·)中增加了新的项，E(·)的值仍然会增加。因此，当匹配矩阵M包含了所有正确匹配的直线编组对而剔除了所有错误匹配的直线编组对，则E(·)的值将会达到一个全局最大值，这时的匹配矩阵M是最优的，同时由最优匹配矩阵计算的变换模型参数向量P也是最优的，即

<M_opt，P_opt>＝arg max(E(·))    (17)

其中M_opt，P_opt分别表示最优的匹配矩阵和最优的变换模型参数向量，arg max(·)表示对·求最大值。这样，直线编组的匹配问题就转变成为一个能量函数最优化问题。

第六步，采用两步迭代优化算法对(17)式描述的优化问题进行迭代优化，得到M_opt和P_opt。包括以下几个步骤：

步骤1：初始化变换模型参数向量P

采用一个序贯组合特征相似性和空间关系一致性的方法来初始化变换模型参数向量P，方法为：首先根据特征相似性测度选择求解变换模型需要的最少数目的3-5倍的直线编组对作为候选匹配直线编组对(如对于仿射变换模型，求解变换模型参数向量最少需要3对匹配直线编组对，那么就选择9-15对的直线编组对)，然后根据空间关系一致性测度利用穷举法或RANSAC算法获得求解变换模型参数向量P需要的最少数目的最优匹配直线编组对，利用它们的提取交点求出变换模型参数向量P的初值P₀。

步骤2：初始化控制参数β＝β₀

β₀、β_r和β_max分别是控制参数β的初始值、增长率和最大值，一般取β₀＝0.5，β_r取1.05□1.1，β_max＝10。

步骤3：已知变换模型参数向量P，采用软指派方法计算匹配矩阵M

已知变换模型参数向量P，空间关系一致性测度矩阵S由公式(15)事先求出，这样公式(16)就变为

$E\left(M\right)=\underset{i=1}{\overset{N_1^{\&prime;}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}{R}_{\mathrm{ij}}{B}_{\mathrm{ij}}{S}_{\mathrm{ij}}=\underset{i=1}{\overset{N_1^{\&prime;}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}{W}_{\mathrm{ij}}---\left(18\right)$

其中W_ij＝R_ijB_ijS_ij是一个常数。

求解(17)式的最大困难就是它必须满足(1)、(2)式的双向约束。为了将不等约束变为相等约束，将匹配矩阵M增加一行和一列，即当参考图像中的直线编组z₁ ⁱ不和待配准图像中的任何一个直线编组匹配时， $M_{i(N_2^{\&prime;}+1)}=1;$ 当待配准图像中的直线编组z₂ ^j不和参考图像中的任何一个直线编组匹配时， $M_{(N_1^{\&prime;}+1)j}=1.$ 这样，(2)式可变为

$\left\{\begin{array}{c}\underset{j=1}{\overset{N_2^{\&prime;}+1}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}=1(1\&le;i\&le;N_1^{\&prime;})\\ \underset{i=1}{\overset{N_1^{\&prime;}+1}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}=1(1\&le;j\&le;N_2^{\&prime;})\end{array}\right.---\left(19\right)$

显而易见，由(19)式约束的(17)式就是一个典型的分配问题。在最近的几十年里，学者们提出了大量的算法解决这一优化问题，如松弛标记算法、线性规划算法、模拟退火算法、遗传算法、神经网络算法以及渐进的软指派算法等。本发明采用Steven提出的渐进的软指派方法(具体参见S.Gold，A.Rangarajan.A Graduated Assignment Algorithm for Graph Matching，IEEETrans.on Pattern Analysis and Machine Intelligence，vol.18(4)，pp：377-387，1996.)来获得匹配矩阵M，具体步骤为：

(1)M_ij＝βW_ij；

(2)行归一化匹配矩阵M，即 $M_{\mathrm{ij}}=M_{\mathrm{ij}}/\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}};$

(3)列归一化匹配矩阵M，即 $M_{\mathrm{ij}}=M_{\mathrm{ij}}/\underset{i=1}{\overset{N_i^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}};$

(4)当收敛时，算法结束；否则转(2)。

匹配矩阵M的收敛准则是

$\underset{i=1}{\overset{N_1^{\&prime;}+1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}+1}{\mathrm{\&Sigma;}}}|M_{\mathrm{ij}}^0-M_{\mathrm{ij}}^1|<{\mathrm{\&epsiv;}}_0---\left(20\right)$

其中M_ij ⁰、M_ij ¹分别表示前一次和当前获得的匹配矩阵M中元素M_ij的值，ε₀是一个预设门限，一般取0.01-0.05。

步骤4：已知匹配矩阵M，采用多元变量变尺度算法计算变换模型参数向量P，然后让β＝β_rβ。

已知匹配矩阵M，空间关系一致性测度矩阵S的每个元素S_ij就可表示成变换模型参数向量P的函数，这样公式(16)就变为

$E\left(P\right)=\underset{i=1}{\overset{N_1^{\&prime;}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}{R}_{\mathrm{ij}}{B}_{\mathrm{ij}}{S}_{\mathrm{ij}}=\underset{i=1}{\overset{N_1^{\&prime;}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{S}_{\mathrm{ij}}\left(P\right)U_{\mathrm{ij}}---\left(21\right)$

其中U_ij＝M_ijR_ijB_ij是一个常数。显然，(17)描述的优化问题就是一个多维变量的最大化问题，目前也有很多算法可以解决它，本发明采用性能较好的多元变量变尺度算法(具体参见W.Press，S.Teukolsky，W.Vetterling and B.Flannery.Numerical Recipes in C，the Art of Scientific Computing，Second Edition，Cambridge University Press，1992.)来获得变换模型参数向量P。

步骤5：若控制参数β＞β_max，判断变换模型参数向量P是否收敛，如果收敛，算法结束，如果不收敛，返回步骤2；若控制参数β≤β_max，则转步骤3。

变换模型参数向量P的收敛准则为

‖P⁰-P¹‖＜ε₁    (22)

其中P⁰、P¹分别表示前一次和当前获得的变换模型参数向量，ε₁是一个预设门限，一般取0.01-0.05。

采用本发明可以达到以下技术效果：

(1)本发明将直线编组作为配准基元，克服了已有方法直接将直线特征作为配准基元，在特征提取不完整时造成获得的同名像点不准确的不足；

(2)本发明在直线特征匹配过程中，考虑了特征提取结果的可靠性，保证越是可靠的特征越是优先用于匹配，克服了已有方法同等看待每一个特征，当把提取不可靠的特征作为匹配特征用于配准时造成配准结果精度不高甚至配准失败的缺点；

(3)本发明有机组合了直线编组的可靠性、空间关系一致性和相似性建立能量函数，克服了基于空间关系方法、基于特征相似性方法以及序贯组合空间关系和相似性方法存在的容易造成误匹配、计算复杂的缺陷。

附图说明

图1是本发明总体流程图；

图2是本发明第二步中直线特征及其可靠性示意图；

图3是本发明第二步中直线编组及其可靠性示意图；

图4是本发明第六步中两步迭代优化算法流程图；

图5是本发明第六步中软指派算法流程图。

具体实施方式

图1是本发明总体流程图，包括以下步骤：

1、分别从参考图像和待配准图像中提取直线编组，作为配准基元。

2、计算参考图像中的直线编组和待配准图像中的直线编组之间的可靠性测度，得到可靠性测度矩阵R。

3、计算参考图像中的直线编组和待配准图像中的直线编组之间的相似性测度，得到相似性测度矩阵B。

4、计算参考图像中的直线编组和待配准图像中的直线编组之间的空间关系一致性测度，得到空间关系一致性测度矩阵S。

5、建立(16)式描述的能量函数，将直线特征匹配问题转化为(17)式描述的最优化问题。

6、采用两步迭代优化算法对(17)式描述的优化问题进行迭代优化，得到最优的匹配矩阵M_opt和变换模型参数向量P_opt。

图2是本发明第二步中直线特征及其可靠性示意图。对于参考图像或待配准图像中任意一条直线特征l_H ^m(H＝1，m＝1，2，…，N₁；H＝2，m＝1，2，…，N₂)，如图所示，假设过l_H ^m的中点M_H ^m的与l_H ^m夹角均为Δθ_H ^m的两条虚线是l_H ^m可能偏离的最大方向，l_H ^m的真实位置只可能落在两条虚线之间。

图3是本发明第二步中直线编组及其可靠性示意图。对于参考图像或待配准图像中任意一个直线编组z_H ^k，z_H ^k由直线特征l_H ^k1，l_H ^k2组成，如图所示，过l_H ^k1的中点M_H ^k1的与l_H ^k1夹角均为Δθ_H ^k1的两条虚线是l_H ^k1可能发生的最大方向偏移，过l_H ^k2的中点M_H ^k2的与l_H ^k2夹角均为Δθ_H ^k2的两条虚线是l_H ^k2可能发生的最大方向偏移，四条虚线交于A，B，C和D四点。这样，直线特征l_H ^k1和l_H ^k2的真实位置只能落在过它们各自中点的两条虚线之间，因此直线编组z_H ^k真实交点Q_H ^k必将落在由A，B，C和D四点围成的四边形r_ABCD中，且可能落在r_ABCD内的任何位置上。

图4是本发明第六步中的两步迭代优化算法的流程图，包括五个步骤：

1、初始化变换模型参数向量P＝P₀。

2、初始化控制参数β＝β₀。

3、已知变换模型参数向量P，采用软指派方法计算匹配矩阵M。

4、已知匹配矩阵M，采用多元变量变尺度算法计算变换模型参数向量P，β＝β_rβ。

5、若控制参数β＞β_max，判断变换模型参数向量P是否收敛，如果收敛，算法结束，如果不收敛，返回步骤2；若控制参数β≤β_max，则转步骤3。

图5是本发明第六步中的软指派算法流程图，包括四个步骤：

1、M_ij＝βW_ij。

2、行归一化匹配矩阵M，即 $M_{\mathrm{ij}}=M_{\mathrm{ij}}/\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}.$

3、列归一化匹配矩阵M，即 $M_{\mathrm{ij}}=M_{\mathrm{ij}}/\underset{i=1}{\overset{N_i^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}.$

4、当匹配矩阵M收敛时，算法结束；否则转2。

Claims (3)

1.一种基于直线特征图像配准中的特征匹配方法，其特征在于包括以下步骤：

第一步，分别从参考图像和待配准图像中提取直线编组作为配准基元，方法是：采用稳健的直线提取方法从参考图像I₁(x₁，y₁)提取N₁条直线特征，从待配准图像I₂(x₂，y₂)中提取N₂条直线特征，分别记为 $L_1=\{l_1^1,l_1^2,...,l_1^{N_1}\}$ 和 $L_2=\{l_2^1,l_2^2,...,l_2^{N_2}\},$ N₁、N₂均为正整数，(x₁，y₁)、(x₂，y₂)分别表示参考图像和待配准图像中的像素点坐标，由L₁或L₂中的任意两条直线形成的一个直线对叫直线编组，用参数 $V=[l_H^{k1},l_H^{k2},\mathrm{lengt}{h}_H^{k1},\mathrm{lengt}{h}_H^{k2},S_H^{k1},S_H^{k2},E_H^{k1},E_H^{k2},M_H^{k1},M_H^{k2},O_H^k,x_H^k,y_H^k]$ 来描述任意一个直线编组z_H ^k(H＝1，k＝1，2，…，N′₁；H＝2，k＝1，2，…，N′₂)的属性，其中l_H ^k1，l_H ^k2表示组成直线编组z_H ^k的两条直线特征，length_H ^k1，length_H ^k2，S_H ^k1，S_H ^k2，E_H ^k1，E_H ^k2，M_H ^k1，M_H ^k2分别表示它们的长度、四个端点和两个中点，O_H ^k是直线编组的提取交点，x_H ^k，y_H ^k是交点坐标；定义一个大小为N′₁×N′₂的匹配矩阵M来表示参考图像和待配准图像中直线编组之间的匹配关系，即

其中N′₁、N′₂均为正整数，N′₁为参考图像中提取的直线编组个数，N′₂为待配准图像中提取的直线编组个数，参考图像中提取的直线编组 $Z_1=\{z_1^1,z_1^2,...,z_1^{N_1^{\&prime;}}\},$ 待配准图像中提取的直线编组 $Z_2=\{z_2^1,z_2^2,...,z_2^{N_2^{\&prime;}}\},$ 1≤i≤N′，1≤j≤N′；匹配矩阵M满足一一匹配准则，即

$\left\{\begin{array}{c}\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}\&le;1(1\&le;i\&le;N_1^{\&prime;})\\ \underset{i=1}{\overset{N_1^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}\&le;1(1\&le;j\&le;N_2^{\&prime;})\end{array}\right.---\left(2\right)$

第二步，计算参考图像和待配准图像中直线编组之间的可靠性测度，得到可靠性测度矩阵R，R是一个大小为N′₁×N′₂的矩阵，它的每一项R_ij表示直线编组对LP_ij的可靠性测度，LP_ij为Z₁中的任意一个直线编组z₁ ⁱ与Z₂的任意一个直线编组z₂ ^j组成的直线编组对，方法是：采用D-S证据理论方法，对c₁ ⁱ和c₂ ^j进行融合计算，得到R_ij，即

$R_{\mathrm{ij}}=\frac{c_1^i\&times;c_2^j}{c_1^i\&times;c_2^j+(1-c_1^i)(1-c_2^j)}---\left(8\right)$

其中c₁ ⁱ为直线编组z₁ ⁱ的可靠性，c₂ ^j为直线编组z₂ ^j的可靠性，它们用下式计算

$c_H^k=\exp [-\frac{{\left({\stackrel{\&OverBar;}{d}}_H^k\right)}^2}{2{{\mathrm{\&sigma;}}_2}^2}]---\left(7\right)$

其中σ₂是高斯函数的方差，取2□4；H＝1时，k＝1，2，…，N′₁，H＝2时，k＝1，2，…N′₂，d_H ^k；d_H ^k为d_H ^k的期望，d_H ^k表示z_H ^k的真实交点Q_H ^k和提取交点Q_H ^k之间的距离，有

${\stackrel{\&OverBar;}{d}}_H^k=\frac{1}{S_{r_{\mathrm{ABCD}}}}{\&Integral;}_{r_{\mathrm{ABCD}}}\left|\right|O_H^k-Q_H^k\left|\right|d{Q}_H^k---\left(6\right)$

其中，z_H ^k由直线特征l_H ^k1，l_H ^k2组成，过l_H ^k1的中点M_H ^k1的与l_H ^k1夹角均为Δθ_H ^k1的两条虚线是l_H ^k1可能发生的最大方向偏移，过l_H ^k2的中点M_H ^k2的与l_H ^k2夹角均为Δθ_H ^k2的两条虚线是l_H ^k2可能发生的最大方向偏移，四条虚线交于A，B，C和D四点，Δθ_H ^k1和Δθ_H ^k2是直线特征l_H ^k1和l_H ^k2的方向偏移，它们的大小由公式(3)确定，r_ABCD为A，B，C和D四点围成的四边形，

表示r_ABCD的面积，其大小根据过l_H ^k1的中点M_H ^k1和过l_H ^k2的中点M_H ^k2的四条虚线方程求出，四条虚线方程通过Δθ_H ^k1和Δθ_H ^k2以及提取的直线特征l_H ^k1、l_H ^k2获得，

$\mathrm{\&Delta;}{\mathrm{\&theta;}}_H^m=\mathrm{\&alpha;exp}[-\frac{{\left({\mathrm{length}}_H^m\right)}^2}{2{{\mathrm{\&sigma;}}_1}^2}]---\left(3\right)$

其中α是一个预设门限，用来限定方向偏移Δθ_H ^m的最大范围，取α＝0.175弧度，即直线特征的方向偏移不超过10度；length_H ^m是l_H ^m的长度；σ₁是高斯函数的方差，它决定多长的直线特征将不发生方向偏移，取20□40，具体取决于图像的分辨率，图像分辨率越高，σ₁越靠近40；

第三步，计算参考图像和待配准图像中直线编组之间的相似性测度，得到相似性测度矩阵B，B是一个大小为N′₁×N′₂的矩阵，它的每一项B_ij表示参考图像和待配准图像中任意两个直线编组z₁ ⁱ和z₂ ^j的相似性测度，相似性测度表示两个直线编组之间的相似程度，相似性测度B_ij定义为

$B_{\mathrm{ij}}=\left|\right|\frac{{\mathrm{length}}_1^{i1}-{\mathrm{length}}_1^{i2}}{{\mathrm{length}}_1^{i1}+{\mathrm{length}}_1^{i2}}|-|\frac{{\mathrm{length}}_2^{j1}-{\mathrm{length}}_2^{j2}}{{\mathrm{length}}_2^{j1}+{\mathrm{length}}_2^{j2}}\left|\right|---\left(9\right)$

对(9)式进行变换，即

B_ij＝F₁(B_ij)(10)

其中F₁(·)是一个单调递减函数，采用高斯函数，即

$F_1\left(x\right)=\exp (-\frac{x^2}{{2\mathrm{\&sigma;}}_3^2})---\left(11\right)$

其中σ₃是高斯函数的方差，在实际中一般取σ₃＝1，将B_ij代入(11)式，得

$B_{\mathrm{ij}}=\exp (-\frac{{B_{\mathrm{ij}}}^2}{2{\mathrm{\&sigma;}}_3^2})---\left(12\right)$

第四步，计算参考图像和待配准图像中直线编组之间的空间关系一致性测度，得到空间关系一致性测度矩阵S，S也是大小N′₁×N′₂的矩阵，它的每一项S_ij表示两幅图像中任意两个直线编组z₁ ⁱ和z₂ ^j之间的空间关系一致性测度，S_ij通过下式计算

$S_{\mathrm{ij}}=\sqrt{{(x_1^i-T_x\left(x_2^j\right))}^2+{(y_1^i-T_y\left(y_2^j\right))}^2}---\left(13\right)$

其中T_x(·)、T_y(·)分别表示参考图像和待配准图像之间的变换模型在x方向和y方向上的变换函数，参考图像和待配准图像之间的变换模型用变换模型参数向量P来描述，P与采用的变换模型类型有关；采用(11)式的高斯函数对(13)式做变换，即

S_ij＝F₁(S_ij)(14)

将S_ij代入(14)式，得

$S_{\mathrm{ij}}=\exp (-\frac{{S_{\mathrm{ij}}}^2}{2{\mathrm{\&sigma;}}_3^2})---\left(15\right);$

第五步，建立能量函数：

$E(M,P)=\underset{i=1}{\overset{N_1^{\&prime;}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}{R}_{\mathrm{ij}}{B}_{\mathrm{ij}}{S}_{\mathrm{ij}}---\left(16\right)$

当匹配矩阵M包含了所有正确匹配的直线编组对而剔除了所有错误匹配的直线编组对，E(·)的值将会达到一个全局最大值，这时的匹配矩阵M是最优的，同时由最优匹配矩阵计算的变换模型参数向量P也是最优的，即

<M_opt，P_opt>＝argmax(E(·))(17)

其中M_opt，P_opt分别表示最优的匹配矩阵和最优的变换模型参数向量，arg max(·)表示对·求最大值，这样，直线编组的匹配问题转变成为一个能量函数最优化问题；

第六步，采用两步迭代优化算法对(17)式描述的优化问题进行迭代优化，得到M_opt和P_opt：

步骤1：采用一个序贯组合特征相似性和空间关系一致性的方法来初始化变换模型参数向量P：首先根据特征相似性测度选择求解变换模型需要的最少数目的3-5倍的直线编组对作为候选匹配直线编组对，然后根据空间关系一致性测度利用穷举法或RANSAC算法获得求解变换模型参数向量P需要的最少数目的最优匹配直线编组对，利用它们的提取交点求出变换模型参数向量P的初值P₀；

步骤2：初始化控制参数β＝β₀，β₀、β_r和β_max分别是控制参数β的初始值、增长率和最大值，取β₀＝0.5，β_r取1.05-1.1，β_max＝10；

步骤3：已知变换模型参数向量P，采用软指派方法计算匹配矩阵M；

步骤4：已知匹配矩阵M，空间关系一致性测度矩阵S的每个元素S_ij就可表示成变换模型参数向量P的函数，公式(16)变为

$E\left(P\right)=\underset{i=1}{\overset{N_1^{\&prime;}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}}{R}_{\mathrm{ij}}{B}_{\mathrm{ij}}{S}_{\mathrm{ij}}=\underset{i=1}{\overset{N_1^{\&prime;}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{S}_{\mathrm{ij}}\left(P\right)U_{\mathrm{ij}}---\left(21\right)$

其中U_ij＝M_ijR_ijB_ij是一个常数，采用多元变量变尺度算法计算变换模型参数向量P，然后β＝β_rβ；

步骤5：若控制参数β＞β_max，判断变换模型参数向量P是否收敛，如果收敛，算法结束，如果不收敛，返回步骤2；若控制参数β≤β_max，则转步骤3；变换模型参数向量P的收敛准则为

||P⁰-P¹||＜ε₁(22)

其中P⁰、P¹分别表示前一次和当前获得的变换模型参数向量，ε₁是一个预设门限，取0.01-0.05。

2.如权利要求1所述的基于直线特征图像配准中的特征匹配方法，其特征在于从参考图像和待配准图像中提取直线编组时只选择那些满足约束条件的直线编组作为配准基元，约束条件是：直线编组的提取交点到组成它的两条直线之间的Hausdroff距离之和小于一个预设门限T，T取10-30个象素。

3.如权利要求1所述的基于直线特征图像配准中的特征匹配方法，其特征在于所述软指派方法为：

(1)M_ij＝βW_ij；

(2)行归一化匹配矩阵M，即 $M_{\mathrm{ij}}=M_{\mathrm{ij}}/\underset{j=1}{\overset{N_2^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}};$

(3)列归一化匹配矩阵M，即 $M_{\mathrm{ij}}=M_{\mathrm{ij}}/\underset{i=1}{\overset{N_i^{\&prime;}}{\mathrm{\&Sigma;}}}{M}_{\mathrm{ij}};$

(4)当收敛时，算法结束，否则转(2)，匹配矩阵M的收敛准则是

$\underset{i=1}{\overset{N_1^{\&prime;}+1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_2^{\&prime;}+1}{\mathrm{\&Sigma;}}}|M_{\mathrm{ij}}^0-M_{\mathrm{ij}}^1|<{\mathrm{\&epsiv;}}_0---\left(20\right)$

其中M_ij ⁰、M_ij ¹分别表示前一次和当前获得的匹配矩阵M中元素M_ij的值，ε₀是一个预设门限，取0.01-0.05。

Images