Delaunay算法的实现与应用.doc

1、分类号：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 在计算机辅助几何设计中 ,它可以简化几何模型和

3、实体模型的数据结构及内部表示;在计算机图形显示中, 它可以加速光线跟踪算法和光照模型的处理以及图形图像插值应用;在地理学中, 它可以快速的生成地理模型。Delaunay 三角剖分甚至在无线 Ad Hoc 通信网络中，可以用于构建发现移动节点间通信路径的在线路由算法。目前，在计算流体力学、统计学、气象学、固体物理学、计算几何学、图形图像学、三维仿真等多个邻域，都能见到它的踪影。我们可以举一些简单的应用：近年来，随着网上购物越来越走入人们的生活，我们可以把某些商品建成可视又可“触” 的模型，人们在选购的时候不仅可以看到商品的样子，还可以通过鼠标和键盘对商品模型进行各种操作，查看它的弹性、柔软度等性

4、质，从而用更直观、更便捷的方法来了解商品的性能，而不是对着大篇幅的数字指标发呆。又如，我们可以建立仿真的地理模型，来做大量的地震、暴雨 等的灾难模拟实验，为灾难预测提供理论依据，从而节省大量开支。所以，我相信三角网格化的研究工作一定有着良好的发展和应用前景。1.4 本课题的研究方法三角网格化主要有两种准则:一种称为 Delaunay 三角剖分,即在生成的三角形网格中,各三角形的最小内角和为最大;另一种是在生成的三角网格中 ,所有三角形的边长和最小.其中, Delaunay 三角剖分是目前研究应用最广的一种剖分方法.本课题的研究方法主要是以 Delaunay 三角网的两个重要性质（空外接圆性质和

5、最大最小角度性质）以及 Delaunay 三角网的基本原理为基础，参照传统算法思路，在构建三角网的过程中，改进算法的实现方法，数据结构，以达到提高效率的目的。2 Delaunay 方法的基本原理2.1 Voronoi 图与 Delaunay 三角网的基本概念Voronoi图（以下简称 V图） ，又叫泰森多边形或 Dirichlet图，它是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成。N 个在平面上有区别的点，按照第 3 页 共 19 页最邻近原则划分平面；每个点与它的最近邻区域相关联。Delaunay 三角形是由与相邻 Voronoi多边形共享一条边的相关点连接而成的三角形。Dela

6、unay 三角形的外接圆圆心是与三角形相关的 Voronoi多边形的一个顶点。Voronoi 三角形是 Delaunay图的偶图。为清晰 V图及 Delaunay三角网的概念，如图 1，设为欧几里德平面上的一个点集，则每一个 对应区域niP1 iP称为点 的 Voronoi区域，各点的 Voronoi区域ijnjiPV,: i共同组成 V图（图中用黑色连线） 。当 为非共线点集时，若其对应 V图中n1有两个点 ， 的 Voronoi区域有公共边，就连接这两个点，以此类推遍历ij这 n个点，可以得到一个连接点集 的唯一确定的网格，就是 Delaunay三ni1角网格（图中用红色连线） 。2.2

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

8、依据。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 之间，并将

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

10、所有三角形，构成插入区域。（2）删除插入区域内的三角形公共边，形成由影响三角形顶点构成的多边形。（3）将插入点与多边形所有顶点相连，构成新的 Delaunay 三角形。（4）重复（1） 、 （2） 、 （3） ，直到所有非凸壳离散点都插入完为止。完成这一步后，就完成了 Delaunay 三角网的构建。3 三角剖分改进算法3.1 算法基本流程总的算法流程图:读入离散点坐标建立凸包凸包分割初始三角形加入其余离散点显示结束图 2 算法 1 流程图算法中的第四部分预定义过程(加入其余离散点)实际包括三个步骤：首先要确定插入点所在的三角形，然后确定插入点的影响域,最后重构影响域内的三角网。此步骤是影响效

11、率的主要步骤，在程序运行的时候是被重复执行最多的步骤，所以算法的时间复杂性主要取决于它们的执行效率。并且，随着离散点开始第 5 页 共 19 页的增加三角形数以点数约 2倍的速度增加，如何降低以上步骤中的查找和定位算法的复杂性及缩短算法的执行时间成为提高算法效率的关键。3.2 Graham 扫描法求凸包凸包的概念：点集 Q的凸包(convex hull)是指一个最小凸多边形，满足 Q中的点或者在多边形边上或者在其内。下图中由红色线段表示的多边形就是点集 Q=p0,p1,.p12的凸包。 图 3 凸包形成示意图算法伪代码描述:对于一个有三个或以上点的点集 Q，Graham 扫描法的过程如下： 令

