CN110211129A - 基于区域分割的低覆盖点云配准算法 - Google Patents

Abstract

本发明公开了一种基于区域分割的低覆盖点云配准算法，涉及三维重建技术领域。本发明包括如下步骤：采用体积积分不变量判断其点的凹凸性，并提取凹的和凸的特征点集；对特征点构建邻接矩阵实现特征点向无向途中节点的转换，并将点的对应关系转换成无向图中的连线；利用MPPCA分析法提取特征点之间的几何关系，并构造邻接矩阵；通过SVD的ICP算法对分割区域进行配准。本发明通过体积积分不变量计算点云上点的凹凸性，并提取凹凸特征点集；再采用基于混合流形谱聚类的分割算法对特征点集进行区域分割，并利用SVD的ICP算法对区域进行配准，无需迭代即可计算刚体变换的最佳旋转矩阵，利用区域分割提高了配准区域覆盖率、配准精度和配准速度。

Description

基于区域分割的低覆盖点云配准算法

技术领域

本发明属于三维重建技术领域，特别是涉及一种基于区域分割的低覆盖点云配准算法。

背景技术

点云配准是点云数据处理的一个重要研究内容，是后续三维点云重建的技术基础。由于受环境和三维数据采集设备分辨率等多种因素的影响，在对实物数字化时采集到的每个点云数据往往只覆盖了实物模型表面的部分数据信息，而且还存在旋转和平移等刚体变换问题。因此，为了获取物体表面的完整点云数据信息，必须通过点云配准将不同角度的点云数据变换到同一坐标系下。

点云配准的算法有很较多，主要有重心重合法、标定点法以及特征提取法等。重心重合法通过将两个点云的重心重合来实现点云配准；标定点法是一种通过在物体表面贴上标定点并利用标定点进行点云配准的方法；特征提取法则是利用物体表面上的一些几何特征，如特征点、曲率、表面纹理以及边缘等，对特征进行提取和描述，通过确定特征之间的对应关系来实现点云配准。

目前应用最广泛的点云配准算法是由Besl等提出的迭代最近点(iterativeclosestpoint,ICP)算法。ICP算法计算简便直观，配准精度高，但对待配准点云的初始位置的依赖性较强，并要求两个点云之间存在一定的包含关系，因此在实际配准中容易陷入局部最优，导致配准效率不高。为了克服ICP算法的局限性，国内外学者对ICP算法进行了改进研究，相继提出PICP(probabilityICP)算法、IAICP(intensity-assistedICP)算法、MICP(modifiedICP)算法以及CICP(clusterICP)算法等多种改进的ICP算法。但是这些配准算法并不具有普适性，特别是对低覆盖率点云的配准效果不佳。

针对上述点云配准的时间复杂度高、收敛速度缓慢、对应点匹配易错等缺点，本文提出一种基于区域分割的低覆盖率点云配准算法。该算法无需迭代即可计算两个特征点集的刚体变换，具有较高的配准精度和配准速度。

发明内容

本发明的目的在于提供一种基于区域分割的低覆盖点云配准算法，通过体积积分不变量计算点云上点的凹凸性，并提取凹凸特征点集；再采用基于混合流形谱聚类的分割算法对特征点集进行区域分割，并利用SVD的ICP算法对区域进行配准，解决了现有的配准算法不具有普适性，特别是对低覆盖率点云的配准效果不佳的问题。

为解决上述技术问题，本发明是通过以下技术方案实现的：

本发明为一种基于区域分割的低覆盖点云配准算法，包括如下步骤：

步骤1：提取凹、凸特征点：采用体积积分不变量判断其点的凹凸性，并提取凹的和凸的特征点集；

步骤2：特征点区域分割：对特征点构建邻接矩阵实现特征点向无向途中节点的转换，并将点的对应关系转换成无向图中的连线；

步骤3：特征点分析：利用MPPCA分析法提取特征点之间的几何关系，并构造邻接矩阵；

步骤4：区域配准：通过SVD的ICP算法对分割区域进行配准；点云配准算法的步骤如下：

