CN103985155A - 基于映射法的散乱点云Delaunay三角剖分曲面重构方法 - Google Patents

Abstract

本发明涉及一种基于映射法的散乱点云Delaunay三角剖分曲面重构方法，属于计算机图形学、虚拟现实技术领域。具体操作步骤为：①获取目标的原始点云数据。②获取原始点云数据中每个点的K阶邻域和单位法向量。③将点云数据进行分片。④参数化已分片的点云到二维平面。⑤在二维平面内对点云进行Delaunay三角剖分并映射回对应的三维空间。⑥优化初始三角网格模型。与已有技术相比较，本发明方法能够在进行大规模点云数据网格建模时，在保证三角网格的质量的同时，能够较快的实现散乱点云的三角网格化，对于海量点云效果更佳。

Description

基于映射法的散乱点云Delaunay三角剖分曲面重构方法

技术领域

本发明涉及一种基于映射法的散乱点云德劳内(Delaunay)三角剖分曲面重构方法，属于计算机图形学、虚拟现实技术领域。

背景技术

在逆向工程中，获取目标的原始点云数据是逆向工程的第一步，点云数据的获取主要是通过三维扫描仪。当前的三维扫描仪主要可以分为两类：普通光学扫描仪(结构光三维扫描仪)和激光三维扫描仪。普通光学扫描仪比较适用于对小型物体进行三维精确建模，具有很高的扫描精度，应用在工业生产中对零件和产品的尺寸进行测量。激光三维扫描仪获得的点云数据除包含数据点坐标外，通常还包含反射强度、法线方向等信息。逆向工程中最为关键的技术在于利用获取的目标原始点云数据构造出被扫描物体的三维模型，即对获得的目标原始点云数据进行处理并生成三维模型。

曲面重构技术是反向工程中的最重要的一步，也是基于点云三维重建中的重要环节。

现有的点云曲面重建算法的可以分为以下几类：参数曲面方法、隐式曲面方法、基于学习的方法、Delaunay三角剖分法等。

参数曲面重建主要应用在CAD领域，典型的重建算法有B-spline曲面、NURBS曲面以及Bezier曲面。基于隐式曲面的方法采用一组隐式曲线方程拟合原始点云，然后在零值面上抽取三角网格。使用最多的为径向基函数(RadialBasis Function，RBF)和移动最小二乘法，典型的代表算法有Hoppe提出的方法。基于学习的曲面重建方法将统计学和机器学习的思想用于散乱点云的曲面重建，输入为通过设备得到的点云数据，输出为逼近的曲面模型，Yumer和Kara使用神经网络重建网格模型。三角网格由其简单实用的曲面表达形式和易于渲染、添加纹理等特性得到了广泛的应用，Delaunay三角剖分的方法用各种不同的方法去实现Voronoi图和Delaunay三角化。代表算法有Amenta等人提出的Crust算法和Bernardini等人提出滚球法。在Crust算法的基础上，Amenta、Choi和Cazals提出了Cocone算法以及Amenta、Choi和Kolluri提出了Power Crust算法。作为常用的曲面重建方法之一的三角剖分技术，既是热点又是难点问题。一般情况下，在众多三角剖分算法中，Delaunay三角剖分构建的网格质量是最优的。

目前，散乱点云的Delaunay三角剖分算法主要包括两类，即空间直接三角剖分算法和投影域三角剖分算法(映射法)。空间直接三角剖分算法是指直接在三维空间内进行三角剖分，直接在三维空间内将散乱点云连接成一个优化的空间三角网格，这类算法虽然能够得到较好的曲面重建效果，但在计算空间拓扑结构和连接三角形计算上较为复杂，效率低下，其算法本身存在很多局限性。投影域三角剖分算法是先将空间的点云数据投影到参数域，在参数域内应用Delaunay三角网格生成算法对其进行网格化，最后将计算结果返回到三维空间区域，由于参数域(如平面区域)的三角网格化技术比较成熟且时间复杂度较低，该方法能够处理大规模的点云数据。但当参数化后的点云数据有重叠区域时，需要对点云数据进行分块和连接操作，增加了算法的复杂程度，由于在这种方法中对三维点云进行了降维操作，会在一定程度导致三维空间内的距离信息的丢失。

