CN102201119A - 一种基于控制点无偏变换的图像配准方法及*** - Google Patents

Abstract

本发明适用于图像处理领域，提供了一种图像配准方法，包括以下步骤：根据待配准数字图像中的当前控制点对集合构造当前变换函数，并根据该变换函数确定新增源控制点以及新增目标控制点的初始位置；根据变换的无偏性和局部配准度量的综合度量，在新增目标控制点的局部区域内确定与新增源控制点对应的目标控制点；将新增源控制点和新增目标控制点暂时加入到当前控制点对集合中，并计算全局图像的配准度量；根据图像配准度量的变化量确定该新增控制点对加入到当前控制点对集合中的概率值，若概率值大于随机产生的概率参考值，则将新增控制点对添加到当前控制点对集合中，重复上述过程直至图像配准精度达到要求或没有新增控制点。本发明有效地提高标志点对应关系的精度，同时得到无偏性较好的变换函数。

Description

一种基于控制点无偏变换的图像配准方法及***

技术领域

本发明属于图像处理技术领域，尤其涉及一种基于控制点无偏变换的图像配准方法及***。

背景技术

图像配准的目的是寻找两幅图像之间的变换，使得图像中相应目标的位置对准。基于标志点对应关系的配准方法是图像弹性配准中一类重要的方法，但由于其变换函数采用参数模型，不能保证是可逆的，因此不适用于大形变配准问题。

为了使基于标志点对应关系的配准方法中的变换函数是可逆的，常用的方法是构造目标函数，通过对变换函数增加约束条件，可以得到可逆的变换函数。常用的方法如对B样条构造的变换函数系数进行约束，以确保变换函数可逆。如果变换函数的Jacobian行列式值恒为正，则变换函数可以保证是可逆的。这些方法构造的变换函数仅仅是可逆的，不能保证是无偏的。如果正向和反向变换函数的Jacobian值在1附近对称分布，则称为无偏变换。无偏变换对于图像配准非常重要，它意味着无论图像A向图像B配准，或是图像B向图像A配准都可以使同一目标得到配准。目前关于无偏变换的研究比较少，Gholipour定义的对称代价函数，使源图像和目标图像互换后的变换函数互逆，但该方法采用FFD形变模型，计算量较大。Johnson构造了标志点的一致性配准方法，该方法可以使标志点在正、反向的变换达到一致，但会导致标志点对应关系发生较大误差，还需要引入一致性强度配准进一步调整。Leow定义了变换函数与理想变换之间的对称KL距离，并定义了对称目标函数寻找无偏变换函数，但Leow采用的是粘性流体配准模型，计算量仍然很大。

在保证标志点对应关系前提下，构造基于基函数扩展的无偏变换是非常吸引人的，因为这类变换函数可以保证特征点之间的对应关系，同时具有明确的解析式，易于分析，计算简单，计算量小。但存在以下困难：(1)寻找精确的标志点对应关系比较困难。标志点对应关系往往受局部极值的影响而出现误差，对配准结果产生严重影响；(2)构造无偏变换比较困难，因为基函数扩展的变换函数本身不能保证可逆性，不适用于大形变配准问题。因此很少有文献讨论如何构造基于标志点对应关系的可逆变换，Miller在保证标志点对应关系的前提下构造了微分同胚函数，但该方法采用了光流模型，需要求解偏微分方程，计算复杂，而研究基于标志点对应关系的无偏变换就更少了。

构造基于标志点对应关系的无偏变换需要解决两个问题：(1)寻找对应的标志点。目前已经提出了许多确定标志点对应关系的方法，具有代表行性的是基于信息论的方法，这些方法都需要在搜索最优的匹配点，如果能提供较好的初始搜索位置，将对提高匹配精度大有帮助。(2)在保证标志点对应关系的前提下，构造基于基函数扩展的无偏变换。如果要保证标志点对应关系，基函数扩展构造的变换函数形式是确定的，其无偏性无法保证，只有通过调整标志点分布的位置，或是在适当的位置添加标志点，可以调整变换函数的无偏性。

发明内容

本发明所要解决的技术问题在于提供一种基于控制点无偏变换的图像配准方法，旨在解决现有技术应用基于标志点的图像配准插值时存在的标志点对应关系难以确定，同时变换函数无偏性较差的问题。

本发明是这样实现的，一种基于控制点无偏变换的图像配准方法，所述方法包括下述步骤：

