Delaunay三角剖分程序文档

算法简介

本次实验我们实现了简单凸多边形的Delaunay三角剖分的分治算法，算法的复杂度为O(NlogN)。算法输入为一系列点的集合，而输出为由这些点所构成的简单凸多边形对应的Delaunay三角剖分。

传统的简单凸多边形的Delaunay三角剖分算法有随机增量方法、扫描线方法和分治方法。这三个算法的理论复杂度均为O(NlogN)，然而，在实际应用做，效率最好的还是分治算法。

算法原理说明

Delaunay三角剖分分治算法说明

整体流程

算法的整体流程如下：

点集初始处理包括两个步骤。首先，对所有点进行按x坐标进行排序（如果x坐标相同，则按y坐标排序），然后，去掉重复点，即x坐标和y坐标均相同的点。

分治算法中，采取的办法是不断地对点集从中间进行剖分，知道剖分到每个点集只含有2-3个点时，则直接加入边，2个点时构成线段，而3个点则构成一个三角形。之后不断地合并这些剖分的点集，合并操作需要先求解上下公切线操作，之后再将两个Delaunay三角剖分的点集合并为一个Delaunay三角剖分点集。通过这样一次次的合并，最后可以得到最初的整个点集的Delaunay三角剖分。

最后，我们将已经得到的结果返回，在屏幕上显示。

数据结构

在程序中我们用来表示点和边为修改后的DCEL（Doubly-Connected Edge List）结构。在点中记录点的坐标信息，以及以他为起点的一条边的指针，即：

class Vertex

{

public:

//坐标信息

double x;

double y;

//该边指向从当前顶点出发的任意一条半边

HalfEdge* pEdge;

};

在边中记录该边对应的起始点在点数组中的下标，对边指针，按逆时针方向以同一点为起始点下一条边的指针和前一条边的指针，即：

class HalfEdge

{

public:

//半边起点在vector里的编号

//Vertex* ori;

int oriNum;

//孪生边

HalfEdge* twin;

//沿着顺时针方向遇到的第一条同起点半边

HalfEdge* pre;

//沿着逆时针方向遇到的第一条同起点半边

HalfEdge* next;

}

具体结构如下图所示：

下切线和上切线算法

算法流程说明：

给定两个已有的三角剖分，寻找这两个点集的下公切线：

1)      X为左侧点集的最右点，Y为右侧点集的最左点

2)      Z为Y沿右侧点集凸包逆时针方向的连接第一个点，Z’为X沿左侧点集凸包顺时针方向的连接第一个点

3)      判断Z’是否在YX连线的右侧且Z是否在XY连线的左侧。如果这一条件同时满足，则XY为下切线，返回。否则若Z’不满足条件，则令Z’等于Z’ 沿左侧点集凸包顺时针方向的连接第一个点，若Z不满足条件，则令Z等于Z沿右侧点集凸包逆时针方向的连接第一个点，继续执行第4步

4)      判断Z’是否在YX连线的左侧，如果在，则令Z’等于Z’ 沿左侧点集凸包顺时针方向的连接第一个点，重复这一过程，直到Z’不在YX连线的左侧为止

5)      判断Z是否在XY连线的右侧，如果在，则令Z等于Z沿右侧点集凸包逆时针方向的连接第一个点，重复这一过程，直到Z不在在XY连线的右侧为止

6)      连接Z和Z’，将X设为新的Z’点，将Y点设为新的Z点，并且令Z’为X沿左侧点集凸包顺时针方向的连接第一个点，Z为Y沿右侧点集凸包逆时针方向的连接第一个点，返回第3步继续执行

算法的一些过程如下图所示：

最初的最左点和最右点分别为Y和X，通过上述算法过程，在第一轮循环结束后，找到了新的X(New X)和Y(New Y)，此时通过判断，此时的XY线即为下公切线。

同理，上公切线线也可以采用同样算法求解，只不过每回判断的是Z’

合并算法

首先，在进行合并算法之前，我们需要知道两个带合并点集的上公切线和下公切线。同时，定义以下几个操作：

PRED(A, B)返回以A为起始点，从AB以A为轴顺时针旋转遇到的第一条边的终点；

SUCC(A, B)返回以A为起始点，从AB以A为轴沿逆时针旋转遇到的第一条边的终点；

INSERT(A, B)将AB边插入到对应的DCEL结构中；

DELETE(A, B)将AB边从对应的DCEL结构中删除；

QTEST(H, I, J, K)判断K点是否在H、I、J三点的外接圆中，如果不在返回TRUE，否则返回FALSE。

合并两个点集的Delaunay三角剖分算法流程如下（这里“=”表示赋值操作）：

