Delaunay算法的实现与应用.docx

1、Delaunay算法的实现与应用分类号：TP311.1 U D C：D10621-408-(2007)5855-0密 级：公 开 编 号：2003031303信息工程学院学位论文Delaunay算法的实现与应用论文作者： 申请学位专业：计算机科学与技术申请学位类别：工学学士指导教师（职称）： 论文提交日期： Delaunay算法的实现与应用摘 要数字地形模型是针对地形地貌的一种数字建模,这种建模的结果通常就是一个数字高程模型（DEM）。不规则三角网（TIN）模型是DEM中存储和表示非规则数据的理想模型，它既减少规则网格方法造成的数据冗余，同时在计算效率方面又优于纯粹基于等高线的方法，所以寻求一

2、种好的TIN算法更能快速逼真的显示与模拟出地貌三维信息。在所有可能的三角网中，狄洛尼（Delaunay）三角网在地形拟合方面表现最为出色，因此常常用于TIN的生成。依据Delaunay三角剖分准则，直接以边为基础向一侧推进,而不是以凸包为基础向推进,从而极提高了Delaunay三角网推进的速度。仿真实验表明，改进后的算法效率有了显著的提高。关键词：数字地形模型；数字高程模型；不规则三角网；Delaunay三角网Delaunay Triangulation Algorithm Realization & ApplicationAbstractDigital Elevation Model(DEM

3、) is a digital modeling process which aims at terrain and physiognomy. Irregular triangulation TIN is the best model when DEM data are stored and expressed. Besides reducing the redundancy of the data caused by regular raster model, it also presents the method purely based on contour lines in calcul

4、ate efficiency. So a well developed arithmetic can show and simulated 3-Dimension information of terrain and geomorphology more quickly and vividly. Among all the available ones, Dlaunay triangulation is the best to simulate the terrain. And so it is used to create TIN usually. According to the anal

5、yse rule, the edges were used as the base when going forward ,other than Vononoi figure as the base. Consequently, the speed of constructing Delaunay triangle was greatly improved . The result of simulating shows that the efficiency of mended algorithm is evidently enhanced.Key words: Digital Elevat

6、ion Model；Digital Terrain Model；Triangulated Irregular Network；Triangulated Delaunay Network论文总页数：19页1 引言 11.1 课题背景 11.2 国外研究现状 11.3 本课题研究的意义 11.4 本课题的研究方法 22 Delaunay方法的基本原理 22.1 Voronoi图与Delaunay三角网的基本概念 22.2 Delaunay的重要性质 32.3 传统Delaunay生成步骤 33 三角剖分改进法 43.1 算法基本流程 43.2 Graham扫描法求凸包 53.3 详细算法描述 53

7、.4 程序运行结果 74 Super三角改进算法 84.1 算法基本流程 84.2 Super三角形的生成 94.3 详细算法描述 94.4 程序运行结果 104.5 面向对象计算机的实现 114.6 测试结果与算法分析 125 Delaunay算法的应用 135.1 插值基本原理 135.2 笔者源程序 145.3 基于网格插值的等值线生成 15结 论 16参考文献 16致 18声 明 191引言1.1课题背景三角网格化问题可以追朔到1907年,G.Voronoi首先提出了此问题.后来Delaunay在1932年首次提出了解决这一问题的方法.近年来,平面任意点集的三角网格化(triangul

8、ation)问题一直是人们密切关注的问题.真三维的地理信息系统的实现仍然存在诸多尚未解决的技术难题。首先，空间三维数据的采集，其成本相当昂贵；其次，空间数据量大，种类多，结构复杂；第三，三维空间的点、线、面和体之间的拓扑关系复杂，技术尚不成熟；第四，空间分析困难。因此，在地理信息的三维可视化（特别是地形三维可视化）的研究中，通常采用2.5维的GIS可视化的方法来实现地理信息的三维可视化。而该方法主要又是以高质量的数字高程模型（DEM）和高逼真度的三维显示技术为基础，其中DEM的质量，对地形三维可视化的效果有着不容忽视的影响；而影响DTM质量的关键是生成DEM的算法。所以，采用一种实用性高、精度

