基于GIS的高质量约束Delaunay三角网格剖分

第26卷 第5期2010年9月地理与地理信息科学
Geog ra phy and Geo-Infor matio n Science V ol.26 N o.5September 2010
收稿日期:2010-04-21; 修订日期:2010-06-10
基金项目:国家自然科学基金重点项目(50839001);国家自然科学基金项目(50874021、50779006);辽宁省高等学校科研项目计划(L20100321) 作者简介:赵晓东(1969-),男,博士,教授,从事GIS 在采矿、岩土工程应用和水动模型耦合的研究。

E-mail:xdong.zhao@
基于GIS 的高质量约束Delaunay 三角网格剖分
赵晓东1
,晏小宝1
,沈永明2
,王 亮
2
(1.大连大学院士创业园中日地层环境科学研究中心,辽宁大连116622;2.大连理工大学海岸和近海工程国家重点实验室,辽宁大连116023)
摘要:在分析现有非结构化网格剖分算法的基础上,提出了一种GIS 支持下的改进分治算法实现约束Delaunay 三角网格剖分。

该方法利用了GIS 的空间拓扑关系对算法输入数据进行预处理,基于三角形的统一数据结构实现了网格细化,对输出剖分网格进行准确的拓扑和约束条件的检查,并基于推进阵面算法思想,结合空间邻近拓扑关系实现了三角剖分节点和网格的重新编号,方便了实际问题中开边界条件的赋值,提高了计算效率。

实例应用表明,该方法大大简化了数值模型非结构化网格剖分的前处理过程,集成了几种综合算法的优点,在保证原分治算法时间复杂度的基础上,提高了约束条件下Delauna y 三角网格生成的质量。

关键词:网格剖分;GIS;约束Delaunay 三角剖分
中图分类号:P208 文献标识码:A 文章编号:1672-0504(2010)05-0024-05
0 引言
在应用有限元、有限差分和有限体积法对力学问题进行数值计算的前处理中,网格自动剖分占有重要的位置。

目前,实现网格剖分的算法很多,而且不断有新的算法推出,其相关研究领域不仅针对数值模拟计算,对GIS 数据表达、地学分析、计算机视觉、表面目标重构等众多领域也是一项重要的应用技术[1-3]。

网格剖分可分为结构化网格(Structured Gr id)和非结构化网格(Unstructured Grid)两种。

非结构化网格易于控制网格单元大小、形状及网格点位置,可以实现合理分布网格的密度,提高计算精度,因此具有比结构化网格更大的灵活性和对复杂边界更强的适应性。

生成二维非结构化网格的常用方法[4]有四叉树法(Quadtr ee)、Delaunay 三角剖分法(Delaunay T riangulation,DT)和推进阵面法(Ad -v ancing Fr ont M ethod,AFM)以及几种方法的综合和改进。

由于Delaunay 三角网是Vor ono i 图的直线对偶图,具有空外接圆和最大最小角特性,还可以尽可能地避免病态三角形的出现,实现约束条件下的三角剖分,因此,在数值计算非结构化网格剖分和GIS 三角网(T IN)的数据格式表达中广泛使用。

目前常见的构建Delaunay 三角网的算法有:分治算法(Div ide -and -conquer )、逐点插入算法(Incre -m ental Insertio n)、生长算法和扫描线算法(Sw eep -line)[3,5-7],其中以分治算法效率最高,时间复杂度
为O(N logN)[5]。

由于实际的不同需要,除几何约束条件外,对三角剖分的角度、大小和节点物理量的控制也都有相应约束,目前的剖分算法尚不能满足以上全部约束条件进行剖分。

本文基于改进的分治
算法,提供了梯度变化和针对实际问题的单元和节点重编号,并利用GIS 强大的空间数据处理和分析功能对算法做数据前后处理,保证了高质量约束Delaunay 三角网格剖分的生成。

1 约束Delaunay 三角化方法
分治算法[5]的基本思路是把输入点集分割为数个较小的点集,在各个子点集内生成小三角网,再逐级合并的一个递归过程。