步骤A，根据待配准数字图像中的当前控制点对集合构造当前变换函数；所述当前控制点对集合包含初始源控制点集合和初始目标控制点集合；

步骤B，根据步骤A中构造得到的当前变换函数在当前源控制点周围选择新增源控制点，并将新增源控制点的映射位置作为新增目标控制点的初始位置；

步骤C，若步骤B中能找到新增源控制点，则根据变换的无偏性和局部配准度量，调整新增目标控制点；

步骤D，将新增源控制点和新增目标控制点暂时加入到当前控制点对集合中，并计算新变换函数的对称KL距离，以及图像全局配准度量；

步骤E，根据图像全局配准度量的变化量，以及新变换函数的对称KL距离，确定新增控制点对加入到当前控制点对集合中的概率值；所述新增控制点对包含相对应的新增源控制点和新增目标控制点；

步骤F，若步骤E中计算得到概率值大于随机产生的概率参考值，则将新增控制点对添加到当前控制点对集合中；

重复步骤A至步骤F。

本发明还提供了一种基于控制点无偏变换的图像配准***，包括；

当前变换函数构造单元，用于根据待配准数字图像中的当前控制点对集合构造当前变换函数；所述当前控制点对集合包含初始源控制点集合和初始目标控制点集合；

新增源控制点选择单元，用于根据所述当前变换函数构造单元构造得到的当前变换函数在当前源控制点周围选择新增源控制点，并将新增源控制点的映射位置作为新增目标控制点的初始位置；

新增目标控制点确定单元，若所述新增源控制点选择单元能找到新增源控制点，则用于根据变换的无偏性和局部配准度量，调整新增目标控制点；

概率计算单元，用于将新增源控制点和新增目标控制点暂时加入到当前控制点对集合中，并计算新变换函数的对称KL距离，以及图像全局配准度量；并根据图像全局配准度量的变化量，以及新变换函数的对称KL距离，确定新增控制点对加入到当前控制点对集合中的概率值；所述新增控制点对保包含相对应的新增源控制点和新增目标控制点；

当前控制点对集合更新单元，若所述概率计算单元中计算得到概率值大于随机产生的概率参考值，则将新增控制点对添加到当前控制点对集合中。

本发明针对标志点对应关系，提出了一种基于无偏变换的渐进式图像配准方法，该方法只需要很少的几个初始标志点，通过迭代不断加入新的标志点逐渐提高配准精度。该方法计算简单，不但可以用于解决小形变配准问题，也适用于大形变图像弹性配准问题。

附图说明

图1是本发明实施例提供的图像配准方法的流程图；

图2是本发明实施例提供的获取新增目标点初始位置的流程图；

图3是本发明实施例提供的确定新增对应控制点对的流程图；

图4是本发明实施例提供的以模拟退火概率模型接收新增控制点对的流程图；

图5是本发明实施例提供的基于控制点无偏变换的图像配准***的逻辑结构原理图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

在本发明实施例中，根据数量很少的初始源控制点和目标控制点集合，构造初始变换函数。在初始源控制点附近确定新增源控制点，根据初始变换函数获取新增目标控制点的初始位置，并确定新增目标控制点的搜索范围。在新增目标控制点的搜索范围中，根据变换函数的无偏性和局部配准度量，获取新增目标控制点的准确位置；将新增源控制点和新增目标控制点暂时添加到当前控制点对集合中，计算在该扩展控制点对集合下的配准图像，并计算图像的配准度量。根据该配准度量和当前变换函数的无偏性确定模拟退火概率，以概率接收该标志点对加入到当前控制点对集合中。重复上述过程，直至图像配准精度达到要求为止。

图1示出了本发明实施例提供的图像配准的流程，详述如下：

在步骤S101中，根据当前控制点对集合，构造当前变换函数。

在本发明实施例中，待插值处理的图像为数字图像，初始控制点对已知且数量很少(一般少于10对)。假设初始源控制点集合为P＝[p₁，p₂，…，p_N]^T，初始目标控制点集合为Q＝[q₁，q₂，…，q_N]^T，构造的初始函数为h(x)，

$h\left(x\right)=x+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i}U\left(\right||x-p_i|\left|\right)---\left(1\right)$

其中x∈R²是图像中的任意点，R为实数集合，p_i，i＝1，…，N是源控制点，U(r)是径向基函数，系数w_i，i＝1，…，N由方程(2)求解得。

KW＝Q-P        (2)

