一种面向保积仿射变换的图像配准方法

摘要

本发明涉及一种面向保积仿射变换的图像配准方法，包括：读取待配准具有保积仿射几何扭曲的图像对，从所读取的图像对中提取并建立特征匹配对；其中，所述图像对中的参考图像为主图像，待配准图像为辅图像；基于特征匹配对，利用迭代估计法反演出图像对间的保积仿射扭曲参数；其中，所述保积仿射扭曲参数包括：保积仿射矩阵参数a、b、c、d以及两个图像方向上的偏移量tx和ty；利用保积仿射扭曲函数，对辅图像进行插值处理，使其与主图像在几何上实现对准。

说明书

技术领域

本发明涉及计算机图形图像处理领域，特别涉及一种面向保积仿射变换的图像 配准方法。

背景技术

仿射变换是计算机图像处理领域一种最常用的六自由度刚性几何变换。不同于 一般的四自由度相似变换，仿射变换可将一个矩形图像块变换为一个任意的平行四 边形图像块，这种几何扭曲对于光学图像、遥感图像以及医学影像等都具有广泛的 适用性。因此，面向仿射变换和相似变换扭曲的图像配准问题在过去一段时间内受 到了大量关注，许多成熟的算法得以提出，并已在不同领域得到了普遍应用。

近期随着机器视觉和遥感技术的深入发展，图像处理领域也出现了一些新的研 究动态，即关于五自由度图像几何变换和配准的研究。迄今共有两种不同的五自由 度几何变换被提出，分别是弱仿射变换(参见参考文献1：“D.Li and Y.Zhang,"A novel  approach for the registration of weak affine images,"Pattern Recognition Letters,vol.33, no.12,pp.1647-1655,2012”)和保积仿射变换(参见参考文献2：“J.Flusser and B. Zitova,"A comment on'a novel approach for the registration of weak affine images'," Pattern Recognition Letters,vol.34,no.12,pp.1381-1385,2013”)。两者分别通过对几 何切变和图像块面积的约束，去掉了原始仿射变换中的一个自由度。从参数估计角 度看，这种操作减少了待估参数个数，然而却极大地增加了估计难度，因为此时我 们面临的是一个约束优化问题，其无法像仿射变换和相似变换那样允许我们直接得 到一个严格的解析解。对于这类五自由度扭曲，常用的参数反演方法将不适用且不 精确，因此亟需发展出新的估计方法。

基于一个分步优化方法，面向弱仿射变换的参数反演算法近期被发展出，实验 表明其具有精确的配准效果(参见参考文献3：“D.Li and Y.Zhang,"A novel approach  for the registration of weak affine images,"Pattern Recognition Letters,vol.33,no.12,pp. 1647-1655,2012”)。然而对于面向保积仿射变换的图像配准，尽管Flusser和Zitova 近期给出了解决该问题的初步设想(参见参考文献4：“J.Flusser and B.Zitova,"A  comment on'a novel approach for the registration of weak affine images',"Pattern Recognition Letters,vol.34,no.12,pp.1381-1385,2013”)，但却没能给出确切的实现 方案，并且Flusser和Zitova的设想存在明显的局限性。

发明内容

本发明的目的在于克服现有面向保积仿射变换图像配准方法存在局限性的缺 陷，从而提供一种精确、有效的图像配准方法。

为了实现上述目的，本发明提供了一种面向保积仿射变换的图像配准方法，包 括：

步骤1)、读取待配准具有保积仿射几何扭曲的图像对，从所读取的图像对中提 取并建立特征匹配对；其中，所述图像对中的参考图像为主图像，待配准图像为辅 图像；

步骤2)、基于步骤1)所得到的特征匹配对，利用迭代估计法反演出图像对间 的保积仿射扭曲参数；其中，所述保积仿射扭曲参数包括：保积仿射矩阵参数a、b、 c、d以及两个图像方向上的偏移量t_x和t_y；

步骤3)、利用步骤2)得到的保积仿射扭曲函数，对辅图像进行插值处理，使 其与主图像在几何上实现对准。

上述技术方案中，在步骤1)中，利用ASIFT算子在主图像与辅图像间建立特 征匹配对，每一对特征匹配的表达式为其中，

$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right),s.t.|\mathrm{ad}-\mathrm{bc}|=1.$

上述技术方案中，所述步骤2)进一步包括：

