一种Delaunay三角剖分的改进算法

余代俊，蒲朝旭，朱逍贤

(1. 成都理工大学 现代工程测量技术及应用研究所，四川 成都 610059; 2. 四川科技职业学院 土木与建筑工程学院，四川 成都 611745)

一、引 言

不规则三角网(TIN)是二维平面空间对任意离散点实行GIS数据表达、管理、可视化，以及地学分析、DEM分析、计算机视觉等方面的一项不可或缺的应用技术与手段[1]。Delaunay三角网即D-TIN，它能够很好地满足TIN的3点基本要求，即唯一性、最大最小角特性、空圆特性，能够对给定区域点集进行最佳的三角剖分[2]。

正因为TIN模型具有众多优点，同时其数据结构比较简单，因此它在许多领域得到了广泛应用。但随着计算机技术和地理信息技术的发展，如何快速高效地处理大量数据和对数据进行高精度模拟与仿真就成了该方面研究的主要内容[3]。

本文就如何改进和优化Delaunay三角网的构建算法来满足海量数据的处理需求，以及如何将构建算法利用较为简洁的程序进行实现，进行了深入研究和详细探讨。

二、算法设计

笔者结合已有三角网的生成算法[4-12]，取长补短，提出了基于凸包实现的Delaunay三角网剖分算法。

1. 数据存储结构

一个好的数据结构直接影响着D-TIN的生成效率，D-TIN的数据结构主要包括点、线、面之间的组织关系和拓扑关系。

本文所使用的数据结构如下

struct Point3d∥点结构

{

double X;∥X坐标

double Y;∥Y坐标

double Z;∥Z坐标

}

struct Edge∥边结构

{

Point3d StartPt;∥起点

Point3d EndPt;∥终点

Triangle LeftTri;∥左三角形

Triangle RightTri;∥右三角形

}

struct Triangle∥三角形结构

{

Edge edge1;∥边1

Edge edge2;∥边2

Edge edge3;∥边3

}

基于上述结构的设计，三角形的各个顶点和边均采用逆时针的方式进行存储，在程序实现过程中采用链表和字典(字典是采用键-值的存储结构，能够通过键直接取出值)等数据存储结构，能够快速实现数据的存取及快速定位。

2. 算法实现

该算法主要分3步完成：首先，对输入数据进行预处理，生成凸包；其次，在凸包的基础之上采用逐点插入法生成三角网；最后对三角网采用LOP优化算法生成最终的三角网。

(1) 生成凸包

凸包的生成是基于Akl-Toussaint启发式函数[13]和斜率判别法来建立的，如图1所示。

1) 构建初始凸包：将输入点集中maxx、maxy、minx和miny的4个点或是与该点集最小外围边框的西南角和东北角最近的4个点作为初始凸包，如图1中的P11、P14、P6和P2点组成的虚多边形，并且按照逆时针的方式将此4点存入凸包链表中，如果存在多个minx最小的点，则取其中y坐标最小的点，其余类似情况依此处理。

2) 清理点：建立初始凸包之后，该初始凸包的4个点肯定位于最终凸包点集中，而该4点所组成的四边形中所包含的点则肯定不会出现在最终凸包点集中，为此使用Akl-Toussaint启发式函数将位于初始凸包中的所有点全部删除，将剩余的点用于凸包的构建计算。该Akl-Toussaint启发式函数的性能为O(n)，对于大量的随机点而言，该法能够排除掉几乎一半的点，能够明显提高后续建立凸包的性能。

图1中的初始凸包由P11、P14、P6和P2组成，将该4点组成的四边形内所包含的点全部剔除，即剔除6个点(P8、P9、P5、P4、P10、P13)，只留下剩余的4个点(P1、P3、P7、P12)来判断是否为凸包点，明显提高了程序的运行效率。

3) 排序点：对进行清理之后的点集按照X坐标升序为主、Y坐标升序为辅的方式进行排序，以便用于后续凸包的生成。对于点集的排序可以采用LINQ方式进行，该方法是所有排序算法中最高效的。

4) 修改凸包：取出初始凸包链表中相邻的两个点，检查位于其右侧斜率最小的点，该点亦即是其右侧距离该两点连线距离最大的点(可以采用三角形面积法判断)。若该点存在，则将该点插入到上述两点之间；否则不插入，进行下一条边的判断。循环执行上述过程，直至所有边均没有插入点则完成凸包的搜索。在进行点搜索的时候，每添加一个点，也可以添加Akl-Toussaint启发式函数的判断，将位于新加入点与该边组成的三角形中的点剔除掉。

5) 生成凸包：将修改之后的凸包代替初始凸包则完成了点集凸包的搜索。

在构建凸包过程中，其时间复杂度主要是修改凸包所花费的时间，其时间复杂度和Akl-Toussaint启发式函数的性能一致，均为O(n)。

