CN102034275A - 面向大尺度变形的非刚性注册方法 - Google Patents

Abstract

本发明公开了一种面向大尺度变形的非刚性注册方法。该方法的操作步骤如下：对扫描数据进行特征点提取；建立两个模型特征点之间的对应关系；全局刚性对准；全局变形优化。本发明基于显示计算对应关系与全局非线性优化的结合，通过建立高可信度的对应关系引导非刚性变形，最终完成大尺度变形下的非刚性注册。它不仅能有效地处理小尺度情况下的非刚性注册，对大尺度变形下的注册也有很好的适用性，鲁棒性强。

Description

面向大尺度变形的非刚性注册方法

技术领域

本发明涉及计算机图形学中动态三维模型重建的非刚性注册方法，特别是针对大尺度物体变形的非刚性注册。

背景技术

动态三维几何模型是随时间变化(运动或非刚性变形)的动态对象。如人体、动物、机器人、车辆、布料、植物、流体等。其研究动机可追溯到20世纪70年代，研究人员曾试图建立动态对象的逼真表示，从而改善动画制作效果。但有别于传统的离散帧动画，时空表面是三维几何模型随时间连续演化的结果。与离散帧动画相比，四维几何模型具有更加细致的动画效果和更高的灵活性。目前，其应用领域远超出了计算机动画和影视制作，可应用于虚拟现实、交互式游戏、模拟训练和机器人控制等方面。

随着扫描技术的发展与计算机性能的提高，数字几何模型获取技术可处理动态对象，它从实物出发，由不同时刻的采样数据生成不断变化的三维点云序列。新发明的深度视频摄像设备更能以20帧/秒的速率采样场景，实时记录观察对象的深度图像连续序列。这就为动态三维几何模型重建及后继编辑提供了研究条件。

注册(registration)可以理解为三维模型(或两帧深度图像)的匹配。注册是动态三维模型重建工作中的关键步骤，通过注册，可以积累扫描得到的各帧深度图像的信息，重建出完整的物体模型(可以称之为模板)，继而再通过模板拟合各帧的数据，重现真实物体的连续动作姿态，即动态三维几何模型。根据动态对象的物理属性，注册可分为刚性和非刚性注册两类。刚性注册基于欧氏全等变换(移动+旋转)合并(merge)二维或三维刚性模型，非刚性注册往往采用迭代最近点(iterative closest point，ICP)或其它变异算法完成求解。

最初的时候，科研工作者采用人工标记的方法，即预先标定部分对应关系来完成注册。后来，科研人员为了减弱约束条件，进行了不懈的努力，发展了各种各样的新方法，这些方法大致可以归为三类：

1.基于隐式曲面的注册方法

Mitra等人在2007年将模型分解为很小的近似刚体子块集合，使用四维表面的运动学属性来完成隐式注册；后来Lee和Lai进一步完成了MPU表面的三维非刚性注册；

2.非刚性ICP方法

①Amberg等人在2007年提出了具有非刚性优化步骤的ICP算法(N-ICP 1)，无需手工选择对应性。N-ICP 1在实现过程中采用点对点度量准则惩罚(penalize)了对应性的滑动；

②Pauly等人在2005年将收敛的刚性ICP算法扩展到非刚性条件下(N-ICP 2)，并集成使用了点对点和点对面两种度量准则。

③Li等人在2008年提出了一种非刚性ICP注册的全局对应性算法。通过优化Summer等人的嵌入式变形方法，通过求解源和目标扫描点云间以置信权值标示的对应性，同时确定非重叠(non-overlapping)区域和对准源点云和目标表面的变形(warping)域。

3.基于先验模板的注册方法

De Aguiar等人在2008年通过预先构造模板，对多视角扫描得到的动态数据进行重建，利用Laplacian的体和表面变形技术实现了物体连续动作的捕捉。

Hao Li等人在2009年进一步改进了空间嵌入变形方法，利用动态Graph的建立，使模板更有效地拟合各帧的扫描数据。

除了以上三类典型的方法，还有科研人员将视频处理中的光流方法用于动态三维模型重建。但是以上这些注册方法基于的前提假设都是两帧之间的物体变形尺度较小，两个物体空间位置接近。一旦变形尺度较大，以上方法将不再适用。虽然Huang等人在2008年提出了在保长前提下的非刚性注册，可以用来处理大尺度变形，不过算法在建立对应关系时对空间位置有初始要求，从而限制了算法的适用范围。

发明内容

