三角网数字地面模型快速构建算法研究_刘学军

第13卷　第2期2000年4月

文章编号:100127372(2000)0220031206

中　国　公　路　学　报

ChinaJournalofHighwayandTransport

Vol113　No12Apr.2000

三角网数字地面模型快速构建算法研究

刘学军,符锌砂,赵建三

(长沙交通学院路桥系,湖南长沙　410076)

摘　要:系统地研究了三角网数字地面模型TIN构建中的几个关键问题,提出了动态创建和维护三角形拓扑关系的方法,建立了快速确定点在三角形中的算法原理及空外接圆判断法则的简易表达式,所设计的逐点插入算法有着较高的执行效率,算法复杂度与点数几乎成线性关系。关键词:数字地面模型;三角不规则网;算法;拓扑关系中图分类号:U41216　　　文献标识码:A

AStudyofalgorithmforfastcontriangulation(LIUXue2jun,2san

(DepartmentofHighwayandBonsUniversity,Changsha410076,China)

Abstract:Thisproceduresthatholduptheefficiencyofconstructing

triangulatiTIN).AwayandalgorithmfordynamicallyrenewingthetopologicalrelationsinTpresented.ThepaperalsostudiesthewayoffindingtrianglethatcontainagivenpointinTINandtheformulaofemptycircum2circletest.TheincrementalinsertionalgorithmofconstructingTIN,hasatimecomplexitythatisaboutlineartothenumberofpointsset.

Keywords:digitalterrainmodel(DTM);topologicalrelation

triangulationirregularnet(TIN);algorithm;

三角网数字地面模型(TriangulationIrregularNet简写为TIN)作为地表(地貌和地物)的数字化表现手段和分析工具,以其几何结构良好、数据结构简单、地表重构精度高及对不规则区域和数据点分布密度适应能力强等特点,并在地学领域如公(铁)路勘测设计一体化、地理信息系统等得到了广泛的应用。其核心技术—散点的三角剖分算法,近十几年来一直是众多学者研究和关注的焦点。迄今为止出现了不少成熟算法,如以BOWYER和GREEN等人为代表的Voronoi图法[9,11]、SHAMES及HOEY的分割—合并算法[1,4]、LAWSON的逐点插入算

[9]法[2,3]、WATSON的空外接圆算法及GREEN和

[10,11]

等。其中较有代表SIBSON的三角网生成算法

性的是LAWSON的逐点插入算法和WATSON的空外接圆算法。

编程易LAWSON的逐点插入算法原理简单、实现。研究结果表明[6],逐点插入算法占用内存资源较少,但时间复杂度差。SHANES和HOEY已证明[4],对N个数据点建立Delaunay三角剖分网,至少需0(NlgN)的时间。笔者对逐点插入算法进行了详细的研究,分析了该算法中效率的制约因素,提出了动态建立和更新三角网拓扑关系的算法,建立了快速三角形定位方法原理及运行次数更少的空外接圆判断公式。基于本文原理所编制的逐点插入算

收稿日期:1999207206

基金项目:交通部“九五”科技攻关项目(95205201207)

作者简介:刘学军(19652),男,陕西合阳人,长沙交通学院副教授,工学博士研究生.

　32　　　　　　　　　　　　　　　中　国　公　路　学　报　　　　　　　　　　　　　　2000年法程序,有着较高的执行效率。借助于Windows的文件映射技术,能一次性处理上千万个点。试验结果表明,时间复杂度与数据量基本成线性关系。

1　逐点插入算法介绍及分析

逐点插入法基本思想是:在一已存在的三角网中插入一点,该点与其所在的三角形三顶点形成新的三个三角形,然后用对角线交换法来优化新形成的三角形,从而保证所建三角网为Delaunay三角网。

逐点插入算法的基本步骤为:

(1)建立包含所有数据点的初始多边形,该多边形可为一个正三角形或两个直角三角形组成的矩形。

(2)从数据域中取出一点P,做如下工作:①找出包括点P的三角形T,设T的三顶点为V1、V2、V3;②P与T三顶点相连,形成三个新的三角形T1、T2、T3;

③对所有新形成的三角形,化。

(3)LOP局部算法(LocalOptimizationProcedureLOP)保证了所生成的三角形为Delaunay三角形。它利用了Delaunay三角形的空

图2　LOP优化过程

形数与点数的关系为:三角形数=边界点数+2×内点数-2)。无论是建立格网索引还是全局搜索,基于点在多边形中判断的过程是一个很费时的过程。该问题也是TIN数模应用中经常碰到的问题。112　LOP局部优化过程(拓扑更新)