使用映射法进行散乱点云的Delaunay三角网格化，点云分片和点云参数化是关键的部分。点云分片的一个重要原则就是每一个分片点云进行映射后能够不产生重叠点云数据，此外，点云分片还要求能够产生较少的分片，减少后期拼接工作量。从目前国内外研究及应用来看，尽管对点云分片已经进行了大量的、面向不同应用问题的研究，但是还没有一种适合所有应用的分片算法。绝大多数算法都是针对具体问题提出的。目前的分片算法往往需要大量的人工干预，降低了三角剖分算法的效率。将三维点云映射到二维平面是无网格参数化的一种，目前大多数研究都是针对于网格参数化的，无网格参数化的研究工作较少。研究的目标为能使映射到二维平面的点云并最大程度地保持三维空间的位置距离信息。Floater将重心坐标映射的三角网格的平面参数化方法应用到点云，提出了一种基于点的邻域的凸线性组合的无网格参数化方法。Tewari等利用离散l-形式参数化方法建立亏格为l的点云的周期参数化。2010年，Zhang Lei提出一种能够较好保持三维点云拓扑信息的方法，采用从局部到整体的思想。该方法在进行整体的刚性变换时计算复杂。Er Li等采用基于KNN图论的方法得到方向场引导参数化过程。由于在二维平面内进行Delaunay三角剖分由于涉及到点云的分片和拼接操作，往往需要的人工干预才能够完成点云数据的重构。Bradley对多视点的图像进行双目立体匹配得到点云，然后根据其特有特征在低维进行网格化。Zhang Lei将其无网格参数化方法应用于点云的重建，但没有对大规模点云分片和拼接作介绍。

本发明涉及到的已有技术为：

1.八叉树(Octree)结构：是平面的四叉树结构向三位空间的拓展，通过对点云所在的三维空间进行递归的划分而得到。建立八叉树时要首先建立一个包含所有点云数据的最小的立方体包围盒，把它作为八叉树结构的根节点，代表整个空间；再将该立方体包围盒分割成大小相同的8个子立方体，每个子立方体均视为根节点的字节点；如此递归分割，直至最小子立方体的边长等于给定的点距。

2.K阶邻域：是指点云数据中距离某一点最近的K个点。

3.Delaunay三角网：是指具有空圆特性和最大化最小角特性的三角网。空圆特性是指，在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在。最大化最小角特性是指在散点集可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小内角最大。

发明内容

本发明的目的是为了克服上述已有技术的缺陷，解决大规模点云数据重建效率较低、人工干预度高等问题，提出一种基于映射法的散乱点云Delaunay三角剖分曲面重构方法。

本发明方法是通过下述技术方案实现的。

本发明的一种基于映射法的散乱点云Delaunay三角剖分曲面重构方法，其具体步骤如下：

步骤一、获取目标的原始点云数据。

设定点云数据坐标系：建立空间直角坐标系，以平向右方向为X轴正方向，竖直向上方向为Z轴正方向，垂直于X轴和Z轴所确定的平面的轴为Y轴。它们的正方向符合右手规则。

所述原始点云数据包含三维坐标信息，还可能包含每个点的法向量信息。

步骤二、获取原始点云数据中每个点的K阶邻域和单位法向量。

在步骤一操作的基础上，获取原始点云数据中每个点的K阶邻域和单位法向量，具体操作步骤为：

步骤2.1：获取原始点云数据中每个点的K阶邻域，具体操作步骤为：

步骤2.1.1：采用八叉树结构处理步骤一得到的点云数据，得到八叉树点云数据；

步骤2.1.2：对于步骤2.1.1得到的八叉树点云数据中的每一点(用符号P表示)，通过搜索点P所在的叶子节点和与该节点相邻的叶子节点中最近的k个点得到点P的K阶邻域。其中，k为正整数，4≤k≤100。

步骤2.2：在步骤2.1操作的基础上，获取八叉树点云数据中每个点的单位法向量，具体操作分为2种情况：

情况1：如果八叉树点云数据中每个点包含法向量信息，对每个点的法向量进行单位化得到其单位法向量。

情况2：如果八叉树点云数据中每个点不包含法向量信息，则采用主元分析法(PCA)计算点云数据的单位法向量。其步骤如下：

步骤2.2.1：对于点云中任意一点P，通过公式(1)使用最小二乘法为点P和点P的K阶邻域计算出一个局部平面(用符号S表示)。

$S(n,d)={\arg \min }_{(n,d)}\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{(n\·P_i-d)}^2---\left(1\right)$

其中，n为平面S的单位法向量；d为点P到坐标原点的距离，P_i为点P的k个邻近点。

步骤2.2.2：找到点P以及点P的k个邻近点的质心(用符号P＇表示)，并通过公式(2)获得半正定协方差矩阵(用符号M表示)；然后对半正定协方差矩阵M进行特征值分解，将半正定协方差矩阵M的最小特征值对应的特征向量作为点P的法向量。

