CN105427308A - 一种稀疏和稠密特征匹配结合的图像配准方法 - Google Patents

Abstract

本发明提供了一种稀疏和稠密特征匹配结合的图像配准方法，将稀疏和稠密特征匹配相结合，得到一种新的数学模型，新模型包含两个变量：非刚性几何变换和离散位移流场，前者适用于稀疏匹配流，我们引入了局部线性约束以调整变换，使得该问题是适定的，后者适用于稠密匹配流，采取与SIFT流相似的模型，同时采用置信传播算法优化求解，能够针对包含未知非刚性运动的遥感图像，获得精确的像素间对比的匹配。

Description

一种稀疏和稠密特征匹配结合的图像配准方法

技术领域

本发明涉及一种稀疏和稠密特征匹配结合的图像配准方法，属于图像处理技术领域。

背景技术

图像配准是遥感领域中的一个根本性和挑战性问题，也是很多广泛应用(包含地形重建、环境监测、变化探测、图像镶嵌、图像融合及地图更新等)中的一个先决条件。

图像配准旨在建立不同时间不同视角或通过不同探测器获取同一场景的两幅图像之间像素的对应关系。配准问题可依据具体应用和数据形式分为刚性或非刚性。刚性配准(仅涉及少量参数)相对容易，且已得到了广泛研究。相比而言，非刚性配准较为困难，由于潜在非刚性变换模型通常无法事先获知且较为复杂，难以建模。尽管如此，非刚性配准对遥感图像而言仍非常重要，因为遥感图像通常存在一些因地貌变化或成像视角变换导致的局部非刚性形变，而这些形变无法基于简单刚性模型实现“精确配准”。

一种早期广泛使用的非刚性图像配准算法为光流法。其直接通过最小化像素间灰度的差异来计算一个全局稠密的位移场，通常适用于对两幅非常相似的图像进行配准，例如视频序列中的相邻帧图像。光流法的典型假设包含光照恒定，以及像素位移场应具有分段光滑性。然而，由于光照变化、视角变换和噪声干扰，像素灰度值往往变得不太可靠。近年来，Liu等人提出了可容忍高度类内变化的SIFT流配准算法。与光流法匹配像素灰度值不同的是，SIFT流算法匹配图像中稠密采样像素的SIFT描述子。该算法在复杂场景情况下展现出令人满意的配准效果，但对大尺度缩放和旋转的鲁棒性依然较差。

然而，各种局部鲁棒特征的涌现与发展为图像配准带来了新的解决途径。这类方法通过匹配从图像中提取出的局部特征信息，估计稀疏特征位置之间的空间变换函数并用其来拟合图像对之间的真实几何变换关系，从而实现对原始图像的配准。这类方法中采用的特征可由不同层次的简单几何实体来描述，包括点、线段、轮廓和区域等。通常，较高层次的特征很难被准确地提取；点特征是最简单的特征形式，同时也是最常见的特征形式，由于高层次特征例如线和轮廓可被描述为点的集合。从这层意义上来说，配准问题可简化为找出两个提取的特征点集之间的对应关系及潜在的空间变换。这种基于特征的方法对典型表观变化和场景位移具有较强的鲁棒性，且在正确实施的情况下，拥有更快的速度。他们在刚性场景(例如：在图像拼接领域)和轻度非刚性场景(例如：在医学成像领域)匹配中取得了很大成就，但在处理高度非刚性(例如：在地形重建领域)情形时无法取得满意效果。其本质在于这里配准是通过对稀疏特征匹配进行插值得到，而不是通过计算所有像素之间的对应关系得到，当真实匹配为非刚性且变换模型未知时就会带来问题。

发明内容

为了解决现有技术的不足，本发明提供了一种稀疏和稠密特征匹配结合的图像配准方法，能够针对包含未知非刚性运动的遥感图像，获得精确的像素间对比的匹配。

本发明为解决其技术问题所采用的技术方案是：提供了一种稀疏和稠密特征匹配结合的图像配准方法，包括以下步骤：

(1)建立基于局部线性约束的稀疏特征匹配能量函数：

