Delaunay三角剖分算法介绍

一.什么是Delaunay三角剖分

从事数值计算相关领域的读者，相信或多或少都听说过“三角剖分”这个概念。在诸如有限元仿真，光线追踪渲染等计算当中，都需要把几何模型转化为三角网格数据，即“三角网格生成”。在这个过程中，三角剖分是十分重要的一步。下面我们就具体了解一下三角剖分的概念。

三角剖分:

顾名思义，三角剖分(triangulation)就是对给定的平面点集，生成三角形集合的过程。考虑平面点集P={ P1...Pn },我们希望得到三角形集合T=在{ t1...tm },满足：

a）所有三角形的端点恰好构成集合P。

b）任意两个三角形的边不相交（要么重合，要么没有交点）。

c）所有三角形的合集构成P的凸包（convex hull)。

注：有时我们也称三角剖分为“三角化”，二者表示同一概念。而后者有时也指从多边形或者多面体生成三角网格的过程，即“三角网格生成”。为避免混淆，本文中不使用“三角化”这一词语，而专用“三角剖分”和“三角网格生成”来指代两个不同的过程。

一般来说给定一个点集，往往存在不止一个三角剖分。以下图为例，我们只要把第一个三角剖分中的两个相邻三角形的公共边做一次翻转，就可以得到一个新的三角剖分。

图一：对于给定点集的两种三角剖分

由于给定点集的三角剖分不唯一，我们希望从中挑选出一个“最优”的三角剖分。那么何为最优呢？这涉及到对三角形“质量”（quality）的评定。一般在数值计算中，我们不希望三角形过于狭长，也就是说越接近等边三角形越好。以下是几种常见的质量评定标准：

a) 最小角(minimum angle)：即所有三角形的内角当中最小的角。

b) 纵横比(aspect ratio)：三角形最短边与最长边的比例。

c) 半径比(radius ratio)：三角形内接圆半径的两倍与外接圆半径的比例。

Delaunay三角剖分

所有三角形的外接圆均满足空圆性质的三角剖分，称为一个Delaunay三角剖分。

空圆性质

即一个三角形（或边）的外接圆范围内（边界除外），不包含点集P中的任何顶点。

对于上面的例子来说，左图的三角剖分是Delaunay的，任何一个三角形的外接圆内部都不包含点集中的顶点；而右图的三角剖分不是Delaunay的，因为下面的两个三角形的外接圆内部都包含了顶点。

图二：左图的剖分满足空圆性质，而右图的剖分不满足

二.Delaunay三角剖分的相关理论

这部分简单介绍一下Delaunay三角剖分的相关理论，是为了更好的理解Delaunay三角剖分的最优性质和构造算法。对理论不感兴趣的小伙伴可以直接跳到结论部分。先介绍一下Delaunay三角形和（局部）Delaunay边的概念。

a) 对三角剖分中的一个三角形，如果其外接圆满足空圆性质，称为一个Delaunay三角形。

b) 对三角剖分中的一条边，若存在一个外接圆满足空圆性质，称为一条Delaunay边。（一条边的外接圆可以是经过这条边两个顶点的任意圆，不唯一。）

c）对三角剖分中的一条属于两个三角形 t1,t2 的边e，若存在一个外接圆不包含三角形 t1,t2中的任何顶点称为局部Delaunay边。若一条边只属于一个三角形，也属于局部Delaunay边。

由定义可知，Delaunay边一定是局部Delaunay的，但反过来未必。以下图为例，左侧三角剖分中的蓝色边是局部Delaunay边，因为可以找到一个外接圆不包含A,B,C,D中任何一个点。但它并非一条Delaunay边，因为它的外接圆总是包含其他两个顶点。而右侧三角剖分中的蓝色边不是局部Delaunay边，因为它的任意一个外接圆都至少包含A,B,C,D中的一个点。

3.左图的边BD是局部Delaunay边，而右图的边AC不是局部Delaunay边。

读者们可以再去回看一下图2，其中左侧三角剖分中的蓝色边既是Delaunay边也是局部Delaunay边，而右侧三角剖分中的蓝色边既非Delaunay也非局部Delaunay。可见判断“局部”二字意味着它只与相邻的两个三角形有关。

以下结论被称为Delaunay引理

对一个三角剖分T，以下三个命题相互等价

a) T中所有三角形均为Delaunay三角形。（Delaunay三角剖分的定义）

b) T中所有三角形的边均为Delaunay边。

c) T中所有三角形的边均为局部Delaunay边。

此处略去证明，感兴趣的读者可以阅读Jonathan Schewchuk关于Delaunay triangulation的讲义。

