CN104143210B - 一种多尺度法向特征点云配准方法 - Google Patents

Abstract

本发明涉及一种多尺度法向特征点云配准方法。读入由点云获取设备采集的目标点云和源点云两视角点云；计算点的三个尺度半径邻域的曲率，从目标点云和源点云中抽取关键点；计算每一个关键点在不同尺度半径邻域的法向量角度偏差和曲率，构成关键点的特征描述子，得到目标点云关键点特征向量集和源点云关键点特征向量集；根据关键点特征描述子的相似程度，确定目标点云关键点和源点云关键点的对应关系；去除错误的对应关系，获得精确对应关系；对精确对应关系使用聚类分析方法进行精简，得到分布均匀的对应关系；对最终的对应关系进行奇异值分解得到刚体变换矩阵。本发明特征计算数据量少，计算效率和计算精度高，为实时高精度点云配准奠定了基础。

Description

一种多尺度法向特征点云配准方法

技术领域

本发明涉及一种多尺度法向特征点云配准方法。

背景技术

三维点云数据配准技术一直是逆向工程、计算机视觉、曲面质量检测及摄影测量学等领域的研究热点与难点。随着点云数据获取方法越来越简单，成本也越来越低，点云处理技术在文物保护、考古研究、3D影视、医学、工业等领域得到不断的实践和发展。在三维测量过程中,受被测物体空间位置、几何外形和测量方式的影响，单视角测量受到限制，设备仪器需从不同视角对物体定位测量,并且要将多视角测量的结果整合到同一坐标下,点云配准方法的选择尤为重要。主流的配准方法是通过计算点的特征(比如FPFH,Fast PointFeature Histogram)来寻找对应点，从而获得多视角点云之间的旋转平移关系；还有Besl等人提出的迭代最近点(ICP，Iterative Closest Point)以及在此基础上的改进算法，通过迭代采样的方式计算使点之间的欧氏距离达到最小的旋转平移关系。FPFH特征需要计算任意点以及它的k紧邻点之间两两相互连接的角度特征，计算量比较大，计算效率较低。ICP对原始点云初始位置要求比较高。并且主流的配准方法在抗噪声性能，配准精度和计算效率方面都存在不足。

发明内容

本发明目的在于提供一种多尺度法向特征点云配准方法，配准精度和计算效率高。

实现本发明目的技术方案：

一种多尺度法向特征点云配准方法，其特征在于:

步骤1：读入由点云获取设备采集的目标点云和源点云两视角点云；

步骤2：计算点的三个尺度半径邻域的曲率，根据目标函数，从目标点云和源点云中抽取关键点；

步骤3：计算每一个关键点在不同尺度半径邻域的法向量角度偏差和曲率，将其作为特征分量，构成关键点的特征描述子，从而得到目标点云关键点特征向量集和源点云关键点特征向量集；

步骤4：根据关键点特征描述子的相似程度，初步确定目标点云关键点和源点云关键点的对应关系；

步骤5：去除错误的对应关系，获得精确对应关系；

步骤6：对获得的精确对应关系使用聚类分析方法进行精简，得到分布均匀的对应关系；

步骤7：对最终的对应关系进行奇异值分解得到刚体变换矩阵。

步骤2中，具体包括以下步骤，

步骤2.1：根据不同点云数据选取适当半径r，作为任意一点p的邻域计算半径，计算点云每一个点p半径为r时邻域的3×3的协方差矩阵E，并计算特征值λ₁₁、λ₁₂、λ₁₃，

E·v_j＝λ_j·v_j,j∈{1,2,3}

式中，p_i是p半径为r邻域内的点，k是点p半径为r邻域内点的数量,表示这个邻域的质心，λ_j、v₁分别是协方差矩阵的特征向量和对应的特征值，j＝1,2,3；

步骤2.2：按照步骤2.1计算点p半径为2r时邻域的协方差矩阵，获得特征值λ₂₁、λ₂₂、λ₂₃；

