不规则三角网(TIN)生成的算法 - 图文

第五章 不规则三角网（TIN）生成的算法

在第四章，基于三角网和格网的建模方法使用较多，被认为是两种基本的建模方法。三角网被视为最基本的一种网络，它既可适应规则分布数据，也可适应不规则分布数据，即可通过对三角网的内插生成规则格网网络，也可根据三角网直接建立连续或光滑表面模型。在第四章中同时也介绍了Delaunay三角网的基本概念及其产生原理，并将三角网构网算法归纳为两大类：即静态三角网和动态三角网。由于增量式动态构网方法在形成Delaunay三角网的同时具有很高的计算效率而被普遍采用。本章主要介绍静态方法中典型的三角网生长算法和动态方法中的数据点逐点插入算法；同时，还将给出考虑地形特征线和其他约束线段的插入算法。而其他非Delaunay三角网算法如辐射扫描法Radial Sweep Algorigthm（Mirante & Weingarten, 1982）等本文将不再介绍。

5.1三角网生长法

5.1.1递归生长法

递归生长算法的基本过程为如图5.1.1所示：

2

1

3

2

1

3

（a）形成第一个三角形 (b) 扩展生成第二个和第三个三角形

图5.1.1 递归生长法构建Delaunay三角网

（1）在所有数据中取任意一点1（一般从几何中心附近开始），查找

1

距离此点最近的点2，相连后作为初始基线1-2；

（2）在初始基线右边应用Delaunay法则搜寻第三点3，形成第一个

Delaunay三角形；

（3）并以此三角形的两条新边（2-3，3-1）作为新的初始基线； （4）重复步骤（2）和（3）直至所有数据点处理完毕。

该算法主要的工作是在大量数据点中搜寻给定基线符合要求的邻域点。一种比较简单的搜索方法是通过计算三角形外接圆的圆心和半径来完成对邻域点的搜索。为减少搜索时间，还可以预先将数据按X或Y坐标分块并进行排序。使用外接圆的搜索方法限定了基线的待选邻域点，因而降低了用于搜寻Delaunay三角网的计算时间。如果引入约束线段，则在确定第三点时还要判断形成的三角形边是否与约束线段交叉。

5.1.2凸闭包收缩法

与递归生长法相反，凸闭包搜索法的基本思想是首先找到包含数据区域的最小凸多边形，并从该多边形开始从外向里逐层形成三角形网络。平面点凸闭包的定义是包含这些平面点的最小凸多边形。在凸闭包中，连接任意两点的线段必须完全位于多边形内。凸闭包是数据点的自然极限边界，相当于包围数据点的最短路径。显然，凸闭包是数据集标准Delaunay三角网的一部分。计算凸闭包算法步骤包括：

（1）搜寻分别对应x-y，x+y最大值及x-y，x+y最小值的各二个点。

这些点为凸闭包的顶点，且总是位于数据集的四个角上，如图5.1.2（a）中的点7，9，12，6所示；

（2）将这些点以逆时针方向存储于循环链表中；

（3）对链表中的点I及其后续点J搜索线段IJ及其右边的所有点，计

算对IJ有最大偏移量的点K作为IJ之间新的凸闭包顶点，如点11对边7-9。

（4）重复（2）-（3）直至找不到新的顶点为止。

（a）初始边界7，9，12，6；（b）搜索凸闭包顶点11，5，4；（c）凸闭包

2

图5.1.2 凸闭包的计算(引自 Tsai，1993)

一旦提取出数据区域的凸闭包，就可以从其中的一条边开始逐层构建三角网，具体算法如下：

(a)第一个三角形 (b)第一层三角形

图5.1.3 凸闭包收缩法形成三角网

(1) 将凸多边形按逆时针顺序存入链表结构，左下角点附近的顶点排

第一；

(2) 选择第一个点作为起点，与其相邻点的连线作为第一条基边，如

图5.1.3（a）中的9-5；

(3) 从数据点中寻找与基边左最邻近的点8作为三角形的顶点。这样

便形成了第一个Delaunay三角形；

(4) 将起点9与顶点8的连线换作基边，重复(3) 即可形成第二个三

角形；

(5) 重复第(4)步，直到三角形的顶点为另一个边界点11。这样，借助

于一个起点9 便形成了一层Delaunay三角形；

