05点集三角刨分(Delaunary Triangulation)

计算几何-笔记目录

本章主体为点集三角刨分。即给定一个平面点集。求一个建立在点上的最大边集，且边没有交叉。

• Point Set Triangulation
  • 介绍点集三角刨分的定义背景。并且基于凸包的知识直接给出了一个构造方案，给出了上下界。
• Delaunary Triangulation
  • 介绍DelaunaryTriangulation概念即其相关内容，后面称之为DT。先从前面介绍的Voronoi图入手。Voronoid图的对偶图即是一个三角刨分。
• Properties
  • 介绍Delaunary三角刨分的特性。
• Proximity Graph
  • 在本节介绍了两类图，其都为DT的子图。实际上DT是很多重要图的超图。
• Euclidean Minimum Spanning Tree
  • 介绍了EMST问题，即给定平面一堆点集，构造平面上的欧式距离最小生成树结构。使用基于DT的结构即可快速解决。
• Euclidean Traveling Salesman Problem
  • 欧几里得距离旅行商问题。这是一个经典NP-hard问题。本节基于EMST简单给出了一个算法结构。
• Minimum Weighted Triangulation
  • 介绍了最小权重的三角刨分概念，即MWT刨分。说明了其与DT的不同，说明其难度为NP-Hard。
• Construction
  • 主要讲解如何构造DT刨分。是从三角刨分形成的角度来说明的。最后介绍了几个构造策略，引入了随机化算法结构。
• RIC With Example
  • 主要介绍RIC算法。即Randomized Incremental Construction。每一次递增过程都随机添加操作。主要分为两步Insert和Update。Insert插入一个点，Update更新数据结构。
• Randomized Incremental Construction
  • 深入RIC算法，讲述整体实现。细化给出两个关键步骤的数据结构以及操作，圆内测试(In-Circle Test)和点定位(Point Location)。
• RIC Analysis
  • 对RIC算法进行一个时间分析。RIC算法的主要消耗在边交换(edge changes)与重新分桶(rebucketing)上面。

---

点集三角刨分介绍

这里的研究对象是平面上一系列的点。我么要对这些点进行处理，处理就是给这些点持续引入一些边，称之为点集$P$的对角线。直到其构成一个边集，覆盖整个点所形成的凸包，具体来说如下图：

确切来说，对于一个点集合$P$的三角刨分(Point Set Triangulation)，是一组不冲突(non-conflicting)对角线集合，其将$P$所形成的凸包$conv(P)$划分成一系列三角形(2-simplices)。

首先我么可以得到一个点集三角刨分的一些特性。

• 对于点集三角刨分来说，其刨分不唯一，显然交换一个四边形内对角线即为成立的不同刨分。
• 三角刨分用最多$3(n-1)-h$条对角线。
• 三角刨分形成$2(n-1)-h$个三角形。

数学上来说，一个点集三角刨分可能的刨分形式会多达$\Omega{}(72^{n/2})$个可能三角刨分，呈一个指数级增长。

而对于构造一个点集三角刨分，显然是存在一个$O(nlogn)$复杂度的算法。这个方式就可以来自于凸包章节中的zig-zag的合并算法。

假设左右两个已经三角刨分好的凸包。那么通过如图的方式，已知可以将两个凸包合并成一个。同时关注这个合并过程，可以看到其zig-zag过程，其连线的方式恰好构成了一系列的三角形。而形成的这一系列三角形恰好就是填充在两个凸包之间的三角刨分方式。所以合并好的凸包，也是已经三角刨分好的结构。

那么这个问题的下界是多少呢？这个时候就是使用规约的时候了。其实有个最直接从凸包到三角刨分的规约。已知凸包问题的下界是$\Omega(nlogn)$。现在进行这样一个规约：

• 对于一个凸包输入，直接转换为点集三角刨分输入。
• 对于点集三角刨分的输出，可以看出来，其也给出了凸包的边界。通过DCEL结构，遍历最外围的边即可。

