《计算几何》学习笔记

前言

这是学堂在线邓老师讲的《计算几何》课程的学习笔记。实际上很多细节，比如点重合、边平行等特殊情况，在课程中没有进行考虑，课程讲的算法只考虑一般性的情况。

（最后两章Geometric Range Search和Windowing Query我还没看）

脚踏实地，勤做笔记，时时温习~

Convex Hull 凸包

课程中讲到的平面凸包算法：

• Extreme Points 极点判断：如果某个顶点在某个三角形内部，那么它不是极点。找出所有的极点构成凸包。$O(n^4)$

  • To-Left测试，判断顶点是否在直线（线段）左边，可以用叉乘实现

• Extreme Edges 极边判断：所有的点都会落在极边的同一侧。找出所有的极边，取顶点构成凸包。$O(n^3)$

• Incremental Construction：对每个新增加的顶点，遍历已有凸包可以判定其是否在凸包内部，如果在外部还可以同时找到两条切线，从而构造处新的凸包。$O(n^2)$

• Jarvis March / Gift Wrapping：首先寻找一个极点（例如 lowest-then-leftmost）$o$ 加入集合，然后对于每个新加入集合的顶点$k$，寻找“最右侧”的顶点$s$，如下图所示，想象从 $k$ 发出的射线逆时针旋转，第一个碰到的顶点就是 $s$，这个过程可以通过遍历顶点做To-Left测试实现。$O(n^2)$

  • 
  • 算法过程类似于“礼品包装”。该算法是Output Sensitive输出敏感的，算法复杂度和输出大小有关。

• Graham Scan：首先寻找LTL点（lowest-then-leftmost）（其实任意一个顶点都可以）作为原点，然后剩下的点按照到原点的极角排序，把原点和极角最小的点（构成极边）压入栈$S$中，然后按照极角顺序遍历剩下的点，假设栈$S$顶部和第二个位置的元素分别为$S_0$和$S_1$，我们判定当前遍历的点是否在$S_1-S_0$左侧，如果是则将此点压入栈顶，如果不是则弹出栈顶。$O(nlogn)$

  • 
  • 正确性证明：可以使用数学归纳法证明对于前$k$个点得到的栈$S$就是前$k$个点的凸包。
  • 简化：可以按照$x$坐标排序，假设原点在$(0,-\infty )$进行扫描可以得到上凸包，假设原点在$(0,\infty )$扫描可得下凸包。

• 分治：首先将点集按照x-axis排序，每次divide的时候从中间将点集分为左右两部分，可以说明这两部分的凸包也会被分开，不会有重合。那么问题就转换为，有了左右两个不相交的凸包，如何合并为一个凸包？合并方法如下图所示，本质上只需要找到上下“切线”即可。

  

  这种算法称为Zig-Zag，做法是先找到连线$rl$，其中$r$是左凸包最右边的点，$l$是右凸包最左边的点，然后比如要找上切线，我们就不断地判断$rl$的端点是否可以往上移，能够往上移的条件是，上面那个点在$rl$的左边（也就是To-Left测试为真），直到不能移动的时候，就找到了上切线。下切线也是同理。寻找的过程传神地描述为“Zig-Zag”。分治算法总体复杂度为$O(nlogn)$。

Geometric Intersection 几何求交

给定两个几何物体，求交问题从易到难分为四个层面：1）判定是否相交；2）相交集合数目；3）求出所有交点；4）构造相交区域几何体。

Preliminary 准备工作

EU（Element Uniqueness）元素唯一性问题：给定一组实数，判定是否存在相同的实数。

• 可以通过排序进行验证，最优复杂度为$O(nlogn)$。

Min-Gap问题：给定一组实数，求出相邻两点间距的最小值。

• 可以通过规约到EU问题，证明其最优复杂度也为$O(nlogn)$。

Max-Gap问题：给定一组实数，求出相邻两点间距的最大值。

• 表面上和Min-Gap问题对称，实际上其最优复杂度可以达到$O(n)$。做法是首先找到最小值和最大值，然后将其区间均匀分为$n$段，得到$n$个桶，将每个元素计算放入对应的桶中，然后过滤出所有有元素的桶，在每个桶中计算出最小值和最大值，问题就等价于求解前一个桶的最大值和后一个桶的最小值的间距，求出这个间距的最大值。
• 
• 课程虽然没有详细解释，但是实际上这个方法的巧妙点在于，最大间距对应的两个相邻的点一定不会存在于同一个桶内，这样就免去了排序的过程，只需要计算桶内的最值即可。因此，只要桶的数目不小于$n$，都可以求解。

IEU问题：EU问题的整数版本，其最优复杂度不变，仍然是$O(nlogn)$。

Segment Intersection 线段求交

Interval Intersection区间求交问题不赘述了。

线段求交问题：平面中给定$n$个线段，检测是否相交（SID，Segment Intersection Detection），或者找出所有相交的线段对（SIR，Segment Intersection Reporting）。

通过将问题规约到EU问题可以确定线段求交问题的时间复杂度下界为$\Omega (nlogn+I)$，其中$I$为交点数目。

课程介绍的求解方法为BO算法（Bentley-Ottmann Algorithm），该算法的思想是，用一条竖直的扫描线（sweepline）从左到右扫描，并处理所有事件event，这里的event对应扫描线遇到线段左端点、右端点或交点。

