PaperRead - Sweep-line algorithm for constrained Delaunay triangulation
号称当时最快的Constrained Delaunay Triangulation Algorithm。
[1] Floriani L D , Puppo E . An on-line algorithm for constrained Delaunay triangulation[J]. Cvgip Graphical Models & Image Processing, 1992, 54(4):290-300.
1. Introduction
从零开始构建三角形，最常用的方法是Delaunay triangulation （DT）。Delaunay三角剖分是一种优化三角剖分，它优化所有三角形的最小内角，保证三角形形状良好。
输出顶点和限制边，在这个基础上进行三角剖分。因此，会使用弱化的Delaunay property。
在图2a中，三角形t2的外接圆包含点P1和P3。P1违反了Delaunay特性，因此交换t1和t2三角形的公共边，形成图2b。这个过程被实现为一个称为合法化的递归过程，稍后将简要讨论。
当t2合法化了之后，我们发现P3位于t2的外接圆内，但是位于受限边e背面。因此，从t2是看不到P3的。受约束的Delaunay三角剖分削弱了该特性，该特性仅对可见顶点有效。
由于三角化的时候，会涉及到很多顶点和限制边，因此需要耗时低的算法。通常CDT算法分为两大类：
（1）第一大类是，将点和边分开来处理。首先针对点进行Delaunay三角化，然后将边插入到DT中。哪些包含插入的边的三角形会被新的三角形替换，这个过程是完全满足weaker Delaunay特性。这些实现中，针对构建第一步DT的时候会采用不同的策略，包括：divide-and-conquer（CHEW, L.P., 1989, Constrained Delaunay triangulations. Algorithmica, 4(1), pp. 97–108），sweep-line（FORTUNE, S.J., 1987, A sweepline algorithm for Voronoi diagrams. Algorithmica, 2, pp.153–174.），incremental algorithms（GUIBAS, L., KNUTH, D. and SHARIR, M., 1992, Randomized incremental construction of Delaunay and Voronoi diagrams. Algorithmica, 7, pp. 381–413），high-dimensional embedding（EDELSBRUNNER, H. and SEIDEL, R., 1986, Voronoi diagrams and arrangements. Discrete and Computational Geometry, 1(1), pp. 25–44.）；
必须确保能够快速定位被插入边穿透的三角形。通过简单的遍历，利用三角形的邻接连接来定位穿透的三角形，以考虑它们相对于约束边的位置。在开始遍历之前，我们需要定位到第一个三角形，通常需要一个高效的数据结构用于获取三角形中的点的位置。还需要考虑到内存需求。然而，当提供合适的数据结构时，插入的边的动态（增量）操作是可能的，这对于许多应用是需要的。
（2）如果数据流不会增加（不是增量的），则同时处理点和边的算法更具有吸引力。通常这样的算法，更快，内存消耗更小，对数据分布的敏感度更低。不幸的是对应的实现会更加复杂。
本篇文章提出了一种新的基于sweep-line的CDT算法。它同时对点和边进行处理。但是，任何特定边的插入都会推迟，直到它被完全扫描。这样，确定被边缘穿透的三角形是高效、简单的，并且不需要任何额外的搜索数据结构。该方法基于sweep-line DT算法（Zˇ ALIK, B., 2005, An efficient sweep-line Delaunay triangulation algorithm. Computer-AidedDesign, 37(10), pp. 1027–1038.）。
2. Background
通过sweep-line可以处理很多几何相关的问题（DE BERG, M., VAN KREVELD, M., OVERMARS, M. and SCHWARZKOPF, O., 1997, Computational Geometry, Algorithms and Applications (Berlin: Springer-Verlag).）。它包含了两个步骤。第一步，对几何元素进行排序（我们的场景中是点和限制边）。第二步，扫描先开始从负无穷到正无穷滑动，在执行所需几何任务和更新数据结构的事件点停止。通常情况下，线已经扫描过的半平面的问题是完全解决了的，而在为扫描过的半平面是未解决的。不幸的是，DT中存在场景不满足上述的情况。Namely, the circumcircle of a triangle whose upper vertex has just been hit by the sweep-line, exceeds the sweep-line and some of the points before the sweep-line could be located in the circle. ？？？？TODO。很长时间认为，线扫描算法不能用于DT。在1987年，Fortune使用了一种非常聪明的方法实现了基于线扫描的方法（FORTUNE, S.J., 1987, A sweepline algorithm for Voronoi diagrams. Algorithmica, 2, pp.153–174.）。
Fortune考虑问题处于三维空间中，针对每个被三角化的点，他放置一个圆锥体并观察圆锥体之间的交点。这些相交的地方，以抛物线的形式存在。如果将这些抛物线投影回平面，能够获取Voronoi图（也就是也可以获取DT）。这个算法非常高效，但是实现很复杂。另一个基于线扫描的算法，结合了线扫描和Lawson's legalization（LAWSON, C.L., 1977, Software for C 1 surface interpolation. In Mathematical Software III, J.R. Rice (Ed.), pp. 161–194 (New York, NY: Academic Press).）。Lawson证明了，任何三角化都可以通过在所有的三角形对上应用empty circle，转变成DT。如果不满足空圆特定，两个三角形的公共边会发生交换，如Fig3所示。
CDT添加了一些限制，正规化也会在两个三角形共享一条受限边的时候停下。扫描线DT算法与Lawson的合法化相结合，显得非常有效，实现起来也相对简单。
3. The algorithm
3.1 Sweep-line DT algorithm using an advancing front
在figure4的例子中，sweep-line已经经过了10个点。他们被两条多线段包围。
(i) 下边界决定了凸包的一部分（图中以点线虚线表示）；
(ii) 上边界折线（图中虚线表示）被认为是advancing front；
所有位于下边界和上边界之间的点通过空圆特性进行三角化。当sweep-line遇到下一个点（如例子中的vertex P11），那么算法将执行以下步骤：
(i) fig4a中所示，将点垂直投影到advancing front。
(ii) 根据投影到的边，构建新的三角形P9P11P8，并且检测空圆特性。如果不满足，主要执行规范化操作（fig4c）。
(iii) 更新advancing front（fig4d）
完整的算法实现，包括特殊情况的处理参见：Zˇ ALIK, B., 2005, An efficient sweep-line Delaunay triangulation algorithm. Computer-AidedDesign, 37(10), pp. 1027–1038.
3.2 CDT algorithm
整个算法的实现包括三个部分
- initialization
- sweeping
- finalization
3.3 Initialization
首先，对所有的输入点按照y轴的进行排序，不管它们是否定义了边。如果顶点有相同的y坐标，那么按照x进行排序。每个点存储了它是不是一条或多条边的upper ending point的信息。如图5所示，P5和P9点，存储了e2和e5边的信息。P7与e1，e3，e4关联。如果边是水平的，那么关联x较小的端点，作为upper point。
第一个三角形在排序和分析顶点之后创建。Zalik（2005）中提到了不同的场景。这使实现变得复杂，是可能出现编程错误的关键因素。针对CDT场景，甚至需要考虑更多的场景。基于这个原因，这里采用了另一种方法，遵循所谓的大三角形或超三角形的想法（被用于 incremental insertion algorithms，（GUIBAS, L., KNUTH, D. and SHARIR, M., 1992, Randomized incremental construction of Delaunay and Voronoi diagrams. Algorithmica, 7, pp. 381–413））。在我们的场景中，只需要两个虚拟点，标记为(P_{-1})和(P_{-2})。他们位于所有顶点的下方，同时他们的x轴的值在输入点最小和最大值的范围之外。公式如下：
初始化的位置类似如下图所示：
对于虚拟点x的位置，只要满足在范围外即可。公式中的(alpha)是任意正数（在我们的实现中采用的是0.3）。
现在剩下的初始化任务都很简单。第一个点(P_0 in V)和虚拟点构建了初始化三角形，如图6。advancing front的初始状态也得到了初始化。它由两条线段组成(P_{-1}P_0)和(P_0P_{-2})。虚拟点的引入，同样简化了sweeping阶段，在后面会提到。在对最后的finalization阶段，所有和虚拟点相关的三角形都会被删除掉。
3.4 Sweeping
遵从sweep-line paradigm，算法在每个点停止，并检测它的状态。这里会有两种情况：
(i) 一个点是孤立的或它是下边缘点（我们将这种场景命名为point event）；
(ii) 一个点是至少一条线段的upper ending point（这种场景称为 edge event）；
显然，点插入是在每个点执行的，而不管其类型如何，但只有当整个边被扫掠时，才执行边插入。
3.4.1 Point event
当扫描线碰到(P_i)，在advancing front边界上的垂直的投影点为(P_Q)，包含该投影点的线段为(P_LP_R)（左边和右边）。Advancing front以结合hash表（用于加速几何检索）双向列表的数据结构实现。虚拟点的引入，只需要考虑两种情况：
(i) middle case。投影点(P_Q)严格的位于advancing front线段(P_LP_R,P_Q e{P_L,P_R})，如图7a所示。然后构建出三角形( riangle P_RP_iP_L)，然后进行规范化操作，并更新advancing front（图7c）；
(ii) left case。投影点(P_Q)恰好是(P_L)。生成两个新的三角形，规范化，更新advancing front。见图8.
但是，形成三角形的过程在这个阶段仍在继续。(P_i)继续和advancing front可见的点形成三角形。这个想法的直接应用会产生许多形状不好的三角形，这些三角形需要在以后通过边交换合法化。为了避免这种情况，提出了两种启发式方法：
(i) 需要对新添加的三角形的边和advancing face的边之间的角度进行判断。如果角度小于(pi/2)，则生成三角形，否则中断三角形的添加（见图9）
(ii) 控制advancing front line segments之间的夹角，减小前表面的起伏（Zalik 2005中这个部分称为basins）。如果角度大于(3pi/4)，basin被三角形填满。TODO???此处的角度是怎么计算的？
每构建一个新的三角形，马上进行规范化，因此Delaunay特性一致是满足的。
3.4.2 Edge event
当sweep-line遇到边的ending point（y比较大的点）的时候才会执行边的插入。实际上，正如前一小节所解释的，首先结束点包含在CDT中。与uppder ending point关联的所有线段以任意顺序一个个的被处理。算法以三个步骤工作：
(i) 定位第一个相交三角形；
(ii) 删除相交的三角形；
(iii) 对删除后的空区域进行三角化。
边e穿过的第一个三角形的位置非常简单。点(P_i)周围的其中一个三角形。三角形可以很快的通过叉乘找到。见图11.
针对第二步，删除相交的三角形，通过三角形邻接链接遍历三角剖分来执行（称作：triangle traversal）。这个方法和Anglada中提到的方法一致。（ANGLADA, M.V., 1997, An improved incremental algorithm for constructing restrictedDelaunay triangulation. Computers and Graphics, 21(2), pp. 215–223）。大致上，在每个访问的三角形中，算法通过确定其相对于约束边的位置来选择下一个三角形进行遍历。每一个遍历到的三角形在之后会被删除，它的顶点存储在独立的列表，分别为above edge e，和below edge e。这样能够获取伪多边形(Pi_U)和(Pi_L)如图12b。
在我们的算法里面，根据edge e相对于advancing face的位置，需要考虑三种情况。第一个情况已经考虑过了，即edge e完全位于advancing front下面。第二个情况是，最简单的情况，当e不和任何三角形相交，如图13。这个时候只会存在lower pseudo-polygon (Pi_L)，它的顶点是从Pi开始沿着advancing front表面直到边的ending point（这个过程被称为advancing front traversal）。
第三个情况是最复杂的情况，edge e边和advancing front表面有一个或一个以上的交点（图14）。但是，这种情况，也可以通过简单的方式解决。算法在上述两种traversal之间进行变换：当traversing advancing front，所考虑的advancing front point一侧是关于约束边计算的。如果位于边的上面，那么将traversal切换成triangle traversal，反之亦然，当点位于边的下面的时候，从triangle切换成advancing front traversal。
如果边e与三角形的一条边完全重合，则可以绕过这三个步骤。三角形边标记为固定，不能交换。
3.5 Finalization
当所有的点和边都被扫过后。三角化依然还没有完成，还需要做额外的两件事：
(i) 移除所有的和虚拟点相关的三角形；
(ii) 添加边界三角形（边界三角形的边应该形成V的凸包）
两个任务时同时解决的。算法从advancing front point (P_{-1})开始。取三个后继的advancing front points，然后计算三角形的带符号的面积。如果是正，决定这个三角形为缺失的三角形，然后进行生成和规范化（图15a）。当抵达第二个虚拟点的时候，这个算法会做细微的更改，因为这个时候我们没有advancing front了（图15b）。算法沿着(P_{-2})到(P_{-1})的虚拟三角形。删除了所考虑的人工三角形，并与以前类似，计算了后三个点的符号面积。在这种情况下，如果有符号的区域是负数，则添加一个新的三角形并使其合法化。
4. Results
本文的算法（SL）和不同的实现进行了比较，包括（Triangle的包里面实现有）：
(i) Lee and Schachter提出的分而治之算法（LS-D&C）（LEE, D.T. and SCHACHTER, B.J., 1980, Two algorithms for constructing a Delaunay triangulation. International Journal of Computer and Information Sciences, 9(3), pp. 219–242.）；
(ii) sweep-line algorithm by Fortune（F-SL）（FORTUNE, S.J., 1987, A sweepline algorithm for Voronoi diagrams. Algorithmica, 2, pp. 153–174.）；
(iii) randomized incremental algorithm by Mucke （M-INC）（MU¨ CKE, E.P., SIAS, I. and ZHU, B., 1996, Fast randomized point location without preprocessing in two- and three-dimensional Delaunay triangulatons. In Proceedings of the 12th annual symposium on computational geometry, pp. 274–283 (Philadelphia, PA: ACM Press).）
