非刚性点集配准(non-rigid point set registration)是将某一个点集(源点集)与其发生非刚性变形后的点集(目标点集)进行匹配的过程.该技术在立体匹配(stereo matching)^{[1, 2]}、基于内容图像检索(content-based image retrieval)^[3]、图像配准(image registration)^[4]、形状识别(shape recognition)、超分辨率问题(super-resolution)^[5]等计算机视觉任务中是一项重要的研究问题, 它同时还被广泛应用在医学图像处理^[6]、模式识别^[7]、计算机图形学、地理信息系统等领域.在已有的非刚性点集配准方法中, 大体可以分为两个研究方向:(1)基于迭代或非迭代的算法; (2)基于学习或非学习的算法.由于本文算法主要涉及基于迭代的非刚性点集配准问题, 所以我们将基于第1类研究方向对当前的非刚性点集配准算法进行介绍和讨论.

基于非迭代的非刚性点集配准算法通常使用某种高级结构特征(high level structural features)来对两个点集之间的对应关系进行仅有一次的相似性评估, 然后直接确定两个点集之间的对应关系.在基于非迭代的配准模型中, 直线(line)^[8]、曲线(curve)^[9]、表面结构(surface)^[10]、Shape context^{[11, 12]}和图论(graphs)^{[13, 14]}等特征被用于两个点集之间相似性的评估.在非迭代算法中, Shape context和Graphs是最受欢迎的两种特征描述法, 前者的核心思想是通过最小化两个点集之间的分布差异来确定它们之间的对应关系^{[11, 12]}, 后者则是通过最小化两个点集之间的拓扑结构差异来确定它们之间的对应关系^{[13, 14]}.最近, 一部分研究人员^[15-19]在传统的基于Graphs特征算法的基础上加入了学习要素, 通过在配准前使用适当的学习样本进行学习, 来优化算法中的参数设置, 从而提高了算法的配准精度.但是这类算法使用的Shape context或Graphs特征在点集的相邻点较为接近的情况下会变得非常相似, 以至于这类算法并不能达到较好的配准效果^[20-22].

基于迭代的非刚性配准算法通常由两个相互交替进行的过程组成:对应关系评估(correspondence estimation)和空间变换更新(transformation updating).与基于非迭代的算法相比, 基于迭代的算法有如下优势:它们可以在迭代过程中逐步调整源点集的初始几何形状和空间位置, 使得源点集的几何形状和空间位置变得与目标点集越来越接近, 从而使得通过几何结构特征评估点集之间的对应关系变得更加容易.TPS-RPM^[23]算法是基于迭代的算法中最著名的, 它使用点集到点集的距离、Softassign^{[24, 25]}和退火算法^[26]来对点集之间的对应关系进行评估和控制薄板样条函数(thin platespline, 简称TPS)的更新.Myronenko等人^[27]在TPS-RPM算法框架基础上提出了在空间变换更新中增加运动一致性约束条件(motion coherence constraint)来提高配准过程中空间变换的稳定性, 点集之间的对应关系则是利用最大似然估计法(maximum likelihoodestimation)来进行评估.然后, Myronenko等人^[28]在之前的文献[27]的基础上发表了著名的CPD算法(coherent points drift, 简称CPD), 他们改良了空间变换模型, 使之既可以适用于刚性和非刚性的点集配准问题, 并可以在牺牲些许配准精度的前提下, 通过使用快速高斯变换(fast Gaussian transform)^[29]和矩阵低秩逼近(low-rank matrix approximation)^[30]技术减少计算量来提升算法的配准速度.在这之后, Jian等人^[21]提出了一种基于高斯混合模型(Gaussian mixture model)的非刚性点集配准算法(被命名为GMMREG算法).该算法不直接在几何空间中配准两个点集, 而是通过将两个点集转变成两个高斯混合模型来进行对应关系的评估, 并基于最小化两个高斯混合模型的L2距离进行空间变换的更新.最近, 国内的Ma等人^[31]提出了一种基于Shape context特征的对应关系评估并采用L2最小化评估(L2- minimizing estimate)^[32]进行空间变换更新, Wang等人^[33]通过使用不对称高斯模型捕捉空间点集的不对称分布, 并用其作为特征描述进行点集的非刚性配准.

TPS-RPM、CPD、GMMREG、Ma等人和Wang等人分别使用TPS、运动一致性理论和L2距离最小化约束条件进行空间变换的更新, 分别使用欧式距离、高斯混合模型、Shape context和不对称高斯模型来进行对应关系的评估.然而, 基于单一的全局特征的对应关系评估, 使这些算法无法清晰地描述出点集中的局部特征, 这使得它们评估出的对应关系在全局特征较为相似的情况下不够可靠.

在本文中, 我们提出了一种基于双特征高斯混合模型和双约束空间变换的非刚性点集配准算法.这种算法通过分析双特征高斯混合模型的分布情况来评估两个点集之间的对应关系, 并基于高斯径向基函数对源点集进行全局与局部结构双约束的空间变换更新.与当前典型的其他基于迭代的算法相比, 本文算法主要包含了3个创新之处.

(1) 我们定义了一个双特征描述子, 并以此改进高斯混合模型.通过对高斯混合模型的改进, 使其可以基于两个特征进行对应关系的评估, 并使用确定性退火技术, 合理地调节两个特征之间的比重.

(2) 我们定义了一个局部特征描述子, 该描述子利用局部相邻点的和向量, 能够敏感地描述局部特征结构.

(3) 我们定义了一个局部结构约束项, 此约束项通过最大化空间变换前后的局部结构相似性, 维护点集的局部结构在空间变换更新过程中的稳定性.

本文第1节将对基于双特征高斯混合模型和双约束空间变换的非刚性点集配准算法进行详细的介绍并阐述其核心思想.第2节对本文算法的相关算法进行讨论.第3节将呈现本文算法在6种不同配准模式下的性能表现, 并与8种当前典型的算法进行比较实验.第4节对本文进行总结, 并提出我们对本文算法的改进设想.

1 方法