步骤41：构建测量模式矩阵并通过平移将点集的重心对其；

步骤42：将SVD用于测量模式矩阵，即并计算旋转矩阵

步骤43：创建旋转模式并根据的最近点构建

步骤44：计算配准误差若ε^(j)小于给定阈值则说明两个点云区域{A_p,B_q}配准成功，否则区域对{A_p,B_q}为不可配准的区域；

步骤45：通过穷举法将区域集合{A₁,A₂,...,A_s}和{B₁,B₂,...,B_t}中区域进行两两配准即可实现初始点云M和点云D的最终配准。

优选地，所述步骤1中，体积积分不变量V_r(p)的计算公式为：

V_r(p)＝∫_p+rQg(x)w(p-x)dx (1)；

式(1)中，g(x)表示为特征函数，w(p-x)表示权重函数，p+rQ表示为半径是r、中心是点p的球；式(1)中特征函数换成I_D(x)，并令w(p-x)＝1，则体积积分不变量V_r(p)的计算公式为：

V_r(p)＝∫_p+rQI_D(x)dx (2)；

式(2)中，当x在曲面外侧，I_D(x)＝1；当x在曲面内侧，则I_D(x)＝0。

优选地，所述步骤3中，利用EM算法得到最优描述点云分布的M个线性主成分分析器和每个点x_i对应于M个线性主成分分析器的概率密度分布，对于第m个线性主成分分析器，任意以个点x_i对应于一个降维后的点y_i，计算式为：

x_i＝V_my_i+μ_m+ε_m (4)；

式(4)中，V_m表示第m个主成分分析器的主子空间，μ_m表示分布于第m个分析器的所有点均值，ε_m表示噪声；

当y_i和ε_m均服从高斯分布，则x_i关于m的条件概率为：

式(5)中，Z＝x_i-μ_m，目标函数F_MPPCA为第m个主成份分析器上任意点x_i的最大似然估计，计算式为：

式(6)中，π_m为各个主成份分析器的混合比例，π_m≥0，

本发明具有以下有益效果：

本发明通过体积积分不变量计算点云上点的凹凸性，并提取凹凸特征点集；再采用基于混合流形谱聚类的分割算法对特征点集进行区域分割，并利用SVD的ICP算法对区域进行配准，提高了配准区域覆盖率、配准精度和配准速度。

当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有优点。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的一种基于区域分割的低覆盖点云配准算法；

图2为体积积分不变量示意图；

图3为待配准的bunny点云图；

图4为待配准的dragon点云图；

图5为bunny和dragon的公共点云的配准结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

请参阅图1所示，本发明为一种基于区域分割的低覆盖点云配准算法，包括如下步骤：

步骤1：提取凹、凸特征点：采用体积积分不变量判断其点的凹凸性，并提取凹的和凸的特征点集；

步骤2：特征点区域分割：对特征点构建邻接矩阵实现特征点向无向途中节点的转换，并将点的对应关系转换成无向图中的连线；

步骤3：特征点分析：利用MPPCA分析法提取特征点之间的几何关系，并构造邻接矩阵；

步骤4：区域配准：通过SVD的ICP算法对分割区域进行配准；点云配准算法的步骤如下：

步骤41：构建测量模式矩阵并通过平移将点集的重心对其；

步骤42：将SVD用于测量模式矩阵，即并计算旋转矩阵

步骤43：创建旋转模式并根据的最近点构建

步骤44：计算配准误差若ε^(j)小于给定阈值则说明两个点云区域{A_p,B_q}配准成功，否则区域对{A_p,B_q}为不可配准的区域；

步骤45：通过穷举法将区域集合{A₁,A₂,...,A_s}和{B₁,B₂,...,B_t}中区域进行两两配准即可实现初始点云M和点云D的最终配准。

其中，步骤1中，体积积分不变量V_r(p)的计算公式为：

V_r(p)＝∫_p+rQg(x)w(p-x)dx (1)；

