三种TIN的生成方法

生长算法

生长算法有两种:递归生长和凸闭包收缩算法。通俗来讲就是前者是从点集中心向外扩展,后者是由外向内收敛。

以递归生长算法为例,其基本步骤是:

  1. 在数据点中任取一点 P1P_1 (一般是位于数据点的几何中心附近),并寻找距离此点最近的点 P2P_2 。两者相连形成初始基线 P1P2P_1P_2 ,利用 Delaunay 三角网剖分准则,在点集中寻找第三点 P3P_3 ,形成第一个三角形 P1P2P3\triangle P_1P_2P_3

  2. 以初始三角形的三条边为基线,利用 Delaunay 三角网剖分准则,分别寻找能与该三条基线形成 Delaunay 三角形的 P4,P5,P6P_4,P_5,P_6 点,构造新的三角形。

  3. 重复步骤2,直到所有数据点处理完毕。

    Delaunay 三角网剖分准则

    • 满足空外接圆特性
    • 满足最大最小角特性

三角网生长算法逻辑简单,易于计算机编程实现,其最大的问题就在于生成过程中搜寻数据点时间效率很低,因此这部分的选择的算法主要是对搜寻数据点进行的优化[1]^{[1]}。由于时间有限,只能列出一篇主要论文。

利用点角改进 Delaunay 三角网生长算法

  1. 点角定义

    Delaunay 三角网中,任一点的点角为该点所构成的三角形中该点的角度之和。如下图所示,在 Delaunay 三角网中,点 P4P_4 所在的三角形有 P4P6P5, P4P3P6, P4P2P3, P4P1P2, P4P5P1\triangle P_4P_6P_5,\ \triangle P_4P_3P_6,\ \triangle P_4P_2P_3,\ \triangle P_4P_1P_2,\ \triangle P_4P_5P_1 ,易知 P4P_4 的点角为:

P4=1+2+3+4+5=360°\begin{equation} \angle P_4=\angle 1+\angle 2+\angle 3+\angle 4+\angle 5=360° \end{equation}

显然,除了 Delaunay 三角网中凸包点(最外层多边形构成的点)外,其他任何点的点角都为 360°。为了保证统一且算法更易实现,规定凸包点的点角也为 360°,如下图中三角网点 P2P_2 的点角为 P2=6+7+8=360°\angle P_2= \angle6+\angle7+\angle8=360°。如此,在构造 Delaunay 三角形之前需先寻找离散点的凸包,计算边界点的点角。

图1

  1. 计算凸壳

    常规的三角网生长算法不需要计算凸壳,这里主要是为查找凸包点并计算每个凸包点的外角,将其赋值为该点的点角,并为确定初始基线边做准备。

    使用格雷厄姆方法来求点集的凸壳。

    • 设点集中纵坐标最小的点为 P1P_1

    • P1P_1 和点集中其他各点用有向线段连接,并计算这些线段与水平线的夹角;

    • 线段按夹角大小排序(夹角相同时,按点到 P1P_1 的距离排序),得到一序列 P1P2,,PnP_1,P_2,\cdots,P_n

    • 依次连接所有点,便得到一多边形。显然,P1P2,,PnP_1,P_2,\cdots,P_n 是凸壳边界上的点(见下图)。依据“凸多边形的各顶点必在该多边形的任意一条边同侧”的原则,删去边界序列中的非凸点,从而得到凸壳点集,并以顺时针连接各凸壳点形成有向凸多边形(见下图)。

      图2

  2. CCW 方位测试

    CCW 方位测试(简称 CCW 检测)主要用来查找符合条件的点构成 Delaunay 三角形,可以判断出点与方向线的相互关系,CCW 方位测试如图下所示。

    图3

CCW(c,a,b)=axay1 bxby1 cxcy1=axcxaycy bxcxbycy\begin{equation} CCW(c,a,b)= \begin{vmatrix} a_x&a_y&1 \\\ b_x&b_y&1 \\\ c_x&c_y&1 \end{vmatrix}= \begin{vmatrix} a_x-c_x&a_y-c_y \\\ b_x-c_x&b_y-c_y \end{vmatrix} \end{equation}

