专利名称：：基于直线特征图像配准中的特征匹配方法
技术领域：
：本发明属于图像处理
技术领域：
的图像配准方法，尤其是图像配准过程中基于两幅图像中提取的直线特征进行匹配的方法。
背景技术：
：图像配准是各种遥感图像处理软件、医学图像处理软件中的重要模块，是遥感图像融合、环境监视、变化检测、图像拼接、天气预测以及医学图像处理应用中的关键一步，目的是将不同时期、不同视点或不同传感器获得的同一地域或物体的两幅或多幅图像进行几何校正，使它们达到空间位置上的一致性。多幅图像之间的配准问题往往可以分解成多个两幅图像之间的配准问题，因此，图像配准问题主要都是研究两幅图像之间的几何校正问题，其中一幅图像称为参考图像，另一幅图像称为待配准图像。由于图像配准的重要性，近三十年来，人们对其进行了广泛深入的研究，取得了很多重要的成果，提出了大量的方法。已有的图像配准方法主要分为两大类基于区域的方法和基于特征的方法。基于区域的方法是从待配准图像中取一小窗口的图像块与参考图像中同样大小的图像块进行统计地相关比较，把相似性测度(一般采用灰度差的平方和、相关系数等，也可能是变换域的一些属性，如傅立叶变换、小波变换等)最大的两个窗口的中心作为一对同名像点，然后利用获得的所有同名像点和最小二乘算法估计变换模型参数向量，完成配准。由于采用的相似性测度一般与图像辐射特性和变换模型有关，因而基于区域的方法不适合辐射特性和几何形变较大的两幅图像之间的配准。基于特征的方法通常包括三个步骤首先分别从参考图像和待配准图像中提取特征，这些特征一般是空域特征，例如点、边缘、轮廓、直线段、区域等；然后利用提取特征的属性信息进行特征匹配，由匹配的特征对获得同名像点；最后根据这些同名像点和最小二乘算法估计变换模型参数向量，完成配准。跟基于区域的方法相比，基于特征的方法对图像辐射失真不太敏感，能适应于不同传感器图像之间的配准，已成为图像配准技术研究的重点和热点。基于特征的图像配准方法的一个关键和难点是如何寻找参考图像和待配准图像中的对应特征，即特征匹配问题。现有文献中已提出多种策略来解决这一问题，主要分为利用空间关系、利用特征相似性、利用序贯组合空间关系和特征相似性这三类方法。利用空间关系的方法是首先从两幅图像提取的每个特征中抽取出一个或多个控制点(如直线的中点、区域的重心、轮廓上的突出点等)，然后寻找两幅图像中各个控制点在空间关系达到全局一致性时的最佳变换模型参数向量。这类方法的最大优点是一般都能获得较好的配准结果，不过它的最大问题是计算复杂，同时要求正确的匹配特征数必须大于任何一类局外情况的特征数目，才能获得比较可靠的结果，而且随着正确匹配特征数量的减少，算法的复杂度将会显著增加，配准的成功率却大为下降。利用特征相似性的方法是首先采用一组参数描述两幅图像中的每个特征，然后利用这些参数计算不同特征之间的相似性测度，根据相似性测度最大的原则寻找匹配特征对。特征描述参数一般要满足不变性、唯一性、稳定性以及独立性等几个条件，常用的特征描述参数有相关系数、矩不变量、形状矩阵、链码等。这类方法的特点是计算简单，但当图像中存在多个形状相同或相近的特征时，容易造成特征之间的误匹配，使得配准失败。正因为单独利用空间关系或特征相似性进行特征匹配都有局限性，人们自然想到组合这两种匹配准则，这就是利用序贯组合空间关系和特征相似性的方法，这种方法首先利用特征相似性剔除一些不太可能匹配的特征组合，以减少基于空间关系方法的搜索空间，然后利用空间关系得到最终的特征匹配结果。这类方法结合了前两类方法的优点，在配准速度和精度的兼顾性上优于前两类方法，但是这种序贯组合方式实际上仍是一种基于空间关系的方法，不能克服基于空间关系的固有缺陷即随着两幅图像中正确匹配特征数量的减少，算法的复杂度将会显著增加、配准的成功率却大为下降。此外，现有的绝大多数基于特征的图像配准方法在特征匹配阶段都等同地看待每一个待匹配特征，没有考虑特征提取结果的可靠性。由于图像噪声和特征提取算法等方面的原因，特征提取结果或多或少存在一些不确定性(如位置提取不准确，特征提取不完整等)。这样，如果将一些不可靠的特征作为匹配特征来估计变换模型参数向量，就会导致获得的配准结果精度不高，甚至配准失败。基于特征的图像配准方法的关键是从参考图像和待配准图像中提取的点、线、面特征中获得同名像点。点特征本身可用作同名像点，因而它成为基于特征的图像配准方法中最常用的一种特征，但当图像分辨率较低时，点特征往往定位不准确，不适合高精度的图像配准，且点特征含有的信息量有限，匹配也比较困难。相比而言，线特征被认为更适合用于图像配准，一方面它含有比较丰富的信息量，利于匹配；另一方面，人造目标的图像中往往含有大量的线特征。面特征尽管也含有丰富的信息量，但图像中并不是总含有面特征，而且面特征通常可以用一系列的线特征表示。因此，线特征是基于特征的图像配准方法中的首选特征。现有的大部分基于直线特征的配准方法往往都是直接将直线特征作为配准基元进行匹配，然后选用匹配直线特征对(参考图像和待配准图像中对应的两条直线组成)的中点作为同名像点，当直线特征提取不完整(即直线的端点位置提取不准确)时，这样获得的同名像点往往是不正确的。
发明内容本发明要解决的技术问题是提出一种新的直线特征匹配方法，克服已有方法计算复杂、容易造成误匹配以及同名像点获得不准确的缺陷，解决基于直线特征图像配准中的特征匹配问题。本发明的技术方案是-第一步，分别从参考图像和待配准图像中提取直线编组，作为配准基元。一般从客观景物得到的图像都是二维的，可以用一个二维数组表示。参考图像和待配准图像分别记为A(^:^)和/2(^h)，它们分别表'示两幅图像的灰度值，其中0c,，乂)、(A,h)分别表示两幅图像中的像素点坐标。采用稳健的直线提取方法从参考图像提取m条直线特征，从待配准图像中提取7V2条直线特征，分别记为a-{Auv'^r'^f^^/"2,其中m、^均为正整数(具体直线提取方法参见文贡坚，王润生.一种稳健的直线提取算法.软件学报，12(11):1660-1667，2001.)。如
背景技术：
中所述，直接利用直线特征进行匹配，当直线特征提取不完整时获得的同名像点往往不准确。为了克服这一缺点，本发明利用直线编组(由a或a中的任意两条直线形成的一个直线对)作为配准基元进行匹配。如果直接将参考图像和待配准图像中提取的任意两条直线特征都组成一个直线编组，获得的直线编组数量太大，增加了后续匹配的计算量。为此，只选择那些满足一定约束条件的直线编组作为配准基元，约束条件是直线编组的提取交点到组成它的两条直线之间的Hausdroff距离之和小于一个预设门限r，r取10[]30个象素。假设参考图像中提取了M'个直线编组，待配准图像中提取了A^个直线编组，分别记为Z,-^^,…,d和Z2-W，Z22,…,《2)，其中7V;、A^均为正整数。定义一个大小为A^^的匹配矩阵M来表示参考图像和待配准图像中直线编组之间的匹配关系，即由于实际中参考图像中的一个直线编组只能与待配准图像中的一个直线编组相匹配，因此匹配矩阵M应满足一一匹配准则，即因此解决基于直线特征图像配准中的特征匹配问题就变为如何得到匹配矩阵M。用参数V4/仏/^,/e"g《,/e"g《X'X2,《，《,《',M^,《，4,乂]来描述参考图像和待配准图像中的任意一个直线编组4(/f"，"l，2,…,A^;/f-2,"l,2，…,AO的属性，其中《,《2表示组成直线编组4的两条直线特征，/e"g《,/e"g《,C《2,《，《,《',M〖2分别表示它们的长度、四个端点和两个中点，《是直线编组的提取交点，4,^是提取交点坐第二步，计算参考图像和待配准图像中直线编组之间的可靠性测度，得到可靠性测度矩阵R。直线编组的可靠性由组成它的两条直线特征的可靠性决定。正如
背景技术：
中提到的，实际中由于图像噪声和特征提取算法本身的原因，直线特征的提取位置与它的实际位置可能不一致，这种不一致性主要体现在两个方面端点和方向。在方向一定的情况下，两条直线特征的端点位置对直线编组交点的提取位置没有影响，因此，只须计算直线特征的方向偏移对直线编组交点提取位置的影响。对于A或丄2中任意一条直线特征=1，m=1，2,…，M;7/=2，m=1,2,…,Ag，假设过^的中点M-的与/卩夹角均为A^;的两条虚线是/;;可能偏离的最大方向，A^称为方向偏移。这样，G的真实位置只可能落在两条虚线之间。一般情况下，直线特征长度越长，其提取结果越稳健，即方向偏移A^越小，当直线特征达到一定长度时，其提取结果不发生方向偏移，即A^-0;另一方面，直线特征的方向偏移A^往往会小于一定的范围。因此，可将直线特征/;;的方向偏移A^定义为A[證p[-^f^](3)其中a是一个预设门限，用来限定方向偏移A^的最大范围，一般取"=0.175弧度，即直线特征的方向偏移一般不超过10度；/ew^^是G的长度；C7,是高斯函数的方差，它决定多长的直线特征将不发生方向偏移，一般取20D40，具体取决于图像的分辨率，图像分辨率越高，^越靠近40。对于参考图像或待配准图像中任意一个直线编组4，4由直线特征《,/^组成，过《的中点M^的与《夹角均为A^的两条虚线是/可能发生的最大方向偏移，过/〖2的中点M〖2的与^夹角均为A《2的两条虚线是/〖2可能发生的最大方向偏移，四条虚线交于4AC和D四点。A^和A《2是直线特征《和《的方向偏移，它们的大小由公式(3)确定。这样，直线特征/^和/^的真实位置只能落在过它们各自中点的两条虚线之间，因此直线编组《真实交点《必将落在由AAC和D四点围成的四边形。，中，且可能落在。，内的任何位置上。因此，4的真实交点《是一个随机变量，记它的概率密度函数为/(&)。用《表示《的真实交点《和提取交点《之间的距离，则4也是一随机变量，记《的期望表示为《，则《=g)=—&||)=l。K-^lk必)赠(4)其中￡(*)表示对*求期望，l卜l表示欧式距离，//=1时，"1,2,…,7V;,//=2时，*=1,2,z^的真实交点&可以是高斯分布、均匀分布等，本发明采用均匀分布，即/必)=^~(5)其中&,。表示^cd的面积，其大小可根据过^的中点M〖1和过/〖2的中点M-2的四条虚线方程求出，四条虚线方程通过A《1和A《2以及提取的直线特征e、《2获得。这样，(4)式可写成将直线编组4的可靠性记为4，它与《有关，《越小，^越可靠，即4越大。《采用下式求得4=exp[—^(7)其中^是高斯函数的方差，实际中一般取2口4。直线编组的可靠性测度矩阵R是一个大小为A^W的矩阵，它的每一项R,表示直线编组对Z4的可靠性测度，i4为Z,中的任意一个直线编组z;与Z2的任意一个直线编组^组成的直线编组对。一般情况下，如果组成直线编组对丄A的两个直线编组都是可靠的，即直线编组z;的可靠性《和直线编组^的可靠性《都比较大，则丄A也是可靠的，即Ry也比较大。因此，可以采用D一s证据理论方法，对c;和《进行融合计算，得到r,，即qx《可见，可靠性测度矩阵R可以事先计算得到，是一个常数矩阵。(8)第三步，计算参考图像和待配准图像中直线编组之间的相似性测度，得到相似性测度矩阵B。相似性测度矩阵B是一个大小为A7xiV;的矩阵，它的每一项B,表示参考图像和待配准图像中任意两个直线编组z,'和^的相似性测度。相似性测度表示两个直线编组之间的相似程度，一般特征编组越相似，它们就越可能匹配。相似性测度8,,定义为<formula>formulaseeoriginaldocumentpage15</formula>(9)由(9)式可知，B,越小表示直线编组《和^越相似，即它们匹配的可能性也越大。为了与直线编组之间的可靠性测度保持一致，即直线编组之间的相似性测度越大，两个直线编组越可能匹配，并且要求它的值落在o-i之间，因此需要对(9)式进行变换，即By=W)(10)其中巧(，)是一个单调递减函数，采用高斯函数，即<formula>formulaseeoriginaldocumentpage15</formula>其中A是高斯函数的方差，在实际中一般取A^。将B,代入(11)式，得<formula>formulaseeoriginaldocumentpage16</formula>(12)可见，相似性测度矩阵B也可以事先计算得到，也是一个常数矩阵。第四步，计算参考图像和待配准图像中直线编组之间的空间关系一致性测度，得到空间关系一致性测度矩阵s。直线编组之间的空间关系一致性就是指参考图像中的直线编组与从待配准图像中变换而来的直线编组之间的几何临近程度。本发明采用直线编组的提取交点在变换模型约束下的欧式距离来度量它们之间的空间关系一致性。空间关系一致性测度矩阵s也是大小7v;x^的矩阵，它的每一项s&表示参考图像和待配准图像中任意两个直线编组z;和s之间的空间关系一致性测度。s,通过下式计算<formula>formulaseeoriginaldocumentpage16</formula>(13)其中t;w、？;("分别表示参考图像和待配准图像之间的变换模型在;c方向和y方向上的变换函数。参考图像和待配准图像之间的变换模型通常可以用一组变换模型参数向量p来描述，p与采用的变换模型类型有关，例如，当采用相似变换模型时，P是4维向量；当采用仿射变换模型时，P是6维向量；当采用投影变换模型时，P是8维向量。为了与直线编组之间的可靠性测度、相似性测度保持一致，即空间关系一致性测度越大，两个直线编组越可能匹配，需要对(13)式做一变换，本发明也采用(ll)式的高斯函数，即S「A(S,)(14)将S,代入(14)式，得<formula>formulaseeoriginaldocumentpage17</formula>空间关系一致性测度矩阵s与7;(，)、z;("有关，因此与变换模型参数向量p有关。第五步，建立能量函数参考图像和待配准图像中的任意两个直线编组如果是正确匹配的，按照上面的计算方法，则它们之间的可靠性测度、相似性测度和空间关系一致性测度都应该比较大。为此，定义如下能量函数如果匹配矩阵M中包含错误的匹配直线编组对，则利用M获得的同名像点和最小二乘算法计算的变换模型参数向量P将是不准确的甚至是错误的，这时大部分实际匹配的直线编组对之间的空间关系一致性测度将会明显的下降，导致通过该M计算的^—比通过不包含错误的匹配直线编组对的M计算的^:(，)要小。另一方面，如果在M中加入正确的匹配直线编组对，则利用M获得的同名像点和最小二乘算法计算的变换模型参数向量P变化很小，因此，尽管一些正确匹配直线编组对的空间关系一致性测度可能会有所减小，但由于在￡(*)中增加了新的项，￡—)的值仍然会增加。因此，当匹配矩阵M包含了所有正确匹配的直线编组对而剔除了所有错误匹配的直线编组对，则￡(*)的值将会达到一个全局最大值，这时的匹配矩阵M是最优的，同时由最优匹配矩阵计算的变换模型参数向量P也是最优的，即〈Mop,,P。p,〉=argmax(￡())(17)其中M。一P^分别表示最优的匹配矩阵和最优的变换模型参数向量，argmax(，)表示对，求最大值。这样，直线编组的匹配问题就转变成为一个能量函数最优化问题。第六步，采用两步迭代优化算法对(17)式描述的优化问题进行迭代优化，得到M—和P。p,。包括以下几个步骤步骤1:初始化变换模型参数向量P采用一个序贯组合特征相似性和空间关系一致性的方法来初始化变换模型参数向量P，方法为首先根据特征相似性测度选择求解变换模型需要的最少数目的3—5倍的直线编组对作为候选匹配直线编组对(如对于仿射变换模型，求解变换模型参数向量最少需要3对匹配直线编组对，那么就选择9一15对的直线编组对)，然后根据空间关系一致性测度利用穷举法或RANSAC算法获得求解变换模型参数向量P需要的最少数目的最优匹配直线编组对，利用它们的提取交点求出变换模型参数向量P的初值P。。步骤2:初始化控制参数々=々。A、A和A^分别是控制参数"的初始值、增长率和最大值，一般取A=0.5，A取1.05。1.1，〃max=10。步骤3:已知变换模型参数向量P，采用软指派方法计算匹配矩阵M已知变换模型参数向量P，空间关系一致性测度矩阵S由公式(15)事先求出，这样公式(16)就变为￡(M)=iiMaba=iiMyw,(is)其中W^-R^^,是一个常数。求解(17)式的最大困难就是它必须满足(1)、(2)式的双向约束。为了将不等约束变为相等约束，将匹配矩阵M增加一行和一列，即当参考图像中的直线编组z;不和待配准图像中的任何一个直线编组匹配时，M=1;当待配准图像中的直线编组^不和参考图像中的任何一个直线编组匹配时，M(和)广l。这样，(2)式可变为<formula>formulaseeoriginaldocumentpage19</formula>(19)显而易见，由(19)式约束的(17)式就是一个典型的分配问题。在最近的几十年里，学者们提出了大量的算法解决这一优化问题，如松弛标记算法、线性规划算法、模拟退火算法、遗传算法、神经网络算法以及渐进的软指派算法等。本发明采用Steven提出的渐进的软指派方法(具体参见S.Gold，A.Rangarajan.AGraduatedAssignmentAlgorithmforGraphMatching,IEEETrans,onPatternAnalysisandMachineIntelligence,vol.18(4)，pp:377—387,1996.)来获得匹配矩阵M，具体步骤为(1)My=;^;(2)行归一化匹配矩阵M，即M,-M々/fx;(3)列归一化匹配矩阵M，即M,-Myi;M,;(4)当收敛时，算法结束；否则转(2)。匹配矩阵M的收敛准则是<formula>formulaseeoriginaldocumentpage19</formula>(20)/=1声l其中M/、1VV分别表示前一次和当前获得的匹配矩阵M中元素M,.的值，c。是一个预设门限，一般取0.01-0.05。步骤4:已知匹配矩阵M，采用多元变量变尺度算法计算变换模型参数向量P，然后让/=/^。已知匹配矩阵M，空间关系一致性测度矩阵S的每个元素S,就可表示成变换模型参数向量P的函数，这样公式(16)就变为<formula>formulaseeoriginaldocumentpage20</formula>(21)其中仏-M,R^々是一个常数。显然，(17)描述的优化问题就是一个多维变量的最大化问题，目前也有很多算法可以解决它，本发明采用性能较好的多元变量变尺度算法(具体参见W.Press,S.Teukolsky,W.VetterlingandB.Flannery.NumericalRecipesinC，theArtofScientificComputing,SecondEdition,CambridgeUniversityPress,1992.)来获得变换模型参数向量P。步骤5:若控制参数">/_，判断变换模型参数向量P是否收敛，如果收敛，算法结束，如果不收敛，返回步骤2;若控制参数^/_，贝U转步骤3。变换模型参数向量p的收敛准则为|p。-p'l",(22)其中P。、P'分别表示前一次和当前获得的变换模型参数向量，^是一个预设门限，一般取0.01-0.05。采用本发明可以达到以下技术效果(1)本发明将直线编组作为配准基元，克服了已有方法直接将直线特征作为配准基元，在特征提取不完整时造成获得的同名像点不准确的不足；(2)本发明在直线特征匹配过程中，考虑了特征提取结果的可靠性，保证越是可靠的特征越是优先用于匹配，克服了已有方法同等看待每一个特征，当把提取不可靠的特征作为匹配特征用于配准时造成配准结果精度不高甚至配准失败的缺点；(3)本发明有机组合了直线编组的可靠性、空间关系一致性和相似性建立能量函数，克服了基于空间关系方法、基于特征相似性方法以及序贯组合空间关系和相似性方法存在的容易造成误匹配、计算复杂的缺陷。图l是本发明总体流程图；图2是本发明第二步中直线特征及其可靠性示意图；图3是本发明第二步中直线编组及其可靠性示意图；图4是本发明第六步中两步迭代优化算法流程图；图5是本发明第六步中软指派算法流程图。具体实施方式图1是本发明总体流程图，包括以下步骤1、分别从参考图像和待配准图像中提取直线编组，作为配准基元。2、计算参考图像中的直线编组和待配准图像中的直线编组之间的可靠性测度，得到可靠性测度矩阵R。3、计算参考图像中的直线编组和待配准图像中的直线编组之间的相似性测度，得到相似性测度矩阵B。4、计算参考图像中的直线编组和待配准图像中的直线编组之间的空间关系一致性测度，得到空间关系一致性测度矩阵S。5、建立(16)式描述的能量函数，将直线特征匹配问题转化为(17)式描述的最优化问题。6、采用两步迭代优化算法对(17)式描述的优化问题进行迭代优化，得到最优的匹配矩阵M。p,和变换模型参数向量P。p,。图2是本发明第二步中直线特征及其可靠性示意图。对于参考图像或待配准图像中任意一条直线特征/-(//=1,附=1,2广.,7^;//=2,附=1,2广.^2)，如图所示，假设过/;;的中点M^;的与/:夹角均为A^的两条虚线是/;;可能偏离的最大方向，的真实位置只可能落在两条虚线之间。图3是本发明第二步中直线编组及其可靠性示意图。对于参考图像或待配准图像中任意一个直线编组《，^由直线特征《,/^组成，如图所示，过《的中点M^的与《夹角均为A《的两条虚线是^可能发生的最大方向偏移，过^的中点M〖2的与^夹角均为A《2的两条虚线是^可能发生的最大方向偏移，四条虚线交于4AC和Z)四点。这样，直线特征C和《的真实位置只能落在过它们各自中点的两条虚线之间，因此直线编组4真实交点《必将落在由X,5,C和D四点围成的四边形。，中，且可能落在内的任何位置上。图4是本发明第六步中的两步迭代优化算法的流程图，包括五个步骤1、初始化变换模型参数向量P-P。。2、初始化控制参数"-A)。3、已知变换模型参数向量P，采用软指派方法计算匹配矩阵M。4、已知匹配矩阵M，采用多元变量变尺度算法计算变换模型参数向量P，5、若控制参数">;9_,判断变换模型参数向量P是否收敛，如果收敛，算法结束，如果不收敛，返回步骤2;若控制参数々^^_，则转步骤3。图5是本发明第六步中的软指派算法流程图，包括四个步骤1、M广;^.。2、行归一化匹配矩阵M，即M,-M々/艺M3、列归一化匹配矩阵M，即M々-M,/IX。4、当匹配矩阵M收敛时，算法结束；否则转2。权利要求1.一种基于直线特征图像配准中的特征匹配方法，其特征在于包括以下步骤第一步，分别从参考图像和待配准图像中提取直线编组作为配准基元，方法是采用稳健的直线提取方法从参考图像I1(x1，y1)提取N1条直线特征，从待配准图像I2(x2，y2)中提取N2条直线特征，分别记为<mathsid="math0001"num="0001"><math><![CDATA[<mrow><msub><mi>L</mi><mn>1</mn></msub><mo>=</mo><mo>{</mo><msubsup><mi>l</mi><mn>1</mn><mn>1</mn></msubsup><mo>,</mo><msubsup><mi>l</mi><mn>1</mn><mn>2</mn></msubsup><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><msubsup><mi>l</mi><mn>1</mn><msub><mi>N</mi><mn>1</mn></msub></msubsup><mo>}</mo></mrow>]]></math>id="icf0001"file="A2008100315750002C1.tif"wi="29"he="4"top="73"left="149"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/></maths>和<mathsid="math0002"num="0002"><math><![CDATA[<mrow><msub><mi>L</mi><mn>2</mn></msub><mo>=</mo><mo>{</mo><msubsup><mi>l</mi><mn>2</mn><mn>1</mn></msubsup><mo>,</mo><msubsup><mi>l</mi><mn>2</mn><mn>2</mn></msubsup><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><msubsup><mi>l</mi><mn>2</mn><msub><mi>N</mi><mn>2</mn></msub></msubsup><mo>}</mo><mo>,</mo></mrow>]]></math>id="icf0002"file="A2008100315750002C2.tif"wi="33"he="5"top="84"left="22"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/></maths>N1、N2均为正整数，(x1，y1)、(x2，y2)分别表示参考图像和待配准图像中的像素点坐标，由L1或L2中的任意两条直线形成的一个直线对叫直线编组，用参数<mathsid="math0003"num="0003"><math><![CDATA[<mrow><mi>V</mi><mo>=</mo><mo>[</mo><msubsup><mi>l</mi><mi>H</mi><mrow><mi>k</mi><mn>1</mn></mrow></msubsup><mo>,</mo><msubsup><mi>l</mi><mi>H</mi><mrow><mi>k</mi><mn>2</mn></mrow></msubsup><mo>,</mo><mi>lengt</mi><msubsup><mi>h</mi><mi>H</mi><mrow><mi>k</mi><mn>1</mn></mrow></msubsup><mo>,</mo><mi>lengt</mi><msubsup><mi>h</mi><mi>H</mi><mrow><mi>k</mi><mn>2</mn></mrow></msubsup><mo>,</mo><msubsup><mi>S</mi><mi>H</mi><mrow><mi>k</mi><mn>1</mn></mrow></msubsup><mo>,</mo><msubsup><mi>S</mi><mi>H</mi><mrow><mi>k</mi><mn>2</mn></mrow></msubsup><mo>,</mo><msubsup><mi>E</mi><mi>H</mi><mrow><mi>k</mi><mn>1</mn></mrow></msubsup><mo>,</mo><msubsup><mi>E</mi><mi>H</mi><mrow><mi>k</mi><mn>2</mn></mrow></msubsup><mo>,</mo><msubsup><mi>M</mi><mi>H</mi><mrow><mi>k</mi><mn>1</mn></mrow></msubsup><mo>,</mo><msubsup><mi>M</mi><mi>H</mi><mrow><mi>k</mi><mn>2</mn></mrow></msubsup><mo>,</mo><msubsup><mi>O</mi><mi>H</mi><mi>k</mi></msubsup><mo>,</mo><msubsup><mi>x</mi><mi>H</mi><mi>k</mi></msubsup><mo>,</mo><msubsup><mi>y</mi><mi>H</mi><mi>k</mi></msubsup><mo>]</mo></mrow>]]></math>id="icf0003"file="A2008100315750002C3.tif"wi="118"he="5"top="106"left="66"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/></maths>来描述任意一个直线编组zHk(H＝1，k＝1，2，…，N′1；H＝2，k＝1，2，…，N′2)的属性，其中lHk1，lHk2表示组成直线编组zHk的两条直线特征，lengthHk1，lengthHk2，SHk1，SHk2，EHk1，EHk2，MHk1，MHk2分别表示它们的长度、四个端点和两个中点，OHk是直线编组的提取交点，xHk，yHk是交点坐标；定义一个大小为N′1×N′2的匹配矩阵M来表示参考图像和待配准图像中直线编组之间的匹配关系，即其中N′1、N′2均为正整数，N′1为参考图像中提取的直线编组个数，N′2为待配准图像中提取的直线编组个数，参考图像中提取的直线编组<mathsid="math0004"num="0004"><math><![CDATA[<mrow><msub><mi>Z</mi><mn>1</mn></msub><mo>=</mo><mo>{</mo><msubsup><mi>z</mi><mn>1</mn><mn>1</mn></msubsup><mo>,</mo><msubsup><mi>z</mi><mn>1</mn><mn>2</mn></msubsup><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><msubsup><mi>z</mi><mn>1</mn><msubsup><mi>N</mi><mn>1</mn><mo>′</mo></msubsup></msubsup><mo>}</mo><mo>,</mo></mrow>]]></math>id="icf0005"file="A2008100315750002C5.tif"wi="34"he="5"top="200"left="149"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/></maths>待配准图像中提取的直线编组<mathsid="math0005"num="0005"><math><![CDATA[<mrow><msub><mi>Z</mi><mn>2</mn></msub><mo>=</mo><mo>{</mo><msubsup><mi>z</mi><mn>2</mn><mn>1</mn></msubsup><mo>,</mo><msubsup><mi>z</mi><mn>2</mn><mn>2</mn></msubsup><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><msubsup><mi>z</mi><mn>2</mn><msubsup><mi>N</mi><mn>2</mn><mo>′</mo></msubsup></msubsup><mo>}</mo><mo>,</mo></mrow>]]></math>id="icf0006"file="A2008100315750002C6.tif"wi="35"he="5"top="210"left="87"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/></maths>1≤i≤N′，1≤j≤N′；匹配矩阵M满足一一匹配准则，即<mathsid="math0006"num="0006"><math><![CDATA[<mrow><mfencedopen='{'close=''><mtable><mtr><mtd><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><msubsup><mi>N</mi><mn>2</mn><mo>′</mo></msubsup></munderover><msub><mi>M</mi><mi>ij</mi></msub><mo>≤</mo><mn>1</mn><mrow><mo>(</mo><mn>1</mn><mo>≤</mo><mi>i</mi><mo>≤</mo><msubsup><mi>N</mi><mn>1</mn><mo>′</mo></msubsup><mo>)</mo></mrow></mtd></mtr><mtr><mtd><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><msubsup><mi>N</mi><mn>1</mn><mo>′</mo></msubsup></munderover><msub><mi>M</mi><mi>ij</mi></msub><mo>≤</mo><mn>1</mn><mrow><mo>(</mo><mn>1</mn><mo>≤</mo><mi>j</mi><mo>≤</mo><msubsup><mi>N</mi><mn>2</mn><mo>′</mo></msubsup><mo>)</mo></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow>]]></math></maths>第二步，计算参考图像和待配准图像中直线编组之间的可靠性测度，得到可靠性测度矩阵R，R是一个大小为N′1×N′2的矩阵，它的每一项Rij表示直线编组对LPij的可靠性测度，LPij为Z1中的任意一个直线编组z1i与Z2的任意一个直线编组z2j组成的直线编组对，方法是采用D-S证据理论方法，对c1i和c2j进行融合计算，得到Rij，即<mathsid="math0007"num="0007"><math><![CDATA[<mrow><msub><mi>R</mi><mi>ij</mi></msub><mo>=</mo><mfrac><mrow><msubsup><mi>c</mi><mn>1</mn><mi>i</mi></msubsup><mo>×</mo><msubsup><mi>c</mi><mn>2</mn><mi>j</mi></msubsup></mrow><mrow><msubsup><mi>c</mi><mn>1</mn><mi>i</mi></msubsup><mo>×</mo><msubsup><mi>c</mi><mn>2</mn><mi>j</mi></msubsup><mo>+</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msubsup><mi>c</mi><mn>1</mn><mi>i</mi></msubsup><mo>)</mo></mrow><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msubsup><mi>c</mi><mn>2</mn><mi>j</mi></msubsup><mo>)</mo></mrow></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中c1i为直线编组z1i的可靠性，c2j为直线编组z2j的可靠性，它们用下式计算<mathsid="math0008"num="0008"><math><![CDATA[<mrow><msubsup><mi>c</mi><mi>H</mi><mi>k</mi></msubsup><mo>=</mo><mi>exp</mi><mo>[</mo><mo>-</mo><mfrac><msup><mrow><mo>(</mo><msubsup><mover><mi>d</mi><mo>‾</mo></mover><mi>H</mi><mi>k</mi></msubsup><mo>)</mo></mrow><mn>2</mn></msup><mrow><mn>2</mn><msup><msub><mi>σ</mi><mn>2</mn></msub><mn>2</mn></msup></mrow></mfrac><mo>]</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中σ2是高斯函数的方差，取2□4；H＝1时，k＝1，2，…，N′1，H＝2时，k＝1，2，…N′2，<overscore>d</overscore>Hk；dHk为dHk的期望，dHk表示zHk的真实交点QHk和提取交点QHk之间的距离，有<mathsid="math0009"num="0009"><math><![CDATA[<mrow><msubsup><mover><mi>d</mi><mo>‾</mo></mover><mi>H</mi><mi>k</mi></msubsup><mo>=</mo><mfrac><mn>1</mn><msub><mi>S</mi><msub><mi>r</mi><mi>ABCD</mi></msub></msub></mfrac><msub><mo>∫</mo><msub><mi>r</mi><mi>ABCD</mi></msub></msub><mo>|</mo><mo>|</mo><msubsup><mi>O</mi><mi>H</mi><mi>k</mi></msubsup><mo>-</mo><msubsup><mi>Q</mi><mi>H</mi><mi>k</mi></msubsup><mo>|</mo><mo>|</mo><mi>d</mi><msubsup><mi>Q</mi><mi>H</mi><mi>k</mi></msubsup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中，zHk由直线特征lHk1，lHk2组成，过lHk1的中点MHk1的与lHk1夹角均为ΔθHk1的两条虚线是lHk1可能发生的最大方向偏移，过lHk2的中点MHk2的与lHk2夹角均为ΔθHk2的两条虚线是lHk2可能发生的最大方向偏移，四条虚线交于A，B，C和D四点，ΔθHk1和ΔθHk2是直线特征lHk1和lHk2的方向偏移，它们的大小由公式(3)确定，rABCD为A，B，C和D四点围成的四边形，id="icf0011"file="A2008100315750003C4.tif"wi="9"he="6"top="223"left="90"img-content="drawing"img-format="tif"orientation="portrait"inline="yes"/>表示rABCD的面积，其大小根据过lHk1的中点MHk1和过lHk2的中点MHk2的四条虚线方程求出，四条虚线方程通过ΔθHk1和ΔθHk2以及提取的直线特征lHk1、lHk2获得，<mathsid="math0010"num="0010"><math><![CDATA[<mrow><mi>Δ</mi><msubsup><mi>θ</mi><mi>H</mi><mi>m</mi></msubsup><mo>=</mo><mi>αexp</mi><mo>[</mo><mo>-</mo><mfrac><msup><mrow><mo>(</mo><msubsup><mi>length</mi><mi>H</mi><mi>m</mi></msubsup><mo>)</mo></mrow><mn>2</mn></msup><mrow><mn>2</mn><msup><msub><mi>σ</mi><mn>1</mn></msub><mn>2</mn></msup></mrow></mfrac><mo>]</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中α是一个预设门限，用来限定方向偏移ΔθHm的最大范围，取α＝0.175弧度，即直线特征的方向偏移不超过10度；lengthHm是lHm的长度；σ1是高斯函数的方差，它决定多长的直线特征将不发生方向偏移，取20□40，具体取决于图像的分辨率，图像分辨率越高，σ1越靠近40；第三步，计算参考图像和待配准图像中直线编组之间的相似性测度，得到相似性测度矩阵B，B是一个大小为N′1×N′2的矩阵，它的每一项Bij表示参考图像和待配准图像中任意两个直线编组z1i和z2j的相似性测度，相似性测度表示两个直线编组之间的相似程度，相似性测度Bij定义为<mathsid="math0011"num="0011"><math><![CDATA[<mrow><msub><mi>B</mi><mi>ij</mi></msub><mo>=</mo><mo>|</mo><mo>|</mo><mfrac><mrow><msubsup><mi>length</mi><mn>1</mn><mrow><mi>i</mi><mn>1</mn></mrow></msubsup><mo>-</mo><msubsup><mi>length</mi><mn>1</mn><mrow><mi>i</mi><mn>2</mn></mrow></msubsup></mrow><mrow><msubsup><mi>length</mi><mn>1</mn><mrow><mi>i</mi><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mi>length</mi><mn>1</mn><mrow><mi>i</mi><mn>2</mn></mrow></msubsup></mrow></mfrac><mo>|</mo><mo>-</mo><mo>|</mo><mfrac><mrow><msubsup><mi>length</mi><mn>2</mn><mrow><mi>j</mi><mn>1</mn></mrow></msubsup><mo>-</mo><msubsup><mi>length</mi><mn>2</mn><mrow><mi>j</mi><mn>2</mn></mrow></msubsup></mrow><mrow><msubsup><mi>length</mi><mn>2</mn><mrow><mi>j</mi><mn>1</mn></mrow></msubsup><mo>+</mo><msubsup><mi>length</mi><mn>2</mn><mrow><mi>j</mi><mn>2</mn></mrow></msubsup></mrow></mfrac><mo>|</mo><mo>|</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mrow>]]></math></maths>对(9)式进行变换，即Bij＝F1(Bij)(10)其中F1(·)是一个单调递减函数，采用高斯函数，即<mathsid="math0012"num="0012"><math><![CDATA[<mrow><msub><mi>F</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><mfrac><msup><mi>x</mi><mn>2</mn></msup><msubsup><mrow><mn>2</mn><mi>σ</mi></mrow><mn>3</mn><mn>2</mn></msubsup></mfrac><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>11</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中σ3是高斯函数的方差，在实际中一般取σ3＝1，将Bij代入(11)式，得<mathsid="math0013"num="0013"><math><![CDATA[<mrow><msub><mi>B</mi><mi>ij</mi></msub><mo>=</mo><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><mfrac><msup><msub><mi>B</mi><mi>ij</mi></msub><mn>2</mn></msup><mrow><mn>2</mn><msubsup><mi>σ</mi><mn>3</mn><mn>2</mn></msubsup></mrow></mfrac><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mrow>]]></math></maths>第四步，计算参考图像和待配准图像中直线编组之间的空间关系一致性测度，得到空间关系一致性测度矩阵S，S也是大小N′1×N′2的矩阵，它的每一项Sij表示两幅图像中任意两个直线编组z1i和z2j之间的空间关系一致性测度，Sij通过下式计算<mathsid="math0014"num="0014"><math><![CDATA[<mrow><msub><mi>S</mi><mi>ij</mi></msub><mo>=</mo><msqrt><msup><mrow><mo>(</mo><msubsup><mi>x</mi><mn>1</mn><mi>i</mi></msubsup><mo>-</mo><msub><mi>T</mi><mi>x</mi></msub><mrow><mo>(</mo><msubsup><mi>x</mi><mn>2</mn><mi>j</mi></msubsup><mo>)</mo></mrow><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo><msup><mrow><mo>(</mo><msubsup><mi>y</mi><mn>1</mn><mi>i</mi></msubsup><mo>-</mo><msub><mi>T</mi><mi>y</mi></msub><mrow><mo>(</mo><msubsup><mi>y</mi><mn>2</mn><mi>j</mi></msubsup><mo>)</mo></mrow><mo>)</mo></mrow><mn>2</mn></msup></msqrt><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中Tx(·)、Ty(·)分别表示参考图像和待配准图像之间的变换模型在x方向和y方向上的变换函数，参考图像和待配准图像之间的变换模型用变换模型参数向量P来描述，P与采用的变换模型类型有关；采用(11)式的高斯函数对(13)式做变换，即Sij＝F1(Sij)(14)将Sij代入(14)式，得<mathsid="math0015"num="0015"><math><![CDATA[<mrow><msub><mi>S</mi><mi>ij</mi></msub><mo>=</mo><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><mfrac><msup><msub><mi>S</mi><mi>ij</mi></msub><mn>2</mn></msup><mrow><mn>2</mn><msubsup><mi>σ</mi><mn>3</mn><mn>2</mn></msubsup></mrow></mfrac><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>15</mn><mo>)</mo></mrow><mo>;</mo></mrow>]]></math></maths>第五步，建立能量函数<mathsid="math0016"num="0016"><math><![CDATA[<mrow><mi>E</mi><mrow><mo>(</mo><mi>M</mi><mo>,</mo><mi>P</mi><mo>)</mo></mrow><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><msubsup><mi>N</mi><mn>1</mn><mo>′</mo></msubsup></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><msubsup><mi>N</mi><mn>2</mn><mo>′</mo></msubsup></munderover><msub><mi>M</mi><mi>ij</mi></msub><msub><mi>R</mi><mi>ij</mi></msub><msub><mi>B</mi><mi>ij</mi></msub><msub><mi>S</mi><mi>ij</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mrow>]]></math></maths>当匹配矩阵M包含了所有正确匹配的直线编组对而剔除了所有错误匹配的直线编组对，E(·)的值将会达到一个全局最大值，这时的匹配矩阵M是最优的，同时由最优匹配矩阵计算的变换模型参数向量P也是最优的，即<Mopt，Popt>＝argmax(E(·))(17)其中Mopt，Popt分别表示最优的匹配矩阵和最优的变换模型参数向量，argmax(·)表示对·求最大值，这样，直线编组的匹配问题转变成为一个能量函数最优化问题；第六步，采用两步迭代优化算法对(17)式描述的优化问题进行迭代优化，得到Mopt和Popt步骤1采用一个序贯组合特征相似性和空间关系一致性的方法来初始化变换模型参数向量P首先根据特征相似性测度选择求解变换模型需要的最少数目的3-5倍的直线编组对作为候选匹配直线编组对，然后根据空间关系一致性测度利用穷举法或RANSAC算法获得求解变换模型参数向量P需要的最少数目的最优匹配直线编组对，利用它们的提取交点求出变换模型参数向量P的初值P0；步骤2初始化控制参数β＝β0，β0、βr和βmax分别是控制参数β的初始值、增长率和最大值，取β0＝0.5，βr取1.05-1.1，βmax＝10；步骤3已知变换模型参数向量P，采用软指派方法计算匹配矩阵M；步骤4已知匹配矩阵M，空间关系一致性测度矩阵S的每个元素Sij就可表示成变换模型参数向量P的函数，公式(16)变为<mathsid="math0017"num="0017"><math><![CDATA[<mrow><mi>E</mi><mrow><mo>(</mo><mi>P</mi><mo>)</mo></mrow><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><msubsup><mi>N</mi><mn>1</mn><mo>′</mo></msubsup></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><msubsup><mi>N</mi><mn>2</mn><mo>′</mo></msubsup></munderover><msub><mi>M</mi><mi>ij</mi></msub><msub><mi>R</mi><mi>ij</mi></msub><msub><mi>B</mi><mi>ij</mi></msub><msub><mi>S</mi><mi>ij</mi></msub><mo>=</mo><munderover><mi>Σ</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><msubsup><mi>N</mi><mn>1</mn><mo>′</mo></msubsup></munderover><munderover><mi>Σ</mi><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><msubsup><mi>N</mi><mn>2</mn><mo>′</mo></msubsup></munderover><msub><mi>S</mi><mi>ij</mi></msub><mrow><mo>(</mo><mi>P</mi><mo>)</mo></mrow><msub><mi>U</mi><mi>ij</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>21</mn><mo>)</mo></mrow></mrow>]]></math></maths>其中Uij＝MijRijBij是一个常数，采用多元变量变尺度算法计算变换模型参数向量P，然后β＝βrβ；步骤5若控制参数β＞βmax，判断变换模型参数向量P是否收敛，如果收敛，算法结束，如果不收敛，返回步骤2；若控制参数β≤βmax，则转步骤3；变换模型参数向量P的收敛准则为||P0-P1||＜ε1(22)其中P0、P1分别表示前一次和当前获得的变换模型参数向量，ε1是一个预设门限，取0.01-0.05。2.如权利要求1所述的基于直线特征图像配准中的特征匹配方法，其特征在于从参考图像和待配准图像中提取直线编组时只选择那些满足约束条件的直线编组作为配准基元，约束条件是直线编组的提取交点到组成它的两条直线之间的Hausdroff距离之和小于一个预设门限r，r取10-30个象素。3.如权利要求1所述的基于直线特征图像配准中的特征匹配方法，其特征在于所述软指派方法为(2)行归一化匹配矩阵M，即M广M,./f]M"/一i(3)列归一化匹配矩阵M，即My-M,/fx;/'=1(4)当收敛时，算法结束，否则转(2)，匹配矩阵M的收敛准则是SZK-M;l"o(20)其中My°、My'分别表示前一次和当前获得的匹配矩阵M中元素M々的值，s。是一个预设门限，取0.01-0.05。全文摘要本发明公开了一种基于直线特征图像配准中的特征匹配方法，目的是解决基于直线特征图像配准中的特征匹配问题。技术方案是首先从参考图像和待配准图像中提取直线编组作为配准基元，定义匹配矩阵M；接着计算这两幅图像中直线编组之间的可靠性测度、相似性测度、空间关系一致性测度，得到可靠性测度矩阵R、相似性测度矩阵B、空间关系一致性测度矩阵S；然后建立以M为自变量的能量函数，将直线编组的匹配问题转变成能量函数最优化问题；最后采用两步迭代优化算法对能量函数最优化问题进行迭代优化。本发明在特征匹配过程中有机组合了直线编组的可靠性、空间关系一致性和相似性，克服了现有方法存在的计算复杂、容易造成误匹配以及同名像点获取不准确的缺陷。文档编号G06T7/00GK101315698SQ20081003157公开日2008年12月3日申请日期2008年6月25日优先权日2008年6月25日发明者吕金建,文贡坚申请人:中国人民解放军国防科学技术大学