BO算法要求数据不存在以下特殊情况：

1. 线段端点和交点有相同的x坐标；
2. 线段端点在另一条线段上；
3. 三个或以上的线段相交于同一点；

BO算法用一个优先队列$PQ$（横坐标排序）来存储event，初始时刻$PQ$里面存储了所有线段的端点；用一个平衡二叉树$BBST$来存储当前扫描线和哪些线段相交，$BBST$里面的value为与线段交点的纵坐标。

遇到的三种event以及对应的处理方式如上图所示：

• 对于左端点event：将线段插入$BBST$；检查该线段上下相邻线段与其的交点，将交点插入$PQ$；
• 对于右端点event：将线段从$BBST$删除；检查该线段上下相邻线段的交点，将交点插入$PQ$；
• 对于交点event：在$BBST$中交换交点对应的两条相邻线段；分别检查两条线段与新的相邻线段的交点，将交点插入$PQ$。

事实上，在扫描线抵达交点前的某个时刻，该交点对应的两条线段在$BBST$中总是相邻的，因此它总是会被加入$PQ$中，这就保证了BO算法的正确性。

$PQ$中会有$(2n+I)$个元素，每一步处理event的复杂度为$logn$，因此BO算法的复杂度为$O((n+I)logn)$。最坏情况下BO算法比暴力两两枚举求交还要慢。通常情况下交点较少，效率还是挺高的。

我还有些不懂的地方，BBST中的value是扫描线与线段的交点纵坐标，那么该如何动态更新这个value呢？我找到一个代码实现的例子，它的操作是遇到左端点和交点的event，直接更新整个BBST所有线段的value。但我寻思这么做不就相当于直接重建了BBST，复杂度直接就来到了$n^{2}logn$。

Convex Polygon Intersection 凸多边形求交

约定如果一个polygon在另一个的内部也视为相交。

Intersection Detection 相交检测

使用Dobkin and Kirkpatrick’s Algorithm进行凸多边形的相交检测。

首先将凸多边形进行单调划分（Monotone Partitioning），也就是找到最高和最低的点，划分成左右两条（y-axis）单调链，并且在末端水平延长至无穷远处。

对于凸多边形$P$，其左右链分别为$P_L$和$P_R$，那么多边形$P$和$Q$相交的充要条件就是：$P_L\bigcap Q_R$和$P_R\bigcap Q_L$都不为空。因此原本区域相交的问题就转变为线段链的相交问题。

那么如何求解线段链的相交问题呢？课程讲述的方法是分治，如下图所示

在分治过程中，找到两条链的中间的线段，根据这两个线段的位置关系分类讨论，至少可以排除掉某条链的一半。课程里面关于如何舍弃一半的链讲得不太清晰，也仅仅举了三个例子，其实它舍弃的原因无非就是“如果这半条链和另一条链有交点，那么另外半条链一定有，因此舍弃掉这半条链”，当然这也只是我的理解，我没有去详细地查阅资料看看到底都有哪些情况。总之最后可以证明这种分治的方法检验凸多边形是否相交的时间复杂度为$O(log(m+n))$。

Intersection Construction 交集构造

课程介绍了两种方法来构造两个凸多边形的交集。

第一种方法是边追赶（Edge Chasing），同时遍历两个凸多边形，假设两个多边形当前遍历的边分别为$e$和$f$，那么关键问题就是下一次迭代要前进哪一条边，老师直接放出了判断条件：

这种方法虽然精巧但是很繁琐，课程没有详细讲，具体的方法阐述还是得去找资料，关键词：O’Rourke, Edge Chasing

边追赶算法复杂度为$O(m+n)$ 。

第二种方法是平面扫描（Plane Sweeping），用一条竖直的扫描线从左到右扫描，由于凸多边形的性质，扫描线最多与4条线段相交，因此其每一步的维护都是常数复杂度的。课程也没有详细讲。

Halfplane Intersection Construction 半平面交集构造

给定若干个半平面，求其相交的凸多边形。

课程给出的方法是分治法，首先先说明了对于凸多边形交集构造的算法可以推广到不闭合的凸多边形的情况，然后分治的方法就很简单了，把半平面集平均分成两部分，分别递归求解后，再用凸多边形交集构造算法进行合并。这种分治法的复杂度是$O(nlogn)$。

以前搞ACM的时候还学过一种类似Graham Scan的算法，时间复杂度也是$O(nlogn)$。

Triangulation 三角剖分

Art Gallery Problem 画廊问题

给定一个多边形，问在多边形内部最少放置多少个哨兵（guard），使得可以监控到多边形所有地方？

前人经过证明，不论是在任意地方or顶点or边上放置哨兵作为条件，这个问题都是NP-hard的，甚至于对于正交多边形（每条边都是垂直或水平的），这个问题也是NP-hard的。

不过有一些定理可以约束其数目：

• 对于任意$n$边形，最多需要的哨兵数目为$\lfloor \frac{n}{3}\rfloor$（足够且可达）；