步骤2-1)、根据步骤1)所得到的特征匹配对计算A参数和B参数，计算公式 如下：

$\left(\begin{array}{c}{A}_1=\frac{1}{N}\underset{i=1}{\overset{n}{\Sigma }}x,A_2=\frac{1}{N}\underset{i=1}{\overset{n}{\Sigma }}y,A_3=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{x}^{\prime },A_4=\frac{1}{4}\underset{i=1}{\overset{N}{\Sigma }}{y}^{\prime }\\ B_1=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{x}^2,B_2=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{y}^2,B_3=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\mathrm{xy},B_4=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{xx}}^{\prime }\\ B_5=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{yx}}^{\prime },B_6=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{xy}}^{\prime },B_7=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{yy}}^{\prime },B_8=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{x}^{\prime 2},B_9=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{y}^{\prime 2}\end{array}\right);$

其中，N为特征匹配对的数目；

步骤2-2)、根据A参数和B参数，计算C参数和D参数，计算公式如下：

$\left(\begin{array}{c}{C}_1=B_1-A_1^2,C_2=B_2-A_2^2,C_3=B_3-A_1{A}_2\\ D_1=A_1{A}_3-B_4,D_2=A_2{A}_3-B_5,D_3=A_1{A}_4-B_6,D_4=A_2{A}_4-B_7\end{array}\right);$

步骤2-3)、根据步骤2-1)和步骤2-2)所得到的A参数、B参数、C参数和D 参数，利用迭代估计法反演保积仿射扭曲参数a、b、c、d、t_x和t_y。

上述技术方案中，所述步骤2-3)包括：

采用内循环自动参数更新法计算出內自新最优迭代次数Iter_in以及内自新最小保 积误差E_in；然后在外循环自动参数更新法中采用内循环自动参数更新法所得到的计 算结果，计算出最优迭代次数Iter_opt与最优κ参数κ_opt；最后根据最优迭代次数Iter_opt与最优κ参数κ_opt采用最优自动参数更新法计算出保积仿射参数a、b、c、d、t_x和t_y； 其中，

内循环自动参数更新法输入：內自新最大迭代次数Iter_max；参数κ；参数F；参 数a、b、c和d；

内循环自动参数更新法输出：內自新最优迭代次数Iter_in；内自新最小保积误差 E_in；

内循环自动参数更新法执行步骤如下：

初始化内自新次数Iter＝1；

步骤101：若Iter≤Iter_max，更新参数a、b、c和d，利用更新后的参数a、b、c 和d计算保积误差E，利用更新后的参数a、b、c和d更新参数F；否则，跳转至步 骤103；其中，

更新参数a、b、c和d的公式为：

$\left(\begin{array}{c}a=\frac{-C_2{D}_1+C_3{D}_2+D_{4}F}{C_1{C}_2-C_3^2-F^2},b=\frac{-C_1{D}_2+C_3{D}_1-D_{3}F}{C_1{C}_2-C_3^2-F^2}\\ c=\frac{-C_2{D}_3+C_3{D}_4-D_{2}F}{C_1{C}_2-C_3^2-F^2},d=\frac{-C_1{D}_4+C_3{D}_3+D_{1}F}{C_1{C}_2-C_3^2-F^2}\end{array}\right);$

计算保积误差E的计算公式为：

E＝((ad-bc)²-1)²；

更新参数F的公式为：

F＝2×10^κ(ad-bc)((ad-bc)²-1)；

步骤102：更新迭代次数Iter＝Iter+1，跳转至步骤101；

步骤103：将Iter_max次迭代得到的E中最小值作为E_in，其对应的迭代次数作为 Iter_in；

外循环自动参数更新法输入：参数κ最大取值κ_max和更新步长κ_Δ；

外循环自动参数更新法输出：最优迭代次数Iter_opt；最优κ参数κ_opt；

外循环自动参数更新法执行步骤如下：

初始化內自新最大迭代次数Iter_max；初始化参数F＝0；初始化参数κ＝κ_min，κ_min<κ_max；

步骤201：执行內自新，获得內自新最优迭代次数Iter_in及最小保积误差E_in；

步骤202：更新参数κ：κ＝κ+κ_Δ，利用更新后的参数κ更新参数F，若κ≤κ_max， 跳转至步骤201；否则，将所有外自新迭代得到的E_in中最小值对应的迭代次数作为 最优迭代次数Iter_opt，对应的κ参数作为最优κ参数κ_opt；其中，