9、较高、生成速度快、使用方便的DEM的生成算法十分必要。不规则三角网（TIN-Triangulated Irregular Network）是一种表示数字高程模型的方法，它既减少规则网格方法造成的数据冗余，同时在计算效率方面又优于纯粹基于等高线的方法。1.2国外研究现状自1934年,俄国数学家Delaunay提出三角形最小角最大的三角化准则,并证明在四点或四点以上共圆条件下的平面散乱点存在的三角化方式后,近年来,有很多学者都致力于三角化理论及各种应用的研究, 比如: Lawson等人提出三角化的最大角最小化原则,使得三角化局部更为均匀。Amenta等人提出的外壳(C rust)算法即基于计算几何

10、中voronoi图和Delaunay三角化的全新的曲面重建算法。Hoppe提出的累进网格 ProgressiveMesh简称PM)表示方法,率先讨论了多分辨率流形三角网格累进传输问题。Aumann用三角网格近似给定的曲面，把三角网格的面积作为原曲面的近似，进而给出了斜直纹曲面近似保面积展开的一个简单算法。Azevedo针对Delaunay三角网的特点,提出的使三角网最大外接圆最小的优化方法,Rajan则将这种方法应用到了三维空间。国方面，周培德教授提出的最小权三角划分方法(三角网格划分中各三角形边的权值之和为最小值的方法)对三角网格进行了优化等等。由此可见，目前国际上关于三角化问题及其应用的研

11、究很多,而国这方面研究还为数不多,因此开展这方面的研究很有必要。1.3本课题研究的意义三角网格化问题在许多领域有广泛应用.例如,在有限元分析中,它是有限元网格生成的有效方法之一;在计算机辅助几何设计中,它可以简化几何模型和实体模型的数据结构及部表示;在计算机图形显示中,它可以加速光线跟踪算法和光照模型的处理以及图形图像插值应用;在地理学中,它可以快速的生成地理模型。Delaunay三角剖分甚至在无线Ad Hoc通信网络中，可以用于构建发现移动节点间通信路径的在线路由算法。目前，在计算流体力学、统计学、气象学、固体物理学、计算几何学、图形图像学、三维仿真等多个邻域，都能见到它的踪影。我们可以举一

12、些简单的应用：近年来，随着网上购物越来越走入人们的生活，我们可以把某些商品建成可视又可“触”的模型，人们在选购的时候不仅可以看到商品的样子，还可以通过鼠标和键盘对商品模型进行各种操作，查看它的弹性、柔软度等性质，从而用更直观、更便捷的方法来了解商品的性能，而不是对着大篇幅的数字指标发呆。又如，我们可以建立仿真的地理模型，来做大量的地震、暴雨 等的灾难模拟实验，为灾难预测提供理论依据，从而节省大量开支。所以，我相信三角网格化的研究工作一定有着良好的发展和应用前景。1.4本课题的研究方法三角网格化主要有两种准则:一种称为 Delaunay三角剖分,即在生成的三角形网格中,各三角形的最小角和为最大;

13、另一种是在生成的三角网格中,所有三角形的边长和最小.其中, Delaunay三角剖分是目前研究应用最广的一种剖分方法.本课题的研究方法主要是以Delaunay三角网的两个重要性质（空外接圆性质和最大最小角度性质）以及Delaunay三角网的基本原理为基础，参照传统算法思路，在构建三角网的过程中，改进算法的实现方法，数据结构，以达到提高效率的目的。2Delaunay方法的基本原理2.1Voronoi图与Delaunay三角网的基本概念Voronoi图（以下简称V图），又叫泰森多边形或Dirichlet图，它是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成。N个在平面上有区别的点，按照最邻

14、近原则划分平面；每个点与它的最近邻区域相关联。Delaunay三角形是由与相邻Voronoi多边形共享一条边的相关点连接而成的三角形。Delaunay三角形的外接圆圆心是与三角形相关的Voronoi多边形的一个顶点。Voronoi三角形是Delaunay图的偶图。为清晰V图及Delaunay三角网的概念，如图1，设为欧几里德平面上的一个点集，则每一个对应区域称为点的Voronoi区域，各点的Voronoi区域共同组成V图（图中用黑色连线）。当为非共线点集时，若其对应V图中有两个点， 的Voronoi区域有公共边，就连接这两个点，以此类推遍历这n个点，可以得到一个连接点集的唯一确定的网格，就是D