K是N×N矩阵，K_ij＝U(||p_i-p_j||)，W＝[w₁，w₂，…，w_N]^T。初始控制点对数量很少，一般取小于10对。这些控制点对可以通过人工方式设定，最好均匀分布在整个图像区域，其对应关系一般比较准确。

在步骤S102中，根据步骤S101中构造得到的当前变换函数在当前源控制点周围选择新增源控制点，确定新增目标控制点的初始位置。

本实施例中，根据初始控制点对集合可以构造基于径向基函数的当前变换函数，在当前源控制点附近确定新增源控制点，并根据当前变换函数确定新增目标控制点的初始位置。

当前变换函数为h(x)，其Jacobian行列式值为Dh(x)。在点x处的Jacobian行列式的计算方法如下：

$\mathrm{Dh}\left(x\right)=\left|\begin{array}{cc}\frac{{\∂h}_x\left(x\right)}{\∂x}& \frac{{\∂h}_x\left(x\right)}{\∂y}\\ \frac{{\∂h}_y\left(x\right)}{\∂x}& \frac{{\∂h}_y\left(x\right)}{\∂y}\end{array}\right|=\frac{{\∂h}_x}{\∂x}\·\frac{{\∂h}_y}{\∂y}-\frac{\∂h_x}{\∂y}\·\frac{{\∂h}_y}{\∂x}---\left(3\right)$

选择新增源控制点p_a为

$p_a=\underset{x}{\arg }\max |\mathrm{Dh}\left(x\right)-1|,$ 且distT_min＜||x-p_i||＜distT_max，x∈E， $\∀p_i\∈P$ (4)

满足式(3)的新增源控制点p_a具有以下特点：(1)该点的Jacobian行列式值与1相差较大，即在该点附近变换函数的无偏性较差；(2)该点与当前标志点的距离不能太远(距离小于最大距离阈值distT_max)，即在当前标志点附近选择，同时也不能距离过远(距离大于最小距离阈值distT_min)，以免初始位置误差过大，影响搜索的有效性；(3)该点落在图像边缘上，E是图像边缘点集合。

目标控制点q_a的初始位置为：

$q_a^0=h\left(p_a\right).$

在步骤S102中，如果找不到合适的新增源控制点p_a，则算法结束。否则，继续步骤S103。

在步骤S103中，根据变换的无偏性和局部配准度量，调整新增目标控制点。具体实现流程如图2所示。

在步骤S201中，在以新增目标控制点的初始位置为中心的局部区域中，选取落在图像边缘上的点作为备选目标控制点，将新增源控制点和备选目标控制点临时添加到控制点集合，计算临时变换函数，具体步骤如S101，这里不再赘述。

在步骤S202中，计算临时变换函数的对称KL距离。

假设h′是将临时点对(p_a，x)添加到集合(P，Q)后得到的临时变换函数。令P→Q的正向映射为h′，反向映射为h′^-1。记h′和h′^-1的Jacobian行列式分别为|Dh′(x)|和|Dh′^-1(x)|。定义三个概率密度函数：

${\mathrm{pdf}}_h\left(\mathrm{\ξ}\right)=\left|{\mathrm{Dh}}^{\′}\left(\mathrm{\ξ}\right)\right|,{\mathrm{pdf}}_{h^{-1}}\left(\mathrm{\ξ}\right)=\left|{\mathrm{Dh}}^{\′-1}\left(\mathrm{\ξ}\right)\right|,{\mathrm{pdf}}_{\mathrm{id}}\left(\mathrm{\ξ}\right)=1---\left(5\right)$

其中pdf_id是单位映射(id(x)＝x)的概率密度函数。KL距离可以描述两个概率密度函数之间的相似性，KL(pdf_h′，pdf_id)是正向变换与单位变换之间的KL距离，

是反向变换与单位变换之间的KL距离，定义变换h′的对称KL距离sKL(h′)为

$\mathrm{sKL}\left(h^{\′}\right)=\mathrm{KL}({\mathrm{pdf}}_{h^{\′}},{\mathrm{pdf}}_{\mathrm{id}})+\mathrm{KL}({\mathrm{pdf}}_{h^{\′-1}},{\mathrm{pdf}}_{\mathrm{id}})$

$=\∫\left(\right|{\mathrm{Dh}}^{\′}\left(x\right)-1\left|\right)\log \left|{\mathrm{Dh}}^{\′}\left(x\right)\right|\mathrm{dx}---\left(6\right)$

$=\∫\left(\right|{\mathrm{Dh}}^{\′-1}\left(x\right)-1)\log \left|{\mathrm{Dh}}^{\′-1}\left(x\right)\right|\mathrm{dx}$