更新参数F的公式为：

F＝2×10^κ(ad-bc)((ad-bc)²-1)；

最优自动参数更新法输入：参数κ；参数F；最优自新最大迭代次数Iter_opt；

最优自动参数更新法输出：保积仿射参数a、b、c、d、t_x和t_y；

最优自动参数更新法执行步骤如下：

初始化参数F＝0；初始化最优自新次数Iter＝1；设置κ参数为κ_opt；

步骤301：若Iter≤Iter_opt，计算参数a、b、c和d，利用更新后的参数a、b、c 和d更新参数F；否则，跳转至步骤303；其中，

计算参数a、b、c和d的公式为：

$\left(\begin{array}{c}a=\frac{-C_2{D}_1+C_3{D}_2+D_{4}F}{C_1{C}_2-C_3^2-F^2},b=\frac{-C_1{D}_2+C_3{D}_1-D_{3}F}{C_1{C}_2-C_3^2-F^2}\\ c=\frac{-C_2{D}_3+C_3{D}_4-D_{2}F}{C_1{C}_2-C_3^2-F^2},d=\frac{-C_1{D}_4+C_3{D}_3+D_{1}F}{C_1{C}_2-C_3^2-F^2}\end{array}\right);$

更新参数F的公式为：

F＝2×10^κ(ad-bc)((ad-bc)²-1)；

步骤302：更新迭代次数Iter＝Iter+1，跳转至步骤301；

步骤303：计算参数t_x和t_y；其计算公式为：

t_x＝A₃-A₁a-A₂b,t_y＝A₄-A₁c-A₂d。

本发明的优点在于：

本发明的方法能够精确有效地实现对具有保积仿射扭曲图像对的几何配准。

附图说明

图1是本发明的面向保积仿射变换的图像配准方法的流程图；

图2是待配准保积仿射扭曲图像对的示意图；其中左图为主图像，右图为辅图 像；

图3是从图2提取的有效特征匹配对的示意图；

图4是本发明方法中所涉及的内循环自动参数更新法的流程图；

图5是本发明方法中所涉及的外循环自动参数更新法的流程图；

图6是本发明方法中所涉及的最优自动参数更新法的流程图；

图7是最终配准结果的示意图；其中左图为原始主图像，右图为对准后的辅图 像。

具体实施方式

现结合附图对本发明作进一步的描述。

本发明的面向保积仿射变换的图像配准方法利用拉格朗日乘子法将保积仿射估 计的约束优化问题转化为一个迭代优化问题，提出一个全新的面向保积仿射变换的 目标函数，推导出严格的参数反演公式，设计出高效的迭代估计方法，从而精确有 效地实现了对扭曲图像对的几何配准。

本发明的面向保积仿射变换的图像配准方法在实现几何配准时，需要在待配准 的图像对中提取和构建具有一定规模、一定匹配精度、一定几何定位精度和具有一 定几何不变性的特征匹配对。一定规模是指提取的特征匹配对数目应不小于待求解 模型的自由度。对于保积仿射变换来说，这要求至少提取出五对有效特征匹配，否 则较少的观测样本无法支持我们得到唯一解。一定匹配精度是指允许得到的特征匹 配对中存在错匹配，但错匹配数目有一定限制。对于常用的稳健LMS和LTS估计， 其要求有效匹配数目至少等于错匹配个数，对于RANSAC估计，其要求有效匹配数 目最少应与待估模型的自由度相当。本发明不考虑存在错匹配情形，因为只要将提 出的参数反演算法嵌入到LMS、LTS或RANSAC等稳健估计中，该问题即可得到有 效解决。在常用的稳健特征提取算子中，提取的特征有几何位置和描述符两方面描 述信息。几何位置的定位精度直接影响参数估计精度，因为参数估计直接使用匹配 对的位置信息。故提取特征的定位精度应至少优于最终图像的配准精度，即要得到 亚像素级配准精度，那么特征匹配的几何定位也需优于亚像素。描述符通常是一串 具有特定几何不变性的矢量，其记录了特征周围的纹理信息，可用于特征匹配的构 建。描述符的几何不变性要不低于待估模型，这才能保证得到的特征匹配能与待估 模型相吻合。在本发明中，描述符几何不变性要不低于保积仿射变换。