在本部分中, 我们首先介绍了基于高斯混合模型进行非刚性点集配准的算法, 然后提出本文算法对上述算法的改进方案, 以此突显本文算法的创新之处, 并在最后详细介绍了本文算法进行非刚性点集配准的完整过程及其可行性的验证.假设A=(a₁, a₂, …, a_X)^T和B=(b₁, b₂, …, b_Y)^T是两组需要进行非刚性配准的点集, 其中, A是源点集, B是目标点集.

1.1 基于高斯混合模型的非刚性点集配准

高斯混合模型在非刚性点集配准中的核心思想是:将源点集内的每个点都当作一个高斯模型的中心, 把整个源点集构造为一个高斯混合模型, 并通过调整高斯混合模型的参数, 使目标点集在高斯混合模型中的概率密度最小.如此反复迭代, 直到达到最大迭代次数或源点集与目标点集足够接近.具体的配准过程可分为两部分.

a) 构建高斯混合模型并评估两个点集之间的对应关系

通过把源点集的每一个点都当作一个高斯模型(Gaussian model, 简称GM)的中心, 使整个源点集被构造为一个高斯混合模型(Gaussian mixture model, 简称GMM).在本文中, 我们用点集A来构造高斯混合模型.于是, 点集B在高斯混合模型中的概率密度函数(probability density function, 简称PDF):

$ \mathcal{P}\mathcal{D}\mathcal{F}(B|\mathcal{G}) = \prod\limits_{j = 1}^Y {\mathcal{P}\mathcal{D}\mathcal{F}({b_j}|\mathcal{G})} = \prod\limits_{j = 1}^Y {\sum\limits_{i = 1}^{X + 1} {P(x)\mathcal{P}\mathcal{D}\mathcal{F}({b_j}|{\varrho _i})} } \\ = \prod\limits_{j = 1}^Y {\left( {(1 - \omega )\sum\limits_{i = 1}^X {P(x)\mathcal{P}\mathcal{D}\mathcal{F}({b_j}|{\varrho _i})} + \omega \frac{1}{Y}} \right)} $

(1)

这里, 我们使用了X+1个高斯模型.其中的前X个高斯模型中, 每一个高斯模型$ {\varrho _i}$均代表了点集A中的一个点a_i.第X+1个高斯模型被用于消除冗余点的影响, 并用ω来约束这个高斯模型, 且满足0＜ω＜1.点b_j在高斯模型$ {\varrho _i}$中的概率密度函数:

$ \mathcal{P}\mathcal{D}\mathcal{F}({b_j}|{\varrho _i}) = \frac{1}{{{{(2\pi {\sigma ^2})}^{D/2}}}}\exp \left( { - \frac{{||{b_j} - \mathcal{C}({\varrho _i},{\theta _i})|{|^2}}}{{2{\sigma ^2}}}} \right) $

(2)

其中, $ {\cal C}$($ {\varrho _i}$, θ_i)是第x个高斯模型的参数取值为θ_i时的中心坐标, 即进行了空间变换更新之后的点a_i坐标; D是进行配准的点集的维度.然后, 就可以基于贝叶斯法则产生高斯模型的后验概率, 其中, 高斯模型的参数的取值等于前一次迭代的参数:

$ {\mathcal{P}^{old}}({\varrho _i}|{b_j}) = \frac{{\exp \left( { - \frac{1}{{2{\sigma ^2}}}||{b_j} - \mathcal{C}({\varrho _i},{\theta _i})|{|^2}} \right)}}{{\sum\limits_{n = 1}^X {\exp \left( { - \frac{1}{{2{\sigma ^2}}}||{b_j} - \mathcal{C}({\varrho _n},{\theta _n})|{|^2}} \right)} + \frac{\omega }{{1 - \omega }}{{(2\pi {\sigma ^2})}^{D/2}}\frac{X}{Y}}} $

(3)

由此得到的一个X×Y矩阵$ {{\cal P}_o}$就是源点集A与目标点集B之间的一组模糊对应关系.

b) 空间变换更新

步骤a)确定了两个点集之间的对应关系矩阵, 在这一步开始对源点集A进行空间变换更新.由于源点集被构造为了一个高斯混合模型, 因此只要基于对应关系矩阵找到高斯混合模型的一组最优参数, 就可以完成对A的空间变换更新.由此可以看出, 非刚性点集配准的空间变换更新被转变为了一个参数的优化过程.我们通过最小化公式(1)的负对数似然函数(negative log-likelihood function)来寻找高斯混合模型最优的参数Θ和σ².公式(1)的负对数似然函数如下:

$ E = - \log \prod\limits_{j = 1}^Y {\sum\limits_{i = 1}^{X + 1} {P(x)\mathcal{P}\mathcal{D}\mathcal{F}({b_j}|{\varrho _i})} } = - \sum\limits_{j = 1}^Y {\log \sum\limits_{i = 1}^{X + 1} {P(x)\mathcal{P}\mathcal{D}\mathcal{F}({b_j}|{\varrho _i})} } $

(4)

我们采用期望最大化算法(expectation maximization, 简称EM)来解决这个优化问题, EM算法分为求解期望的E步骤和最大化期望的M步骤.

● 通过E步骤, 可以得到如下能量方程:

$ Q(\mathit{\Theta } ,{\sigma ^2}) = \sum\limits_{j = 1}^Y {\sum\limits_{i = 1}^X {{\mathcal{P}^{old}}({\varrho _i}|{b_j})} } \frac{{||{b_j} - \mathcal{C}t({\varrho _i},{\theta _i})|{|^2}}}{{2{\sigma ^2}}} + \frac{{{N_p}D}}{2}\log {\sigma ^2} $

(5)

其中, $ {N_p} = \sum\limits_{j = 1}^Y {\sum\limits_{i = 1}^X {{\mathcal{P}^{old}}({\varrho _i}|{b_j})} }$.

● 然后进行EM算法的M步骤, 即通过对能量方程式(5)求偏导数来计算最优参数Θ和σ².