步骤2.3：按照步骤2.1计算点p半径为3r时邻域的协方差矩阵，获得特征值λ₃₁、λ₃₂、λ₃₃；

步骤2.4:由以上所得的特征值估算曲率，曲率c_i由如下公式所得：

由此，得到点云数据中每个点的三个半径的曲率c₁,c₂,c₃。选择满足如下关系的点为特征点，

其中，ε为关键点选择阈值。

步骤3中，每一个关键点在原点云中分别以r、2r、3r、4r、5r为半径，按照步骤2.1计算特征值，按照步骤2.4计算曲率；将最小特征值对应的特征向量作为法向量；以不同邻域求得的法向量之间的角度变化作为特征，可得到10个角度偏差作为特征分量α_ij，如下：

以n₁作为标准，其他法向量与其做差，可得到另外四个向量，计算这四个向量两两之间角度，可得到6个偏差作为特征分量β_ij，如下：

特征分量α_ij、特征分量β_ij、曲率共同组成每个关键点的特征分量。

步骤4中，设定阈值，寻找目标点云中关键点的特征描述子在源点云中最小和次小欧氏距离的点并计算最小距离与次小距离的比值，小于阈值则认定为对应关系；

目标点云中关键点特征描述子：R_i＝(r_i1 r_i2 ... r_i21)

源点云中关键点特征描述子：S_j＝(s_j1 s_j2 ... s_j21)

任意两个特征描述子的相似性度量：

被认定的对应关系，需满足：

S_m是源点云中与R_i相似性度量最小的特征描述子，S_n是源点云中与R_i相似性度量次小的特征描述子，是最小与次小距离比值阈值，▽表示相似性度量阈值；

满足条件的特征向量R_i和S_m对应的关键点被认为是一个对应的关键点对。

步骤5中，通过如下方法实现去除错误的对应关系，获得精确对应关系，

从已获得的对应关系样本集中随机地选择出一个大小为n的样本子集S，并根据样本子集评估总样本集参数模型，计算总样本集中的所有对应关系的源点经参数模型变换后与目标点的距离偏差，如果偏差小于设定的阈值，此对应关系属于模型内样本，否则即为模型外样本，保存所有内样本；再次随机选择样本子集S，重复上述过程；每一次重复都保存内样本，并统计内样本数量，直到迭代次数达到设定值，迭代结束；之后，找到内样本个数最大的参数模型，并令其为最终参数模型，认为此时的内样本为精确的对应关系。

步骤6中，通过如下方法实现对获得的精确对应关系进行精简，得到分布均匀的对应关系；

把对应关系中所有源点分离出来组成源点集，任意选曲一个点作为聚类中心，获取半径为r的邻域内所有点，这个聚类中心和其邻域内所有点被认为属于一类；选取这一类中的点所在的对应关系中相似程性度量最小的作为这一类的对应关系，在去除这一类所有点的剩余点集中再次随机选取一个点作为聚类中心，重复上述过程，直到源点集中的点都被分类；每一类中选取的对应关系被认为最终的对应关系。

步骤7中，对最终的对应关系进行奇异值分解得到刚体变换矩阵。

对应点协方差矩阵并奇异值分解，获取R(旋转矩阵)、T(平移矩阵)。源点集M和目标点集N的质心分别为

k是对应关系的数量，p_mi和p_ni分别是是对应关系中的源点和目标点。

构造矩阵

其中，E为3×3矩阵。解E＝UΛV^T，得X＝UV^T,即得：

R＝X,T＝C_n-RC_m

矩阵TR就是两幅点云之间的刚体变换矩阵。通过刚体变换矩阵就可以将多视角点云数据共同部分重叠在一起，从而得到完整的物体模型。

本发明具有的有益效果：

