CN110288640A - 基于凸密度极值的点云配准方法 - Google Patents

Abstract

本发明公开了基于凸密度极值的点云配准方法，该方法利用凸密度极值约束牙颌点云重叠区域的对应估计，再采用迭代最优变化进行牙颌点云配准。本发明相较于现有算法，利用重叠区域凸状特征，利用这个特征提出了凸密度极值来约束对应估计，最后利用迭代最优变换利用重叠区域进行牙颌点云配准。本发明的方法重叠率高，且配准精度有所提升，并大大提升了配准效果和配准效率。

Description

基于凸密度极值的点云配准方法

技术领域

本发明涉及计算机视觉技术领域，具体涉及基于凸密度极值的点云配准方法。

背景技术

随着计算机图形图像技术的发展，对牙颌进行三维数字化建模和分析，可使口腔医疗诊断更加直观精确。Hiew等研究得出在三维空间里可实现对牙颌模型多角度和无阻挡地检查咬合关系，还能在对数字模型进行截面观测，直观和准确地完成的诸如对面积、体积的测量，可确保前期的信息更加全面、准确。牙颌咬合分析扩展到三维空间上是口腔正畸牙颌咬合分析的重要趋势，能够更灵活准确地将牙颌分析研究提升到三维立体空间，进而转变为动态分析，以达到更好地分析效果。正畸医生可以更好地观测牙齿和口腔内部构造，从不同角度对咬合关系进行评估，还能根据不同指标进行准确的参数测量。研究在点云空间上的咬合分析具有重要价值和必要性。

然而，在三维空间上直观分析咬合关系的挑战在于不同治疗时期牙颌点云的配准。这是因为治疗前后牙齿和牙龈的变化明显，数据源具备高噪声、低重叠率等特点，现有技术为低重叠点云问题提供了不同的解决方案。

1、基于特征或关键点

大多数方法都基于提取特征或关键点的方法寻找对应关系。这些特征包括几何特征，或更高级别的形状描述符和特殊图案。例如，在Wehr，A等提出了使用了分段和曲线，而Dold，C等提出了使用局部平面，球体和圆柱体作为对应关系的候选。4PCS快速鲁棒匹配算法，从全局配准的角度出发，把配准问题看成一个搜索平行四边形共面四点的问题，具有抗噪和鲁棒性高的特点。

基于3D表面特征的算法涉及提取局部特征，获得两点云中特征之间的匹配，最后使用RANSAC或其他稳健估计器估计相对姿态。这些算法容易受到大部分不正确的特征匹配和重复的场景元素的影响。

然而，上述使用局部形状描述符的方法高度依赖于特征，虽然提高了效率，也扩大了配准应用场景，但是，信息的减少会导致不准确和模糊；并且当局部表面不具备特定特征或者特征不明显，由于不具备正确信息的特征描述子的使用导致匹配错误率高，而且描述子主要利用点的相邻区域或欧几里球域，这些区域对网格细分也很敏感。

2、基于点云统计特性

第一类方法为局部优化技术，然而其仅在具有小的相对变化和大的重叠率的情况下才能发挥作用，而对于重叠率很低的数据，配准数据源之间的相似性不是很明显，点对点的对应关系不一定存在，点云可能只有部分重叠，而底层对象本身通常是非凸的，导致潜在的大量对齐出现局部最小值的陷阱。

第二类方法，依赖于两点云的统计特性。基于正态分布的配准方法，建立数据的统计模型，使用标准优化技术来确定两个点云间的最优的匹配，由于利用点云密度的统计函数而不是点集，精度不高；因为其在配准过程中不利用对应点的特征计算和匹配，所以时间比其他方法快。

3、基于重叠区域定位

基于重叠区域的方法，从数据源特点出发，重点放在提取重叠区域上，去除数据中噪声的影响，找出待配准点云之间最为接近的区域用于对齐。

