三维Delaunay算法

将Delaunay剖分算法推广到三维具有重要意义。三维Delaunay剖分构成的Delaunay四面体是进一步构成任意块体的中间工具。三维空间的四面体对于三维射线追踪非常方便。三维Delaunay算法基本原理与二维Delaunay算法十分相似，但编制起来更为复杂。它数庆也要用一个n×8维数组记录四面体构成点和相邻四面体信息，它大致也分为以下几步：

（1）判断哪些四面体的外接球包含新加入的点；

（2）将这些四面体汇总到一块，形成一个凸多面体；

（3）找到这多面体外表面，用一个二维三角形相邻关系数组记录下来；

（4）由多面体表面的三角形与新加入的点构制新四面体，用一个三维四面体数据结构数组，存贮新形成的四面体信息；

（5）用新四面体信息去更新原来总的四面体数据结构信息。

其中第（3）步是比二维Delaunay剖分复杂。待修改的四面体的外表面不再能够用环表示，而要用一个二维Delaunay三角形数据结构来表示。这种三角形的相邻关系要从原四面体的数据结构关系中去寻找。A、C面是多面体的外表面。B是其内表面须寻找A的一个棱的相邻三角形B，薯扮握从B找到下一个四面体，再从四面体上找到与B三角形棱相邻的面C。

在三维Delaunay算法中，在通过外接球找到所有待修改的四面体后，将其排成一个待修改四面体数据结构缺历数组，将待修改四面体的所有面全部排列起来，其中包含各个四面体的相互公共面和这些四面体组成的多面体的外表面，构成一个关于面的数据结构数组。从这个数据结构中删去相互相邻的三角形面，这样就构成了外表面的数据结构数组。

三维Delaunay剖分时，我们总利用三角形网格中的数量关系，来检查所形成的四面体的正确性。如果待删四面体数目为T，外表面三角形的数目为F，则

地质模型计算机辅助设计原理与应用

从上文中介绍的对点集的Delaunay三角剖分方法可以空游败知道，点集的Delaunay的三角剖分方法是最简单、最基础的，所以，可以考虑在点集的Delaunay三角剖分方法的基础上获得区域的Delaunay三角剖分，实际上，区域的Delaunay三角剖分方法正是这样做的。

为了利用点集的Delaunay三角剖分方法，首先应该将剖分对象从区域转换到点集，而为了实现上述目的，需要经过两个过程:布点、离散边界。布点即是向区域内插入一些合适的点，离散边界即是将区域边界分割成若干小线段。在进行上述处理后，就将需要剖分的对象从一个区域转变成一个点集，点集中的点为向区域内插入的点和边界离散后形成小线段的端点然后采用点集的Delaunay三角剖分方法，如Lawson算法或Bowyer-Wat-son算法，对从区域转换过来的点集进行Delaunay三角剖分 *** 作最后需要删除那些在区域之外的三角形，因为对一个点集进行Delaunay三角剖分 *** 作，最后获得的三角剖分的范围是该点集的凸壳，而点集的凸壳绝大多数情况下与需要剖分的区域不一致，通常情况下是凸壳的范围大于需要剖分的区域，因此那些属于凸壳但不属于需要剖分的区域的三角形需要删除。

区域的Delaunay三角剖分方法可以概括为3个步骤:

第一步:布点并离散边界，将剖分对象从区域转换为点集。根据剖分规模或想要得到的三角形单元的大概边长，向区域内插入一系列的内部点，并将内外边界离散成一系列的小线段。图3.14(a)为某剖分区域，图3.14(b)为向区域内布点的结果，图3.14(c)则时在布点之后进行离散边界的结果。

图3.14 二维区域布点与边界离散

第二步:对点集进行Delaunay三角剖分 *** 作。将离散后的内外边界的端点和根据需要布设的内部点组成输入点集，采用逐点插入法(Lawson算法或Bowyer-Watson算法)将上述点集进行Delaunay三角剖分，详细步骤请参见本章中点集的Delauany三角剖分方法，最后点集的Delaunay三角剖分如图3.15所示。剖分区域之外的三角形只能是由同一条边界上线段的端点组成，不可能由内部点或不同边界(如2个顶点位于外边界，1个顶点位于内边界)上端点组成区域之外的三角形，所以，判断一个三角形是否多磨模余，首先判断该三角形的3个点是否位于同一条边界，若是，该三角形有可能为多余三角形，但不一定，如图3.16中由点9、10和11组成的三角形位于区域内部，但由点21、22和23组成的三角形却在区域外因此，三角形的3个顶点位于同一条边界上只是必要条件，还需要进一步判断。