(6) 适当修改边界点序列，依次选取前一层三角网的顶点作为新起点，

重复前面的处理，便可建立起连续的一层一层的三角网。 该方法同样可以考虑约束线段。但随着数据点分布密度的不同，实际情况往往比较复杂。比如边界收缩后一个完整的区域可能会分解成若干个相互独立的子区域。当数据量较大时如何提高顶点选择的效率是该方法的关键。

5.2数据逐点插入法

5.1节介绍的三角网生长算法最大的问题是计算的时间复杂性，由于每个三角形的形成都涉及所有待处理的点，且难于通过简单的分块或排序予

3

以彻底解决。数据点越多，问题越突出。本节将要介绍的数据逐点插入法在很大程度上克服了数据选择问题。其具体算法如下（见图5.2.1）：

（1） 首先提取整个数据区域的最小外界矩形范围，并以此作为最简

单的凸闭包。

（2） 按一定规则将数据区域的矩形范围进行格网划分，为了取得比

较理想的综合效率，可以限定每个格网单元平均拥有的数据点数。

（3） 根据数据点的（x,y）坐标建立分块索引的线性链表。 （4） 剖分数据区域的凸闭包形成两个超三角形，所有的数据点都一

定在这两个三角形范围内。

（5） 按照（3）建立的数据链表顺序往（4）的三角形中插入数据点。

首先找到包含数据点的三角形，进而连接该点与三角形的三个顶点，简单剖分该三角形为三个新的三角形。

（6） 根据Delaunay三角形的空圆特性，分别调整新生成的三个三

角形及其相邻的三角形。对相邻的三角形两两进行检测，如果其中一个三角形的外接圆中包含有另一个三角形除公共顶点外的第三个顶点，则交换公共边。

（7） 重复（5）-（6）直至所有的数据点都被插入到三角网中。

（a）第一分块数据插入后 (b) 第二分块数据插入后 (c)全部三角形

图5.2.1 逐点插入法构建Delaunay三角网

可见，由于步骤（3）的处理，保证相邻的数据点渐次插入，并通过搜寻加入点的影响三角网（Influence Triangulation），现存的三角网在局部范

4

围内得到了动态更新。从而大大提高了寻找包含数据点的三角形的效率。

5.3带约束条件的Delaunay三角网

当不相交的地形特征线、特殊的范围边界线等被作为预先定义的限制条件作用于TIN的生成当中时，必须考虑带约束条件的Delaunay三角网。最简单的处理方法是所谓的“加密法”，即通过加密约束线段上的数据点，将约束数据转换为普通数据，从而按标准Delaunay三角形剖分即可。尽管该方法加大了数据量并改变了原始数据集，但由于简单易行、稳定可靠，在许多情况下可以很好地满足需要。该方法唯一的问题在于如何恰当地确定特征线上加密数据点之间的距离，一般取平均数据点间距的一半或更小即可。以下内容主要介绍直接处理约束线段的算法。

5.3.1带约束条件的Delaunay三角网的定义

定义 1：给定一个d维欧基里德空间E和一个N点 mi 集 M。那么，关联的 Voronoi图 (又称 Thiessen多边形 )为覆盖E的一个凸多边形序列 (V(m1 )，V(m2 )，…，V(m N))，其中，V(mi)包括E中所有以M中的 mi 为最近点的点，即 V(mi)=p∈E∶ Vj，1≤j≤N，d(p，mi)≤ d(p，mj)，d表示欧基里德距离。Voronoi图的几何对偶 (dual)，即把点 mi 联结起来而得到的邻接格网称为M的Delaunay三角网。显然，Delaunay三角网的元素之并等于M的凸包之内部。Delaunay三角网自然推广到输入数据不仅包括点集 M，还包括不相交叉的直线段集L。在计算几何里，这类问题称约束Delaunay三角网（Constrained Delaunay Triangles，简称 CDT）问题。对地形数据来说，L即地形特征线段集（朱庆，陈楚江，1998）。

定义 2：令单点集M和线段端点集E之并为V（V=M∪E），如果在V的每个Delaunay三角形的外接圆范围内不包含任何与三角形的顶点均通视的其它点，而点Pi与Pj（Pi，Pj∈V）当且仅当连线PiPj 不与L中的任何约束线段相交叉（除在端点处外）时才互相通视，那么称这个 Delaunay三角网为V由L约束的 Delaunay三角网（朱庆，陈楚江，1998）。

