基于GMaps的Delaunay三角网格剖分——逐点插入法

A.逐点插入法基本原理

逐点插入法的优势不仅在于其理论上比较简单，容易理解，更重要的是该方法对于点的插入是动态的、局部的，正是其局部性，所以是高效的，特别适合于可能动态插入钻孔数据的地质建模方法。

其基本步骤为：①定义一个包含所有钻孔孔口坐标数据点的超级三角形(或者矩形，由两个三角形构成)；②遍历待插入点，对于某一个待插入点P，找到其所在的三角单元T；③把P与T的3个顶点相连，生成3个新三角形；④用局部最优方法(LOP)优化三角网，使其符合Delaunay空圆准则。

B.几何查询与拓扑修改

点插入过程中，需要大量的几何查询和拓扑修改。这些 *** 作主要分为两类：①查询。从三角单元中提取信息，比如在结点、边和三角形之间检查邻接关系，提取三角形单元进行可视化显示；②修改。改变TIN中的拓扑和几何信息，如在局部优化时进行边交换 *** 作和基于点插入的细分 *** 作。

在进行点插入时，首先是几何上的查询，即判断待插入点P(x，y)是否位于某个三角形内。传统上，对点是否在三角形内的这一看似简单的问题，是通过严格的几何计算来确定。这里，充分发挥GMaps拓扑数据模型的优势，判断方法如下，遍历三角形的三条半边E1、E2、E3(图3.4)，若点P不在半边E1所在的半平面(half plane)H1，2左侧，则直接返回“false”，进入下一个三角形；否则继续。利用2-Orbit(即 )由E1得到E2，判断点P是否位于其所在的半平面H2，3左侧。只有点P全部在半平面H1，2、H2，3、H3，1的左侧，点P位于三角形内。即：

P∈H1，2∩H2，3∩H3，1 　(3.2)

图3.4 点与三角形的关系判断

若点P不在当前三角形内，可以通过1-Orbit( )获取相邻的三角形继续进行判断，直到找到点P所在的三角形，整个过程如图3.5所示。实践表明，在整个查找过程中，对于任何初始拓扑单元d，查找路径总是以近似直线(图3.5所示的虚箭头)的“速度”趋近于目标的三角单元。因此，只需要有限的若干步骤即可完成，大大提高了算法的执行效率，也体现GMaps的优势。

在找到P所在的三角形Ti后，用点P将Ti一劈为三(图3.6)。首先，新建六条半边：V1P、V2P、V3P和PV1、PV2、PV3；然后，需要修改三角单元V1V2V3的拓扑关系：将V1V2的后继半边修改为V2P，并设V2P的后继边为PV1，PV1的后继边为V1V2；这样，PV1V2构成一个新的拓扑三角单元。同理，对V2V3和V3V1做拓扑修改，得到另外两个新的拓扑三角单元：PV2V3、PV3V1。

经过点劈裂后形成三个新的三角单元，接下来需要对每个三角单元进行局部优化(图3.7)。如图3.7a所示，问题将转化为拓扑单元V1V2V3的外接圆是否包含相邻的拓扑单元V2V4V3的顶点V4。如图3.7b所示，判断点V4是否在V2V4V3的外接圆内，只需考虑α+β的大小。若α+β＞π，则V4在V2V4V3的外接圆内，需要执行边交换。根据四边形的性质，同时有α+β＜2π：

图3.5 初始三角网和待插入点P的查找过程

图3.6 点劈裂三角形

图3.7 三角化过程中进行边交换

数字地下空间与工程三维地质建模及应用研究

因此，只需判断sinαcosβ+cosαsinβ＜0即可。如图3.7b所示，计算过程采用如下方法：

sinα=crossProduct(d1，d2)

sinβ=crossProduct(d3，d4)

cosα=scalarProduct(d1，d2)

cosβ=scalarProduct(d3，d4)

if((sinαcosβ+cosαsinβ)＜0)

return ok

else

return false

其中，crossProduct表示向量叉积，scalarProduct表示向量点积。

