CN107886529B - 一种用于三维重建的点云配准方法 - Google Patents

Abstract

本发明公开了一种用于三维重建的点云配准方法包括：S101、获得多视角三维点云，取两视角三维点云分别作为源点集和目标点集；S102、构建KD‑tree；S103、求取源点集中各点的法向量；S104、计算法向量夹角平均值；S105、将源点集中的点分级并设置最大分辨率，初始时分辨率为1；S106、计算当前分辨率下源点集中每一级的采样比例并提取采样点；S107、基于本发明中的匹配度求采样点在目标点集的匹配点；S108、使用四元数法计算旋转矩阵和平移矩阵；S109、对源点集变换，得到新的源点集；S110、重复步骤S107到S109直到目标函数最小；S111、若满足预设条件，则结束，否则将当前分辨率加1，并返回步骤S106。本发明利用低分辨率匹配点对迅速完成配准，高分辨率匹配点对提高精度，并利用匹配度查找匹配点，大大提高了大规模点云的配准速度及精度。

Description

一种用于三维重建的点云配准方法

技术领域

本发明属于图像处理技术领域，具体涉及一种用于三维重建的点云配准方法。

背景技术

近年来，随着三维扫描技术的日益成熟与计算机视觉技术的飞速发展，三维重建技术正广泛应用于3D打印、机器视觉、数字考古、医疗发展等诸多领域，受到了各大研究院校与科研学者的高度关注。

由于三维扫描设备有限的扫描视角，每次测量物体时只能得到物体的一部分点云数据，在实际应用过程中通常从不同视角对物体进行多次扫描后配准来实现三维物体的重现，因此，在三维物体重建过程中点云配准显得尤为重要。

当前，使用的自动配准算法主要包括基于曲面的配准、基于几何特征的配准与迭代配准算法。基于曲面的配准算法用切平面代替真正的曲面使得在迭代过程中可以使迭代误差快速收敛，但需求解非线性最小二乘的问题，时间复杂度比较高，同时，当物体表面曲率变化明显时，算法收敛性得不到保证。由于曲率具有平移、旋转不变性，在基于几何特征的配准算法中大多学者选择基于曲率特征的配准算法，该算法可以在一定程度上保证算法的配准精度，但易导致存在多个相似点对的情况，增加算法的复杂度，且对于局部特征不明显的物体不太适用。目前最经典的迭代配准算法就是Besl和Chen于1991年提出的迭代最近点(ICP)算法，该算法通过每次迭代都在目标点集中选择最近的点作为控制点来估算变换矩阵参数，直到目标函数值不变或小于条件预设阈值为止。但该算法收敛速度慢、计算量大、鲁棒性很差，同时对点云间位姿要求较高，容易陷入局部最优。

为了更快速的获得精度更高的配准结果，国内外学者对传统ICP算法做了一系列的改进与优化。如点到投影的ICP算法，提出将待配准点云中的点沿着目标点云视点方向穿过，选取与目标点云中的点相交的点作为匹配点进行配准，但该算法应用范围不广且配准精度不高；如结合法向量和曲率特征的ICP改进算法，首先利用点云法向量夹角初步提取关键点，并使用主曲率约束选取初始点集，然后使用点间距离和高斯曲率来确认精确匹配点，但该方法预设阈值较多，较难确定，且适合特征突出的点云匹配，一旦匹配点选取错误，容易造成配准失败；如结合法向量和直方图的精确配准算法，提出法向量直方图的概念，并将点的法向量直方图特征量作为点云的特征描述子来确定匹配点对，该算法避免了直接使用法向量进行匹配带来的歧义性问题，但计算量很大。

综上所述，现有的基于ICP的改进算法有如下不足之处：

(1)速度快，但精度较低；关键特征点选取越少，计算量越小，其收敛速度就越快，但精度相应较低。

(2)精度高，但速度较慢；关键特征点选取越多，精度相应提高，但由于计算量增加，速度较慢。

(3)算法的初始位置要求较高，否则易陷入局部最优。

发明内容

针对现有技术中存在的上述不足，本发明提供了一种速度较快、精度较高的三维重建的点云配准方法。