本发明要解决的技术问题是提供适用范围广，鲁棒的非刚性注册方法，使得不仅能有效地完成小尺度变形下的注册，更能克服大尺度变形带来的挑战和困难，处理大尺度变形下的注册问题。

本发明的技术方案是：

第一步，在大尺度变形模型中，对源模型和目标模型的扫描数据进行特征点提取。采用基于多尺度(multi-scale)的slippage特征分析，为两帧模型提取一系列显著的特征点。如果一个点的邻域相对于自身是不可滑动的(no slippable motion)，则确定此点为显著的特征点，具有明显的几何特征。

对于模型表面上的一点p，其坐标和法向分别表示为p_i＝(p_ix，p_iy，p_iz)^T和n_i＝(n_ix，n_iy，n_iz)^T。刚性运动M由两部分组成：r＝(r_x，r_y，r_z)^T表示一个3x1的旋转向量(绕x，y，z坐标轴)，t＝(t_x，t_y，t_z)^T表示一个3x1的平移向量。如果在M的作用下，p的邻域内的每个点的瞬时运动都与表面相切，即满足：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left((r_i\×p_i+t_i)n_i\right)}^2=0$

则称M相对于p的邻域是一个滑动运动。(1)可认为是一个最小二乘问题进行求解：

$\underset{r_i,t_i}{\min }\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left((r_i\×p_i+t_i)n_i\right)}^2$

上式等同于求解线性方程Cx＝0，C表示目标函数二阶偏导的协方差矩阵：

$C=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left[\begin{array}{cccccc}{c}_{\mathrm{ix}}{c}_{\mathrm{ix}}& c_{\mathrm{ix}}{c}_{\mathrm{iy}}& c_{\mathrm{ix}}{c}_{\mathrm{iz}}& c_{\mathrm{ix}}{n}_{\mathrm{ix}}& c_{\mathrm{ix}}{n}_{\mathrm{iy}}& c_{\mathrm{ix}}{n}_{\mathrm{iz}}\\ c_{\mathrm{iy}}{c}_{\mathrm{ix}}& c_{\mathrm{iy}}{c}_{\mathrm{iy}}& c_{\mathrm{iy}}{c}_{\mathrm{iz}}& c_{\mathrm{iy}}{n}_{\mathrm{ix}}& c_{\mathrm{iy}}{n}_{\mathrm{iy}}& c_{\mathrm{iy}}{n}_{\mathrm{iz}}\\ c_{\mathrm{iz}}{c}_{\mathrm{ix}}& c_{\mathrm{iz}}{c}_{\mathrm{iy}}& c_{\mathrm{iz}}{c}_{\mathrm{iz}}& c_{\mathrm{iz}}{n}_{\mathrm{ix}}& c_{\mathrm{iz}}{n}_{\mathrm{iy}}& c_{\mathrm{iz}}{n}_{\mathrm{iz}}\\ n_{\mathrm{ix}}{c}_{\mathrm{ix}}& n_{\mathrm{ix}}{c}_{\mathrm{iy}}& n_{\mathrm{ix}}{c}_{\mathrm{iz}}& n_{\mathrm{ix}}{n}_{\mathrm{ix}}& n_{\mathrm{ix}}{n}_{\mathrm{iy}}& n_{\mathrm{ix}}{n}_{\mathrm{iz}}\\ n_{\mathrm{iy}}{c}_{\mathrm{ix}}& n_{\mathrm{iy}}{c}_{\mathrm{iy}}& n_{\mathrm{iy}}{c}_{\mathrm{iz}}& n_{\mathrm{iy}}{n}_{\mathrm{ix}}& n_{\mathrm{iy}}{n}_{\mathrm{iy}}& n_{\mathrm{iy}}{n}_{\mathrm{iz}}\\ n_{\mathrm{iz}}{c}_{\mathrm{ix}}& n_{\mathrm{iz}}{c}_{\mathrm{iy}}& n_{\mathrm{iz}}{c}_{\mathrm{iz}}& n_{\mathrm{iz}}{n}_{\mathrm{ix}}& n_{\mathrm{iz}}{n}_{\mathrm{iy}}& n_{\mathrm{iz}}{n}_{\mathrm{iz}}\end{array}\right]$

其中c_ik＝(p_i×n_i)_k，k＝x，y，z。

