带断层约束的Delaunay三角剖分混合算法＊

张群会，解子毅

（西安科技大学计算机科学与技术学院，陕西 西安 710054）

0 引言

目前比较流行的构建数据高程模型(DEM)的方法主要有2种:一种是基于不规则三角形(Triangular Irregular Networks，TIN)网格法，另一种是四边形(GRID)网格法。和GRID相比，TIN能够精确地表达网格边界和断层，是一种比较理想的三维层面建模方法，在地质层面可视化方面得到了广泛应用。Delaunay三角剖分法是建立TIN模型最常用的方法。由于Delaunay三角化满足最小角最大准则和外接圆不包含其他点准则，总是能尽可能避免狭长的三角形，自动向等边三角形逼近，具有网格形态优美等特点。因此利用Delaunay三角剖分来建立带断层约束的地质模型具有十分重要的价值。

现阶段的 CDT生成算法较多，主要以2步法［1－2］为主，其基本思想是第一步构建 DT网格；第二步将约束线段强行嵌入到DT网格中进行整体剖分。传统CDT两步法在约束线段嵌入方面的方法很多，比较流行的有3种:一种［1］是利用生长法和分治法［3］的思想，首先删除与约束线段相交的线段，形成影响域。然后以约束线段为基边把影响域分成左右2部分分别进行剖分。由于约束线段的长度没有进行细分，构网时在断层边界处可能出现狭长三角形。解决上述问题可以采用卢朝阳的对半划分增量型附加点插入算法［4］，通过对约束线段上插入数据点来优化三角网格。第一种算法简单但效率不高，网格形态也不优美。第二种是李立新优化的多对角线交换算法［5］。该算法的思想是对影响区域的三角形组成的凸四边形进行对角线交换，通过多次交换对角线可实现在任意影响域中强行嵌入约束线段。但该算法执行效率较低、编程难度大，并且也会出现狭长三角形的情况。第三种是基于复连通区域的三角剖分算法［6］。该算法首先要搜索到与约束线段相交或靠近约束线段的三角形，并把这些三角形合并成一个单连通区域，即找到包含约束线段的最小外环。然后求取该区域的外边界，将最小外环和其外边界组成复连通区域，最后将复连通区域三角剖分转化成了简单多边形的三角剖分。该方法虽然能够有效的解决约束线段的嵌入问题，但是由于其步骤较多，构网效率不高。而文中的三角剖分混合算法可以较好的解决上述问题。

1 DT的基本性质

1.1 空外接圆性质

在由点集P所形成的DT网格中，其每个三角形的外接圆均不包含P中的其他任意点。在二维空间中，给定点集P，若P中任意4个点不共圆，则Delaunay三角化是唯一的，否则Delaunay三角化不唯一。对于Delaunay三角化，局部优化可以保证全局最优，如图1(a)。

1.2 最小角度最大性质

在由点集P所能形成的三角网中，DT网格中三角形的最小角度是最大的。即在三角形网格中，对任意相邻的、构成严格凸四边形的2个三角形来说，其6个内角中的最小值大于交换四边形的对角线所形成的另外2个三角形的6个内角中的最小值，如图1(b)。

图1 Delaunay三角化的特点Fig.1 Characteristics of the Delaunay triangulation

2 数据结构

在插入点及约束边的嵌入过程中，需要由三角网拓扑信息快速搜索到插入点和约束线段所影响的三角形，并不断的更新三角网的拓扑信息。因此，保证三角网拓扑信息的正确性和有效性对于建立CDT是至关重要的。文中将三角网内查询的拓扑信息主要定义为点、边、三角形3个方面，其各自的结构如下:

3 带断层约束的三角剖分混合算法

传统CDT剖分算法所耗费的构网时间主要用在断层数据插入和约束线段嵌入这2方面上。该算法由于要不断的判断插入点是否包含在三角形的外接圆内，因此需要耗费大量的时间来计算三角形的外接圆。而插入点所影响的三角形一般都是它周围邻近的三角形，没有必要对所有的三角形进行处理。针对上述问题，文中提出的混合算法在插点方面对传统算法进行了优化:首先以插入点作为中心，然后在以搜索半径为圆的范围内寻找受插入点影响的三角形(只要搜索半径设置合理，就一定能够找到所有受插入点影响的三角形)，并只对搜索到的三角形进行处理，从而减少了处理过程，大大提高了构网效率。而在约束线段嵌入方面，文中算法首先利用均分法对约束线段进行细分，解决了约束线段过长而产生狭长三角形的问题，并且还能够使约束线段尽量包含在DT网格中，减少后续过程对约束线段的嵌入处理。