$M=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}(\mathrm{Pi}-P^{\′}){(\mathrm{Pi}-P^{\′})}^T/k---\left(2\right)$

步骤2.2.3：对步骤2.2.2得到的点P的法向量进行单位化得到点P的单位法向量。

步骤2.2.4：由于步骤2.2.3计算出的单位法向量方向与真实的单位法向量可能相反，因此需要判断单位法向量方向是否需要进行调整，如果单位法向量方向与真实的单位法向量相反则进行调整。

判断单位法向量方向是否需要进行调整以及调整方法为：

在八叉树点云数据的所有点中，找到其Z坐标值最大的点(用符号A表示)，然后选取点A的邻域中任意一点(用符号N表示)，计算向量NA与点A的单位法向量(用符号na表示)之间的夹角(用符号β(NA,na)表示)，如果β(NA,na)>π/2，则将点A的单位法向量变换方向。以点A为基准点调整其k个邻域点的单位法向量(用n₁表示)方向，计算单位法向量na与单位法向量n₁之间的夹角(用符号β(na,n₁)表示)，如果β(na,n₁)>π/2，则改变该邻域点的单位法向量方向；此后再依次以点A的邻域点为基准点，调整点A的邻域点的邻域点的单位法向量，重复此过程，直到所有单位法向量都调整完毕。

步骤三、将点云数据进行分片。

在步骤二计算点云数据单位法向量的基础上，对点云数据进行分片。具体操作步骤如下：

步骤3.1：搜索种子点。

所述搜索种子点的方法为：在步骤二得到的点云数据中，从未分片的点中随机选取一个点(用符号O＇表示)，将点O＇的k个邻域点的单位法向量投影到点O＇与点O＇的单位法向量确定的平面上，形成k条射线，对每一条射线(用符号L表示)，以射线L的起点为起点，经过点O＇的形成一个新的射线，如果该射线与射线L的夹角全部大于π或者全部都小于π，则点O＇为种子点，否则，点O＇不是种子点。如果点O＇不是种子点，重复该步骤，直到找到一个种子点(用符号O表示)。

步骤3.2：将步骤3.1得到的种子点O加入到待拓展队列(用符号Queue表示)中，将该种子点为点云分片(用符号s表示)的第一个点，同时将种子点标记为已分片；

步骤3.3：若待拓展队列Queue不为空，从队列中取出一个点(用符号Q表示)，遍历点Q的k个邻域点，判断点Q的当前邻域点的法向量与种子点O的法向量夹角(用符号α表示)是否低于一个阈值(用符号θ表示)，其中，θ为人为设定值，如果α≤θ成立，则将该邻域点加入待拓展队列Queue和当前的点云分片s中，同时标记该邻域点为已分片。重复此过程，直到待拓展队列Queue为空时，点云的一个分片的分割完成。

步骤3.4：重复执行步骤3.1至步骤3.3，直到所有的种子点都标记为已分片。

步骤3.5：经过步骤3.1至步骤3.4的操作，如果存在少数点没有被划分到任意一个分片当中，或者分片中的点的数量少于阈值(用符号σ表示)，其中，σ为人为设定值，5≤σ≤100。对于这样的点，将其划分到含有最多该点的领域点的分片中。

步骤3.6：对每一个分片的边缘进行拓展。其具体操作方法为：遍历分片中全部点，查看每个点的邻域点的所属分片，若邻域点的分片和该点的分片不同，则将该点加入到邻域点所属的分片中。

经过步骤三的操作，即可完成全部点云数据的分片。

步骤四、参数化已分片的点云到二维平面。

在步骤三对点云数据分片的基础上，进行分片点云的无网格参数化，将空间的三维点云影射到二维平面空间。其具体操作步骤如下：

步骤4.1：从步骤三得到的分片中，取一个未参数化的分片，将该分片中的种子点O加入待参数化队列(用符号Queue2表示)。

步骤4.2：若待参数化队列Queue2不为空，从待参数化队列Queue2中取出一个点(用符号B表示)，将点B标记为已调整，将点B直接投影射到以点O和点O的法向量确定的平面上，得到投影点(用符号B′表示)。对于点B的k个邻域点中的每一点(用符号D表示)，将点D直接投影射到以点O和点O的法向量确定的平面上，得到投影点(用符号D′表示)。沿着B′D′的方向将B′D′延伸，得到一条长度为|BD|的线段，该线段记为B′D″，点B′和点D″为线段B′D″的两个端点。点D″为点D经过一次调整后得到的映射点。将点D″加入到点D的调整点集合(用符号coSet表示)中。如果点D为未调整，则将点D加入待参数化队列Queue2中。

