点云数据三角化

1、第1章、1.1空间曲面上散乱数据三角剖分的概念目前有多种方法可以获得物理模型的形状信息。在制造工业中最常使用的是 坐标测量机(Coordinate Measuring Machine,CMM)。坐标测量机能精确测量物体表 面上点的位置，但其测量速度较慢，当测量点数较多时，效率很低，一般用在对 精度要求较高的场合，如检查零件的形状精度、位置精度等。当需要大量获取零 件表面的数据点时，一般使用激光扫描仪(Laser Scanner)。激光扫描仪能在相对 较短的时间内得到大量零件表面的数据点。另外一种在医学上常用的测量设备是 计算机断层扫描仪(Computerized Tomography,CT)。

2、CT得到的是物体的轮廓线， 数据点呈层状分布，每一层代表物体的一个剖面。这些测量设备得到的数据点形式各不相同，虽然在局部上某些数据点具有有 组织的状态，如激光扫描仪和CT所得的数据点呈现层状的特点，但在全局上基 本均表现出散乱的特点。所谓散乱数据的三角剖分就是给定一组散乱数据点，将各数据点之间以三角形相互连接，形成一张三角网格。其实质是以三角网格反映数据点与其邻近 点间的拓扑连接关系。而正确的拓扑连接关系将有效揭示散乱数据集所蕴涵的原 始物体表面的形状和拓扑结构。1.2空间曲面上散乱数据三角剖分的研究意义及应用范围空间曲面上散乱数据的三角剖分是构造散乱数据插值曲面的必不可少的前 置处理步骤，也

3、是最重要最关键的一步，基于散乱数据点三角剖分构造散乱数据 插值曲面的过程如图1所示：图1.1构造散乱数据插值曲面的步骤空间曲面上散乱数据的三角剖分是在对测量数据点必要的处理之后进行的， 它是构造散乱数据插值曲面的前置处理步骤。平面域内的散乱数据的三角剖分研 究己经经历了相当长的时间，相关理论与算法己经相当成熟，特别是Delaunay 三角剖分及其优化准则等研究成果使得平面内的散乱数据点的三角剖分已经不 再困难。但在把平面内的算法推广到空间曲面上时，由于空间曲面散乱数据点之 间拓扑关系的复杂性，对其直接剖分的理论和算法尚不完善。在处理空间曲面上 数据点时，一般的算法都是基于平面凸域的或者是已知各

4、种约束条件的情况，对 于与多值曲面对应的空间曲面上散乱数据点的三角剖分算法研究的较少，且在算 法效率和剖分效果上还远远不能让人满意。散乱数据曲面重构的难点在于如何在数据点集中快速自动得到邻近点间正 确的拓扑连接关系。目前的测量设备能够在短时间内得到数万乃至几十万数据 点，所能体现的曲面形状信息越来越精细和复杂，因此对曲面构造的效率和效果 提出了较高的要求。研究散乱数据特别是大规模散乱数据的三角剖分，对于迅速 构建数据点之间的拓扑连接关系，进而构造插值曲面具有十分重要的意义。散乱数据的三角剖分不但是构造散乱数据插值曲面的重要基础，对三维数据 场的可视化、快速原型制造等新技术的研究也有巨大的推动作

5、用。因此广泛应用 于测量造型、计算机视觉、医学、气像、勘探、环保等领域中。本文所要研究三维散乱数据的直接剖分方法旨在探索解决与复杂曲面对应 的散乱数据点的三角化方法，快速准确地生成优化的三角网格，为构造散乱数据 插值曲面做好准备。因此本文关于空间曲面上散乱数据三角剖分方法的研究对散 乱数据插值曲面的构造以及逆向工程曲面重构方法的研究有着重要的现实意义， 对三维数据场的可视化、基于CT图像的体数据的三维模型重建等应用研究也有 一定的借鉴与参考价值。第2章22.1引言空间曲面上散乱数据的三角剖分是科学计算与分析中的一种重要方法，在计 算几何散乱数据插值曲面构造和三维数据场可视化中得到了广泛应用。对