为解决上述技术问题，本发明 采用了如下的技术方案：

一种用于三维重建的点云配准方法，包括以下步骤：

S101、用三维扫描仪从不同视角扫描物体表面，获得至少两个视角三维点云，取两个视角三维点云分别作为源点集P和目标点集Q；

S102、对P和Q构建KD-tree；

S103、基于KD-tree求取P中每个点的k_1个近邻点，并计算P中每个点的法向量；

S104、计算P中每个点的法向量夹角平均值；

S105、基于P中每个点的法向量夹角平均值对P中的点分级并设置最大分辨率，初始时分辨率为1；

S106、计算当前分辨率下P中每一级的采样比例并提取采样点；

S107、基于KD-tree求P中各采样点在Q中的k_3个近邻点，并基于本发明中的匹配度求各采样点与其k_3个近邻点的匹配值，计算匹配度时若某点无曲率信息，则求取其曲率信息并存储，取匹配值最小的点作为该采样点的匹配点；

S108、根据获得的匹配点对，使用四元素法计算旋转矩阵和平移矩阵；

S109、使用获得的旋转矩阵和平移矩阵对源点集P进行变换，得到转换后的点云集P¹，并将P¹作为新的源点集P；

S110、重复步骤S107到S109直到使目标函数最小；

S111、若分辨率满足预设条件或目标函数的均方根误差满足预设条件，则结束，否则将当前分辨率加1，并返回步骤S106。

优选地，步骤S104包括：

对于源点集P中任意点p_i，使用KD-tree求取p_i的k_2个近邻点，计算P中点的法向量夹角平均值M_i，其中M_i为：

上式中

为点p_i的法向量，

为p_i近邻点的法向量，M_i值的大小反映p_i局部区域的弯曲程度，M_i值越大，p_i局部区域弯曲程度越大，M_i值越小，p_i局部区域越平坦，由此计算出P中每个点的法向量夹角平均值。

优选地，步骤S105包括：

基于源点集P中每个点的法向量夹角平均值将P中的点分为m级，则第l级点的取值范围为：min_M+(l-1)*G≤M_l＜min_M+l*G，1≤l≤m；

其中min_M为源点集P中所有点的法向量夹角平均值的最小值，max_M为P 中所有点的法向量夹角平均值的最大值，G是每级的间隔；

设最大分辨率为n，设初始时分辨率为1。

优选地，步骤S106包括：

当前分辨率为g时，1≤g≤n，第l级采样点的采样比例为：

其中count_m为第m级总点数，count_l为第l级总点数，fix为向零取整；

计算当前分辨率为g时第l级采样点数C_l,g＝count_l·R_l,g，若当前分辨率 g＝1，对第i级随机采样C_l,g个点，否则，对第l级随机采样C_l,g-C_l,g-1个点。

优选地，步骤S107包括：

对源点集P中的任意采样点p_i，使用KD-tree求取点p_i在目标点集Q中的 k_3个近邻点；

分别计算p_i和其k_3个近邻点的匹配度W(p_i,q_j)：

上式中p_i1、p_i2、p_i3、p_i4分别为点p_i的主曲率k₁、k₂、高斯曲率K、平均曲率H， q_j1、q_j2、q_j3、q_j4分别为近邻点q_j的主曲率k₁、k₂、高斯曲率K、平均曲率H；

计算匹配度时若某点无曲率信息，则使用KD-tree求取其在对应点集的 k_4个近邻点，计算曲率信息并存储；选取W(p_i,q_j)值最小的点作为p_i的匹配点。

优选地，所述预设条件为当前分辨率等于最大分辨率，或者，所述预设条件为目标函数的均方根误差满足预设阈值或迭代至收敛。

优选地，目标函数为

式中Q_i为源点集P 中的点在目标点集Q中的匹配点，P_i为源点集P中的点，N为匹配点对总数，通过迭代使目标函数的均方根误差满足预设阈值或迭代至收敛，则目标函数最小。