子点集最终只有两点或者三点,形成一条边、共线边或者三角形。

算法必须将点集连接为凸区域(凸壳),对凹区域和多连通区域会产生域外三角形,无法保证约束边的存在。

约束Delaunay 三角化的分治算法的关键过程及实现方法如下。

1.1 点集划分
原分治算法只有Y 轴方向的分割,这里改进为交替分割[8],即点域X 轴方向的长度大于Y 轴方向的长度,则以X 轴方向对半分割点集,否则以Y 轴方向对半分割点集。

该过程为递归调用的分割过程。

1.2 约束特征线重构
在Delaunay 三角剖分上,采用逐次加入约束特
征线的方法进行约束重构,可以证明该重构方法是收敛的[4],具体步骤如下:1)对于任一条约束直线,查找与之相交的三角形,即该约束线穿越的所有三角形组成的影响域;2)删除该影响域内的所有三角形,形成一个多边形空腔,称为恢复域;3)以该约束直线为界将恢复域一分为二,形成两个多边形子域;4)对这两个多边形子域进行Delaunay 三角剖分,得到包含约束直线的影响域内的三角剖分;5)重复上述步骤,直至所有的约束线都被加入。

1.3 基于三角形的统一数据结构
Delaunay 三角剖分涉及三角形和线段两类几何数据,如果将线段表示为退化的三角形势必增加算法特殊情况判断的冗余。

如图1a 所示,利用/影子0层表示完整的三角形结构,凸壳上的每一线段有一个顶点在无穷远处(算法中采用空指针)的/影子0三角形,特征线段有两个/影子0三角形。

图1b 中下端虚线三角形为子网合并中新建的/影子0三角形,阴影的三角形采用局部LOP 优化算法[7],保证三角形的Delaunay
特性。

图1 分治算法的三角形数据结构和子网合并[9]
Fig.1 Data s tructure of triangle and halfway through a merge
step of the dividing and conqu ering algorithm
1.4 网格细化修正
剖分细化是在网格中插入顶点细化网格,使得
网格最大和最小角满足约束。

本文采用Shew chuk 的细化算法[9,10],它是对Ruppert 和Chew 提出算法[11,12]的综合改进算法,顶点插入使用LOP 算法[7],保证Delaunay 特性,遵循两个规则:
规则1:一条线段的直径圆是唯一且最小的包含该线段的圆。

如图2所示,如果有一点位于直径园内,则称该线段被侵入。

任何被侵入线段在其中点
位置插入顶点进行分割。

图2 线段被递归细分直至不被侵入为止[9]
Fig.2 Segments are split recursively until no
segments are encroached
规则2:一个三角形的外接圆是唯一经过三角形3个顶点的圆。

如果三角形的一个角度过小或面积过大的情况下满足了约束条件,则称该三角形为坏三角形。

如图3所示,以外心顶点(外接圆的中心)作为插入点进行分割。

图3 每个坏三角形线在外心顶点插入分割[9]
Fig.3 Each bad triangle is split by inserting a vertex
at its circumcenter
剖分细化具体算法如下:1)取约束线段集合中的每一条约束线段,检查它是否被其他点侵入,是则转第2步;否则转第3步。

2)按规则1将约束线段分割,将新的约束线段加入约束线段集合,转第1步。

3)取剖分三角形集合的三角形进行检查,全部合格,转第5步;否则,对于不满足约束条件的三角形按规则2插入点。

4)如果插入点不侵入约束线段,将该点插入网格中,转第3步;否则,转第2步。

5)算法收敛停止。

2 GIS 支持下约束Delaunay 三角网格剖分
高质量Delaunay 三角网格意味着必须满足一定的约束条件,而约束条件是根据工程需要或者数值模拟的条件给定的,如流体动力学问题一般在边
界、结构物和扩散源附近采取网格加密方法,这些都要对输入点集进行预处理;输出的剖分三角形内角不能太大也不能太小,面积梯度变化适中。