步骤4.3：重复步骤4.2直到待参数化队列Queue2为空。

步骤4.4：依次遍历步骤4.1中所述分片中的点，对每一点计算其调整点集合coSet中所有点的坐标平均值，将此均值坐标作为参数化结果。

经过步骤4.1至步骤4.4的操作，完成了一个分片的无网格参数化。

步骤4.5：重复执行步骤4.1至步骤4.4，直到所有的分片都完成无网格参数化。

步骤五、在二维平面内对点云进行Delaunay三角剖分并映射回对应的三维空间。

在步骤四操作的基础上，对全部无网格参数化结果的点云分片进行德劳内(Delaunay)三角网的构建。具体操作步骤为：

步骤5.1：采用德劳内三角网格的方法重构无网格参数化结果的点云分片，德劳内三角网所具有的空圆特性和最大最小角特性，保证了德劳内三角网中不会出现过于狭长的三角形，使得三角网的构建更合理与准确，从而具有极大的应用价值。

步骤5.2：经过步骤5.1对所有的点云分片进行德劳内三角网构建完成后，将全部分片的构网结果返回到三维空间，得到原始点云数据的初始三角网格模型。

步骤六、优化初始三角网格模型。

在步骤五初始三角网格模型的基础上，进行网格重叠和空洞的优化。具体操作步骤如下：

步骤6.1：统计步骤五得到的初始三角网格模型中，每条边使用的次数(用符号count表示)。若某一条边使用次数count＝3，则表示共同使用该边(用符号EF表示)的有3个三角形，其中点E和点F为边EF的两个端点。计算使用此边的三角形所在平面之间的夹角，选择其中夹角最小的两个三角形，将这两个三角形分别称为第一三角形和第二三角形，将另外一个三角形称为第三三角形。第一三角形、第二三角形和第三三角形的另外一个顶点分别用符号G、H和I表示。选取边EF的中点(用符号J表示)，如果∠GJI>∠HJI，则删除第二三角形；否则，删除第一三角形。重复该步骤，直到每条边的使用次数count均不大于2。

步骤6.2：统计步骤6.1得到的三角网格模型中，每条边使用的次数count。若某一条边使用次数count＝1，将此边称为边界边(用符号RT表示)，其中点R和点T为边RT的两个端点。将使用边界边RT的三角形的另外一个顶点用符号V表示。将点R和点T称为孔洞点。搜索点R和点T的邻域点中的孔洞点，并组成点集(用符号SET表示)。选取边RT的中点(用符号U表示)。如果点集SET不为空，则遍历点集SET中的点(用符号W表示)，选取∠VUW中的最大值对应的点(用符号W_max表示)。将点W_max、R和T组成一个新的三角形加入到三角网络模型中，同时更新ΔW_maxRT的三条边的使用次数count。如果点集SET为空，将边界边RT的使用次数count值更新为2。重复该步骤，直到每条边的使用次数count均等于2。

经过上述步骤的操作，得到原始点云数据的最终三角网格模型。

有益效果

与已有技术相比较，本发明方法能够在进行大规模点云数据网格建模时，在保证三角网格的质量的同时，能够较快的实现散乱点云的三角网格化，对于海量点云效果更佳。

附图说明

图1为本发明具体实施方式中基于映射法的散乱点云Delaunay三角剖分曲面重构方法的整体框架设计流程图；

图2为本发明具体实施方式中模具的原始点云数据示意图；

图3为本发明具体实施方式中对图2中的点云数据进行分片的结果示意图；

图4为本发明具体实施方式中图2中的原始点云数据的初始三角网格模型示意图；

图5为本发明具体实施方式中图2中的原始点云数据的最终三角网格模型示意图；

图6为本发明具体实施方式中兔子的原始点云数据示意图；

图7为本发明具体实施方式中猫的原始点云数据示意图；

图8为本发明具体实施方式中汽车配件的原始点云数据示意图；

图9为本发明具体实施方式中人脸的原始点云数据示意图。

具体实施方式

根据上述技术方案，下面结合附图和实施实例对本发明进行详细说明。

本实施例中，原始点云数据如图2所示，对图2由50548个点组成，每一个点只包含其三维坐标数值。对如图2所示的点云数据采用本发明提出的基于映射法的散乱点云Delaunay三角剖分曲面重构方法进行三角曲面重构，其整体框架设计流程如图1所示，其具体步骤如下：

步骤一、获取目标的原始点云数据。

设定点云数据坐标系：建立空间直角坐标系，以平向右方向为X轴正方向，竖直向上方向为Z轴正方向，垂直于X轴和Z轴所确定的平面的轴为Y轴。它们的正方向符合右手规则。