综上所述，本发明公开了一种用于三维重建的点云配准方法，包括以下步骤：S101、用三维扫描仪从不同视角扫描物体表面，获得至少两个视角三维点云，取两个视角三维点云分别作为源点集P和目标点集Q；S102、对P和Q构建KD-tree；S103、基于KD-tree求取P中每个点的k_1个近邻点，并计算P中每个点的法向量；S104、计算P中每个点的法向量夹角平均值；S105、基于P中每个点的法向量夹角平均值对P中的点分级并设置最大分辨率，初始时分辨率为1；S106、计算当前分辨率下P中每一级的采样比例并提取采样点；S107、基于KD-tree求P中各采样点在Q中的k_3个近邻点，并基于本发明中的匹配度求各采样点与其k_3个近邻点的匹配值，计算匹配度时若某点无曲率信息，则求取其曲率信息并存储，取匹配值最小的点作为该采样点的匹配点；S108、根据获得的匹配点对，使用四元素法计算旋转矩阵和平移矩阵；S109、使用获得的旋转矩阵和平移矩阵对源点集P进行变换，得到转换后的点云集P¹，并将 P¹作为新的源点集P；S110、重复步骤S107到S109直到使目标函数最小；S111、若分辨率满足预设条件或目标函数的均方根误差满足预设条件，则结束，否则将当前分辨率加1，并返回步骤S106。

同现有的改进的一些ICP算法相比较，本发明具有以下有益技术效果：

本发明首先引入多分辨率概念，利用低分辨率匹配点对迅速完成配准，高分辨率匹配点对提高配准精度。其次引入匹配度概念并改进计算方法，从而寻找更合适的匹配点，同时根据计算需要来计算相应点的曲率信息，以进一步提升配准效率。同时引入Logistic模型，曲线为S型曲线，具备开始增长迅速，后来增长缓慢的特点，符合点云多分辨率配准思想：先采样少数点云进行迅速配准，然后大幅增加点云采样，提升配准精度，最后小幅增加采样，使点云整体收敛，各分辨率间采样呈S型增长。

附图说明

图1为本发明公开一种用于三维重建的点云配准方法的流程图；

图2为配准前的Bunny点云模型原始图；

图3为传统ICP算法对Bunny点云配准的结果图；

图4为现有算法对Bunny点云配准的结果图；

图5为本发明公开一种用于三维重建的点云配准方法对Bunny点云配准的结果图；

图6为配准前的Dragon点云模型原始图；

图7为传统ICP算法对Dragon点云配准的结果图；

图8为现有算法对Dragon点云配准的结果图；

图9为本发明公开一种用于三维重建的点云配准方法对Dragon点云配准的结果图；

图10为本发明公开一种用于三维重建的点云配准方法相对传统ICP算法和现有算法所提高百分比的折线图。

具体实施方式

下面结合附图对本发明作进一步的详细说明。

如图1所示，本发明公开了一种用于三维重建的点云配准方法，包括以下步骤：

S101、用三维扫描仪从不同视角扫描物体表面，获得至少两个视角三维点云，取两个视角三维点云分别作为源点集P和目标点集Q；

S102、对P和Q构建KD-tree；

其中，本发明中点云KD-tree的构建步骤如下：

1)分别求点云中所有点在x,y,z维度上的方差；

2)将方差最大的维度设为split域；

3)将所有点按split域的值进行排序，取中值点作为根节点；

4)将split域的值小于中值点split域的值的点分配到左子空间，反之分配到右子空间；

5)对左子空间和右子空间分别重复步骤1)到4),直到只剩一个数据点。

S103、基于KD-tree求取P中每个点的k_1个近邻点，并计算P中每个点的法向量；

其中，本发明使用主成分分析法(PCA)来求取法向量，估计点云中p点的法向量的问题近似于估计该点切面法线的问题，即可转化为估计p点k_1邻域最小二乘拟合平面的问题，具体步骤如下：

1)假设该最小二乘平面为

点s(s₁,s₂,s₃)为该平面上一点，

为该平面的法线。则点p的邻域点p_i(p_ix,p_iy,p_iz)到

的距离

为：

s₁,s₂,s₃分别表示点s上的x,y,z坐标；a,b,c分别表示

的x,y,z坐标。

2)最小二乘平面

要求满足式

最小，