文中分3步来实现带断层约束的三角剖分混合算法。第一步首先建立一个包含所有数据点的矩形框，用逐点插入法［7－8］将数据点逐一插入到矩形框中进行剖分，然后利用重心布点法［9］对剖分好的网格进行加密，最后通过局部优化法(LOP)［10］对剖分好的三角形进行优化，形成无约束的DT网格；第二步首先对断层数据进行细分，并以细分好的数据作为中心，在一定的范围内对DT网格进行局部剖分，直到所有断层点嵌入到DT网格中。然后连接相邻的断层点，形成断层线段，判断其是否包含在网格中，如果不包含，采用三角网生长法和分治法将断层线段嵌入到网格中，直到所有断层线段处理完毕，最后用局部优化法优化网格；第三步将由断层线段围成的断层区域中的三角形删除，形成带断层约束的CDT网格。

3.1 构造无约束的DT网格

1)设边界点数组 P=(p1，p2，…，pn)，n 为边界数据点个数，三角形数组 T=(t1，t2，tm)，m 为三角形个数。构造初始矩形框，用逐点插入法在矩形框中插入数据点pj，找到T中外接圆含pj的三角形ti，并将ti从T删除，将ti的相邻边删除形成多边形空腔，把pj与空腔多边形各顶点相连，生成新三角形存入到T中。

2)用重心布点法对网格进行加密，并用局部优化法(LOP)进行优化。重心布点法的思想是首先设置三角形的最大高度为H，并通过在三角形重心位置布点，使所生成的三角网相对均匀，且三角形高度均不超过H.

3.2 构造带断层约束的CDT网格

1)利用细分算法把由断层数据构成的断层边进行细分，将细分结果有序的存储在数组F中。具体步骤如下

步骤1 从P中取出相邻的2点pi和pi+1，i=(1，2，…，n －1)，计算出 pi，pi+1的间距 di，并设置一定的步长Step，如果Step＜d转到步骤2，否则继续步骤1；

步骤2 设D=d/Step，将线段pipi+1均分为D+1段，将分段点按顺序存放在F中；

2)将细分后的断层数据点作为新点逐一插入到DT网格中进局部剖分，使断层点嵌入到DT网格中。具体步骤如下

步骤1 根据三角形的最大高度H来设置搜索半径r；

步骤2 从F中任取一点fi∈F，计算fi到ti∈T 3个顶点的距离，取其中最短距离作为di，如果di＜r，将ti存入到临时数组V中。

步骤3 判断fi是否在V中三角形的外接圆内，如果在，将三角形进行标记；

步骤4 删除标记的三角形邻边，形成多边形空腔，将fi与多边形各顶点进行V连接形成新的三角形存入到T中，清空V.迭代2，3，4步骤，直到所有断层数据点处理完毕。

图2为剖分过程图(以f1∈F为例):从图2(a)可知f1C为搜索半径r，C为圆上的任意一点，f1到 t6，t7，t8，t9的距离为 d1=f1A，f1到 t1，t2，t3，t4，t5的距离为 d2=f1B，且 d1＜ r，d2＜ r因此 t1，t2，t3，t4，t5，t6，t7，t8∈V 为 f1邻近三角形；从图 2(b)可知，f1在 t4，t5，t6的外接圆内，将 t4，t5，t6从 T 中删除；从图2(c)可知，删除 t4，t5，t6的相邻边所形成多边形空腔为Q={ABCDE}，将Q的顶点与f1进行连接，形成新的三角形 t10，t11，t12，t13，t14存入到 T中，完成该点局部剖分。

图2 空外接圆法剖分过程图Fig.2 Empty circum circle subdivision process diagram

3)将F中的第一个元素存入F，实现首尾相连，并采用生长法和分治法将断层线段嵌入到DT网格中。具体步骤如下:

步骤1 取F中相邻的两点形成约束线段fifi+1，i≤n －1，判断 fi，fi+1是否包含在 DT 网格中，如果包含继续步骤1，否则转到步骤2；

步骤2 找到与fifi+1的相交边，将相交边删除形成多边形空腔，并以fifi+1为基边将多边形空腔分成左右2个区域；

步骤3 利用三角网生长法对QL和QR2个区域分别进行剖分(以QL为例)。

1)建立一堆栈S，并把fifi+1边压入S；

2)在S中弹出一边为fifi+1作为基边，在QL顶点集合中找出距fifi+1最近的点A进行连接得到fi+1A，将三角形fifi+1A存入T中；

3)如果fiA或fi+1A是QL的一条边界边，则不入栈，反之把fiA或fi+1A压入堆栈S中；

4)如果S不空，则返2，否则退出；

5)对QL和QR剖分后形成的三角形进行局部优化(LOP)处理；

步骤4 重复上述步骤直到所有约束线段处理完毕。删除F中的最后一个元素；

从图3(a)中可以看出 fifi+1与 AD，BD，BC相交，删除线段 AD，BD，BC形成多边形空腔 Q={fiABfi+1CD}.图3(b)中约束线段fifi+1将多边形空腔分成QL={fiABfi+1}和QR={fifi+1CD}2个区域。图3(c)中QL和QR分别进行了剖分，最终实现了约束线段的嵌入。

4)将断层区域中的三角形删除，最终形成带断层约束的CDT网格。具体步骤如下