老师这里则给出了一个更底层的规约，即$Sorting\leq_N{}Tri$。方式如下图：

• 对于sorting的输入，通过如图的方式，在线性时间内，映射的圆弧上。同时添加一个左侧远点$s$后，作为三角刨分输入。
• 对于三角刨分的输出，其应该有两部分。一个是圆弧点之间的连线，还有一部分是$s$和原来圆弧上的点之间的连线。此时关注跟$s$相连的线，可以注意到这个顺序就是点之间的排序。而这个顺序，通过DCEL结构是容易得到的。

Delaunary Triangulation

Delaunary Triangulation(简称DT)是一种很特别的三角刨分。实际上这种三角刨分的对偶图就是这些点形成的Voronoi图。这种关系可以通过下图表示：

对于一般情况来说，对Voronoi图来说，里面每个Voronoi Vertex点的度都是3。通过对偶图的操作，就可以得出，此时对偶图中的每个面都只有3条边，而这说明正好将对偶图刨分成三角形。但是注意是这是一般情况，对于Voronoi图来说，其也存在度数大于3的情况，这里可以参看CGAA中的描述来详细了解。

总之通过这样一种方式可以看出来，对于三角刨分来说，求出Voronoi也基本上可以求出其三角刨分，这也说明三角刨分的复杂度不超过$O(nlogn)$。而对于DCEL结构来说，其进行对偶图转化是非常方便的操作。只要把其中对应的数据结构，进行对调即可。

现在来看DT有哪些性质，实际上作为Voronoi图的对偶图，其性质也是关联化的，假设一个点集为$S$，其Delaunary Triangulation为$DT(S)$，则$DT(S)$有如下性质：

• Empty Circumcircle：$DT(S)$中每个三角性三点所张成的空圆必定是空的。这个来源就是Voronoi图，而且这个空圆的圆心就是Voronoi Vertex。
• Empty Circle：除了上述之外，还有很多空圆存在于$DT(S)$之中。这些圆只经过两个点，并且以一条边为弦。换句话说，这个边之所以被采用，是因为有一个极大空圆刚好过这条边的两个端点。
• Nearest Neighbors：对于一个点$p$，如果有一个点$q$，使得以$pq$为直径的圆为空的。那么这个边$pq$一定会被DT三角刨分采用。这个是来自于Empty Circle的推论。所以最近邻的$q$一定会被采用。

复杂度分析：

对于二维来说，其边数面数，前面已经看过，都是线性量级。但是对于高维来说则不一定了。对于d维空间$DT(S)$可能包含$O(n^{\lfloor{}n/2\rfloor})$个三角形。

点集上的一些应用

首先看两类DT图的子图来了解图之间的近邻性以及关系

Gabriel Graph

先来看GabrielGraph定义：

直观来讲，就是如果一条边$pq$被$GG(S)$选中，那么以$pq$为直径的圆是空的。然后反观公式，满足$|pq|^2=|pr|^2+|rq|^2$的轨迹形成的是一个圆，如果落在圆内则有$|pq|^2>|pr|^2+|rq|^2$，进而$pq$不会被采用。

如果$pq$被采用，那么空圆存在，也满足之前$DT$的采用条件，所以显然有$GG(S)\subseteq{}DT(S)$

除此之外也可以通过Voronoi图的视角来出发，对于Voronoi图来说，如果$Cell(p),Cell(q)$之间存在边界$e$。那么$pq$一定满足GabrielGraph定义，进而会被$GG(S)$采用。同时这也给出了一个Voronoi图来计算$GG(S)$的方法，遍历所有边，取其两侧的Site来构造即可。

Relative Neighborhood Graph

直观来讲。$RNG(S)$来挑选边的原则是，以$pq$为半径，$p,q$为圆心的两个圆有一个相交的纺锤型区域，在这个区域内没有其他点。回到公式上来说，对于$S$中一点$r$，先取其到$p,q$上距离的最大值。然后再取这些距离的最近值即$pq$。