所述原始点云数据只包含三维坐标信息。

步骤二、获取原始点云数据中每个点的K阶邻域和单位法向量。

在步骤一操作的基础上，获取原始点云数据中每个点的K阶邻域和单位法向量，具体操作步骤为：

步骤2.1：获取原始点云数据中每个点的K阶邻域，具体操作步骤为：

步骤2.1.1：采用八叉树结构处理步骤一得到的点云数据，得到八叉树点云数据；

步骤2.1.2：对于步骤2.1.1得到的八叉树点云数据中的每一点P，通过搜索点P所在的叶子节点和与该节点相邻的叶子节点中最近的k个点得到点P的K阶邻域。其中，k＝10。

步骤2.2：在步骤2.1操作的基础上，获取八叉树点云数据中每个点的单位法向量，本点云数据中每个点不包含法向量信息，采用主元分析法(PCA)计算点云数据的单位法向量。其步骤如下：

步骤2.2.1：对于点云中任意一点P，通过公式(1)使用最小二乘法为点P和点P的K阶邻域计算出一个局部平面S。

步骤2.2.2：找到点P以及点P的10个邻近点的质心P＇，并通过公式(2)获得半正定协方差矩阵M；然后对半正定协方差矩阵M进行特征值分解，将半正定协方差矩阵M的最小特征值对应的特征向量作为点P的法向量。

步骤2.2.3：对步骤2.2.2得到的点P的法向量进行单位化得到点P的单位法向量。将点P的法向量与点P法向量模长的比值作为点P的单位法向量。

步骤2.2.4：由于步骤2.2.3计算出的单位法向量方向与真实的单位法向量可能相反，因此需要判断单位法向量方向是否需要进行调整，如果单位法向量方向与真实的单位法向量相反则进行调整。

判断单位法向量方向是否需要进行调整以及调整方法为：

在八叉树点云数据的所有点中，找到其Z坐标值最大的点A，然后选取点A的邻域中任意一点N，计算向量NA与点A的单位法向量na之间的夹角β(NA,na)，如果β(NA,na)>π/2，则将点A的单位法向量变换方向。以点A为基准点调整其10个邻域点的单位法向量n₁方向，计算单位法向量na与单位法向量n₁之间的夹角β(na,n₁)，如果β(na,n₁)>π/2，则改变该邻域点的单位法向量方向；此后再依次以点A的邻域点为基准点，调整点A的邻域点的邻域点的单位法向量，重复此过程，直到所有单位法向量都调整完毕。

步骤三、将点云数据进行分片。

在步骤二计算点云数据单位法向量的基础上，对点云数据进行分片。具体操作步骤如下：

步骤3.1：搜索种子点。

所述搜索种子点的方法为：在步骤二得到的点云数据中，从未分片的点中随机选取一个点O＇，将点O＇的10个邻域点的单位法向量投影到点O＇与点O＇的单位法向量确定的平面上，形成10条射线，对每一条射线L，以射线L的起点为起点，经过点O＇的形成一个新的射线，如果该射线与射线L的夹角全部大于π或者全部都小于π，则点O＇为种子点，否则，点O＇不是种子点。如果点O＇不是种子点，重复该步骤，直到找到一个种子点O。

步骤3.2：将步骤3.1得到的种子点O加入到待拓展队列Queue中，将该种子点为点云分片s的第一个点，同时将种子点标记为已分片；

步骤3.3：若待拓展队列Queue不为空，从队列中取出一个点Q，遍历点Q的k个邻域点，判断点Q的当前邻域点的法向量与种子点O的法向量夹角α是否低于3/2π。如果α≤3/2π成立，则将该邻域点加入待拓展队列Queue和当前的点云分片s中，同时标记该邻域点为已分片。重复此过程，直到待拓展队列Queue为空时，点云的一个分片的分割完成。

步骤3.4：重复执行步骤3.1至步骤3.3，直到所有的种子点都标记为已分片。

步骤3.5：经过步骤3.1至步骤3.4的操作，如果存在少数点没有被划分到任意一个分片当中，或者分片中的点的数量少于20。对于这样的点，将其划分到含有最多该点的领域点的分片中。

步骤3.6：对每一个分片的边缘进行拓展。其具体操作方法为：遍历分片中全部点，查看每个点的邻域点的所属分片，若邻域点的分片和该点的分片不同，则将该点加入到邻域点所属的分片中。

经过步骤三的操作，即可完成全部点云数据的分片，如图3所示。

步骤四、参数化已分片的点云到二维平面。

在步骤三对点云数据分片的基础上，进行分片点云的无网格参数化，将空间的三维点云影射到二维平面空间。其具体操作步骤如下：