为了方便对能量方程求导, 我们将其重新写为矩阵形式.在这之前, 我们首先基于高斯径向基函数(Gaussian radial basis function, 简称GRBF)对源点集A的空间变换做了如下定义:

$ {\cal C}({\cal G},\mathit{\Theta } ) = A + GW $

(6)

其中, G是由GRBF:g_ij=exp(2ξ^(-2)‖a_i-a_j‖²)计算得到的X×X内核矩阵, W是高斯内核的X×D权重矩阵.因此, 我们将要进行优化的参数由Θ和σ²变为了W和σ².通过公式(6)将公式(5)改写为

$ Q(W,{\sigma ^2}) = \sum\limits_{j = 1}^Y {\sum\limits_{i = 1}^X {{\mathcal{P}^{old}}({\varrho _i}|{b_j})} } \frac{{||{b_j} - ({a_i} + {g_{i, \cdot }}W)|{|^2}}}{{2{\sigma ^2}}} + \frac{{{N_p}D}}{2}\log {\sigma ^2} $

(7)

然后, 将能量方程式(7)写为矩阵形式:

$ \begin{gathered} Q = \frac{1}{{2{\sigma ^2}}}(Trace({B^T}{\mathcal{P}_B}B) - 2Trace({A^T}{\mathcal{P}_o}B) - 2Trace({W^T}G{\mathcal{P}_o}B) + Trace({A^T}{\mathcal{P}_A}A) + \hfill \\ {\text{ }}2Trace({W^T}G{\mathcal{P}_A}A) + Trace({W^T}G{\mathcal{P}_A}GW)) + \frac{{{N_p}D}}{2}\log {\sigma ^2} \hfill \\ \end{gathered} $

(8)

其中, Trace(·)是某个矩阵的迹, $ {{\cal P}_A}$(X×X)是由列向量$ {{\cal P}_o}$1所形成的对角矩阵, $ {{\mathcal{P}}_{B}}$(Y×Y)是由列向量$\mathcal{P}_o^T$1所形成的对角矩阵, 1代表的是一个元素全是1的列向量.

最后, 下述这组偏导数方程组的解就是进行空间变换更新所需要的最优参数:

$ \frac{{\partial Q}}{{\partial W}} = \frac{1}{{{\sigma ^2}}}G{\mathcal{P}_A}A - \frac{1}{{{\sigma ^2}}}G{\mathcal{P}_o}B + \frac{1}{{{\sigma ^2}}}G{\mathcal{P}_A}GW = 0 $

(9)

$\left. \begin{align} &\frac{\partial Q}{\partial {{\sigma }^{2}}}=-\frac{1}{{{\sigma }^{4}}}(Trace({{B}^{T}}{{\mathcal{P}}_{B}}B)-2Trace(A{{\mathcal{P}}_{o}}B)-2Trace({{W}^{T}}G{{\mathcal{P}}_{o}}B)+\\ &Trace({{A}^{T}}{{\mathcal{P}}_{A}}A)+ \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2Trace({{W}^{T}}G{{\mathcal{P}}_{A}}A)+Trace({{W}^{T}}G{{\mathcal{P}}_{A}}GW))+\frac{{{N}_{p}}D}{2{{\sigma }^{2}}}=0 \\ \end{align} \right\} $

(10)

至此就完成了基于高斯混合模型的非刚性点集配准过程中的一次迭代过程, 只要反复迭代进行a)步骤和b)步骤, 直到达到最大迭代次数或公式(5)收敛, 就完成了对两个点集的非刚性配准.

本文算法在此非刚性配准算法的两个步骤均做出了改进.

a) 定义了一种双特征描述子和局部特征描述子来改进评估对应关系时仅使用了单一的全局特征, 并在后文对改进方法进行了详细介绍;

b) 在能量方程式(5)上添加了一个全局结构约束项和一个局部结构约束项, 并在后文对局部结构约束项做了详细介绍.

1.2 对应关系评估的改进

在这一节, 我们使用双特征描述子代替原始高斯混合模型上的单一特征描述, 使其同时包含两种特征, 并通过确定性退火技术来平衡这两种特征之间的关系, 以改进对应关系评估方法.我们分为两节对改进方法进行了详细介绍.

1.2.1 双特征高斯混合模型

我们首先定义了双特征描述子, 具体如下:

$ F(\cdot )=\mathcal{F}{{(\cdot )}_{1}}+\alpha \mathcal{F}{{(\cdot )}_{2}} $

(11)

其中, $ {{\mathcal{F}}_{1}}$为第1种特征的描述子; $ {{\mathcal{F}}_{2}}$为第2种特征的描述子; $ {{\mathcal{F}}_{2}}$的权重参数α被用来平衡两个特征之间的关系, 定义如下:

$ \alpha = \exp \left( { - \frac{t}{{{c_1}}}} \right),{c_1} \in \mathbb{R} $

(12)

其中, t是当前迭代次数.由此可以看出, 随着迭代次数的增加, α会逐渐地趋于0.即随着配准的进行, 第2种特征在双特征描述子中产生的影响会逐渐减小.常数c₁被用来控制减小的速度.在本文中, 我们使用全局特征作为第1种特征, 局部特征作为第2种特征.在配准前期, 目标点集由于产生较大级别非刚性形变, 影响了全局特征描述的稳定性, 本文算法通过混合使用全局特征与稳定性较强的局部特征以解决此问题; 随着配准的进行, 源点集与目标点集的空间位置逐渐变得相似, 全局特征描述的稳定性逐渐恢复.本文通过确定性退火技术逐渐降低混合特征中局部特征的权重, 使全局特征在配准后期占主要地位.

为了使高斯混合模型可以包含两种特征, 我们使用双特征描述子把普通高斯模型的概率密度函数改写为如下形式:

$ \left. \begin{align} &\mathcal{P}\mathcal{D}{{\mathcal{F}}_{\varrho }}=\omega \exp \frac{{{[(\mathcal{F}{{(B)}_{1}}+\alpha \mathcal{F}{{(B)}_{2}})-(\mathcal{F}{{(A)}_{1}}+\alpha \mathcal{F}{{(A)}_{2}})]}^{2}}}{2{{\sigma }^{2}}} \\ &\ \ \ \ \ \ \ \ \ \ \ =\omega \exp \left( \frac{{{[(\mathcal{F}{{(B)}_{1}}-\mathcal{F}{{(A)}_{1}})+\alpha (\mathcal{F}{{(B)}_{2}}-\mathcal{F}{{(A)}_{2}})]}^{2}}}{2{{\sigma }^{2}}} \right) \\ &\ \ \ \ \ \ \ \ \ \ \ =\omega \exp \left( \frac{{{({{\mathit{\Psi } }_{1}}+\alpha {{\mathit{\Psi } }_{2}})}^{2}}}{2{{\sigma }^{2}}} \right) \\ \end{align} \right\} $

(13)

其中, Ψ₁表示两个点集在第1种特征上的X×Y差异矩阵, Ψ₂表示两点集在第2种特征上的X×Y差异矩阵.

通过公式(13), 把公式(3)改写为如下形式:

$ {\mathcal{P}^{old}}({\varrho _i}|{b_j}) = \frac{{\exp \left( { - \frac{1}{{2{\sigma ^2}}}{{({\psi _{1ij}} + \alpha {\psi _{2ij}})}^2}} \right)}}{{\sum\limits_{n = 1}^X {\exp \left( { - \frac{1}{{2{\sigma ^2}}}{{({\psi _{1nj}} + \alpha {\psi _{2nj}})}^2}} \right)} + \frac{\omega }{{1 - \omega }}{{(2\pi {\sigma ^2})}^{D/2}}\frac{X}{Y}}} $

(14)

至此, 高斯混合模型经过双特征描述子的改进, 基于两种特征加强了它对点集特征的描述.我们为经过改进的新型高斯混合模型命名为双特征高斯混合模型(daul-feature Gaussian mixture model).

本文中, 全局特征描述采用了原始高斯混合模型的全局特征(global feature), 即点的空间坐标.为了在评估对应关系时能正确评判两个点集之间全局特征的相似性, 使两个模型的空间结构最大程度地重叠, 并加强与局部特征的混合效果, 我们将两个点的欧氏距离作为全局特征差异.令Ψ₁=Ψ_G, 则矩阵中每一个元素的定义如下.

$ {{\mathit{\psi }}_{1ij}}={{\mathit{\psi }}_{Gij}}=||{{b}_{j}}-\mathcal{C}({{\varrho }_{i}},{{\theta }_{i}})|| $

(15)

为了使点集的局部特征(local feature)作为第2种特征, 我们定义了一个局部特征描述子, 令Ψ₂=Ψ_L, 并在下一小节进行了详细介绍.

1.2.2 局部特征描述子

在本小节中, 我们基于向量对局部特征描述子进行了定义.假设有一个包含了V个点的点集P, 点q_ik是点p_i的第k近的相邻点, 则点集P的第i个点的局部特征描述如下:

$ f{(P)_i} = \sum\limits_{k = 1}^K {{\mu _{ik}}{v_{ik}}} $

(16)

其中, K是相邻点的数量, v_ik是点p_i到点q_ik之间的向量, μ_ik是控制对应向量对局部特征描述贡献程度的权重参数.由于向量特征既包含了欧式距离还包含了角度的信息, 因此, 描述子采用和向量对p_i的局部特征进行描述.

由公式(16)的定义可以看出, 此描述子对局部特征的描述主要受到两个因素的影响:1)相邻点个数K的取值; 2)权重参数μ_ik的取值.

● 首先, 由于点p_i的局部特征描述主要受到它周围点的影响, 故K的取值应当顾及周围大部分会产生影响的相邻点.例如在二维的情况下, 最少在4个方向(上、下、左、右)的点可能存在较大的影响, 故进行二维的非刚性点集配准时, K应取一个大于4的值.由于本文算法可以应用在大于等于二维的任何维度, 因此K的取值应根据进行配准时所在的维度来确定.在本文实验中, K在二维情况下取值为5, 三维情况下取值为7.

● 其次, 考虑到在某种情况下, 向量v_ik的长度可能会由于过长而对局部特征描述子造成超出预期的影响, 故我们用欧几里得范数对μ_ik进行了定义, 使v_ik对局部特征描述子的影响满足长度越短影响越大, 具体定义如下:

$ {\mu _{ik}} = \frac{1}{{\sqrt {2 \pi} }}\exp \left( { - \frac{{||{p_i} - {q_{ik}}|{|^2}}}{2}} \right) $

(17)

至此, 完成了对局部特征描述子的定义.其核心思想是:某个点的局部特征可以由该点到一组指定数量的相邻点之间的向量之和来表示, 且这些向量均会受到他们本身长度的影响.

因此, 第2种特征差异矩阵Ψ_L的每个元素可以基于公式(16)写成如下形式:

$ {{\mathit{\psi }}_{Lij}}=||f{{\left( B \right)}_{j}}-f{{(\mathcal{C}(\mathcal{G},\mathit{\Theta } ))}_{i}}|{{|}^{2}} $

(18)

于是, 将全局特征当作第1种特征, 局部特征当作第2种特征, 公式(14)即可写为

$ {\mathcal{P}^{old}}({\varrho _i}|{b_j}) = \frac{{\exp \left( { - \frac{1}{{2{\sigma ^2}}}{{({\psi _{Gij}} + \alpha {\psi _{Lij}})}^2}} \right)}}{{\sum\limits_{n = 1}^X {\exp \left( { - \frac{1}{{2{\sigma ^2}}}{{({\psi _{Gnj}} + \alpha {\psi _{Lnj}})}^2}} \right)} + \frac{\omega }{{1 - \omega }}{{(2\pi {\sigma ^2})}^{D/2}}\frac{X}{Y}}} $

(19)

至此就完成了对对应关系评估的全部改进.根据改进之后的公式(19), 可以得到新的模糊对应关系矩阵$ {{\cal P}_o}$.

1.3 空间变换更新的改进方法