6、空间曲 面上散乱数据的三角剖分方法的研究不论是在二维平面区域还是在三维空间区 域上都己经有了很多成果，尤其是对二维平面散乱数据三角剖分的研究，其理论 和算法已经比较成熟，相对来说对空间曲面乱数据的三角剖分，特别是曲面形状 比较复杂、散乱数据点的数目很大时，目前的算法在适应性、执行效率等方面还 有待进一步提高。2.2问题描述及分类问题2.1:给定一组散乱数据点Vi(i=1,2,?，n)，如何将各数据点之间以三角 形相互连接，形成一张三角网格，并使网格质量较优。该问题的解是散乱点集Vi 的一个三角剖分T，其实质是以三角网格M反映各个数据点与其邻近点之间的 拓扑连接关系，从而揭示数据点之间的内在本质

7、联系。三角剖分所涉及的问题在 实际应用中根据数据点位置的不同有三种情况:二维，三维实体和空间曲面。根 据这种情况，空间曲面上散乱数据的三角剖分可分为对空间曲面上散乱数据投影 域的剖分和在空间中直接剖分两种类型。空间曲面上散乱数据的投影域包括平面 域和球面域。直接三角剖分方法研究如何直接将空间曲面上散乱数据点在空间中 连接成一个优化的三角网格。2.3三角剖分2.3.1定义定义2.1：给定Ed中k个不相同的点P , P , PP点集 123kP= a p + a p + a p +a p ( a gr, a 河，a + a + a + a =1) (3.1)112233kk k1123k是由p ,

8、 p2，p生成的凸集。定义2.2：给定一个点集的任意子集L, L的凸壳是包含L的最小的子集。定义2.3:给定平面内的顶点的集合Vi(i=1,2,n),用不相交的直线段连接vi和vj, 1Wi,jWn,i勺.使得n个点的凸壳内每一个区域是一个三角形。图2.2平面点集的一种三角剖分一般来说，对顶点集合Vi的三角剖分不是唯一的，有多种剖分结果都能满 足上述定义。当然，其中只有部分结果的三角剖分网格的形态较优，能够满足际 应用的需求。在许多应用中，最好是三角形尽可能是“等边”的，或边总长度最 小，当三角剖分边的总长度减至最小时，则称为最小权三角剖分。2.3.2三角剖分基本理论Voronoi 图与 De

9、launay 三角剖分平面三角剖分的实质是以三角形反映数据点与其邻近点之间的拓扑连接关 系。若能首先找出平面上一点的所有邻近点，则问题2.1在平面上的情况就好办 多了，为此需要解决下面的问题：问题2:给定平面中n个点的集合S=P1,P2,P3,P4对于平面中比其它点更 接近于P的点(x,y)轨迹是什么？图2.3定义2.4:给定平面中n个点的集合s,V(pi)在平面S中比其它点更接近于的 Pi点的轨迹是n-1个半平面的交，它是一个不多于n-1条边的凸多边形域，称 V(pi)为关联于P的Voronoi多边形或关联于P的Voronoi域。图2.4表示关联于 pi的一个Voronoi多边形，它是一个五

10、边形，n=6。图2.4对于S中的每一个点都可以作一个Voronoi多边形，这样的n个Voronoi多 边形组成的图成为Voronoi图，记为Vor(S)，如图2.5所示。该图中的顶点分别 称为Voronoi顶点和Voronoi边。Vor(S)的边是某点对的垂直平分线上的一条线段 或者半直线，从而为该点对所在的两个多边形域所共有。Vor(S)中有的多边形域 是无界的。图2.5在叙述Voronoi图的性质时作假设:原来集合S中没有四个点是共圆的。若此 假设不成立，则需要在证明中加入一些无关紧要的，但却是冗长的陈述。在叙述Voronoi图的性质时作假设“原来集合S中没有四个点是共圆的。若 此假设不成

