delaunay三角剖分（区域的Delaunay三角剖分）

本文目录

• 区域的Delaunay三角剖分
• 带权限定Delaunay三角化的算法步骤及实现
• MATLAB中Delaunay三角剖分怎么求各个离散点法矢，并且在图上画出法矢信息
• Delaunay三角剖分算法的准则特性
• matlab 有个函数英文看不懂，名字叫delaunay，功能是用三角形划分区域
• 基于GMaps的Delaunay三角网格剖分——逐点插入法
• 半球体在MATLAB中属于空间散乱点还是三维凸体，怎么对他进行三角剖分呢
• 点集的Delaunay三角剖分方法

区域的Delaunay三角剖分

从上文中介绍的对点集的Delaunay三角剖分方法可以知道，点集的Delaunay的三角剖分方法是最简单、最基础的，所以，可以考虑在点集的Delaunay三角剖分方法的基础上获得区域的Delaunay三角剖分，实际上，区域的Delaunay三角剖分方法正是这样做的。

为了利用点集的Delaunay三角剖分方法，首先应该将剖分对象从区域转换到点集，而为了实现上述目的，需要经过两个过程:布点、离散边界。布点即是向区域内插入一些合适的点，离散边界即是将区域边界分割成若干小线段。在进行上述处理后，就将需要剖分的对象从一个区域转变成一个点集，点集中的点为向区域内插入的点和边界离散后形成小线段的端点;然后采用点集的Delaunay三角剖分方法，如Lawson算法或Bowyer-Wat-son算法，对从区域转换过来的点集进行Delaunay三角剖分操作;最后需要删除那些在区域之外的三角形，因为对一个点集进行Delaunay三角剖分操作，最后获得的三角剖分的范围是该点集的凸壳，而点集的凸壳绝大多数情况下与需要剖分的区域不一致，通常情况下是凸壳的范围大于需要剖分的区域，因此那些属于凸壳但不属于需要剖分的区域的三角形需要删除。

区域的Delaunay三角剖分方法可以概括为3个步骤:

第一步:布点并离散边界，将剖分对象从区域转换为点集。根据剖分规模或想要得到的三角形单元的大概边长，向区域内插入一系列的内部点，并将内外边界离散成一系列的小线段。图3.14(a)为某剖分区域，图3.14(b)为向区域内布点的结果，图3.14(c)则时在布点之后进行离散边界的结果。

图3.14 二维区域布点与边界离散

第二步:对点集进行Delaunay三角剖分操作。将离散后的内外边界的端点和根据需要布设的内部点组成输入点集，采用逐点插入法(Lawson算法或Bowyer-Watson算法)将上述点集进行Delaunay三角剖分，详细步骤请参见本章中点集的Delauany三角剖分方法，最后点集的Delaunay三角剖分如图3.15所示。剖分区域之外的三角形只能是由同一条边界上线段的端点组成，不可能由内部点或不同边界(如2个顶点位于外边界，1个顶点位于内边界)上端点组成区域之外的三角形，所以，判断一个三角形是否多余，首先判断该三角形的3个点是否位于同一条边界，若是，该三角形有可能为多余三角形，但不一定，如图3.16中由点9、10和11组成的三角形位于区域内部，但由点21、22和23组成的三角形却在区域外;因此，三角形的3个顶点位于同一条边界上只是必要条件，还需要进一步判断。

图3.15 点集的Delaunay三角剖分(未进行边界处理)

第三步:边界处理。根据区域边界，删除多余的三角形，得到区域的Delaunay三角剖分。

完整的边界处理方法如下:

(1)将区域的每条边界离散后按照固定顺序排列，边界上端点的ID也依次排列。需要特别注意的是，外边界必须按照逆时针排列，外边界上的端点编号也依次由小到大排列，如图3.16所示，从点0到点36依次排列;所有内边界离散后按照顺时针排列，同样，每条内边界上端点编号依次从小到大排列。

(2)依次判断点集的Delaunay三角剖分中每个三角形的三个顶点是否位于同一边界，若某个三角形的3个顶点位于同一边界上，将3个顶点由小到大排序，并按照排序后点号组成一个新的三角形，计算该新三角形的面积，若面积＞0，则表示该三角形为区域内部三角形，保留;若面积＜0，则表示该三角形为区域外三角形，删除。

