PaperRead - Sweep-line algorithm for constrained Delaunay triangulation-编程思维

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

整个算法的实现包括三个部分

  1. initialization
  2. sweeping
  3. 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轴的值在输入点最小和最大值的范围之外。公式如下:

\[P_{-1} = (x_{min} - \delta_x, y_{min} - \delta_y) \\ P_{-2} = (x_{max} + \delta_x, y_{min} - \delta_y) \\ \delta_x = \alpha * (x_{max} - x_{min}) \\ \delta_y = \alpha * (y_{max} - y_{min}) \\ \alpha >0 \]

初始化的位置类似如下图所示:

对于虚拟点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\ne\{P_L,P_R\}\),如图7a所示。然后构建出三角形\(\triangle 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)。如果角度大于\(3\pi/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).

版权声明:本文版权归作者所有,遵循 CC 4.0 BY-SA 许可协议, 转载请注明原文链接
https://www.cnblogs.com/grass-and-moon/p/14872636.html