根据公式 (2)(2) 计算出来的结果的正负性 (>0, <0, =0)(>0,\ <0,\ =0),可判断出点 ccab\overrightarrow{ab} 的左侧、右侧或直线上。

该论文在构造 Delaunay 三角形时,只有 ab\overrightarrow{ab} 右侧的点才能与 ab\overrightarrow{ab} 构成三角形。

  1. 构建三角网

    • 确定初始基线边,如下图 b 所示,子块的凸壳形成后,以纵坐标最小的点为初始基线边的第一点,从图中易知 P1P11\overrightarrow {P_1 P_{11}} 便是初始基线边,且凸多边形每条有向边都是待扩展的基线边。

    • 构建首个三角形,从所有离散点中查找能与 P1P11\overrightarrow {P_1 P_{11}} 构三角形最优的点,形成首个三角形存入 Triangle 中。如下图 b,首个三角形为 P1P10P11\triangle P_1P_{10}P_{11},同时,计算该三角形每个点的点角,其中 P10P_{10} 的点角为 P10=1\angle P_{10}= \angle 1P1P_1 的点角为 P1=2+3\angle P1= \angle 2+\angle 3P11P_{11} 的点角 P11=4+5\angle P_{11}= \angle 4+\angle 5。首个三角形同时也形成了三条有向边,将 P1P10, P10P11\overrightarrow {P_{1} P_{10}},\ \overrightarrow {P_{10} P_{11}} 存入 Edge(待扩展基线边)中, P11P1\overrightarrow {P_{11} P_{1}}P1P11\overrightarrow {P_1 P_{11}} 则不再需要扩展,从 Edge 中删除 P1P11\overrightarrow {P_1 P_{11}} 的记录( P11P1\overrightarrow {P_{11} P_{1}} 直接不需存)。

    • 扩展三角网,按顺序从 Edge(待扩展基线边)中取出一条边,从点角不为360°的点集中查找最优点构建 Delaunay 三角形。每扩展一个三角形,需要从 Edge 中删除该有向边,同时计算三个角的点角,生成两条新的有向边。判断 Edge 中有无与新生成的边方向相反且点相同的边,若有,则从 Edge 中删除该边,否则存入新生成的有向边。当 Edge 中无元素时子三角网构建完成,此时点集中所有点的点角为360°。

      图4

个人看法

在另一个生长算法的论文[3]^{[3]}中提到了几个概念,与上述论文有一定相通之处。

排除封闭点

在寻找第三点构网的过程中会不断产生封闭点。以 P1P1 为起始点经过3次遍历后的结果如下图所示,此时点 P3P3 已被三角形所包围,称之为封闭点。封闭点将不再参与下一次遍历边集寻找第三点的过程,因此第 4 次遍历前可将点 P3P3 从备选点集中移除。通过动态地排除封闭点,可将第三点的备选范围不断缩小,以减少寻点时间。

图5

这里的封闭点,其实和点角的概念相近,只不过点角封闭点一个具体判断方法。

由于以前参加编程小组学习建TIN的时候,就是顺着相似的思路写了一个简陋的算法,看到点角就有种亲切感,所以生长算法篇幅就比较多。

逐点插入法

逐点插入法的基本步骤[1]^{[1]}

  • 定义一个包含所有数据点的初始多边形;
  • 在初始多边形中建立初始三角网,然后迭代以下步骤,直至所有数据点都被处理;
  • 插入一个数据点 PP ,在三角网中找出包含 PP 的三角形 tt ,把 PPtt 的三个顶点相连,生成三个新的三角形;
  • LOP 算法优化三角网。

可以看出,逐点插入算法的思路非常简单,先在包含所有数据点的一个多边形中建立初始三角网,然后将余下的点逐一插入,用 LOP 算法确保其成为 Delaunay 三角网。各种实现方法的差别在于其初始多边形的不同以及建立初始三角网的方法不同。

一种改进的 Delaunay 三角剖分快速实现算法