在步骤S203中，根据临时变换函数h′计算新增源控制点p_a所处局部区域Loc(p_a)和新增目标控制点q_a所处局部区域Loc(q_a)的配准度量。

$\mathrm{MSD}(\mathrm{Loc}\left(q_a\right),\mathrm{Loc}\left(p_a\right))=\underset{(i^{\′},j^{\′})\∈\mathrm{Loc}\left(p_a\right)}{\underset{(i,j)\∈\mathrm{Loc}\left(q_a\right)}{\mathrm{\Σ}}}|f(i,j)-f(i^{\′},j^{\′})|---\left(7\right)$

其中，(i，j)表示参考图像中的点，(i′，j′)表示待配准图像中的点，f(i，j)是参考图像在点(i，j)处的强度值，f(i′，j′)是待配准图像在点(i′，j)处的强度值。

在步骤S204中，根据临时变换函数h′的对称KL距离和图像局部配准度量，确定最优的目标控制点作为新增源控制点的对应点。

在局部区域

中寻找p_a的对应点q_a。最优目标控制点q_a为

其中h′是将临时点对(p_a，x)添加到集合(P，Q)后得到的变换函数，Sοh′是经过h′变换后的图像，MSD是均方误差，sKL是对称KL距离。

在步骤S103中，将新增源控制点和新增目标控制点暂时加入到当前控制点对集合中，并计算新变换函数的对称KL距离，以及全局配准度量。具体流程如图3所示。

在步骤S301中，将新增源控制点和新增目标控制点临时添加到控制点集合，计算临时变换函数，具体步骤如S101，这里不再赘述。

在步骤S302中，计算临时变换函数的对称KL距离，具体步骤如S202，这里不再赘述。

在步骤S303中，根据临时变换函数计算配准图像的全局配准度量。

假设将点对(p_a，q_a)添加到集合(P，Q)后得到变换函数h_a，计算形变后图像与参考图像之间的互信息MI(Sοh_a，R)作为全局配准度量。

在步骤S105，根据图像配准度量的变化量，以及新变换函数的对称KL距离，确定新增控制点对加入到当前控制点对集合中的概率值。

点对(p_a，q_a)添加到集合(P，Q)中的模拟退火概率P为

(9)