本发明根据不同尺度半径所得到的曲率找到曲面变化程度比较大的点作为关键点，这些点特征比较明显，可辨识性强，能很好地表现点邻域曲面变化；同时，以关键点在原点云中的邻域计算法向量以及法向量之间的角度偏差，计算过程简单、计算数据量小。通过对欧氏距离阈值、最小和次小欧氏距离比阈值法计算的初步对应关系进行两次优化，去除了错误对应关系，并使得对应关系分布均匀，避免由于对应关系集中导致的局部配准精度高而整体精度差，使得整体配准精度得到提高。本发明点云配准方法特征计算数据量少，计算效率和计算精度高，并且能够不依赖点云初始位置、抗噪声能力强，无需二次配准，为实时高精度点云配准奠定了基础。

附图说明

图1是三组点云数据原始位置图；

图2是三种点云关键点分布图；

图3是特征向量计算图；

图4是对应关系聚类选择图；

图5是三组点云数据配准结果图；

图6是添加噪声后的配准结果图；

图7是配准过程及结果参数表。

具体实施方式

使用bunny、dragon和happy_backdrop三组点云数据给出实施方式，图1是三组点云数据初始位置，同一组点云数据中含有两个视角的点云数据。

步骤1:输入两幅点云(源点云，目标点云)。读入三维点云获取设备，如激光三角测距扫描仪采集的点云数据。点云数据必须包含从两个不同视角，并且具有重叠部分的两幅点云数据。

步骤2:计算点的三个尺度半径邻域的曲率，根据目标函数，从目标点云和源点云中抽取关键点；

步骤2.1：根据不同点云数据选取适当半径r，作为任意一点p的邻域计算半径。本发明中设定r为根据点云中点之间的最小距离2～4倍。计算点云每一个点p半径为r时邻域的3×3的协方差矩阵E，并计算特征值λ₁₁、λ₁₂、λ₁₃，

E·v_j＝λ_j·v_j,j∈{1,2,3} (2)

此处，p_i是p半径为r邻域内的点，k是点p半径为r邻域内点的数量,表示这个邻域的质心，λ_j、v₁分别是协方差矩阵的特征向量和对应的特征值，j＝1,2,3；

步骤2.2：按照步骤2.1计算点p半径为2r时邻域的协方差矩阵，获得特征值λ₂₁、λ₂₂、λ₂₃。

步骤2.3:按照步骤2.1计算点p半径为3r时邻域的协方差矩阵，获得特征值λ₃₁、λ₃₂、λ₃₃。

步骤2.4:由以上所得的特征值估算曲率，曲率c_i由如下公式所得：

由此，得到点云数据中每个点的三个半径的曲率c₁,c₂,c₃。选择满足如下关系的点为特征点。

其中，ε为关键点选择阈值。对于bunny、dragon和happy_backdro三组点云数据，阈值ε分别为0.1、0.15、0.13，见图7。图2展示了三组点云数据其中一个视角点云关键点分布情况，由图可以看出被选作关键点的点都是邻域曲面变化剧烈的点。

步骤3:计算每一个关键点在不同尺度半径邻域的法向量角度偏差和曲率，将其作为特征分量，构成关键点的特征描述子，从而得到目标点云关键点特征向量集和源点云关键点特征向量集；

每一个关键点在原点云中分别以r、2r、3r、4r、5r为半径，按照公式(1)计算协方差矩阵，求得特征值和特征向量。最小特征值对应的特征向量可以近似为法向量，按照公式(3)计算曲率。这样每一个关键点拥有五个不同法向量和五个不同的曲率。

在图3中每两个法向量之间存在角度偏差，我们以不同邻域求得的法向量之间的角度变化作为特征。可得到10个角度偏差作为特征分量α_ij，如下：

图3中n₂-n₁与n₃-n₁之间存在角度偏差β₂₃，这个偏差表明多个法向量之间三维空间里的变化情况。以n₁作为标准，其他法向量与其做差，可得到另外四个向量。计算这四个向量两两之间角度，因此，可得到6个偏差作为特征分量β_ij，如下：

曲率变化也是点的基本特征之一，每一个点计算了五个曲率，这样对于每一个关键点，我们就有了拥有21个特征分量的特征描述子。

步骤4:根据关键点特征描述子的相似程度，初步确定目标点云关键点和源点云关键点的对应关系；