(1-1)采用特征探测法推断得到待配准的两个图像的假定匹配集x_n与y_n分别表示待配准的两幅图像中特征点空间位置的二位列向量；所述假定匹配集S包含错误匹配和正确匹配，其中正确匹配根据带匹配的两个图像之间的几何变换Z确定，即如果(x_n,y_n)是一个正确的匹配，则y_n＝Z(x_n)是一个正确的匹配；

(1-2)初始位置加位移函数v将定义变换Z：Ζ(x)＝x+v(x)，其中v在函数空间H中模拟，所述函数空间H为向量值的再生核希尔伯特空间，通过矩阵值的核Γ:R²×R²→R^2×2与对角线高斯核定义得到，其中R为实数空间，I为单位矩阵，x_i和x_j为再生核希尔伯特空间实数空间中两个特征点的空间位置，β为高斯窗参数，参数β为设置值，其范围为0.01～1，则变换Z通过以下形式表示：

$Z\left(x\right)=x+v\left(x\right)=x+\underset{n=1}{\overset{N}{\Σ}}\mathrm{\Γ}(x,x_n)c_n...\left(1\right)$

其中，x表示任意一个图像点的空间位置，c_n表示一个2×1维待求解的向量系数；

(1-3)利用对角元素{p_n}组成的对角矩阵P指出匹配可信度，对角矩阵P中的任意元素p_n∈[0,1]为0到1之间的实数，当p_n＝1时表示(x_n,y_n)为一个正确的匹配，当p_n＝0时表示(x_n,y_n)为一个错误的匹配；

(1-4)创建大小为N×N的权值矩阵W，且使当x_j不属于x_i的Q邻域时，W中的元素W_ij＝0，x_i的Q邻域指欧氏距离最近的Q个元素；在约束下通过代价函数最小化基于局部线性约束的重构误差，并采用最小二乘求解出W，得到以下能量函数：

${\mathrm{\ϵ}}_1(P,Z)=\underset{n=1}{\overset{N}{\Σ}}{p}_n\left|\right|y_n-Z\left(x_n\right)|{|}^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n+\mathrm{\λ}\underset{i=1}{\overset{N}{\Σ}}{p}_i||Z\left(x_i\right)-\underset{j=1}{\overset{N}{\Σ}}{W}_{ij}Z\left(x_j\right)|{|}^2...\left(2\right)$

其中，参数T₁表示确定性退火的温度，该参数用于结合退火速率参数r用于逐步求解能量函数，η表示对错误匹配的惩罚，λ表示正则化参数基于局部线性约束的正则化参数，参数T₁、η和λ均为设置值，T₁初始值的范围为0.05～5，η的范围为0.01～1，λ的范围为100～10000；

(2)建立基于SIFT流的能量函数：

用m表示图像像素的网格坐标，u(m)表示位移向量，第一幅图像中的点m对应第二幅图像中的点m+u(m)，s₁(x)和s₂(x)分别为两幅图像逐像素采样的SIFT特征，集合e包含四邻域系统中所有空间邻域，则基于稠密像素SIFT配准的能量函数为：

${\mathrm{\ϵ}}_2\left(u\right)=\underset{m}{\Σ}min\left(\right||s_1\left(m\right)-s_2(m+u(m))|{|}_1,t)+\underset{m}{\Σ}\mathrm{\γ}\left|\right|u\left(m\right)|{|}_1+\underset{(m,q)\∈e}{\Σ}\underset{i=1}{\overset{2}{\Σ}}min(\mathrm{\α}|u^i\left(m\right)-u^i\left(q\right)|,d)...\left(3\right)$

其中，t和d为两个截断L1范数的阈值，γ表示保持位移场平滑性的权重，α表示保持位移场连续性的权重，m和q表示图像像素的网格坐标，uⁱ(m)或uⁱ(q)表示对应位移向量的第i个分量，这里t设为s₁(x)与s₂(x)差值的中位数；d、γ和α均为设置值，d的范围为4～400，γ的范围为0.05～0.5，α的范围为0.2～20；