图3.16 边界处理

图3.17 二维区域的Delaunay三角剖分

如由点0、1和36组成的三角形，排序后依次为0、1、36，新的三角形面积＞0，表示该三角形为区域内部三角形，保留;由点21、22和23组成的三角形按照逆时针规则正常时排序为23、22、21，排序后依次为21、22、23，新的三角形面积＜0，表示该三角形为区域外部三角形，删除。

边界处理后，得到最终的区域的Delaunay三角剖分，如图3.17所示。

布点和离散边界的代码详见推进波前法代码中的函数CreateInnerPoint()和Create-BoundaryPoint();点集的Delauany三角剖分方法的完整代码详见Bowyer-Watson算法的代码。以下为进行边界处理的代码，函数DeletedInvalidTriangle()用于删除剖分区域之外的三角形，函数AreaOfTrgl()用于计算三角形的面积。

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

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

带权限定Delaunay三角化的算法步骤及实现

1.二维的算法步骤及实现

带权的限定Delaunay三角剖分（Weighted CDT）的算法的输入是一个包含限定线段和限定点的平面直线图（planar straight line graph，简称PSLG），算法的输出是与限定条件（限定点和限定线段）一致的一个三角形集合。

算法4.7二维的带权限定Delaunay三角剖分

ConformingWeightedDelaunayTriangulation（PSLG）

{

规范化算法

调用算法WeightAssignment2d对PSLG中的点赋权值

建立初始大三角形

调用二维带权Delaunay空洞算法（算法4.2）生成受限点集的带权Delaunay三角剖分；

调用二维恢复受限边的算法（算法4.5）生成边界一致的带权Delaunay三角网格

删除边界外的多余的三角形，得到边界一致的带权限定Delaunay三角剖分

}

2.三维的算法步骤及实现

加权的限定Delaunay剖分的算法步骤：

输入：一个分段线性复合形（a piecewise linear complex）

输出：满足受限条件的具有带权Delaunay性质的四面体集合

算法三维的带权限定Delaunay四面体剖分算法：

ConformingWeightedDelaunayTetrahedralization（PLC）

{

调用算法4.4 WeightAssignment3d对规范化后的点和边赋权值；

建立初始大四面体

for所有的PLC中的平面

将该平面通过坐标变换转换到XOY平面上；

调用算法WeightedConstrainDelaunayTriangulation进行二维的带权限定Delaunay三角剖分

将得到的二维带权限定Delaunay三角网格坐标变换到空间中其原来位置

endfor

调用算法三维带权Delaunay空洞算法生成三维点集的带权Delaunay四面体剖分

调用算法RecoverConformSegments和算法RecoverConformFacets恢复受限边和受限面

删除边界外的多余的四面体，得到边界一致的带权限定Delaunay四面体剖分

}

MATLAB中Delaunay三角剖分怎么求各个离散点法矢，并且在图上画出法矢信息

N=10;
=meshgrid(1:N);
z=peaks(N); %这里用matlab自带peaks函数产生的数据代替你的点云数据
x=x(:);
y=y(:);
z=z(:);
tri = delaunay(x,y); %三角网格划分
tr = TriRep(tri,); %三角类表示每个三角形面元
P = incenters(tr); %求每个三角面元的内接圆圆心
fn = faceNormals(tr); %求每个三角面元的法向量
trisurf(tri,x,y,z, ...
’FaceColor’, ’cyan’, ’faceAlpha’, 0.8); %画出曲面
axis equal;
hold on;
quiver3(P(:,1),P(:,2),P(:,3), ...
fn(:,1),fn(:,2),fn(:,3),1, ’color’,’r’); %画出法向量
hold off;

Delaunay三角剖分算法的准则特性