表示点p的 k_1个邻域点质心。

令：

则可求得：

即：

3)由式

可知

在

上。约束

则

可通过最小化J求解，故问题可转化为求解协方差矩阵C的特征值和特征矢量,其最小特征值所对应的特征向量即为p的法向量，其中协方差矩阵C：

上式中δ_i为p_i的权重，取值为1，此时计算出的法向量方向无法确定；

4)对

的方向进行调整，调整后即为所求p的法向量,设置视点方向

使得

的方向和

的方向一致，即：

S104、计算P中每个点的法向量夹角平均值；

S105、基于P中每个点的法向量夹角平均值对P中的点分级并设置最大分辨率，初始时分辨率为1；

最大分辨率跟配准精度有关，配准精度要求越高最大分辨率就设置得越大。

S106、计算当前分辨率下P中每一级的采样比例并提取采样点；

S107、基于KD-tree求P中各采样点在Q中的k_3个近邻点，并基于本发明中的匹配度求各采样点与其k_3个近邻点的匹配值，计算匹配度时若某点无曲率信息，则求取其曲率信息并存储，取匹配值最小的点作为该采样点的匹配点；

其中，点云曲率估计的方法如下：

对于某点h_i，利用构建的KD-tree求取其在对应点集中的k_4近邻点，记h_i的k_4个近邻点为Nbhd(h_i)，设拟合的二次曲面η方程为：

z＝f(x,y)＝ax²+by²+cx+dy+e，则在Nbhd(h_i)中，η满足使得下式取值最小：

分别对式ε²中各系数求导，并使其值为0，然后联立式ε²即可求解各系数，从而得到拟合二次曲面方程z＝f(x,y)。所以曲面的参数方程为：

P(x,y)＝(x y z,(x,y,f(x,y)))；

设

则可得：

由曲面第一和第二基本形式可得：

从而可得高斯曲率K、平均曲率H、主曲率k₁和主曲率k₂为：

其中，本发明使用曲率来描述点云中各点的局部特征，依据曲率信息来引导点云进行配准。但为进一步降低配准计算过程的时间复杂度，只有需要使用点p_i和点q_i的曲率时才计算其曲率并存储，若其曲率存在，则不进行计算；

S108、根据获得的匹配点对，使用四元素法计算旋转矩阵和平移矩阵；

匹配点对指的是求得的采样点及其匹配点，根据匹配点对求旋转矩阵及评议矩阵的方法为现有技术，因此不再赘述。

S109、使用获得的旋转矩阵和平移矩阵对源点集P进行变换，得到转换后的点云集P¹，并将P¹作为新的源点集P；

根据求取的旋转矩阵R和平移矩阵T，利用刚体变换公式P¹＝R*P+T得到转换后的点云集P¹，并将P¹作为新的源点集P。

S110、重复步骤S107到S109直到使目标函数最小；

目标函数为

式中Q_i为源点集P中的点在目标点集Q中的匹配点，P_i为源点集P中的点，N为匹配点对总数，通过迭代使目标函数的均方根误差满足预设阈值或迭代至收敛，则目标函数最小。

S111、若分辨率满足预设条件或目标函数的均方根误差满足预设条件，则结束，否则将当前分辨率加1，并返回步骤S106。

具体实施时，步骤S104包括：

对于源点集P中任意点p_i，使用KD-tree求取p_i的k_2个近邻点，计算P中点的法向量夹角平均值M_i，其中M_i为：

上式中

为点p_i的法向量，

为p_i近邻点的法向量。

具体实施时，步骤S105包括：

基于源点集P中每个点的法向量夹角平均值将P中的点分为m级，则第l级点的取值范围为：min_M+(l-1)*G≤M_l＜min_M+l*G，1≤l≤m

其中min_M为源点集P中所有点的法向量夹角平均值的最小值，max_M为P 中所有点的法向量夹角平均值的最大值，G是每级的间隔；

设最大分辨率为n，设初始时分辨率为1。

具体实施时，步骤S106包括：

当前分辨率为g(1≤g≤n)时，第l(1≤l≤m)级采样点的采样比例为：