12、 p0为 Q中 Y-X坐标排序下最小的点 设为对其余点按以 p0为中心的极角逆时针排序所得的点集（如果有多个点有相同的极角，除了距 p0最远的点外全部移除压 p0进栈 S压 p1进栈 S压 p2进栈 Sfor i 3 to mdo while 由 S的栈顶元素的下一个元素、S 的栈顶元素以及 pi构成的折线段不拐向左侧对 S弹栈压 pi进栈 Sreturn S;此过程执行后，栈 S由底至顶的元素就是 Q的凸包顶点按逆时针排列的点序列。需要注意的是，我们对点按极角逆时针排序时，并不需要真正求出极角，只需要求出任意两点的次序就可以了。而这个步骤可以用矢量叉积性质实现。 3.3 详细算法描述算法基于

13、上述的传统构建算法，但仅有两步：第 6 页 共 19 页第一步：（1） 在离散点集中寻找一纵坐标最小的点 A。（2） 以点 A为起点，寻找两个点 B、D，使得向量 AB与横坐标轴夹角最小，向量 AD与横坐标轴夹角最大。若点 A、B、D 共线，将原 B点标记为 A，寻找点 D，使得向量 AD与直线 AB夹角最大；寻找点 C使得向量 BC与线段 AB夹角最小。否则，若 A、B、D 不共线，则寻找点 C使得向量 BC与线段 AB夹角最小。这样，所有点都在逆时针旋转的折线 DABC的左侧。（3） 上面一步生成的点 C、D 如果为同一点，则ABC（或ABD）即为包含所有不规则点的 Delaunay三角形

14、，生成凸包的过程结束跳过一下各步；否则，继续第四步。（4） 寻找一个距离直线 AB最远的点 E，过 E作直线 AB的平行线，与射线 AD、BC 分别交于点 D、C ，将点 D、C重新标记为点 D、C，则凸四边形DABC包括所有不规则点（5） 判断ABC 的外接圆是否包含点 D，若是则求ABC 的外接圆半径R，在射线 BC上距点 C2.1R远处取一点 D，并将点 D重新标记为点 D，则凸四边形 DCBA即为所求得的初始凸包。若ABC 的外接圆不包含点 D，则凸四边形 DCBA即为所求得的初始凸包。第二步：不规则点内插在原有 Delaunay三角形网的基础上，将其余离散点依次内插，形成Delaun

15、ay三角网。该算法原理遵循传统（TSAI）离散点内插算法。内插的计算机实现：将第一步中形成的初始 Delaunay三角形放入 Delaunay三角形链表。并建立临时三角形链表，放置新生成的三角形，初始值为空。（1） 若没到点集链表的尾端，按顺序取不规则点中的一个点 O；将临时三角形链表清空。（2） 若 Delaunay三角形链表不为空，从该链表中按顺序取一个，称为当前三角形。若为空，转（5） 。（3） 若当前三角形的外接圆不包含点 O，转（2） ；否则，将当前三角形与临时链表中的所有三角形依次比较，将临时三角形链表中的与当前三角形有公共边的全部删除；并且将当前三角形中的公共边标记出，若当前三角

16、形的公共边数达到 3，转（4）（4） 将点 O与当前三角形的非公共边连接成新的三角形，放入临时三角形链表的尾部，同时从 Delaunay三角形链表中删除当前三角形转（2） 。（5） 判断临时三角形链表是否为空，若否，将临时三角形中的所有三角第 7 页 共 19 页形全部放入 Delaunay三角形链表。3.4 程序运行结果以某地形数据为例，程序运行结果如下：（1） 所有不规则点：图 4 所有不规则点（2） 包含所有不规则点的凸包和不规则点集：（上下两条直线向左交于一点 D）图 5 凸包和不规则点集（3） 初始的两个 Delaunay 三角形：（上下两条直线向左交于一点 D）图 6 初始的两个

17、Delaunay 三角网（4） 包含辅助三角形的所有 Delaunay 三角形：（上下两条直线向左交于一点 D）第 8 页 共 19 页图 7 包含辅助三角形的 Delaunay 三角网（5） 去掉辅助点后的 Delaunay 三角形（即所有不规则点构成的Delaunay 三角形网）图 8 最终 Delaunay 三角形网4 Super三角改进算法4.1 算法基本流程总的算法流程图:读入离散点坐标排序离散点构建 Super 三角形加入离散点显示结束图 9 算法 2 流程图开始第 9 页 共 19 页和算法 1一样，算法 2中的第四部分预定义过程(加入离散点)也是整个算法效率的关键，通过第二部分

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

19、(1) 分别用 xMax、xMin 表示数组 vertices中最大、最小的 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 可以表示为：A ),0.3/(yMindyxinB aC )5.0