对于这个定理的证明（不考虑空洞，即便似乎有空洞这个定理也成立）非常简单，首先将这个$n$边形三角剖分，然后证明了可以只用3种颜色对顶点染色，使得所有边的两端颜色不同（这个证明过程是：给多边形构造对偶图，这个对偶图必然是树形结构，遍历树形结构时，第一个节点对应三角形的3个顶点，之后的节点都会删去并增加1个顶点，因此可以成功染色），然后说明只需要在某种颜色的顶点放置哨兵，就可以覆盖所有的三角形，因此可以覆盖整个多边形。

• 对于正交多边形（Orthogonal Polygon），最多需要的哨兵数目为$\lfloor \frac{n}{4}\rfloor$（足够且可达）；

Triangulation Existence 三角剖分存在性

三角剖分只对**简单多边形（Simple Polygon）**有效，所谓简单，表示只有相邻的边相交，并构成封闭路径。

三角剖分的适用对象也可推广到内部含有空洞的简单多边形。

如何证明不含空洞的简单多边形的三角剖分的存在性呢？

首先定义一下耳朵（Ear）和嘴巴（Mouth）（起名鬼才）：

前人（G.H.Meisters）已经证明，任何非平凡的简单多边形一定存在至少两个耳朵，也就是双耳定理（Two-Ear Theorem）。因此我们使用“切耳法”就可以进行三角剖分了。

**如何证明含有空洞的简单多边形的三角剖分的存在性呢？**可以用数学归纳法：

首先对于顶点数目$n=3$时，存在三角剖分；

然后对于顶点数目$n>3$的简单多边形$P$，做出归纳假设：对于小于（空洞越少越小，空洞数目相同时，顶点数越少越小）$P$的简单多边形都存在三角剖分。然后证明$P$存在三角剖分：

我们断言多边形$P$必定在外边界上含有一个凸顶点（与相邻两边夹角小于180°），事实上lowest-then-leftmost那个顶点就是凸顶点。

假设我们找到的凸顶点为$J$，与其相邻的两个顶点为$I$和$K$，就得到了可能为内对角线（Internal Diagonal）（整条线完全在多边形内部）的$IK$，于是接下来有两种可能性：

• 如果$IK$是内对角线，不做处理；
• 如果$IK$不是内对角线，我们在$\Delta IJK$内部寻找距离$IK$最远的顶点$M$（必然存在这样的点），那么$JM$就是一条内对角线。

沿着内对角线切开，我们要么会得到两个顶点数更少的简单多边形，要么会使得原来的多边形空洞数目减一，这两种情况都使得多边形”变小了“，因此得证！

Triangulating Monotone Polygons 单调多边形的三角剖分

如果一条链沿着某条直线的投影与原本链上节点的顺序相同，那么称之为单调链；如果某个多边形可以分解为互补的两条链，这两条链沿着同一个方向（一般情况取y轴）单调，那么称之为单调多边形。

对于不带洞的简单多边形$P$，判定其单调性可以在$O(n)$时间内完成，甚至可以同时求出所有单调的轴。

单调多边形的三角剖分算法如下：

（注意这里的多边形默认都是没有空洞的情况，并且默认单调轴是y轴）

首先根据顶点的y轴坐标从上到下合并（merge）排序，按照顺序得到事件队列。处理事件时我们用栈$S$存放待处理的节点。这里先总结出四条性质：

1. $S$中的节点按照高度顺序已排好序；
2. $S$中所有节点在同一条单调链上；
3. $S$中任意3个相连的节点的夹角是凹角（reflex angle）；
4. 下一个事件节点一定和栈顶或者栈底相连。

这四条性质是这个算法执行过程中满足的一组不变性。

算法的初始时刻，“栈的初始状态”是压入了最上面的两个节点。在算法的执行过程中，也就是处理每个事件对应的顶点$c$的时候，可能会遇到3种情况，分别讨论如下：

• same side + reflex，同侧且凹角：

  

  将$c$推入栈顶；

• same side + convex，同侧且凸角：

  

  一直将栈顶、次栈顶和$c$组成的三角形进行切分，然后弹出栈顶，直到变成凹角的情况或者栈中元素不多于1，最后将$c$推入栈顶；

• opposite side，异侧：

  

  一直将栈顶、次栈顶和$c$组成的三角形进行切分，直到栈中元素为1，然后弹出最后一个元素，最后将原始栈顶和$c$推入栈顶；可以看到在这种情况下，最后又回到了和“栈的初始状态”等价的状态。

整个算法的复杂度为$O(n)$。

Others 其它

课程还讲述了如何对简单多边形进行单调性剖分，使其分解为若干个单调多边形，从而完整地解释了三角剖分的算法，总体的时间复杂度为$O(nlogn)$。单调性剖分也是运用扫描线的策略，记录扫描线所在水平直线上存在的梯形和对应的helper，然后在遇到上凸顶点时在相应梯形进行构建内对角线。（这部分暂时没有详细记录了，因为懒）

最后还简单提了一下多面体的四面体剖分，说明了四面体剖分的数目不确定性，以及存在不可剖分的四面体（例如 Schönhardt polyhedron）。

另外还举了一个例子（Seidel’s pylyhedron）：

来说明在三维情况下的画廊问题，可能会在多面体内部存在一个点，在这个点看不到任何顶点，因此就算在所有顶点处放置哨兵，也无法覆盖整个多面体。