图3 影响域的形成和剖分过程Fig.3 Formation of influence domain and subdivision process

步骤1 首先利用F中的断层数据点建立断层多边形 Q=(f1，f2，…，fn)；

步骤2 计算T中三角形ti的重心gi，判断gi是否在Q内，如果在就从T中删除ti，重复步骤2，直到所有三角形处理完毕。

4 实验结果与分析

文中算法在Microsoft Visual C++2010平台上编写测试通过，实现了带断层约束的Delaunay三角剖分。对Release版本进行了时间性能测试，测试环境为:操作系统Windows 7旗舰版(64位)；处理器:AMD A8－4500 M with Radeon(tm)HD Graphics 44 CPU 1.9 GHz；内存:4 GB.图4为文中算法构造出的(a)不带断层的DT网格和(b)带断层约束的CDT网格，其时间性能见表1.

表1 算法执行效率测试结果Tab.1 Test results of algorithm performance

图4 混合算法剖分效果图Fig.4 Mixed algorithm subdivision renderings

从图4可以看出，混合CDT剖分算法所构造的DT网格和CDT网格形态都比较均匀优美，没有出现狭长三角形的情况，并且断层边能够完好的嵌入到网格中。图4(b)中的黑点为断层数据点。通过分析表1可得，在剖分相同数量的三角形情况下，文中算法在构建DT和CDT网格的时间效率上都明显好于传统算法。

5 结语

文中通过对几种传统算法进行分析和研究，逐一找出了各算法在构网过程中耗时最多的算法步骤，并在深入了解其核心思想的基础上，通过实验对比，判断出几种传统算法的优缺点，并根据研究结果将各算法的优点进行结合和优化，提出了一种带断层约束的Delaunay三角剖分混合算法。该算法通过对断层数据的插入和约束线段的嵌入过程进行优化，快速实现了带断层约束的CDT网格。最后通过实验结果的分析和对比，说明文中算法在构网效率和构网质量上都优于传统算法，在数据建模方面具有较好的应用性和实用价值。

References

［1］刘学军，龚健雅.约束数据域的Delaunay三角剖分与修改算法［J］.测绘学报，2001，30(1):82 －88.LIU Xue-jun，GONG Jian-ya.Constrained Data Delaunay triangulation and modification algorithms［J］.Acta Geodaetica et Cartographica Sinica，2001，30(1):82 －88.

［2］刘少华，程朋根，史文中.约束Delaunay三角网生成算法研究［J］.测绘通报，2004，26(1):4 －7.LIU Shao-hua，CHENG Peng-gen，SHI Wen-zhong.Constrained delaunay triangulation algorithm research［J］.Bulletin of Surveying and Mapping，2004，26(1):4 －7.

［3］谢增广.平面点集带Delaunay三角剖分的分治算法［J］.计算机工程与设计，2012，33(7):70 －73.XIE Zeng-guang.Partition of planar point set with delaunay triangulation subdivision algorithm［J］.Computer Engineering and Design.2012，33(7):70 －73.

［4］卢朝阳，吴成柯，周幸妮.满足全局Delaunay特性的带特征约束的散乱数据最优三角剖分［J］.计算机学报，1997，20(2):118 －124.LU Chao-yang，WU Cheng-ke，ZHOU Xing-ni.Meet the global feature of constraint delaunay properties of approximate optimal triangle subdivision［J］.Chinese Journal of Computers，1997，20(2):118 －124.

［5］李立新，谭建荣.约束Delaunay三角剖分中强行嵌入约束边的多对角线交换算法［J］.计算机学报，1999，22(10):1114－1118.LI Li-xin，TAN Jian-rong.Constraint delaunay triangulation subdivision forcibly embedded constraints in multiple diagonal exchange algorithm［J］.Chinese Journal of Computers，1999，22(10):1114 －1118.

［6］刘少兵.带断层地震数据的Delaunay三角剖分算法［J］.煤田地质与勘探，2008，36(6):71 －73.LIU Shao-bing.With seismic data of fault delaunay triangular subdivision algorithm［J］.Coal Geology and Exploration，2008，36(6):71 －73.

［7］Watson D F.Computing the n-dimension delaunay tessellation with application to voronoi polylopes［J］.Computer Journal，1981，24(2):167 －172.

［8］周雪梅，黎应飞.基于Bowyer-Watson三角网生成算法的研究［J］.计算机工程与应用，2013，49(6):198－201.ZHOU Xue-mei，LI Ying-fei.Based bowyer-watson triangulation generation algorithm research［J］.Computer Engineering and Applications，2013，49(6):198 －201.

［9］吴 芬.平面域Delaunay三角剖分新加密算法［J］.计算机与现代化，2007(7):19－22.WU Fen.Plane domain delaunay triangular subdivision new encryption algorithm［J］.Computer and Modernization，2007(7):19 －22.

［10］Lawson C L.Software for C1surface interpolation mathematical software Ⅲ［M］.NewYork:Academic Press，1977.