5.3.2带约束条件的Delaunay法则

带约束条件的三角网仍然满足Delaunay法则，但其局部等角特性有较小的改变。当需要考虑约束条件时，可视图有助于重新定义Delaunay法则和Lawson LOP交换原则。对数据点及作为约束条件的断裂线，可视图由互相可视的任意两点连接而成。在可视图中，除在断裂线的端点处外，连

5

接线与任一断裂线都不相交（图5.3.1）。由此Delaunay法则及Lawson LOP交换可以重新定义为：

带约束条件的Delaunay法则：只有当三角形外接圆内不包含任何其它点，且其三个顶点相互通视时，此三角形才是一个带约束条件的Delaunay三角形。

带约束条件的DelaunayLawson LOP交换：只有在带约束条件的Delaunay法则满足的条件下，由两相邻三角形组成的凸四边形的局部最佳对角线（Locally Optimal Diagonal）才被选取。

图5.3.1 9个点与两条约束线段的通视图(引自 Tsai，1993)

5.3.3顾及线段约束的三角网生成算法

考虑线段约束可以在形成Delaunay三角形的同时进行，如根据带约束条件的Delaunay法则建立静态三角网的生长算法就是如此。而采用更多的方法是在动态生成三角网的基础上，采用两步法实现CDT的建立。所谓两步法即分以下两步完成：

（1） 将所有数据包括约束线段上的数据点，建立标准的Delaunay

三角网。

（2） 嵌入线段约束，根据对角线交换法LOP调整每条线段影响区

域内的所有三角形。

在作为约束条件的地形特征信息存在时，当标准Delaunay三角网建立起来后，便可加入预先给定的约束线段以完成带约束条件的Delaunay三角网的构建。如图5.3.2所示，下面步骤用于完成约束线段的插入：

（1） 在三角网中插入一约束线段；

（2） 确定边界与约束线段相交的三角形，如果两个这样的三角形有

公共边，则将此公共边删除，最后形成约束线段的影响多边形； （3） 将影响多边形其它各顶点与约束线段的起始节点相连；

6

（4） 应用带约束条件的LOP交换，更新影响多边形内的三角网，使

约束边成为三角网中的一边；

（5） 重复步骤1-4，直至所有约束线段都加入三角网中。

(a)插入线段ab，搜索其影响多边形；（b）连接节点a与影响多边形的所有顶点；（c）应用带约束条件的Lawson LOP交换对三角网进行优化；（d）带约束线段ab的三角网

图5.3.2 约束线段ab插入到已有Delaunay三角网的过程(引自Tsai，1993)

5.3.4从等高线生成三角网

等高线是一种特殊的特征线，等高线也可以作为约束线段。从等高线生成三角网一般有三种算法：等高线离散点直接生成TIN方法、将等高线作为特征线的方法、自动增加特征点及优化TIN的方法。等高线离散点直接生成TIN方法直接将等高线上的点离散化，然后采用上面所讲的从不规则点生成TIN的方法。但是由于这种算法只独立地考虑了数据中的每一个点，而并未考虑等高线数据的特殊结构，所以会导致很坏的结果，如出现三角形的三个顶点都位于同一条等高线上（即所谓的“平三角形”）或者三角形某一边穿过了等高线这样的情况（图5.3.3）。这些情形按TIN的特性都是不允许的。因此，在实际应用中，这种算法很少直接使用。通常将等高线作为特征线来构建三角网。

7

（a） 三角形与等高线相交；（b）三角形的三个顶点都位于同一条等高线上

图5.3.3 对等高线进行不合理三角化的例子

将等高线作为特征线生成三角网一般有两种算法：将等高线作为特征线的方法、自动增加特征点及优化TIN的方法。

将每一条等高线当作断裂线或结构线时，对三角形而言，至多只能从同一等高线取两个点。图5.3.4显示了一个考虑等高线特性的Delaunay三角网。

图5.3.4 将等高线当作断裂线以建立三角网

