清华大学-计算几何-点集三角化
1、Point Set Triangulation
- 目标：对一系列的点进行三角化(网格化)
当我们有一系列点集时，我们需要对这些点集数据进行结构化。对点集进行结构化的一个重要的方法就是对这些点进行三角化。
这里的对角线分为两类，如上图所示，深绿色的外面的边对应于原来的凸包，另外一些内部的浅绿色的线则是凸包的内对角线。
所谓的点集的三角化，其实就是对这个点集所对应的凸包的区域进行剖分。
- 对于点集的三角化并不唯一
如果一个点集的点数目是n，如果凸包的size记做h，那么任意一个极大的对角线的集合都是由3(n-1)-h条对角线组成；这些对角线剖分出来的三角形的数目为2(n-1)-h。
每新增一个点，大致就增加三条对角线，两个三角形。
边翻转（Edge Flipping）：如上面不唯一的剖分所示，从左边的竖着的那条红色的内对角线变成右边横着的内对角线的这个操作称作边翻转。每进行一次边翻转就会得到一个不同的剖分方案。
- 剖分算法的界
如图所示，若我们已经有两个做好剖分的子集，那么根据第一章凸包算法找上下边界的Zigzag方法，我们就可以得到两个子集连接后的三角化，时间复杂度为O(nlogn)，那么上界为O(nlogn)。
下界：将sorting问题规约到trangulation问题。下界为O(nlogn)
将sorting的一系列输入映射到弧上的点，然后找到输入的最大值和最小值，再在外面引入一个点构成三角化的问题输入。对于这个点与弧上的点的这些对角线就构成了三角化，然后可以通过这个点的DCEL结构枚举各条边及其对应的邻居们，这样就得到了这些点在圆周上的排列的次序，也就得到了sorting的输出。这样就将sorting问题规约到了剖分问题。
因此点集三角化的下界是O(nlogn)
2、Delaunay Triangulation（简称DT）
Delaunay三角化定义：平面上的点集P是一种三角化，使得P中没有点严格处于剖分后中任意一个三角形外接圆的内部。
对偶图：对于一个Voronoi图，若任何两个site之间有一条非空边界，那么这两个site之间连接一条边。
由此得到的对偶图，就是一个三角化。这是由于对于一个一般的Voronoi图（即没有四点共圆的情况出现），它的每一个Vertex的度数都恰好是3，也就是说每一个Vertex都恰好和三条边关联，那么与同一个Vertex关联的三条边就会围成一个三角形。这就说明在对偶图中每一个face都是由三条边围成的三角形。
Delaunay三角化的复杂度是O(nlogn)。因为对偶变换是线性时间，从数学上将只需要把线性个顶点，线性个边对应过去即可，而Voronoi图的构建是O(nlogn)的时间，因此总体的剖分算法的复杂度也就是O(nlogn)
3、Properties
Delaunay三角化与Voronoi图是对偶关系，它们的性质也是有所联系。
（1）、空圆性（Empty Circle）
①对于Delaunay剖分中的任何一张face的外接圆必然是空的。如下图所示的七个圆都是空的。
②在Delaunay三角化中，每一条边之所以会成为一条边，本质上是存在一个空圆以它为弦。
这个是它的充分条件。
左图是剖分图，右图是相应的Voronoi图。
（2）、最近邻性
任何一条连接于最近邻之间的边都会被Delaunay剖分所采用，因为这里头会存在一个以这条边为直径（弦）的空圆。
而按照这种最近邻关系所生成的图又被称为Nearest Neighbor Graph（最近邻图），它是对应的Delaunay剖分的一个子图。NNG是一个有向图，由上图可知最近邻关系不是对等的。
（3）、复杂度
在二维平面中，每增加一个点三角形的数目都会大概增加2，边数增加3。可以说在二维上的Delaunay剖分中是一个线性规模的数据结构。但在三维的情况下这两个指标最多会达到平房的量级，更高维的空间的一般结论也会达到2^d量级。
4、Proximity Graph
- Gabriel Graph（简称GG）
定义：对于给定点集S，若线pq是属于Gabriel Graph的一条边当且仅当
|pq|²=min{|pr|²+|rq|² | r∈S}
几何解释：以pq为直径的圆是空的。如图所示
其他等价定义：若边pq被Gabriel Graph所采用，它的充要条件是①其他的点r与p,q两点形成的角不是钝角②在同样这个点集所对应的Voronoi图中，在p和q之间必然有一段非空的公共边界，且pq的连线必然会和实际上存在但是看不见的那条公共边界存在一个实在的交点。
对于第②点的解释：在Voronoi图中，对应的p和q是site。如果确实p和q中间的连边会和p和q这两个Cell的公共边界相交于一点，那一定是中点x。这个点是Voronoi edge上的一个点，那么以它为中心会存在一个空圆且经过p和q；反过来如果p和q之间的这条边的确会被Gabriel Graph采用并且相应的有一个空圆，那么也说明这个圆心是具有Voronoi edge的必要条件的一个点，即这个圆心就在某一条Voronoi edge上，且这条edge就是介于p和q所对应的Cell之间的那一段公共边界上。
结论：Gabriel Graph中的任何一条边也就属于同样的一个点集所对应的Voronoi图中的一条边，因此Gabriel Graph必然是Delaunay剖分的一个子图。
构造Gabriel Graph的复杂度：在拥有Voronoi图的基础上可以花费O(n)时间来得到它，而构造Voronoi图需要O(nlogn)，因此累计是O(nlogn)的时间即可构造Gabriel Graph。
- Relative Neighborhood Graph（相对邻居图，简称RNG）
定义：|pq| = min{max(|pr|, |rq|) | r∈S}
几何解释：若线pq属于RNG，那么以pq为半径，分别以p和q为圆心作出的两个圆的公共部分必定是空的。如图所示的两个圆的公共的梭形部分。
性质：RNG一定是Gabriel Graph的子图。由上图所示，显然以pq为直径的圆一定是空的。
构造RNG的复杂度：O(nlogn)
- 构造GG或RNG的下界：O(nlogn)
证明：考虑一维的情况，将这个问题规约到sorting。对于sorting输入的一系列点，其实就相当于一维空间中GG和RNG的输入。如图所示，根据GG的定义，相邻两个点的连线必然会被GG所采用。那么构造出来的GG的点的排列是按照紧邻的规则来确定，而这个规则恰好就是对应的排序的结果，因此可以规约到sorting问题上。RNG图亦然。
5、Euclidean Minimum Spanning Tree（欧氏最小生成树，简称EMST）
定义：最小生成树两个节点之间的权重就是它们之间的欧氏距离。
如图所示，给了一系列的点用适当的方式连接起来，使得这个连接的成本达到欧氏意义上的最小。那么这个图实际上是一个隐式的图，因为所有点对之间的距离都是通过欧氏距离隐式的给出。这样，我们的计算输入是一幅完全图。
由于有了欧氏距离这个条件，EMST是我们前面提到的RNG的子图。
证明：思路是任何一条不属于RNG的边也必不属于EMST。
如图所示，图中的pq必不被RNG所采用，假设它被EMST采用，那么将pq剪除后会把EMST图分为不连通的两个部分。而r作为EMST中的一个点，它至少会和p或q中的一个点连通。不失一般性，假设r与p连通，当我们将pq剪除后再连接rq，那么我们会得到一棵新的树，这棵树与之前的EMST的差别就在于pq和rq的权重不同。由于权重的定义是欧氏距离，那么从几何上可知rq的权重必然小于pq的权重。这说明这棵新的树的权重比之前的最小生成树的权重还要小，矛盾。因此不被RNG采用的边一定不被EMST采用。
在有上述结论EMST是RNG的子图后，那么EMST也必然是DT的子图。那对于EMST的构造就可以先花费O(nlogn)的时间构造DT，然后再套用Kruskal或Prim算法即可。
总览：
6、Euclidean Traveling Salesman Problem（简称ETSP）
定义：给出的所有点之间的距离采用欧氏距离，在这样的意义下找到一条环路使得距离和最小。
已证明：ETSP依然是NP-hard问题，但有一些简单的近似方法。我们可以在O(n)时间内从一个已经构造好的EMST中得到一个ETSP的近似解，且近似程度最多不过最优解的两倍。
算法：首先构造点集对应的EMST，然后利用DCEL中类似的方法构造一个初始的环路（如图所示绿色部分）。然后从一个特定的点出发如最左侧的点，沿着初始环路往前走，能往前走就走直到它试图回到以前已经访问过的一个点的时候，就跳过这个点并尝试下一个点，最终会回到起点，并且总的路程（tour）会构成一个环路（如图所示的蓝色虚线）。
在没有优化的情况下，这条tour的总权重|W| = 2 * |EMST| < 2 * |ETST|
7、Minumum Weighted Trangulation（最小权三角化，简称MWT）
定义：对于给定的一系列点集，找到一个三角化，使得它的权重之和最小。
MWT与DT并不是同一码事，如下图，MWT会选择竖直的那条内对角线，而DT会选择水平的内对角线。
MWT问题的是NP-Hard问题。
8、Construction
- 基本操作：edge flip（边翻转）
- 前置工作：
对于每一个三角化所对应的角的数目为6n-3h-6个。
一个指标：将一个三角化中的所有角按照从小到大排列，构成一个向量用来度量和评价这个三角化。如下图所示
若采用AC作为内对角线剖分，则角度序列是 <15,30,60,75,75,105>
若采用BD作为内对角线剖分，则角度序列是<15,30,45,60,75,135>
用字典序比较这两个角度序列，我们可以说AC＞BD
Delaunay三角化存在一个趋势：尽可能使每一个三角形都匀称均衡。如下图所示，上图采用的是DT，下图则是其他方案。因为要保证空圆性，上图剖分的三角形要比下图更为均衡。
从数学上讲，Delaunay三角化会尽可能使得刚才所说的角度序列极大化，也可以说是尽可能使小的角更大，又被称为Maximizing the minimum angle.
解释：由上图所示，我们来考虑ad,ab,bc和cd对应的四个角，由几何知识得：θab＜φab，θbc＜φbc，θcd＜φcd，θda＜φda。因此由DT得到的角相对而言比其他方法得到的要大。
问题：我们总共有6对角，为何只考虑其中的四对角？
个人思考：虽然四边形内角和是360度固定，有角变大自然就有角变小，但是因为我们在进行比较的时候是从小到大按字典序比较，只要偏小的角度比其他方法得到的小角度要大，那么后面较大的角度比其他方法得到的大角度小也无所谓了，毕竟前面就已经比较出大小了。
- 策略总结
（1）确定性算法
①对偶：通过Voronoi图的对偶图可以在O(nlogn)+O(n)=O(nlogn)的时间内得到。
②分而治之
③平面扫描
（2）随机化算法：蛮力算法+随机化处理
9、RIC With Example
RIC：Randomized Incremental Construction
思想：在一个已经成形的现有框架下再增加一个点进行构造。
（1）、点定位
对于给定的一个新的点，判断它处于哪个face中。
（2）、In-Circle test
对于新引入的这个点，判断它的哪个三角化是属于Delaunay剖分的。
当如上图的右图的剖分不属于Delaunay剖分时，将该边进行边翻转（edge flip）翻到左边的那种情况。但当进行了这个操作后，要对那些粉红色的边进行判断是否还是属于Delaunay剖分，这些边又被称为嫌疑边。
（3）、Frontier
对嫌疑边做空圆检测，若不符合，那么继而进行边翻转，再检查其余的嫌疑边。这些嫌疑边实际上是首尾相连构成一个环的，有一个中心G，是一个星形多边形。这个星形多边形又称为Frontier Polygon。
（4）、Convergence
对于这个多边形的角度序列，每做一次边翻转，角度序列都会有所递增，而这个角度序列是有一个上界的，因此这个操作存在一个上界。当抵达上界时算法结束。
10、Randomized Incremental Construction
- 基本操作
①、查找新插入的点在哪个三角形的内部
②、把新引入的点在几何结构上缝合到其中去
③、边翻转
- 递归版本：伪代码如下
1 | Insert (p) |
- 迭代版本：伪代码如下
1 | swapTest(p, a, b, c) |
- In-Circle Test
判断点是否在三角形外接圆内部。
若InCircle(a,b,c,p)>0，则说明这个点落在圆内部；
若InCircle(a,b,c,p)<0，则说明这个点落在圆外部；
若InCircle(a,b,c,p)=0，则说明这个点落在圆周上。
时间复杂度O(1)。
- Point Location
首先一开始的时候将所有的未插入的点进行归类，每一类用一个桶结构（bucket）来进行存储，具体的分类就是按照它们在已有的三角形中分别归属于哪个三角形。
如上图所示，abc同属于一个三角形内，那么这三个点就放在一个桶里。
在进行边翻转的时候会导致三角形结构的变更。如下图，一开始是水平切分成两个三角形，那么abc是一个桶，de是另一个桶；若发生边翻转，那么这两个桶消亡，重新生成新的两个桶，一个桶放置bce，一个桶放置ad。
11、RIC Analysis
- 复杂度分析：主要包含边的变化的复杂度和桶结构修正的复杂度。
后向分析法：在对于已知的结果情况下对前面的时间成本进行分析，站在一个过去完成时，随机的在整个计算过程中的某一个时刻叫停，然后反问这个算法在刚刚过去已经发生的这样一步迭代过程消耗了多少时间。
需要具备的条件：①最终的结果必须是唯一确定的②在中间过程中的任何一个时刻的结果都是确定的③前面的任何一个操作都是概率均等的。
- 边变化：对于每一次顶点的插入，边的变化是常数O(1)的时间。
使用后向分析的条件对应：①假设所有的点都在一般性位置没有退化，那么最终结果是确定的②在中间的任何一个时刻的k个点的临时局部的Delaunay剖分都是确定的③前面k个点的插入从概率上来说每个点是最后插入的概率都是相等的。虽然不知道前k个点究竟是哪个点作为最后一个点插入，但是只需要把每个点所作为最后一个元素插入所对应的时间成本平均起来就可以得到这样一个期望值。
我们将假设是最后一个插入的点称为点p，分析每一个点作为这个点p相应的需要做多少次边翻转以及相应的边操作。
假设这个点G是作为最后一个点即p点被插入，它会引起的修改操作是：①初始插入：即引入三条新增的边；②：边翻转，每一次边翻转会生成一条新的边。对于这个图而言，在第二步的时候会有两次边翻转，所以总共有3+2次边操作。
分析：因为所有的点都可能是作为最后一个点被插入，理论上而言我们应该把所有点的可能统计出来再求平均。具体的来讲，假设上图的G是作为最后一个点被插入，如果我们能把它精确所需要的时间准确的估计出来就足够了，若这个方法能通用就可以施加到所有的点上去。因此现在我们将目光就放在一个点上。
分析这个点G，在每经过一次边的翻转时，都会有一条老的边死掉并且引入一条新的边，同时会有两条新的边成为需要检查的嫌疑边，因此对于这个Frontier而言每做过一次flip后size都是-1+2=1，这说明这个Frontier Polygon忠实的记录了我们此前进行过多少次操作。巧合的是如果算上初始生成的三条边这个数恰好相当于Frontier Polygon的size。注意到点G引起的edge change恰好是5次，而这个Frontier Polygon的size也是5，因此这个polygon最后的size就度量了我们此前如果G是最后一个点被插入的话它所花费的时间。而前面我们说过这个Frontier Polygon是一个星形多边形，G就是它的核，我们就可以直接得出结论说这个点G的度数就度量了花费的时间。
结论：任何一个点作为最后一个点被插入的成本就等于它在当下这张图中所具有的度数。那么平均时间就是当下所有点的度数的平均值。
由欧拉公式知，从渐进的意义上来讲一个平面图的边数大概是顶点数的3倍。而每一条边贡献2度，因此总共不过6n度数，分配到n个点每个点的平均度数不会超过6。也就是说在最后一步插入的过程中，从平均意义上讲所需要执行的边操作无非是O(6)即O(1)，也就是之前所说的常数时间。
- 桶变化（rebucketing）
思路：考虑某一个点，在整个算法的执行过程中这个点从期望的角度讲会参与rebucketing多少次？
结论：O(logn)
假设点pi是第i个点，即截止到它已经插入了i个点，那么剩下n-i个点待插入。而在剩下的n-i个点钟每一个点从期望的意义上讲从概率上都是3/i的概率会参与刚刚过去的最后那次迭代中参与了rebucketing，它们所归属的那个三角形以及相应的bucket有变化。
对于上图，假设在最后的那步迭代中q曾经参与一次rebucketing，它的必要条件是在当下q所属的那个三角形必然是在刚刚最后那步迭代中最新创建的。可以再将这个条件转化为另一个必要条件：如果这个三角形确实是刚才第i次迭代中创建的，那么这个三角形肯定会与最新插入的顶点相关联，也就是说新插入的这个点必然是这个三角形的三个顶点之一，即要么a要么b要么c作为最后一个点引入，才会引起q的rebucketing。因此这个概率是3/i。
总体期望：
因此最后结果是O(nlogn)。
