基于Voronoi最小邻近点集的Delaunay三角化方法

孟宪海，成文迪，徐 博，杨 钦

（1.北京航空航天大学计算机学院，北京 100191；2.中国石油天然气勘探开发公司，北京 100034）

1 绪 论

网格生成，即是将复杂几何结构离散成一组简单几何单元的过程，生成三角形或四面体网格的过程又称三角化。其中，Delaunay三角形/四面体网格具有很高的质量，对复杂区域有很好的逼近性，应用最为广泛。

Delaunay三角网格在1934年由俄国数学家Delaunay提出[1]，在此类网格中，任何一个三角形/四面体网格单元的外接圆/球中，均不包含网格中的其他顶点。Delaunay三角网格以其优良的几何性质，在工程中得到了十分广泛的应用。

点集的Voronoi图，由每个点的Voronoi单元组成，点a的Voronoi单元，定义为空间中所有到a点的欧氏距离不大于到点集中其他任何点的欧氏距离的所有点的集合。点集的Delaunay三角化结果，和点集的Voronoi图互相对偶。即如果点集中任意两点的Voronoi单元之间存在公共面，则这两点的连线即为Delaunay三角化中的一条边。二者是彼此等价的。

传统的点集Delaunay三角形/四面体网格化算法，可以生成高质量的网格，还能处理一定的限定条件[2-3]。但近年来，由于高精度科学计算，数值仿真等领域的需求，需要处理更大数据量、生成更高精度的网格，还对网格的生成次序、形态等有特殊的要求，这些因素促使人们不断寻求新的网格生成思路。

经典的Delaunay三角化生成算法已经较为成熟。其中以B/W增量算法应用最为广泛[4-5]。在用B/W算法进行点集的Delaunay三角化时，首先用少量点生成初始的粗网格；之后，每插入一个新点时，检查哪些三角网格单元的外接球包含这个点，把他们删去，形成一个空洞；最后将新点与组成空洞的每个顶点相连，形成新的网格。不断重复这个过程直到所有点都被加入到网格之中，算法终止[6]。

在这种思路下，每插入一个点，都要在现有的整体网格中搜索外接球包含它的网格单元，每一个新网格的生成都与当前整体网格相关，算法具有全局性、连续性。虽然效率高，但当存在某些限定条件、或者无法得到全局信息时，该算法不太适用。但实际上，Delaunay三角网格是具有局部性特点的，本文即是充分利用这个性质，提出一种新的Delaunay三角化策略。

2 Delaunay三角网格的局部性原理

2.1 Delaunay三角网格的局部性特点

传统的B/W增量算法[4-5]是基于整个点集的。但实际上，一个点的三角网格，只与其邻近有限距离内的点相关。我们可以充分利用这个性质，来设计新的生成策略。

2.2 最小Voronoi邻近点集

点集的Delaunay三角网格中，与某个点相关的网格单元，只由邻近局部范围内的有限点决定。我们把这些点，称为该生长点的最小Voronoi邻近点集。

定义1 设S为空间点集，点o∈S，点集M为S的子集。任取p∈S，若p的Voronoi单元和o的Voronoi单元相邻（Vor(p) ∩Vor(o)≠空集），则p∈M；若p的Voronoi单元和o的Voronoi单元不相邻（Vor(p) ∩Vor(o)＝空集），则p∉M。称M为生长点o在点集S中的最小Voronoi邻近点集（Minimum Voronoi Neighbors，简称MVN）。

对于点集中的每个点，它作为生长点时的Voronoi单元的形状只与其最小Voronoi邻近点集相关，而且最小Voronoi 邻近点集中点的个数必然是有限个，分布在该点的局部有限距离内。因此，只要找到其MVN点集，就可以在当前生长点的局部构建Voronoi单元了。只要构造出整个点集的Voronoi图，就可以通过对偶得到点集的最终Delaunay三角网格。

为了从点集中找到生长点o的MVN，我们给出如下准则：

设S为空间点集，点o∈S，M为S的一个子集，T为M中的点与o点组成的Delaunay三角形/四面体网格单元集合。如果M满足：

1）T中的每个三角形/四面体网格单元都以o点为一个顶点。

2）T中的三角形/四面体网格单元组成一个凸包，如果o点不在凸壳上，则T围绕o点构成星形三角形/四面体集合。如果o在凸壳上，则凸壳边界包含以o点为顶点的网格单元。

3）T中任意网格单元的外接圆/球中均不包含S中的其他点，则M为o点的MVN。

上面给出的三条规则，可以作为o点的最小Voronoi邻近点集判定准则。

3 Delaunay局部网格生成算法

3.1 算法流程

3.1.1 最小Voronoi邻近点集搜索策略