11、立，则需要在证明中加入一些无关紧要的，但却是冗长的陈述。图2.6图2.6定理2.1:对于S的Voronoi图每一个顶点v，圆C(v)不包含S的其它点。由 于定理2.1： Vor(S)又称为最近点意义下的Voronoi图。定理2.2：在S中，Pi的每一个最近点确定Voronoi多边形v(pi)的一条边。忙甲什缈*P.图2.7图2.7 P的每一个最近邻近点确定v(pi)定义2.5:Voronoi图的直线对偶是由S的每个点对之间加入一个直线段以获 得嵌入平面的图，产生的图是原来n个点上的一个图。如图2.5中的虚线所示。对偶的重要性主要归于下面的Delaunay定理。定理2.3:Voronoi图的直线

12、对偶是S的一个三角剖分。图2.5中的虚线即是S 的一个三角剖分。Voronoi图的这些性质可以用来快速构造Voronoi图，且可应用它解决最邻 近点问题和三角剖分问题。定义2.6:在一个三角剖分中，若每一个三角形的外接圆为空(即在外接圆中 不含有其它点)，则该三角剖分称为Delaunay三角剖分。对于给定的一组平面数据点集S，可以多种不同的三角剖分，其中Delaunay 三角剖分是最优的，二维Delaunay三角剖分由满足最小内角最大准则的三角形 组成。下一节将介绍详细三角剖分准则。2.3.2.2三角剖分优化准则在三角剖分过程中，往往用一种比较简单的方法构造散乱数据点的初始三角 网格，然后对其

13、进行优化，以获得较优化的三角形网格。优化的方法和效果取决 于所采用的优化准则。平面三角剖分常采用的优化准则有Thiessen区域准则、最 小内角最大准则、圆准则。Sibson证明了这三个准则的等价性，并指出符合这三个准则的三角剖分只有 一个，艮口 Delaunay 三角化。今 Thiessen 区域准则(Thiessen regioncriterion)Thiessen 区域指Voronoi分割后形成的多边形区域。如果两个区域具有非零长度的公共线 段，则称这两个区域的生成点为Thiessen强邻接点(strongThiessen neighbours);如 果它们的公共部分仅是一个点，则称为T

14、hiessen弱邻接点(weak Thiessen neighbours)o 一个严格凸的四边形至多有一对相对顶点是Thiessen强邻接点。 Thiessen区域准则指对一个严格凸的四边形三角化时，将Thiessen强邻接点相连， 若两对顶点都是Thiessen弱邻接点，则任选一对相连，这样构造的三角化是最优 的，见图2.8:图2.8 Thiessen区域准则最小内角最大准则在平面内，对一个严格凸的四边形进行三角化时，有两种选择，最小内角最大准则就是要保证对角线两侧两个三角形中的最小内角最大。如图2.9所示：图2.9最小内角最大准则圆准则严格凸的四边形中的三个点确定一个圆，如果第四个顶点落在

15、圆内，则将第 四个顶点与其相对的顶点相连，否则将另两个相对顶点相连，如图210所示：图2.10圆准则平面三角剖分优化准则理论为平面内散乱点集的快速三角剖分提供了判断 标准，也为三维空间散乱点集的三角剖分优化标准提供了借鉴依据。2.3.3经典三角化算法三角化算法虽然很多，但大多算法生成的是Delaunay三角网格，即以空外 接圆(球)准则为优化准则。根据实现方法的不同，Delaunay三角化方法。主要有三类：换边法(Swapping edges):首先构造非优化的初始三角化，然后对2个共边三 角形形成的凸四边形迭代换边优化。以Lawson为代表的对角线交换算法属于换 边法，换边法适用于二维Del