对C进行特征值分解C＝XΛX^T，与零特征值相对应的特征向量即表示了M相对于p的邻域的一个滑动运动。在实际情况中，C很有可能是一个满秩矩阵，所以只要特征值足够小，就可认为对应了一个滑动运动。假定|λ₁|≤|λ₂|≤|λ₃|≤|λ₄|≤|λ₅|≤|λ₆|是矩阵C的六个特征值，如果只考虑C特征值中最小值λ₁来进行滑动性的判别，这种分析在多尺度邻域的划分下将是不稳定，应该同时考虑最大特征值λ₆，并定义点p的slippage特征指示符φ＝|λ₁|/|λ₆|(任何条件下，最大特征值λ₆总是一个非零值)。采用指示符φ对p的邻域的滑动性进行状态估计，当φ的值接近0时，即表示p是可滑动的。

通过在多尺度的邻域下寻找slippage特征指示符φ的局部极大值，可以得到一系列顶点，这些顶点的相应邻域相对其自身认为是不可滑动的，即为所求的特征点。

第二步，建立两个模型特征点之间的对应关系。基于多尺度下的高斯曲率定义每个特征点的柱状图描述符，并根据柱状图描述符对源模型和目标模型上的特征点进行两两匹配，建立初始的候选对应关系集合。但是集合中不仅包括正确的对应关系，也不可避免地存在很多错误的对应关系，需要采取有效的策略剔除掉错误的对应关系。借鉴PROSAC(PROgressiveSAmple Consensus)思想，筛选候选对应关系集合，建立高可信度的显式对应关系。

首先，利用全局的谱裁剪算法，根据测地距为所有的对应关系构建相容矩阵M

c₀表示进行相容性判别的阈值，在本文中取c₀＝0.7，c_i，j表示基于模型表面的测地距定义的两组对应关系的相容性。对于两组对应关系：(s_i，t_i)和(s_j，t_j)，s_i，s_j表示源模型S上的顶点，t_i，t_j表示目标模型T上的顶点，如果同一模型上的两个顶点的测地距离近似相等，d_g(s_i，s_j)≈d_g(t_i，t_j)，其中d_g(s_i，s_j)、(t_i，t_j)表示两点之间的测地距离，则认为这两组对应关系是相容的，即采用测地距离对相容性进行度量。相容矩阵M的最大特征值的特征向量[w₁，w₂，w₃，...]^T的每个元素即为每组对应关系定义了一个权值。选取高权值的对应关系构成采样集合H。其次，对于H中的每三组对应关系，检查它们的相容性。如果这三组对应关系彼此相容，则将它们加入到相容集合K中(集合K初始为空)。然后依据每组对应关系的权值，从高到低依次向K中加入新的对应关系。每次加入的对应关系时，进行相容性的判别：如果待加入的对应关系与K中的每组对应关系的相容性度量大于等于0.7，则将其加入集合K，否则舍弃。这样，对于H中的每三组对应关系都能得到与之对应的集合K。集合K内的对应关系都是彼此相容的，K中的对应关系组数表明了该集合对于全局对应关系的相容能力，相容能力越强，可信度越高。最终，选取包含元素最多的集合作为最终所求高可信度的对应关系集合。

第三步，全局刚性对准。基于建立的稀疏对应关系，完成源模型和目标模型间的整体刚性对准。刚性对准为接下来的全局变形优化提供了很好的初始条件，加速了大尺度变形情况下全局变形优化过程的收敛。因为极分解不能区分映像和旋转，故本文采用奇异值分解(SVD)的方法。基于奇异值分解得到的负奇异值表示映像，正奇异值表示旋转。利用正奇异值，求解出正确的刚性旋转变换。

第四步，全局变形优化。在预处理阶段建立的稀疏对应关系的基础上，采用全局变形优化方法完成非刚性注册。与既有工作不同，本发明利用高可信度的显式对应关系以提高非刚性注册算法的正确性和效率，解决既有全局优化方法的局部最优解的不足。

步骤一，非刚性变形框架。采用空间嵌入变形框架，利用表面均匀采样得到若干变形控制节点，对源模型进行离散化。每个控制节点的局部作用空间定义了一个3x 3的矩阵R和一个3x 1的平移向量t，表示一个仿射变换。所有节点的仿射变换及对其邻近空间的作用构成了全局的非刚性变形。基于空间变换的局部一致性的考虑，指定源模型上的每个顶点只受其最近的k个节点的影响(本文中取k＝4)。如果两个节点的作用区域包含相同的顶点，则建立它们的邻接关系，即这两个节点互为相邻节点。在整个优化过程中，使用两个能量函数项来控制源模型的非刚性变形。

第一个能量函数项E_rigid，度量矩阵R与真正意义上的刚性旋转的差异，以保持模型表面的几何信息。这种约束是非线性的。E_rigid将每个节点的仿射变换的旋转度量误差相加：