(3)将基于局部线性约束的稀疏特征匹配能量函数和基于SIFT的能量函数整合，得到稀疏和稠密特征匹配结合的图像配准能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_3\left(P,v,u\right)=\underset{n=1}{\overset{N}{\Σ}}{p}_n\left|\right|y_n-x_n+v\left(x_n\right)|{|}^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n\\ +\mathrm{\λ}\underset{i=1}{\overset{N}{\Σ}}{p}_i\left|\right|v\left(x_i\right)-\underset{j=1}{\overset{N}{\Σ}}{W}_{ij}v\left(x_j\right)|{|}^2\\ +\frac{\mathrm{\δ}}{L}\underset{m}{\Σ}\min \left(\left|\right|s_1\left(m\right)-s_2\left(m+u\left(m\right)\right)|{|}_1,t\right)+\frac{\mathrm{\δ}}{L}\underset{m}{\Σ}\mathrm{\γ}\left|\right|v\left(m\right)-u\left(m\right)|{|}_1\\ +\frac{\mathrm{\δ}}{L}\underset{\left(m,q\right)\∈e}{\Σ}\underset{i=1}{\overset{2}{\Σ}}\min \left(\mathrm{\α}\left|u^i\left(m\right)-u^i\left(q\right)\right|,d\right)\end{array}\mathrm{...}\left(4\right)$

其中，L表示评价像素的总数，δ表示设置的控制稀疏和稠密匹配平衡的正数参数，ε表示所有四邻域的集合；

(4)解答能量函数，执行图像配准：

(4-1)首先提取稀疏和稠密特征匹配结合的图像配准能量函数ε₃(P,v,u)中与v相关的项，得到以下能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_4\left(P,C\right)=\left|\right|P^{1/2}\left(Y-X-KC\right)|{|}_Z^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n\\ +\mathrm{\λ}\left|\right|P^{1/2}(I-W)KC|{|}_Z^2+\frac{\mathrm{\δ}\mathrm{\γ}}{L}||VC-U|{|}_Z^2\end{array}\mathrm{...}\left(5\right)$

其中，X＝(x₁,...,x_N)^T，Y＝(y₁,...,y_N)^T，K∈R^N×N且C表示格拉姆矩阵，为位移函数v的系数矩阵，V∈R^L×N且U＝(u₁,...,u_L)^T表示尺寸为L×2的流场，为弗罗贝尼斯范数；

(4-2)通过确定性退火步骤解答位移函数v：

(4-2-1)初始化参数T₁、r和λ，设置迭代次数；

(4-2-2)初始化P＝I，C＝0；

(4-2-3)通过计算式(5)的极值完成用当前位移函数v更新匹配可信度，得到匹配可信度的闭合形式为：

$p_n=e^{-\frac{\left|\right|y_n-x_n-v\left(x_n\right)|{|}^2-\mathrm{\η}+T_1}{T}}...\left(6\right)$

通过公式(6)更新p_n；

(4-2-4)通过以下线性方程式更新C：

$\[KPK+\mathrm{\λ}K{(I-W)}^{T}P(I-W)K+\frac{\mathrm{\δ}\mathrm{\γ}}{L}{V}^{T}V\]C=KPY-KPX+\frac{\mathrm{\δ}\mathrm{\γ}}{L}{V}^{T}U...\left(7\right)$

(4-2-5)重复步骤(4-2-3)和(4-2-4)直到...，结束交替更新；

(4-2-6)减少T₁和λ，若未达到预设的迭代次数或能量函数ε₄(P,C)未收敛，则返回步骤(4-2-3)；否则进入步骤(4-2-7)；

(4-2-7)结束确定性退火步骤；

(4-4)提取稀疏和稠密特征匹配结合的图像配准能量函数ε₃(P,v,u)中与u相关的项，得到以下能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_5\left(u\right)=\underset{m}{\Σ}\min \left(\left|\right|s_1\left(m\right)-s_2\left(m+u\left(m\right)\right)|{|}_1,t\right)\\ +\underset{\left(m,q\right)\∈\mathrm{\ϵ}}{\Σ}\underset{i=1}{\overset{2}{\Σ}}\min \left(\mathrm{\α}\left|u^i\left(m\right)-u^i\left(q\right)\right|,d\right)+\underset{m}{\Σ}\mathrm{\γ}\left|\right|v\left(m\right)-u\left(m\right)|{|}^2\end{array}\mathrm{...}\left(8\right)$

(4-5)利用SIFT流算法修改公式(8)中从Σ_mγ||u(m)||₁到Σ_mγ||v(m)-u(m)||²的小位移项以解答u；

(4-6)通过u执行图像配准。

步骤(1-2)所述的高斯窗参数β设置为0.1。