15、elaunay三角网格（图中用红色连线）。2.2Delaunay的重要性质空外接圆性质：在由点集V生成的Delaunay三角网中,每个三角形的外接圆均不包含该点集的其他任意点。最大最小角度性质：在由点集V生成的Delaunay三角网中，所有三角形中的最小角度是最大的，即在生成的三角形网格中,各三角形的最小角和为最大。唯一性：不论从区域何处开始构网，最终都将得到一致的结果。由于以上特性,决定了Delaunay三角网具有极大的应用价值。Miles证明了Delaunay三角网是“好的”三角网；Lingas进一步论证了“在一般情况下，Delauany三角网是最优的。”同时以上特性也成为建立Delaun

16、ay三角网的重要算法依据。2.3传统Delaunay生成步骤Tsai(1993)提出的在6 维欧拉空间中构造Delaunay 三角形的通用算法：凸包插值算法。该算法主要包括3 步：第一步：生成包含所有离散点的多边形（1）求出点集中满足min(x-y)、min(x+y)、max(x-y)、max(x+y)的四个点，并按逆时针方向组成一个点的链表。这四个点是离散点中与包含离散点的外接矩形的四个角点最近的点。这四个点构成的多边形作为初始凸包。（2）对于每个凸包上的点I，设它的后续点为J，计算矢量线段IJ 右侧的所有点到IJ的距离，求出距离最大的点K。（3）将K 插入I、J 之间，并将K 赋给J。（4

17、）重复（2）、（3）步，直到点集中没有在线段IJ 右侧的点为止。（5）将J赋给I,J 取其后续点，重复（2）、（3）、（4）步。（6）当凸包中任意相邻两点连线的右侧不存在离散点时，结束点集凸包求取过程。第二步：环切边界法凸包三角剖分在凸包链表中每次寻找一个由相邻两条凸包边组成的三角形，在该三角形的部和边界上都不包含凸包上的任何其它点。将这个点去掉后得到新的凸包链表。重复这个过程，直到凸包链表中只剩三个离散点为止。将凸包链表中的最后三个离散点构成一个三角形，结束凸包三角剖分过程。 第三步：离散点插（1）找出外接圆包含待插入点的所有三角形，构成插入区域。（2）删除插入区域的三角形公共边，形成由影响

18、三角形顶点构成的多边形。（3）将插入点与多边形所有顶点相连，构成新的Delaunay 三角形。（4）重复（1）、（2）、（3），直到所有非凸壳离散点都插入完为止。完成这一步后，就完成了Delaunay 三角网的构建。3三角剖分改进算法3.1算法基本流程总的算法流程图:图2 算法1流程图算法中的第四部分预定义过程(加入其余离散点)实际包括三个步骤：首先要确定插入点所在的三角形，然后确定插入点的影响域,最后重构影响域的三角网。此步骤是影响效率的主要步骤，在程序运行的时候是被重复执行最多的步骤，所以算法的时间复杂性主要取决于它们的执行效率。并且，随着离散点的增加三角形数以点数约2倍的速度增加，如何降

19、低以上步骤中的查找和定位算法的复杂性及缩短算法的执行时间成为提高算法效率的关键。3.2Graham扫描法求凸包凸包的概念：点集Q的凸包(convex hull)是指一个最小凸多边形，满足Q中的点或者在多边形边上或者在其。下图中由红色线段表示的多边形就是点集Q=p0,p1,.p12的凸包。 图3 凸包形成示意图算法伪代码描述:对于一个有三个或以上点的点集Q，Graham扫描法的过程如下： 令p0为Q中Y-X坐标排序下最小的点设为对其余点按以p0为中心的极角逆时针排序所得的点集（如果有多个点有相同的极角，除了距p0最远的点外全部移除压p0进栈S压p1进栈S压p2进栈Sfor i 3 to mdo