为了寻找点集中某个点的最小Voronoi邻近点集，需要对点集中的点进行遍历，并依照空球准则对每个点进行处理。处理过程实际上是一个局部加点过程。

在B/W增量算法中，用于初始化的粗网格是覆盖全体点集的，因此，每一次被加入的点，一定落在现有的某个三角形/四面体网格单元内。但是在本文提出的局部遍历过程中，初始网格并不是覆盖全体的，在现有网格单元达到闭合之前，我们需要额外维护一个结构来记录现有网格的边界。

二维情况下，需要维护的边界记录比较简单。如图1(1)所示，在加入初始的两个点之后，我们得到一个Delaunay三角网格单元，当前网格的边界就是由生长点和两个加入点组成的两条直线，平面被分成4个区域。可以简单的把它们记录为“左点”和“右点”。

在局部网格对生长点达到闭合之前，每次加入新点，需要首先判断它与边界记录的关系，即判断它落在如图所示的A、B、C、D四个区域中的哪一个里：

1）落在A区域：判断干涉情况加点，不更新左/右点，如图1(2)。

2）落在B/C区域：判断干涉情况加点，同时更新左/右点，如图1(3)。

3）落在D区域：判断干涉情况加点，达到闭合，如图1(4)。

达到闭合之后，舍弃左点和右点的记录。

图1 二维条件下的局部加点过程示意

三维情况要比二维要复杂些。首先，三维情况下的边界记录是面，空间被划分成若干区域。其次，边界记录的数目在加点过程中可能变化。

图2展示了三维情况下维护边界记录的过程。在图2(1)中，生长点和初始加入的3个点构成了第一个Delaunay四面体网格单元。此时，边界由它的3个包含生长点的面确定。在加入第四个点时，需要首先判断它落在八个区域中的哪一个，从而采取不同的策略。

随着加点过程的继续，边界记录中的记录数可能发生变化，考虑图2(2)，新点D对于当前边界记录(OAB,OBC,OCA)的位置向量是(-,+,+)，D点的加入会导致OAB所在边界面被两个新的边界面OAD和ODB所取代，边界记录的长度变成4。上述过程一直继续，直到出现图2(5)，这时新加入的点P位于所有边界面的负侧，加入P点之后，当前局部Delaunay四面体网格对生长点o达到了闭合，如图2(6)。此后的加点过程可以用统一的方式来处理。

图2 三维情况下的局部加点过程示意

在当前局部Delaunay四面体网格对生长点o闭合前，维护一个边界面的列表。在每次加入一个新点时，首先计算它和边界列表中每个面的位置关系，得到它的位置向量。对于位置向量中的每个元素：

1）如果为"-"，处理下一个元素。

2）如果为"+"，删除该元素对应的位置向量，并形成两个新的候选边界面。考察该元素在位置向量中的前一个元素，如果为"+"，删除前一个候选面；同时考察该元素在位置向量中的前一个元素，如果为"-"，删除后一个候选面。

三维情况下判断干涉加点的过程也与二维情况类似，可以把二维情况下的局部加点算法当成三维情况下边界记录链表长度固定的一个特例。

给定一个点集，我们可以对其中的每个点，利用上述搜索算法，生成其局部Delaunay网格。这个算法称为Delaunay局部网格生成算法。

3.1.2 引入背景索引网格的辅助

如果处理每个点时，都要遍历整个点集，效率将会十分低下。但Delaunay网格的局部性特点，保证了每一次遍历都会在生长点局部的有限范围内中止。当两个点同时插入一个Delaunay三角网格时，如果它们距离较远，彼此一定互不影响，则称它们为非冲突点。反之，称它们为冲突点。有如下定理[7-8]：

定理1 如果r是网格中所有单元外接圆半径的上界，则距离大于4r的两个点一定是非冲突点。

对于给定的点集，r的距离是一定的。因此，Voronoi邻近点集一定分布在生长点周围的有限区域内。但r的具体值很难预先得到，我们可以采取间接的手段，建立一个二级索引数据结构，人为将点集按照位置关系划分为许多单元块，称为背景索引网格。背景索引网格应该尽可能简单、规则。这里选择最简单的矩形/立方体。

首先将点集中的每个点布置到其所对应的三角网格单元中。之后，对当前生长点o所在的背景网格单元及其一阶邻接单元进行检索，而后再逐级扩展至更高阶邻接单元上的网格单元。当某个网格单元与当前生长点的局部Delaunay三角形/四面体集合中所有三角形/四面体的外接圆/球的交为空，则不再对该网格单元的邻接单元进行检索。由于最小Voronoi邻近点集中的点是有限个，遍历过程会在有限步骤内终止，如图3所示。