戴玉成等将配准分成两个阶段：粗配准阶段应用SVD方法，根据特征点对找出重叠区域；精配准阶段采用ICP算法对重叠区域求出全局最优解，然而ICP方法简单且效率高，但是当数据源重叠率不高的时候，欧式距离最近的点作为对应点的方式并不能保证正确的对应关系，即该方法的有效实施需要对应关系的估计建立在具有良好重叠率数据源的前提下。孙楠提出的方案将点云相对关系信息作为已知量，通过投影法找到重叠区域，从而降低了误匹配率。张晓娟等设计了一种扫描方法来对线的曲率进行分析，同时引入遗传算法来提取点云重叠区域。沈萦华等使用k-means聚类算法来分块，根据质心组成的近似三角形关系，并基于八叉树的区域生长算法得到重叠区域。然而，上述方法在配准精度、重叠率确定、效率以及配准效果方向均存在不足。

发明内容

为了解决现有技术存在的低重叠率、噪声大引起配准精度和配准效果不足，以及算法复杂，耗时长，效率低等技术问题，本发明提供了解决上述问题的基于凸密度极值的点云配准方法。

本发明通过下述技术方案实现：

基于凸密度极值的点云配准方法，该方法包括以下步骤：

步骤一，获取同一病例不同时期的牙颌点云和牙颌点云的重叠区域；

步骤二，基于重叠区域的凸状特征，利用凸密度极值约束对应估计，得到重叠区域对应关系；

步骤三，基于上述对应估计得到的对应关系，利用迭代最优变换优化对牙颌点云进行配准。

优选的，所述步骤二具体包括：从重叠区域中寻找凸密度最大的区域得到WALA嵴，利用WALA嵴约束对应估计。

优选的，所述寻找凸密度最大的区域得到WALA嵴具体包括：

步骤1.1，由牙颌点云的重叠区域获取牙弓线邻域点集；

步骤1.2，计算牙弓线邻域点集内的点的平均曲率；

步骤1.3，判断该点是否为凸状特征，如果是，则执行步骤1.4，否则返回步骤1.2；

步骤1.4，计算得到不同距离位置下的凸密度值；

步骤1.5，判断该距离处的凸密度值是否为极值，如果是，则该点为WALA嵴点；否则返回步骤1.2；

步骤1.6，重复步骤1.2至步骤1.5，直到遍历完毕牙弓线邻域点集内的所有点，得到WALA嵴点集。

优选的，利用WALA嵴约束对应估计包括：

步骤2.1，根据最邻近搜索法确定待配准点云之间的对应点；

步骤2.2，将得到的WALA嵴点集作为约束条件；

步骤2.3，从对应点中剔除不符合约束条件的对应点对，得到最终的对应点集，即为重叠区域的对应关系。

优选的，所述步骤三还包括：

步骤3.1，基于步骤二得到的重叠区域的对应关系，利用迭代最优变换完成重叠区域的配准；

步骤3.2，采用上述步骤3.1重叠区域配准中所利用的旋转矩阵R和平移矩阵T完成整体牙颌点云的配准。

优选的，所述步骤3.1还包括：

步骤3.1.1，基于凸密度极值约束后，得到S_Q中与S_P存在对应关系的点集S_Qnearst,记作S_P＝{sp_i|i＝1,2...m},S_Qnearst＝{sq_i|i＝1,2...m}，设置目标函数E(R,T)：

其中，S_Q与S_P分别为待配准的两个牙颌点云P,Q的重叠区域；

步骤3.1.2，利用相关变化矩阵求解方法和基于最小二乘优化算法计算旋转矩阵R和平移矩阵T，使得上述目标函数最小；

步骤3.1.3，对sp_i使用上一步求得的旋转矩阵R和平移矩阵T进行旋转和平移变换，得到sq_i新的对应点集sp_i′＝R·sp_i+T；

步骤3.1.4，计算sq_i和对应点集sp_i'的平均距离d：

步骤3.1.5，返回步骤3.1.1进行迭代计算，直到迭代次数达到设定的最大值或d小于设定阈值。