基于对特征匹配对的上述要求，在本申请中采用ASIFT算子在待配准的两幅图 像之间建立特征匹配对，该ASIFT算子的定位精度为亚像素级，对一般仿射变换具 有几何不变性。

在成功构建好特征匹配点对后，即可开启参数估计流程。假设为任意一对 待配准保积仿射扭曲图像，其中I为主图像，I'为辅图像，若为建立 的N组特征匹配之一，保积仿射变换可写为：

$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=A\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)---(1-1)$

其中，A为2×2保积仿射矩阵，其由参数a、b、c和d构成，t_x和t_y为两个图像 方向上的偏移量。矩阵A满足约束：

|det(A)|＝|ad-bc|＝1          (1-2)

图像配准的目的是基于构建的特征匹配对反演保积仿射矩阵参数a、b、c、d以 及两个图像方向上的偏移量t_x和t_y。由上式可看到，此时面临的是一个约束优化问题。 利用拉格朗日乘子法，上述问题可转化为对下述优化问题的求解：

$O=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\left|\right|{\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)-A\left(\begin{array}{c}x\\ y\end{array}\right)-\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)\left|\right|}_2^2+\lambda {\left(\right|\det A|-1)}^2,\lambda \ge 0---\left(2\right)$

式中，λ为拉格朗日乘因子。为了便于计算，Flusser和Zitova建议将上述问题 进一步放松为：

$O=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\left|\right|{\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)-A\left(\begin{array}{c}x\\ y\end{array}\right)-\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)\left|\right|}_2^2+\lambda {\left(\right|\det A|-1)}^2,\lambda \ge 0---\left(3\right)$

相比于式(2)，式(3)明显忽视了det A<0的情况，这限制该估计的适用范围， 即使如此，Flusser和Zitova也没有给出确切的求解式(3)的方法，因此尚缺乏一种 全面有效的面向保积仿射变换的参数估计方法。本发明中，申请人提出了一个新的 优化目标问题：

$O=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\left|\right|{\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)-A\left(\begin{array}{c}x\\ y\end{array}\right)-\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)\left|\right|}_2^2+\lambda {({(\det \mathrm{><mo>)</mo>}}^2-1)}^2,\lambda \ge 0---\left(4\right)$

式(4)与式(2)完全等价，解决了式(3)的不足。进一步带入参数a、b、c 和d有：

$O=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}({[x^{\prime }-(\mathrm{ax}+\mathrm{by}+t_x)]}^2+{[y^{\prime }-(\mathrm{cx}+\mathrm{dy}+t_y)]}^2)+\lambda {({(\mathrm{ad}-\mathrm{bc})}^2-1)}^2---\left(5\right)$

展开、整理并化简，有如下表达式：

$\left(\begin{array}{c}O=B_8+B_9+(a^2+c^2)B_1+(b^2+d^2)B_2+2(\mathrm{ab}+\mathrm{cd})B_3+(t_x^2+t_y^2)+2({\mathrm{at}}_x+{\mathrm{ct}}_y)A_1\\ +2({\mathrm{bt}}_x+d{t}_y)A_2-2(t_x{A}_3+t_y{A}_4+a{B}_4+b{B}_5+c{B}_5+c{B}_6+{\mathrm{dB}}_7)+\lambda {({(\mathrm{ad}-\mathrm{cb})}^2-1)}^2\end{array}\right)---\left(6\right)$

其中：

$\left(\begin{array}{c}{A}_1=\frac{1}{N}\underset{i=1}{\overset{n}{\Sigma }}x,A_2=\frac{1}{N}\underset{i=1}{\overset{n}{\Sigma }}y,A_3=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{x}^{\prime },A_4=\frac{1}{4}\underset{i=1}{\overset{N}{\Sigma }}{y}^{\prime }\\ B_1=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{x}^2,B_2=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{y}^2,B_3=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\mathrm{xy},B_4=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{xx}}^{\prime }\\ B_5=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{yx}}^{\prime },B_6=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{xy}}^{\prime },B_7=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{yy}}^{\prime },B_8=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{x}^{\prime 2},B_9=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{y}^{\prime 2}\end{array}\right);$

式(6)两边分别对参数a、b、c、d、t_x和t_y取偏导，并令其为0，有：