经过上述计算过程，若需要边交换，则如图3.7c所示，修改V3V1的后继半边为V1V4，V1V4的后继半边为V4V3，V4V3的后继半边为V3V1。同理对V1V2也做相应的拓扑修改。整个过程始终维持拓扑单元d的拓扑关系的一致性，结果如图3.7d所示。

C.点的动态插入

之所以采用逐点插入法构建Delaunay三角网，是因为其具有“开放性”。即当有新的数据点加入时，只需在局部进行修改或者扩充，而不需要改变整个体系结构。能够便于空间数据的动态插入与局部更新。

新插入数据点P，依据插入点所在位置的不同分为两种情况：

(1)新插入点P在原有Delaunay三角网边界之内，此时只需要判断新点位于哪一个三角形内，然后按照上节所述方法通过几何查询和拓扑修改即可完成。

(2)新插入点P在原有Delaunay三角网边界之外，此时只需要将新数据点追加到离散点数据库的最后。从Delaunay三角网边界为凸包的角度考虑，先找出已有三角网的边界边，可以通过α2(d)=d来获取边界边所形成的凸包多边形。然后，找到与待插入点P最近的边界点Q，可以证明待插入点与该边界点形成一条Delaunay边。在边界点集中，以点Q为参照向两个方向扩展形成新的边界，直到最后所形成的边界满足凸包原则即可。最后，对新形成的三角形进行局部优化调整，可以保证最终形成的三角网为Delaunay三角化。

D.点的动态删除

点删除是点插入的逆过程，需要使点删除之后的三角网仍然满足Delaunay法则，且要动态更新点、线和三角形之间的拓扑关系。关于点删除算法国内外研究相对较少，Chew(1986)最早提出从D-TIN中删除点的算法，但比较复杂，没有得到实际应用；此后，朱庆等(1998)提出了类D-TIN中点删除算法，Devillers(1999)提出基于凸耳权的点删除算法——凸耳消元法(ear elimination，EE)，实现了从二维D-TIN中删除单一点的删除算法；Marc(1997)在规则三角剖分(regular triangulation，RT)中也涉及对点的删除算法，Mostafavi et al.(2003)改进了凸耳消元法。

所谓凸耳，是指在D-TIN中删除点P时，从点P的影响域中逆时针选取3个相邻的点组成一个三角形，若其相对于点P是凸的，且其外接圆内部包含其他点，则该三角形称为一个凸耳。Devillers提出的凸耳权值点删除算法的基本原理是：查找以被删除点P为顶点的所有三角形，组成任意多边形H={ q0，q1，…，qk-1，qk=q0}。按照逆时针顺序形成该多边形的凸耳权值队列。从队列中选择权值最小的凸耳ear={qi，qi+1，qi+2}，则qiqi+2为D-TIN中的一条边，凸耳权计算公式为(Devillers，1999)：

数字地下空间与工程三维地质建模及应用研究

由于凸耳定义严格，动态编辑过程中凸耳队列的更新维护较复杂，可以以特征三角形来取代凸耳。所谓特征三角形，是指沿多边形边界(按逆时针方向)连接任意相邻3点组成的三角形，其权值定义如下：若三角形的矢量面积为负，其权值为无穷大，否则，其权值用式(3.4)计算。

删除点的步骤如下：

第一步：查找被删除点P相关联的所有顶点和边，利用0-Orbit( )即可完成(图3.8a)。

第二步：形成带权值的特征三角形队列。根据第一步找到的多边形顶点，按逆时针依次计算每一个三角形的凸耳权值。

第三步：查找符合Delaunay特征的D-TIN边。从特征三角形队列中删除权值最小的△tri=qiqi+1qi+2{ }，则qiqi+2为D-TIN的一条边。修改此特征三角形的前后数据元素，并计算新的特征三角形权值。