步骤4.1：从步骤三得到的分片中，取一个未参数化的分片，将该分片中的种子点O加入待参数化队列Queue2。

步骤4.2：若待参数化队列Queue2不为空，从待参数化队列Queue2中取出一个点B，将点B标记为已调整，将点B直接投影射到以点O和点O的法向量确定的平面上，得到投影点B′。对于点B的10个邻域点中的每一点D，将点D直接投影射到以点O和点O的法向量确定的平面上，得到投影点D′。沿着B′D′的方向将B′D′延伸，得到一条长度为|BD|的线段，该线段记为B′D″，点B′和点D″为线段B′D″的两个端点。点D″为点D经过一次调整后得到的映射点。将点D″加入到点D的调整点集合coSet中。如果点D为未调整，则将点D加入待参数化队列Queue2中。

步骤4.3：重复步骤4.2直到待参数化队列Queue2为空。

步骤4.4：依次遍历步骤4.1中所述分片中的点，对每一点计算其调整点集合coSet中所有点的坐标平均值，将此均值坐标作为参数化结果。

经过步骤4.1至步骤4.4的操作，完成了一个分片的无网格参数化。

步骤4.5：重复执行步骤4.1至步骤4.4，直到所有的分片都完成无网格参数化。

步骤五、在二维平面内对点云进行Delaunay三角剖分并映射回对应的三维空间。

在步骤四操作的基础上，对全部无网格参数化结果的点云分片进行德劳内(Delaunay)三角网的构建。具体操作步骤为：

步骤5.1：采用德劳内三角网格的方法重构无网格参数化结果的点云分片，德劳内三角网所具有的空圆特性和最大最小角特性，保证了德劳内三角网中不会出现过于狭长的三角形，使得三角网的构建更合理与准确，从而具有极大的应用价值。

步骤5.2：经过步骤5.1对所有的点云分片进行德劳内三角网构建完成后，将全部分片的构网结果返回到三维空间，得到原始点云数据的初始三角网格模型，如图4所示。

步骤六、优化初始三角网格模型。

在步骤五初始三角网格模型的基础上，进行网格重叠和空洞的优化。具体操作步骤如下：

步骤6.1：统计步骤五得到的初始三角网格模型中，每条边使用的次数count。若某一条边使用次数count＝3，则表示共同使用该边EF的有3个三角形，其中点E和点F为边EF的两个端点。计算使用此边的三角形所在平面之间的夹角，选择其中夹角最小的两个三角形，将这两个三角形分别称为第一三角形和第二三角形，将另外一个三角形称为第三三角形。第一三角形、第二三角形和第三三角形的另外一个顶点分别用符号G、H和I表示。选取边EF的中点J，如果∠GJI>∠HJI，则删除第二三角形；否则，删除第一三角形。重复该步骤，直到每条边的使用次数count均不大于2。

步骤6.2：统计步骤6.1得到的三角网格模型中，每条边使用的次数count。若某一条边使用次数count＝1，将此边称为边界边RT，其中点R和点T为边RT的两个端点。将使用边界边RT的三角形的另外一个顶点用符号V表示。将点R和点T称为孔洞点。搜索点R和点T的邻域点中的孔洞点，并组成点集SET。选取边RT的中点U。如果点集SET不为空，则遍历点集SET中的点W，选取∠VUW中的最大值对应的点W_max。将点W_max、R和T组成一个新的三角形加入到三角网络模型中，同时更新ΔW_maxRT的三条边的使用次数count。如果点集SET为空，将边界边RT的使用次数count值更新为2。重复该步骤，直到每条边的使用次数count均等于2。

经过上述步骤的操作，得到原始点云数据的最终三角网格模型，如图5所示，其三角网格的质量较高。

为了说明本发明的效率，采用本发明发明对兔子(bunny)、猫(cat)、汽车配件(car)和人脸(head-half)进行曲面重建，同时采用经典的crust算法和TCocone算法对相同的原始点云数据进行处理，耗时情况如表1所示。

表1三种不同方法的耗时对比情况

通过对比，可以得出本发明方法在保证三角网格的质量的同时，能够较快的实现散乱点云的三角网格化。

虽然结合附图描述了本发明的实施方式，但是对于本领域技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些也应视为属于本发明的保护范围。

Claims (2)

1.一种基于映射法的散乱点云Delaunay三角剖分曲面重构方法，其特征在于：其具体步骤如下：

步骤一、获取目标的原始点云数据；

设定点云数据坐标系：建立空间直角坐标系，以平向右方向为X轴正方向，竖直向上方向为Z轴正方向，垂直于X轴和Z轴所确定的平面的轴为Y轴；它们的正方向符合右手规则；