Voronoi Diagram 维诺图

直观上来看，Voronoi图是对空间中的一种划分，例如在平面上指定若干个顶点，划分为Voronoi图之后，每个四面体内部的点离相应指定的点最近。

这有许多应用意义，例如图中红色的顶点为食堂，不论你在哪里，你想要去最近的食堂，就可以参照Voronoi图去找到最近的；又例如红色顶点为火源，火势蔓延速度固定，那么绿色的线就是火势交界面。

Terminologies 术语

令 S = { p 1 , . . . , p n } S=\{p_1,...,p_n\} S={p1​,...,pn​}是欧式空间中的一组顶点，每个顶点称为一个site，每个site $p_k$对应的一个空间区域cell记作$Cell(p_k)$，其具有的意义是：在$Cell(p_k)$中的任意一点离$p_k$的距离都要小于到其它site的距离。

对于某个site $p_k$来说，它的cell实际上就是它与其它所有顶点的中垂线所组成的它那一侧的半平面的交集，因此必然是凸的。

Voronoi图中的线段、射线或直线称为Voronoi边（Voronoi edge），边末端的顶点称为Voronoi点（Voronoi vertex）。

Properties 性质

• 任何一个site都有一个非空的cell；
• 空圆性（Empty Disks）：平面上任取一点$x$为圆心，不断放大其半径，直到接触到最近的site，这个site必然是$x$所在的cell的site，这个圆内部一定没有其它的site；（如果$x$在Voronoi边或者Voronoi点上，则会同时接触到多个site）
• 共圆性（Concyclic）：任取一点$x$，其所有最近的site都共圆；

实际上这些性质根据定义来看都很显然的。

另外给出了Voronoi图占用空间的大小：对于$n$个sites的Voronoi图，点数目为$O(2n-4)$，边数目为$O(3n-6)$。证明过程如下：

首先把无界的边统一连接到一个新建的顶点上，这样得到的图是一个平面图。然后假设$v$为顶点数目，$e$为边的数目，$f$为面的数目（也就是cell的数目），$c$为连通块数目，$D$为所有顶点的度数和，就可以得到如下方程组：
$$\begin{array}{rl}planarity& :3v\le D=2e\\ Eule{r}^{\prime }sformula& :v-e+f-c=1\end{array}$$
上面那个公式是因为，每条边的度数都是2，而每个顶点的度数至少是3；下面的公式是著名的平面图的欧拉公式。上式中$f=n$，$c=1$，联立求解就可以得证。

Fary’s Theorem：任何一个可平面化的图（边可能是曲线），都可以等价转换为一个直线平面图。

DCEL 双向链接边列表

DCEL是一种用于存储平面图的结构，便于高效的查询操作。它的核心思想是将每一条边分解为相反方向的两条有向边，并且保证每个面周围的边都是逆时针方向的。DCEL用三张表存储顶点、边和面的信息，分别存储的内容如下：

• 顶点（Vertex）：

  

• 边（Half-Edge）：

• 面（Face）：

根据DCEL的数据存储结构，可以很容易的实现如下功能：

• 根据某一个face遍历其环绕的边；
• 根据某一个顶点遍历其环绕的边；

Hardness 复杂度分析

课程通过规约的方法说明了Voronoi图的时间复杂度下界是$nlogn$。

然后介绍了一些思想：在凸包的求解过程中，如果输入的顶点是sorted的，例如沿着x轴排好序，那么Graham Scan的复杂度就是$O(n)$；但是就算输入的顶点有序，Voronoi图的复杂度仍然不会降低。原因是经过仿射变换后，凸包问题的解不会改变，而Voronoi图会。从本质上讲，Voronoi图依赖于距离，因此在仿射变换这种改变距离的变换下，Voronoi图不能保持不变性。

Construction 构造方法

课程讲述了4中构造Voronoi图的方法。

Naive Construction

根据Voronoi图的性质，对于每个site，我们求出其与其它所有sites的平分线，得到了若干个半平面，那么这个site的cell就是这些半平面的交。时间复杂度为$O(n^{2}logn)$。

Incremental Construction（Bowyer-Watson）

这种构造方法是已有了一个Voronoi图的DCEL结构，现在新增加一个site $p_{n+1}$，该如何进行切分得到新的Voronoi图。

方法的思想就是枚举$p_k$，计算$p_{n+1}$和$p_k$的平分线，如果这个平分线和$Cell(p_k)$相交，那么就找到了$Cell(p_{n+1})$的第一条边界。根据平分线和$Cell(p_k)$的交点，从相交的边可以通过查询对偶边的方式“翻墙”到相邻的cell，并且可以确信，这个相邻的cell也必将被切分，因为这个相邻的cell已经有一个交点了，从这个交点出发沿着平分线可以找到下一个交点……如此下去就可以把$Cell(p_{n+1})$构造出来。

课程中讲述该算法的复杂度时说到，求平分线与凸多边形的交点的复杂度是$O(logn)$，最多可能切分$n$个cells，因此新增加一个site的复杂度为$O(nlogn)$。但是根据我的理解，计算平分线和边的交点的总计算次数，应该是不会超过边的总数的2倍的，而Voronoi图边的数目是线性的，因此就算增加新的site时对每个可能的cell一个一个遍历所有的边的时间也是线性的，所以新增加一个site的复杂度是$O(n)$。