第四步：对角线交换及D-TIN的局部更新。根据拓扑关系和点与边的对应关系在D-TIN中查找以qi、qi+1、qi+2、p为顶点的两个三角形。交换这两个三角形所组成的凸四边形的对角线，这样以被删除点为顶点的三角形个数少1。交换对角线后，生成1条新边和2个新三角形，并维护三角形、边之间的拓扑关系。如图3.8所示，特征三角形队列中权值最小的三角形为△q1q2q3，以p、q1、q2、q3为顶点的两个三角形为F、F2，交换这两个三角形形成的凸四边形的对角线为E0={q1q3}。

图3.8 对角线交换与局部更新

第五步：重复第三、第四步，直到多边形顶点个数为3，转入第六步。

第六步：点删除与局部更新。

经过以上多次对角线交换，以被删除点P为顶点的三角形只剩下3个，合并这3个三角形，删除2个三角形元素、3条边和数据点P，更新三类元素间的拓扑关系，即完成了Delaunay点的动态删除。

3.2.1.1 基本理论

B.Delaunay于1934年提出了Delaunay三角网格的概念，它是Voronoi图(简称V图)的几何对偶图，具有严格的数学定义和完备的理论基础。

图3.1 Voronoi图(虚线)及对应的Delaunay三角剖分(实线)

3.2.1.1.1 Voronoi图

假设V={v1，v2，…，vN}，N≥3是欧几里得平面上的一个点集，并且这些点不共线，四点不共圆。用d(vi，vj)表示点vi与vj间的欧几里得距离。

设x为平面上的点，则:

区域V(i)={x∈E2d(x，vi)≤d(x，vj)，j=1，2，…，N，j≠i}称为Voronoi多边形，也称为该点的邻域。点集中所有点的Voronoi多边形组成Voronoi图，如图3.1所示。

平面上的Voronoi图可以看做是点集V中的每个点作为生长核，以相同的速率向外扩张，直到彼此相遇为止而在平面上形成的图形。除最外层的点形成开放的区域外，其余每个点都形成一个凸多边形。

3.2.1.1.2 Delaunay三角剖分

Delaunay三角形网格为V图的几何对偶图。在二维平面中，点集中若无四点共圆，则该点集V图中每个顶点恰好是3个边的公共顶点，并且是3个Voronoi多边形的公共顶点上述3个Voronoi多边形所对应的点集中的点连成的三角形称为与该Voronoi顶点对应的Delaunay三角形，如图3.1所示。如果一个二维点集中有四点共圆的情况，此时，这些点对应的Voronoi多边形共用一个Voronoi顶点，这个公共的Voronoi顶点对应多于3个Voronoi多边形，也就是对应于点集中多于3个的点这些点可以连成多于一个的三角形。此时，可以任意将上述几个点形成的凸包划分为若干三角形，这些三角形也称为和这个Voronoi顶点对应的Delaunay三角形。

所有与Voronoi顶点对应的Delaunay三角形就构成了Delaunay三角剖分。当无退化情况(四点共圆)出现时，点集的Delaunay三角剖分是唯一的。

3.2.1.1.3 Delaunay三角剖分的特性

Delaunay三角剖分具有两个重要特性:

(1)最小角最大化特性:即要求三角形的最小内角尽量最大，具体地说是指在两个相邻的三角形构成凸四边形的对角线，在相互交换后，6个内角的最小角不再增大，并且使三角形尽量接近等边。

(2)空外接圆特性:即三角形的外接圆中不包含其他三角形的顶点(任意四点不能共圆)，该特性保证了最邻近的点构成三角形，使三角形的边长之和尽量最小。

3.2.1.2 常用算法

Delaunay三角剖分方法是目前最流行的通用的全自动网格生成方法之一。比较有效的Delaunay三角剖分算法有分治算法、逐点插入法和三角网生长法等(Tsai，1993)，其中逐点插入法由于其算法的简洁性且易于实现，因而获得广泛的应用。其主要思路是先构建一个包含点集或区域的初始网格，再依次向初始网格中插入点，最后形成Delaunay三角剖分。