对以往提出的逐点插入算法的研究发现,影响逐点插入法时间效率的主要因素为点的插入次序。针对这一特点,对逐点插入法进行改进。

  1. 基本概念

    • 极大三角形

      在平面内,加入三个点形成一个三角形,使其能够包含点集中的所有点。这个极大三角形就作为三角剖分的辅助初始三角形。

    • 影响域

      对每个插入点,外接圆包含该点的所有三角形所构成的多边形区域称为该点的影响域。

  2. 数据结构

    算法的执行效率与数据结构之间有着很密切的关系。该论文中使用的数据结构简要列出如下:

    • 有序顶点类

      保存散乱点集和排序后的顶点序列。

    • 有向边类

      保存三角形的边以及与该边相邻的三角形信息。

    • 三角形类

    • 保存生成的各个三角形的顶点及边信息

  3. 算法的主要步骤

    • 将点集中的点在平面内以 XX 轴为主,YY 轴为辅进行升序排序;

    • 虚构一个极大三角形成为最初的三角形,使其包含点集中的所有点;

    • 插入第一个点,这个点肯定在初始三角形内,连接这个点和极大三角形的各个顶点,形成最初的三角网;

    • 依次插入排好序的点到已生成的三角网中,对插入点找到该点的影响域,连接该插入点和它影响域的各个顶点,即重构了 Delaunay 三角网;

    • 删除所有以极大三角形顶点为顶点的三角形,即多余三角形,最后得到的三角网就是所求的 Delaunay 三角网。

  4. 新插入一个顶点的过程

    • 在已存在的 Delaunay 内新插入一个点,首先定位该点所在的三角形,由于对点集中的所有点进行了有序化,所有定位三角形时不像传统方法那样需要遍历当前所有三角形;
    • 其次找出有哪些三角形的外接圆包含该插入点,即确定出该点的影响域,删除影响域内的三角形,形成一个空腔,连接新插入点与空腔的各个顶点,形成多个三角形;
    • 最后把新生成的三角形放回三角网中。

个人看法

  1. 在上述论文的描述中,虽然并没有明确地使用 LOP 算法。但是个人认为实际上在新插入一个顶点的过程中,“寻找外接圆包含插入点的三角形,并删去这些三角形,再连接新插入点”的操作,就等效于使用了 LOP 算法。
  2. 无论是在生长算法,还是逐点插入法中,排序都是很重要的一种优化思路,如论文 [3,4] 都提到了排序的必要性。并且能连接成三角形的点基本上距离很近。结合无人机摄影中学到的 KD-树的概念,我们是否能采用类似的树结构进行优化呢?在论文 [5] 中也确实有四叉树划分点云数据文件进行并行计算动态调度的思路。

分割-合并算法

将分治算法思想应用于生成 Delaunay 三角网,递归地分割点集,直至子集中只包含三个点而形成三角形,然后自下而上地逐级合并生成最终的三角网。其基本步骤如下[1]

  • 把点集 VV 以横坐标为主,纵坐标为辅按升序排序,然后递归地执行以下步骤;
  • 把点集 VV 分为近似相等的两个子集 VLV_LVRV_R
  • VLV_LVRV_R 中生成三角网;
  • LOP 算法优化所生成的三角网,使之成为 Delaunay 三角网;
  • 找出连接 VLV_LVRV_R 中两个凸壳的底线和顶线;
  • 由底线至顶线合并 VLV_LVRV_R 中两个三角网。

以上步骤显示,这种算法的基本思路是使问题简化,把点集划分到足够小,使其易于生成三角网,然后把子集中的三角网合并生成最终的三角网,用 LOP 算法保证其成为 Delaunay 三角网。