(2) 建立三角网

1) 预处理：将所有凸包点从原始点集中剔除，所有凸包点已经是点集的最外围点，不再用于插入构网。

2) 构建初始三角网：从剔除凸包点之后的点集中任意取出一点P9，如图2所示。将该插入点与凸包的所有顶点连接，将生成的三角形存入三角形链表中，将各边存入边链表中。

图2 Delaunay三角网建立过程示意图

3) 修改三角网：将剩余的点依次插入到该初始三角网中，每插入一个点需要判断该点位于哪一个三角形内，同时对该三角形进行重构，而不用重构整个三角网，能够大大提高构网效率。判断点位于哪一个三角形时需要遍历三角形链表，在遍历比较时可以进行排斥试验，如果该点在该三角形的最小外围边框(BoundingBox，矩形边与WCS的X、Y、Z轴平行)之外，则肯定不在该三角形内；如果在该矩形框内，则再使用内角和法、同向法及重心法等方法进行判断，也可以使用射线法进行判断等，但进行排斥试验之后能够节约很多程序运行时间。

在建立初始三角网之后，任意取一点，此处以点P13作为示例(如图2所示)，插入到初始三角网中，经过计算可知点P13位于△P9P12P14中，将点P13与该三角形的3个顶点分别连接，形成新的3个三角形，并加入到三角形存储链表中，新生成的3条边存入边链表中，并且更新原始边链表中△P9P12P14的3条边的属性，如图2中的虚线所示，同时在三角形链表中删除原始的△P9P12P14[14]。

4) 生成三角网：所有点插入完成之后，将最终的三角形链表和边链表代替初始三角网链表，则得到最终的三角网。

该部分程序的算法在当数据点排列比较规则、分散比较均匀的情况下，其时间复杂度为O(n)，在最坏的情况下，其时间复杂度为O(n2)，其时间复杂度平均保持在O(nlgn)。

(3) LOP优化

LOP优化即局部优化方法，是Lawson于1977年提出的，它的核心思想是通过交换凸四边形的对角线来获得等角性更好的TIN。

1) 预处理：将凸包边从边链表中剔除，凸包边是最外围的边界，只包含在一个三角形中，故不再用于优化处理。

2) 改进LOP优化方法：LOP优化的目的是满足空圆特性和最大最小角特性，即保证最邻近的点构成三角形且三角形尽量接近等边。文献[15]中提出了通过比较凸四边形中对角线的长短来进行三角形优化处理的方法，该方法选取最短的一条对角线作为该凸四边形的划分线，但是该法在某些情况下无法处理，如图3所示。

图3 改进LOP优化算法示意图

图3中，LP2P4=51.71，LP1P3=54.02，如果按照文献[15]的说法，取最短的边作为该四边形的分割线，则应该取图3中的虚线(对角线)，而实际上应该取图3中边P1P3作为该四边形的分界线。从图中可以看出，∠P1P4P3为55.96°，∠P1P2P3为102.93°，∠P2P3P4为147.44°，∠P2P1P4为53.67°，如果将∠P1P4P3与∠P1P2P3相加，将∠P2P3P4与∠P2P1P4相加，并将它们各自的和分别与120°相减可以得出：∠P1P4P3+∠P1P2P3-120°<∠P2P3P4+∠P2P1P4-120°。

之所以要与120°相减是因为D-TIN需要满足最大最小角特性，即三角形尽量接近等边，等边三角形中各角均为60°，同时也很好地满足空圆特性。

从上述讨论可以看出，与120°相差最小的那一对角所对的对角线恰好为所选择的对角线，如图3中的边P1P3，该方法同样适用于其他凸多边形。

基于上述计算和检验笔者提出了角度判别对角线的方法，该方法能够很好地处理LOP优化，同时该方法的时间复杂度仅为O(n)，能够很好地完成三角网的局部优化。

3) LOP优化：基于上面所提出的角度判别对角线的方法，使用改进后的LOP优化方法进行三角网优化的处理过程如下：

从边链表中任意取出一条边，判断该边的左三角形和右三角形组成的四边形是否为凸四边形，如果是凹四边形则不进行后续操作而继续取下一条边；若是凸四边形则需要按照介绍的方法进行LOP优化处理，如果交换了对角线，则需要将交换前的两个三角形从三角形链表中删除，同时加入新生成的三角形，还需要将交换前的边从边链表中删除，同时加入新生成的边到边链表中。

由于设计边链表时存储了边的左右三角形，因此进行LOP优化时可以直接取出该两个三角形进行优化，省去了对左右三角形的搜索时间，同时三角形与三角形之间通过边进行链接，各个三角形之间形成有向链接，能够很方便地进行等值线生成等后续处理。

依次循环检查各条边，当边链表循环完毕则LOP优化处理结束，此时输出的三角形链表和边链表为最终Delaunay三角网构网结果。