因此,要实现本文提出的非结构网格生成算法,需要提取计算区域边界点的坐标,包括外边界和内边界上各节点光滑、间距的处理、内部约束特征线的设置等。

如图4所示,这些预处理都可以在GIS 的支持下完成,同时可以实现空间数据的管理,保证数据的拓扑准确性和剖分三角网格的正确性。

2.1 边界点的输入
剖分算法输入的为点集的坐标值,非实际的拓扑数据,而且边界点组织直接影响到剖分输出的结果,需要在输入前进行必要的预处理。

如图4所示,按照拓扑关系,将输入数据划分为Poly gon 、Poly line 和Point 3种矢量数据,其涵盖了流体动力学问题中的海岸边界、岛屿、开边界和水深等高线等所有输入
数据。

其中,海岸边界采用Do ug las -Peucker 算法做简化处理后,根据实际计算需要,对边界和岛屿边界
页
25第第5期
赵晓东等:基于G IS 的高质量约束Delaunay 三角网格剖分
线上节点进行光滑处理。

图4GIS下Delaunay三角化网格剖分流程Fig.4Flow diagram illustrating Delaunary triangulation
mesh generation in G IS 2.2Delaunay三角网格剖分
采用本文改进的分治算法实现约束Delaunay 三角网格剖分,并将算法集成在A rcGIS平台中。

如图5所示,直接输入经过预处理的边界控制线Poly-line、边界控制区域Polyg on和控制点Point,可直接输出Delaunay三角网格剖分网格及节点。

2.3三角剖分网格拓扑检查及编辑
如果初始输入的约束之间没有形成尖锐角度,算法能够生成最小角度30b的Delaunay剖分网格,并具有较好的梯度。

但在初始约束不当和特殊约束条件下,算法还不能保证满足所有约束条件,因此对算法生成的剖分三角形的拓扑检查是必要的。

如图6所示,对局部没有满足约束条件的三角网格,可通过GIS下的拓扑编辑修正节点以满足其约束条件。

2.4剖分三角节点及网格的重编号
算法生成的Delaunay三角网格是没有规律的,给实际数值模拟带来不便,如水动数值模型的边界条件是从开边界开始设定的,因此需要对剖分网格节点及网格重新编号。

这样不仅可以和数值模型开边界的编号完全耦合,还使得相邻网格之间和相邻节点间的编号跨度减小,从而减少计算中调用各网格信息的时间花费,提高计算效率。

重编号算法采用A FM法的思想,以开边界为起始边界,按照空间邻接关系向前推进,直到遍历所有剖分节点,
具体步
图5G IS下的Delaunay三角网格剖分
Fig.5Delaunay triangulation mesh generation in
GIS
图6三角剖分网格拓扑检查
Fig.6T opological ch eck on triangulation meshes
骤如下[13]:1)如图7所示,给定起始边界线(这里为
开边界线)、三角剖分Polyg on和剖分节点Point数
据,设定排序方法;2)根据起始边界线空间检索与
之相交的节点Point,检索结果按照设定方法排序;
3)根据检索到的节点Po int空间检索与之相交的三
角剖分Poly gon,以其中心点为基点对未排序三角剖
分Poly gon按照设定方法排序;4)根据检索到的三
角剖分Poly gon空间检索与之相交剖分节点Point,
判断是否有未排序节点,如果没有,算法结束,否则,
按照设定方法排序后返回第3步,直到算法收敛。

页
26
第地理与地理信息科学第26卷
图7 三角剖分网格重编号
Fig.7 Re -numbering of triangulation meshes
3 算例及应用
采用本文算法分别使用随机离散点集和实例计
算海域进行了测试和应用,测试硬件为Intel Core Duo CPU 1.8GH z,内存2GB,操作系统Window s XP SP3,算法时间不包括输入/输出时间。