后面讲到的更快的算法都是一次性构建$n$个sites的Voronoi图，因此这个增量法还是有意义的。

Divide-And-Conquer

分治构建算法的策略就是首先将所有site按照x-axis排序，然后每次从中间切开，分别构建之后再合并，因此该算法的核心就在于，给定左右两个DCEL结构的Voronoi图，如何合并为一个。

假设最后合并完成后，新增加的一条连续的**轮廓线（Contour）**如上图中左边的绿色线所示，那么我们可以得到轮廓线的一些性质：

• 轮廓线只会改变与其相邻的cell；
• 轮廓线上的点到左右两边的点集的距离相等，这个距离指的是离点集中的site的最短距离；
• 轮廓线必然沿着y-axis是单调的；
• 轮廓线的长度不会超过$m+n$，其中$m$和$n$分别是左右site的数目；

然后我们考察轮廓线最上面的一条边，这条边必然是左右两个凸包的公切线的中垂线，如下图所示：

首先这条最上面的边一定是无限向上延伸的，假设我们考察这条边无穷远的顶点，它到左右两边的距离相等，因此必然有上面说的结论。轮廓下的最下面的边也是同理。

接下来的计算轮廓线的过程如下：

它的过程就是，首先找到左右的上公切点$l$和$r$，以及轮廓线最上面的一条边$b$，这个时间复杂度是$O(n)$的。然后计算$b$和$Cell(l)$以及$Cell(r)$的交点，假设先和$b$相交的是$Cell(l)$的边界，那么根据相交的那条边，我们可以将$l$更新为与其相邻的下一个site，那么现在$b$就应该更新为新的$l$和原本的$r$的平分线，之后再次计算交点，判断谁先相交，一直向下延伸即可。

（可以结合下面的图理解，每一段黄色的轮廓线$b$都要与左右两个cell分别求交，先相交的更新到下一个site）

这个算法流程还是很清晰的，但是有一些细节需要注意，如下图所示：

图中$L$为左侧的site，$1,2,3$为右侧的site，在某些特殊情况，可能左侧的某个site会对应右侧的非常多的site，如果每次求交都要从$0$号黄色点遍历一遍$Cell(L)$的边，时间复杂度可能会退化到$O(n^2)$。幸运的是黄色的轮廓线同时会构成$Cell(L)$的边，而$Cell(L)$是凸的，因此轮廓线与$Cell(L)$的交点是沿着$Cell(L)$单调的，因此可以保证时间复杂度为$O(n)$。

综上所述，merge两个Voronoi图的复杂度是$O(n)$的，因此用分治法构建Voronoi图的复杂度为$O(nlogn)$。

Plane-Sweep（Fortune）

也可以使用扫描线算法来构造Voronoi图。一般意义的扫描线存在局限性，对于Voronoi图的构造来说，位于扫描线下方的site会影响已经扫描区域的构造结果，因此这个构造算法还要维护扫描线上方的BL（Beach Line，可以形象的理解为海岸线），如上图所示，扫描线和BL将整个平面分为了frozen、unfrozen和unscanned三个部分，分别对应了颜色由深到浅的三个部分。最上面的frozen部分表示形状已经固定，中间的unfrozen部分还会受最下面的unscanned部分的影响。

BL实际上是已经扫描过的sites对应抛物线（Parabola）的下包络面（Lower Envelope），这个抛物线是以site为焦点、扫描线为准线定义出来的。在BL上相邻两段抛物线的交点称为Break Point，可以很容易地证实，Break Point到相邻两个site的距离相等，恰好满足Voronoi图的定义，因此Break Point的移动轨迹恰好可以“勾勒”出Voronoi图的一条边。

虽然从理解上来说扫面线是连续地往下扫的，但是计算机只能处理离散事件，因此我们要找出会影响Voronoi图构造的事件类型。课程总结出了两种事件类型，并介绍了分别应该如何处理：

• Circle Event：对应的情况是“Voronoi图中两条边合为一条边”。这个事件发生在BL上相邻的3个site三点共圆时，因此在扫描线到达这个圆的最下方时，我们来处理这个事件。也就是说每当BL结构发生改变时，我们都要对于改变的site，考察其附近的相邻3个site，将增加的Circle Event加入事件队列中，同时也要从事件队列中删掉变为frozen的site对应的Circle Event。

  处理Circle Event的方法如上图(左)所示。

• Site Event：对应的情况是“Voronoi图中一条边分为两条边”。这个事件发生在当扫面线扫到site的时候，此时需要新增加一段抛物线，并增加一段“悬挂的”Voronoi边。具体的处理方法如上图(右)所示。

扫描线算法的复杂度是$O(nlogn)$。

Delaunay Triangulation 完美三角剖分

与之前的三角剖分不同的是，本节研究内容是对点集的三角剖分。

Point Set Triangulation 点集三角剖分

所谓对点集的三角剖分，就是找到一个极大的对角线集合，对角线之间相互不相交，并且把这个点集的凸包三角剖分。

对于顶点数目为$n$、凸包顶点数目为$h$的点集，其三角剖分获得的对角线数目为$3(n-1)-h$，三角形数目为$2(n-1)-h$。