从图形上来看，一个边$pq$属于$RNG(S)$，那么$pq$为直径的圆为空，所以$pq$必然属于$GG(S)$。进而有如下公式：

$$RNG(S)\subseteq{}GG(S)\subseteq{}DT(S)$$

现在考虑这两类图的求解复杂度，其思路任然是规约：

直接考虑1维情况，规约如下：

• 对于一个Sorting输入可以直接作为$GG(S)$和$RNG(S)$的输入。
• 对于$GG(S)$或者$RNG(S)$的输出，我么根据定义可以直到，其采用的边一定是直线上相邻两个点之间的边。而通过边的顺序，我么就可以得到排序的结果。

Euclidean Minimum Spanning Tree

最小生成树都很熟悉。这里考虑的是欧几里得平面上的最小生成树EMST，或者说此时点间的距离，就取其欧式距离。注意这里的距离是隐式给出的，而且这样形成的一个计算用底图是一个完全图，其边数可能高达$n^2$条。

此时如果我们关注经典算法都不可避免的会去遍历所有边，这个稠密图导致所有算法效率直接$n^2$起步：

• Kruskal对边距离要预排序，自然效率为$O(eloge)=O(2n^2logn)$
• Prim算法，多叉堆优化，那么这个分叉数跟点度数对应，自然效率有$O(e)=O(n^2)$

但是这个欧式图也带来了一定的新特性，例如三角不等式。这个三角不等式的运用就是对于某些点之间边会自然的不被$EMST(S)$采用。这里的说明思路式证明$EMST(S)\subseteq{}RNG(S)$。思路如下图：

对于一条边如果不被$RPG(S)$采用，必然也不会被$EMST(S)$采用。反证如下：

如果$pq$不被采用，说明存在一个点$r$在梭形区域中。如果此时$pq$仍然被$EMST(S)$采用，这一定是一个树的一边。移除该边后，必然会分成两个不连通的联通分支。不失一般性的$r$一定跟其中一个相连，例如左边。那么此时我么可以发现，如果把$rq$连接起来，仍然是一个树，但是树的距离值则小于原先的生成树。所以此时我们知道$EMST(S)\subseteq{}RNG(S)\subseteq{}DT(S)$。

所以我么可以这么来处理这个问题，先把$DT(S)$计算出来，然后再套用之前的算法，即可更快的算出最小生成树。

Euclidean Traveling Salesman

旅行商问题是一个经典问题。这里考虑欧几里得旅行商问题ESTP。所有点之间的距离由欧式距离隐式给出。很可惜的是，ESTP任然是NP-Hard问题。但是我们可以找到一些近似算法来处理这个问题。

这里老师介绍的就是一个简明的近似算法。基于前面的EMST问题可以在$O(n)$时间内，构造出一个长度最多不超过最佳答案2倍长度的解。

这个近似过程如下图：

• 首先生成EMST。
• 然后我们要得到环路，操作就是将EMST中树拆分成两个方向的方向边。可以设想，此时不是在线上走，而是围绕线走一圈。
• 然后我们随便选一个点，然后按照上面有方向的环路前进。对于一个已经遍历过的点，我们只要跳过就可以了，直接走向下一个点即可了。

这样一个结果，确实可以得到我们TSP问题的一个解。而对于这个过程得到的回路权重我们显然有如下结果：

$$|W|=2|EMST|<2|ESTP|$$

因为我们将EMST拆分成了两条回路所以$|W|=2|EMST|$，而EMST严格小于ESTP。因为ESTP是一个环路，移除一条边后就是一颗树，而树中权重最小得就是最小生成树EMST。所以不等式得证。

Minimum Weighted Triangulation

即最小权重三角刨分，找到一个三角刨分，其所有边得权重是最小的。这个问题至今仍是一个开放问题。可以说明其仍然是一个NP-Hard问题。如下图：

构造问题

前面已经看到，对于点集三角刨分来说，其刨分方案存在无数种。但是其包含的三角形个数是给定的，只有$2(n-1)-h$个。首先便是建立在不变量上一个指标来度量三角刨分效果。

