摘 要: 平面点集的三角剖分在数值分析以及图形领域,都是极为重要的一项预处理技术。作为一种广泛应用的三角剖分技术,Delaunay三角剖分通过最大化最小角确保接近与规则的三角网和唯一性。本文通过概述 Delaunay 三角剖分的原理,实现了一种增量的 Delaunay 三角剖分构造算法。实验在真实人脸特征点数据和模拟数据上进行,并分别在不同数据规模下进行测试,结果表明了实现算法的有效性。
关键词: Delaunay 三角剖分;Voronoi 图;Delaunay 图; 三角剖分
引言
在数学和计算几何中,对于一个给定一般位置的离散点集 $P$,其 Delaunay 三角剖分记为 $DT(P)$,满足 $P$ 中的任何点都不在 $DT(P)$ 中的任何三角形的外接圆内。Delaunay 三角剖分是三角剖分中所有三角形最小角度的最大化,倾向于避免狭长的三角形出现。由于 Boris Delaunay 从 1934 年开始研究这个课题,因此这种三角剖分是以 Boris Delaunay 的名字进行命名。
对于同一直线上的一组点,不存在 Delaunay 三角剖分;同一圆上的四个或多个点,存在不唯一的 Delaunay 三角剖分。通过考虑外切球,Delaunay 三角剖分的概念可以扩展到三维和更高的维度上,对欧几里得距离以外的度量进行推广是可能的。然而,在这些情况下,Delaunay 三角剖分不能保证存在或唯一。
本文实现 Delaunay 三角剖分算法,实验验证了算法的有效性,并将用在人脸特征点数据的三角剖分上。
Delaunay 三角剖分
本节中,首先由 Voronoi 图的对偶图引入 Delaunay 图的定义,继而引出 Delaunay 三角剖分的定义以及探讨如何构建 Delaunay 三角剖分。
Voronoi 图
平面离散点集 $P$ 的 Voronoi 图是平面的一个子区域划分,其中包含 $n$ 个子区域,分别对应于 $P$ 中的各个基点 $p \in P$ 所对应的子区域,由平面上以 $p$ 为最近基点的所有点组成。
$P$ 的 Voronoi 图,记作 $Vor(P)$。与基点 $p$ 相对应的子区域,称为 $p$ 的 Voronoi 单元,记作 $V(p)$。如 所示,本文主要研究 Voronoi 图的对偶图,记作 $G$。其中,对应于每一个 Voronoi 单元,各有一个节点,若两个单元之间公用一条边,则在这两个单元各自对应的节点之间,联接一条弧。因此,对应于 $Vor(P)$ 中的每一条边,$G$ 中都有一条弧与之对应。正如 所示,在 $G$ 的所有有界面与 $Vor(P)$ 中的所有顶点之间,存在一个一一对应关系。
Delaunay 图
考察 $G$ 的如下直线嵌入。其中,与 Voronoi 单元 $V(p)$ 对应的节点,用点 $p$ 来实现;而联接于 $V(p)$ 和 $V(q)$ 之间的弧,则用线段 $\overline{pq}$ 来实现(如 下图 所示)。这一直线嵌入,称 作 $P$ 的 Delaunay 图,记作 $DG(P)$。 点集的 Delaunay 图总是一个平面图,该直线嵌入中的任何两条边都不会相互跨越。
定理 1. 任何平面点集的 Delaunay 图都是一个平面图。
证明 1. 为证明这一结论,需要利用 Voronoi 图的一个特性,此处将借助 Delaunay 图的概念,进行具体证明。
如 所示,$\overline{p_ip_j}$ 是 Delaunay 图 $DG(P)$ 中的一条边,当且仅当存在这样一个闭圆盘 $C_{ij}$:它的边界经过 $p_i$ 和 $p_j$,而且其中不包含来自 $P$ 的任何其它基点。
由顶点 $p_i$、$p_j$ 以及 $C_{ij}$ 的中心确定的那个三角形,记作 $t_{ij}$。可以看出,$t_{ij}$ 的联接于 $p_i$ 和圆心 $C_{ij}$ 之间那条边,必然完全落在 $V(p_i)$ 内;$p_j$ 也有类似的性质。现在,任取 $DG(P)$ 中的另一条边 $\overline{p_kp_l}$,参照 $C_{ij}$ 和 $t_{ij}$ 的定义方法,也可以定义出圆盘 $C_{kl}$ 以及三角形 $t_{kl}$。
若本定理不成立,即 $\overline{p_ip_j}$ 与 $\overline{p_kp_l}$ 相交。既然 $p_k$ 和 $p_l$ 必然都落在圆盘 $C_{ij}$ 之外,故必然也落在三角形 $t_{ij}$ 之外。这就意味着,在围成 $t_{ij}$ 的三条边中,与 $C_{ij}$ 的圆心相关联的那两条边必有其一与 $\overline{p_kp_l}$ 相交。同理,在围成 $t_{kl}$ 的三条边中,与 $C_{kl}$ 的圆心相关联的那两条边也必有其一与 $\overline{p_ip_j}$ 相交。于是,若考虑在 $t_{ij}$ 的边界上、与 $C_{ij}$ 的圆心相关联的那两条边,以及在 $t_{kl}$ 的边界上、 与 $C_{kl}$ 的圆心相关联的那两条边,则前两条边中的某一条必然与后两条边中的某一条相交。然而,这种情况是不可能的,因为这两条边必然分别完全落在两个不同的 Voronoi 单元之内。 ◻
平面点集 $P$ 的 Delaunay 图,是 Voronoi 图的对偶图的一个直线嵌入。正如我们在此前已经注意到的,$Vor(P)$ 中的每个顶点,都分别对应于 Delaunay 图中的某张面。而在 Delaunay 图中围成每张面的各边,分别对应于 Voronoi 图中与某个 Voronoi 顶点相关联的各条 Voronoi 边(如 所示)。具体而言,在 Vor(P) 中,若基点 $p_1, p_2, p_3, \cdots, p_k$ 各自对应的(共 $k$ 个)Voronoi 单元有一个共同的顶点 $v$,则在 $DG(P)$ 中,与 $v$ 相对应的面 $f$ 的各个顶点,必然就是 $p_1, p_2, p_3, \cdots, p_k$。在这种情况下,点 $p_1, p_2, p_3, \cdots, p_k$ 必然散落在某个以 $v$ 为中心的圆周上。因此可以看出,$f$ 不仅是一个 $k$-多边形,而且它必然还是凸的。
如果 $P$ 中各点是随机分布的,任何四点恰好共圆的可能性就会很小。任何集合,只要其中没有任何四点共圆,我们就称它是处于一般性位置的。若 $P$ 的确处于 一般性位置,则在其对应的 Voronoi 图中,每个顶点的度数必然都是 3。于是,$DG(P)$中的每一张有界面都必然是三角形。我们之所以经常将 $DG(P)$ 称作 Delaunay 三角剖分(Delaunay triangulation),原因正在于此。然而,在此我们还是应该更为谨慎一些,姑且将 $DG(P)$ 称作 $P$ 的 Delaunay 图(Delaunay graph)。至于 Delaunay 三角剖分,我们对其有另一番定义:以 Delaunay 图为基础,通过引入联边 而得到的一个三角剖分。既然 $DG(P)$ 中的每张面都是一个凸集,这样一个三角剖 分就可以很容易地得到。需要注意的是,$P$ 的 Delaunay 三角剖分是唯一确定的,当且仅当 $DG(P)$ 本身已经是一个三角剖分,换言之,$P$ 是处于一般性位置的。
可以借助 Delaunay 图的概念,对关于 Voronoi 图有如下定理:
定理 2. 设 $P$ 为任一平面点集,则在 $P$ 的 Delaunay 图中:
- 三个点 $p_i, p_j, p_k \in P$ 同为某张面的顶点,当且仅当 $p_i, p_j, p_k$ 外接圆的内部不含 $P$ 中的任何点;
- 两个点 $p_i, p_j \in P$ 同时与某条边相关联,当且仅当存在一个闭圆盘 $C$,除了 $p_i$ 和 $p_j$ 落在其边界上之外,该圆盘不包含 $P$ 中其它的任何点。
证明 2. 证明略,详见参考文献 [@cgaa_book] 定理 7.4。 ◻
同时,可以对 定理 1 重新表述为:
定理 3. 设 $P$ 为平面上的任一点集,而 $T$ 为 $P$ 的任一三角剖分。则 $T$ 是 $P$ 的 Delaunay 三角剖分,当且仅当在 $T$ 中每个三角形的外接圆的内部,都不包含 $P$ 中的任何点。
证明 3 证明同 定理 1。 ◻
合法三角剖分
针对 $P$ 的任一三角剖分 $T$,我们来考察其中的某一条边 $e = \overline{p_ip_j}$。如果边 $e$ 不属于 $T$ 中那张无界面的边界,它必然会同时与两个三角形 $\triangle p_ip_jp_k$ 和 $\triangle p_ip_jp_l$ 相关联。如果这两个三角形合起来构成一个凸四边形,那么只要将 $\overline{p_ip_j}$ 从 $T$ 中删去,代之以 $p_kp_l$,我们就可以得到另一个三角剖分 $T’$。如下图所示,这一操作称作边翻转(edge flip)。
对比 $T$ 与 $T’$ 的角度向量,共有六处不同:$A(T)$ 中的六个角度 ${\alpha_1, \alpha_2, \cdots, \alpha_6}$,在 $A(T’)$ 中被换成了 ${\alpha’_1, \alpha’_2, \cdots, \alpha’_6}$。如果
\[\begin{equation*} \min_{1 \leq i \leq 6} \alpha_i < \min_{1 \leq i \leq 6} \alpha'_i, \end{equation*}\]则将 $e = \overline{p_ip_j}$ 称作一条非法边(illegal edge)。换而言之,只要在对某条边进行边翻转操作之后,我们能够使局部的(即对应的六个角度中的)最小角增大,它就必然是一条非法边。由非法边的定义,可以立即得出如下观察结论。
推论 1. 设 $e$ 为三角剖分 $T$ 中的一条非法边。在 $T$ 中对 $e$ 进行边翻转操作之后,设新的三角剖分为 $T’$,则必有 $A(T’) > A(T)$。
不含任何非法边的三角剖分,称作合法三角剖分(legal triangulation)。由上面的观察结论可知,角度最优的三角剖分必然也是合法三角剖分。
为了得到好的三角剖分,也就是要使其对应的角度向量尽可能地大。由合法三角剖分,引入对 Delaunay 三角剖分的角度向量的考察。
定理 4. 设 $P$ 为平面上的任一点集。则 $T$ 是 $P$ 的一个合法三角剖分,当且仅当 $T$ 是 $P$ 的 Delaunay 三角剖分。
证明 4. 证明略,详见参考文献 [@cgaa_book; @dengjunhui] 定理 9.8。 ◻
既然任一角度最优的三角剖分都必合法,故由 可得出推论:$P$ 的任一角度最优的三 角剖分,必是P的一个 Delaunay 三角剖分。当P处于一般性位置时,合法三角剖分是唯一存在的。它就是唯一的那个角度最优的三角剖分,即与 Delaunay 图完全吻合的那个唯一的 Delaunay 三角剖分。
可以证明,在 $P$ 的所有三角剖分中,Delaunay 三角剖分使最小角达到最大。$P$ 的任一角度最优的三角剖分,必是 $P$ 的一个 Delaunay 三角剖分。
构造 Delaunay 三角剖分
$P$ 的 Delaunay 三角剖分的确是一种适宜的三角剖分。其原因在于,Delaunay 三角剖分可使其中的最小角最大化。本节中,主要介绍采用随机增量式算法来直接计算 Delaunay 三角剖分。
首先用一个足够大的三角形将整个点集 $P$ 包围起来,需要引入两个辅助点 $p_{-1}$ 和 $p_{-2}$,它们与 $P$ 中的最高点联合构成的三角形,将包含所有的点。于是,我们需要构造 ${p_{-1}, p_{-2}} \cup P$ 的三角剖分,在得到的三角剖分中只需要删除 $p_{-1}$ 和 $p_{-2}$ 以及与它们关联的各边即为 $P$ 的三角剖分。为此,我们所选择的 $p_{-1}$ 和 $p_{-2}$ 必须相距足够远,才不致于对 $P$ 的 Delaunay 三角剖分中的任何三角形有所影响。尤其必须保证的一点是,它们不能落在 $P$ 中任何三点的外接圆内。
按随机次序逐一引入各点,整个过程中,都要维护并更新一个与当前点集对应的 Delaunay 三角剖分。考虑引入点 $p_r$ 时的情况,如图所示。首先,要在当前的三角剖分中,确定 $p_r$落在哪个三角形内。然后,将 $p_r$ 与该三角形的三个顶点分别联接起来,生 成三条边。倘若 $p_r$ 碰巧落在三角剖分的某条边 $e$ 上,就需要找到与 $e$ 关联的那两个三角形,然后将 $p_r$ 与对顶的那两个顶点分别联接起来,生成两条边。
这样,就得到了一个三角剖分,但它不一定是一个 Delaunay
三角剖分。因为,在引入点 $p_r$ 之后,原来的某些边可能不再合法。为消除这些不合法性,需要针对每一条可能的非法边,调用一次子函数 LEGALIZEEDGE
。这个子函数通过边翻转操作,将所有的非法边转换为合法边。为便于分析,不妨假定集合
$P$ 由 $(n+1)$ 个点组成。
令 $p_0$ 为 $P$ 中字典序最高的点 在 $\mathcal{R}^2$ 中选取点 $p_{-1}$ 和 $p_{-2}$ ,使 $P$ 完全包含于三角形 $\triangle p_0p_{-1}p_{-2}$ 之中将 $T$ 初始化为单独的一个三角形 $\triangle p_0p_{-1}p_{-2}$ 随机地选取 $P/{p_0}$ 中各点的一个次序: $p_1, p_2, \cdots, p_n$ 将点 $p_{-1}, p_{-2}$ 以及与之关联的所有边从 $T$ 中剔除掉。
由 定理 2 可知,一个三角剖分是 Delaunay 三角剖分,当且仅当其中的所有边都是合法的。按照算法 LEGALTRIANGULATION
的原则,不断地对非法边实施翻转操作,直到重新回到一个合法三角剖分。由于原来的任何一条合法边 $\overline{p_ip_j}$,只有在与其相关联的三角形发生变化时,才有可能会成为一条非法边。因此,我们只需检查新生成的那些三角形的各边。这项工作是由子程序 LEGALIZEEDGE
来完成的,它会对有关的各边进行检查,若有必要,则进行边翻转。每翻转一条边之后,可能又会进而使得其它的某些边变得非法。对于所有可能的新非法边,LEGALIZEEDGE
都要递归地调用自己,逐一进行核查。
其中, 第 2 行要检查某条边的合法性。从 LEGALIZEEDGE
的代码可以清楚地看出,由于 $p_r$ 的插入而新生出来的每一条边, 都必然与 $p_r$ 相关联。这一点也可以从中看出:在原有的某些三角形被销毁之后,新生出来的三角形用灰色表示。任何一条边若(从合法)变成非法,则与之相关联的(至多两个)三角形中必有其一发生了变化。所有可能变为非法的边,都必然会接受该算法的检查。也就是说,该算法是正确的。需要指出的是:该算法也不致于陷入无限的死循环。这是因为每经过一次翻转,三角剖分的角度向量总是会单调地增长。
实验与结果
本节中,主要通过使用 Python 程序实现 Delaunay 三角剖分算法,通过仿真实验验证了算法的有效性。
实验环境
实验在配置为 Intel(R) Core(TM) 7-10510U CPU @ 1.80GHz.2.30 GHz 和 16 GB memory 的笔记本上使用 Python 3.7 进行。
结果与分析
为了方便显示, 展示了 24 个点的 Delaunay 三角剖分构造与 Voronoi 图生成步骤的可视化结果。关于完整的生成步骤视频,详见支撑材料。
Delaunay 三角剖分可视化
首先,随机生成包括若干点的集合 $P$,采用 依次迭代,最终生成 $P$ 的一个三角剖分。为了可视化方便,设置点集 $P$ 的大小为 24,生成的三角剖分如 所示。
给出了人脸部特征点的三角剖分实例, 是 Trump 脸部原图, 是对其脸部 68 个特征点的 Delaunay 三角剖分。
Voronoi 图可视化
在 Delaunay 三角剖分的基础上,构造 Voronoi 图。作每一个三角面的外接圆(如 ),这些外接圆的圆心即为 Voronoi 图各边的端点(如 )。连接同一 Delaunay 边为弦的两个外接圆圆心形成 Voronoi 图的各边,分别将 Voronoi 图各个子区域涂色如 和 所示。
算法运行时间
为了统计数据规模对算法运行时间的影响,分别生成 10、100、1000 的基点进行实验,Delaunay 三角剖分算法运行时间实验结果如下表所示。由此可见,随着数据规模的增加,算法运行时间对数线性增加。可以证明,的时间复杂度为 $\mathcal{O}(n \log n)$。
1 |
|
另外, 也报告了不同规模基点 Delaney 三角剖分生成三角网格的数目。可以证明, 由 算法1 生成的三角形,总数目的期望值不超过 $9n + 1$。
小结
本文概述了 Delaunay 三角剖分的原理,详细阐述了构建 Delaunay 三角剖分的算法。分别在模拟数据和真实数据上进行了实验验证,实验结果表明算法的有效性。
参考资料
[1] M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars, Computational Geometry: Algorithms and Applications, 2008. [2] 邓俊辉, 计算几何算法与应用 (中文版), 2011.
[3] https://github.com/jmespadero/pyDelaunay2D
[4] https://blog.csdn.net/weixin_42512684/article/details/106650061
附录
主要程序
1 |
|
人脸特征点提取
1 |
|
可视化程序
1 |
|