基于Delaunay三角剖分和区域生长的散乱点云重构

王先泽，李忠科，马亚奇，张晓娟

(第二炮兵工程大学 401教研室，西安 710025)

近年来，对无规则点云表示的模型的曲面重构已经在逆向工程、计算机视觉、计算机图形学和虚拟现实中得到了广泛的应用。目前解决重构问题的方法可以划分为3种:隐式曲面法、基于Delaunay三角剖分法和区域生长法。

隐式曲面法首先从点云数据中定义一个带符号的距离函数，然后利用该函数的零值集来近似生成三角网格曲面。Hoppe等［1］提出的重构方法就是基于隐式曲面法。但是该方法不能插值点云中的点，在CAD系统中的应用受到了限制。基于Delaunay三角剖分法首先通过Delaunay三角剖分或是Voronoi图从散乱点云数据中建立几何邻近结构，然后从结构中提取出一些面片来近似表示曲面。该方法很好地建立了散乱点云内部的几何信息，具有很强的鲁棒性和系统性。但是由于使用了中间数据结构，计算量大，时效性比较差。Edelsbrunner等［2］提出的著名的alpha shapes算法和Dey等［3］提出的co-cone算法都属于基于Delaunay三角剖分的算法。区域生长算法先初始化一个三角形作为初始的区域，然后在区域的边界迭代加入三角形生成曲面。比较有代表性的区域生长算法有 Bernardini等［4］提出的 ball-pivoting算法和最近 Cohen-Steiner［5］提出的基于 Delaunay的贪婪算法。区域生长算法的优点是效率高，缺点是网格的构建质量过分依赖种子面片的选取以及自定义的参数，在采样密度不均匀的情况下容易留下孔洞，需要后续的孔洞修补算法以构建封闭的曲面。

本文提出了一种基于Delaunay三角剖分和区域生长相结合的散乱点云重构方法，首先计算采样点集的Delaunay三角剖分，得到点集的Delaunay剖分结果，然后用区域生长算法从中“生长”出三角网格曲面。生长的过程是先初始化一个三角形作为初始区域，然后仅在区域的边界边上迭代添加新的三角面片，这样在迭代的过程中生成的区域会不断扩展。图1为实验模型horse的区域生长过程的展示。与传统的区域增长算法相比，本文的算法在区域生长的过程中使用自定义的可接受度因子作为判断添加新三角面片的标准，在一定程度上避免了狭长的三角形，提高了网格质量。区域生长完成后，对其中的孔洞进行后处理，可使生成的曲面更加完整。

图1 重构模型horse的区域生长过程

1 基本概念

首先，对本文所用到的一些主要符号作一下简要说明:P表示所要重构的散乱点云数据集;S表示重构出的三角网格曲面;∂S表示区域生长过程中边界边集;t表示三角形;v表示顶点;e表示边;tseed表示选取出的种子面片。

定义1:设三角形t的外接圆半径为rt，他是指通过三角形顶点的不包含其他采样点的最小圆的半径。

在Voronoi图中，每个Voronoi点恰好是三条Voronoi边的交点［6］。也就是说，Voronoi点是由点集P中3点形成的三角形的外接圆的圆心，因此，三角形t的外接圆半径rt也可以定义为t的任意一个顶点到其对偶Voronoi边的距离。在本算法的计算过程中所用的就是这个定义。

在区域生长算法中，迭代添加三角面片时通常使用的原则是添加外接圆半径rt最小的三角形［5］，这样对于一些狭长但外界圆半径小的三角形也被添加进去，影响生成网格的质量。在本算法中，提出了一种新的解决方案，自定义一个参数因子可接受度Acceptance(t)作为区域生长中迭代添加三角面片的标准。

定义2:在区域生长中，一个三角形t作为迭代添加的三角形的可接受度Acceptance(t)为

其中:remin表示所有与边e相关的三角形的外接圆中的最小半径;aud(t)是指三角形的最小内角与π/3的比值。对散乱点云进行三角剖分的时候，理想状况就是三角网格的各个单元都是等边三角形，但这一目标通常难以实现，通常只能是三角形的最小内角最大，这也是判断网格质量的一个重要标准。三角形t的最小内角可以通过aud(t)反映出来。aud(t)越大表示三角形3条边越均匀，则t的最小内角越大，网格质量也越好。当aud(t)等于1时，表示三角形为等边三角形。

权重α和β的值可以调整。α越大，表示候选三角形倾向于外接圆半径最小的三角形;β大则表示对三角网格的质量要求较高，但α和β的取值应该满足α+β=1。通过对α和β的取值，当2个三角形的外接圆半径相差不大时，可以选取到能使网格质量更好的三角形。Acceptance(t)越大，表示该三角形下一次被添加的可能性越高。