这里取的就是所有三角形中所有角度，从小到大得一个排序，来作为度量。并且可以定义在这个度量上的一个字典序来定义度量的一个大小关系。

假设所有角度一共有$6n-3h-6$个，构造一个向量$a=a_1,\cdots,a_{6n-2h-6}$。这个向量显然是跟每一种三角刨分关联的。然后定义其上的大小关系如下$a

而Delaunary三角刨分，则倾向于让这个值更大，或者换句话说，DT刨分，会避免非常尖锐角的出现。这个特性来自于DT刨分的空圆特性(Empty Circumcircle)，因为空圆，说明圆内没有其他点，其他点都在圆外，这说明当前任意三角形边所对的角，都比跟圆外的点所成的角大。

对于构造来说，从前面来看也已经有很多确定算法方式了，例如：

• 通过Voronoi图的对偶图来构造
  • 花费$O(nlogn)+O(n)$时间
• 分而治之
  • L. Guibas and J. Stolfi, 1985
  • R. A. Dwyer, 1987
• 扫面线的方式
  • S. Fortune, 1987

而这里老师则讲述了一个随机化算法，可以做到期望意义上的最优。开启随机算法的篇章！

RIC

RIC即Randomized Incremental Construction

这个算法已经知道了40多年了，但是它的期望复杂度直到1992年才被发现在$O(nlogn)$量级。

整个过程大致可以分为四步来走：

• Point Location：给定一个点，确定其落在Subdivision中的哪个face之中。然后可以把所在的face进行快速三角刨分。
• In-Circle Test：判定新划分的三角形是否属于DT三角刨分。方法就是判定新三角形三顶点构成的圆内是否有其他点。这个可以通过InCircleTest来完成。
• Edge Flipping：如果发现圆内有其他点。则进行边对换(Edge Flipping)操作。每次Edge Flipping都会产生新的三角形刨分，进而有新的检测边。
• Frontier：可以发现，上面每次检测边翻转，都会围绕新加入点操作。这样会形成一条以新加入点为核心的星状多边形。而上面要进行检测的嫌疑边，就是这个多边形的边界。

形式化描述算法步骤如下：

对于一个点$p$：

• 找到包含$p$的三角形$T(a,b,c)$。
• 连接$pa,pb,pc$三条边。
• 然后对$pab，pbc，pca$三个三角形进行交换检测。这实际上是考察对应$ab,bc,ca$三条边是否需要交换。
• 这种考察最关键的点就是$ab$对面那个点$x$，如果$x$存在则进行InCircleTest，否则直接返回，$ab$不需要交换。
  • 如果InCircleTest为true，则交换双边。并且标记新的待检测三角形$pax,pxb$
  • 如果InCircleTest为false。则直接返回。

上图是一个递归版本(Recursive Implementation)。对应递归都可以写成一个迭代版本(Iterative Implementation)。具体思路就是用一个数据结构取存储要检测的边，每次遇到可疑边就加入队列。实际上前面可以看到，对于检测是否合法来说，就是依次检测这些可能不正确的边外接圆是否有点。

现在来看这一系列基础操作的时间消耗：

• Finding：找到点$p$所在三角形。
• Inserting：连接点$p$与三角形三条边。
• Flipping：翻转边。

这些操作在DCEL的基础上，都可以认为是期望$O(1)$的。对于点定位来说，实际可以通过持续维护一个桶结构来加快算法。其结构就是存储每个未插入的点归属于哪个三角形来进行。而在整个过程，会动态修改三角形，所以这个时候要重新把这些点来分类的新三角形对应的桶中。而InCircleTest可以根据补充计算看出来，也是$O(1)$的。

所以整体来讲，对于这个算法复杂度主要有两个消耗点：

• Edge Change：平均而言，插入一个点所期望要进行的边交换次数。
• Rebucketings：上面过程中，将三角形中所有点，重新分桶的时间消耗。

Edge Change