优选的，所述步骤3.2还包括：

步骤3.2.1，待配准的两个牙颌点云P,Q利用最邻近点的对应估计方法找出Q中与P存在对应关系的点集Q_nearst，记作P＝{p_i|i＝1,2...m},Q_nearst＝{q_i|i＝1,2...m}；

步骤3.2.2，对p_i使用重叠区域对齐求得的旋转矩阵R和平移矩阵T进行旋转和平移变换，得到q_i新的对应点集p_i′＝R·p_i+T；从而获得了P旋转平移后新的矩阵P'＝{p_i'|i＝1,2...m}，得到配准后的点云P',Q。

本发明具有如下的优点和有益效果：

相较于现有方法，利用重叠区域凸状特征，利用这个特征提出了凸密度极值来约束对应估计，最后利用迭代最优变换利用重叠区域进行牙颌点云配准。本发明的方法重叠率高，且配准精度有所提升，并大大提升了配准效果和配准效率。

附图说明

此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不构成对本发明实施例的限定。在附图中：

图1为本发明的方法流程示意图。

图2为本发明的点云凸密度实验图。

图3为本发明的WALA嵴提取流程示意图。

图4为第一实施例牙颌配准效果对比图。其中，(a)表示采用ICP算法的牙颌配准效果图，(b)表示采用本发明的牙颌配准效果图。

图5为第二实施例牙颌配准效果对比图。其中，(a)表示采用ICP算法的牙颌配准效果图，(b)表示采用本发明的牙颌配准效果图。

图6为对应估计效果对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。

实施例1

本实施例提出了基于凸密度极值的点云配准方法，该方法利用凸密度极值约束牙颌点云重叠区域的对应估计，再采用迭代最优变化进行牙颌点云配准，如图1所示。本实施例提出的点云配准方法在重叠率和配准精度、配准效果等方面相较于传统算法都有明显提升。

本实施例对重叠区域特征定义如下：

空间凹凸性：通过研究发现牙齿与牙龈的边缘正是符合“人类视觉沿着主曲率的负最小值描绘了对象区域的边界”这一理论。选取病例上颌点云，计算每个点的平均曲率。可以发现牙齿与牙龈的边界线上的点集的平均曲率的方向为负方向，数值为负。

因此可以得出，平均曲率为负数的点集合，在三维视觉中呈现出凹陷的特点；同理，平均曲率为正数的点集合，在三维视觉中呈现出凸状的特点，研究发现WALA嵴的平均曲率的方向为正方向，数值为正，这正与其凸状特征相对应。从而平均曲率的大小反应了凹凸的程度，方向反应属于凹状特征还是凸状特征。所以可以利用平均曲率的方向和大小定义重叠区域的凸状特征。

点云平均曲率：本实施例利用Pottmann的基于主分量分析的积分不变量理论来得到点云表面的曲率特征。对于点集A，主分量分析主要是对重心b和协方差矩阵J(A)的分析。主曲率和主分量是可以利用协方差来对应起来的，这是因为某点所在曲面和协方差矩阵的两个特征向量方向是大致相切，这两个特征向量对应的特征值绝对值是较大的。具体来说主曲率指的是曲率的主分量，而主曲率方向用协方差矩阵的特征向量e_i(i＝1,2,3...)表示。曲率的数值则由特征值λ_i(i＝1,2,3...)来表示。

曲率邻域的选择通常有球体邻域，球面邻域，和曲面邻域。本实施例选择球体邻域，是希望最大程度上降低牙颌点云的高噪声和低重叠对结果的影响。

在p点处选一个半径值r生成包围球，设内部点个数为n_B。将包围球得到的每个点p_i投影到p点及其邻近点确定的MLS移动最小二乘平面，得到投影点p_i'，当即该点与其投影点连线与法矢方向n(x)的夹角小于90度，那么定义p_i是我们需要的球体邻域内部，反之，即在外部。选择所有被定义在内部的点，设为点集B＝{p_in}，n_Bin为内部的点的总数，那么邻域的体积可以表示为式(1)：