$E_{\mathrm{rigid}}=\underset{i}{\mathrm{\Σ}}\mathrm{Rot}\left(R_i\right)$

其中， $\mathrm{Rot}\left(R\right)={\left(c_1^{\′}{c}_2\right)}^2+{\left(c_1^{\′}{c}_3\right)}^2+{\left(c_2^{\′}{c}_3\right)}^2+(1-c_1^{\′}{c}_1)+{(1-c_2^{\′}{c}_2)}^2+(1-c_3^{\′}{c_3}^2)$

c₁，c₂，c₃表示R_i的列向量。上式表示每个R_i需要满足以下条件：每个列向量为单位长度，同时两两正交。

因为相近节点的仿射变换存在相交重叠的作用区域，所以必须保证这些仿射变换彼此之间的一致性。第二个能量函数项E_smooth对这种一致性进行度量，以保证变形的平滑：

$E_{\mathrm{smooth}}=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}\underset{k\∈N\left(j\right)}{\mathrm{\Σ}}{\left|\right|R_j(g_k-g_j)+g_j+t_j-(g_k+t_k)\left|\right|}_2^2$

其中g_j表示节点的空间坐标，N(j)表示节点j的相邻节点。

步骤二，对对应关系进行加权处理。在整个优化的过程中，对应关系作为位置约束引导非刚性变形，构成全局变形优化过程中的第三个能量函数项：

$E_{\mathrm{cor}}=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\left|\right|{v^{\′}}_i-v_i\left|\right|}_2^2$

v′表示源模型S变形后的顶点坐标，v_i表示目标模型上与之对应的空间位置。该对应关系既包括了第三节基于几何特征建立的对应关系，也包括基于迭代最近点[18]为每个节点计算得到的对应关系。即

$E_{\mathrm{cor}}^{\mathrm{inthal}/\mathrm{ICP}}=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{w}_i{\left|\right|{v^{\′}}_i-v_i\left|\right|}_2^2$

预处理阶段建立的显式对应关系，而

则是根据通过最近点算法为每个节点建立的临时对应关系。

为了防止变形过程中产生失真，为每个对应关系都设置一个置信度量w。这是因为不论是基于几何特征还是迭代最近点算法建立的对应关系都不是完全正确的，错误的对应关系将会对全局变形产生影响。所以，采用这种置信度策略可以提高整个优化过程对错误的对应关系的容错能力。由于加入了置信度量，全局变形优化中定义了一个新的能量函数项E_conf：

$E_{\mathrm{conf}}=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\left|\right|1-w_i\left|\right|}_2^2$

若w_i的值接近1，则表示该对应关系对全局非刚性变形具有高可信度的引导；相反，w_i接近0则隐含的表示了该对应关系会引起E_rigid和E_smooth的显著增加，在变形过程中可以忽略该对应关系。能量函数项E_conf促使w_i向1逼近，达到源模型与目标模型的最优化匹配。

步骤三，对全局优化函数进行数值求解。综合上述各个子项，全局变形优化的能量函数E_total定义为：

$E_{\mathrm{total}}={\mathrm{\α}}_{\mathrm{rigid}}{E}_{\mathrm{rigid}}+{\mathrm{\α}}_{\mathrm{smooth}}{E}_{\mathrm{smooth}}+{\mathrm{\α}}_{\mathrm{inthal}}{E}_{\mathrm{cor}}^{\mathrm{inthal}}+{\mathrm{\α}}_{\mathrm{ICP}}{E}_{\mathrm{cor}}^{\mathrm{ICP}}+{\mathrm{\α}}_{\mathrm{conf}}{E}_{\mathrm{conf}}$

对应每个节点g_j，相应的仿射变换的R中有9个变量，t中有3个变量，而每组对应关系都有1个相应的置信变量，所以整个优化过程中一共有13n+m个未知变量。n表示变形控制节点的数目，m表示对应关系数目。

总体能量函数E_total的最优化是一个非线性问题，对应着全局非刚性变形相关参数的求解。采用高斯-牛顿方法进行迭代求解。首先设置相关参数的初始值：每个节点对应的仿射变换中R＝I，t＝0；相应能量项的权重为a_rigid＝1000，a_smooth＝100，a_conf＝100，a_inthal＝10。在每次迭代中，动态更新控制权重的值以避免局部最优解，即a_rigid、a_smooth、a_conf和a_inthal的值更新为原值的1/2，值得注意的是，a_conf的更新会在a_conf＜1的时候停止，a_inthal的更新会在a_inthal＜0.1的时候停止。α_ICP在整个优化过程中保持固定值0.1。同一般的非线性优化问题一样，全局变形优化需要经过多次迭代求解才能达到最终的非刚性注册。本方法设定整个优化过程在a_rigid＜0.01或是迭代次数大于50的时候停止。