边交换次数可以通过backward analysis得出。所谓backward即从最后落定的结果来看待整个过程的消耗。这个想法非常的贝叶斯。

我们从结果来看，给定点集$S$，最后会得出其DT刨分确定唯一的结果$DT(S)$。并且在每一步插入点的过程中，这个局部的点集也是有唯一结果的。最后在这一步中，对于已经插入的任何一个点都有可能是最后一个加入该结构的点。实际上，这个时候我们可以看到这个结果形成了，但是我们并不知道最后一个点是哪一个，所以我们只能假设均匀随机的插入每一个点。

现在来看随机的每一步插入点之后，这个点带来了多少新边生成或者修改。可以看到对于每个点来说，此时操作的个数等于这个点的度数。所以这个时候考虑其平均操作次数，就是每个点的平均度。而对于平面图$DT(S)$来说其边数不超过$O(3n)$条，所以平局度数不超过$O(6)$。

Rebucketings

直接来算每一步中Rebucketing的点数并不是很好计算。这里变换一个视角，我们关注在某一个特定的点上，我们来看一个特定点$p$会被Rebucketing多少次。

现在来看插入第$i$个点之后，剩余点被Rebucketing的概率是多少。这也是个backward analysis方式。有一个明锐的观察可以看到，如果$q$被Rebucketing，那么$a,b,c$一定是刚刚插入的顶点。所以$q$被Rebucketing的概率为$3/i$。

那么对于Rebucketing的期望，就是每一步中被Rebucketing的点的概率次数和。所以其期望值有：

$$\sum_{i=1}^n\frac{3}{i}(n-i)\leq\sum_{i=1}^n\frac{3}{i}n\leq{}3nlogn$$

所以结合这两步的消耗我们可以得出最后的时间复杂度度是$O(nlogn)$

补充

InCircleTest

在RIC章节部分，老师给出了一个简短的InCircleTest计算方法。类似ToLeft测试一般，只用计算出符号就可以得出点$p$是否在三角形$a,b,c$的外界圆内，但是没有给出证明来源，这里做一个简短补充。

计算公式如下：

$$ \begin{vmatrix} a_x & a_y & a_x^2+a_y^2 & 1 \\\ b_x & b_y & b_x^2+b_y^2 & 1 \\\ c_x & c_y & c_x^2+c_y^2 & 1 \\\ p_x & p_y & p_x^2+p_y^2 & 1 \\\ \end{vmatrix} $$

这里其实计算的思路是，线性化和升维。对于一个二维圆公式我们有如下结构：

$$ (x-t_x)^2+(y-t_y)^2=r^2 \\\ x^2-2xt_x+t_x^2+y^2-2yt_y+t_y^2=r^2 $$

这里我们可以认为点$t$就是这个圆的圆心，显然该公式是不能表示为$t_x,t_y$的线性化的，因为存在$t_x^2,t_y^2$两个项目。但是我们可以引入一个新的变量$R=r^2-t_x^2-t_y^2$则原公式可以做如下变形：

$$ 2xt_x+2yt_y+R=x^2+y^2\\\ \begin{bmatrix} 2x & 2y & 1 \end{bmatrix} \begin{bmatrix} t_x \\\ t_y \\\ R \end{bmatrix}=x^2+y^2 $$

这里注意到对于平面三点共圆相当于三个点都在圆上，三个点$a,b,c$都满足方程。所以有如下等式：

$$ \begin{bmatrix} a_x & a_y & 1 \\\ b_x & b_y & 1 \\\ c_x & c_y & 1 \\\ \end{bmatrix} \begin{bmatrix} 2t_x \\\ 2t_y \\\ R \end{bmatrix}= \begin{bmatrix} a_x^2+ a_y^2 \\\ b_x^2+ b_y^2 \\\ c_x^2+ c_y^2 \\\ \end{bmatrix} $$

这个线性方程组可以通过Cramer’s Rule来求解出$t_x,t_y,R$，而$R$实际上包含了半径数值。

同时这个方法可以扩展到任意维度，三维就是要四个点，三个坐标。