在这一节, 为了使本文算法可以在进行空间变换更新时维护点集的全局结构和局部结构稳定, 我们通过添加全局结构约束项和局部结构约束项重新设计了用于空间变换更新的能量方程式(7), 并对局部结构约束项的定义进行了具体描述.我们分为两小节详细介绍了本文算法在空间变换更新上的改进.

1.3.1 全局结构约束项

为了使本文算法在进行空间变换更新时可以维护点集的全局结构稳定, 我们在能量方程上加入了一个全局结构约束项——非刚性变换的正则化运算符γ($ {\cal C}$($ {\cal G}$, Θ)).为了方便对空间变换更新过程的改进, 我们给出了此约束项的具体矩阵形式:

$ R = \gamma (\mathcal{C}(\mathcal{G},\mathit{\Theta } )) = \frac{\lambda }{2}Trace({W^T}GW) $

(20)

其中, 常数λ是控制全局结构约束项约束强度的权重参数.

由此, 可以将公式(7)改进为如下形式:

$ {Q_G}(W,{\sigma ^2}) = Q(W,{\sigma ^2}) + \frac{\lambda }{2}Trace({W^T}GW) $

(21)

然后, 进行参数计算的方程组(9)和公式(10)被重新写为如下形式:

$ \frac{{\partial {Q_G}}}{{\partial W}} = \frac{{\partial Q}}{{\partial W}} + \frac{\partial }{{\partial W}}\left[ {\frac{\lambda }{2}Trace({W^T}GW)} \right] = 0 $

(22)

$ \frac{{\partial {Q_G}}}{{\partial {\sigma ^2}}} = \frac{{\partial Q}}{{\partial {\sigma ^2}}} + \frac{\partial }{{\partial {\sigma ^2}}}\left[ {\frac{\lambda }{2}Trace({W^T}GW)} \right] = 0 $

(23)

1.3.2 局部结构约束项

为了使本文算法在进行空间变换更新时可以维护点集的局部结构稳定, 我们首先基于第1.2.2节的局部特征描述子公式(16), 对局部结构约束项进行了定义:

$ L = \eta \sum\limits_{i = 1}^X {||f{{(A)}_{Li}} - f{{(A + GW)}_{Li}}|{|^2}} $

(24)

其中, η是控制局部结构约束项约束强度的权重参数.此约束项的核心思想是:通过判断点集A在空间变换前后的局部相似性, 来约束点集A在进行空间变换时的局部形变.

由于在此约束项中, 权重参数η的取值对约束强度占主要影响, 而在不同情况下和不同配准进度下, 本文算法对局部结构约束强度的需求是不同的, 因此, η的取值策略对本文算法在非刚性点集配准上的性能同样会产生一定影响.为了使本文算法的性能最大化, 我们针对大部分情况对局部结构约束项的约束强度在不同的配准进度下的需求进行了如下设计:在配准一开始, 按照预先设定的η的最大值对源点集的空间变换进行局部结构约束; 随着配准的进行, 逐渐减弱局部结构约束项的约束强度.因此, 采用确定性退火技术对η进行了定义:

$ \eta = m\exp \left( { - \frac{{t - 1}}{{{c_2}}}} \right),{\text{ }}{c_2} \in \mathbb{R} $

(25)

其中, m是预先设定的η的最大值, t是当前迭代次数, 常数c₂被用来控制约束强度变化的速度.

将局部结构约束项加在已经可以控制全局结构稳定的能量方程式(21)上, 使本文算法在进行空间变换更新时, 可以同时维护全局和局部结构稳定.在这之前, 为方便改进的过程, 我们给出了局部结构约束项的矩阵形式:

$ L=\eta Trace({{\mathcal{D}}^{T}}\mathcal{D}) $

(26)

其中,

$ \mathcal{D}=UA-U\left( A+GW \right)=-UGW $

(27)

这里, U=H(A)-KI, I是一个单位矩阵, K是相邻点的数量.假设有一个包含了V个点的点集P, 且点q_ik是点p_i第k近的相邻点, 则我们将H(P)定义为一个V×V的矩阵, 矩阵内的元素定义如下:

$ {h_{ij}} = \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{\sqrt {2 \pi} }}\exp \left( { - \frac{{||{p_j} - {p_i}|{|^2}}}{2}} \right),{\text{ }}{p_j} \in \{ {q_{ik}}\} _{k = 1}^K} \\ {0, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {p_j} \notin \{ {q_{ik}}\} _{k = 1}^K} \end{array}} \right. $

(28)

其中, K是相邻点的数量.可以看出, 这里的H(P)表示的是由点p_i的K个最近相邻点的权值所组成的稀疏矩阵.为了使矩阵W的取值满足点a_i+g_{i, ·}W的K个最近相邻点与a_i的K个最近相邻点相同, 在公式(27)中我们使用点集A的稀疏矩阵来计算点集A+GW的局部特征.

由此, 可以将能量方程式(21)改写为

$ {{Q}_{GL}}(W,{{\sigma }^{2}})={{Q}_{G}}(W,{{\sigma }^{2}})+\eta Trace({{\mathcal{D}}^{T}}\mathcal{D}) $

(29)

然后, 进行参数计算的方程组式(22)和公式(23)被重新写为

$ \frac{{\partial {Q_{GL}}}}{{\partial W}} = \frac{{\partial {Q_G}}}{{\partial W}} + \frac{\partial }{{\partial W}}[\eta Trace({\mathcal{D}^T}\mathcal{D})] = \frac{{\partial Q}}{{\partial W}} + \frac{{\partial R}}{{\partial W}} + \frac{{\partial L}}{{\partial W}} = 0 $

(30)

$ \frac{{\partial {Q_{GL}}}}{{\partial {\sigma ^2}}} = \frac{{\partial {Q_G}}}{{\partial {\sigma ^2}}} + \frac{\partial }{{\partial {\sigma ^2}}}[\eta Trace({\mathcal{D}^T}\mathcal{D})] = \frac{{\partial Q}}{{\partial {\sigma ^2}}} + \frac{{\partial R}}{{\partial {\sigma ^2}}} + \frac{{\partial L}}{{\partial {\sigma ^2}}} = 0 $