与现有技术相比，采用本发明可以达到以下技术效果：

本发明基于显示计算对应关系与全局非线性优化的结合，通过建立高可信度的对应关系引导非刚性变形，最终完成大尺度变形下的非刚性注册。它不仅能有效地处理小尺度情况下的非刚性注册，对大尺度变形下的注册也有很好的适用性，鲁棒性强。

附图说明

图1是本发明总体算法流程图；

图2是本发明第二步对应关系的加权示意图。

具体实施方式

本发明提出的面向大尺度变形的非刚性注册方法，图1是本发明总体算法流程图，具体过程如下：

1、在大尺度变形模型中，对源模型和目标模型的扫描数据进行特征点提取，采用基于多尺度(multi-scale)的slippage特征分析，为两帧模型提取一系列显著的特征点，如果一个点的邻域相对于自身是不可滑动的(no slippable motion)，则确定此点为显著的特征点，具有明显的几何特征；

2、建立两个模型特征点之间的对应关系，筛选候选对应关系集合，建立高可信度的显式对应关系；

3、全局刚性对准，基于建立的稀疏对应关系，完成源模型和目标模型间的整体刚性对准；

4、全局变形优化，在预处理阶段建立的稀疏对应关系的基础上，采用全局变形优化方法完成非刚性注册。

Claims (5)

1.面向大尺度变形的非刚性注册方法，其特征在于该方法的具体过程为：

(1)、在大尺度变形模型中，对源模型和目标模型的扫描数据进行特征点提取，采用基于多尺度(multi-scale)的slippage特征分析，为两帧模型提取一系列显著的特征点，如果一个点的邻域相对于自身是不可滑动的(no slippable motion)，则确定此点为显著的特征点，具有明显的几何特征；

(2)、建立两个模型特征点之间的对应关系，筛选候选对应关系集合，建立高可信度的显式对应关系；

(3)、全局刚性对准，基于建立的稀疏对应关系，完成源模型和目标模型间的整体刚性对准；

(4)、全局变形优化，在预处理阶段建立的稀疏对应关系的基础上，采用全局变形优化方法完成非刚性注册。

2.根据权利要求1所述的面向大尺度变形的非刚性注册方法，其特征在于：通过在多尺度的邻域下寻找slippage特征指示符φ的局部极大值，可以得到一系列顶点，这些顶点的相应邻域相对其自身认为是不可滑动的，即为所求的特征点。

3.根据权利要求1所述的面向大尺度变形的非刚性注册方法，其特征在于筛选候选对应关系集合，建立高可信度的显式对应关系，具体步骤如下：

首先，利用全局的谱裁剪算法，根据测地距为所有的对应关系构建相容矩阵M，相容矩阵M的最大特征值的特征向量[w₁，w₂，w₃，...]^T的每个元素即为每组对应关系定义了一个权值，选取高权值的对应关系构成采样集合H；

其次，对于H中的每三组对应关系，检查它们的相容性，如果这三组对应关系彼此相容，则将它们加入到相容集合K中，集合K初始为空，然后依据每组对应关系的权值，从高到低依次向K中加入新的对应关系，每次加入的对应关系时，进行相容性的判别：如果待加入的对应关系与K中的每组对应关系的相容性度量大于等于0.7，则将其加入集合K，否则舍弃；这样，对于H中的每三组对应关系都能得到与之对应的集合K；

最终，选取包含元素最多的集合作为最终所求高可信度的对应关系集合。

4.根据权利要求1所述的面向大尺度变形的非刚性注册方法，其特征在于刚性对准采用奇异值分解(SVD)的方法，基于奇异值分解得到的负奇异值表示映像，正奇异值表示旋转，利用正奇异值，求解出正确的刚性旋转变换。

5.根据权利要求1所述的面向大尺度变形的非刚性注册方法，其特征在于利用高可信度的显式对应关系以提高非刚性注册算法的正确性和效率，解决既有全局优化方法的局部最优解的不足，具体步骤如下：

步骤一，非刚性变形框架，采用空间嵌入变形框架，利用表面均匀采样得到若干变形控制节点，对源模型进行离散化；

步骤二，对对应关系进行加权处理；

步骤三，对全局优化函数进行数值求解。

Images