对于顶点数目$n\gg 1$的点集，其三角剖分方案数目上界为$\Omega (72^{\frac{n}{2}})=\Omega (8.48528^n)$。

构造点集三角剖分的方法是分治法，和凸包算法的分治法很相似，如下图所示：

同样按照x-axis排序后进行分治，在合并的时候通过ZigZag方法往上和往下遍历，中间那条线的轨迹构成新增加的对角线。

Delaunay Triangulation（DT）

Delaunay三角剖分的定义是Voronoi图的对偶图。

对偶图的构造：对于给定平面图$G$，其对偶图$G^{\ast }$的构造方法如下：

• 在$G$的每个face任取一点作为$G^{\ast }$的顶点（这些face应当铺满整个平面，也就是说如果$G$是封闭的，那么外面无限大的地方也算一个face）；
• 对于相邻边为$e$的两个face，这两个face对应的两个顶点连边$e^{\ast }$，$e$和$e^{\ast }$相交；
• 当$e$只为一个face的边界时，这个face对应的顶点构造一个自环，自环和$e$相交。

一个通过Voronoi图的对偶图得到Delaunay三角剖分的例子：

Delaunay三角剖分的时间复杂度为$O(nlogn)$。

需要注意的是，上述通过对偶图构造三角剖分的方法，必须要求原来的Voronoi图中每个顶点的度都是3，如果有特殊情况还需要一些额外处理，课程中没有讲述处理方法。

DT Properties DT的性质

由Voronoi图的性质可以推出Delaunay三角剖分的一些性质，假设三角剖分得到图$G$，那么：

• Empty Circle 空圆性：$G$的每个面的外接圆必然是空圆，即内部没有别的顶点，另外外接圆的圆心是原本Voronoi图的顶点；
• Nearest Neighbor 最近邻：对于site $p$，如果其最近邻（最近的site）为$q$，那么$pq$必然在DT中；由此可知，对于点集$S$，其最近邻图（Nearest Neighbor Graph，NNG）必然是DT图的子图；
• Complexity 复杂度：$d$维空间的DT含有的三角形数目大致是$O(n^{\lfloor \frac{d}{2}\rfloor })$ 。与凸包一样。

对于点集$S$，$|S|=n$，$|CH(S)|=h$，那么其三角剖分含有的夹角数目为$6n-3h-6$。把所有夹角拿出来从小到大排序，就得到某个三角剖分的夹角序列。

Delaunay Triangulation这种三角剖分总是会趋向于避免“狭长三角形”的出现，可以证明，在$S$的所有三角剖分中，$DT(S)$的夹角序列最大，也可以描述为$DT(S)$总是最大化最小的角。

课程讲述了一种简单情况的样例，用到的分析方法主要是圆周角的“等弧对等角”性质。

Proximity Graph

Porximity的意思是接下来介绍的图反映了顶点之间的“邻近性”。

Grabriel Graph

顶点集$S$的Grabriel图的定义为，其中任意一条边$pq$，都满足$|pq{|}^2=\min (|pr{|}^2+|rq{|}^2),r\in S$。其等价的可以描述为：

• 以$pq$为直径的圆中间不存在任何其它顶点；（如果存在，那么以它作为$r$，平方和会更小）
• 若以$S$构造Voronoi图，$Cell(p)$和$Cell(q)$的边界为$e$，那么$pq$与$e$相交；（因为相交，那么以交点为圆心，到$p$和$q$的距离都为最短，自然等价于上面一条描述）

只要有了Voronoi图，就可以以此为基础，以$O(n)$的时间构造出Grabriel图。构造方法就是遍历Voronoi图的每一条边，判断以边为公共边界的相邻site的连线是否与这条边相交，如果相交那么连线属于Grabriel图的边集。

Relative Neighborhood Graph

顶点集$S$的Relative Neighborhood图的定义为，其中任意一条边$pq$，都满足$|pq|=\min (|pr|+|rq|),r\in S$。其等价的可以描述为：

• 以$|pq|$为公共半径，分别以$q$和$p$为圆心的圆的交集部分$lune(p,q)$内部不存在其它顶点；

因为RNG要求更大部分为空，因此它是Grabriel Graph的子图。

只要有了Delaunay Triangulation，就可以在$O(n)$时间内构造出RNG，课程里没有讲方法，但实际上方法很简单，只要在DT中遍历每一个顶点，对于每个顶点遍历其相连的边，找到最短的那条加入RNG即可。

Euclidean Minimum Spanning Tree 欧几里得最小生成树

给定一个点集$S$，其欧几里得最小生成树（EMST）定义为，任意两点之间连边，边权重为欧氏距离，然后得到的最小生成树。

求最小生成树的经典算法Kruskal和Prim似乎最少都要$O(n^2)$或$O(n^{2}logn)$才能实现。

课程首先证明了$EMST(S)\subseteq RNG(S)$，证明方法如下：

首先RNG是之前说到的Relative Neighborhood Graph，然后证明上述包含关系等价于证明，对于任意不属于$RNG(S)$的边，也必然不属于$EMST(S)$。如上图(左)所示，$pq\notin RNG(S)$，因为在相交梭形里面有个顶点$r$，那么假设$pq\in EMST(S)$，我们把$pq$这条边剪掉，整个生成树就被分为左右两部分。不失一般性假设$r$在左边，那么$r$和$p$是联通的，我们新建一条边$rq$，可以得到一个新的生成树，这个生成树的权值更小。因此可得$pq\notin EMST(S)$。