1)      初始化操作：令BT为对应的下公切线，UT为上公切线，L为BT的左端点，R为BT的右端点，A和B为两个标志值

2)      判断BT是否等于UT，若等于，则合并过程结束，返回结果，否则继续进行后面的操作

3)      令A=B=FALSE，INSERT(L, R)，并令R1=PRED(R, L)，L1=SUCC(L, R)

4)      判断R1是否在LR连线的左侧，如果在，则执行下面这一循环：令R2=PRED(R, R1)，判断QTEST(R1, L, R, R2)是否为TRUE，若为TRUE，则退出循环；否则，DELETE(R, R1)，令R1=R2，R2=PRED(R, R1)，然后重复执行这一循环。如果R1不在LR连线的左侧，则令A=TRUE

5)      判断L1是否在RL连线的右侧，如果在，则执行下面这一循环：令L2=SUCC(L, L1)，判断QTEST(L, R, L1, L2)是否为TRUE，若为TRUE，则退出循环；否则，DELETE(L, L1)，令L1=L2，L2=SUCC(L, L1)，然后重复执行这一循环。如果L1不在RL连线的右侧，则令B=TRUE

6)      如果A为TRUE，令L=L1，否则，如果B为TRUE，令R=R1。如果A和B均不为TRUE，则判断QTEST(L, R, R1, L1)的值，为TRUE，则令R=R1，否则，令L=L1

7)      令BT为LR的连线，返回第2步

算法中间步骤如下图所示：

判断右侧点集是否符合算法要求。

判断左侧点集是否符合算法要求。

判断应该挪动哪一个端点，以作为下一次循环的LR

连接LR，作为新的BT，返回算法第2步，重复上述过程，知道BT与UT相等，合并操作完成。

判定点是否在指定三点外接圆内算法

在上述算法中，需要判断一个点是否在其他指定三点对应的外接圆内，我们采用如下算法进行判断：

给定三角形的三个点A(x1,y1)，B(x2,y2)，C(x3,y3)，将其投影到三维抛物面z = x2+y2上后可确定一平面F，如果给定点P(x0,y0)在抛物面上的投影在F上方，则P在A,B,C所确定的三角形外接圆外，P在F下方则表示P在三角形ABC外接圆内，P在F上则表示P刚好位于三角形ABC的外接圆上。

根据此原理，将A,B,C，P分别转换成空间点D(x1,y1,x12+y12)，E(x2,y2,x22+y22)，F(x3,y3,x32+y32)和点G(x0,y0,x02+y02)，然后求向量DE和向量DF的叉积向量v，v向量垂直于DEF三点所确定的平面。如果v向量与z轴反向，则将v调整为与z轴相对同向，即向量v与空间向量(0,0,1)点积应该为正。接下来做向量v与向量DG的点积，其结果大于0，小于0，等于0分别对应点G位于平面DEF上方，下方和刚好位于平面DEF上。由此可以确定点P是否在ABC所确定的外接圆内。

复杂度分析

查找下切线和上切线的复杂度显而易见为O(n)。

对于合并操作，由算法可以看出，因为每条被删除的边之后都不会再被添加，所以删除边的操作数为O(n)，同理，插入边的次数同样也是O(n)级的，另外，判断是否在圆内以及判断处于左侧右侧这些操作的时间都是常数量级的。所以，加上寻找上下公切线的算法，一次合并操作的复杂度也是O(n)。令t(N)表示构造n个点的Delaunay三角剖分所需要的时间，则有如下公式：

从而可以得到算法的复杂度为O(NlogN)

穿刺算法说明

穿刺算法流程

穿刺算法是用来判定剖分结果是否正确的手段。根据用户输入，可以得到穿刺的起点和终点。由于三角剖分的边界为输入点集的凸包，首先沿着凸包的边界检查是否有边与穿刺相交，根据相交边的个数确定穿刺与剖分的位置关系，具体分为三种情况：起点和终点都在剖分内，起点和终点都在剖分外，起点和终点有且仅有一个在凸包内。

接下来沿着起点到终点的方向找到剖分中第一条与穿刺相交的边startEdge和最后一条与穿刺相交的边endEdge，并将这些边加入到结果集合中，接下来沿着穿刺方向从startEdge开始找下一条与穿刺相交的边，具体的找法是通过startEdge围成的三角形来找，找到三角形中下一条与穿刺相交的边edge，取edge的孪生边作为新的startEdge，如次进行下去直到找到的边为endEdge，每找到一条边就将其加入到结果集合。

算法结束后，通过顺序遍历结果集合中的半边，找到半边围城的三角形，我们就可以顺次得到从穿刺起点到穿刺终点依次经过的三角形。