LOP局部优化过程是基于具有公共边的两个

三角形进行的,LOP时,要快速,虽说可以通过拓,但由于逐点插入算,因而如何在动态过程中创建和,直接影响到算法的执行效率。113　数值计算问题

对每一次LOP过程都需要进行空外接圆检测,该过程在整个算法过程中的执行频率极高,故空外接圆检测方法的好坏对程序效率影响很大。空外接圆检测过程是一个数值分析与计算过程,因而应在计算稳定可靠的前提下,尽量减少计算次数和较费时的函数计算,从而提高执行效率。

外接圆特性,即一个三角形为Delaunay三角形,该三角形外接圆中不含有其余任何数据点(图1)。方法为:对个有公共边的两个三角形组成的四边形进行判断,如果其中一个三角形的外接圆中包含有另一三角形除公顶点外的第三顶点,则交换公共边,图2表示了该过程

。

[1]

2　点在三角形中查找(定位问题)

在三角剖分和应用中,经常碰到的问题是定位一个点在哪个三角形中,一般做法是扫描整个或局部三角形网格,利用点在多边形中原理判断计算。当三角形数目较大时,这是一个很费时的过程,如果建立了三角形网络的拓扑关系,则可利用三角形面积坐标和拓扑关系来解决这一问题。211　三角形面积坐标

如图3所示,三角形△V1V2V3的三顶点坐标为V1(x1,y1)、V2(x2,y2)、V3(x3,y3),任给点为P(x,

图1　Delaunay三角形空外接圆特性

y),按有限元理论,P在三角形△V1V2V

3

内的面积

(1)

对上述算法的研究分析可知,影响算法执行效率主要有以下几个因素:

111　点在三角形中的查找(定位问题)

坐标L1、L2、L3定义为

L1=A1 A;L2=A2 A;L3=A3 A

式中:A为△V1V2V

3

面积;A1、A2、A

3

分别为

算法中每插入一点,首先要定位该点在哪个三