其中count_m为第m级总点数，count_l为第l级总点数，fix为向零取整；

计算当前分辨率为g(1≤g≤n)时第l(1≤l≤m)级采样点数C_l,g＝count_l·R_l,g，若当前分辨率g＝1，对第i级随机采样C_l,g个点，否则，对第l级随机采样 C_l,g-C_l,g-1个点。

具体实施时，步骤S107包括：

对源点集P中的任意采样点p_i，使用KD-tree求取点p_i在目标点集Q中的 k_3个近邻点；

分别计算p_i和其k_3个近邻点的匹配度W(p_i,q_j)：

上式中p_i1、p_i2、p_i3、p_i4分别为点p_i的主曲率k₁、k₂、高斯曲率K、平均曲率H， q_j1、q_j2、q_j3、q_j4分别为近邻点q_j的主曲率k₁、k₂、高斯曲率K、平均曲率H；

计算匹配度时若某点无曲率信息，则使用KD-tree求取其在对应点集的 k_4个近邻点，计算曲率信息并存储；选取W(p_i,q_j)值最小的点作为p_i的匹配点。

在本发明中，曲率信息包括高斯曲率、平均曲率和主曲率。

具体实施时，所述预设条件为当前分辨率等于最大分辨率，或者，所述预设条件为目标函数的均方根误差满足预设阈值或迭代至收敛。

本发明采用斯坦福大学开放点云数据库的Bunny模型来验证效果，对Bunny 进行配准实验时，估计点云法向量、曲率以及选取目标点集配准点的k值均设为k＝8，分级m＝3，最大分辨率n＝6。

图2为Bunny(0度和45度)配准前点云模型图，图3、图4、图5分别为在传统ICP算法、现有算法和本文算法下对Bunny点云的配准结果图，表一为在不同算法下的耗时与配准误差的对比：

表1 Bunny点云模型下不同算法的配准结果

为进一步验证本文算法的有效性，使用曲面更为复杂的Dragon点云数据集进行实验，算法中估计点云法向量、曲率以及选取目标点集配准点的k值均设为k＝8，分级m＝5，最大分辨率n＝8。

图6为Dragon(0度和45度)配准前点云模型图，图7、图8、图9分别为在传统ICP算法、现有算法和本文算法下对Dragon点云的配准结果图，表二为在不同算法下的耗时与配准误差的对比：

表2 Dragon点云模型下不同算法的配准结果

为验证本文算法的执行效率，分别对不同规模的点云进行配准，其耗时如表3所示：

表3 不同规模点云耗时

根据表3，分别计算本文算法相对传统ICP算法和现有算法所提高的百分比并绘制折线图，所绘折线图如图10所示，其中横坐标为点云中的点的数量，纵坐标为提高的速度的百分比，Compared with the existing ICP表示在不同点云规模下本发明相对于现有算法速度提高的百分比的曲线，Compared with the traditional ICP表示在不同点云规模下本发明相对于传统算法速度提高的百分比的曲线；

从以上实验可以看出，无论是对Bunny点云模型还是对Dragon点云模型进行配准，传统ICP算法虽然精度较高但比较耗时；而现有算法根据点云特征精简了配准点，在查询邻域点上使用KD-tree加速，并且根据所提出的方法筛选匹配点对，一定程度上降低了迭代次数并提高了匹配点对的匹配精度，在精度和速度上都有一定提升；本发明提出的多分辨率匹配点ICP算法，根据点云特征，选取少数匹配点对快速配准，然后增加匹配点对，提高配准精度，在查询邻域点上使用KD-tree加速，在选择匹配点对上，使用所提出的匹配度选择匹配点，增加匹配点对的匹配精度，无论在精度还是速度上，较现有算法都有进一步的提升。根据本发明算法相对传统ICP算法和现有算法所提高的速度百分比所绘制的折线图(如图10)所示，可以看出，本发明算法相比传统ICP算法，速度提升了77％以上，相比现有算法，速度提升了62％以上，且随着点云规模的增大，有递增的趋势。

上仅是本发明优选的实施方式，需指出是，对于本领域技术人员在不脱离本技术方案的前提下，还可以作出若干变形和改进，上述变形和改进的技术方案应同样视为落入本申请要求保护的范围。