2 基于Delaunay和区域生长的散乱点云重构算法

2.1 Delaunay三角剖分

近年来，在曲面重建算法研究领域中，Delaunay三角剖分和Voronoi图已经得到了广泛的应用。一个三维的Delaunay三角剖分后的结果就是将采样点形成一个四面体的凸包，每个四面体包含4个要素:顶点，边，三角形和四面体。在所形成的四面体中，其外接球内都不包含其他的采样点，这也是使用Delaunay三角剖分对散乱点云进行重构的一个重要优点［7］。

本文提出的算法是直接使用CGAL算法库来计算点云的三维Delaunay三角剖分DT，数据结构也是采用CGAL提供的Delaunay三角剖分的数据结构。使用CGAL进行Delaunay三角剖分采用了逐点插入的算法策略，包括点的定位和正确性检查等过程。点的定位有2种方法:一种是通过查找与定位点最近的初始点，再由初始点沿一直线搜索定位;另一种是使用一种层次的数据结构，随机采样，完成定位。正确性检查首先检查三角化数据结构的一致性，每个面与其邻面具有3个共点，并具有连续编号。总的点线面数量也要做检查。进而更改每一个体元的朝向，以及顶点凸集的正确性。三角剖分完成后，可以很容易地得到任意的边或是顶点相邻的四面体和三角形，这为获得后面区域生长中所必需的几何信息提供了捷径。

2.2 区域生长算法

区域生长方法通常首先从点云数据中选定一个初始三角片，并作为网格曲面的种子面加入三角曲面集，3条边分别加入边界集。然后按照一定的规则从边界处向三角曲面集中添加新的三角片元素，并实时刷新边界集，通过曲面网格由局部到整体的动态增长最终生成一张完整的三角网格曲面。

2.2.1 种子面片的选取

区域生长的第1步就是选取种子三角面片tseed。首先就是从三角剖分后的点集中选出z坐标值最大的点vi，然后计算所有与其相关的三角形的外接圆半径rt和角均匀度aud(t)，得到三角形的可接受度Acceptance(t)，选取可接受度最大的三角形作为种子面片，加入到曲面S中。同时，也将该三角形的3条边e加入到边界边集∂S中，作为曲面的边界边。另外，种子三角面片tseed的法向量应与z轴正方向一致，也就是说，tseed的法线方向与z轴正方向的夹角必须小于π/2。

2.2.2 有效三角形的几何和拓扑约束

在本算法中，要求所构建出的曲面S必须为封闭的或带边界的二维流型，这也是我们希望得到的好的构建曲面。在这里，采用文献［8］中对几何一致性的证明，通过4种拓扑结构得到想要构建的流型曲面。假设e为构建曲面S边界边，t为e所在的三角形，v为t中与e所对的点，则这4种拓扑关系分别为面延伸(v∉S)、填充洞(v∈∂S，并且v在∂S中的2个邻点同时是e的2个端点)、填充耳(v∈∂S，并且v在∂S中的1个邻点是e的端点)和粘贴面(v∈∂S，但是v在∂S中的2个邻点都不是e的端点)。4种拓扑结构如图2所示，满足这4种结构之一的都认为t是边e的有效三角形，反之则是无效三角形。

在这4种拓扑结构中，面延伸、填充洞和填充耳都不会出现非流型的点或边，但粘贴面操作后会出现网格仅通过1点相连的情况，在点云噪声过大且密度较低的区域执行此操作时可能造成平面不可定向并引起错误，通常是通过增加附加三角形t'来避免这种情况。添加三角形时得满足2个条件:一是t'中必须有1条边的端点作为顶点，并且在∂S中与v相关的边作为对边;二是加入t和t'后不能导致曲面不可定向。

2.2.3 候选三角形的选取

种子面片选择好之后，下一步就是按照一定的标准选择一个新的三角形加入到构建曲面S中。这可以分为2步来完成:首先从∂S中的每1条边界边e的所有有效三角形中选出候选三角形，然后从所有的候选三角形中选择出要添加的三角形。

图2 4种拓扑结构

在从边界边的有效三角形中选择候选三角形时，选取的是可接受度最大的三角形，但是由于有片状四面体的存在，有可能出现生长方向往回长的情况，从而引起错误。解决这一问题的方法是对三角形与构建曲面间的二面角加以限制。假设与边界边e∈∂S相关的三角形为t，包含在曲面S中的与e相关的三角形为tb，θt表示t与tb之间法线的夹角的绝对值，αsliver是一个常数，则边界边e的候选三角形ce可以表示为