(31)

1.4 基于双特征高斯混合模型和双约束空间变换的非刚性点集配准

结合第1.2节给出的经过改进的对应关系评估方法和第1.3节给出的经过改进的空间变换更新, 在本节, 我们给出了本文算法进行非刚性点集配准的完整过程, 并在本节最后给出了本文算法的伪代码.

使用本文算法进行非刚性点集配准的过程可以具体分为反复迭代进行的两部分.

● 用公式(1)在源点集上构建高斯混合模型, 并基于双特征描述子, 用公式(19)计算出源点集和目标点集之间的模糊对应关系矩阵.

● 通过方程组(30)和公式(31)计算出本文算法对源点集进行空间变换时所需要的参数W和下一次迭代所需要的参数σ², 然后用公式(6)对源点集进行空间变换更新.

由于在介绍本文算法的双约束空间变换更新时并未验证它的可行性, 因此我们在本节给出了方程组(30)和公式(31)的解, 以说明本文算法在进行双约束空间变换更新时, 用于计算最优参数的能量方程式(2)是可行的.

由于

$ \frac{{\partial R}}{{\partial W}} = \lambda GW,\frac{{\partial L}}{{\partial W}} = 2\eta {G^T}{U^T}UGW,\frac{{\partial R}}{{\partial {\sigma ^2}}} = 0,\frac{{\partial L}}{{\partial {\sigma ^2}}} = 0 $

(32)

因此, 与公式(9)和公式(10)相结合, 可以得出方程组(30)和公式(31)的解:

$ W={{(G{{P}_{A}}G+{{\sigma }^{2}}\lambda G+2{{\sigma }^{2}}\eta {{G}^{T}}{{U}^{T}}UG)}^{-1}}(G{{P}_{o}}B-G{{P}_{A}}A) $

(33)

$ \left. \begin{gathered} {\sigma ^2} = \frac{2}{{{N_p}D}}(Trace({B^T}{\mathcal{P}_B}B) - 2Trace(A{\mathcal{P}_o}B) - 2Trace({W^T}G{\mathcal{P}_o}B) + \hfill \\ \ \ \ \ \ \ Trace({A^T}{\mathcal{P}_A}A) + 2Trace({W^T}G{\mathcal{P}_A}A) + Trace({W^T}G{\mathcal{P}_A}GW)) \hfill \\ \end{gathered} \right\} $

(34)

算法1给出了本文算法的伪代码.

算法1. 基于双特征高斯混合模型和双约束空间变换的非刚性点集配准算法.

输入:两个点集A和B.

参数设置.初始化K, m, ω, c₁, c₂和λ.

开始.迭代优化双特征高斯混合模型计划.

步骤1. 使用公式(19)评估点集A与点集B之间的模糊对应关系矩阵.

步骤2. 使用公式(33)和公式(34)更新参数W和σ².

步骤3. 使用公式(6)对点集A进行空间变换更新得到新的点集A.

结束.直到达到最大迭代次数或公式(29)收敛.

输出.配准后的点集A.

2 相关研究

当前, 基于迭代的算法中, 与本文算法较为相似的主要有7种:TPS-RPM^[23]、CPD^[28]、GMMREG^[21]、Ma- L2E^[31]、Wang等人^[33]、Yang等人^[22]和Ma等人的PRGLS^[34], 我们给出了其中相似度最高的Yang等人和PRGLS与本文算法的不同之处.

(1) Yang等人使用全局和局部混合距离(global and local mixture distance, 简称GLMD)特征对两个点集进行对应关系的评估, 并通过控制TPS函数的更新对源点集进行空间变换更新.本文算法使用双特征高斯混合模型对两个点集进行对应关系的评估, 并基于GRBF进行维护了源点集的全局和局部结构稳定的空间变换更新.虽然Yang等人和本文算法都是基于混合全局与局部双特征对两点集之间的对应关系进行评估, 但本文算法对源点集的空间变换更新是在有双约束的前提下进行的, 而Yang等人对源点集的空间变换更新是在只维护了全局结构稳定的前提下进行的.

(2) PRGLS是基于全局特征对两点集之间的对应关系进行评估, 并在维护了源点集的全局和局部结构稳定的前提下进行空间变换更新.虽然PRGLS和本文算法在进行空间变换更新时均维护了点集的全局和局部结构, 但本文算法在进行两点集之间的对应关系评估时, 同时使用了全局与局部特征, 而PRGLS仅使用了单一的全局特征.

由此可以看出本文算法与这两种相关算法之间的不同.为了测试本文算法的性能, 我们在本文的实验部分分别与这些相似算法进行了比较实验.

3 实验

我们用Matlab 2016a实现了本文算法, 并通过下述6种配准模式的实验测试了本文算法的各项性能.

1) 二维轮廓配准(2D synthetic point set);

2) 三维轮廓配准(3D synthetic point set);

3) 序列图像(CMU hotel sequence);

4) 遥感图像配准(remote sensing registration);

5) IMM人脸数据配准(IMM face database);

6) 真实图像特征点配准(Pascal 2007 Challenge datasets).

与下列当前较为流行的8种相关算法进行了性能比较实验:(1) TPS-RPM^[23], CPD^[28], GMMREG^[21], Ma- L2E^[31], PRGLS^[34]与Wang等人^[33]这5种基于迭代的算法; (2) Caetano等人^[15]与Leordeanu等人^[18]这两种基于Graph的学习算法; (3) Zhou等人^[14]这一种基于Graph的非学习算法.

3.1 实验设计

Fish 1、Fish 2、Chinese character和3D face是非刚性点集配准算法在轮廓点集配准测试中普遍使用的几个经典点集, 它们分别来自TPS-RPM^[23]和CPD^[28].本文实验首先使用这4个点集作为源点集, 并使用下面列出的人工合成方法来创建发生了非刚性形变的目标点集(由于本研究主要解决包含高级别形变、噪音、旋转的非刚性点集配准问题, 本文算法并非通用普适, 因此我们未进行经过冗余点处理(outlier)与局部结构残缺(occlusion)的配准模式实验).