其中V_b是半径为r的包围球的体积。接下来进行主分量分析，根据式(2)可以得到质心b和协方差矩阵J(A)的表达式：

协方差矩阵的特征向量e_i(i＝1,2,3...)对应的特征值λ_i(i＝1,2,3...)对应了该方向上的曲率。取较大的前两个特征值，根据式(3)再由平均曲率K是主曲率k₁,k₂的平均值，进而得到平均曲率：

本实施例的点云配准方法具体包括以下步骤：

步骤一、基于凸密度极值的对应估计

因为凸状特征由平均曲率表示，那么凸密度就可以用平均曲率密度来表示，之所以要定义凸密度是因为它能反应凸状特征的变化，从而得到重叠区域中凸状特征最为明显的地方。而重叠区域的凸状特征在WALA嵴区域最为明显，通过寻找凸密度最大的区域得到WALA嵴，从而得到重叠区域中凸状特征最为明显的区域。然后利用WALA嵴约束对应估计，从而使得待配准点云之间的对应关系更为准确。

1.1凸密度极值

WALA嵴为前牙弓下方牙龈区域中具备最多凸状特点的连续区域，对应的点云特征描述就是前牙弓远离牙齿方向的连续邻域点集的一个子集，其具备凸密度最大的特点。

首先对本实施例涉及到的概念进行介绍：

(1)点云的连续定义：用空间点云相关技术来表达，就是如果某点集及其邻域具备连续特征，则可以认为这样的点集及其邻域能够表征三维表面连续的特点。其中离散点云的连续特征定义，设三维空间里两个点p_i＝(x_i1,x_i2,x_i3),q_j＝(y_j1,y_j2,y_j3)，如式(4)所示，其欧式距离小于或等于三维激光扫描精度下的最小距离ε，这样就可以认为三维空间里这两点是连续的，是可以近似表征三维表面的连续特点。

并且提取的前牙弓线的点云点集满足点云的连续定义，其邻域的点集也满足这样的定义，所以前牙弓的点云点集及远离牙齿方向的邻域点集能够满足WALA嵴定义中关于连续牙龈区域的要求。

(2)远离牙齿方向的邻域：设Q＝{q_j|j＝1,2...m},P＝{p_i|i＝1,2...n}，表示牙弓线点云Q共有m个点，牙龈点云P共有n个点。由式(5)可以得到每个p_i都有一个d_inearst和q_inearst与之对应，分别表示该点与牙弓线的最短距离和牙弓线上的对应点。

d_inearst＝d₁,(s.t.d_j≤d_j+1)

q_inearst＝q₁ (5)

其中d_inearst表示p_i计算它到点云Q每个点之间的距离中最近的一个距离。同理可以得到点云Q中每个点q_j的对应点p_jnearst。

类似骨架线区分颊侧牙龈线和舌侧牙龈线的方法，由骨架线曲线拟合作为分界线可以区分P点云的颊侧和舌侧，从而P的颊侧处于牙齿点云的下方，所以牙弓线点云Q在P的颊侧邻域是远离牙齿方向的。

(3)点云的凸状特征：由上述重叠区域特征定义中的三维空间中一个点如果平均曲率方向为正，则可以认为该点具备凸状的特点，其平均曲率的大小反应凸状特点的明显与否，即平均曲率如果为正方向，数值越大则凸状特征越明显。

我们令λ_i＝d_inearst，另外设每个p_i的平均曲率为K_i，则每个牙龈点云上的点都有两个特征数值λ_i,K_i，分别表示了该点距离牙弓线的距离特征和该点的凹凸特征。即向着远离牙齿方向，距离前牙弓线的距离为λ_i的点，如果K_i＞0，则表示该处具有凸状特征。而确定属于哪种凹凸特性，可以通过定义w_i变量来确定，如式(6)：