角形中,随着点数增加,三角形数量成倍增加(

三角

△V2V3P、△V3V1P、△V1V2P的面积。

面积坐标(L1、L2、L3)与直角坐标的关系为

第2期　　　　　　　　刘学军,等:三角网数字地面模型快速构建算法研究　　　　　　　　　　　　33

L1=[(x3-x2)(y3-y2)-(x-x2)(y-y2)] AL2=[(x1-x3)(y1-y3)-(x-x3)(y-y3)] AL3=[(x2-x1)(y2-y1)-(x-x1)(y-y1)

]A

L1=(x3-x2)(y3-y2)-(x-x2)(y-y2)

(2)

212　点在三角形中的判断

参看图3,并设在三角形顶点逆时针排列,则当P在三角形中时,必有其所有的面积坐标大于零,即

L1

>0;L2>0;L3>0　　而若P不在三角形中,则至少一个面积坐标小于零。如图4所示,P在三角形之外,且为V2V3的外侧,则有

L1<0;L2>0;L3>0

(3)　L2=(x1-x3)(y1-y3)-(x-x3)(y-y3)

L3=(x2-x1)(y2-y1)-(x-x1)(y-y1　　另外,该算法还取决于初始三角形的位置,这可通过建立三角形索引来解决。

3　动态三角形拓扑关系的创建与维护

事实上,无论是构网过程中,还是在以后的应用中,三角形拓扑关系都扮演着相当重要的角色,而在构网过程中,要解决的问题是如何动态地创建和更新维护三角形的拓扑关系。

以下为叙述方便,给出几个定义及约定。

定义一:三角形数据结构定义如下(仅为本文定义)

typedefstructtagT{longA[3:三角形顶点按逆时针排列,三角形边的

　图3　三角形面积坐标　　　　　图4　P　　事实上,向。,利用面积坐标的这一特性,。图5。

0、1、2,即V[0]V[1]为第0号边,V[1]V[2]为第1号边,V[2]V[3]为第2号边,并边号可循环,即当边号i=i-1≤0时,i=2,当i=i+1≥2时,i=0。

约定二:三角形表示方法,任意三角形i,其三顶点表示为i.V[j](j=0,1,2),相邻三角形为i.A[j](j=0,1,2)。

定义二:基三角形,当插入一点P时,包含P的三角形为基三角形,用bt表示。

定义三:分裂三角形,由P点与基三角形顶点

图5　点在三角形中查找

连接而成的三角形,用st表示。

定义四:拓扑三角形,任意三角形三条边的相邻三个三角形为该三角形的拓扑三角形,用tt表示。

以上三种三角形存在着如下的数量关系:当P点落在三角形中时,基三角形数=1,分裂三角形数=3,基三角形的拓扑三角形数=3;当P落在三角形某一边上时,基三角形数=2,分裂三角形数=4,基三角形的拓扑三角形数=3或4。图6表示了上述三种三角形。

约定三:对分裂三角形进行编号和组成时,由基三角形的第一顶点开始,按逆时针方向进行,并依次称为第一、第二、第三分裂三角形。

定义五:第四顶点,任何一分裂三角形与其拓扑三角形组成具有公共边的四边形,在拓扑三角形中除公共边顶点之外的另一顶点,用Q表示,当前插

设初始三角形号为△1,计算P在△1的面积坐标,其符号情况为(上标为三角形号):

L1<0;L2>0;L3>0

1

1

1

取小于零的面积坐标对应边(图中为23边)的邻面三角形△2,则有

L>0;L<0;L>0

21

22

23

计算P在△5中的面积坐标

L>0;L>0;L>0

51

52

53

即P的三个面积坐标值都大于零,故P在△5中。

实际工作中,仅关心的是面积坐标的正负情况,故当三角形顶点以逆时针方向排列时,式(2)的分母

A恒大于零,因而式(2)可简化为

　34　　　　　　　　　　　　　　　中　国　公　路　学　报　　　　　　　　　　　　　　2000年

当P落在某一三角形的边上时,基三角形为

bt1、bt2,增加两临时记录tmp1和tmp2来记录bt1和

则第一、第二分裂三角形编号为bt1、bt2信息(图8)。

第四分裂三角形编号为nt+1、bt2,第三、nt+2,依约定三,则可组成分裂三角形及其拓扑三角形如表2

所示。

图6　三种三角形定义

入点用P表示。

311　插点时拓扑组成与更新

设已存在的三角网中的三角形数为nt,则在当前三角网中插入一点P时,无论P是落在三角形中还是落在三角形边上(区域边界除外),三角形数目总是增加两个,因而分裂三角形及其拓扑组成有如下规律:

当P落在三角形中时,

设基三角形号为bt(图7),另增加一临时记录tmp,用以记录基三角形bt未修改前的信息。把第一分裂三角形st1编号为原基三角形号bt,第二、第三分裂三角形st2、st3编号为nt+1、nt+2,依约定三和图7,其拓扑三角形如表1所示。

图8　P表2V]

112=bt2st3=nt+1st4=nt+2

p1.V[0]tmp1.V[1]tmp2.V[0]tmp2.V[1]

V[1]tmp1.V[1]tmp1.V[2]tmp2.V[1]tmp2.V[2]

V[2]PPPP

A[0]tmp1.A[0]tmp1.A[1]tmp2.A[0]tmp2.A[1]

A[1]A[2]bt2nt+1nt+2bt1

nt+2bt1bt2nt+1

拓扑组成

原基三角形的拓扑三角形的拓扑信息修改同点在三角形中。

312　LOP过程中的拓扑更新

如图9所示,以第一分裂三角形为例说明该过程。按311中所述,st1、st2结构如表3所示,设立临时变量tmp1和tmp2用以记录优化前st1、tt1的信息。经过LOP优化后,st1、tt1的结构组成如表4所示。

表3　优化前st1与st2的结构

三角形编号

st1tt1

V[0]

图7　P点落在三角形中

表1　分裂三角形和拓扑三角形的组成

三角形编号

V[0]

st1=btst2=nt+1st3=nt+2

tmp.V[0]tmp.V[1]tmp.V[2]

V[1]tmp.V[1]tmp.V[2]tmp.V[0]

V[2]PPP

A[0]tmp.A[0]tmp.A[1]tmp.A[2]

A[1]A[2]nt+2nt+1nt+2bt

btnt+1

顶点组成拓扑组成

顶点组成

V[1]

V[2]P

A[0]tt1at1

拓扑组成

A[1]st2st1

A[2]st3at2

24

33

2

最后修改原基三角形的拓扑三角形的拓扑信息。由于此时原基三角形信息已改变,故有关原基三角形信息应由变量tmp中取得,修改方法为:

若(tmp.A[0]).A[j]=bt,则(tmp.A[0]).A[j]=bt.(j=0,1,2);

若(tmp.A[1]).A[j]=bt,则(tmp.A[1]).A[j]=nt+1.(j=0,1,2);

若(tmp.A[2]).A[j]=bt,则(tmp.A[2]).A[j]=nt+2.(j=0,1,2)。

表4　优化后st1与tt1的结构

三角形编号

st1tt1

V[0]

顶点组成

V[1]

V[2]P

A[0]at1st1

拓扑组成

A[1]st2st3

A[2]tt1at2

44

3

P

2

下面由表3和表4的演变过程来分析其拓扑变化规律。

在LOP优化前后,三角形顶点变化及拓扑变化

第2期　　　　　　　　刘学军,等:三角网数字地面模型快速构建算法研究　　　　　　　　　　　　35

处的拓扑三角形拓扑关系。依312节方法有:

如果(tm

p1.A[is-1]).A[j]=tt1,(j=0,1,2),则

(tmp1.A[is-1]).A[j]=st1.;

如果(tmp2.A[it-1]).A[j]=st1,(j=0,1,2),则

(tmp2.A[it-1]).A[j]=tt1.

图9　LOP拓扑更新过程

与公共边在两三角形中的位置序号密切相关。为此,首先求出公共边23在st1、tt1中的位置序号,设为is及it,本例中is=0,it=1;由于边的位置序号由0开始,故其值也可看作是顶点数组V或拓扑数组A的下标。

31211　三角形组成

若以公共边位置序号作为顶点数组V的下标,则交换后三角形顶点变化正好处于该下标处。对分裂三角形st1,是用第四顶点Q(4)值交换is位置的原有值,即st1.V[is]=Q。而拓扑三角形tt1则是用加入点P取代位置it的原有值,即tt1.V[it]=P。31212　拓扑变化

若把公共边位置序号kA的下标,k-1和k处,4　空外接圆检测公式(数值计算问题)

在Delaunay三角网中,每一个三角形都要经过空外接圆检测。在算法中,这一过程是恒定的,它具有累计性。当数据较大,它在整个程序执行中所占用的CPU时间不容忽视。目前常见的做法是计算分裂三角形的外接圆圆心及半径,然后利用第四顶点到圆心距离和外接圆的半径关系进行判定。这一过程中要多次执行三角函数、开方、除法、平方等运算,,则是比较低的。。

,设分裂三角形st为△V1V2V3,其V2V3P,第四顶点为P。点P与△V1V2V3外接圆关系为

分裂三角形0)

st1.A[is-]st1.A[2]=tt1

st1.A[is]=at1=tmp2.A[it-

1]

拓扑三角形tt1(it=1)

tt1.A[it-1]=st1

tt1.A[it]=tmp1.A[is-

1]=st3

图10　空外接圆的简化处理

31213　拓扑三角形的拓扑变化

>0

-=0Β′

<0

表5列出了LOP前后st1、tt1拓扑三角形的拓扑信息变化情况。

表5　st1与tt1拓扑三角形的拓扑信息变化

LOP前

A[0]

st1.A[0]=tt1st1.A[1]=st2st1.A[2]=st3tt1.A[0]=at1tt1.A[1]=st1tt1.A[2]=at2

at1tt2A[1]st1st3st2st2

A[2]at2st1tt3

A[0]tst1tt2(P在圆外)(P在圆上)(P在圆内)

LOP后

A[1]st3st3st2st2

A[2]at2st1tt3

考虑Β′=Β,由于Β′+Α=180°,Β′=180°-Α

故有　-tgΑ=tgΒ

则有空外接圆检测函数f

>1　(P在圆外)

f=-tgΑ tg=1　(P在圆上)

<1　(P在圆内)

…

tt1tt1

…

st3

…

tt1tt1

…

st3

…………

依表中的变化不难看出,拓扑三角形的拓扑信

息变化之处也是在公共边k的前一位置和本身处。由于本身处的两个三角形即为LOP三角形,其相应拓扑关系在LOP时已改变,故仅要改变的是k-1

由三角函数关系和图10不难看出f的坐标公式

(4)f=-(△x13△x12-△y13△y12)(△yp3△xp2-△xp3△yp2)

△xab=xa-xb;△yab=ya-yb　(a=1,2,3,p;b=1,2,3,p)

上式f的执行效率比计算圆半径、圆心坐标及距离的效率高得多。

　36　　　　　　　　　　　　　　　中　国　公　路　学　报　　　　　　　　　　　　　　2000年为保证算法的稳定性,应注意当P落在三角形边上的情况,即当(△yp3△xp2-△yp2△xp3)=0时可不交换对角线。

号95205201207)有着运行速度快、对数据量无限制、整体建立地表模型的特点。参考文献:

[1]　LEEDT,SCHACHTERBJ.Twoalgorithmsforcon2

structingadelaunaytriangulation[J].242.

[2]　LAWSONCL.SoftwareforC1SurfaceInterpolation

[M].MathematicalSoftware .J.Rice,Ed.NewYork:AcademicPress,1977.161—194.

[3]　SLOANSW.Afastalgorithmforconstructingdelau2

nay

triangulation

in

theplane[J].Advanced

EngineeringSoftware,1987,9(1):34—55.

[4]　SHAMOSMI,HONYD.Closest2pointproblems[J].

Proceedingsofthe16thSymposiumontheons,1975:151—162.[,,.三角网的生成算法

Int.J.of

ComputerandInformationScience,1980,9(3):219—

5　算法分析和应用

现对上述算法原理做一简要的分析。

在LOP过程中,其主要工作是交换对角线及空外接圆检测,而在数值计算时,无论采用何种检测方法,对每一三角形而言,其运算工作量相同,因而这一过程是一常量。同时在维护拓扑关系时,动态维护算法无须进行全局或局部的三角形比较,故该过程

是一线性执行过程,其时间复杂度为0(N)。

每插入一点,其所在三角形的定位工作是比较费时的,理论上可以证明,笔者所给查找方法其平均

4

),而最坏情况为0(N2),因时间复杂度约为0(N1 而对所有点,其复杂度为0(N

5 4

),最坏为0(N2)。设

若对数据点或三角形建立了索引关系,则其时间复杂度还可减少,故可认为该部分的变化亦为常量变化。当数据量比较大时,的,即具有0(N)的复杂度。

的C++程序量,因而采用了Ws的文件映射技术调度和管理数据,即所有数据不进入内存,故比起数据在内存时,其运算要稍慢些,但依表6还是不难看出:基于上述原理的程序,其执行效率还是较高的。

表6　程序执行效率

点　　数

26030515651614771896966

J].,1999,28(1):28—35.

],,周建华.CAD软件设计[M].北京:北