图3.15 点集的Delaunay三角剖分(未进行边界处理)

第三步:边界处理。根据区域边界，删除多余的三角形，得到区域的Delaunay三角剖分。

完整的边界处理方法如下:

(1)将区域的每条边界离散后按照固定顺序排列，边界上端点的ID也依次排列。需要特别注意的是，外边界必须按照逆时针排列，外边界上的端点编号也依次由小到大排列，如图3.16所示，从点0到点36依次排列所有内边界离散后按照顺时针排列，同样，每条内边界上端点编号依次从小到大排列。

(2)依次判断点集的Delaunay三角剖分中每个三角形的三个顶点是否位于同一边界，若某个三角形的3个顶点位于同一边界上，将3个顶点由小到大排序，并按照排序后点号组成一个新的三角形，计算该新三角形的面积，若面积＞0，则表示该三角形为区域内部三角形，保留若面积＜0，则表示该三角形为区域外三角形，删除。

图3.16 边界处理

图3.17 二维区域的Delaunay三角剖分

如由点0、1和36组成的三角形，排序后依次为0、1、36，新的三角形面积＞0，表示该三角形为区域内部三角形，保留由点21、22和23组成的三角形按照逆时针规则正常时排序为23、22、21，排序后依次为21、22、23，新的三角形面积＜0，表示该三角形为区域外部三角形，删除。

边界处理后，得到最终的区域的Delaunay三角剖分，如图斗颤3.17所示。

布点和离散边界的代码详见推进波前法代码中的函数CreateInnerPoint()和Create-BoundaryPoint()点集的Delauany三角剖分方法的完整代码详见Bowyer-Watson算法的代码。以下为进行边界处理的代码，函数DeletedInvalidTriangle()用于删除剖分区域之外的三角形，函数AreaOfTrgl()用于计算三角形的面积。

三维地质建模方法及程序实现

三维地质建模方法及程序实现

随机增量法较为简单，遵循增量法的一贯思路，即按照随机的顺序依次插入点集中的点，在整个过程中都要 维护并更新一个与当前点集对应的Delaunay三角剖分盯宏慎。考虑插入vi点的情况，由于前面插入所有的点v1,v2,...,vi-1构成的DT(v1, v2,...,vi-1)已经是Delaunay三角剖分，只需要考虑插入vi点后凯敬引起的变化，并作调整使得DT(v1,v2,...,vi-1) U vi成为新的Delaunay三角剖分DT(v1,v2,...,vi)。

(插绝庆入调整过程：首先确定vi落在哪个三角形中（或边上），然后将vi与三角形三个顶点连接起来构成三个三角形（或与共边的两个三角形的对顶点连接起来构 成四个三角形），由于新生成的边以及原来的边可能不是或不再是Delaunay边，故进行边翻转来调整使之都成为Delaunay边，从而得出DT (v1,v2,...,vi)。)

其Pseudocode如下：

Algorithm IncrementalDelaunay(V)

Input: 由n个点组成的二维点集V

Output: Delaunay三角剖分DT

1.add a appropriate triangle boudingbox to contain V ( such as: we can use triangle abc, a=(0, 3M), b=(-3M,-3M), c=(3M, 0), M=Max({|x1|,|x2|,|x3|,...} U {|y1|,|y2|,|y3|,...}))

2.initialize DT(a,b,c) as triangle abc

3.for i <- 1 to n

do (Insert(DT(a,b,c,v1,v2,...,vi-1), vi))

4.remove the boundingbox and relative triangle which cotains any vertex of triangle abc from DT(a,b,c,v1,v2,...,vn) and return DT(v1,v2,...,vn).

Algorithm Insert(DT(a,b,c,v1,v2,...,vi-1), vi)

1.find the triangle vavbvc which contains vi // FindTriangle()

2.if (vi located at the interior of vavbvc)

3.then add triangle vavbvi, vbvcvi and vcvavi into DT // UpdateDT()

FlipTest(DT, va, vb, vi)

FlipTest(DT, vb, vc, vi)

FlipTest(DT, vc, va, vi)

4.else if (vi located at one edge (E.g. edge vavb) of vavbvc)

5.then add triangle vavivc, vivbvc, vavdvi and vivdvb into DT (here, d is the third vertex of triangle which contains edge vavb) // UpdateDT()

FlipTest(DT, va, vd, vi)

FlipTest(DT, vc, va, vi)

FlipTest(DT, vd, vb, vi)