(4)点云的凸密度：随着远离牙齿方向，距离的增加，牙弓线点云及其邻域内点的集合越来越大，统计集合内点的正平均曲率密度，可以表示随着距离的增加，点云的凸状特征密度的变化。对{λ_i,K_i}表示的点p_i根据λ_i进行升序排序得到新的集合P满足P＝{p_i|λ₁≤λ₂≤...≤λ_n,i＝1，2...n}。

则根据式(7)点云的凸密度β描述为：

其中N表示符合λ_i＜θ条件下i的最大取值，即计算凸密度的点集内的元素个数；θ表示确定距离位置的参数，通过改变θ，可以得到不同距离位置下的凸密度数值。

如图2所示的点云凸密度实验图，横坐标是距离λ_i，纵坐标是凸密度β，可以发现随着距离的增加，点云的凸密度特性出现先减小后增大的现象，出现一个凸密度峰值后便逐步减小到一个稳定值。

结合实验现象做出如下分析：

(1)减小阶段：这是因为牙弓线附近的点大多处于牙齿和牙龈的交界区域，由“人类视觉沿着主曲率的负最小值描绘了对象区域的边界”，可知这些区域曲率为负值，使得凸密度计算公式中的w_i大多数为0，所以凸密度很小；

(2)增大阶段：随着距离的增加，点集合逐渐并入了WALA嵴区域的点，这些区域的点的凸状特征明显，使得凸密度很大。

(3)稳定阶段：距离进一步加大后，逐渐远离了WALA嵴区域，新并入的点并不具备明显的凸状特征，甚至部分新并入的点不具备凸状特征，使得凸密度下降。

综上所述，重叠区域的特征为凸状特征，通过实验分析凸密度，可以发现凸密度能够反映凸状特征的变化。

根据凸密度公式，当距离参数θ使得凸密度β(θ)取得最大值的时候，N满足WALA嵴点集的条件，记为N_s，也就是S＝{p_i|λ_i≤λ_i+1,i＝1，2...N_s-1}为WALA嵴点集，如式(8)：

WALA嵴区域提取流程如图3所示。

1.2对应估计约束

对应估计仅仅采用最邻近点算法存在缺陷，因为并不是所有最近点都是正确的对应关系，这些误点会降低配准的准确性。而对应估计约束就是根据某种约束条件对已有的对应关系进行分析，不符合约束条件的点对需要删除，对符合约束条件的点对进行保留。一个好的约束条件将使得对应估计更为精确，提升配准算法的抗噪能力。

本实施例提出利用凸密度极值作为约束条件是指，如果对应关系能最大程度上的使得WALA嵴作为对应点集合，那么根据WALA嵴的生物恒定性，治疗前后这部分区域的变化最小，对这种对应关系进行刚性变换，将在最大程度上保证配准达到期望的结果。

设待配准的两个点云的重叠区域S_P,S_Q，基于凸密度极值的对应估计约束算法如下：

(1)计算待配准点云之间的最邻近点。假设将S_P的点集合数量m作为对应点对的个数，根据最邻近搜索找到在S_Q中与之对应的点，将这些点集设为S_Qnearst，每一个对应点对设为(x_i,y_i),i＝1，2，...m。

(2)计算凸密度极值得到WALA嵴点集。根据凸密度公式，当距离参数θ使得凸密度β(θ)取得最大值的时候，得到WALA嵴点集，设S_P的WALA嵴点集为W_P。设S_Q的WALA嵴点集为W_Q。

(3)对应估计约束。从对应点(x_i,y_i),i＝1，2，...m剔除不符合约束条件的对应点对。其中约束条件是如果x_i∈W_P，那么检查y_i是否属于W_Q，当时候，删除这对对应点对(x_i,y_i)；当y_i∈W_Q时候，保留这对对应点对(x_i,y_i)。假设经过约束条件之后约束点对数量m'，则得到最终的对应点对(x_i',y_i'),i＝1，2，...m'。

基于凸密度极值的对应估计约束，剔除了错误对应点对，保证了WALA嵴点集为W_P和W_Q一定成为对应点集，基于WALA嵴的生物恒定性，这样的对应估计更为准确。