采用逐点插入法建立Delaunay三角网的算法思想最初是由Lawson于1977年提出的(Lawson，1977)，Bowyer和Watson等先后对该算法进行了发展和完善(Bowyer，1981Watson，1981)。目前涌现出的大量逐点插入法中，主要为以Lawson算法代表的对角线交换算法和以Bowyer-Watson算法代表的空外接圆法。

3.2.1.2.1 Lawson算法

Lawson算法的主要思想是将要插入的数据点逐一插入到一个已存在的Delaunay三角网内，然后再用局部优化算法(Local Optimization Procedure，LOP)优化使其满足Delau-nay三角网的要求，其主要步骤如下:

图3.2 包含点集的三角形

第一步:构建一个三角形或多个三角形(通常情况下为一个三角形，该三角形也称超级三角形)，使点集中所有点都落在该三角形内，如图3.2所示。另外，也可以用一个矩形包围所有点，再将矩形对角线相连生成一对三角形作为初始网格。

第二步:将包含点集的三角形作为第一个三角形单元加入初始三角形网格中。

第三步:依次插入一个点P，并在三角网中找出包含P的三角形T，此时存在3种情况:

(1)P在三角形T内部，即P不落在三角形T的边或顶点上，把P与T的3个顶点相连，生成3个新的三角形，将新生成的三角形加入三角网，删除三角形T，如图3.3(a)所示。

(2)P在三角形T的一条边上，但不在顶点上，把P与T的3个顶点中与P相对的顶点相连，生成两个新的三角形，将新生成的三角形加入三角网，删除三角形T，如图3.4(b)所示。

(3)P在三角形T顶点上，不需处理，进行下一步 *** 作。实际上，若点集中不存在相互重合的点，则P落在三角形T顶点上这种情况不会发生。

图3.3 向三角形网格中插入点P

上一步处理中，直接将三角形的顶点与P相连生成新的三角形单元不一定满足Delau-nay三角剖分的要求，因此，需要采用LOP算法对局部三角形进行优化。该算法的主要 *** 作是进行边交换使新生成的三角形单元及其相邻的单元满足空外接圆特性。这样一个重要的过程其实非常简单，它就是运用Delaunay三角剖分的空外接圆特性对由两个有共用边的三角形组成的四边形进行判断，如果其中一个三角形的外接圆包含第四个顶点，则将这个四边形的对角线交换，如图3.4所示，在交换对角线后两个新三角形都满足空外接圆特性。

当两个有共用边的三角形具有相同的外接圆时，即存在四点共圆情况，此时，按照空外接圆特性，可以交换对角线也可以不交换对角线，那么这时候就要按照最小角最大化特性进行处理。具体做法是，先计算不进行边交换前三角形对中的6个内角的最小值，再计算进行边交换后6个内角的最小值，判断边交换前、后最小角是否变大，若是，交换，否则，不交换。上述 *** 作就是为了尽量满足最小角最大化特性。

图3.4 四点不共圆时的边交换

图3.5 点集的Delaunay三角网格

第四步:在所有点被插入之后，从网格中删除多余的三角形单元，最后得到的网格就是点集的Delaunay三角网格，如图3.5所示。

按照上述算法，编程环境为VC++，Lawson算法的代码如下，其中参数CSurf*surf表示网格平面，网格结点和单元分别存储在成员pNodes和pTrgls中。函数IsPointInTrian-gle(x，y，tx［0］，ty［0］，tx［1］，ty［1］，tx［2］，ty［2］)用于判断点(x，y)是否落在由三点(tx［0］，ty［0］)、(tx［1］，ty［1］)和(tx［2］，ty［2］)组成的三角形中函数LOP(CTrgl*ta，CTrgl*tb)采用LOP方法优化三角形对ta和tb。

三维地质建模方法及程序实现

三维地质建模方法及程序实现

三维地质建模方法及程序实现

3.2.1.2.2 Bowyer-Watson算法