采用链表式的方式管理需要遍历的背景索引网格单元，可以保证遍历顺序是以节点为中心环绕进行，逐级向外的，能使得当前局部网格尽早达到闭合。

图3 背景索引网格辅助下的一个点的MVN

3.2 算法描述

Delaunay局部网格生成算法分为两个步骤，一是背景辅助网格的生成，二是针对每个点的局部网格生成。生成背景辅助网格的过程是串行的，对输入点的分布进行分析，创建背景索引网格链表，把点集中的点布置到对应的背景网格单元中，并建立背景网格单元的邻接关系。这里不再赘述。

针对每个点的每个局部网格生成，以三维情况为例，描述如下：

算法3.2每个点的最小邻近点集生成算法

SearchMVNwithBackRect(growPoint,growBox,backGridSet)；

算法输入：矩形背景网格集合backGridSet，生长点growPoint，生长点所在的背景网格单元growBox；

算法输出：生长点growPoint的最小邻近点集；

1）定义backRectList为需要遍历的背景网格单元集合队列，并定义整型动态数组backRectFlags，存放已经加入遍历队列的背景网格单元在背景网格集合backGridSet中的下标。首先将growBox加入到backRectList，并把growRect的下标加入backRectFlags；

2）设i=0，对backRectList中的第i个单元t，对t中的每个点轮流执行加点算法addPointToTetraList。通过返回值表明该点是否被加入当前局部网格；

3）如果t中任意一点被加入当前局部网格，则逐一考察t的邻接单元。在立方体背景网格的条件下，通过计算下标的方式即可得到t的邻接单元；

4）查找每个邻接单元的下标是否存在于backRectFlags中，如果存在，说明该单元已经在遍历队列中。否则，把该单元加入backRectList，下标加入backRectFlags；

5）如果t中没有点被加入当前局部网格，则判断当前局部网格中所有三角形的外接圆是否与t相交，如果存在任意外接圆与t相交，则再次执行步骤3)、4)，逐一考察t的邻接单元，查找其下标是否存在于backRectFlags中，如果存在，说明该单元已经在遍历队列中。否则，把该单元加入backRectList，下标加入backRectFlags;

6）i++，如果i小于backRectList的当前长度，则处理backRectList中的下一个t。

其中UpdateTriToDTSet为局部Delaunay三角形集合三角形更新算法，UpdatePtToDTSet为局部Delaunay三角形集合点更新算法，均指按照Delaunay准则将三角形newT和点p加入到核心点o的邻近Delaunay三角形集合中，不再赘述。

最终结果网格中的每个单元，都在以不同顶点为生长点时重复生成过。在实际计算中，我们可以给每个点设置一个“是否进行了局部网格生成”的标志变量。避免重复加入。

3.3 算法执行效率分析

假设输入点的个数为n，背景索引网格的生成过程较为简单，在对输入点集进行一次遍历后即可完成，时间复杂度为O(n)。

每个点的最小Voronoi邻近点集生成算法的过程则较为复杂。在算法中，每个点的最小Voronoi邻近点集数目有限，且遍历过程都是集中在点集附近的，所以对每个点来说，生成其局部网格的时间复杂度介于O(n)到O(1)之间，我们可以根据网格中的数量关系进行简单估计。

先考虑二维情况下，根据欧拉定理，对于三角网格中的数量关系，可以得到下述结论[7]：

定理2 平面三角网格中，三角网格单元数T、边数E和点数n是同一数量级，即O(T)=O(E)=O(n)。

对于一个点来说，以它为生长点的三角网格单元数是O(1)级的，根据经验，在点集分布较为均匀时，所需遍历到的点的数目一般是最小Voronoi邻近点集大小的5倍左右，由于背景索引网格的存在，这个基数不会太大。我们可以把二维Delaunay三角网格局部生成算法的时间复杂度归纳为O(kn)，其中k是个较大的常数。

在三维情况下，每个点最小Voronoi邻近点集的基数较大。从欧拉定理可知[7]：

定理3 三维四面体网格的边数E、面数F和四面体数T的下限是O(n)，上限是O(n2)。

在三维情况下，对于一个点来说，以它为生长点的四面体网格单元数是下限是O(1)，上限是O(n)。整个算法的下限是O(n)，上限是O(n2)。

不过在大多数实际问题中，三维四面体网格的数目并不会达到上限，该算法在一般情况下仍然具有可行性和实用价值。

4 算法的应用

4.1 地质断裂层面上的PEBI网格生成

在油藏数值模拟计算中，常需要精确描述地质层面的网格，PEBI网格是较好的选择[9]。PEBI网格在定义上等价于Voronoi图，生成PEBI网格的问题也就是确定生长点集并生成Voronoi图的问题。