16、aunay三角化，对于三维Delaunay三角化，则需要 对共面四面体进行换面优化。加点法(Adding points):从一个三角形开始，每次加一个点，保证每一步得到 的当前三角化是局部优化的。以英国Bath大学数学分校Bowyer,Green,Sibson为 代表的计算Dirichlet图的方法属于加点法，是较早成名的算法之一;以澳大利亚 悉尼大学地学系Watson为代表的空外接球法亦属于加点法。加点法算法简明， 是目前应用最多的算法。分割占有法(Devide and conquer):将数据域递归细分为子 块，每块实现局部优化三角化，然后合并。对应于上述三类算法各有一些著名的 算法。Bo

17、wyer 算法该算法基于Dirichlet图的构造，适用于n维空间中的离散点集的网格剖分， 是英国Bath大学数学分校的A.Bowyer在总结该校Robin Sibson教授、 PeterJ.Green博士七十年代所做工作的基础上，于1981年提出的。算法思路Bowyer算法是针对Voronoi图的实现，首先开始于一个由n+1个数据点形 成的Delaunay单纯体，这样将得到一个包含一个真实顶点、而其它顶点为无穷 远点0的Voronoi图:如图2-17所示，往己有数据形成的Voronoi图中加入新点Q， 从包含新点的Voronoi多边形开始，利用Voronoi多边形相邻关系，找到最近点， 构造