步骤(1-4)中，T₁＝0.5，η＝0.1，λ＝1000。

步骤(2)中，d＝40，γ＝0.005，α＝2。

本发明基于其技术方案所具有的有益效果在于：

(1)本发明的稀疏和稠密特征匹配结合的图像配准方法，针对遥感图像存在地形起伏导致非刚性形变的问题，对点匹配进行局部线性约束，能够在图像变换后保护特征集中的局部结构，从而提高图像配准精度；

(2)本发明的稀疏和稠密特征匹配结合的图像配准方法，将稀疏和稠密特征匹配相结合，得到一种新的数学模型，新模型包含两个变量：非刚性几何变换和离散位移流场，前者适用于稀疏匹配流，我们引入了局部线性约束以调整变换，使得该问题是适定的，后者适用于稠密匹配流，采取与SIFT流相似的模型，同时采用置信传播算法优化求解；

(3)本发明的稀疏和稠密特征匹配结合的图像配准方法，在理想情况下，这两个变量是一致的，但同时求解这两个变量非常困难，我们采用迭代的策略固定一个变量求解另一变量，使得求解变得相对容易，并且同时可以避免陷入不理想的局部最优解；

(4)本发明的稀疏和稠密特征匹配结合的图像配准方法，同时具有基于稀疏特征匹配的高效的优点，以及基于稠密特征匹配的精确的优点。

具体实施方式

下面结合实施例对本发明作进一步说明。

本发明提供了一种稀疏和稠密特征匹配结合的图像配准方法，包括以下步骤：

(1)建立基于局部线性约束的稀疏特征匹配能量函数：

(1-1)采用特征探测法推断得到待配准的两个图像的假定匹配集x_n与y_n分别表示待配准的两幅图像中特征点空间位置的二位列向量；所述假定匹配集S包含错误匹配和正确匹配，其中正确匹配根据带匹配的两个图像之间的几何变换Z确定，即如果(x_n,y_n)是一个正确的匹配，则y_n＝Z(x_n)是一个正确的匹配；

(1-2)初始位置加位移函数v将定义变换Z：Ζ(x)＝x+v(x)，其中v在函数空间H中模拟，所述函数空间H为向量值的再生核希尔伯特空间，通过矩阵值的核Γ:R²×R²→R^2×2与对角线高斯核定义得到，其中R为实数空间，I为单位矩阵，x_i和x_j为再生核希尔伯特空间实数空间中两个特征点的空间位置，β为高斯窗参数，参数β为设置值，其范围为0.01～1，则变换Z通过以下形式表示：

$Z\left(x\right)=x+v\left(x\right)=x+\underset{n=1}{\overset{N}{\Σ}}\mathrm{\Γ}(x,x_n)c_n...\left(1\right)$

其中，x表示任意一个图像点的空间位置，c_n表示一个2×1维待求解的向量系数；

(1-3)利用对角元素{p_n}组成的对角矩阵P指出匹配可信度，对角矩阵P中的任意元素p_n∈[0,1]为0到1之间的实数，当p_n＝1时表示(x_n,y_n)为一个正确的匹配，当p_n＝0时表示(x_n,y_n)为一个错误的匹配；

(1-4)创建大小为N×N的权值矩阵W，且使当x_j不属于x_i的Q邻域时，W中的元素W_ij＝0，x_i的Q邻域指欧氏距离最近的Q个元素；在约束下通过代价函数最小化基于局部线性约束的重构误差，并采用最小二乘求解出W，得到以下能量函数：

${\mathrm{\ϵ}}_1(P,Z)=\underset{n=1}{\overset{N}{\Σ}}{p}_n\left|\right|y_n-Z\left(x_n\right)|{|}^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n+\mathrm{\λ}\underset{i=1}{\overset{N}{\Σ}}{p}_i||Z\left(x_i\right)-\underset{j=1}{\overset{N}{\Σ}}{W}_{ij}Z\left(x_j\right)|{|}^2...\left(2\right)$

其中，参数T₁表示确定性退火的温度，该参数用于结合退火速率参数r用于逐步求解能量函数，即，T₁ ^new＝T₁ ^old·r，η表示对错误匹配的惩罚，λ表示正则化参数基于局部线性约束的正则化参数，参数T₁、η和λ均为设置值，T₁初始值的范围为0.05～5，η的范围为0.01～1，λ的范围为100～10000；