Claims (3)

1.一种用于三维重建的点云配准方法，其特征在于，包括以下步骤：

S101、用三维扫描仪从不同视角扫描物体表面，获得至少两个视角三维点云，取两个视角三维点云分别作为源点集P和目标点集Q；

S102、对P和Q构建KD-tree；

S103、基于KD-tree求取P中每个点的k_1个近邻点，并计算P中每个点的法向量；

S104、计算P中每个点的法向量夹角平均值；

S105、基于P中每个点的法向量夹角平均值对P中的点分级并设置最大分辨率，初始时分辨率为1；

S106、计算当前分辨率下P中每一级的采样比例并提取采样点；

S107、基于KD-tree求P中各采样点在Q中的k_3个近邻点，并基于匹配度求各采样点与其k_3个近邻点的匹配值，计算匹配度时若某点无曲率信息，则求取其曲率信息并存储，取匹配值最小的点作为该采样点的匹配点；

S108、根据获得的匹配点对，使用四元素法计算旋转矩阵和平移矩阵；

S109、使用获得的旋转矩阵和平移矩阵对源点集P进行变换，得到转换后的点云集P¹，并将P¹作为新的源点集P；

S110、重复步骤S107到S109直到使目标函数最小；

S111、若分辨率满足预设条件或目标函数的均方根误差满足预设条件，则结束，否则将当前分辨率加1，并返回步骤S106；

步骤S104包括：

对于源点集P中任意点p_i，使用KD-tree求取p_i的k_2个近邻点，计算P中点p_i的法向量夹角平均值M_i，其中M_i为：

上式中

为点p_i的法向量，

为p_i近邻点的法向量，M_i值的大小反映p_i局部区域的弯曲程度，M_i值越大，p_i局部区域弯曲程度越大，M_i值越小，p_i局部区域越平坦；由此计算出P中每个点的法向量夹角平均值；

步骤S105包括：

基于源点集P中每个点的法向量夹角平均值将P中的点分为m级，则第l级点的取值范围为：min_M+(l-1)*G≤M_l＜min_M+l*G，1≤l≤m；

其中min_M为源点集P中所有点的法向量夹角平均值的最小值，max_M为P中所有点的法向量夹角平均值的最大值，G是每级的间隔；

设最大分辨率为n，设初始时分辨率为1；

步骤S106包括：

当前分辨率为g时，1≤g≤n，第l级采样点的采样比例为：

其中count_m为第m级总点数，count_l为第l级总点数，fix为向零取整；

计算当前分辨率为g时第l级采样点数C_l,g＝count_l·R_l,g，若当前分辨率g＝1，对第i级随机采样C_l,g个点，否则，对第l级随机采样C_l,g-C_l,g-1个点；

步骤S107包括：

对源点集P中的任意采样点p_i，使用KD-tree求取点p_i在目标点集Q中的k_3个近邻点；

分别计算p_i和其k_3个近邻点的匹配度W(p_i,q_j)：

上式中p_i1、p_i2、p_i3、p_i4分别为点p_i的主曲率k₁、k₂、高斯曲率K、平均曲率H，q_j1、q_j2、q_j3、q_j4分别为近邻点q_j的主曲率k₁、k₂、高斯曲率K、平均曲率H；

计算匹配度时若某点无曲率信息，则使用KD-tree求取其在对应点集的k_4个近邻点，计算曲率信息并存储；选取W(p_i,q_j)值最小的点作为p_i的匹配点；由此选取出源点集P中各个采样点的匹配点。

2.如权利要求1所述的用于三维重建的点云配准方法，其特征在于，所述预设条件为当前分辨率等于最大分辨率，或者，所述预设条件为目标函数的均方根误差满足预设阈值或迭代至收敛。

3.如权利要求1所述的用于三维重建的点云配准方法，其特征在于，目标函数为

式中Q_i为源点集P中的点在目标点集Q中的匹配点，P_i为源点集P中的点，N为匹配点对总数，通过迭代使目标函数的均方根误差满足预设阈值或迭代至收敛，则目标函数最小。

Images