有了这个结论之后，对于点集$S$，我们就可以在$O(nlogn)$时间内先求出其Delaunay Triangulation $DT(S)$，这样图的边数目就降到了$O(n)$，然后再用Kruskal或者Prim算法求最小生成树即可。总体时间复杂度为$O(nlogn)$。

包含关系

综上，我们总结一下目前讲到的图的包含关系：
$$NNG(S)\subseteq EMST(S)\subseteq RNG(S)\subseteq GG(S)\subseteq DT(S)$$

---

（课程还讲了两个实际问题，一个是Euclidean Traveling Salesman Problem，就是旅行商问题并且距离为欧几里得距离；另一个是Minimum Weighted Triangulation，最小权重三角剖分，要求使得三角剖分后边权和最小。两个问题都是NP-hard的，但是ETSP问题可以用DT逼近，使得得到的解不超过最小距离的两倍。）

---

Randomized Incremental Construction 随机化递增构造

构造Delaunay Triangulation的方法有很多，最直接的就是先构造Voronoi Graph，然后得到对偶图，其它的方法也有分治法、扫描线法等等。本小节介绍一种随机化+暴力的方法，Randomized Incremental Construction，简记为RIC，这种方法较为简单。

RIC的递归形式的算法如下：

对于当前的DT，新增加一个site就调用Insert(p)，它首先找到site $p$所在的三角形$\Delta abc$，然后将$p$与$abc$分别连边，然后针对$p$相对的三条边，分别做sTest检测，这个检测的方案如上图(右)所示，假设检测的边为$ab$，那么首先找到$ab$相对的另一个site $x$，如果不存在则返回，如果存在则检测$x$是否在$pab$的外接圆内部，如果在外接圆内部则说明这条原本的边$ab$不合法，就把$ab$删去，同时新增一条边$px$（课程把这个操作称为flip），并且此时新增加两条待检测的边$xa$和$xb$，再次递归做sTest检测。

上面的算法描述有些地方需要详细解释。首先是如何寻找$p$所在的三角形，解决办法是一直维护每个三角形所对应的bucket（桶），那么新增加一个顶点时，我们就需要将原本的一个桶划分为3个桶；flip某条边的时候，就应当把两个桶的元素重新分配。