20、*3,5.*)dxaxi 经过以上三步 Super三角形便构建完成。4.3 详细算法描述第一步：排序，即将所有离散点按 x方向从小到大排列，并用 vertices数组存储。第二步：按照上面所述的 Super三角形的生成方法求出 Super三角形。第三步：不规则点依次内插第 10 页 共 19 页 将第二步生成的 super三角形放入 Delaunay临时链表 workset。 建立存放边集的临时链表 edgeset，初始值为空。 输出链表 output初始值为空。(1) 若没有取到 vertices的最后一个点，从 vertices中顺序取一个点，称为当前点 verticesi，并将边集链表

21、edgeset清空。(2) 将临时链表 workset中所有“完成”的三角形加入到输出链表output中，并将它们从 workset中移除。 “完成”的三角形是指外接圆完全在当前点 verticesi左边的三角形，也就是组成最终 Delaunay三角网的三角形。(3) 若临时链表 workset不为空，从 workset中顺序取一个三角形，称为当前三角形。如果当前三角形的外接圆包含当前点 verticesi，则删除当前三角形，但保留它的边，若边集链表 edgeset中没有相同的边，将它（它们）加入边集链表 edgeset。(4) 将当前点与边集链表 edgeset中的所有边构成三角形，加入临时

22、链表 workset，转（1） 。(5) 将临时链表 workset中的三角形全部加到输出链表 output。经过以上步骤，输出链表 Output中则包含了所有的 Delaunay三角形，将它们输出则形成清晰的三角网图。4.4 程序运行结果以某地形数据为例，程序运行结果如下：（1） 所有不规则离散点：图 10 所有不规则点 （2） 包含所有不规则点的 Super 三角形和不规则点集：第 11 页 共 19 页图 11 不规则点和 Super 三角形 （3） 最终形成的 Delaunay 三角网：图 12 最终形成的 Delaunay 三角形网 （4） 10000 个随机点形成的 Delauna

23、y 三角网:图 13 一万个随机点形成的 Delaunay 三角形网4.5 面向对象计算机的实现程序大量的涉及到点、边、三角形等元素的运算，它们的设计应该合理和高效，在程序中由于需要实现两个算法，所以对类的设计综合了两个算法的要求，其中有些成员变量和成员函数在一个算法中有用而另一个算法并不会用到，以下为它们的结构设计（部分成员变量和成员函数）：class Pointpublic:void SetX(double PointX);void SetY(double PointY);void SetZ(double PointZ);double GetX();double GetY();double