(2)建立基于SIFT流的能量函数：

用m表示图像像素的网格坐标，u(m)表示位移向量，第一幅图像中的点m对应第二幅图像中的点m+u(m)，s₁(x)和s₂(x)分别为两幅图像逐像素采样的SIFT特征，集合e包含四邻域系统中所有空间邻域，则基于稠密像素SIFT配准的能量函数为：

${\mathrm{\ϵ}}_2\left(u\right)=\underset{m}{\Σ}min\left(\right||s_1\left(m\right)-s_2(m+u(m))|{|}_1,t)+\underset{m}{\Σ}\mathrm{\γ}\left|\right|u\left(m\right)|{|}_1+\underset{(m,q)\∈e}{\Σ}\underset{i=1}{\overset{2}{\Σ}}min(\mathrm{\α}|u^i\left(m\right)-u^i\left(q\right)|,d)...\left(3\right)$

其中，t和d为两个截断L1范数的阈值，即设定阈值来限定L1范数的最大值，γ表示保持位移场平滑性的权重，α表示保持位移场连续性的权重，m和q表示图像像素的网格坐标，uⁱ(m)或uⁱ(q)表示对应位移向量的第i个分量，这里t设为s₁(x)与s₂(x)差值的中位数；d、γ和α均为设置值，d的范围为4～400，γ的范围为0.05～0.5，α的范围为0.2～20；

(3)将基于局部线性约束的稀疏特征匹配能量函数和基于SIFT的能量函数整合，得到稀疏和稠密特征匹配结合的图像配准能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_3\left(P,v,u\right)=\underset{n=1}{\overset{N}{\Σ}}{p}_n\left|\right|y_n-x_n+v\left(x_n\right)|{|}^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n\\ +\mathrm{\λ}\underset{i=1}{\overset{N}{\Σ}}{p}_i\left|\right|v\left(x_i\right)-\underset{j=1}{\overset{N}{\Σ}}{W}_{ij}v\left(x_j\right)|{|}^2\\ +\frac{\mathrm{\δ}}{L}\underset{m}{\Σ}\min \left(\left|\right|s_1\left(m\right)-s_2\left(m+u\left(m\right)\right)|{|}_1,t\right)+\frac{\mathrm{\δ}}{L}\underset{m}{\Σ}\mathrm{\γ}\left|\right|v\left(m\right)-u\left(m\right)|{|}_1\\ +\frac{\mathrm{\δ}}{L}\underset{\left(m,q\right)\∈e}{\Σ}\underset{i=1}{\overset{2}{\Σ}}\min \left(\mathrm{\α}\left|u^i\left(m\right)-u^i\left(q\right)\right|,d\right)\end{array}\mathrm{...}\left(4\right)$

其中，L表示评价像素的总数，δ表示设置的控制稀疏和稠密匹配平衡的正数参数，ε表示所有四邻域的集合；

(4)解答能量函数，执行图像配准：

(4-1)首先提取稀疏和稠密特征匹配结合的图像配准能量函数ε₃(P,v,u)中与v相关的项，得到以下能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_4\left(P,C\right)=\left|\right|P^{1/2}\left(Y-X-KC\right)|{|}_Z^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n\\ +\mathrm{\λ}\left|\right|P^{1/2}(I-W)KC|{|}_Z^2+\frac{\mathrm{\δ}\mathrm{\γ}}{L}||VC-U|{|}_Z^2\end{array}\mathrm{...}\left(5\right)$

其中，X＝(x₁,...,x_N)^T，Y＝(y₁,...,y_N)^T，K∈R^N×N且C表示格拉姆矩阵，为位移函数v的系数矩阵，V∈R^L×N且U＝(u₁,...,u_L)^T表示尺寸为L×2的流场，为弗罗贝尼斯范数；

(4-2)通过确定性退火步骤解答位移函数v：

(4-2-1)初始化参数T₁、r和λ，设置迭代次数；

(4-2-2)初始化P＝I，C＝0；

(4-2-3)通过计算式(5)的极值完成用当前位移函数v更新匹配可信度，得到匹配可信度的闭合形式为：