Bowyer-Watson算法也是一种先形成初始网格，再逐步插入数据点进行细化的算法。与Lawson算法不同的是，Bowyer-Watson算法插入点时需要判断网格中三角形单元的外接圆是否包含该结点，而不是判断该三角形单元本身是否包含该结点另外一个区别是Bowyer-Watson算法在每次插入一个点之后不需要像Lawson算法那样用LOP算法优化三角网。

Bowyer-Watson算法的主要步骤如下:

第一步:建立一个包含整个点集的初始网格。通常情况下初始网格为一个三角形或由一个矩形连接对角线得到的一对三角形。

第二步:依次将数据点插入到最近更新后得到的网格中，形成新的网格并更新，分为三小步:

(1)插入一个新点 P 到现有的 Delaunay 三角网格中，如图 3.6(a)所示。

(2)寻找并删除所有外接圆包含 P 点的三角单元，形成一个 Delaunay 空腔。图3.6(a)中的阴影三角形为外接圆包含 P 点的三角单元图 3.6(b)所示为形成的一个Delaunay 空腔。

(3)连接 P 点和 Delaunay 空腔边界上所有各点，形成新的 Delaunay 三角网格，如图3.6(c)所示。

图3.6 向三角形网格中插入点

第三步:删除多余的三角形，即如果某个三角形中的某个结点不属于原始的点集，则删除该三角形最后结果即为点集的Delaunay三角剖分网格。

按照上述算法，编程环境为VC++，Bowyer-Watson算法的完整代码如下，其中参数CSurf*surf表示网格平面，网格结点和单元分别存储在成员pNodes和pTrgls中。函数TriangleInCircle()的作用是判断点(xp，yp)是否落在由三点(x1，y1)、(x2，y2)和(x3，y3)组成的三角形的外接圆内。

三维地质建模方法及程序实现

三维地质建模方法及程序实现

三维地质建模方法及程序实现

三维地质建模方法及程序实现

采用上述算法和代码，对一个包含46个点的点集进行Delaunay剖分，图3.7(a)所示为输入的点集，得到拥有74个三角形单元的网格，如图3.7(b)所示。

图3.7 Bowyer-Watson算法剖分实例

具体什么不太懂，不过我收藏了一个以前别人讲这方面的东西，你看看，希望对你有帮助。

Delaunay三角网是俄国数学家B.Delaunay于1934年发现的。关于Delaunay三角网构建的研究有许多，但由于本课题具有数据量大的特征，不宜直接沿用已有构建方法，笔者针对本课题数据特征，研究获得了适应本课题，速度较快的构建方法。Delaunay三角网有一个特性，每个三角网形成的外接圆都不包含其他参考点。利用这一个性质，我们可以直接构成Delaunay三角网：

一、建立第一个三角形

1、判断用来建立TIN的总脚点数，小于3则报错退出。

2、第一点的选择：

链表的第一节点，命名为Pt1

3、第二点的选择：

A.非Pt1点； B.距Pt1最近

命名为Pt2

4、第三点的选择

A.非Pt1，Pt2点；

B.与Pt1，Pt2点组成的三角形的外接圆内无其他节点；

C.与Pt1，Pt2组成的三角形中的角Pt1Pt3Pt2最大。

命名为Pt3

5、生成三边，加入边表

6、生成第一个三角形，组建三角形表

二、扩展TIN

1、从边表头取一边，要求：该边flag标志为假（只在一个三角形中）

2、从点链表中搜索一点，要求：

A、与边中的Pixel3在边的异侧；

B、该点与边组成的三角形的外接圆内无其他点

C、满足上面两条件的点中角Pt1Pt3Pt2最大的点为Pt3。

3、判断新生成的边，如果边表中没有，则加入边表尾，设定标志flag为假，如果有，则设定该边flag为真。

4、将生成的三角形假如三角形表

5、设定选中的边的标志flag为真

6、转至1，直至边表边的标志flag全部为真。

数据结构：

struct Pixel //脚点数据

{

double x,y,z,g

bool flag

}

struct List //数据链表

{

Pixel *pixel

List *next

}

struct Line //三角形边