24、 GetZ();private:第 12 页 共 19 页double x;double y;double z;；class Triangle public:Point GetPointA();Point GetPointB();Point GetPointC();void SetPointA(Point PointA);void SetPointB(Point PointB);void SetPointC(Point PointC);BOOL IsInExternalCircle(Point AScatteredPoint);TriangleBOOL operator = (const Tr

25、iangleprivate:Point m_PointC;Point m_PointB;Point m_PointA;4.6 测试结果与算法分析实验采用随机生成的数据。当点数从 10000到 100000变化时，所耗费的时间如下图所示（1.73GHz，内存 1G） ，从图中可以看出，构建三角网所需时间和点数几乎成线形增长。将数据做相关分析，回归方程为 Y = -12.9428 + 7.329945X，复相关系数高达 0.97215678，可见线性相关程度是很高的。所以从统计意义上讲，此算法还是比较成功的。第 13 页 共 19 页计 算 耗 时 统 计010203040506070801 2

26、3 4 5 6 7 8 9 10点 集 个 数 (万 个 )所耗时间(秒)图 14 耗时统计图更详细的数据请看下表：表 1 构网耗时详细数据表构网点数（个）构网时间（秒）三角形个数（个）10000 1.731 1987520000 4.68 3960630000 8.471 5917340000 13.058 7860450000 18.611 9787460000 25.225 11696070000 33.416 13601480000 43.353 15493190000 55.442 173641100000 69.732 1922395 Delaunay 算法的应用5.1 插值基本原

27、理其实插值的原理非常简单，如下图所示：第 14 页 共 19 页ABCMM图 15 三角插值原理示意图若 M（x, y）所在的三角形为ABC，三顶点的坐标为（x1, y1, z1） ， （x2, y2, z2） ， （x3, y3, z3） ，则由 A，B，C 确定的平面方程为或者 133221zyxz 0131313222 zyx令;121zy13zyx则 21312 2131321 )()(. yxxzxMA、B、C 三点确定一个平面，在平面上任意一点 M的值，可由这个平面方程得出。在 Delanunay三角网内的点都可以通过它所在三角形进行插值，若不在三角网内则可以通过延伸三角形同样进行

28、插值。这样就可以完成整个范围内的插值了。5.2 笔者源程序已知三点求在这三点确定的平面上的任意一点 M的值double CalculatePointHeight(Point M,Point M1,Point M2,Point M3)double a23;a00 = M2.x - M1.x;a01 = M2.y - M1.y;第 15 页 共 19 页a02 = M2.z - M1.z;a10 = M3.x - M1.x;a11 = M3.y - M1.y;a12 = M3.z - M1.z;double A , B , C , D;/A B C平面法向量A = a01 * a12 - a02

29、* a11;B = a02 * a10 - a00 * a12;C = a00 * a11 - a01 * a10;D = A * M1.x + B * M1.y + C * M1.z;double x = (D - A * M.x - B * M.y) / C;return x;5.3 基于网格插值的等值线生成等值线是某地地理现象数值相等的各点的连线。等值线图是用布满一定区域内的若干条等值线表示某地理现象的数量分布状况，由于等值线标有数值，而且数值间隔是相同的，所以可以根据等值线的数值大小、疏密程度、排列的方向、形状变化等反映出该地理事物变化的急缓、递变的方向及分布特点等。这里只简单描述算法

30、基本原理: 等值线跟踪是在已知格网点的基础上，内插出等值线点，然后跟踪等值点，形成封闭或非封闭等值线。非封闭等值线的端点必然落在边界线上，而边界线必然是封闭的，通过跟踪边界线，再把非封闭等值线端点插入边界线，通过边界线建立等值线间的拓扑关系，在此基础上跟踪出封闭多边形，实现等值线的填充。该算法分为以下几大步骤：边界线的跟踪，非封闭等值线端点插入边界线中并排序，非封闭等值线建立拓扑关系，跟踪封闭多边形并排序，等值线的充填。以下为程序运行效果图:第 16 页 共 19 页图 16 插值效果与等值线图结 论Delaunay 三角剖分是国际三大流行的全自动网格生成方法之一，它的最主要特点之一就是它自动

31、避免了生成小内角的长薄单元。经过几十年不懈的研究，在二维欧氏空间上的 Delaunay 自动网格生成技术已经趋于完善，理论的完备性和实践的多样化均给这一有限元网格剖分方法以新的活力。时至今日，如何评价一种 Delaunay 网格剖分程序的优劣，该程序的时间、空间复杂度以及强壮性仍是人们热切关注的话题。同时，针对 Delaunay 三角网格的不同评价标准，人们仍会有设计出许多新的方法来更好的适应实际的需求。另外，在向三维欧氏空间拓展的形势下，Delaunay 三角剖分还有极大的发展空间。就本设计的以上两种算法示例方法而言仍有许多可以改进的地方，如点集加密技术的讨论、优化方法设计的时间与空间复杂性

32、研究、针对各类实际数据的模拟和剖分状态比较、与其他网格剖分技术的对比研究等等，这里我所给出的仅是一种 Delaunay 三角剖分的尝试，还有许多后续工作可以继续讨论。总之，Delaunay 三角网格自动剖分技术作为一种成熟的有限元网格剖分方法，在实践中有着极其广泛的应用，因此对它的研究也是具有实际价值、意义深远的科研工作。参考文献1 胡金星,马照亭,吴焕萍,潘懋.基于格网划分的海量数据Delaunay三角剖分J.测绘学报,2004,33(2):163-167。http:/www.biyezuopin.cc2 李志林,朱庆.数字高程模型M.武汉:武汉大学出版社,2001。第 17 页 共 19

33、页3 吴立新,史文中.地理信息系统原理与算法M.北京:科学出版社,2003。4 杨杰.基于凹凸顶点判定单多边形的三角剖分J.小型微型计算机系统,2000,21(9):974-975。http:/5 周培德.计算机几何算法分析与设计M.清华大学出版社、广西科学技术出版社,2000。http:/www.17chaogu8.wang6 候俊杰.深入浅出WINDOWS MFC程序设计M.华中理工大学出版社,1998。7 杨钦.限定Delaunay三角网格剖分技术M.电子工业出版社,2005。8 邬吉明,沈隆钧,张景琳.Delaunay三角网格的一种快速生成法J.数值计算与计算机应用,2001(7):2

34、67-275。9 胡恩球,张新访等.有限元网格生成方法发展综述J.计算机辅助设计与图形学学报，1997,9(4):378-383。第 18 页 共 19 页致 谢本文是在老师的热情关心和指导下完成的，他渊博的知识和严谨的治学作风使我受益匪浅，对顺利完成本课题起到了极大的作用。在此向他表示我最衷心的感谢！感谢金虎老师在课题研究初期给予的帮助，在金虎老师的帮助下课题才得以很好的展开，有了一个很好的开端！还要感谢四川省地质调查院主任，有了他的细心提点，课题中的一些重要算法实现才得以顺利进行。在论文完成过程中，本人还得到了其他老师和许多同学的热心帮助，本人向他们表示深深的谢意！最后向在百忙之中评审本文的各位专家、老师表示衷心的感谢！