步骤二、利用迭代最优变换进行牙颌点云配准

本实施例提出的算法是利用凸密度极值约束对应估计之后配准，根据WALA嵴的相关理论，重叠区域能够具有治疗前后变化小的特点，而凸状特征使得对应估计更为准确。将重叠区域配准得到的旋转平移矩阵用于整个牙颌点云的配准，在配准精度尤其是重叠率上效果明显。

2.1重叠区域配准

设待配准的两个点云P,Q，它们的重叠区域，分别令作S_P,S_Q。将他们对齐就是确定这两个重叠区域的对应估计和迭代最优化变换。

基于凸密度极值的对应估计约束后，得到S_Q中与S_P存在对应关系的点集S_Qnearst，记作S_P＝{sp_i|i＝1,2...m},S_Qnearst＝{sq_i|i＝1,2...m}，其中sp_i,sq_i属于一组基于欧式距离最近的对应点。

然后进行迭代最优化变换，其目的就是通过不断迭代来最优化一个目标函数，从而确定一个旋转矩阵R和平移矩阵T，从而使得待配准的两个点云的平均欧式距离最近，就认为完成了对齐。这个目标函数也就是误差函数.

具体算法流程如下：

(1)S_P＝{sp_i|i＝1,2...m},S_Qnearst＝{sq_i|i＝1,2...m}，也就是本实施例凸密度极值约束的对应估计得到的对应关系。式(9)定义误差函数E(R,T)为：

(2)利用相关变换矩阵求解方法和基于最小二乘的优化算法计算旋转矩阵R和平移矩阵T，使得误差函数最小；

(3)对sp_i使用上一步求得的旋转矩阵R和平移矩阵T进行旋转和平移变换，得到sq_i新的对应点集sp_i′＝R·sp_i+T；

(4)计算sq_i和对应点集sp_i'的平均距离d：

(5)迭代计算需要回到第(1)个步骤，直到迭代次数达到设定的最大值或d小于设定阈。

关于参数的选择：

算法的参数：一个是最近点作为对应点需要确定多小距离以内下两点最近；另外一个是迭代的最大迭代次数；还有一个目标函数结果的最小阈值。这些参数的选择，与实际的工程应用相关，例如应用场景是激光扫描的牙颌点云配准，数据精度是0.0001，而最小距离参数可以根据实际数据来确定，但不能小于0.0001；而迭代的最大迭代次数和目标函数结果的最小阈值都属于迭代终止条件，可以根据数据和实验来确定。

2.2牙颌点云配准

整体牙颌的配准方法为将上述2.1中重叠区域对齐所利用到的旋转矩阵R和平移矩阵T作为整体牙颌对齐的刚性变换矩阵。

(1)设待配准的点云P,Q同样利用最邻近点的对应估计方法可以找出Q中与P存在对应关系的点集：Q_nearst，记作：P＝{p_i|i＝1,2...m},Q_nearst＝{q_i|i＝1,2...m}，其中p_i,q_i属于一组基于欧式距离最近的对应点。

(2)对p_i使用重叠区域对齐求得的旋转矩阵R和平移矩阵T进行旋转和平移变换，得到q_i新的对应点集p_i′＝R·p_i+T；从而获得了P旋转平移后新的矩阵P'＝{p_i'|i＝1,2...m}，得到配准后的点云P',Q。

实施例2

本实施例采用上述实施例1提出的基于凸密度极值的点云配准方法与采用传统ICP配准方法在重叠率和配准精度进行实验分析。本实施例中重叠率采用最大公共点集LCP来衡量，配准精度采用均方根误差RMS进行衡量。本实施例给出的实验结果均基于Matlab实验环境，在CPU Intel(R)Core(TM)i3-6100@3.70GHz和内存8G的PC机上运行，操作系统为Windows64位操作系统。