Delaunay 三角网的生成算法研究

  1. 其基本步骤为:

    • 把点集 VV 以横坐标为主,纵坐标为辅按升序排序,然后递归地执行以下步骤;

    • VV 中数据量大于一给定值,把 VV 分为近似相等的两个子集 VLV_LVRV_R

    • VLV_LVRV_R 中用合成算法生成三角网;

    • 找出连接 VLV_LVRV_R 中两个凸壳的底线和顶线;

    • 由底线至顶线合并 VLV_LVRV_R 中两个三角网

    • 否则生成基本三角网

  2. 基本三角网生成过程

    • 生成凸壳过程

      使用 Larkin 改进和完善了的凸壳生成算法。

      把具有 x,y,x+y,xyx,y,x+y,x-y 的极大值和极小值的点相连作为初始凸壳,然后用一个递归过程 convex(I,J)convex(I,J) 把位于其中相邻两点间凸壳上的其余点找出来。为提高搜索效率,先把数据分块。 convex(I,J)convex(I,J) 搜索包含线段 IJ\overrightarrow{IJ} (从 IIJJ 逆时针方向)的矩形内的数据块。找出位于 IJ\overrightarrow{IJ} 右侧与 IJ\overrightarrow{IJ} 距离最大的点。该点就是要找的凸壳上的点。如果最大距离为零,则要判断它是否位于 I,JI,J 之间。如是,它也是凸壳上的点。

    • 初始三角网生成过程

    • 点插入过程

    • 局部优化过程(LOP

      局部优化过程(LOP)是所有生成 Delaunay 三角网的算法都要用到的关键过程。理论上说,不论用何种方法生成的三角网,只要用 LOP 过程进行处理,就能使它变为 Delaunay 三角网。这样一个重要的过程其实非常简单,它就是运用 Delaunay 三角网的性质,对由两个有公共边的三角形组成的四边形进行判断。如果其中一个三角形的外接圆包含第四个顶点,则将这个四边形的对角线交换。

      图6

  3. 寻找连接两个凸壳的底线过程

    分治算法分两步进行,首先找出连接左右两个凸外界的顶线和底线,然后由底至顶合并两个子三角网。这两个过程都要用到函数 pred(vi,vj)pred(v_i,v_j)succ(vi,vj)succ(v_i,v_j) 。函数 pred(vi,vj)pred(v_i,v_j) 找出在以 viv_i 为共同端点的一簇线段中,与 vivjv_iv_j 顺时针方向相邻线段的另一端点,succ(vi,vj)succ(v_i,v_j) 找出与 vivjv_iv_j 逆时针方向相邻线段的另一端点。

    通过这两个函数,可以找到连接 HL,HRH_L,H_R 两个凸壳的底线和顶线。

  4. 子三角网归并过程

    子三角网的合并过程也是一个迭代过程,它从连接 HL,HRH_L,H_R 的底线开始,在两个子三角网 TL,TRT_L,T_R 中寻找与底线组成 Delaunay 三角形的第三点 L1,R1L_1,R_1 ,选其中外接圆半径小的一个插入到最最终的三角网中。新生成的连接左右两个子三角网的边成为新的底线,逐步上推直到顶线结束,如下图所示。

    图7

个人看法

逐点插入算法实现比较简单,占用内存较小,但它的时间复杂度高;分治算法时间复杂度方面最好,但由于递归执行,内存占用大。为了均衡各种算法的不足,有不少新算法都是围绕着将多种算法融合的思路提出的,如论文 [6,7],以及上面提到的论文 [1] 。


视频推荐

一个关于Graham算法的参考视频

一个逐点插入法的参考视频(和上面那个视频是同一个作者!)。讲得还不错,可以去看看!里面还提到了 CCW 方位测试

参考文献

[1]武晓波, 王世新, 肖春生. Delaunay三角网的生成算法研究[J]. 测绘学报, 1999, 28(1):28-35.

[2]李建平,徐猛.利用点角改进Delaunay三角网生长算法[J].地理空间信息,2018,16(02):82-84+12.

[3]韦云舰,卢中.Delaunay三角网生长法的优化与应用[J].地理空间信息,2021,19(05):101-103+119+6.

[4]宋晓宇,王守金,王永会. 一种改进的Delaunay三角剖分快速实现算法[C]//.2008’中国信息技术与应用学术论坛论文集(二).,2008:129-131.

[5]LI Jian, LI Deren, SHAO Zhenfeng. A Streaming Data Delaunay Triangulation Algorithm Based on Parallel Computing[J]. Geomatics and Information Science of Wuhan University, 2013, 38(7): 794-798.

[6]芮一康,王结臣.Delaunay三角形构网的分治扫描线算法[J].测绘学报,2007(03):358-362.

[7]刘云,夏兴东,黄北生.基于分治算法与逐点插入法的Delaunay三角网建立算法的改进[J].现代测绘,2010,33(04):14-16.