设定阈值，寻找目标点云中关键点的特征描述子在源点云中最小和次小欧氏距离的点，并计算最小距离与次小距离的比值，小于阈值则认定为对应关系。三组点云最小距离与次小距离的比值阈值和距离阈值如图7所示。

目标点云中关键点特征描述子：

R_i＝(r_i1 r_i2 ... r_i21) (7)

源点云中关键点特征描述子：

S_j＝(s_j1 s_j2 ... s_j21) (8)

任意两个特征描述子的相似性度量：

被认定的对应关系，需满足：

S_m是源点云中与R_i相似性度量最小的特征描述子，S_n是源点云中与R_i相似性度量次小的特征描述子，是最小与次小距离比值阈值，▽表示相似性度量阈值；满足条件的特征向量R_i和S_m对应的关键点被认为是一个对应的关键点对。这样，源点云和目标点云的对应关系被初步确定下来。

步骤5:去除错误的对应关系，获得精确对应关系；

一对对应关系中的两个点，一个称为源点，另一个称为目标点。从已获得的对应关系样本集中随机地选择出一个大小为n(在三维点云中n≥3)的样本子集S，并根据样本子集评估总样本集参数模型(在点云数据中参数模型是指源点与目标点之间的变换矩阵)。计算总样本集中的所有对应关系的源点经参数模型变换后与目标点的距离偏差，如果偏差小于设定的阈值，此对应关系属于模型内样本，否则即为模型外样本，保存所有内样本。再次随机选择样本子集S，重复上述过程。每一次重复都保存内样本，并统计内样本数量，直到迭代次数达到设定值，迭代结束。之后，找到内样本个数最大的参数模型，并令其为最终参数模型，认为此时的内样本为精确的对应关系。

步骤6:对获得的精确对应关系使用聚类分析方法进行精简，得到分布均匀的对应关系；

如图4(a)所示，当较小时，会出现对应关系都集中在局部区域的现象，这会造成局部配准精度比较高而整体配准效果比较差。增大降低初步对应关系选择限制，这样对应关系点云整体范围内分布均匀，如图4(b)。但是对应关系太多，每一个对应关系之间的误差都会对配准精度产生影响。提出聚类的方法使得对应关系即得到精简又分布均匀。

把对应关系中所有源点分离出来组成源点集，任意选曲一个点作为聚类中心，获取半径为r的邻域内所有点。这个聚类中心和其邻域内所有点被认为属于一类。选取这一类中的点所在的对应关系中相似程性度量最小的作为这一类的对应关系。在去除这一类所有点的剩余点集中再次随机选取一个点作为聚类中心，重复上述过程，直到源点集中的点都被分类。每一类中选取的对应关系被认为最终的对应关系。如图4(c),经过聚类分选后对应关系分布比较均匀，并且对应关系数量较少。

步骤7:对最终的对应关系进行奇异值分解得到刚体变换矩阵；

求对应点协方差矩阵并奇异值分解，获取R(旋转矩阵)、T(平移矩阵)。源点集M和目标点集N的质心分别为

k是对应关系的数量，p_mi和p_ni分别是是对应关系中的源点和目标点。

构造矩阵

解E＝UΛV^T，得X＝UV^T,即得：

R＝X,T＝C_n-RC_m (13)

矩阵TR就是两幅点云之间的刚体变换矩阵。通过刚体变换矩阵就可以将多视角点云数据共同部分重叠在一起，从而得到完整的物体模型。

图5是三组点云配准结果。由配准结果可以看出使用本发明设计的方法配准结果精度比较高，图7是配准过程中参数需选择和结果列表，其中显示了配准精度、转换矩阵、配准时间等。Bunny和dragon两组点云数据量相近，用时大约5秒，happy_backdrop数据量比较多，用时20秒左右，计算效率也比较高。图6中，(a)添加10％噪声配准结果,(b)添加30％噪声配准结果,(c)添加50％噪声配准结果，分别展示了点云添加10％、30％、50％噪声的配准结果，噪声对配准结果产生了影响，但是配准结果良好。