表1为几种主要算法在随机离散点集下的测试结果,从时间复杂度看,改进的分治算法基本与原分治算法保持一致,在使用了/影子0三角形的统一数据结构下,时间上明显优于其他几种传统算法。

本文提出的算法在保持原分治算法时间复杂度的前提下,满足了约束条件以适应实际问题计算的要求,其结果令人满意。

表2为日本博多湾(图8)和中国渤海区域的海洋数值模拟中的应用算例,Delaunay 三角剖分能很好地满足特征约束,计算速度快,能满足水动数值模拟非结构网格的计算要求。

表1 随机离散点集算例
Table 1 Examples using random points
s
改进分治传统分治扫描线逐点插入
算例1(1万节点)0.310.42 0.78 1.07算例2(10万节点) 4.25 5.6010.7822.65算例3(100万节点)
55.50
71.50
147.15
515.53
表2 应用算例计算输入参数及输出结果
T able 2 Inp ut parameters and output results of examples in case study
输入节点数输出节点数三角单元数计算时间
(ms )备注
实例1
69192179935251250无单元约束实例269192272137078255有单元约束实例3
4400
26623
15712
123
有单元约束
图8 日本博多湾Delaunay 三角剖分
Fig.8 Delaun ay triangu lation of H akata Bay,Japan
4 结语
本文在GIS 平台下基于改进的分治算法实现了高质量约束Delaunay 三角网格剖分,改进的分治算法保持了原算法的时间复杂度,并利用GIS 空间数据处理功能完成了数据输入和输出剖分网格的准确拓扑检查,提高了非结构化网格的质量和准确度;同时针对实际问题的需要,基于推进阵面算法的思想,结合空间邻近拓扑关系实现了三角剖分节点和网格的重新编号算法。

算例测试和实例应用表明,本文提出的方法集成了几种综合算法的优点,算法收敛速度快,网格生成质量高,满足实际问题的约束条件,应用效果较好。