所述原始点云数据包含三维坐标信息，还可能包含每个点的法向量信息；

步骤二、获取原始点云数据中每个点的K阶邻域和单位法向量；

在步骤一操作的基础上，获取原始点云数据中每个点的K阶邻域和单位法向量，具体操作步骤为：

步骤2.1：获取原始点云数据中每个点的K阶邻域；

步骤2.2：在步骤2.1操作的基础上，获取八叉树点云数据中每个点的单位法向量，具体操作分为2种情况：

情况1：如果八叉树点云数据中每个点包含法向量信息，对每个点的法向量进行单位化得到其单位法向量；

情况2：如果八叉树点云数据中每个点不包含法向量信息，则采用主元分析法计算点云数据的单位法向量；其步骤如下：

步骤2.2.1：对于点云中任意一点P，通过公式(1)使用最小二乘法为点P和点P的K阶邻域计算出一个局部平面S；

$S(n,d)={\arg \min }_{(n,d)}\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{(n\·P_i-d)}^2---\left(1\right)$

其中，n为平面S的单位法向量；d为点P到坐标原点的距离，P_i为点P的k个邻近点，k为正整数，4≤k≤100；

步骤2.2.2：找到点P以及点P的k个邻近点的质心P＇，并通过公式(2)获得半正定协方差矩阵M；然后对半正定协方差矩阵M进行特征值分解，将半正定协方差矩阵M的最小特征值对应的特征向量作为点P的法向量；

$M=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}(\mathrm{Pi}-P^{\′}){(\mathrm{Pi}-P^{\′})}^T/k---\left(2\right)$

步骤2.2.3：对步骤2.2.2得到的点P的法向量进行单位化得到点P的单位法向量；

步骤2.2.4：由于步骤2.2.3计算出的单位法向量方向与真实的单位法向量可能相反，因此需要判断单位法向量方向是否需要进行调整，如果单位法向量方向与真实的单位法向量相反则进行调整；

判断单位法向量方向是否需要进行调整以及调整方法为：

在八叉树点云数据的所有点中，找到其Z坐标值最大的点A，然后选取点A的邻域中任意一点N，计算向量NA与点A的单位法向量na之间的夹角β(NA,na)，如果β(NA,na)>π/2，则将点A的单位法向量变换方向；以点A为基准点调整其k个邻域点的单位法向量n₁方向，计算单位法向量na与单位法向量n₁之间的夹角β(na,n₁)，如果β(na,n₁)>π/2，则改变该邻域点的单位法向量方向；此后再依次以点A的邻域点为基准点，调整点A的邻域点的邻域点的单位法向量，重复此过程，直到所有单位法向量都调整完毕；

步骤三、将点云数据进行分片；

在步骤二计算点云数据单位法向量的基础上，对点云数据进行分片；具体操作步骤如下：

步骤3.1：搜索种子点；

所述搜索种子点的方法为：在步骤二得到的点云数据中，从未分片的点中随机选取一个点O＇，将点O＇的k个邻域点的单位法向量投影到点O＇与点O＇的单位法向量确定的平面上，形成k条射线，对每一条射线L，以射线L的起点为起点，经过点O＇的形成一个新的射线，如果该射线与射线L的夹角全部大于π或者全部都小于π，则点O＇为种子点，否则，点O＇不是种子点；如果点O＇不是种子点，重复该步骤，直到找到一个种子点O；

步骤3.2：将步骤3.1得到的种子点O加入到待拓展队列Queue中，将该种子点为点云分片s的第一个点，同时将种子点标记为已分片；

步骤3.3：若待拓展队列Queue不为空，从队列中取出一个点Q，遍历点Q的k个邻域点，判断点Q的当前邻域点的法向量与种子点O的法向量夹角α是否低于一个阈值θ，其中，θ为人为设定值，如果α≤θ成立，则将该邻域点加入待拓展队列Queue和当前的点云分片s中，同时标记该邻域点为已分片；重复此过程，直到待拓展队列Queue为空时，点云的一个分片的分割完成；

步骤3.4：重复执行步骤3.1至步骤3.3，直到所有的种子点都标记为已分片；

步骤3.5：经过步骤3.1至步骤3.4的操作，如果存在少数点没有被划分到任意一个分片当中，或者分片中的点的数量少于阈值σ，其中，σ为人为设定值，5≤σ≤100；对于这样的点，将其划分到含有最多该点的领域点的分片中；