综上可以看出，本发明配准方法无论在精度上，还是在速度上都很理想，在有大量噪声影响下的配准结果良好，是一种非常好的点云快速高精度抗噪声配准方法。

Claims (7)

1.一种多尺度法向特征点云配准方法，其特征在于:

步骤1：读入由点云获取设备采集的目标点云和源点云两视角点云；

步骤2：计算点的三个尺度半径邻域的曲率，根据目标函数，从目标点云和源点云中抽取关键点；

步骤3：计算每一个关键点在不同尺度半径邻域的法向量角度偏差和曲率，将其作为特征分量，构成关键点的特征描述子，从而得到目标点云关键点特征向量集和源点云关键点特征向量集；

步骤4：根据关键点特征描述子的相似程度，初步确定目标点云关键点和源点云关键点的对应关系；

步骤5：去除错误的对应关系，获得精确对应关系；

步骤6：对获得的精确对应关系使用聚类分析方法进行精简，得到分布均匀的对应关系；

步骤7：对最终的对应关系进行奇异值分解得到刚体变换矩阵。

2.根据权利要求1所述的多尺度法向特征点云配准方法，其特征在于：步骤2中，具体包括以下步骤，

步骤2.1：根据不同点云数据选取适当半径r，作为任意一点p的邻域计算半径，计算点云每一个点p半径为r时邻域的3×3的协方差矩阵E，并计算特征值λ₁₁、λ₁₂、λ₁₃，

$E=\frac{1}{k}\underset{i=1}{\overset{k}{\Σ}}(p_i-\stackrel{\‾}{p})\·{(p_i-\stackrel{\‾}{p})}^T$

E·v_j＝λ_j·v_j,j∈{1,2,3}

式中，p_i是p半径为r邻域内的点，k是点p半径为r邻域内点的数量,表示这个邻域的质心，λ_j、v_j分别是协方差矩阵的特征向量和对应的特征值，j＝1,2,3；

步骤2.2：按照步骤2.1计算点p半径为2r时邻域的协方差矩阵，获得特征值λ₂₁、λ₂₂、λ₂₃；

步骤2.3：按照步骤2.1计算点p半径为3r时邻域的协方差矩阵，获得特征值λ₃₁、λ₃₂、λ₃₃；

步骤2.4:由以上所得的特征值估算曲率，曲率c_i由如下公式所得：

$c_i=\frac{{\mathrm{\λ}}_{i1}}{{\mathrm{\λ}}_{i1}+{\mathrm{\λ}}_{i2}+{\mathrm{\λ}}_{i3}},i=1,2,3$

由此，得到点云数据中每个点的三个半径的曲率c₁,c₂,c₃，选择满足如下关系的点为特征点，

$\frac{|c_1-c_2|}{c_1+c_2}\≤\mathrm{\ϵ},\frac{|c_2-c_3|}{c_2+c_3}\≤\mathrm{\ϵ}$

其中，ε为关键点选择阈值。

3.根据权利要求2所述的多尺度法向特征点云配准方法，其特征在于：步骤3中，每一个关键点在原点云中分别以r、2r、3r、4r、5r为半径，按照步骤2.1计算特征值，按照步骤2.4计算曲率；将最小特征值对应的特征向量作为法向量；以不同邻域求得的法向量之间的角度变化作为特征，可得到10个角度偏差作为特征分量α_ij，如下：

${\mathrm{\&alpha;}}_{ij}=\frac{n_i\&CenterDot;n_j}{\left|n_i\right|\&CenterDot;\left|n_j\right|},\left(\begin{array}{cc}\begin{array}{c}i=1,2,3,4\\ j=2,3,4,5\end{array}& i<j\end{array}\right)$

以n₁作为标准，其他法向量与其做差，可得到另外四个向量，计算这四个向量两两之间角度，可得到6个偏差作为特征分量β_ij，如下：