$p_n=e^{-\frac{\left|\right|y_n-x_n-v\left(x_n\right)|{|}^2-\mathrm{\η}+T_1}{T}}...\left(6\right)$

通过公式(6)更新p_n；

(4-2-4)通过以下线性方程式更新C：

$\[KPK+\mathrm{\λ}K{(I-W)}^{T}P(I-W)K+\frac{\mathrm{\δ}\mathrm{\γ}}{L}{V}^{T}V\]C=KPY-KPX+\frac{\mathrm{\δ}\mathrm{\γ}}{L}{V}^{T}U...\left(7\right)$

(4-2-5)重复步骤(4-2-3)和(4-2-4)直到...，结束交替更新；

(4-2-6)减少T₁和λ，若未达到预设的迭代次数或能量函数ε₄(P,C)未收敛，则返回步骤(4-2-3)，迭代次数可设置为500次，未收敛指上次与本次迭代能量函数ε₄(P,C)的变化率大于一定范围，如10^-5；否则进入步骤(4-2-7)；

(4-2-7)结束确定性退火步骤；

(4-4)提取稀疏和稠密特征匹配结合的图像配准能量函数ε₃(P,v,u)中与u相关的项，得到以下能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_5\left(u\right)=\underset{m}{\Σ}\min \left(\left|\right|s_1\left(m\right)-s_2\left(m+u\left(m\right)\right)|{|}_1,t\right)\\ +\underset{\left(m,q\right)\∈\mathrm{\ϵ}}{\Σ}\underset{i=1}{\overset{2}{\Σ}}\min \left(\mathrm{\α}\left|u^i\left(m\right)-u^i\left(q\right)\right|,d\right)+\underset{m}{\Σ}\mathrm{\γ}\left|\right|v\left(m\right)-u\left(m\right)|{|}^2\end{array}\mathrm{...}\left(8\right)$

(4-5)利用SIFT流算法修改公式(8)中从Σ_mγ||u(m)||₁到Σ_mγ||v(m)-u(m)||²的小位移项以解答u；

(4-6)通过u执行图像配准。

本实施例中，步骤(1-2)所述的高斯窗参数β设置为0.1。

步骤(1-4)中，T₁＝0.5，η＝0.1，λ＝1000。

步骤(2)中，d＝40，γ＝0.005，α＝2。

目标函数(4)的最后一项是图像点阵v和u的平均差，可通过降低图像点阵的取样率估算系数矩阵C以实现一个显著加速但却不影响其性能。在评估中统一使用取样间隔为10像素的取样策略。另外，通过步骤(4-2)可初始化独立于U的C。也可采用数据规范化实现和两组稀疏特征点的零平均值和单位方差。

参数设定在公式中主要有4个匹配稀疏特征的参数：K，r，λ和β。参数K控制线性重建最邻近的数字，参数T和r是确定性退火的初始变换和退火速率。参数A控制局部几何约束对变换T的影响。P确定非刚性变换特征点之间相互作用的范围。本实施例设定K＝15，T＝0.5，r＝0.93，λ＝1000，β＝0.1。稠密匹配公式中的参数根据SIFT流算法设定。方程式(8)在第一次迭代后需要根据方程式δγ＝103设置参数δ，其控制稀疏匹配和稠密匹配之间的平衡。

Claims (4)

1.一种稀疏和稠密特征匹配结合的图像配准方法，其特征在于包括以下步骤：

(1)建立基于局部线性约束的稀疏特征匹配能量函数：

(1-1)采用特征探测法推断得到待配准的两个图像的假定匹配集x_n与y_n分别表示待配准的两幅图像中特征点空间位置的二位列向量；所述假定匹配集S包含错误匹配和正确匹配，其中正确匹配根据带匹配的两个图像之间的几何变换Z确定，即如果(x_n,y_n)是一个正确的匹配，则y_n＝Z(x_n)是一个正确的匹配；

(1-2)初始位置加位移函数v将定义变换Z：Ζ(x)＝x+v(x)，其中v在函数空间H中模拟，所述函数空间H为向量值的再生核希尔伯特空间，通过矩阵值的核Γ:R²×R²→R^2×2与对角线高斯核定义得到，其中R为实数空间，I为单位矩阵，x_i和x_j为再生核希尔伯特空间实数空间中两个特征点的空间位置，β为高斯窗参数，参数β为设置值，其范围为0.01～1，则变换Z通过以下形式表示：