自动增加特征点及优化TIN的方法是：仍将等高线离散化建立TIN，但采用增加特征点的方式来消除TIN中的“平三角形”，并使用优化TIN的方式来消除不合理的三角形比如三角形与等高线相交等，另外对TIN中的三角形进行处理以使得TIN更接近理想化的情况。使用手工方式增加特征点线，无论在效率方面，还是在完整性、合理性等方面都是很有限的。因此需要设计一定的算法来自动提取特征点。这些算法的原理大都基于原始等高线的拓扑关系。对TIN进行优化则需对三角形进行扫描判断并以一定的准则进行合理化的处理。

由等高线重建地形的方法中使用骨架线可以保留曲线段之间的拓扑关系。从等高线图生成的Voronoi图上提取骨架线，骨架线可用于提取附加点以消除“平三角形”。附加点的高程可由估算获得。基于该方法可以估计出合理的地形坡度，并且为TIN提取有意义的中间点（Gold，2000）。

从等高线图生成的Voronoi图上提取骨架线的原理如图（5.3.5）所示， 当Delaunay三角形的外接圆不包含Voronoi图的顶点时，Voronoi图顶点在骨架线上（图5.3.5（a））；当Delaunay三角形的外接圆包含Voronoi图的顶点时，Delaunay三角形的边就是边界线（图5.3.5（b））。

8

(a) (b) 图5.3.5 提取骨架线的原理（引自Gold，2000）

图5.3.6 提取骨架线后的等高线图（引自Gold，2000）

提取等高线图的骨架线后（图（5.3.6）），还要估计骨架线上点的高程，其原理如图（5.3.7）所示。假设Zc是有新增点的等高线高程，Zb是相邻等高线的高程，Zi是待估计骨架点的高程，RR是参考圆的半径，Ri是骨架点的半径，则高程Zi可由下式计算：

Zi?Zc?(Ri/RR)?(Zc?Zb)/2 （5.3.1）

9

图5.3.7 骨架点高程估计原理（引自Gold，2000）

图5.3.8 骨架线高程（引自Gold，2000）

图5.3.9（a）中的“平三角形”扭曲了实际地形，而使用增加了骨架点的等高线建立TIN并对TIN进行优化后，对地形表达的效果则要好得多

10

（图5.3.9（b））。

（a）山谷和山顶区域的平三角形

（b）地形地貌的实际表达

图5.3.9 从相同的等高线生成的不同的TIN模型（引自Gold，2000）

5.4基于栅格的三角网生成算法

前面几种方法都是由矢量方式来形成三角网，实际上使用栅格的方式也可建立三角网。在栅格方式下，数学形态学方法是比较好的选择之一（陈晓勇，1991；马飞，1996）。

5.4.1形态变换原理

数学形态学（Mathematic Morphology）是Matheron和Serra于1965

11

年创立的，主要用于研究数字影像形态结构特征与快速并行处理方法。通过对栅格数据形态结构的变换而实现数据的结构分析和特征提取。其中二值形态学（函数值域定义在0或1）是将图形视作集合，通过集合逻辑运算（交、并和补）与集合形态变换（平移、扩张和侵蚀），在结构元作用下转换到新的形态结构。 如果将要建立TIN的区域与一幅数字影像相对应，凡是与数据点对应的像素灰度值为1，其它的像素灰度值为0，则可对这个二值影像进行形态变换建立TIN。运用数学形态学建立TIN，一般要有如下步骤：

用形态学建立TIN，主要是为了确定相邻参考点间的拓扑关系，因而只与点之间的相对距离有关，而与点之间的实际距离无直接关系。因此，为了能快速处理，以参考点间的最小距离为像素大小，将内插区域转化为一幅二值影像，参考点所在的像素灰度值为1，其它像素的灰度值为0。

5.4.2泰森多边形的形成

设X为参考点像素集合，则除去这些参考点后的剩余部分（即X的余集X）的骨架（skeleton），即为建立TIN的泰森多边形。

定义：连续影像A的骨架SK（A）就是A的最大内切圆的圆心集合。所谓最大内切圆是指那些与A的边界至少在两点相切的圆。

利用条件序贯细化形态变换可求得骨架，且能保证A中各分量的拓扑邻接关系。其结果为连续单像元宽度以及各向同性的像元集合。具体算法如下：

c

设Ck为半径为k的栅格圆环，A为影像的一个子集，令

Ak??(A?Ci)?Ak?1?(A?Ck)

i?1k (5.4.1)

A0?A

选用结构元Li（i=1，2，…，8），则