$\left(\begin{array}{c}\frac{\partial O}{\partial a}=2(B_{1}a+B_{3}b+A_1{t}_x-B_4)+2\mathrm{Fd}=0\\ \frac{\partial O}{\partial b}=2(B_{3}a+B_{2}b+A_2{t}_x-B_5)-2\mathrm{Fc}=0\\ \frac{\partial O}{\partial c}=2(B_{1}c+B_{3}d+A_1{t}_y-B_6)-2\mathrm{Fb}=0\\ \frac{\partial O}{\partial d}=2(B_{3}c+B_{2}d+A_2{t}_y-B_7)+2\mathrm{Fa}=0\\ \frac{\partial O}{\partial t_x}=2(A_{1}a+A_{2}b+t_x-A_3)=0\\ \frac{\partial O}{\partial t_y}=2(A_{1}c+A_{2}d+t_y-A_4)=0\end{array}\right)---\left(8\right)$

其中：

F＝2λ(ad-bc)((ad-bc)²-1),λ≥0         (9)

进一步令：

λ＝10^κ        (10)

则式(9)可进一步写为：

F＝2×10^κ(ad-bc)((ad-bc)²-1)        (11)

求解式(8)中最后二式，有：

t_x＝A₃-A₁a-A₂b,t_y＝A₄-A₁c-A₂d       (12)

将式(12)带入式(8)中前四式，进一步化简，有：

$\left(\begin{array}{c}(B_1-A_1^2)a+(B_3-A_1{A}_2)b+(A_1{A}_3-B_4)+\mathrm{Fd}=0\\ (B_3-A_1{A}_2)a(B_2-A_2^2)b+(A_2{A}_3-B_5)-\mathrm{Fc}=0\\ (B_1-A_1^2)c+(B_3-A_1{A}_2)d+(A_1{A}_4-B_6)-\mathrm{Fb}=0\\ (B_3-A_1{A}_2)c+(B_2-A_2^2)d+(A_2{A}_4-B_7)+\mathrm{Fa}=0\end{array}\right)---\left(13\right)$

进一步令：

$\left(\begin{array}{c}{C}_1=B_1-A_1^2,C_2=B_2-A_2^2,C_3=B_3-A_1{A}_2\\ D_1=A_1{A}_3-B_4,D_2=A_2{A}_3-B_5,D_3=A_1{A}_4-B_6,D_4=A_2{A}_4-B_7\end{array}\right)---\left(14\right)$

则式(13)可进一步写为：

$\left(\begin{array}{c}{C}_{1}a+C_{3}b+\mathrm{Fd}+D_1=0\\ C_{3}a+C_{2}b-\mathrm{Fc}+D_2=0\\ C_{1}c+C_{3}d-\mathrm{Fb}+D_3=0\\ C_{3}c+C_{2}d+\mathrm{Fa}+D_4=0\end{array}\right)---\left(15\right)$

求解式(15)有如下公式：

$\left(\begin{array}{c}a=\frac{-C_2{D}_1+C_3{D}_2+D_{4}F}{C_1{C}_2-C_3^2-F^2},b=\frac{-C_1{D}_2+C_3{D}_1-D_{3}F}{C_1{C}_2-C_3^2-F^2}\\ c=\frac{-C_2{D}_3+C_3{D}_4-D_{2}F}{C_1{C}_2-C_3^2-F^2},d=\frac{-C_1{D}_4+C_3{D}_3+D_{1}F}{C_1{C}_2-C_3^2-F^2}\end{array}\right)---\left(16\right)$

式(11)、(12)和(16)构成了为优化求解算法的迭代基础。利用这些公式， 在本申请中可精确估计保积仿射参数a、b、c、d、t_x和t_y，利用这些参数对辅图像进 行插值处理，可将辅图像与主图像在几何上实现精确对准。

以上是对本发明的面向保积仿射变换的图像配准方法中如何选取特征匹配对、 如何估计保积仿射参数的背景说明，下面结合实例对本发明方法的具体实现步骤做 进一步的说明。

参考图1，本发明的面向保积仿射变换的图像配准方法包括以下步骤：

步骤1)、读取待配准具有保积仿射几何扭曲的图像对，提取并建立有效特征匹 配对；

步骤2)、基于得到的特征匹配对，利用设计的迭代估计算法精确反演出图像对 间的保积仿射扭曲参数；

步骤3)、利用得到的保积仿射扭曲函数，对辅图像进行插值处理，使其与主图 像在几何上实现对准。