如何判断顶点在三角形外接圆内部呢？最简单的方法就是计算一个行列式：
$$InCircle(a,b,c,p)=\left|\begin{array}{cccc}{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{array}\right|$$
如果$InCircle(a,b,c,p)>0$，那么点$p$在$\Delta abc$的外接圆内。

总结来说，对于新增加的顶点$p$，我们主要的需要做的事情如下：

• 找到$p$所在的三角形，通过bucket的方法可以在$O(1)$找到；
• sTest检测，通过Backward Analysis（逆向分析），我们可以证明如果顶点$p$是最后插入的顶点，那么它需要sTest的次数等于它的度数，因此插入一个顶点的期望次数为$总度数/顶点数$，是一个常数，因此插入一个顶点的总sTest时间是$O(1)$的；
• 更新bucket，课程证明每一步插入时期望更新bucket的次数是$O(logn)$的；

因此构造Delaunay Triangulation的RIC算法复杂度为$expected-O(nlogn)$，这需要保证输入顶点序列的随机化，才能得到一个期望的复杂度。

Point Location 点定位

我们讨论的是平面点定位问题，也就是对于一个平面的分割，我们给定一个询问点，要求找到这个点在哪个区域内。如下图所示：

一般这种询问都是在线的（Online），要求每次查询的复杂度足够低，一般为$O(logn)$。因此会采用一些算法和数据结构对已有的分割做一些预处理。

接下来介绍几种解决点定位问题的算法：

Slab Method

这种方法首先根据每一个顶点，作一条竖直的直线，这样就可以由$n$个顶点的点集，将平面分割成竖直的$n+1$个slab。

每个slab都是由从上到下有序的梯形（或者三角形）组成，其有序性体现在，我们可以存储slab的线段，对于每个线段可以通过ToLeft来判断点的方位。因为其有序性，我们可以建立搜索树，这样根据y-axis可以在$O(logn)$时间内得知在slab的哪个梯形内。

因此查询某个顶点的过程，就分解为先根据x-axis查询在哪个slab，再根据y-axis查询在哪个梯形。总体查询时间复杂度为$O(logn)$。

但是空间复杂度是非常糟糕的，为$O(n^2)$。

我们可以考虑使用可持久化数据结构（Persistent Structure），可持久化的一些思想我以前在这里总结过。

大致思想就是可以保存历史信息，并且对于树形结构，当我们需要修改一个节点时，可以将那个节点复制出来，然后修改，复制的节点往下的连接关系可以不变，往上的节点也需要复制新的（比如上面的图中就是修改了最下层第二个节点）。使用可持久化数据结构，每次修改一个节点时所耗费的时间和空间都是$O(logn)$的。

再回到上面的Slab Method，我们可以发现，相邻两个slab之间很大程度是一样的，唯一变化的地方就是在节点处删掉左边的线段，新增右边的线段。因此我们可以在可持久化数据结构中实现这一种变化。再考虑到整个图线段数目是线性的，所以在优化过后，空间复杂度和预处理的时间复杂度都是$O(nlogn)$，查询复杂度为$O(logn)$。

Kirkpatrick Structure

根据我听课后的理解，Kirkpatrick结构实际上是一个类似线段树的结构，是一个高度为$O(logn)$的树形结构（严格来说是DAG，有向无环图，但是仍然可以类似树形结构从上往下查找）。如下图所示：

首先Kirkpatrick Sturcture必须满足整个图是一个三角剖分，并且它会在足够远的地方增加3个顶点，使得点集的凸包是一个三角形。然后给三角剖分的每一个三角形做标记，就可以确定它在原来的图中属于哪个区域。

然后它的构造必须要满足，在每一层，选取至少$\frac{n}{K}$个互相独立的顶点（互相独立是指彼此不相连），并且这些互相独立的顶点的度数不超过$D$，这里$K$和$D$是两个常数。然后将这$\frac{n}{K}$个顶点及其连边删去，就得到了$\frac{n}{K}$个空多边形，然后对这些空多边形重新做三角剖分，这样“上层”的三角形就会连接到“下层”不超过$D$个三角形。如此反复直到最后只剩下一个三角形。

观察这个树形结构，最底层有$n$个节点，然后相邻两层中，上层的节点数目是下层的$\frac{K-1}{K}$，因此总的层数为$O(logn)$的；再者，从上层往下搜索时，对于某个三角形其儿子数目不超过$D$个，因此每层搜索时间为$O(1)$，所以每次查询时间为$O(logn)$。另外Kirkpatrick Structure的构造时间复杂度是$O(nlogn)$，而空间复杂度是$O(n)$。

最后就是要证明能够找到这么两个常数$K$和$D$。首先根据欧拉定理可以得到平面图的平均度数$\overline{deg}\le 6$，因此度数至少为$9$的顶点不会超过$\lfloor \frac{2n}{3}\rfloor$个，反过来说，度数不超过$8$的顶点至少有$\lceil \frac{n}{3}\rceil$个。那么我们遍历所有$n$个顶点，每当找到度数不超过$8$的顶点时，我们就将其加入集合$IS$中，然后将这个顶点的相连顶点都标记为不可选取，由于我们标记的顶点最多是选取的顶点的$8$倍，因此我们最终集合$IS$至少可以选取到$\lceil \frac{n}{27}\rceil$个顶点。这样我们就找到了$K=27$，$D=8$。

Kirkpatrick Structure这种结构虽然从理论上来说渐进复杂度很优秀，但实际运用中常系数非常大，难以接受。

Trapezoidal Map 梯形图

Trapezoidal Map（TM，梯形图）方法解决的点定位问题是这样的：给定一个点的坐标，要求输出从这个点往上发出的射线第一个接触到的线段。

这种点定位问题是一个比较普遍的问题，实际上之前讲到的区域定位也可以转换为这一类问题。这个问题还要求线段之间不可相交，但是允许顶点重合（其实也就是平面图的特征）。

TM这种方法和之前讲到的Slab Method在思想上比较相似，都是根据顶点作垂线将平面分割，以得到若干个梯形。但是TM的分割在线段处是不穿透的，如下图所示。另外TM会在外部创建一个bounding box，避免处理无穷远的情况。

另外TM还会创建一个DAG查询结构，在给定了输入点之后，能够快速地查询到它属于哪个梯形。DAG中的数字节点为查询其x-axis在对应顶点的左右；黄色字母节点为查询其在对应线段的上下；蓝色字母节点为终点，对应属于的梯形。

可以看到在上图的例子中，假设我们的输入点在梯形$D$中，那么在查询结构中就会按照$1\to 2\to s\to 3\to t\to D$的顺序找到其属于的区域。

那么如何创建这么一个DAG查询结构（之后简记为SS，Search Structure）呢？

可以使用随机增量构造法（RIC）：首先我们初始化SS，最初始的结构就是只有一个梯形终止节点；然后对于每个新增加的线段$pq$，我们找到$p$所在的梯形T（这一步就是SS的查询功能），根据相应规则updateSS和TM，然后沿着线段$pq$找到与T相邻的下一个梯形，再次update，直到updage完$q$所在的梯形为止。

在上面的算法描述中，有两个实现上需要详细介绍的。其中一个是如何找到在线段$pq$上与梯形T相邻的梯形？我们可以用DCEL结构存储TM，这样寻找相邻以及更新的问题就迎刃而解了。另一个问题是在处理某个梯形T时如何update SS？这需要分三种情况来解决：

• 
• 
• 

至此我们知道了如何根据给定的线段序列构造出TM和SS。需要注意的是，TM与线段序列是无关的，但SS的结构与线段序列（也就是增加线段的先后顺序）是有关的。

如果将线段序列随机化（Randomize），通过一系列的数学证明（我没详细看这一块儿），我们可以得知梯形图的预处理复杂度为$expected-O(nlogn)$，占用空间为$expected-O(n)$，查询时间为$expected-O(logn)$。