数据源为正畸牙科的青少年牙颌点云，利用ICP算法和上述实施例1提出的算法比较治疗前下颌和治疗后下颌的配准效果。

具体实验分析过程如下：

1、衡量指标

重叠率评价指标为最大公共点集LCP。重叠率指治疗前后点云重叠部分点集个数占整体点云点集个数的比例。设待配准的点云为P和Q，Q中与P存在对应关系的点集：Q_nearst，记作：P＝{p_i|i＝1,2...m},Q_nearst＝{q_i|i＝1,2...m}，其中p_i,q_i属于一组基于欧式距离最近的对应点，定义两点重叠为p_i,q_i的欧式距离小于一个阈值α。采用最大公共点集LCP来评估重叠率，先计算每一对对应点p_i,q_i之间的欧式距离，并按照距离大小排序，将欧式距离小于阈值α下的对应点的个数作为重叠部分点集个数，进而得到重叠率。实验设置α＝0.1。

配准精度评价指标为均方根误差RMS来表示，如式(11)。对p_i进行刚性变换，得到q_i新的对应点集p_i′＝R·p_i+T，对于不同的旋转平移矩阵R,T，通过均方根误差函数的计算，评估配准精度，均方根误差数值越小说明对应点集合之间的平均距离越为接近，也就说明配准效果越好，配准精度越高。

2、实验分析

数据来源是某医院儿童口腔科提供的治疗前后下颌点云数据样本，通过基于凸密度极值约束对应估计，并利用迭代变换进行配准，对比治疗前后点云的配准效果。图4为样本1的下颌点云不同算法配准效果的对比图，(a)为ICP算法配准结果，(b)为本实施例算法配准结果。图5为样本2的下颌点云不同算法配准效果的对比图，(a)为ICP算法配准结果，(b)为本实施例算法配准结果。

由图可以发现，上述实施例1提出的方法由于是利用凸密度极值约束对应估计之后配准，根WALA嵴的相关理论，重叠区域能够具有治疗前后变化小的特点，而凸状特征使得对应估计更为准确。将重叠区域配准得到的旋转平移矩阵用于整个牙颌点云的配准，在配准精度尤其是重叠率上效果明显。

本实施例针对配准精度、重叠率、算法耗时的实验的相关参数与结果如表1所示，本实施例算法和ICP算法的对比结果如表2所示。

表1本实施例算法指标表

表2配准算法指标对比表

由上表可以看到，本发明在算法耗时上有比较大的改进，原因在于迭代变化估计的每次迭代都需要重新估计对应关系，本发明的方法的对应估计来自于治疗前后重叠区域的对应估计，对应点对远少于ICP算法选取某一点云所有点进行对应点估计。

重叠率效果明显。重叠率的计算是满足预设间距在阈值α之下的对应点对的数量和总对应点对数量的比值，ICP算法得到的结果中，对应点对满足这样的点对并不多，而本发明提出的算法使得满足这样阈值条件的点对得到很大提升；

在配准精度上，相较于ICP算法，本发明提出的算法也有所改进。

下面分析算法在对应估计方面的差异，以样本1的数据为例，分析不同算法获得到的前5006个对应点，将对应点之间的欧式距离按照等距离分为12组，欧式距离依次增大，横坐标表示距离，纵坐标表示落在不同距离分组里面的对应点的数量，如图6所示(图中对应不同距离的两组直方图，左侧表示采用ICP算法，右侧表示采用本发明提出的算法)。

本发明利用凸密度极值约束对应估计进而改进配准方法，可以发现在前面第1、2、3的距离分组中，本发明算法的对应点对数量多于ICP算法，说明在欧式距离较小的对应点对的数量多于ICP算法。对应点对距离上越接近就说明点云重叠度高，验证了本发明的算法配准后的重叠率高的现象。另外在较大的距离分组，如10、11、12组，本实施例算法得到的对应关系点对并没有落在这些区域，说明本实施例算法配准效果比较好，不存在距离过远的对应关系点对。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.基于凸密度极值的点云配准方法，其特征在于，该方法包括以下步骤：