20、while 由S的栈顶元素的下一个元素、S的栈顶元素以及pi构成的折线段不拐向左侧对S弹栈压pi进栈Sreturn S;此过程执行后，栈S由底至顶的元素就是Q的凸包顶点按逆时针排列的点序列。需要注意的是，我们对点按极角逆时针排序时，并不需要真正求出极角，只需要求出任意两点的次序就可以了。而这个步骤可以用矢量叉积性质实现。 3.3详细算法描述算法基于上述的传统构建算法，但仅有两步：第一步：（1）在离散点集中寻找一纵坐标最小的点A。（2）以点A为起点，寻找两个点B、D，使得向量AB与横坐标轴夹角最小，向量AD与横坐标轴夹角最大。若点A、B、D共线，将原B点标记为A，寻找点D，使得向量AD与直线AB

21、夹角最大；寻找点C使得向量BC与线段AB夹角最小。否则，若A、B、D不共线，则寻找点C使得向量BC与线段AB夹角最小。这样，所有点都在逆时针旋转的折线DABC的左侧。（3）上面一步生成的点C、D如果为同一点，则ABC（或ABD）即为包含所有不规则点的Delaunay三角形，生成凸包的过程结束跳过一下各步；否则，继续第四步。（4）寻找一个距离直线AB最远的点E，过E作直线AB的平行线，与射线AD、BC分别交于点D、C，将点D、C重新标记为点D、C，则凸四边形DABC包括所有不规则点（5）判断ABC的外接圆是否包含点D，若是则求ABC的外接圆半径R，在射线BC上距点C2.1R远处取一点D，并将点D

22、重新标记为点D，则凸四边形DCBA即为所求得的初始凸包。若ABC的外接圆不包含点D，则凸四边形DCBA即为所求得的初始凸包。第二步：不规则点插在原有Delaunay三角形网的基础上，将其余离散点依次插，形成Delaunay三角网。该算法原理遵循传统（TSAI）离散点插算法。插的计算机实现：将第一步中形成的初始Delaunay三角形放入Delaunay三角形链表。并建立临时三角形链表，放置新生成的三角形，初始值为空。（1）若没到点集链表的尾端，按顺序取不规则点中的一个点O；将临时三角形链表清空。（2）若Delaunay三角形链表不为空，从该链表中按顺序取一个，称为当前三角形。若为空，转（5）。（

23、3）若当前三角形的外接圆不包含点O，转（2）；否则，将当前三角形与临时链表中的所有三角形依次比较，将临时三角形链表中的与当前三角形有公共边的全部删除；并且将当前三角形中的公共边标记出，若当前三角形的公共边数达到3，转（4）（4）将点O与当前三角形的非公共边连接成新的三角形，放入临时三角形链表的尾部，同时从Delaunay三角形链表中删除当前三角形转（2）。（5）判断临时三角形链表是否为空，若否，将临时三角形中的所有三角形全部放入Delaunay三角形链表。3.4程序运行结果以某地形数据为例，程序运行结果如下：（1）所有不规则点：图4 所有不规则点（2）包含所有不规则点的凸包和不规则点集：（上下

24、两条直线向左交于一点D）图5 凸包和不规则点集（3）初始的两个Delaunay三角形：（上下两条直线向左交于一点D）图6 初始的两个Delaunay三角网（4）包含辅助三角形的所有Delaunay三角形：（上下两条直线向左交于一点D）图7 包含辅助三角形的Delaunay三角网（5）去掉辅助点后的Delaunay三角形（即所有不规则点构成的Delaunay三角形网）图8 最终Delaunay三角形网4Super三角改进算法4.1算法基本流程总的算法流程图:图9 算法2流程图和算法1一样，算法2中的第四部分预定义过程(加入离散点)也是整个算法效率的关键，通过第二部分预定义过程对离散点按X方向排序

