Delaunay算法的实现与应用

分类号 分类号 TP311 1TP311 1 U U D D C C D10621 408 2007 5855 0D10621 408 2007 5855 0 密密 级 公级 公 开开 编编 号 号 20030313032003031303 成成 都都 信信 息息 工工 程程 学学 院院 学学 位位 论论 文文 DelaunayDelaunay 算法的实现与应用算法的实现与应用 论文作者姓名 论文作者姓名 申请学位专业 申请学位专业 计算机科学与技术计算机科学与技术 申请学位类别 申请学位类别 工学学士工学学士 指指导导教教师师姓姓名名 职职称称 论文提交日期 论文提交日期 Delaunay 算法的实现与应用算法的实现与应用 摘摘 要要 数字地形模型是针对地形地貌的一种数字建模 这种建模的结果通常就是一 个数字高程模型 DEM 不规则三角网 TIN 模型是 DEM 中存储和表示非 规则数据的理想模型 它既减少规则网格方法造成的数据冗余 同时在计算效 率方面又优于纯粹基于等高线的方法 所以寻求一种好的 TIN 算法更能快速逼 真的显示与模拟出地貌三维信息 在所有可能的三角网中 狄洛尼 Delaunay 三角网在地形拟合方面表现最为出色 因此常常用于 TIN 的生成 依据 Delaunay 三角剖分准则 直接以边为基础向一侧推进 而不是以凸包为基 础向内推进 从而极大地提高了 Delaunay 三角网推进的速度 仿真实验表明 改进后的算法效率有了显著的提高 关键词关键词 数字地形模型 数字高程模型 不规则三角网 Delaunay 三角 网 Delaunay Triangulation Algorithm Realization 在计算机辅助几何设计中 它可以简化几何模型和实体 模型的数据结构及内部表示 在计算机图形显示中 它可以加速光线跟踪算法和光 照模型的处理以及图形图像插值应用 在地理学中 它可以快速的生成地理模型 Delaunay 三角剖分甚至在无线 Ad Hoc 通信网络中 可以用于构建发现移动节 点间通信路径的在线路由算法 目前 在计算流体力学 统计学 气象学 固 体物理学 计算几何学 图形图像学 三维仿真等多个邻域 都能见到它的踪 影 我们可以举一些简单的应用 近年来 随着网上购物越来越走入人们的生 活 我们可以把某些商品建成可视又可 触 的模型 人们在选购的时候不仅可 以看到商品的样子 还可以通过鼠标和键盘对商品模型进行各种操作 查看它 的弹性 柔软度等性质 从而用更直观 更便捷的方法来了解商品的性能 而 不是对着大篇幅的数字指标发呆 又如 我们可以建立仿真的地理模型 来做 大量的地震 暴雨 等的灾难模拟实验 为灾难预测提供理论依据 从而节省大 量开支 所以 我相信三角网格化的研究工作一定有着良好的发展和应用前景 1 41 4本课题的研究方法本课题的研究方法 三角网格化主要有两种准则 一种称为 Delaunay 三角剖分 即在生成的三角 形网格中 各三角形的最小内角和为最大 另一种是在生成的三角网格中 所有三 角形的边长和最小 其中 Delaunay 三角剖分是目前研究应用最广的一种剖分方 法 本课题的研究方法主要是以 Delaunay 三角网的两个重要性质 空外接圆性质 和最大最小角度性质 以及 Delaunay 三角网的基本原理为基础 参照传统算法 思路 在构建三角网的过程中 改进算法的实现方法 数据结构 以达到提高 效率的目的 2 2Delaunay 方法的基本原理方法的基本原理 2 12 1Voronoi 图与图与 Delaunay 三角网的基本概念三角网的基本概念 Voronoi 图 以下简称 V 图 又叫泰森多 边形或 Dirichlet 图 它是由一组由连接两邻 点直线的垂直平分线组成的连续多边形组成 N 个在平面上有区别的点 按照 第 2 页 共 19 页 最邻近原则划分平面 每个点与它的最近邻区域相关联 Delaunay 三角形是由 与相邻 Voronoi 多边形共享一条边的相关点连接而成的三角形 Delaunay 三角 形的外接圆圆心是与三角形相关的 Voronoi 多边形的一个顶点 Voronoi 三角 形是 Delaunay 图的偶图 为清晰 V 图及 Delaunay 三角网的概念 如图 1 设 为欧几里德平面上的一个点集 则每一个对应区域 n ii P 1 i P 称为点的 Voronoi 区域 各点的 Voronoi 区域 ijnj jii PPPPPV 1 i P 共同组成 V 图 图中用黑色连线 当为非共线点集时 若其对应 V 图中 n ii P 1 有两个点 的 Voronoi 区域有公共边 就连接这两个点 以此类推遍历 i P j P 这 n 个点 可以得到一个连接点集的唯一确定的网格 就是 Delaunay 三 n ii P 1 角网格 图中用红色连线 2 22 2Delaunay 的重要性质的重要性质 空外接圆性质 在由点集 V 生成的 Delaunay 三角网中 每个三角形的外 接圆均不包含该点集的其他任意点 最大最小角度性质 在由点集 V 生成的 Delaunay 三角网中 所有三角 形中的最小角度是最大的 即在生成的三角形网格中 各三角形的最小内角和为 最大 唯一性 不论从区域何处开始构网 最终都将得到一致的结果 由于以上特性 决定了 Delaunay 三角网具有极大的应用价值 Miles 证明 了 Delaunay 三角网是 好的 三角网 Lingas 进一步论证了 在一般情况下 Delauany 三角网是最优的 同时以上特性也成为建立 Delaunay 三角网的重要 算法依据 2 32 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 重复 2 3 步 直到点集中没有在线段 IJ 右侧的点为止 5 将 J 赋给 I J 取其后续点 重复 2 3 4 步 6 当凸包中任意相邻两点连线的右侧不存在离散点时 结束点集凸包求取 第 3 页 共 19 页 过程 第二步 环切边界法凸包三角剖分 在凸包链表中每次寻找一个由相邻两条凸包边组成的三角形 在该三角形 的内部和边界上都不包含凸包上的任何其它点 将这个点去掉后得到新的凸包 链表 重复这个过程 直到凸包链表中只剩三个离散点为止 将凸包链表中的 最后三个离散点构成一个三角形 结束凸包三角剖分过程 第三步 离散点内插 1 找出外接圆包含待插入点的所有三角形 构成插入区域 2 删除插入区域内的三角形公共边 形成由影响三角形顶点构成的多边 形 3 将插入点与多边形所有顶点相连 构成新的 Delaunay 三角形 4 重复 1 2 3 直到所有非凸壳离散点都插入完为止 完成这 一步后 就完成了 Delaunay 三角网的构建 3 3三角剖分改进算法三角剖分改进算法 3 13 1算法基本流程算法基本流程 总的算法流程图 读入离散点坐标 建立凸包 凸包分割初始三角形 加入其余离散点 显示 结束 图 2 算法 1 流程图 算法中的第四部分预定义过程 加入其余离散点 实际包括三个步骤 首先 要确定插入点所在的三角形 然后确定插入点的影响域 最后重构影响域内的三 角网 此步骤是影响效率的主要步骤 在程序运行的时候是被重复执行最多的 步骤 所以算法的时间复杂性主要取决于它们的执行效率 并且 随着离散点 开始 第 4 页 共 19 页 的增加三角形数以点数约 2 倍的速度增加 如何降低以上步骤中的查找和定位 算法的复杂性及缩短算法的执行时间成为提高算法效率的关键 3 23 2Graham 扫描法求凸包扫描法求凸包 凸包的概念 点集 Q 的凸包 convex hull 是指一个最小凸多边形 满足 Q 中的点或者在多边形边上或者在其内 下图中由红色线段表示的多边形就是点 集 Q p0 p1 p12 的凸包 图 3 凸包形成示意图 算法伪代码描述 对于一个有三个或以上点的点集 Q Graham 扫描法的过程 如下 令 p0 为 Q 中 Y X 坐标排序下最小的点 设为对其余点按以 p0 为中心的极角逆时针排序所得的点集 如果有多个点 有相同的极角 除了距 p0 最远的点外全部移除 压 p0 进栈 S 压 p1 进栈 S 压 p2 进栈 S for i 3 to m do while 由 S 的栈顶元素的下一个元素 S 的栈顶元素以及 pi 构成的折 线段不拐向左侧 对 S 弹栈 压 pi 进栈 S return S 此过程执行后 栈 S 由底至顶的元素就是 Q 的凸包顶点按逆时针排列的点 序列 需要注意的是 我们对点按极角逆时针排序时 并不需要真正求出极角 只需要求出任意两点的次序就可以了 而这个步骤可以用矢量叉积性质实现 3 33 3详细算法描述详细算法描述 算法基于上述的传统构建算法 但仅有两步 第 5 页 共 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 三角形 生成凸包的过程结束跳过一下各步 否 则 继续第四步 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 三角形网的基础上 将其余离散点依次内插 形成 Delaunay 三角网 该算法原理遵循传统 TSAI 离散点内插算法 内插的计算机实现 将第一步中形成的初始 Delaunay 三角形放入 Delaunay 三角形链表 并建 立临时三角形链表 放置新生成的三角形 初始值为空 1 若没到点集链表的尾端 按顺序取不规则点中的一个点 O 将临时 三角形链表清空 2 若 Delaunay 三角形链表不为空 从该链表中按顺序取一个 称为 当前三角形 若为空 转 5 3 若当前三角形的外接圆不包含点 O 转 2 否则 将当前三角形 与临时链表中的所有三角形依次比较 将临时三角形链表中的与当前三角形有 公共边的全部删除 并且将当前三角形中的公共边标记出 若当前三角形的公 共边数达到 3 转 4 4 将点 O 与当前三角形的非公共边连接成新的三角形 放入临时三角 形链表的尾部 同时从 Delaunay 三角形链表中删除当前三角形转 2 5 判断临时三角形链表是否为空 若否 将临时三角形中的所有三角 第 6 页 共 19 页 形全部放入 Delaunay 三角形链表 3 43 4程序运行结果程序运行结果 以某地形数据为例 程序运行结果如下 1 所有不规则点 图 4 所有不规则点 2 包含所有不规则点的凸包和不规则点集 上下两条直线向左交于 一点 D 图 5 凸包和不规则点集 3 初始的两个 Delaunay 三角形 上下两条直线向左交于一点 D 图 6 初始的两个 Delaunay 三角网 4 包含辅助三角形的所有 Delaunay 三角形 上下两条直线向左交于 一点 D 第 7 页 共 19 页 图 7 包含辅助三角形的 Delaunay 三角网 5 去掉辅助点后的 Delaunay 三角形 即所有不规则点构成的 Delaunay 三角形网 图 8 最终 Delaunay 三角形网 4 4SuperSuper 三角改进算法三角改进算法 4 14 1算法基本流程算法基本流程 总的算法流程图 读入离散点坐标 排序离散点 构建 Super 三角形 加入离散点 显示 结束 图 9 算法 2 流程图 开始 第 8 页 共 19 页 和算法 1 一样 算法 2 中的第四部分预定义过程 加入离散点 也是整个算 法效率的关键 通过第二部分预定义过程对离散点按 X 方向排序 在第四部中 离散点依次内插 可以缩小构网区域 只要是最终 Delaunay 三角网中的三角 形 可以一次性输出 减少遍历次数 并且避免了最终三角形的重复运算 大 大地提高了效率 4 24 2Super 三角形的生成三角形的生成 Super 三角形即是一个包含了所有离散点的大三角形 它是根据所有的离 散点运算而来 从效率上考虑 Super 三角形不能过大 过大的话会增加运算 成本 最好能恰好包含所有的离散点 与算法 1 相比较 可以认为这个三角形 就是初始的 Delaunay 三角形 求解的过程比较简单 运用简单的高中知识就可 以解决 以下是生成方法 设 vertices 数组中包含了所有的离散点 1 分别用 xMax xMin 表示数组 vertices 中最大 最小的 x 坐标 yMax yMin 表示最大 最小的 y 坐标 用 dx dy 分别表示 xMax xMin yMax yMin 2 为了显得更加安全 建立 super 三角形时将它构建地稍微大一些 ddx dx 0 01 ddy dy 0 01 xMin ddx xMax ddx dx 2 ddx yMin ddy yMax ddy dy 2 ddy 3 super 三角形的三个顶点 A B C 可以表示为 A 0 3 3 yMindyxMin B 0 3 3 yMindyxMax C 5 0 3 5 0 dxxMaxxMaxxMin 经过以上三步 Super 三角形便构建完成 4 34 3详细算法描述详细算法描述 第一步 排序 即将所有离散点按 x 方向从小到大排列 并用 vertices 数组存储 第二步 按照上面所述的 Super 三角形的生成方法求出 Super 三角形 第三步 不规则点依次内插 第 9 页 共 19 页 将第二步生成的 super 三角形放入 Delaunay 临时链表 workset 建立存放边集的临时链表 edgeset 初始值为空 输出链表 output 初始值为空 1 若没有取到 vertices 的最后一个点 从 vertices 中顺序取一 个点 称为当前点 vertices i 并将边集链表 edgeset 清空 2 将临时链表 workset 中所有 完成 的三角形加入到输出链表 output 中 并将它们从 workset 中移除 完成 的三角形是指外接圆完全在 当前点 vertices i 左边的三角形 也就是组成最终 Delaunay 三角网的三角形 3 若临时链表 workset 不为空 从 workset 中顺序取一个三角形 称 为当前三角形 如果当前三角形的外接圆包含当前点 vertices i 则删除当 前三角形 但保留它的边 若边集链表 edgeset 中没有相同的边 将它 它们 加入边集链表 edgeset 4 将当前点与边集链表 edgeset 中的所有边构成三角形 加入临时链 表 workset 转 1 5 将临时链表 workset 中的三角形全部加到输出链表 output 经过以上步骤 输出链表 Output 中则包含了所有的 Delaunay 三角形 将 它们输出则形成清晰的三角网图 4 44 4程序运行结果程序运行结果 以某地形数据为例 程序运行结果如下 1 所有不规则离散点 图 10 所有不规则点 2 包含所有不规则点的 Super 三角形和不规则点集 第 10 页 共 19 页 图 11 不规则点和 Super 三角形 3 最终形成的 Delaunay 三角网 图 12 最终形成的 Delaunay 三角形网 4 10000 个随机点形成的 Delaunay 三角网 图 13 一万个随机点形成的 Delaunay 三角形网 4 54 5面向对象计算机的实现面向对象计算机的实现 程序大量的涉及到点 边 三角形等元素的运算 它们的设计应该合理和 高效 在程序中由于需要实现两个算法 所以对类的设计综合了两个算法的要 求 其中有些成员变量和成员函数在一个算法中有用而另一个算法并不会用到 以下为它们的结构设计 部分成员变量和成员函数 class Point public void SetX double PointX void SetY double PointY void SetZ double PointZ double GetX double GetY double GetZ private 第 11 页 共 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 Triangle BOOL operator const Triangle private Point m PointC Point m PointB Point m PointA 4 64 6测试结果与算法分析测试结果与算法分析 实验采用随机生成的数据 当点数从 10000 到 100000 变化时 所耗费的时 间如下图所示 1 73GHz 内存 1G 从图中可以看出 构建三角网所需时间和 点数几乎成线形增长 将数据做相关分析 回归方程为 Y 12 9428 7 329945X 复相关系数高达 0 97215678 可见线性相关程度是很高的 所以 从统计意义上讲 此算法还是比较成功的 第 12 页 共 19 页 计算耗时统计 0 10 20 30 40 50 60 70 80 12345678910 点集个数 万个 所耗时间 秒 图 14 耗时统计图 更详细的数据请看下表 表 1 构网耗时详细数据表 构网点数 个 构网时间 秒 三角形个数 个 100001 73119875 200004 6839606 300008 47159173 4000013 05878604 5000018 61197874 6000025 225116960 7000033 416136014 8000043 353154931 9000055 442173641 10000069 732192239 5 5Delaunay 算法的应用算法的应用 5 15 1插值基本原理插值基本原理 其实插值的原理非常简单 如下图所示 第 13 页 共 19 页 A B C M M 图 15 三角插值原理示意图 若 M x y 所在的三角形为 ABC 三顶点的坐标为 x1 y1 z1 x2 y2 z2 x3 y3 z3 则由 A B C 确定的平面方程为 或者 1 1 1 1 1 333 222 111 zyx zyx zyx zyx 0 131313 121212 111 zzyyxx zzyyxx zzyyxx 令 1221 1221 1221 zzz yyy xxx 1331 1331 1331 zzz yyy xxx 则 21313121 213131211213131211 1 yxyx xzxzyyzyzyxx zzM A B C 三点确定一个平面 在平面上任意一点 M 的值 可由这个平面方程 得出 在 Delanunay 三角网内的点都可以通过它所在三角形进行插值 若不 在三角网内则可以通过延伸三角形同样进行插值 这样就可以完成整个范围 内的插值了 5 25 2笔者源程序笔者源程序 已知三点求在这三点确定的平面上的任意一点 M 的值 double CalculatePointHeight Point M Point M1 Point M2 Point M3 double a 2 3 a 0 0 M2 x M1 x a 0 1 M2 y M1 y 第 14 页 共 19 页 a 0 2 M2 z M1 z a 1 0 M3 x M1 x a 1 1 M3 y M1 y a 1 2 M3 z M1 z double A B C D A B C 平面法向量 A a 0 1 a 1 2 a 0 2 a 1 1 B a 0 2 a 1 0 a 0 0 a 1 2 C a 0 0 a 1 1 a 0 1 a 1 0 D A M1 x B M1 y C M1 z double x D A M x B M y C return x 5 35 3基于网格插值的等值线生成基于网格插值的等值线生成 等值线是某地地理现象数值相等的各点的连线 等值线图是用布满一定区 域内的若干条等值线表示某地理现象的数量分布状况 由于等值线标有数值 而且数值间隔是相同的 所以可以根据等值线的数值大小 疏密程度 排列的 方向 形状变化等反映出该地理事物变化的急缓 递变的方向及分布特点等 这里只简单描述算法基本原理 等值线跟踪是在已知格网点的基础上 内插出等值线点 然后跟踪等值点 形成封闭或非封闭等值线 非封闭等值线的端点必然落在边界线上 而边界线 必然是封闭的 通过跟踪边界线 再把非封闭等值线端点插入边界线 通过边 界线建立等值线间的拓扑关系 在此基础上跟踪出封闭多边形 实现等值线的 填充 该算法分为以下几大步骤 边界线的跟踪 非封闭等值线端点插入边界 线中并排序 非封闭等值线建立拓扑关系 跟踪封闭多边形并排序 等值线的 充填 以下为程序运行效果图 第 15 页 共 19 页 图 16 插值效果与等值线图 结结 论论 Delaunay 三角剖分是国际三大流行的全自动网格生成方法之一 它的最主 要特点之一就是它自动避免了生成小内角的长薄单元 经过几十年不懈的研究 在二维欧氏空间上的 Delaunay 自动网格生成技术已经趋于完善 理论的完备性 和实践的多样化均给这一有限元网格剖分方法以新的活力 时至今日 如何评 价一种 Delaunay 网格剖分程序的优劣 该程序的时间 空间复杂度以及强壮性 仍是人们热切关注的话题 同时 针对 D