步骤一，获取同一病例不同时期的牙颌点云和牙颌点云的重叠区域；

步骤二，基于重叠区域的凸状特征，利用凸密度极值约束对应估计，得到重叠区域对应关系；

步骤三，基于上述对应估计得到的对应关系，利用迭代最优变换优化对牙颌点云进行配准。

2.根据权利要求1所述的基于凸密度极值的点云配准方法，其特征在于，所述步骤二具体包括：从重叠区域中寻找凸密度最大的区域得到WALA嵴，利用WALA嵴约束对应估计。

3.根据权利要求2所述的基于凸密度极值的点云配准方法，其特征在于，所述寻找凸密度最大的区域得到WALA嵴具体包括：

步骤1.1，由牙颌点云的重叠区域获取牙弓线邻域点集；

步骤1.2，计算牙弓线邻域点集内的点的平均曲率；

步骤1.3，判断该点是否为凸状特征，如果是，则执行步骤1.4，否则返回步骤1.2；

步骤1.4，计算得到不同距离位置下的凸密度值；

步骤1.5，判断该距离处的凸密度值是否为极值，如果是，则该点为WALA嵴点；否则返回步骤1.2；

步骤1.6，重复步骤1.2至步骤1.5，直到遍历完毕牙弓线邻域点集内的所有点，得到WALA嵴点集。

4.根据权利要求3所述的基于凸密度极值的点云配准方法，其特征在于，利用WALA嵴约束对应估计包括：

步骤2.1，根据最邻近搜索法确定待配准点云之间的对应点；

步骤2.2，将得到的WALA嵴点集作为约束条件；

步骤2.3，从对应点中剔除不符合约束条件的对应点对，得到最终的对应点集，即为重叠区域的对应关系。

5.根据权利要求1-4任一项所述的基于凸密度极值的点云配准方法，其特征在于，所述步骤三还包括：

步骤3.1，基于步骤二得到的重叠区域的对应关系，利用迭代最优变换完成重叠区域的配准；

步骤3.2，采用上述步骤3.1重叠区域配准中所利用的旋转矩阵R和平移矩阵T完成整体牙颌点云的配准。

6.根据权利要求5所述的基于凸密度极值的点云配准方法，其特征在于，所述步骤3.1还包括：

步骤3.1.1，基于凸密度极值约束后，得到S_Q中与S_P存在对应关系的点集S_Qnearst,记作S_P＝{sp_i|i＝1,2...m},S_Qnearst＝{sq_i|i＝1,2...m}，设置目标函数E(R,T)：

其中，S_Q与S_P分别为待配准的两个牙颌点云P,Q的重叠区域；

步骤3.1.2，利用相关变化矩阵求解方法和基于最小二乘优化算法计算旋转矩阵R和平移矩阵T，使得上述目标函数最小；

步骤3.1.3，对sp_i使用上一步求得的旋转矩阵R和平移矩阵T进行旋转和平移变换，得到sq_i新的对应点集sp_i′＝R·sp_i+T；

步骤3.1.4，计算sq_i和对应点集sp_i'的平均距离d：

步骤3.1.5，返回步骤3.1.1进行迭代计算，直到迭代次数达到设定的最大值或d小于设定阈值。

7.根据权利要求5所述的基于凸密度极值的点云配准方法，其特征在于，所述步骤3.2还包括：

步骤3.2.1，待配准的两个牙颌点云P,Q利用最邻近点的对应估计方法找出Q中与P存在对应关系的点集Q_nearst，记作P＝{p_i|i＝1,2...m},Q_nearst＝{q_i|i＝1,2...m}；

步骤3.2.2，对p_i使用重叠区域对齐求得的旋转矩阵R和平移矩阵T进行旋转和平移变换，得到q_i新的对应点集p_i′＝R·p_i+T；从而获得了P旋转平移后新的矩阵P'＝{p_i'|i＝1,2...m}，得到配准后的点云P',Q。