SK(A)?A?{LK}；{Ak)

8 (5.4.2)

即A的骨架由A的条件序贯细化变换生成。迭代的终止条件为

i?1 ?(SK(A)?Li)?（空集合）c (5.4.3)

则以上骨架算法得到SK(X)，即所需要的泰森多边形。

12

5.4.3三角网的形成

若X为参考点集，Pi?X是X的任意一参考点，将与Pi所在的泰森多边形相邻的泰森多边形中的参考点与Pi相连接，就构成了以Pi为顶点的所有的三角形的边。其步骤为:

(i) 将Pi所在多边形扩张至边界（即y的骨架）

(5.4.4)

Di?Pi?{H}；SK(Xc)c

?111???H??111?

?111???则将Pi进行条件序贯扩张，直至充满该泰森多边形，同时不越过多边形的边界。

(ii) 提取与Pi所在的泰森多边形Di相邻的多边形集合。首先作H对

Di的扩张，跨越边界，然后将Di的元素去掉，剩下Di的边界与相邻多

边形的元素，再作条件序贯扩张，条件是不超越边界（即Xc的骨架）。Di相邻多边形的集合Di：

'Di'?[(Di?H)?Dic]?{H}；(SK(xc))c

' (5.4.5)

（iii）提取Di中属于X的点，即提取位于与Pi所在泰森多边形相邻的泰森多边形中的参考点集：

Qi?Di'?X

(5.4.6)

依次连接Pi与Qi中的点，生成TIN相应的边。

对X中的每一点作相同的处理，记录网点邻接以及有关信息并存贮，就构建了三角网数字地面模型TIN。图5.4.1所示为根据等高线数据建立的三角网。

13

图5.4.1 采用数学形态学方法建立Delaunay三角网

应当指出，将形态学的方法用于DEM研究还是近十几年来的事，并且由于它的抽象性和复杂性而不为许多人所知，但使用数学形态学建立TIN，可简化许多矢量方法所考虑的操作，而且如果进行并行处理，则建立TIN的速度会大大提高（陈晓勇，1991；马飞，1996）。另外，使用形态变换还可以很容易处理特征点和特征线约束问题，只需要在建立泰森多边形时加上这些点线即可。

参考文献

陈晓勇，1991，数学形态学与影像分析，测绘出版社，北京，共154页。 马飞，1996，数学形态学在遥感和地理信息系统数据分析与处理中的应用

研究，武汉测绘科技大学博士学位论文，共144页。

吴华意，1999，拟三角网数据结构及其算法研究，武汉测绘科技大学博士

学位论文，共142页。

朱庆，陈楚江，1998，不规则三角网的快速建立及其动态更新，武汉测绘

科技大学学报，23（3）：204-207。

Aumann, G., H. Ebner, and L. Tang, 1991. “Automatic derivation of skeleton

lines from digitized contours”, ISPRS Journal of Photogrammetry and Remote Sensing, Vol. 46:259-268.

Brassel, K. and Reif, D., 1979. Procedure to generate Thiesen polygon.

Geographical Analysis, 11(3): 289-303.

Thibault, D and Gold，C. 1999, Terrain receonstruction from contours by

skeleton retraction. Proceedings of the 2nd International Workshop on Dynamic and Multi-dimensional GIS, Oct. 4-6, 1999. Beijing. 23-27. Li Zhilin, 1990, Sampling Strategy and Accuracy Assessment for Digital Terran

Modelling, Ph.D. Thesis, The University of Glasgow, 298pp.

14

Mirante A. and Weingarten N., 1982, The Radial Sweep Algorithm for

Constructing Triangulated Irregular Networks, IEEE Computer Graphics and Applications, 2(3):11-21。

Petrie, G. and Kennie, T. (eds.), 1990. Terrain Modelling in Surveying and

Civil Engineering, Whittles Publishing, Caitness, England. 351pp. Thibault, D. and Gold, C.M., 1999, Terrain Reconstruction from Contours by

Skeleton Retraction. Proceedings of 2nd International Workshop on Dynamic and Multi-dimensional GIS, October 1999, Beijing, pp23-27。 Tsai, Victor J.D., 1993. Delaunay triangulations in TIN creation :an overview

and a linear-time algorithm, International Journal of GIS, 7(6): 501-524.

15

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