18、新加入点的Voronoi多边形(粗虚线)。图2-L7在已有山”叫皿图中如入新泌堂成其MWQjKi争讪莎的过程算法分析该算法整体效率为O(n(k +1/k),k为空间的维数。基于Voronoi图的Bowyer 算法计算复杂，空间消耗大，算法时空效率不高。Watson 算法该算法是澳大利亚悉尼大学Geology与Geophyics系D.F Watson于1981年 提出的，是用计算机构造品体模型的研究成果。算法思想首先给出由一个或多个外接球不包含其它数据点的单纯体组成的初始网格， 然后往其中加入一个数据点，考察外接球的包含情况，去除包含新点的n维单纯 体，用n+2个点能组合的、外接圆不包含其它点的

19、单纯体取代。如图2-18所示 是Watson算法的具体思路。具体实现时，可一次性全部找出并删除所有包含新 加入点的单纯体，得到一个包含新点的空洞，空洞的边界面与新点相连，得到新 点加入后的Delaunay网格，这样可避免进行新生成单元的外接圆是否包含old points的计算，具体加入一点的算法流程叙述如下。算法流程加入新点，搜索单纯体链表，查找外接球包含新点的单纯体，所有包含 新点的单纯体组成一个多面体。包含新点的单纯体的各个面加入一个临时链表，若一个面加入两次，说 明是两个单纯体共享的面，必然位于多面体的内部，从链表中删除，若出现新点 位于外接球上的退化情形，则抛弃链表和新点，改用其它方法

20、处理。若未出现退化情形，则将新点与多面体的各个面相连，得到新的单纯体。新点加入过程结束。算法分析Watson算法的思路非常简明，易于编程实现。但当出现k+2（k为空间维数） 点位于同一圆球面时，则三角化结果不唯一，这种情形称为退化情形。如图2-19, 为二维空间的退化情形。除了如图2-20所示的规则数据，实际应用中散乱数据 点集很少出现退化情形。但由于计算机的计算精度是有限的，当新点与外接球球 面之间的距离小于预设计算精度，则认为新点位于球面上;这种新点与球面距离 的计算误差可能引起拓扑关系的不一致。2.3.3.3四叉树、八叉树方法四叉树、八叉树本身是数据结构，当用于空间编码时，可进行实体造型

21、，适 当修改即成为网格剖分算法。M A Yerry.M S Shephard于 1983年、1984年发表 13 了四叉树、八叉树在二维、三维网格剖分中的应用，有些文献中将他们的算 法称为Shephard-Yerry算法。Shephard-Yerry算法只适用于域剖分。网格单元可 以是四边形/六面体，也可以是三角形/四面体。这种算法的特点是与实体造型相 结合，自动化程度高，网格密度可调整，剖分速度快，内部单元形状比较好，但 边界单元形状很难保证。在二维空间，以矩形网格为例，边界单元共有16种被 截情形;三维六面体网格单元被边界面截切后共有4096种情形，被切后的单元与 原有网格单元拓扑关系可能

22、不一致，需要进行拓扑变换。此外，这种方法不具备 几何不变性，即剖分对像旋转后，剖分结果发生变化。韩国中央大学YH Jung和K Lee于1993年提出了一种直接基于四面体的八 叉树空间编码法。这种算法一步到位，内部单元为四面体，边界单元易于处理为 四面体，不需要拓扑变换。总之，Shephard-Yerry算法由于与实体造型技术相结 合，是一种很有前途的方法，一旦边界单元得到很好处理，则必在实体网格剖分 方面占据主要地位。2.3.3.4换边换面1977年，美国加利福利亚工学院喷气推进实验室的Charles L.Lawson提出了 基于边交换的二维Delaunay三角化。加拿大埃德蒙顿Albert

23、a大学计算机科学系 Barry Joe分别于1989年、1991年给出了基于局部换面的三维Delaunay三角化 的算法和证明。算法分析换边换面法适用于离散点集剖分和域剖分，算法过程简明，容易实现。但是 三维三角化算法过程中需要检测的面很多，有效控制变换的范围，合理的数据结 构和快速的查询法是提高速度的关键。换边换面法的最大的困难在于如何处理不 可变换的情形，这是本文研究的关键问题。只有解决不可变换的情形，才能得到 delaunay三角剖分网格；否则，结果是非delaunay的。2.3.3.5网格的前言生成法。网格的前沿生成法从提出到现在发展最快，现在有了很多的算法。最早法国 学者S H Lo

24、于1985年在文献中提出了网格前沿法的雏形，只将网格前沿法作为 节点连元的方法，没有与节点生成同时考虑。英国学者J.Peraire,M.Vabdati,K.Morga.等于1987年实现了按方向细化的网格前沿法。此后 研究的各种网格前沿法最大的特征是能够生成复杂形状的非结构网格，其按方向 细化的特点特别适合于三维可压缩流的优化算法。算法思路以剖分域的边界为网格的初始前沿，按-设定的网格单元的形状，尺度等要 求向域内生成节点，连成单元，同时更新网格前沿，如此逐层向剖分域内推进， 直到所有的空间都被剖分。算法分析网格前沿法能够处理比较复杂的对像，其主要特点是提供了控制网格密度和 质量的手段。但是网

25、格前沿法中存在大量的查询操作以及网格前沿面的相交检 测，这些操作是很浪费时间的。很多算法在这两点上作了快速的处理。下图是网 格前沿法生成的三角网格实体。第3章空间曲面上散乱数据点的快速三角剖分算法实现3.1采用的数据结构点类，点节点类，点链表类，边类，边节点类，边链表类（多边形类），三 角形类，三角形节点类，三角形链表类。以下对这几种结构做简单的介绍（只介绍主要的）：点类中的数据成员class Point_Tdouble m_X;/空间点的坐标double m_Y;double m_Z;下面是一些方法点节点类中的数据成员。class PNode_TPont_T m_Point;PNode_T*

26、m_Left;PNode_T*m_Right;下面是一些方法点链表类中的数据成员。class PNList_TPNode_T*m_Head;PNode_T*m_Tail;下面是一些方法略边类中的数据成员。class Edge_TPont_T m_Point1;Pont_T m_Point2;Pont_T m_Point3;/所在三角形第三个点bool m_Used;/表示这条边是活边还是死边/下面是一些方法略边节点类中的数据成员。class ENode_TEdge_T m_Edge;ENode_T*m_Right;ENode_T*m_Left;/下面是一些方法略多边形类中的数据成员。class

27、 Polygone_TENode_T*m_Head;ENode_T*m_Tail;Int m_LivingEdgeNum;/边界前沿中的火边数量。下面是一些方法略三角形类中的数据成员。class Triangle_TPont_T m_Point1;Pont_T m_Point2;Pont_T m_Point3;/所在三角形第三个点/下面是一些方法略三角形节点类中的数据成员class TriNode_TTriangle_T m_Triangle;TriNode_T*m_Right;TriNode_T*m_Left;/下面是一些方法略三角形链表类中的数据成员class TriNList_TTriN

28、ode_T*m_Head;TriNode_T*m_Tail;/下面是一些方法略点链表：存储点云数据的一个单向链表。三角形链表：生成的三角形存储在单链表中。前沿边界（多边形）：多边形采用边的双循环链表的形式存储。3.2三角剖分的过程定义：内边：三角剖分结果中被两个三角形公用的边。外边：三角剖分的结果中只被一个三角形拥有的边。内点：如果一个点的相连边都是内边这个点就是内点。活边：没有经历找点生成新边过程的边。这里要强调，活边经历找点生成新 边过的程，可能找到了匹配点，可能没找到匹配点，可能会产生0条新边，1条 新边或2条新边。死边：经历一次找点生成新边过程的边就是死边。死边可以是边界边，也可 以是

29、外边。外边框：由外边首尾相连构成的空间多边形就是外边框。边界点：在外边框中的顶点。外点：没有被选择过的点。3.2.1算法简单描述：（1）：点云的预处理。（2）：构造一个相对饱满的初始三角形作为种子三角形。（3）：开始循环：如果外边界框中的活边的数量不是零，就继续下面的步骤。 否则算法结束。（4）：对外边框做循环：ENMov=:外边框的头节点。（5）：如果ENMov是活边：选择匹配的点，如果选择到了匹配点就更新点 链表，更新三角形链表，更新外边框链表；如果没选择到匹配点，把ENMov标 志成死边used=1。如果ENMov是死边：什么也不做。（6）： ENMov=：外边框更新前ENMov的下一个

30、节点。如果ENMov是空 的，那么返回（3），否则返回（5）。以下对算法中的各个步骤做详细的解释。3.2.2数据点的预处理（1）：把点云数据文件中的点读到点单链表中。（2）： PNMov=:点链表的头节点。（3）：在链表中删掉PNMov节点后面到PNMov距离小于DIST的所有的节 点。其中DIST是一个可以指定的数，其数值越大剩下的点越少。（4）： PNMov=:点链表中PNMov的下一个节点。（5）： PNMov如果不是链表的尾节点：返回（3）。否则结束。处理后的点 相对少得多，也保证处理后的点云要相对的均匀，最主要的是去掉了曲率过大的 点。3.2.3初始三角形的建立（1）：先得到点云链表

31、中的第一个点作为firstPoint。（2）：然后得到在firstPoint后面，距离firstPoint最近点作为第二个点 secondPoint。（3）：这两个点构成了一条边。在secondPoint后面的点中寻找距离第一条 边的两个端点的距离和最小的点作为第三个点thirdPoint。（4）：如果这个三角形的最小的内角小于30度。那么选择点云链表firstPoint 节点的下一个节点为firstPoint，返回（2）。否则，结束。3.2.4活边寻找最佳匹配点的算法为了加快运算的速度我们采用包围盒算法。一个点要成为当前边的候选点要 具备的必要条件。i:这个点和当前边构成的三角形于当前边所在的三角形钩成的面角要大于给 定的阀值 MINFACEANGLEO卜面介绍空间中两个面的二面角的计算方法。三角形p1p2p3和三角形p1p2p4所成的二面角的大小等于向量p5p4与向量 p6p3所成的角P1p5是p1p4在p1p2上面的投影向量,p1p6是向量p1p3在