${\mathrm{\&beta;}}_{ij}=\frac{(n_i-n_1)\&CenterDot;(n_j-n_1)}{|n_i-n_1|\&CenterDot;|n_j-n_1|},\left(\begin{array}{cc}\begin{array}{c}i=2,3,4\\ j=3,4,5\end{array}& i<j\end{array}\right)$

特征分量α_ij、特征分量β_ij、曲率共同组成每个关键点的特征分量。

4.根据权利要求3所述的多尺度法向特征点云配准方法，其特征在于：步骤4中，设定阈值，寻找目标点云中关键点的特征描述子在源点云中最小和次小欧氏距离的点并计算最小距离与次小距离的比值，小于阈值则认定为对应关系；

目标点云中关键点特征描述子：R_i＝(r_i1r_i2...r_i21)

源点云中关键点特征描述子：S_j＝(s_j1 s_j2 ... s_j21)

任意两个特征描述子的相似性度量：

被认定的对应关系，需满足：

$\frac{d(R_i,S_m)}{d(R_i,S_n)}<\&part;,d(R_i,S_m)<\&dtri;$

S_m是源点云中与R_i相似性度量最小的特征描述子，S_n是源点云中与R_i相似性度量次小的特征描述子，是最小与次小距离比值阈值，▽表示相似性度量阈值；

满足条件的特征向量R_i和S_m对应的关键点被认为是一个对应的关键点对。

5.根据权利要求4所述的多尺度法向特征点云配准方法，其特征在于：步骤5中，通过如下方法实现去除错误的对应关系，获得精确对应关系；

从已获得的对应关系样本集中随机地选择出一个大小为n的样本子集S，并根据样本子集评估总样本集参数模型，计算总样本集中的所有对应关系的源点经参数模型变换后与目标点的距离偏差，如果偏差小于设定的阈值，此对应关系属于模型内样本，否则即为模型外样本，保存所有内样本；再次随机选择样本子集S，重复上述过程；每一次重复都保存内样本，并统计内样本数量，直到迭代次数达到设定值，迭代结束；之后，找到内样本个数最大的参数模型，并令其为最终参数模型，认为此时的内样本为精确的对应关系。

6.根据权利要求5所述的多尺度法向特征点云配准方法，其特征在于：步骤6中，通过如下方法实现对获得的精确对应关系进行精简，得到分布均匀的对应关系；

把对应关系中所有源点分离出来组成源点集，任意选曲一个点作为聚类中心，获取半径为r的邻域内所有点，这个聚类中心和其邻域内所有点被认为属于一类；选取这一类中的点所在的对应关系中相似程性度量最小的作为这一类的对应关系，在去除这一类所有点的剩余点集中再次随机选取一个点作为聚类中心，重复上述过程，直到源点集中的点都被分类；每一类中选取的对应关系被认为最终的对应关系。

7.根据权利要求6所述的多尺度法向特征点云配准方法，其特征在于：步骤7中，通过如下方法实现对最终的对应关系的奇异值分解得到刚体变换矩阵，源点集M和目标点集N的质心分别为

$C_m=\frac{1}{k}\underset{i=1}{\overset{k}{\&Sigma;}}{p}_{mi},C_n=\frac{1}{k}\underset{i=1}{\overset{k}{\&Sigma;}}{p}_{ni}$

k是对应关系的数量，p_mi和p_ni分别是对应关系中的源点和目标点；

构造矩阵

$E=\frac{1}{k}\underset{i=1}{\overset{k}{\&Sigma;}}(p_{mi}-C_m)(p_{ni}-C_n)$

其中，E为3×3矩阵，解E＝UΛV^T，得X＝UV^T,即得：

R＝X,T＝C_n-RC_m

$TR=\left[\begin{array}{cccc}{R}_{11}& R_{12}& R_{13}& T_{11}\\ R_{21}& R_{22}& R_{23}& T_{21}\\ R_{31}& R_{32}& R_{33}& T_{31}\\ 0& 0& 0& 1\end{array}\right]$

由此获得了旋转矩阵R、平移矩阵T，矩阵TR就是两幅点云之间的刚体变换矩阵。