参考文献:
[1] 芮一康,王结臣.Delaun ay 三角形构网的分治扫描线算法[J ].
测绘学报,2007,36(3):358-362.
[2] 胡金星,马照亭,吴焕萍,等.基于网格划分的海量数据Delau -nay 三角剖分[J].测绘学报,2004,33(2):163-167.
[3] 武晓波,王世新,肖春生.Delau nay 三角网的生成算法研究[J ].
测绘学报,1999,28(1):28-35.
[4] 王盛玺,宋松和,邹正平.基于约束Delau nay 三角化的二维非结
构网格生成方法[J].计算物理,2009,26(3):335-348.[5] LEE D T ,SCH ACH TER B J.Tw o algorithms for constructin g
a Delaun ay triangulation [J ].Computer and Information Sc-i ences,1980,9(3):219-242.
[6] FORTUNE S.A sweepli ne algorithm for Voronoi Diagrams[J].Al -gorithmica,1987,2(2):153-174.
[7] L AW SON C L.Softw are for C1surface interp olation[A].RICE
J R.M athem atical S oftw are III [C ].New York:Academic Press,1977.161-194.
[8] DW YE R R A.A Faster divide -and -con qu er algorithm for con -stru cting Delaunay triangulations[J].Algorith mica,1987,2(2):137-151.
[9] S HEW CH UK J R.Delaun ay refinement alg orith ms for triangu -lar mesh generation[J].Computational Geometry,2002,22(1-3):
页
27第第5期
赵晓东等:基于G IS 的高质量约束Delaunay 三角网格剖分
21-74.
[10] SH EWCH UK J R.Triangle:Engineering a 2D quality mesh gener -ator and Delau nay triangulator[A].LIN M C,M ANOCH A D.Applied Compu tational Geometry:T ow ards Geom etric En g-i n eering[C].First ACM W orksh op on Applied Compu tational Geometry,Lecture Notes in Computer Science,Vol.1148.Ber -lin:Springer -Verlag,1996.203-222.
[11] RUPPERT J.A Delaun ay refin ement alg orith m for quality 2-
dim ens ion al mesh gen eration [J].Algorithms ,1995,18(3):548-585.
[12] CHEW L P.Guaranteed -quali ty mesh generation for curved sur -faces[A].Proceedin gs of th e Ninth Annual ACM S ymposium on Computational Geom etry[C].San Diego,CA:ACM Pr ess ,1993.274-280.
[13] 赵晓东,王海龙,沈永明.基于GIS 的二维非结构化剖分网格
优化[J ].地理与地理信息科学,2010,26(2):46-48.
High -Quality C onstrained Delaunay Triangulation Based on GIS
Z HA O Xiao-dong 1,Y AN Xiao-bao 1,SH EN Y ong-ming 2,W A NG L iang 2
(1.China -J ap an Resear ch Center f or Geo -env ironmental S cience ,Pioneer P ar k of A cademician,
D alian Univ er sity ,D alian 116622;2.State K ey L abor ator y o f Coastal and Of f shor e
Engineer ing,D alian Univer sity of T echnology ,D alian 116023,China)
Abstract:T he unstr uctured mesh g eneratio n is one o f the key technical issues in many fields such as mechanical computat ion and numerical simulation.Based on the analysis o f ex isting unstructured mesh generation algo rithms,the impro ved divide -and -con -quer alg or ithm suppor ted by GIS is pro posed to deal w ith co nstr ained Delaunay t riangulation.T his method makes use o f the G IS spatia l to po lo gical relations to handle t he pre -processing of input algo rithm dat a,implements Delaunay refinement using a trian -g le -based data structur e,and checks the output meshes w ith accurate to po lo gy and co nstr aints.Based o n the idea of adv ancing front alg or ithm and spatial to po log ical relatio ns betw een the mesh nodes and tr iangulat ion,the re -numbering algo rit hm is also pr oposed to facilitate the assignment of the o pen bo undar y conditions as w ell as improv ing the computatio nal efficiency.T he a p -plication to the practical simulation show s that the pro po sed method can simplify numerical model pre -pro cessing of unstr uc -tur ed mesh g eneratio n,preserve the advantages of sever al algo rit hms w ith the or iginal running time,and impro ve the quality o f co nstr ained Delaunay t riang ulatio ns.
Key words:mesh generat ion;G IS;constrained Delaunay tr iangulation
(上接第19页)
A Automatic Mosaic Method in Unmanned Aerial Vehicle Images Based on Feature Points
L U Heng,L I Yo ng -shu,HE Jing,CHEN Q iang,REN Z hi-ming
(GI S eng ineer ing Center of Southw est J iaotong Univ er s ity ,Chengd u 610031,China)
Abstract:T he methods w hich curr ent unmanned aeria l vehicle(U A V)phantom fast splicing may select and t heir feasibility wer e elaborat ed,t he steady SIFT alg or ithm was intr oduced into the auto matic mo saic of U A V imag es,and the time co nsumpt ion o f the algo rithm during var ious st ages was analy zed.Considering the feature of U A V to impro ve the algo rithm,w hich can reduce the r ang e for sear ching feature points thr ough estimating the ov erla p o f neighbor ing images,and the most superio r size o f Gaussian kernel adapting in the U AV images was o bta ined thro ug h the ex periment w hen detect ing the criterio n space ext reme po ints,it can over come the tr aditional SIF T alg o rithm c s flaw by using the fix ed size met ho d,no t only reduce the time consum p -t ion,but also o bta in t he character istic po ints as fa r as po ssible.F inally the pr ecise transformat ion matr ix was o bta ined to com -plete imag es mo saic by using the L M method.Ex per imenta l results show that this algo rit hm has go od compatibility to automat ic mosaic of U A V imag es,which can remain the alg or ithm c s robust ,meanwhile,the precision and the efficiency w ere increased.Key words:U AV imag es;automatic mo saic;o verlapping of imag es;size of Gaussian ker nel
页
28第地理与地理信息科学
第26卷。