式(1)中，g(x)表示为特征函数，w(p-x)表示权重函数，p+rQ表示为半径是r、中心是点p的球；式(1)中特征函数换成I_D(x)，并令w(p-x)＝1，则体积积分不变量V_r(p)的计算公式为：

V_r(p)＝∫_p+rQI_D(x)dx (2)；

式(2)中，当x在曲面外侧，I_D(x)＝1；当x在曲面内侧，则I_D(x)＝0；如图2所示，V_r(p)反应了曲面在p点处的凹凸程度，当V_r(p)＝2πr³/3时，即V_r(p)是球体积的一半时，p在半径r领域内为平面顶点，当半径r越大，则V_r(p)越大，且p+rQ内部噪声与V_r(p)无关。

顶点的凹凸程度采用点的体积积分不变量与球体积的比值来定义，定义式为：

由式(3)可知，当charV(p)＝0.5时，点p位于平坦面；当charV(p)＞0.5时，点p位于凸面；当charV(p)＜0.5时，点p位于凹面。

其中，步骤3中，利用EM算法得到最优描述点云分布的M个线性主成分分析器和每个点x_i对应于M个线性主成分分析器的概率密度分布，对于第m个线性主成分分析器，任意以个点x_i对应于一个降维后的点y_i，计算式为：

x_i＝V_my_i+μ_m+ε_m (4)；

式(4)中，V_m表示第m个主成分分析器的主子空间，μ_m表示分布于第m个分析器的所有点均值，ε_m表示噪声；

当y_i和ε_m均服从高斯分布，则x_i关于m的条件概率为：

式(5)中，Z＝x_i-μ_m，目标函数F_MPPCA为第m个主成份分析器上任意点x_i的最大似然估计，计算式为：

式(6)中，π_m为各个主成份分析器的混合比例，π_m≥0，

描述端云的最佳参数通过利用EM算法使目标函数F_MPPCA最大化来求解，将混合比例最大的主成分分析器作为点的所在流形，那么点i和点j之间的相似关系计算式为：

由此根据点之间的相似关系得到点云的邻接矩阵W，然后采用N-cut准则对点进行切割，切割的目标函数F_N-cut为：

式(8)中，A_l(l＝1,2,...,k)为第l个分割区域，为A_l的补集，通过最小化F_N-cut即可实现特征点集的区域分割。

SVD是一种应用非常广泛的矩阵分解方法，它比特征分解等其它方法更加稳定，目前己经在数据压缩、去噪以及配准等方面获得了较为广泛的应用。

设矩阵则存在m阶正交矩阵U和n阶正交矩阵V，使得其中，Σ＝diag(σ₁,σ₂,...,σ_r)，Σ为矩阵A的全部非零奇异值，满足σ₁≥σ₂≥...≥σ_r＞0；U的列向量(左奇异向量)是AA^T特征向量，V的列向量(右奇异向量)是A^TA特征向量。

SVD具有以下性质：

(1)稳定性

定理1设矩阵A,B∈R^m×n，A和B的SVD分别为λ₁≥λ₂≥...≥λ_p和τ₁≥τ₂≥...≥τ_p，其中，p＝min(m,n)，则有|λ_i-τ_i|≤||A-B||。该定理表明当矩阵A有微小扰动时，扰动前后矩阵奇异值的变化不会大于扰动矩阵的2范数。这个性质表明，对于存在噪声干扰等情况的点云，通过SVD后，点云的特征向量不会出现过大的变化。这一性质放宽了对点云预处理的要求，并使点云配准结果的准确性得到了保证；

(2)比例不变性

定理2,若举证A的奇异值为σ₁≥σ₂≥...≥σ_r≥0,举证aA(a≠0)的奇异值为σ₁′≥σ₂′≥...≥σ_s′≥0，则有

|a|(σ₁,σ₂,...,σ_s)＝(σ′₁,σ′₂,...,σ′_s) (9)；

为了消除幅度大小对特征提取的影响，所进行的归一化处理不会从本质改变奇异值的相对大小。

(3)旋转不变性

定理3若P为酉矩阵，则矩阵PA的奇异值与矩阵A的奇异值相同，即