准则:
要满足Delaunay三角剖分的定义，必须符合两个重要的准则：
1、空圆特性：Delaunay三角网是唯一的（任意四点不能共圆），在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在。如下图所示：
2、最大化最小角特性：在散点集可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小角最大。从这个意义上讲，Delaunay三角网是“最接近于规则化的“的三角网。具体的说是指在两个相邻的三角形构成凸四边形的对角线，在相互交换后，六个内角的最小角不再增大。如下图所示：
特性:
以下是Delaunay剖分所具备的优异特性：
1.最接近：以最近的三点形成三角形，且各线段(三角形的边)皆不相交。
2.唯一性：不论从区域何处开始构建，最终都将得到一致的结果。
3.最优性：任意两个相邻三角形形成的凸四边形的对角线如果可以互换的话，那么两个三角形六个内角中最小的角度不会变大。
4.最规则：如果将三角网中的每个三角形的最小角进行升序排列，则Delaunay三角网的排列得到的数值最大。
5.区域性：新增、删除、移动某一个顶点时只会影响临近的三角形。
6.具有凸多边形的外壳：三角网最外层的边界形成一个凸多边形的外壳。

matlab 有个函数英文看不懂，名字叫delaunay，功能是用三角形划分区域

它存的是点的索引，所以3个数能代表3个点。比如(1, 2, 4)就是就是第1点，第2点，第4点组成的三角形。而这些"第几点"是指哪个列表里的？看注释的第一句话：

TRI = delaunay(X,Y) creates a 2-D Delaunay triangulation of the points (X,Y), where X and Y are column-vectors.

翻译一下就是：TRI = delaunay(X,Y) 会创建点集(X, Y)的2D Delaunay三角剖分，其中X和Y分别是两个列向量。

——换句话说，X和Y如果有n行，那么(X,Y)就代表n个点。而之后TRI矩阵里的数就是指向这些点的序号。1就是(X(1), Y(1)), 2就是(X(2), Y(2))...

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

A.逐点插入法基本原理

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

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

B.几何查询与拓扑修改

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

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

P∈H_1，2∩H_2，3∩H_3，1 　　　(3.2)

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

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

在找到P所在的三角形T_i后，用点P将T_i一劈为三(图3.6)。首先，新建六条半边：V₁P、V₂P、V₃P和PV₁、PV₂、PV₃；然后，需要修改三角单元V₁V₂V₃的拓扑关系：将V₁V₂的后继半边修改为V₂P，并设V₂P的后继边为PV₁，PV₁的后继边为V₁V₂；这样，PV₁V₂构成一个新的拓扑三角单元。同理，对V₂V₃和V₃V₁做拓扑修改，得到另外两个新的拓扑三角单元：PV₂V₃、PV₃V₁。

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

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

图3.6 点劈裂三角形

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

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

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

sinα=crossProduct(d₁，d₂)

sinβ=crossProduct(d₃，d₄)

cosα=scalarProduct(d₁，d₂)

cosβ=scalarProduct(d₃，d₄)

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

return ok

else

return false

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

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

C.点的动态插入

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

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

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

(2)新插入点P在原有Delaunay三角网边界之外，此时只需要将新数据点追加到离散点数据库的最后。从Delaunay三角网边界为凸包的角度考虑，先找出已有三角网的边界边，可以通过α₂(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={ q₀，q₁，…，q_k-1，q_k=q₀}。按照逆时针顺序形成该多边形的凸耳权值队列。从队列中选择权值最小的凸耳ear={q_i，q_i+1，q_i+2}，则q_iq_i+2为D-TIN中的一条边，凸耳权计算公式为(Devillers，1999)：

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

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

删除点的步骤如下：

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

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

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

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

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

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

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

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

半球体在MATLAB中属于空间散乱点还是三维凸体，怎么对他进行三角剖分呢

可以。先用用delaunay三角剖分，然后用trimesh命令显示。假设你的三维散点的空间坐标分别存在向量x，y，z（列向量）中。照如下方式操作。tri=delaunay(x,y)%将散点在XoY平面做delaunay三角剖分。trimesh(tri,x,y,z);%显示曲面，利用上一步的三角剖分结果，将空间的三角形画出来。如果效果不好，那是因为你的点云数据不理想，比如不够稠密，不够均匀。可以通过插值解决。

点集的Delaunay三角剖分方法

3.2.1.1 基本理论

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

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

3.2.1.1.1 Voronoi图

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

设x为平面上的点，则:

区域V_(i)={x∈Ed(x，v_i)≤d(x，v_j)，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，1981;Watson，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()的作用是判断点(x_p，y_p)是否落在由三点(x₁，y₁)、(x₂，y₂)和(x₃，y₃)组成的三角形的外接圆内。

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

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

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

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

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

图3.7 Bowyer-Watson算法剖分实例