基于局部Delaunay边的概念，可以定义下面的翻转边算法（flip algorithm）。

定理：对三角剖分T中的一条边e，若它不是局部Delaunay的，则可以被翻转成为一条局部Delaunay边e’。

边的翻转：假设一条边e属于两个三角形 t1,t2 ，这两个三角形的合集所构成的四边形有两条对角线，e为其中一条。现在用另一条对角线e’替换e，得到两个新的三角形 t1′,t2′ ，并用这两个新的三角形替换原三角剖分中的 t1,t2 ,得到新的三角剖分T’。这个操作称为边的翻转。

图2和图3从右到左都是一次翻转边操作。可以看到，每进行一次翻转边操作时，都把三角剖分“局部”地改善了。更具体的讲，有以下结论：

a) 最小角增大：翻转边e后的两个相邻三角形的6个内角中的最小角大于等于原先的两个三角形6个内角中的最小角。

b) 外接圆半径缩小：翻转边e后的两个相邻三角形的外接圆半径最大值小于等于原先的两个三角形外接圆半径的最大值。

c) 四点共圆的情况：若一条边e的两个相邻三角形的顶点A,B,C,D四点共圆，则无论翻转与否，这条边都是局部Delaunay边。

图4.经过一次翻转边操作后，两个相邻三角形的最小角增大了，而外接圆半径缩小了，这说明这次翻转边操作改善了三角剖分的质量。

根据前面介绍的Delaunay引理，对任意一个三角剖分T，只要持续进行上述的翻转边操作，最终总可以转化为一个Delaunay三角剖分。这实际上提供了一种构造Delaunay三角剖分的思路（Lawson算法，详见第三部分）。而根据上面提到的关于翻转边的结论，可以得到下面三条关于Delaunay三角剖分的最优性质。

Delaunay三角剖分的最优性质

a) 最大化最小角性质：

在给定点集P的所有三角剖分当中，Delaunay三角剖分得到的最小角（所有三角形的内角中的最小值）是最大的。

b) 最小化外接圆性质：

在给定点集P的所有三角剖分当中，Delaunay三角剖分得到的最大外接圆半径（所有三角形的外接圆半径中的最大值）是最小的。

c）若点集P中任意四点不共圆，则存在唯一的Delaunay三角剖分T。若点集P中四点A,B,C,D共圆，且△ABC,△BCD属于Delaunay三角剖分T，那么将边BC翻转后得到的三角剖分T’(包含△ABD,△ACD)同样是一个Delaunay三角剖分。

其中最大化最小角性质正是我们在第一部分中提到的，三角形质量的度量标准之一。由于具备这些最优性质，在许多三角剖分的应用场景中，Delaunay三角剖分总是会被优先考虑。

三.Delaunay三角剖分的构造算法

大部分Delaunay三角剖分的构造算法都是通过逐点插入来实现的，这里我们介绍两种： Lawson算法和Bowyer-Watson算法。

前面提到的翻转边算法已经可以帮助我们消除“坏边”（指非局部Delaunay的边）。而每当插入一个新的顶点时，将其与附近的旧顶点连接会加入一些新的边，我们只要用翻转边算法来修正其中的“坏边”，就可以将其重新“修复”为Delaunay三角剖分。顺着这个思路我们就会得到下面的Lawson算法。

Lawson算法：

a) 先计算点集P的包围盒(bounding box)，将包围盒的四个顶点加入P中得到P’。根据包围盒生成两个超三角形(super triangles)，构成初始三角剖分 T0 。由于只包含两个直角三角形， T0 是（包围盒四个顶点的）一个Delaunay三角剖分。

b) 将点集P中的顶点逐一插入现有的三角剖分 Ti 中，并进行如下调整：

➤设插入的顶点v位于三角形t中，将v与三角形的三个顶点连接，使t分裂为3个三角形 t1,t2,t3。

➤分别检查 t1,t2,t3 是否满足空圆性质，若不满足则进行翻转边操作，直到没有坏边为止。此时得到一个包含顶点v的新Delaunay三角剖分。

c) 当最后一个顶点插入到三角剖分中，并且完成所有翻转边操作后，我们得到了点集P’的一个Delaunay三角剖分。现在删除第一步中加入的包围盒的四个顶点，并且去除所有与它们连接的三角形，则剩下的三角形就构成点集P的Delaunay三角剖分T。

在Lawson算法的基础上，可以再做一点改进：在插入新的顶点之后，不要立刻与旧顶点连接生成新的三角形，再去逐一翻转坏边。而是先去寻找哪些边是需要翻转的，直接将这些边删掉形成一个“空穴”(cavity)，然后在空穴当中连接新的边。这就是Bowyer-Watson算法的思路。下面来详细描述一下调整的过程。