现在我们考虑一个新的点$p$，这个点是平面给上任意一个点。并且设$p$到圆心的距离为$u$。我们自然有：

$$ (p_x-t_x)^2+(p_y-t_y)^2=u^2 $$

我们现在考察如下结构

$$ \begin{bmatrix} a_x & a_y & 1 & 0\\\ b_x & b_y & 1 & 0\\\ c_x & c_y & 1 & 0\\\ p_x & p_y & 1 & 1\\\ \end{bmatrix} \begin{bmatrix} 2t_x \\\ 2t_y \\\ R \\\ U \end{bmatrix}= \begin{bmatrix} a_x^2+ a_y^2 \\\ b_x^2+ b_y^2 \\\ c_x^2+ c_y^2 \\\ p_x^2+ p_y^2\\\ \end{bmatrix} ,U=u^2-r^2 $$

可以验证最后一行是成立的。将该方程表示为$AT=B$，$A$为左侧矩阵，$T$为中间向量，$B$为右侧向量。注意根据Cramer’s Rule，我们可以看到U的值可以表示为：

$$ U = \frac{1}{|A|} \begin{vmatrix} a_x & a_y & 1 & a_x^2+a_y^2\\\ b_x & b_y & 1 & b_x^2+b_y^2\\\ c_x & c_y & 1 & c_x^2+c_y^2\\\ p_x & p_y & 1 & p_x^2+p_y^2\\\ \end{vmatrix} $$

注意其中$|A|$为三角形$a,b,c$的按顺序有向面积。所以此时U的符号取决于后面行列式，注意的是$U=u^2-r^2$。刚好是到目标点的距离，其符号说明了在圈外还是圈内。当行列式值大于0说明在圈外，反之说明在圈内。

此外值得注意的一点就是，我们考虑三个三维空间中的点$a,b,c$。三个点如果不共线，则可以确定一个平面这个平面可以表述如下：对于平面内的一点$p=(p_x,p_y,p_z)$有

$$ \vec{ap}\cdot(\vec{ab}\times\vec{ac})=0 $$

语义也非常明确即，面上$p$的点到$a$的向量都与面上的法向量$\vec{ab}\times\vec{ac}$垂直。在三维空间中如果用行列式展开有：

$$ \begin{vmatrix} p_x-a_x & p_y-a_y & p_z-a_z \\\ b_x-a_x & b_y-a_y & b_z-a_z \\\ c_x-a_x & c_y-a_y & b_z-a_z \\\ \end{vmatrix}= \begin{vmatrix} 1 & 0 & 0 & 0 \\\ 0 & p_x-a_x & p_y-a_y & p_z-a_z \\\ 0 & b_x-a_x & b_y-a_y & b_z-a_z \\\ 0 & c_x-a_x & c_y-a_y & b_z-a_z \\\ \end{vmatrix}= \begin{vmatrix} 1 & a_x & a_y & a_z \\\ 1 & p_x & p_y & p_z \\\ 1 & b_x & b_y & b_z \\\ 1 & c_x & c_y & b_z \\\ \end{vmatrix} $$

如果该值大于0说明点$p$在平面法线上方，反之则在下方。对比这个公式跟我们上面的计算公式我们不难发现，对于平面上的点$a,b,c$，可以认为其就是做了一个映射，附加一个z方向值有$(a_x,a_y)\rightarrow(a_x,a_y,a_x^2+a_y^2)$。此时$p$在圈内的问题转化为了在平面上下方的问题。

令$z=x^2+y^2$我们可以看到这个三维图像实际上是一个抛物面。而我们上面所作的映射使得，平面所有点都落在了这个抛物面上面。落在这个抛物面上的点，即可形成一个平面。根据我们上面描述在该平面下方的会落在圆内，而上方的会落在园外。现在我们关注这个平面跟抛物线的截面，该截面会把整个抛物面分成两部分。两部分都可被平面上的点映射过来，可以看到，其实圆内的点就是其中一个有限截面在平面上的投影部分。如下图