式中:t是边e的有效三角形，并且qt＜asliver。如果边界边e的有效三角形中没有满足上面条件的，则ce为空。αsliver的取值并不是严格规定的，在本算法中将其设置为5π/6。在计算可接受度Acceptance(t)时，根据先验知识，取α=0.8，β=0.2。

2.2.4 三角面片的迭代添加

区域生长中每次迭代只能将一个候选三角形添加进构建曲面S中，文献［5］中采用贪婪算法，每次都选择最有可能的候选三角形加入，可以避免前面提出的启发式搜索策略在采样点稀疏时可能失效的缺陷。本算法在选取要添加的三角形时，主要是以三角形与曲面S之间的二面角为标准，在采样密度大的时候，考虑到曲率问题，二面角越小则曲面质量越好。但是，当2个候选三角形与曲面S的二面角都很小或者是有噪声的时候，这个选择标准就不适用。在这里同样需要引入可接受度Acceptance(t)，对于那些二面角θt小于1个常数θ的三角形t来说，就需要通过他们的可接受度来排列其优先级。因此，可以对候选三角形t的优先级P(t)作如下定义

当t为空集时，优先级P(t)为－∞。根据实验，设置θ=π/6。

3 后处理

在实际的操作过程中，由于噪声等因素的干扰，构建出来的曲面S通常会含有一些小的孔洞。本文参照文献［9］中提出的三角网格模型的补洞算法，首先根据网格中的点、边和三角形之间的关系提取孔洞边界，然后采用最小角度法在孔洞中依次填补三角形，接着对一些新增加的高度弯曲的三角形进行细分操作，最后对修补后的孔洞网格进行几何形态调整和曲面光顺化。

4 实验结果

本算法已经用VC++2008在2.8GHz Intel CPU，2G内存的计算机上编译实现。表1为部分模型的实验数据，按类别分为了封闭模型和带边界的模型。对于带边界的模型，采用启发式搜索策略，计算所有三角面片外接圆半径的平均值radius_aver，设置阀值γ，将那些外接圆半径大于γ倍radius_aver的三角面片删除。根据实验分析，γ取值通常在5～100。图3为一些重构出的实验模型效果。

表1 几个实验模型的实验数据

图3 一些实验模型效果

5 结束语

本文提出的基于Delaunay三角剖分和区域生长相结合的对散乱点云进行三角网格重构的方法，是在对散乱点云进行三角剖分的基础上采用区域生长算法“生长”出三角网格曲面。与其他的区域生长算法相比，本算法根据自定义的可接受度因子作为区域生长过程中迭代添加三角形的判断标准，可以在一定程度上避免狭长的三角形，提高网格质量。通过对一些散乱点云数据进行实验，结果表明，该算法是快速、稳定、有效的。但是，该方法还存在一些缺陷，需要自定义的参数很多，没有一个统一的标准，特别是在对Acceptance(t)中α和β的取值上，如果β取值过大，则区域生长的过程中有可能出现添加一些外界圆半径很大的三角形，使区域生长不能完成。

［1］ Hoppe H，Derose T，Duchamp T，et al.Surface reconstruction from unorganized points［C］//SIGGRAPH，［S.l.］:［s.n.］，1992:71 －8.

［2］ Amenta N，Bern M，Kamvysselis M.A new Voronoi-based surface reconstruction algorithm［C］//SIGGRAPH，［S.l.］:［s.n.］，1998:415 －21.

［3］ Dey TK，Giesen J，Leekha N，et ct.Detecting boundaries for surface reconstruction using co-cones［C］//Int J Comput Graph CAD/CAM，［S.l.］:［s.n.］，2001(16):141 －59.

［4］ Bernardini F，Mittleman J，Rushmeier H，et al.The ball-pivoting algorithm for surface reconstruction［J］.IEEE Trans Vis Comput Graph，1999，5(4):349 －59.

［5］ David C S.A greedy Delaunay based surface reconstruction algorithm［R］.Research Report，INRIA，2002.

［6］ 周培德.计算几何－算法与设计与分析［M］.3版.北京:清华大学出版社，2008:127－128.

［7］ Chuan-Chu Kuo，Hong-Tzong Yau.A Delaunay-based region-growing approach to surface reconstruction from unorganized points［C］//Computer-Aided Design，［S.l.］:［s.n.］，2005(37):825 －835.

［8］ Huang J，Menq C H.Combinatorial manifold mesh reconstruction and optimization from unorganized points with arbittary topology［J］.Computer-Aided Design，2002，34(2):149－165.

［9］ 田建磊，刘旭敏，关永.三角网格模型的补洞算法研究［J］.计算机应用，2009，29(8):2035 －2037.

(责任编辑鲁 进)