FlipTest(DT, vb, vc, vi)

6.return DT(a,b,c,v1,v2,...,vi)

Algorithm FlipTest(DT(a,b,c,v1,v2,...,vi), va, vb, vi)

1.find the third vertex (vd) of triangle which contains edge vavb // FindThirdVertex()

2.if(vi is in circumcircle of abd) // InCircle()

3.then remove edge vavb, add new edge vivd into DT // UpdateDT()

FlipTest(DT, va, vd, vi)

FlipTest(DT, vd, vb, vi)

复杂度分析

问题的规模为点集中的点的总个数n（没有重合的点），循环内的基本的 *** 作有：

1.寻找插入点所在的三角形(FindTriangle())

2.测试点是否在外接圆内(InCircle())

3.更新三角网(UpdateDT())

4.寻找共测试边的第三顶点(FindThirdVertex())

考虑最坏情况：

时间复杂度T ＝ T(addboudingbox()) + Sum(T(insert(i), i=1,2,...,n)) + T(removeboundingbox)

因为addboudingbox()和removeboundingbox不随n变化，是个常量。T(addboudingbox()) = O(1), T(removeboundingbox()) = O(1).

T = Sum(T(insert(i), i=1,2,...,n)) + O(1) + O(1).

考虑插入第i个的点的时间:

T(insert(i)) = T(FindTriangle(i)) + T(UpdateDT(i)) + K*T(FlipTest(i)) (K为常数)

故T = Sum(T(FindTriangle(i)), i=1,2,..,n) + Sum(T(UpdateTD(i)), i=1,2,..,n) + K*Sum(T(FlipTest(i)), i=1,2,..,n)

挨个考虑：

FindTriangle(i)是要找出包含第i个点的三角形，由欧拉公式知道，平面图的面数F是O(n), n为顶点数，故此时总的三角形数是O(i)的。所以问题相当于在O(i)个记录中查找目标记录，如果不借助特殊的数据结构，按照一般顺序查找，需要O (i)的时间。

T(FindTriangle(i)) = O(i),故：Sum(T(FindTriangle(i)), i=1,2,..,n) = O(n*n)

UpdateTD(i)是更新三角网数据结构，插入和删除三角形到当前的三角网是个常量 *** 作，因为已经知道插入和删除的位置。

T(UpdateDT(i)) = O(1),故：Sum(T(UpdateTD(i)), i=1,2,..,n) = O(n)

FlipTest(i)比较复杂些，是个递归过程。细分为：T(FlipTest(i)) = T(FindThirdVertex(i)) + T(InCircle(i)) + T(UpdateDT(i)) + 2*T(FlipTest(j))。(这里，用j来区分不同的深度)

因为InCircle(i)是测试点是否在外接圆内，是个常量 *** 作，不随点数变化而变化。故T(InCircle(i)) = O(1), 又T(UpdateDT(i) = O(1)(见上)

FindThirdVertex(i)是查找目标点，在O(i)个记录中查找目标记录，如果不借助特殊的数据结构，按照一般顺序查找需要O(i)的时间。T(FindThirdVertex(i)) = O(i).

剩下的是递归调用FlipTest的过程，不过还好，因为FlipTest的递归深度只有常数级(设为K)。正比于点在三角网中的度数(degree)。故FlipTest(i)最多调用的次数为2*2*,...,2 = K',还是常数。

故T(FlipTest(i)) = O(i) + O(1) + O(1) + K'*(O(i) + O(1) + O(1)) = O(i) + O(1) + O(1) .

Sum(T(FlipTest(i)), i=1,2,..,n) = O(n*n) + O(n) + O(n)

综上，最坏情况下算法总时间复杂度 T = O(n*n) + O(n) + K*(O(n*n) + O(n) + O(n)) + O(1) + O(1) = O(n*n)

其中，关键的 *** 作是FindTriangle()和FindThirdVertex()，都是n*n次的。

在网上很多资料说随机增量法是O(nlogn)的，但是分析下来，却是O(n*n)的。后来看到别人的实 现，原来是用的别的数据结构来存储三角网，减少了FindTriangle()和FindThirdVertex()的复杂度。使得某次查找三角形和共边 三角形的第三顶点能在O(logn)，而非O(n)的时间内实现。这样，总的查找的时间为 O(log1)+O(log2)+,...+O(logn) = O(nlogn)。程序=算法+数据结构，看来一点没错。

比如说用DAG，Quad-edge等，都能达到O(nlogn)的复杂度。

分治法（Divide and Conquer）

据说是现在实际表现最好的。

---