$Z\left(x\right)=x+v\left(x\right)=x+\underset{n=1}{\overset{N}{\Σ}}\mathrm{\Γ}(x,x_n)c_n...\left(1\right)$

其中，x表示任意一个图像点的空间位置，c_n表示一个2×1维待求解的向量系数；

(1-3)利用对角元素{p_n}组成的对角矩阵P指出匹配可信度，对角矩阵P中的任意元素p_n∈[0,1]为0到1之间的实数，当p_n＝1时表示(x_n,y_n)为一个正确的匹配，当p_n＝0时表示(x_n,y_n)为一个错误的匹配；

(1-4)创建大小为N×N的权值矩阵W，且使当x_j不属于x_i的Q邻域时，W中的元素W_ij＝0，x_i的Q邻域指欧氏距离最近的Q个元素；在约束 下通过代价函数最小化基于局部线性约束的重构误差，并采用最小二乘求解出W，得到以下能量函数：

${\mathrm{\ϵ}}_1(P,Z)=\underset{n=1}{\overset{N}{\Σ}}{p}_n\left|\right|y_n-Z\left(x_n\right)|{|}^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n+\mathrm{\λ}\underset{i=1}{\overset{N}{\Σ}}{p}_i||Z\left(x_i\right)-\underset{j=1}{\overset{N}{\Σ}}{W}_{ij}Z\left(x_j\right)|{|}^2...\left(2\right)$

其中，参数T₁表示确定性退火的温度，该参数用于结合退火速率参数r用于逐步求解能量函数，η表示对错误匹配的惩罚，λ表示正则化参数基于局部线性约束的正则化参数，参数T₁、η和λ均为设置值，T₁初始值的范围为0.05～5，η的范围为0.01～1，λ的范围为100～10000；

(2)建立基于SIFT流的能量函数：

用m表示图像像素的网格坐标，u(m)表示位移向量，第一幅图像中的点m对应第二幅图像中的点m+u(m)，s₁(x)和s₂(x)分别为两幅图像逐像素采样的SIFT特征，集合e包含四邻域系统中所有空间邻域，则基于稠密像素SIFT配准的能量函数为：

${\mathrm{\ϵ}}_2\left(u\right)=\underset{m}{\Σ}min\left(\right|\left|s_1\right(m)-s_2(m+u\left(m\right)\left)\right|{|}_1,t)+\underset{m}{\Σ}\mathrm{\γ}\left|\right|u\left(m\right)|{|}_1+\underset{(m,q)\∈e}{\Σ}\underset{i=1}{\overset{2}{\Σ}}min\left(\mathrm{\α}\right|u^i\left(m\right)-u^i\left(q\right)|,d)...\left(3\right)$

其中，t和d为两个截断L1范数的阈值，γ表示保持位移场平滑性的权重，α表示保持位移场连续性的权重，m和q表示图像像素的网格坐标，uⁱ(m)或uⁱ(q)表示对应位移向量的第i个分量，这里t设为s₁(x)与s₂(x)差值的中位数；d、γ和α均为设置值，d的范围为4～400，γ的范围为0.05～0.5，α的范围为0.2～20；

(3)将基于局部线性约束的稀疏特征匹配能量函数和基于SIFT的能量函数整合，得到稀疏和稠密特征匹配结合的图像配准能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_3(P,v,u)=\underset{n=1}{\overset{N}{\Σ}}{p}_n\left|\right|y_n-x_n+v\left(x_n\right)|{|}^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n\\ +\mathrm{\λ}\underset{i=1}{\overset{N}{\Σ}}{p}_i\left|\right|v\left(x_i\right)-\underset{j=1}{\overset{N}{\Σ}}{W}_{ij}v\left(x_j\right)|{|}^2\\ +\frac{\mathrm{\δ}}{L}\underset{m}{\Σ}min\left(\right|\left|s_1\right(m)-s_2(m+u\left(m\right)\left)\right|{|}_1,t)+\frac{\mathrm{\δ}}{L}\underset{m}{\Σ}\mathrm{\γ}\left|\right|v\left(m\right)-u\left(m\right)|{|}_1\\ +\frac{\mathrm{\δ}}{L}\underset{(m,q)\∈e}{\Σ}\underset{i=1}{\overset{2}{\Σ}}min\left(\mathrm{\α}\right|u^i\left(m\right)-u^i\left(q\right)|,d)\end{array}...\left(4\right)$