步骤3.6：对每一个分片的边缘进行拓展；其具体操作方法为：遍历分片中全部点，查看每个点的邻域点的所属分片，若邻域点的分片和该点的分片不同，则将该点加入到邻域点所属的分片中；

经过步骤三的操作，即可完成全部点云数据的分片；

步骤四、参数化已分片的点云到二维平面；

在步骤三对点云数据分片的基础上，进行分片点云的无网格参数化，将空间的三维点云影射到二维平面空间；其具体操作步骤如下：

步骤4.1：从步骤三得到的分片中，取一个未参数化的分片，将该分片中的种子点O加入待参数化队列Queue2；

步骤4.2：若待参数化队列Queue2不为空，从待参数化队列Queue2中取出一个点B，将点B标记为已调整，将点B直接投影射到以点O和点O的法向量确定的平面上，得到投影点B′；对于点B的k个邻域点中的每一点D，将点D直接投影射到以点O和点O的法向量确定的平面上，得到投影点D′；沿着B′D′的方向将B′D′延伸，得到一条长度为|BD|的线段，该线段记为B′D″，点B′和点D″为线段B′D″的两个端点；点D″为点D经过一次调整后得到的映射点；将点D″加入到点D的调整点集合coSet中；如果点D为未调整，则将点D加入待参数化队列Queue2中；

步骤4.3：重复步骤4.2直到待参数化队列Queue2为空；

步骤4.4：依次遍历步骤4.1中所述分片中的点，对每一点计算其调整点集合coSet中所有点的坐标平均值，将此均值坐标作为参数化结果；

经过步骤4.1至步骤4.4的操作，完成了一个分片的无网格参数化；

步骤4.5：重复执行步骤4.1至步骤4.4，直到所有的分片都完成无网格参数化；

步骤五、在二维平面内对点云进行Delaunay三角剖分并映射回对应的三维空间；

在步骤四操作的基础上，对全部无网格参数化结果的点云分片进行德劳内三角网的构建；具体操作步骤为：

步骤5.1：采用德劳内三角网格的方法重构无网格参数化结果的点云分片，德劳内三角网所具有的空圆特性和最大最小角特性，保证了德劳内三角网中不会出现过于狭长的三角形，使得三角网的构建更合理与准确，从而具有极大的应用价值；

步骤5.2：经过步骤5.1对所有的点云分片进行德劳内三角网构建完成后，将全部分片的构网结果返回到三维空间，得到原始点云数据的初始三角网格模型；

步骤六、优化初始三角网格模型；

在步骤五初始三角网格模型的基础上，进行网格重叠和空洞的优化；具体操作步骤如下：

步骤6.1：统计步骤五得到的初始三角网格模型中，每条边使用的次数count；若某一条边使用次数count＝3，则表示共同使用该边EF的有3个三角形，其中点E和点F为边EF的两个端点；计算使用此边的三角形所在平面之间的夹角，选择其中夹角最小的两个三角形，将这两个三角形分别称为第一三角形和第二三角形，将另外一个三角形称为第三三角形；第一三角形、第二三角形和第三三角形的另外一个顶点分别用符号G、H和I表示；选取边EF的中点J，如果∠GJI>∠HJI，则删除第二三角形；否则，删除第一三角形；重复该步骤，直到每条边的使用次数count均不大于2；

步骤6.2：统计步骤6.1得到的三角网格模型中，每条边使用的次数count；若某一条边使用次数count＝1，将此边称为边界边RT，其中点R和点T为边RT的两个端点；将使用边界边RT的三角形的另外一个顶点用符号V表示；将点R和点T称为孔洞点；搜索点R和点T的邻域点中的孔洞点，并组成点集SET；选取边RT的中点U；如果点集SET不为空，则遍历点集SET中的点W，选取∠VUW中的最大值对应的点W_max；将点W_max、R和T组成一个新的三角形加入到三角网络模型中，同时更新ΔW_maxRT的三条边的使用次数count；如果点集SET为空，将边界边RT的使用次数count值更新为2；重复该步骤，直到每条边的使用次数count均等于2；

经过上述步骤的操作，得到原始点云数据的最终三角网格模型。

2.如权利要求1所述的一种基于映射法的散乱点云Delaunay三角剖分曲面重构方法，其特征在于：步骤二中所述获取原始点云数据中每个点的K阶邻域的具体操作步骤为：

步骤2.1.1：采用八叉树结构处理步骤一得到的点云数据，得到八叉树点云数据；

步骤2.1.2：对于步骤2.1.1得到的八叉树点云数据中的每一点P，通过搜索点P所在的叶子节点和与该节点相邻的叶子节点中最近的k个点得到点P的K阶邻域；其中，k为正整数，4≤k≤100。