如上图所示，红色的边即为穿刺找到的边。

三角形测试算法

流程图如下：

界面绘制说明

在这次的实验中，我们在显示结果时并没有采用MFC原有的绘制方法，而是采用了OpenGL来绘制显示区域，从而使得我们的绘制速度有了很大的提升，效果也更好，同时也加入了放大、缩小、平移等操作，使得显示过程稳定效率高，并且与用户的交互性好。

动画演示方法说明

在这次的实验中，我们分别加入了三角剖分和穿刺的动画演示过程，在处理时，我们分别开辟了两个线程，其中一个线程用于后台处理三角剖分，并在加入或删除边的时候加入时间上的延迟，从而慢动作的模拟算法的处理过程；另一个线程则用于不断的重新绘制，使得屏幕能够动态显示算法过程。

通过增加多线程，使得我们的程序在动画演示过程中同样可以进行拖动、放大、缩小等操作，也可以通过点击按钮结束这一动画演示过程。

另外，考虑到可能在绘制过程中会出现边刚刚生成，关系还没有建立就在屏幕上绘制的情况，或者是边刚刚删掉，但是其他边的关联还没有更新的情况，我们在这些基本操作上加了锁，以保证程序的稳定性。

测试结果

测试平台

CPU：Intel Core2 E7200 2.53GHz

内存: 2G DDR2

结果展示

三角剖分：

判断外接圆和三角形：

穿刺：

测试结果

下表是分别随机生成不同点数时，程序完成Delaunay三角剖分所需时间，其中每个数量级别上的点我们分别进行了三次计算，并取对应的平均时间。

点数

时间(ms)

时间(ms)

时间(ms)

平均时间(ms)

1000

16

15

15

15.33333

20000

64

78

78

73.33333

40000

172

188

172

177.3333

60000

266

282

281

276.3333

80000

375

391

375

380.3333

100000

484

484

485

484.3333

200000

1016

1015

1016

1015.667

300000

1578

1578

1578

1578

400000

2140

2125

2141

2135.333

500000

2750

2735

2734

2739.667

600000

3328

3328

3328

3328

700000

3907

3907

3906

3906.667

800000

4500

4500

4500

4500

900000

5125

5125

5125

5125

1000000

5750

5735

5750

5745

1500000

8844

8829

8828

8833.667

2000000

12078

12016

10278

11457.33

2500000

14453

15406

14860

14906.33

结果分析

由上图，红色的线表示按公式0.0003*N*logN得到的曲线，而蓝色的线表示实际中测得的平均运行时间。因此，从图中红色线与蓝色线的接近程度，我们可以知道程序运行时间的变化规律符合O(NlogN)这一趋势，与算法的复杂度为O(NlogN)的理论结果相一致。

总结

问题与解决

在这次的实验中，我们同样也遇到了不少的问题。首先一个就是判断点在外接圆内的问题，最开始，我们采取的是通过求解圆心和距离判断点是否在外接圆内，可是在之后的测试中，我们发现当点的坐标挨得很近时，这时求解出的外接圆的半径相当于无限大，已经超出了double的表示范围，此时的判断已经不可靠了。所以，我们采用了文档前面提到的判断点在外接圆内外的算法解决了这一问题。

另外，在判断ToLeft和ToRight时，我们也需要考虑当点在线上的处理情况，所以，很多时候一些错误也是因为忽视了这些特殊情况，所以在程序中不同地方，我们也都针对这些情况进行了处理，比如有些时候该用ToLeft，而有些时候则改用!ToRight。

最后，在处理动画演示的过程中，我们采取的方法是多线程处理，因为之前没有考虑到给插入边和删除边这些基础操作加锁处理，使得一开始的动画演示总是时不时的出问题，很大的可能是边刚刚生成，关系还没建立好，就在屏幕上显示，导致出现空指针之类的错误。所以我们后来在这些基础的操作上进行了加锁处理，解决了这一问题。

收获

在这次的实验中，我们收获很大。首先是更深刻的理解了Delaunay三角剖分的一些性质，同时对分治算法的一些细节有了自己的一些体会。同时学会了使用DCEL结构来处理一些图的问题。

同时，在编程技巧上，掌握了不少多线程的编程知识，同时也学会了用OpenGL进行绘制的一些基本操作。

参考文献

[1] L Guibas, J Stolfi. Primitives for the manipulation of general subdivisions and the computation of Voronoi. ACM Transactions on Graphics (TOG), 1985

[2] D.T.Lee and B.J.Schachter. Two Algorithms for Constructing a Delaunay Triangulation. International Journal of Computer and Information Sciences, Vol. 9, No. 3. 1980.