25、，在第四部中离散点依次插，可以缩小构网区域，只要是最终Delaunay 三角网中的三角形，可以一次性输出，减少遍历次数，并且避免了最终三角形的重复运算，大提高了效率。4.2Super三角形的生成 Super三角形即是一个包含了所有离散点的大三角形，它是根据所有的离散点运算而来，从效率上考虑，Super三角形不能过大，过大的话会增加运算成本，最好能恰好包含所有的离散点。与算法1相比较，可以认为这个三角形就是初始的Delaunay三角形。求解的过程比较简单，运用简单的高中知识就可以解决，以下是生成方法：设vertices数组中包含了所有的离散点。(1)分别用xMax、xMin表示数组vertice

26、s中最大、最小的x坐标，yMax、yMin表示最大、最小的y坐标。用dx ，dy 分别表示（xMax xMin），（ yMax yMin）。(2)为了显得更加安全，建立super三角形时将它构建地稍微大一些：ddx = dx * 0.01ddy = dy * 0.01xMin -= ddx;xMax += ddxdx += 2 * ddxyMin -= ddyyMax += ddydy += 2 * ddy(3)super三角形的三个顶点A，B，C可以表示为：ABC经过以上三步Super三角形便构建完成。4.3详细算法描述第一步：排序，即将所有离散点按 x方向从小到大排列，并用vertices

27、数组存储。第二步：按照上面所述的Super三角形的生成方法求出Super三角形。第三步：不规则点依次插将第二步生成的super三角形放入Delaunay临时链表workset。建立存放边集的临时链表edgeset，初始值为空。输出链表output初始值为空。(1)若没有取到vertices的最后一个点，从vertices中顺序取一个点，称为当前点verticesi，并将边集链表edgeset清空。(2)将临时链表workset中所有“完成”的三角形加入到输出链表output中，并将它们从workset中移除。“完成”的三角形是指外接圆完全在当前点verticesi左边的三角形，也就是组成最终D

28、elaunay三角网的三角形。(3)若临时链表workset不为空，从workset中顺序取一个三角形，称为当前三角形。如果当前三角形的外接圆包含当前点verticesi，则删除当前三角形，但保留它的边，若边集链表edgeset中没有相同的边，将它（它们）加入边集链表edgeset。(4)将当前点与边集链表edgeset中的所有边构成三角形，加入临时链表workset，转（1）。(5)将临时链表workset中的三角形全部加到输出链表output。经过以上步骤，输出链表Output中则包含了所有的Delaunay三角形，将它们输出则形成清晰的三角网图。4.4程序运行结果以某地形数据为例，程序运

29、行结果如下：（1）所有不规则离散点：图10 所有不规则点 （2）包含所有不规则点的Super三角形和不规则点集：图11 不规则点和Super三角形 （3）最终形成的Delaunay三角网：图12 最终形成的Delaunay三角形网 （4）10000个随机点形成的Delaunay三角网:图13一万个随机点形成的Delaunay三角形网4.5面向对象计算机的实现程序大量的涉及到点、边、三角形等元素的运算，它们的设计应该合理和高效，在程序中由于需要实现两个算法，所以对类的设计综合了两个算法的要求，其中有些成员变量和成员函数在一个算法中有用而另一个算法并不会用到，以下为它们的结构设计（部分成员变量和成

30、员函数）：class Point public: void SetX(double PointX); void SetY(double PointY); void SetZ(double PointZ); double GetX(); double GetY(); double GetZ(); private: double x; double y; double z;；class Triangle public: Point GetPointA(); Point GetPointB(); Point GetPointC(); void SetPointA(Point PointA); voi

31、d SetPointB(Point PointB); void SetPointC(Point PointC); BOOL IsInExternalCircle(Point AScatteredPoint); Triangle& operator = (const Triangle& src); BOOL operator = (const Triangle& src);private: Point m_PointC; Point m_PointB; Point m_PointA;4.6测试结果与算法分析实验采用随机生成的数据。当点数从10000到100000变化时，所耗费的时间如下图所示（1.73GHz，存1G），从图中可以看出，构建三角网所需时间和点数几乎成线形增长。将数据做相关分析，回归方程为Y = -12.9428 + 7.329945X，复相关系数高达0.97215678，可见线性相关程度是很