京航天航空大学出版社,1996.174—183.

[7]　扬　钦,徐永安,陈其明,谭建荣.任意平面域上离散点

的三角化方法[J].软件学报,1998,9(1):241—245.

[8]　胡恩球,等.有限元网格生成方法发展综述[J].计算机

辅助设计与图形学学报,1997,9(4):378—383.

[9]　putingthen2dimensiondelaunay

tessellationwithapplicationtovoronoipolytopes[J].ComputerJournal,1981,24(2):167—172.

[10]　GREENPJ,putingdirichlettesse2

lationsintheplane[J].ComputerJournal,1978,21(2):168—173.

[11]　putingdirichlettesselations[J].Co2

mputerJournal,1981,24(2):162—166.

[12]　TSUNGPAOFANG,LESAPIEGL.Delaunaytri2

angulationusingauniformgrid[J].IEEE,ComputerGraphicsandApplication,1993,(5):34—47.

[13]　ANDREWJH,PETERLL.Onconformingdelau2

naymeshgeneration[J].AdvancesinEngineeringSoftware,1992,(14):129—135.

[14]　LARRYL.Schumaker.TriangulationinCAGD[J].

IEEE,ComputerGraphicsandApplication,1993,(1):47—52.

[15]　刘学军,符锌砂.TIN点单位算法及网形优化[J].中

三角形数

37005700582629542662893

CPU时间 s

12301261600

注:机型为奔腾MMX166,32M内存。

6　结　语

笔者系统地研究了逐点插入算法中的制约因素,提出了动态维护三角网拓扑关系的原理与算法。而正是在该算法维护下的三角网拓扑关系,得以保证了插点算法的线性执行效率。同时出于执行效率的要求,笔者还对在三角形拓扑关系下的三角形定位问题和数值计算问题做了探索研究,得出了有使用价值的结论。基于本文算法而研制的公路数字地

面系统HDTM(交通部九五科技攻关课题,合同编

国公路学报,1997,10(2):24—31.

[16]　符锌砂.基于航测数模的公路测设一体化系统[J].中

国公路学报,1997,10(2):31—38.

[17]　符锌砂.公路计算机辅助设计[M].北京:人民交通出

版社,1998.55—67.

本文来源：https://www.bwwdw.com/article/uak4.html