点云的奇异值特征向量不但具有正交变换、旋转、位移、镜像映射等代数和几何上的不变性，而且具有良好的稳定性和抗噪性，广泛应用于模式识别与图形图像分析中。对点云进行奇异值分解可以得到唯一、稳定的特征描述，降低特征空间的维数，同时提高抵抗干扰和噪声的能力。

在区域配准时，若两个待配准的初始点云为M＝{m_i∈R³|i＝1,2,...,N_M}和D＝{d_i∈R³|i＝1,2,...,N_D},其中，N_M和N_D分别表示点云M和D包含的数目。首先利用第2部分的特征点提取方法获得点云M和D的凹凸特征点集，假设分别为A＝{a_i∈R³|i＝1,2,...,N_A}和B＝{b_i∈R³|j＝1,2,...,N_B}；其中，N_A和N_B分别表示特征点集A和B包含的点的数目。然后利用第3部分的区域分割算法对特征点集A和B进行区域分割，假设分割的区域分别为{A₁,A₂,...,A_s}和{B₁,B₂,...,B_t}，s和t分别表示特征点集A和B划分的区域的数目，由于划分的区域的数目远远少于特征点的数目，因此接下来采用穷举法对每一对区域{A_p,B_q}进行配准，p＝1,2,...,s，q＝1,2,...,t。

对于待配准区域对{A_p,B_q}，点a_pi∈A_p和b_pj∈B_q的相关性可表示为：

b_qj＝Ra_pi+t (11)；

式(11)中，R表示刚体变换的旋转矩阵，t表示刚体变换的平移矢量；

刚体变换(R,t)的求解可以转换为如下的最小化问题：

式(12)中,集合SE(3)是一个3D空间的欧几里得训练组，由于平移矢量t可以通过点集的重心移动来求解，因此式(12)可以写成仅依赖于旋转矩阵R的形式：

式(13)中，R∈SO(3)是一个三维旋转组，A′_p和B′_q分别定义为：

A′_p＝[a′_p1,...,a′_pN]＝A′_p{I_pN-(1/N)11^T} (14)；

B′_q＝[b′_q1,...,b′_qN]＝B′_q{I_qN-(1/N)11^T} (15)；

式(14)、(15)中，I_N是一个单元矩阵，a′_pi和b′_qi是在点集平移时分别从a_pi和b_qj中删减的3D点。

对于式(13)，采用奇异值分解(SDV)法进行求解，于是有：

式(16)中，U∈R^3×3是一个左奇异矢量矩阵，Σ∈R^3×3是一个包含奇异值的对角矩阵，V∈R^3×3是一个右奇异矢量矩阵，于是可以计算出：

S＝diag(11|VU^T|)∈R^3×3 (17)；

式(17)中，S是一个用于在点集中含噪声时避免微小图像匹配的矩阵。

运用式(16)，可以从式(13)中求得旋转矩阵R为：

R＝VSU^T (18)；

本实施例的一个具体应用为：

本文算法在Visual studio2010环境下，Intel Core i73.33GHz的CPU、16GB内存的Windows764位PC机上进行了实现。算法应用到了公共点云数据模型bunny和dragon的配准。实验所用数据均为采用三维激光扫描仪获取的点云数据模型或由三角网格数据模型转换的点云数据模型，模型均满足简单的三维刚体变换。

待配准公共点云数据模型bunny点云和dragon点云分别如图3和图4所示，两组点云的覆盖率均低于60％。采用本文所提配准算法，首先提取点云的凹凸特征点，然后对特征点进行区域分割，最后对区域进行两两配准，从而实现点云模型的配准。公共点云的最终配准结果如图5所示，显然，本文所提配准算法可以实现公共点云的准确配准。

为了进一步验证本文所提配准算法在配准精度和耗时等方面的性能，针对图3和图4的两组公共点云，再分别采用ICP算法、传统算法1和传统算法2进行配准，四种算法的配准结果如下表所示：