在实际问题中，所研究的地质层面经常具有复杂的地质断裂构造，传统的二维限定PEBI网格生成算法，是从整体生成点集的Delaunay图，然后通过对偶的方式获得Voronoi图。这种方法对限定条件的处理复杂、繁琐，通用性差。特别是当其应用于逆断层附近时，由于其该局部在垂直方向上的重叠，从而会造成上下盘层面上的点之间相互干扰，生成错误的Voronoi单元。

图4 地质层面模型上的正断层与逆断层

在这种情况下，本文提出的基于索引的局部Delaunay三角化算法可以派上用场。利用搜索最小Voronoi邻近点集的方法在生长点局部构造Voronoi单元可以有效的避免全局搜索所带来的逆断层上下盘点之间相互干扰的问题，只需合理地设置边界附近的生长点，对断裂边界的恢复问题也能够得到很好的解决。图5和图6显示了利用Delaunay对包含逆断层的地质层面进行PEBI网格生成后的结果。

图5 利用本文算法生成的PEBI网格

图6 图5中逆断层附近的PEBI网格

4.2 Delaunay三角化算法的并行化处理

子步骤之间的独立性，是该算法的另一个优势，使得它十分适合进行共享内存、无需通信的并行化处理。适合在同一台机器上进行多核并行，也可以在支持分布式共享内存的机群上运行。

表1显示了三维条件下算法在不同线程数下并行执行所耗的时间，表格中还包含了本算法与传统的B/W增量算法的执行效率的比较。从结果可以看出，随着点集的规模逐渐增大，和并行度的增加，并行算法的优势逐渐显露，其时间复杂度优于传统的未经优化的B/W增量算法。

表1 三维条件下算法并行执行时间

4.3 无线传感器网络中的定位问题

无线传感器网络（wireless sensor networks,WAN）是一个由大量简单、廉价的传感器集成无线通信接口所组成的网络，广泛应用于环境监测、灾难救助等领域。在无线传感器网络的应用中，传感节点自我定位是一个基础问题[10]。

目前的定位算法通常是依靠传感器收集通信或者跳数信息，来计算节点位置的。但是，在一个大型的无线传感器网络中，收集全局的信息是非常困难的，所以采用分布式的定位算法，而且消息的复杂度不能过高。上文提出的基于最小Voronoi邻近点集的Voronoi图生成方法，用于节点定位是很合适的。在定位某个节点S时，通过广播消息收集周围节点的拓扑信息，生成自己的局部Voronoi邻近点集，每个节点不需要知道全局的节点分布信息也可以正确定位。

5 结论与展望

本文以Delaunay三角网格的局部性原理为理论基础，研究并实现一种基于Voronoi最小邻近点集的Delaunay局部网格生成算法，它是基于节点的，子过程具有局部性和独立性，因此可以用于解决一些传统算法不适用的问题，包括地质断裂层面上的PEBI网格生成，Delaunay网格并行化生成，无线传感器网络定位等。

实现点集的网格生成对于整个研究过程还只是一个开端，后续工作主要在以下两个方面进行。一是改善算法本身，使其能够应用于更多的实际问题。二是进一步研究网格细化、质量改善、简单限定条件的恢复、复杂限定条件恢复等问题。

[1]Delaunay B N.Sur la sphere vide.Bulletin [M].Acade′mie des Sciences URSS.1934: 793-800.

[2]Hazlewood C.Approximating constrained tetrahedrizations [J].Computer Aided Geometric Design, 1993, (1): 67-87.

[3]Weatherill N P, Hassan O.Efficient three dimensional Delannay triangulation with automatic point creation and imposed boundary constraints [J].International Journal for Numerical Methods in Engineering, 1994,37(2): 2005-2039.

[4]Bowyer A.Computing Dirichlet tessellations [J].The Computer Journal, 1981, 24(2): 162-166.

[5]Watson D F.Computing the Delaunay tessellation with applications to Voronoi Polytopes [J].The Computer Journal, 1981, 24(2): 167-172.

[6]Thompson J F.Handbook of grid generation [M].CRC Press LLC.1999: 644-661.

[7]闵卫东, 唐 泽.圣三角网格中的数量关系[J].计算机辅助设计与图形学学报, 1996, 8(2): 81-86.

[8]Toomre A.On the distribution of matter within highly flattened galaxies [J].The Astrophysical Journal, 1963,138: 385-392.

[9]孙 迈.面向油藏数值模拟的PEBI网格生成技术研究[D].北京: 北京航空航天大学, 2007.

[10]王继春, 刘黄生.基于Voronoi图的无需测距的无线传感器网络节点定位算法[J].计算机研究与发展,2008, 45(1): 119-125.