(1) 形变级别(deformation degree):我们设置了8个控制点(三维的情况下是6个控制点)均匀分布在每组轮廓点集边缘(二维:上、下、左、右、左上、右上、左下、右下; 三维:上、下、左、右、前、后).每个控制点拥有上、下、左、右4个方向的自由移动以及0.2的移动步长.8个(或6个)控制点的移动顺序以及方向被随机设定.TPS空间变换被用于使用这8个(或6个)控制点使前述源点集发生形变创建新目标点集.本实验将移动控制点的数量定义为形变级别(二维和三维情况下的最大形变级别分别为8和6).

(2) 噪音比(noise level):我们通过利用均值为0且标准偏差以0.01~0.05间隔为0.01的高斯白噪声(Gaussian whitenoise)创建了5个噪音级别的目标点集.

(3) 旋转角度(rotation degree):我们以15°为间隔, 以-30°~30°对源点集进行人工旋转, 共设置了5个级别的旋转角度(在三维的情况下, 沿Z轴旋转).

其中, 方法(2)、方法(3)在人工产生噪音、旋转之前均会随机发生中等级别的形变(二维情况下为4级, 三维情况下为3级).并且上述的每一种人工合成方法, 在每一个级别都进行了100次的随机实验.

我们采用均方根误差(root mean square error, 简称RMSE)来测量轮廓配准结果的平均误差.

在开始性能测试实验与比较实验之前, 我们首先采用了二维轮廓点集Fish 1与复杂真实三维点集进行了本文算法的性能展示.此部分同时展示了本文算法基于双特征的配准较单特征的优势与本文算法的复杂真实三维点集非刚性配准能力.

3.2 性能展示

在本文中, 我们对与局部特征无关的参数设置如下:K=5(三维时K=7), ω=0.2, λ=8.对于其余参数, 我们进行了A与B两种设置——A:m=0, α=0;B:m=2, c₁=10, c₂=10.可以看出, 在A设置下, 本文算法仅使用了单一的全局特征与全局结构约束; 在B设置下, 本文算法使用了双特征与双约束.

在性能展示的第1部分, 我们采用了二维轮廓点集Fish 1, 并在两种不同移动步长下对目标点集进行了100次随机人工非刚性形变处理, 本文算法分别在单特征单约束参数设置与双特征双约束参数设置下进行了配准实验.测试统计数据(平均误差与标准偏差)展示在图 1(误差线表示了100次随机测试中平均误差的标准偏差值).本文算法在两种参数设置下均给出了正确的配准结果, 双特征双约束参数设置展现出了最优的配准性能.在之后的性能展示实验、性能测试实验与比较实验中, 我们均采取双特征双约束参数设置.