三、算法运行效果

将本文介绍的方法采用C#语言实现，对某一测区10 235个数据点进行构网，先利用Akl-Toussaint启发式函数建立凸包，然后再将凸包以外的数据点进行插入，最后利用角度判别对角线法进行LOP优化以生成Delaunay三角网。图4是原始数据点集，图5是生成的Delaunay三角网的效果图。

图4 原始离散点数据集

图5 三角网生成后的效果

笔者利用笔记本电脑，处理器为i5，采用本文所提出的方法生成上述三角网的时间(包括数据的读取时间)为2.67 s，与分治算法的2.70 s基本齐平。在综合对比了其他数量级的数据生成算法时间之后，得出该算法比较稳定，且唯一性较好，算法的效率大部分情况下与同样唯一性较好的分治算法相当，偶尔比分治算法高效。

在常用的几种三角剖分算法中[16-17]，逐点插入法的时间复杂度平均为O(n2)，三角网生长算法的时间复杂度平均也为O(n2)，分治算法的时间复杂度虽然也可以降到O(n)，但是其子块的递归算法的时间复杂度仍然达到了O(n2)。

本算法从凸包的生成、点的插入、三角网的优化等方面入手，采用了Akl-Toussaint启发式函数、边的有向性存储和角度判别对角线法进行LOP优化等，大大降低了程序的时间复杂度，使时间复杂度平均保持在O(nlgn)，在最好的情况下能够达到O(n)，在最坏情况下也小于O(n2)，明显小于常用三角剖分的时间复杂度。

四、结束语

本文将凸包法与逐点插入法相结合，提出了对于三角剖分算法的改进方法。此法从凸包生成、三角网存储结构、三角形的快速定位，以及LOP优化等方面对算法进行优化，应用Akl-Toussaint启发式函数、三角形边的有向性存储，以及角度判别对角线法等策略，大幅度提高算法的运行效率。从算法分析和程序运行效果可以看出，该算法不仅能够完全满足大数据量的Delaunay三角剖分对高效和高精度的需求，而且其算法比分治算法更加容易理解、容易实现。

参考文献：

[1] 胡鹏,杨传勇,吴艳兰，等.新数字高程模型:理论、方法、标准和应用[M].北京：测绘出版社，2007.

[2] 张宏,温永宁,刘爱利，等.地理信息系统算法基础[M].北京：科学出版社，2006.

[3] 凌海滨,吴兵.改进的自连接Delaunay三角网生成算法[J].计算机应用，1999，19(12):10-12.

[4] LAWSON C L. A Software for C Surface Interpolation[J].Mathematical Software，1977(3): 161-194.

[5] LEE D T, SCHACHTER B J. Two Algorithms for Constructing a Delaunay Triangulation[J]. International Journal of Computer and Information Science, 1980,9(3):219-242.

[6] BOWYER A. Computing Dirichlet Tesselations[J]. Computer Journal,1981(24): 162-166.

[7] MACEDONIA G, PARESCHI M T. An Algorithm for the Triangulation of Arbitrarily Distributed Points: Applications to Volume Estimate and Terrain Fitting[J]. Computers & Geosciences,1991(17): 859-874.

[8] FLORIANI L D， PUPPO E. An On-line Algorithm for Constrained Delaunay Triangulation[J].CVGIP: Graphical Models and Image Processing, 1992, 54(3): 290-300.

[9] 武晓波,王世新,肖春生.Delaunay 三角网的生成算法研究[J].测绘学报,1999，28(1):28-35.

[10] MCCULLAGH M J.Delaunay Triangulation of a Random Data Set for Lirarithmic Mapping [J].The Cartographic Journal,1980(17):93-99.

[11] MIRANTE A,WEINGARTEN N.The Radical Sweep Algorithm for Constructing Triangulated Irregular Networks[J]. IEEE Computer Graphics and Application,1982，2(3):11-21.

[12] 芮一康,王结臣.Delaunay三角形构网的分治扫描线算法[J].测绘学报,2007，36(3):358-362.

[13] HEINEMAN G T, POLLICE G，SELKOW S.算法技术手册[M].杨晨,李明，译.北京：机械工业出版社，2010.

[14] 徐道柱,刘海砚.Delaunay三角网建立的改进算法[J].测绘与空间地理信息,2007，30(1)：38-41.

[15] 魏向辉,夏春林,鲁庆.一种基于凸包的Delaunay三角网算法设计[J].测绘科学,2010，35(5):152-153.

[16] 邵春丽,胡鹏,黄承义,等.Delaunay三角网的算法详述及其应用发展前景[J].测绘科学，2004，29(6):68-71.

[17] 余杰,吕品,郑昌文.Delaunay三角网构建方法比较研究[J].中国图象图形学报,2010，15(8):1158-1167.