$P\_\mathrm{sKL}(p_a,q_a)=\left\{\begin{array}{cc}\exp \left(\frac{\mathrm{sKL}\left(h\right)-\mathrm{sKL}\left(h_a\right)}{T}\right)& \mathrm{sKL}\left(h\right)<\mathrm{sKL}\left(h_a\right)\\ 1& \mathrm{otherwise}\end{array}\right.$

其中MI是变换后的图像与参考图像之间的互信息，h是将点对(p_a，q_a )添加到集合(P，Q)之前的变换函数，Sοh是经过h变换后的图像，h_a是将点对(p_a，q_a)添加到集合(P，Q)之后的变换函数，Sοh_a是经过h_a变换后的图像。与h′(x)|相比，h′(x)|是没有确定对应关系时的临时变换函数，而h_a则是确定对应关系后的变换函数。

在步骤S106中，以概率接收该控制点对(p_a，q_a)加入到当前控制点对集合(P，Q)中，具体流程如图4。

在步骤S401中，利用随机数产生器产生一个[0，1]区间的随机数r。

在步骤S402中，若满足MI(Sοh_a，R)＞MI(Sοh，R)，则标志点对集合更新为P^new＝[P，p_a]，Q^new＝[Q，q_a]，否则，若r＜P(p_a，q_a)，将点对(p_a，q_a)添加到原始点对集合(P，Q)中。

转至步骤S101，重复以上步骤。

本领域普通技术人员可以理解实现上述各实施例提供的方法中的全部或部分步骤可以通过程序指令及相关的硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，该存储介质可以为ROM/RAM、磁盘、光盘等。

图5示出了本发明实施例提供的基于控制点无偏变换的图像配准***的逻辑结构原理，为了便于描述，仅示出了与本实施例相关的部分。此图像配准***可以为内置于图像处理设备中的软件单元、硬件单元或软硬件结合的单元。

参照图5，本发明实施例提供的基于控制点无偏变换的图像配准***包括当前变换函数构造单元51、新增源控制点选择单元52、新增目标控制点确定单元53、概率计算单元54以及当前控制点对集合更新单元55。

其中，当前变换函数构造单元51根据待配准数字图像中的当前控制点对集合构造当前变换函数，所述当前控制点对集合包含初始源控制点集合和初始目标控制点集合。新增源控制点选择单元52根据当前变换函数构造单元51构造得到的当前变换函数在当前源控制点周围选择新增源控制点，并将新增源控制点的映射位置作为新增目标控制点的初始位置。若新增源控制点选择单元52能找到新增源控制点，则由新增目标控制点确定单元53根据变换的无偏性和局部配准度量，调整新增目标控制点。然后，概率计算单元54将新增源控制点和新增目标控制点暂时加入到当前控制点对集合中，并计算新变换函数的对称KL距离，以及图像全局配准度量；概率计算单元54并根据图像全局配准度量的变化量，以及新变换函数的对称KL距离，确定新增控制点对加入到当前控制点对集合中的概率值；所述新增控制点对保包含相对应的新增源控制点和新增目标控制点。最后则由当前控制点对集合更新单元55根据概率计算单元54中计算得到概率值，将新增控制点对添加到当前控制点对集合中。上述各模块重复上述过程，直至图像配准精度达到要求为止。

进一步地，上述当前变换函数构造单元51根据如下公式构造当前变换函数h(x)：

$h\left(x\right)=x+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{i}U\left(\right||x-p_i|\left|\right)$

其中x是图像中的任意点，x∈R²，R为实数集合，p_i，i＝1，…，N是初始源控制点集合，U(r)是径向基函数，系数w_i，i＝1，…，N由方程KW＝Q-P求解得到，其中K是N×N矩阵，K_ij＝U(||p_i-p_j||)，W＝[w₁，w₂，…，w_N]^T，Q为初始目标控制点集合。

进一步地，上述新增源控制点选择单元52根据如下公式选择的新增源控制点：

$p_a=\underset{x}{\arg }\max |\mathrm{Dh}\left(x\right)-1|,$ 且distT_min＜||x-p_i||＜distT_max，x∈E， $\&ForAll;p_i\&Element;P$

其中，x∈R²，R为实数集合，h(x)为当前变换函数，P＝{p_i，i＝1，…，N}为初始源控制点集合，E是图像边缘点集合，Dh(x)为当前变换函数h(x)的Jacobian行列式值，distT_min和distT_max分别是最小距离阈值和最大距离阈值；

上述新增目标控制点确定单元53首先在以新增目标控制点的初始位置为中心的局部区域中，选取落在图像边缘上的点作为备选目标控制点，将新增源控制点和备选目标控制点临时添加到控制点集合，构造临时变换函数h′并计算此临时变换函数h′的对称KL距离sKL(h′)；

sKL(h′)＝∫(|Dh′(x)-1|)log|Dh′(x)|dx

其中，x∈R²，R为实数集合，h′(x)|为临时变换函数，Dh′(x)为临时变换函数h′(x)|的Jacobian行列式值；

上述新增目标控制点确定单元53再根据临时变换函数h′计算新增源控制点所处局部区域和新增目标控制点所处局部区域的配准度量MSD：

$\mathrm{MSD}(\mathrm{Loc}\left(q_a\right),\mathrm{Loc}\left(p_a\right))=\underset{(i^{\&prime;},j^{\&prime;})\&Element;\mathrm{Loc}\left(p_a\right)}{\underset{(i,j)\&Element;\mathrm{Loc}\left(q_a\right)}{\mathrm{\&Sigma;}}}|f(i,j)-f(i^{\&prime;},j^{\&prime;})|$

其中，p_a和Loc(p_a)分别为新增源控制点及其所处局部区域，q_a和Loc(q_a)分别为新增目标控制点及其所处局部区域，(i，j)表示参考图像中的点，(i，j′)表示待配准图像中的点，f(i，j)是参考图像在点(i，j)处的强度值，f(i′，j′)是待配准图像在点(i′，j′)处的强度值；

上述新增目标控制点确定单元53再根据临时变换函数的对称KL距离和所述配准度量MSD，根据如下公式确定最优的新增目标控制点以作为新增源控制点的对应点：

其中h′是将临时点对(p_a，x)添加到当前控制点对集合(P，Q)后得到的变换函数，Sοh′是经过h′变换后的图像，MSD是均方误差，sKL是对称KL距离。

进一步地，上述概率计算单元54根据如下公式确定新增控制点对加入到当前控制点对集合中的概率值：

$P\_\mathrm{sKL}(p_a,q_a)=\left\{\begin{array}{cc}\exp \left(\frac{\mathrm{sKL}\left(h\right)-\mathrm{sKL}\left(h_a\right)}{T}\right)& \mathrm{sKL}\left(h\right)<\mathrm{sKL}\left(h_a\right)\\ 1& \mathrm{otherwise}\end{array}\right.$

其中MI是变换后的图像与参考图像之间的互信息，h是将新增控制点对(p_a，q_a)添加到当前控制点对集合(P，Q)之前的变换函数，Sοh是经过当前变换函数h变换后的图像；h_a是将点对(p_a，q_a)添加到集合(P，Q)之后的变换函数，Sοh_a是经过h_a变换后的图像；

上述当前控制点对集合更新单元55首先利用随机数产生器产生一个[0，1]区间的随机数r；若满足MI(Sοh_a，R)＞MI(Sοh，R)，则将当前控制点对集合更新为P^new＝[P，p_a]，Q^new＝[Q，q_a]；否则，若r＜P(p_a，q_a)，则将点对(p_a，q_a)添加到原始点对集合(P，Q)中；

其中，P为当前源控制点集合，Q为当前目标控制点集合。

上述图像配准***进行图像配准的原理上文中描述的图像配准方法相同，不再一一赘述。

在本发明实施例中，在已知少量控制点对的前提下，根据当前控制点对提供的变换函数，为搜索到对应控制点确定一个有限的、比较精确的局部区域，以利于对应点的搜索。在搜索对应点的过程中，利用了变换无偏性和配准度量的综合度量，可以构造无偏性较好的变换函数。通过不断增加新的控制点对，逐步提高图像配准精度。该方法算法简单，只需要少量的初始控制点，是一种简便、有效的基于控制点对应关系的图像配准方法。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于控制点无偏变换的图像配准方法，其特征在于，所述方法包括下述步骤：

步骤A，根据待配准数字图像中的当前控制点对集合构造当前变换函数；所述当前控制点对集合包含初始源控制点集合和初始目标控制点集合；

步骤B，根据步骤A中构造得到的当前变换函数在当前源控制点周围选择新增源控制点，并将新增源控制点的映射位置作为新增目标控制点的初始位置；

步骤C，若步骤B中能找到新增源控制点，则根据变换的无偏性和局部配准度量，调整新增目标控制点；

步骤D，将新增源控制点和新增目标控制点暂时加入到当前控制点对集合中，并计算新变换函数的对称KL距离，以及图像全局配准度量；

步骤E，根据图像全局配准度量的变化量，以及新变换函数的对称KL距离，确定新增控制点对加入到当前控制点对集合中的概率值；所述新增控制点对包含相对应的新增源控制点和新增目标控制点；

步骤F，若步骤E中计算得到概率值大于随机产生的概率参考值，则将新增控制点对添加到当前控制点对集合中；

重复步骤A至步骤F。

2.如权利要求1所述的图像配准方法，其特征在于，所述步骤A根据如下公式构造当前变换函数h(x)：

$h\left(x\right)=x+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{i}U\left(\right||x-p_i|\left|\right)$

其中x是图像中的任意点，x∈R²，R为实数集合，p_i，i＝1，…，N是初始源控制点集合，U(r)是径向基函数，系数w_i，i＝1，…，N由方程KW＝Q-P求解得到，其中K是N×N矩阵，K_ij＝U(||p_i-p_j)，W＝[w₁，w₂，…，w_N]^T，Q为初始目标控制点集合。

3.如权利要求1所述的图像配准方法，其特征在于，所述步骤B中根据如下公式选择的新增源控制点：

$p_a=\underset{x}{\arg }\max |\mathrm{Dh}\left(x\right)-1|,$ 且distT_min＜||x-p_i||＜distT_max，x∈E， $\&ForAll;p_i\&Element;P$

其中，x∈R²，R为实数集合，h(x)为当前变换函数，P＝{p_i，i＝1，…，N}为初始源控制点集合，E是图像边缘点集合，Dh(x)为当前变换函数h(x)的Jacobian行列式值，distT_min和distT_max分别是最小距离阈值和最大距离阈值。

4.如权利要求1所述的图像配准方法，其特征在于，所述步骤C具体包括：

步骤C1，在以新增目标控制点的初始位置为中心的局部区域中，选取落在图像边缘上的点作为备选目标控制点，将新增源控制点和备选目标控制点临时添加到控制点集合，构造临时变换函数h′并计算此临时变换函数h′的对称KL距离sKL(h′)；

sKL(h′)＝∫(|Dh′(x)-1|)log|Dh′(x)|dx

其中，x∈R²，R为实数集合，h′(x)|为临时变换函数，Dh′(x)为临时变换函数h′(x)|的Jacobian行列式值；

步骤C2，根据临时变换函数h′计算新增源控制点所处局部区域和新增目标控制点所处局部区域的配准度量MSD：

$\mathrm{MSD}(\mathrm{Loc}\left(q_a\right),\mathrm{Loc}\left(p_a\right))=\underset{(i^{\&prime;},j^{\&prime;})\&Element;\mathrm{Loc}\left(p_a\right)}{\underset{(i,j)\&Element;\mathrm{Loc}\left(q_a\right)}{\mathrm{\&Sigma;}}}|f(i,j)-f(i^{\&prime;},j^{\&prime;})|$

其中，p_a和Loc(p_a)分别为新增源控制点及其所处局部区域，q_a和Loc(q_a)分别为新增目标控制点及其所处局部区域，(i，j)表示参考图像中的点，(i′，j′)表示待配准图像中的点，f(i，j)是参考图像在点(i，j)处的强度值，f(i′，j)是待配准图像在点(i′，j′)处的强度值；

步骤C3，根据临时变换函数的对称KL距离和所述配准度量MSD，根据如下公式确定最优的新增目标控制点以作为新增源控制点的对应点：

其中h′是将临时点对(p_a，x)添加到当前控制点对集合(P，Q)后得到的变换函数，Sοh′是经过h′变换后的图像，MSD是均方误差，sKL是对称KL距离。

5.如权利要求1所述的图像配准方法，其特征在于，所述步骤E根据如下公式确定新增控制点对加入到当前控制点对集合中的概率值：

$P\_\mathrm{sKL}(p_a,q_a)=\left\{\begin{array}{cc}\exp \left(\frac{\mathrm{sKL}\left(h\right)-\mathrm{sKL}\left(h_a\right)}{T}\right)& \mathrm{sKL}\left(h\right)<\mathrm{sKL}\left(h_a\right)\\ 1& \mathrm{otherwise}\end{array}\right.$

其中，MI是变换后的图像与参考图像之间的互信息，h是将新增控制点对(p_a，q_a)添加到当前控制点对集合(P，Q)之前的变换函数，Sοh是经过当前变换函数h变换后的图像；h_a是将点对(p_a，q_a)添加到集合(P，Q)之后的变换函数，Sοh_a是经过h_a变换后的图像。

6.如权利要求5所述的图像配准方法，其特征在于，所述步骤F包括：

步骤F1，利用随机数产生器产生一个[0，1]区间的随机数r；

步骤F2，若满足MI(Sοh_a，R)＞MI(Sοh，R)，则当前控制点对集合更新为P^new＝[P，p_a]，Q^new＝[Q，q_a]；否则，若r＜P(p_a，q_a)，则将点对(p_a，q_a)添加到当前控制点对集合(P，Q)中；

其中，P为当前源控制点集合，Q为当前目标控制点集合。

7.一种基于控制点无偏变换的图像配准***，其特征在于，包括：

当前变换函数构造单元，用于根据待配准数字图像中的当前控制点对集合构造当前变换函数；所述当前控制点对集合包含初始源控制点集合和初始目标控制点集合；

新增源控制点选择单元，用于根据所述当前变换函数构造单元构造得到的当前变换函数在当前源控制点周围选择新增源控制点，并将新增源控制点的映射位置作为新增目标控制点的初始位置；

新增目标控制点确定单元，若所述新增源控制点选择单元能找到新增源控制点，则用于根据变换的无偏性和局部配准度量，调整新增目标控制点；

概率计算单元，用于将新增源控制点和新增目标控制点暂时加入到当前控制点对集合中，并计算新变换函数的对称KL距离，以及图像全局配准度量；并根据图像全局配准度量的变化量，以及新变换函数的对称KL距离，确定新增控制点对加入到当前控制点对集合中的概率值；所述新增控制点对保包含相对应的新增源控制点和新增目标控制点；

当前控制点对集合更新单元，若所述概率计算单元中计算得到概率值大于随机产生的概率参考值，则将新增控制点对添加到当前控制点对集合中。

8.如权利要求7所述的基于控制点无偏变换的图像配准***，其特征在于，所述当前变换函数构造单元根据如下公式构造当前变换函数h(x)：

$h\left(x\right)=x+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{i}U\left(\right||x-p_i|\left|\right)$

其中x是图像中的任意点，x∈R²，R为实数集合，p_i，i＝1，…，N是初始源控制点集合，U(r)是径向基函数，系数w_i，i＝1，…，N由方程KW＝Q-P求解得到，其中K是N×N矩阵，K_ij＝U(||p_i-p_j)，W＝[w₁，w₂，…，w_N]^T，Q为初始目标控制点集合。

9.如权利要求7所述的基于控制点无偏变换的图像配准***，其特征在于，所述新增源控制点选择单元根据如下公式选择的新增源控制点：

$p_a=\underset{x}{\arg }\max |\mathrm{Dh}\left(x\right)-1|,$ 且distT_min＜||x-p_i||＜distT_max，x∈E， $\&ForAll;p_i\&Element;P$

其中，x∈R²，R为实数集合，h(x)为当前变换函数，P＝{p_i，i＝1，…，N}为初始源控制点集合，E是图像边缘点集合，Dh(x)为当前变换函数h(x)的Jacobian行列式值，distT_min和distT_max分别是最小距离阈值和最大距离阈值；

所述新增目标控制点确定单元，首先在以新增目标控制点的初始位置为中心的局部区域中，选取落在图像边缘上的点作为备选目标控制点，将新增源控制点和备选目标控制点临时添加到控制点集合，构造临时变换函数h′并计算此临时变换函数h′的对称KL距离sKL(h′)；

sKL(h′)＝∫(Dh′(x)-1|)log|Dh′(x)|dx

其中，x∈R²，R为实数集合，h′(x)|为临时变换函数，Dh′(x)为临时变换函数h′(x)|的Jacobian行列式值；

所述新增目标控制点确定单元再根据临时变换函数h′计算新增源控制点所处局部区域和新增目标控制点所处局部区域的配准度量MSD：

$\mathrm{MSD}(\mathrm{Loc}\left(q_a\right),\mathrm{Loc}\left(p_a\right))=\underset{(i^{\&prime;},j^{\&prime;})\&Element;\mathrm{Loc}\left(p_a\right)}{\underset{(i,j)\&Element;\mathrm{Loc}\left(q_a\right)}{\mathrm{\&Sigma;}}}|f(i,j)-f(i^{\&prime;},j^{\&prime;})|$

其中，p_a和Loc(p_a)分别为新增源控制点及其所处局部区域，q_a和Loc(q_a)分别为新增目标控制点及其所处局部区域，(i，j)表示参考图像中的点，(i′，j′)表示待配准图像中的点，f(i，j)是参考图像在点(i，j)处的强度值，f(i′，j)是待配准图像在点(i，j)处的强度值；

所述新增目标控制点确定单元再根据临时变换函数的对称KL距离和所述配准度量MSD，根据如下公式确定最优的新增目标控制点以作为新增源控制点的对应点：

其中h′是将临时点对(p_a，x)添加到当前控制点对集合(P，Q)后得到的变换函数，Sοh′是经过h′变换后的图像，MSD是均方误差，sKL是对称KL距离。

10.如权利要求7所述的基于控制点无偏变换的图像配准***，其特征在于：

所述概率计算单元根据如下公式确定新增控制点对加入到当前控制点对集合中的概率值：

$P\_\mathrm{sKL}(p_a,q_a)=\left\{\begin{array}{cc}\exp \left(\frac{\mathrm{sKL}\left(h\right)-\mathrm{sKL}\left(h_a\right)}{T}\right)& \mathrm{sKL}\left(h\right)<\mathrm{sKL}\left(h_a\right)\\ 1& \mathrm{otherwise}\end{array}\right.$

其中，MI是变换后的图像与参考图像之间的互信息，h是将新增控制点对(p_a，q_a)添加到当前控制点对集合(P，Q)之前的变换函数，Sοh是经过当前变换函数h变换后的图像；h_a是将点对(p_a，q_a)添加到集合(P，Q)之后的变换函数，Sοh_a是经过h_a变换后的图像；

所述当前控制点对集合更新单元首先利用随机数产生器产生一个[0，1]区间的随机数r；若满足MI (Sοh_a，R)＞MI(Sοh，R)，则将当前控制点对集合更新为P^new＝[P，p_a]，Q^new＝[Q，q_a]；否则，若r＜P(p_a，q_a)，则将点对(p_a，q_a)添加到原始点对集合(P，Q)中；

其中，P为当前源控制点集合，Q为当前目标控制点集合。

Images