在步骤1)中，首先读入两幅待配准保积仿射扭曲图像。在一个实施例中，所读 入的两幅待配准保积仿射扭曲图像如图2所示，其中左图为主图像I(768×770)； 右图为辅图像I'(1055×1124)。I'是将I按照表1实际值所示的保积仿射扭曲参数 变换而成，如果为一对特征匹配，那么：

$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right),s.t.|\mathrm{ad}-\mathrm{bc}|=1---\left(17\right)$

利用ASIFT算子在两幅图像间建立特征匹配共得到了715组 有效匹配(即N＝715)，如图3所示。

在步骤2)中，由步骤1)所得到的有效特征匹配对，利用迭代估计方法精确反 演图像间的保积仿射扭曲参数。该步骤进一步包括：

步骤2-1)、首先按照公式(7)计算A参数和B参数；

步骤2-2)、根据A参数和B参数，参照公式(14)计算C参数和D参数；

步骤2-3)、根据步骤2-1)和步骤2-2)所得到的A参数、B参数、C参数和D 参数，利用迭代方法反演保积仿射扭曲参数a、b、c、d、t_x和t_y。

反演保积仿射扭曲参数a、b、c、d、t_x和t_y的迭代过程由内循环自动参数更新 法(简称內自新)、外循环自动参数更新法(简称外自新)和最优自动参数更新法(简 称最优自新)三部分组成，其实现流程分别如图4、图5和图6所示。

內自新输入：內自新最大迭代次数Iter_max；参数κ；参数F；参数a、b、c和d；

內自新输出：內自新最优迭代次数Iter_in；内自新最小保积误差E_in；

內自新执行步骤如下：

初始化内自新次数Iter＝1；

步骤101：若Iter≤Iter_max，利用式(16)更新参数a、b、c和d，基于式(18) 计算保积误差E，利用式(11)更新参数F；否则，跳转至步骤103；

E＝((ad-bc)²-1)²           (18)

步骤102：更新迭代次数Iter＝Iter+1，跳转至步骤101；

步骤103：将Iter_max次迭代得到的E中最小值作为E_in，其对应的迭代次数作为 Iter_in。

外自新输入：参数κ最大取值κ_max和更新步长κ_Δ；

外自新输出：最优迭代次数Iter_opt；最优κ参数κ_opt；

外自新执行步骤如下：

初始化內自新最大迭代次数Iter_max；初始化参数F＝0；初始化参数κ＝κ_min(κ_min<κ_max)；

步骤201：执行內自新，获得內自新最优迭代次数Iter_in及最小保积误差E_in；

步骤202：更新参数κ：κ＝κ+κ_Δ，利用式(11)更新参数F，若κ≤κ_max，跳转 至步骤201；否则，将所有外自新迭代得到的E_in中最小值对应的迭代次数作为最优 迭代次数Iter_opt，对应的κ参数作为最优κ参数κ_opt。

最优自新输入：参数κ；参数F；最优自新最大迭代次数Iter_opt；

最优自新输出：保积仿射参数a、b、c、d、t_x和t_y；

最优自新执行步骤如下：

初始化参数F＝0；初始化最优自新次数Iter＝1；设置κ参数为κ_opt；

步骤301：若Iter≤Iter_opt，利用式(16)计算参数a、b、c和d，利用式(11) 更新参数F；否则，跳转至步骤303；

步骤302：更新迭代次数Iter＝Iter+1，跳转至步骤301；

步骤303：利用式(12)计算参数t_x和t_y。

由上述迭代过程最终得到的保积仿射参数如表1中估计值所示。通过与这些参 数的实际值相比较可看到，本发明可高精度反演出两幅图像间的保积仿射变换参数， 实际值和估计值间只存在细微的差异，这验证了设计方案的有效性。

在步骤3)中，利用反演出的参数对辅图像进行插值处理，本申请中采用一般的 双线性插值。图7左图所示为原始主图像，右图所示为变换后的辅图像，可看到， 两者在几何上完美地实现了对准，这进一步说明了该方法的有效性。提出的方案可 简单嵌入到现有的稳健估计算法例如LMS、LTS和RANSAC中，使其可有效用于存 在大量错匹配的图像配准情形。

表1

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管 参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明 的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均 应涵盖在本发明的权利要求范围当中。