在性能展示的第2部分, 我们采用了动物狼在两种姿势下共包含3 400个点的表面点集(数据获取自http://www.vovakim.com/), 展示了本文算法在复杂三维点集非刚性配准中的配准性能.图 2展示了配准结果.其中点是源点集, 圈是目标点集.本文算法给出了正确的配准结果.

3.3 二维轮廓配准的实验结果

在本文实验的二维轮廓配准测试部分, 我们将本文算法与TPS-RPM, CPD和PRGLS这3种相关算法进行了比较实验, 并分别给出了本文算法在形变级别为7、噪音比为0.04、旋转角度为30°时的配准实例.

3.3.1 Fish 1

我们将本文算法与相关算法在点集Fish 1的配准实验中的性能测试统计数据展示在了图 3的第1行(误差线表示了100次随机测试中平均误差的标准偏差值.第1行~第3行分别是Fish 1, Fish 2和Chinese character的实验结果).这4种算法在所有实验中均展现出了准确的配准结果.本文算法在所有形变级别的实验中都给出了最优的配准结果.在目标点集含有噪音的配准实验中, PRGLS表现得更好.在目标点集发生了旋转的配准实验中, TPS-RPM表现得更好.图 4给出了本文算法的一个配准实例(其中, 十字是源点集, 圈是目标点集).

3.3.2 Fish 2

我们将本文算法与相关算法在点集Fish 2的配准实验中的性能测试统计数据展示在了图 3的第2行.这4种算法在所有实验中均展现出了准确地配准结果.本文算法在所有形变级别、噪音比和旋转角度的实验中都给出了最优的配准结果.图 5给出了本文算法的一个配准实例(其中, 十字是源点集, 圈是目标点集).

3.3.3 Chinese character

我们将本文算法与相关算法在点集Chinese character的配准实验中的性能测试统计数据展示在了图 3的第3行.这4种算法在所有实验中均展现出了准确地配准结果.本文算法在所有形变级别和旋转角度为-30°~15°的实验中都给出了最优的配准结果.在目标点集含有噪音的配准实验中, PRGLS表现得更好.图 6给出了本文算法的一个配准实例(其中, 十字是源点集, 圈是目标点集).

3.4 三维轮廓配准的实验结果

在本文实验的三维轮廓配准测试部分, 我们测试了本文算法在进行三维轮廓点集配准时的性能, 并与两种相关算法进行了比较实验.我们使用了点集3D face来进行性能测试和比较实验, CPD和GMMREG等算法同样使用了这个三维人造点集来测试性能.我们将本文算法与CPD, GMMREG在配准实验中的性能测试统计数据展示在了图 8(误差线表示了100次随机测试中平均误差的标准偏差值).本文算法和CPD在所有实验中均给出了准确地配准结果.本文算法在所有形变级别、噪音比和旋转角度的实验中都给出了最优的配准结果.图 7给出了本文算法分别在形变级别为5、噪音比为0.04、旋转角度为30°时的一个配准实例(其中, 十字是源点集, 圈是目标点).

3.5 序列图像的配准结果

在本文实验的序列图像配准部分, 我们测试了本文算法在进行序列图像特征点配准时的性能.与人造二维、三维点集相比, 序列图像拥有更少的特征点, 这些点稀疏地分布在图像中.CMU系列的序列图像是目前用于测试基于Graph的学习算法最流行的实验数据.在本文实验中, 我们采用了CMU系列的CMU hotel序列图像, 这组序列图像有101幅图组成, 每幅图拥有30个标记的特征点.本实验使用正确配准点数的百分比(即配准率)进行误差的评估.

本文算法与当今流行的两种基于Graph的学习算法和一种基于迭代的算法分别在这组序列图像中进行了所有配准可能(10 201次配准)的比较实验.我们将实验结果展示在了表 1(表中的S和U分别代表“监督学习”和“无监督学习”, S和U内的数字代表用于训练的图像对数量).与GMMREG、Leordeano等人和Caetano等人相比较, 本文算法展现出了更高的配准率.图 9展示了本文算法的一个配准实例, 使用的是第001号序列图像和101号序列图像(图上的黄色的连线代表正确的匹配, 蓝色的连线代表错误的匹配).

3.6 遥感图像配准

在本文实验的遥感图像配准部分, 我们使用了两对分辨率分别为951×651和785×734的谷歌地球卫星图像, 拍摄地点分别是夏威夷和湄公河.我们用原始的卫星图像(the reference image)和同一地点但发生了视角变化的卫星图像(the sensed image)进行特征点配准, 以此测试我们算法的性能, 并与GMMREG进行了比较实验.

● 特征点选取方法:使用SIFT算法选取卫星图像中的特征点并剔除冗余点, 图 10展示了我们使用的卫星图像, 并用点分别标出了SIFT特征点(左边一列标出的是目标点集, 右边一列标出的是源点集).

● 误差测量方法:我们使用TPS形变函数将发生了视角变化的卫星图像进行变形, 并将变形后的图像(the trasformed image)和原始的图像平均分为5×5个部分, 交替显示在棋盘图像(checkboard)上, 通过观察棋盘图像与原始卫星图像的相似度来评判配准算法的性能, 并从棋盘图像和原始卫星图像上手动选取10对界标点(landmark), 通过计算这10对界标点之间的平均距离作为配准结果的平均误差.我们将本文算法和GMMREG的配准结果展示在了图 11(前两列是本文算法的配准结果, 后两列是GMMREG的配准结果), 将配准结果的平均误差展示在了表 2(我们同时测试了GMMREG在遥感图像配准上的性能来与本文算法进行比较实验).在两个实验中, 本文算法均给出了较为准确的配准结果, 并且展现出了更高的配准精度.

3.7 IMM人脸数据配准

在本文实验的人脸数据配准部分, 我们测试了本文算法在在进行人脸特征点配准时的性能.我们选择了可以公开获取的标准IMM人脸数据库(IMM face database)(数据获取自http://www.imm.dtu.dk/~aam/datasets/datasets.html)^[35]对本文算法进行测试, 此数据库包含了40个人脸的图像组, 并且每个人脸图像组均包含了6种不同表情、姿势和照明度的图像, 每张图像分辨率为640×480.此数据库的图像均用58个特征点把人脸的眉毛、眼睛、鼻子、嘴巴和下巴结构标注了出来.本部分实验的目标是配准同一个人的不同表情、姿势和照明度的人脸图像中的特征点集.我们从此数据库中挑选了前8个人脸, 并对每种人脸的所有可能配准图像进行了配准实验.同时测试了Wang等人和PRGLS与本文算法进行了比较实验, 并使用配准率来衡量配准精度.我们在图 12展示了IMM人脸数据库的两个人脸的12张图像作为示例(每一列均是一个人不同表情、姿势、照明度的人脸图像, 其中的白点代表了人脸结构上的特征点), 并用白点标出了特征点.在表 3中展示了本文算法、Wang等人和PRGLS的比较实验结果(我们同时测试了Wang等人和PRGLS在IMM人脸数据库中的性能来与本文算法进行比较实验), 本文算法展现出了更高的配准精度.图 13展示了本文算法的配准实例, 第1行~第3行分别是IMM人脸数据库中的第2种人脸、第3种人脸和第5种人脸, 其中, 前两列是第1个形态和第4个形态的特征点配准实例, 后两列是第1个形态和第6个形态的特征点配准实例(十字代表源点集, 圈代表目标点集).

3.8 真实图像特征点的配准结果

在本文实验的真实图像配准部分, 我们使用Leordeanu等人^[18]的测试数据测试了本文算法的性能.这组测试数据从Pascal 2007 Challenge数据库中挑选出了30对汽车图像与20对摩托车图像, 并且每对图像包含30到60个特征点.在真实图像特征点配准中, 本文算法与CPD、GMMREG、Zhou等人的FGM^[14]和Leordeanu等人^[18]进行了比较实验, 实验结果展示在表 4.对于FGM和Leordeanu等人, 我们报告了他们公布的实验结果.本文算法给出了最优的配准率, 并在图 14展示了本文算法分别在汽车和摩托车上的两个配准实例(第1行是在Pascal 2007 Challenge汽车数据上的配准实例, 第2行是在Pascal 2007 Challenge摩托车数据上的配准实例.其中, 连线代表正确的匹配).

4 结论

我们已经详细介绍了一种基于双特征高斯混合模型的非刚性点集配准方法:(1)定义了一种基于和向量的局部特征描述子和双特征描述子, 并基于双特征描述子构建了双特征高斯混合模型, 以此进行点集的对应关系评估; (2)基于局部特征描述子定义了可以维护局部结构稳定的局部结构约束项, 与全局结构约束项相结合设计了双约束空间变换更新; (3)使用确定性退火技术来控制局部特征描述子和局部结构约束项在配准的不同时期的影响程度.在本文的实验部分, 我们将本文算法与当前流行的8种算法在人造数据、序列图像和多种不同现实图像中进行了性能测试和比较实验, 并且本文算法在绝大多数的形变和旋转情况下均展现出了最好的配准结果.关于本文算法的改进方案, 我们提出如下两种设想:(1)在双特征描述子的基础上合理地添加新的特征描述, 使对两点集之间的对应关系评估更加准确; (2)使用混合特征描述将高斯模型改进为多维高斯模型.