其中，L表示评价像素的总数，δ表示设置的控制稀疏和稠密匹配平衡的正数参数，ε表示所有四邻域的集合；

(4)解答能量函数，执行图像配准：

(4-1)首先提取稀疏和稠密特征匹配结合的图像配准能量函数ε₃(P,v,u)中与v相关的项，得到以下能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_4(P,C)=\left|\right|P^{1/2}(Y-X-KC)|{|}_Z^2-\mathrm{\η}\underset{n=1}{\overset{N}{\Σ}}{p}_n+T_1\underset{n=1}{\overset{N}{\Σ}}{p}_n{\mathrm{logp}}_n\\ +\mathrm{\λ}\left|\right|P^{1/2}(I-W)KC|{|}_Z^2+\frac{\mathrm{\δ}\mathrm{\γ}}{L}||VC-U|{|}_Z^2\end{array}...\left(5\right)$

其中，X＝(x₁,...,x_N)^T，Y＝(y₁,...,y_N)^T，K∈R^N×N且C表示格拉姆矩阵，为位移函数v的系数矩阵，V∈R^L×N且U＝(u₁,...,u_L)^T表示尺寸为L×2的流场，为弗罗贝尼斯范数；

(4-2)通过确定性退火步骤解答位移函数v：

(4-2-1)初始化参数T₁、r和λ，设置迭代次数；

(4-2-2)初始化P＝I，C＝0；

(4-2-3)通过计算式(5)的极值完成用当前位移函数v更新匹配可信度，得到匹配可信度的闭合形式为：

$p_n=e^{-\frac{\left|\right|y_n-x_n-v\left(x_n\right)|{|}^2-\mathrm{\η}+T_1}{T}}...\left(6\right)$

通过公式(6)更新p_n；

(4-2-4)通过以下线性方程式更新C：

$\[KPK+\mathrm{\λ}K{(I-W)}^{T}P(I-W)K+\frac{\mathrm{\δ}\mathrm{\γ}}{L}{V}^{T}V\]C=KPY-KPX+\frac{\mathrm{\δ}\mathrm{\γ}}{L}{V}^{T}U...\left(7\right)$

(4-2-5)重复步骤(4-2-3)和(4-2-4)直到...，结束交替更新；

(4-2-6)减少T₁和λ，若未达到预设的迭代次数或能量函数ε₄(P,C)未收敛，则返回步骤(4-2-3)；否则进入步骤(4-2-7)；

(4-2-7)结束确定性退火步骤；

(4-4)提取稀疏和稠密特征匹配结合的图像配准能量函数ε₃(P,v,u)中与u相关的项，得到以下能量函数：

$\begin{array}{c}{\mathrm{\ϵ}}_5\left(u\right)=\underset{m}{\Σ}min\left(\right|\left|s_1\right(m)-s_2(m+u\left(m\right)\left)\right|{|}_1,t)\\ +\underset{(m,q)\∈\mathrm{\ϵ}}{\Σ}\underset{i=1}{\overset{2}{\Σ}}min\left(\mathrm{\α}\right|u^i\left(m\right)-u^i\left(q\right)|,d)+\underset{m}{\Σ}\mathrm{\γ}\left|\right|v\left(m\right)-u\left(m\right)|{|}^2\end{array}...\left(8\right)$

(4-5)利用SIFT流算法修改公式(8)中从Σ_mγ||u(m)||₁到Σ_mγ||v(m)-u(m)||²的小位移项以解答u；

(4-6)通过u执行图像配准。

2.根据权利要求1所述的稀疏和稠密特征匹配结合的图像配准方法，其特征在于：步骤(1-2)所述的高斯窗参数β设置为0.1。

3.根据权利要求1所述的稀疏和稠密特征匹配结合的图像配准方法，其特征在于：步骤(1-4)中，T₁＝0.5，r＝0.93，η＝0.1，λ＝1000。

4.根据权利要求1所述的稀疏和稠密特征匹配结合的图像配准方法，其特征在于：步骤(2)中，d＝40，γ＝0.005，α＝2。