从上表可见，本文算法最有最高的配准精度和最快的配准速度。与ICP算法相比，本文算法的配准精度和速度分别提高了约35％和50％；与传统算法1相比，本文算法的配准精度和速度分别提高了约25％和38％；与传统算法2相比，本文算法的配准精度和速度分别提高了23％和35％。这是由于ICP算法对点云的初始相对位置要求较高，且要求待配准点云之间存在包含关系；传统算法1是一种基于迭代因子的点云配准算法，该迭代因子可以提高点云配准的迭代速度，但是不能提高低覆盖率点云的配准精度；传统算法2算法是一种基于ISS特征点的点云配准算法，该算法通过对点云上特征点的特征序列的匹配来实现点云配准，取得了较精确的匹配结果，但是对覆盖率低于60％的点云配准效果不佳。而本文配准算法通过区域分割有效避免了低覆盖区域对点云配准的影响，可以实现点云的快速精确配准。

值得注意的是，上述系统实施例中，所包括的各个单元只是按照功能逻辑进行划分的，但并不局限于上述的划分，只要能够实现相应的功能即可；另外，各功能单元的具体名称也只是为了便于相互区分，并不用于限制本发明的保护范围。

另外，本领域普通技术人员可以理解实现上述各实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，相应的程序可以存储于一计算机可读取存储介质中。

以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。

Claims (3)

1.一种基于区域分割的低覆盖点云配准算法，其特征在于，包括如下步骤：

步骤1：提取凹、凸特征点：采用体积积分不变量判断其点的凹凸性，并提取凹的和凸的特征点集；

步骤2：特征点区域分割：对特征点构建邻接矩阵实现特征点向无向途中节点的转换，并将点的对应关系转换成无向图中的连线；

步骤3：特征点分析：利用MPPCA分析法提取特征点之间的几何关系，并构造邻接矩阵；

步骤4：区域配准：通过SVD的ICP算法对分割区域进行配准；点云配准算法的步骤如下：

步骤41：构建测量模式矩阵并通过平移将点集的重心对其；

步骤42：将SVD用于测量模式矩阵，即并计算旋转矩阵

步骤43：创建旋转模式并根据的最近点构建

步骤44：计算配准误差若ε^(j)小于给定阈值则说明两个点云区域{A_p,B_q}配准成功，否则区域对{A_p,B_q}为不可配准的区域；

步骤45：通过穷举法将区域集合{A₁,A₂,...,A_s}和{B₁,B₂,...,B_t}中区域进行两两配准即可实现初始点云M和点云D的最终配准。

2.根据权利要求1所述的一种基于区域分割的低覆盖点云配准算法，其特征在于，所述步骤1中，体积积分不变量V_r(p)的计算公式为：

V_r(p)＝∫_p+rQg(x)w(p-x)dx (1)；

式(1)中，g(x)表示为特征函数，w(p-x)表示权重函数，p+rQ表示为半径是r、中心是点p的球；式(1)中特征函数换成I_D(x)，并令w(p-x)＝1，则体积积分不变量V_r(p)的计算公式为：

V_r(p)＝∫_p+rQI_D(x)dx (2)；

式(2)中，当x在曲面外侧，I_D(x)＝1；当x在曲面内侧，则I_D(x)＝0。

3.根据权利要求1所述的一种基于区域分割的低覆盖点云配准算法，其特征在于，所述步骤3中，利用EM算法得到最优描述点云分布的M个线性主成分分析器和每个点x_i对应于M个线性主成分分析器的概率密度分布，对于第m个线性主成分分析器，任意以个点x_i对应于一个降维后的点y_i，计算式为：

x_i＝V_my_i+μ_m+ε_m (4)；

式(4)中，V_m表示第m个主成分分析器的主子空间，μ_m表示分布于第m个分析器的所有点均值，ε_m表示噪声；

当y_i和ε_m均服从高斯分布，则x_i关于m的条件概率为：

式(5)中，Z＝x_i-μ_m，目标函数F_MPPCA为第m个主成份分析器上任意点x_i的最大似然估计，计算式为：

式(6)中，π_m为各个主成份分析器的混合比例，π_m≥0，