Bowyer-Watson算法

a) 同Lawson算法

b) 将点集P中的顶点逐一插入现有的三角剖分中，并进行如下调整：

➤在现有三角剖分中，所有外接圆包含顶点v的三角形的合集构成一个“星形多边形”(star shaped polygon)。星形多边形的含义是多边形的任何一个顶点到v的连线都在多边形内部。

➤对于上述星形多边形，将其内部的三角形全部删除，形成一个“空穴”。将空穴边界的顶点与新添加的顶点v连接得到新的三角形，替代剖分中被删除的三角形，此时得到一个包含顶点v的新Delaunay三角剖分。

c) 同Lawson算法

当然，上述算法的正确性是有一些理论来证明的，这里就不详细探讨了，感兴趣的读者可以查看相关讲义。下图展示了Bowyer-Watson算法的调整过程。

图5.左图：原Delaunay三角剖分及新插入的顶点v，中间：所有外接圆包含顶点v的三角形构成一个星形多边形，右图：将星形多边形内部的三角形删除，并将v与空穴顶点重连后得到的新三角剖分，这是一个包含顶点v的Delaunay三角剖分。

由于Delaunay三角剖分的唯一性，两种算法在效果上没什么区别。但Bowyer-Watson算法因为缩短了搜索坏边的过程，效率更高一些，所以在实际应用中我们往往会使用Bowyer-Watson算法。

四.基于Delaunay三角剖分的网格生成算法

三角网格生成(mesh generation)算法有很多种，这里我们介绍基于Delaunay三角剖分的算法（其他还有波前法，八叉树方法等等）。

图6.左图的几何模型由点、线、面等几何元素组成，右图为基于Delaunay三角剖分生成三角网格，全部由三角形组成。

在实际应用中，三角网格生成的对象大多是三维空间中的几何曲面。常见的几何曲面，如平面，圆柱面，球面，B样条曲面，都是在二维空间参数化。因此实际上我们是在二维的参数空间上生成三角网格，最后再根据参数表达式投射到三维空间中。这里我们就简述一下对于给定的二维参数区域，如何基于Delaunay三角剖分去生成一个三角网格。

生成三角网格与构造三角剖分的主要区别在于，点集P并非给定，需要自己去生成。而一旦确定了点集P，剩下的工作就是构造三角剖分了。实际上，生成点集P与构造三角剖分的过程是同时进行的。其过程如下：

Delaunay平面区域网格生成算法

a) 计算区域的包围盒，生成两个超三角形作为初始三角网格。根据边界曲线生成边界点。

b) 将边界点逐一插入到三角网格中。每一步插入后，用Bowyer-Watson算法得到新的Delaunay三角网格。

c) 将第一步插入的辅助顶点删除（同时删除与其连接的三角形），得到关于全部边界点的Delaunay三角网格（亦即边界点的Delaunay三角剖分）。

d) 在区域内部的三角形，根据一定的条件，将重心（centroid）插入三角网格，并根据Bowyer-Watson算法调整，得到新的Delaunay三角剖分。直到不再有需要插入的重心点，此时得到区域的完整Delaunay三角网格。

下面用一个简单的例子来直观地展示这个过程：

图7.从上至下，从左至右依次为：1)参数区域 2)由包围盒生成的初始三角网格及边界曲线上生成的边界点3)将边界点逐一插入并调整后得到的Delaunay三角网格 4)去除辅助顶点及相连的三角形得到的关于边界点的Delaunay三角网格5)加入内部点之后得到的完整Delaunay三角网格。

每个曲面都可以用上述算法在参数区域上生成其Delaunay三角网格，然后投射到三维空间中。然而实际的模型往往包含多个曲面，因此要保证边界的一致性。下面我们介绍完整算法：

三维模型的Delaunay三角网格生成算法

a) 遍历模型的所有边（edge），生成边界点。

b) 遍历模型的所有面（face），收集其边界点在参数区域的坐标。然后在参数区域上按照上面的算法生成Delaunay三角网格。

c) 将每个面在参数区域上生成的内部点和三角形投射到三维空间，添加到三角网格中。

这里要注意的是，每个面上的点需要两份编号，一个局部（local）编号用于参数区域上的三角网格，及一个全局（global）编号用于三维空间的整体三角网格。在网格生成过程中需要时刻保持这两种编号的对应关系。

五.结果展示

以下结果为齿轮和斯坦福兔子的几何模型以及shonMesh中生成的Delaunay三角网格。目前相关功能还在开发中，欢迎对网格生成算法感兴趣的读者与作者联系。感谢阅读！