{

Pixel *pixel1 //三角形边一端点

Pixel *pixel2 //三角形边另一端点

Pixel *pixel3 //三角形边所对顶点

bool flag

}

struct Linelist //三角形边表

{

Line *line

Linelist *next

}

struct Triangle //三角形表

{

Line *line1

Line *line2

Line *line3

Triangle *next

}

以下是我程序中关于建网的部分：

// DelaunayTIN.cpp: implementation of the CDelaunayTIN class.

//

//////////////////////////////////////////////////////////////////////

#include "stdafx.h"

#include "Reconstruction.h"

#include "DelaunayTIN.h"

#include "MyMath.h"

#include <math.h>

#include "ListControl.h"

#ifdef _DEBUG

#undef THIS_FILE

static char THIS_FILE[]=__FILE__

#define new DEBUG_NEW

#endif

//////////////////////////////////////////////////////////////////////

// Construction/Destruction

//////////////////////////////////////////////////////////////////////

CDelaunayTIN::CDelaunayTIN()

{

}

CDelaunayTIN::~CDelaunayTIN()

{

}

/////////////////////////////////////////////////////////////////////////////

//函数名： CreateDelaunayTIN

//编写者： Polaris

//参考资料：

//功能： 用给定的数据链表数据，组建Delaunay不规则三角网

//输入参数：数据链表list区域范围(XMin,YMin),(XMax,YMax)

//输出参数：不规则三角网首三角形地址

//备注：

/////////////////////////////////////////////////////////////////////////////

struct Triangle * CDelaunayTIN::CreateDelaunayTIN(List *list)

{

//组建第一个三角形

CMyMath MyMath

CListControl ListControl

struct List *node

struct Pixel *pt1,*pt2,*pt3

bool flag

struct Triangle *TIN

pt1=list->pixel

pt2=list->next->pixel

node=list->next->next

while(node!=NULL)

{

if(MyMath.Calculate2PtDistanceIn3D

(pt1->x,pt1->y,pt1->z,node->pixel->x,node->pixel->y,node->pixel->z)

<MyMath.Calculate2PtDistanceIn3D

(pt1->x,pt1->y,pt1->z,pt2->x,pt2->y,pt2->z))

{

pt2=node->pixel

}

node=node->next

}

node=list->next

pt3=NULL

while(node!=NULL)

{

if(node->pixel==pt1 || node->pixel==pt2)

{

node=node->next

continue

}

if(pt3==NULL)

{

pt3=node->pixel

}

else

{

if((pow(MyMath.Calculate2PtDistanceIn2D(pt1->x,pt1->y,node->pixel->x,node->pixel->y),2)+pow(MyMath.Calculate2PtDistanceIn2D(pt2->x,pt2->y,node->pixel->x,node->pixel->y),2)-pow(MyMath.Calculate2PtDistanceIn2D(pt1->x,pt1->y,pt2->x,pt2->y),2))/(2*MyMath.Calculate2PtDistanceIn2D(pt1->x,pt1->y,node->pixel->x,node->pixel->y)*MyMath.Calculate2PtDistanceIn2D(pt2->x,pt2->y,node->pixel->x,node->pixel->y))

<(pow(MyMath.Calculate2PtDistanceIn2D(pt1->x,pt1->y,pt3->x,pt3->y),2)+pow(MyMath.Calculate2PtDistanceIn2D(pt2->x,pt2->y,pt3->x,pt3->y),2)-pow(MyMath.Calculate2PtDistanceIn2D(pt1->x,pt1->y,pt2->x,pt2->y),2))/(2*MyMath.Calculate2PtDistanceIn2D(pt1->x,pt1->y,pt3->x,pt3->y)*MyMath.Calculate2PtDistanceIn2D(pt2->x,pt2->y,pt3->x,pt3->y)))

{

pt3=node->pixel

}

}

node=node->next

}

//LineList

Linelist *linehead,*linenode,*linelast

Line *ln1,*ln2,*ln3

linenode=new Linelist

linenode->line=new Line

linenode->line->pixel1=pt1

---