最优三角剖分算法

⑴ 请问什么是双向插入法

逐点插入法进行构网，分析该算法中制约运行速度的因素，在三角网拓扑关系、三角形的查找、LOP算法（Local Optimization Procere）等方面进行了优化处理，使之有了较高的执行效率。
1 数据结构
在采用逐点插入进行Delaunay三角剖分的过程中，存在大量的查询操作，利用数据结构提供三角形之间的拓扑信息，能够有效地提高算法效率。边结构应包含点与三角形的信息，三角形结构应包含点与边的信息，再考虑到构网过程中动态的数据更新，算法中采用了双向链表结构，以方便于剖分过程中新旧边（三角形）的添加、删除工作。
2 算法原理
逐点插入法是Lawson在1977年提出的，该算法思路简单，易于编程实现。基本原理为：对于已有的三角网络，向其中插入一点，该点与包含它的三角形三个顶点相连，形成三个新的三角形，然后逐个对它们进行空外接圆检测，通过交换对角线的方法来保证所形成的三角网为Delaunay三角网。
3 算法基本步骤

a: 点在三角形内，b:点在三角形边上。
在按a进行剖分时，添加了三条新边及三个新三角形NT，删除了旧三角形T。同样，在b的情况中，记录点所在的边E，根据拓扑关系，找出该边的左右相邻三角形T1，T2，添加四条新边和四个新三角形NT，删除T1，T2以及边E。
构网的关键是对新剖分出的三角形用LOP算法进行优化，其基本过程为：新三角形与周围的三角形构成共用同一条对角线的四边形，逐个对四边形中的两个三角形进行空外接圆检测，如果满足空外接圆准则，则略过。如果不满足，则用另一对角线替换现有对角线，在交换对角线后，又会产生相应的新四边形，继续进行空外接圆检测。这种方法可连续进行,直到全部满足空外接准则为止。这个过程是随着数据点P的逐次插入，与三角剖分同时进行的，可以通过递归调用优化程序来实现。

⑵ 点集的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∈E²d(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 Procere，LOP)优化使其满足Delau-nay三角网的要求，其主要步骤如下:

图3.7 Bowyer-Watson算法剖分实例

⑶ 　灵活约束的四面体化

研究几何模型所固有的物理特性的基础步骤是网格化。基于简单形体（即二维的三角形和三维的四面体）的结构化网格对自然物体的建模非常有意义，因为它们对填充由简单边界（即二维情形中的边和三维情形中的三角形）所定义的区域的物性非常灵活。

自从提出了三角剖分一个多边形的O（n）最优算法（Chazelle,1991）后，从计算几何学观点来看，二维情况现在已可以很好地处理了。二维情况的其他重要算法有Chew（1987）的约束Delaunay三角剖分方法和Ruppert（1993）的限制纵横比三角剖分方法，前者给出一种来考虑给定连接关系的点集的三角剖分方法，后者可提供一个具有“好”的元素的网格。

所以，给定一个可在其端点彼此相交的边界集合定义的二维模型，便可能利用三角形来填充所定义的全部区域。图2.13显示一个二维三角剖分模型的例子。

图2.13二维模型的三角形化（Joel Conraud,1995）

●输入角点；——输入边界

如果所有给定的边界都在剖分结果中出现，这样一个模型的三角剖分被称作约束三角剖分；而如果一些边界被分割后出现在剖分结果中，则被称为匹配三角剖分（即在这些边界上增加称为“Steiner点”的一些点，用来分割它们）。为了考虑特定的连接，Chew修改了对Delaunay三角剖分的“空球标准”的定义，在文献（Nackman and Srinivason,1991）中提出了一种在输入的边界上增加一些点，用来实现一种保持Delaunay特点的匹配三角剖分方法。

三维情况。三维情况要难得多，所有在二维中的一系列命题不再有效：

●有n条边的多边形的三角剖分有（n-2）个三角形；但一个多面体的四面体剖分具有不同数目的剖分结果。

●所有的多边形都能被三角剖分；但有些多面体，如果其边界拓扑不改变就不能被四面体填充；换句话说，不能对它们进行约束三维三角剖分，而只能进行匹配三维三角剖分的计算。

●二维点集的Delaunay三角剖分对所有可能的三角剖分的最小角最大化；但三维点集的Delaunay三角剖分做不到。

文献（Ruppert and Seidel,1992）对非凸多面体的四面体化所涉及的困难有一个更详细的描述。因此，如果有人想对一个涉及节点、边和三角形的几何结构四面体化，一般情况下需要增加Steiner点，它们的数目可能是很大的，导致一个复杂得多的四面体化。本节给出一个类似于匹配三维三角剖分的方法，称为灵活约束的四面体化方法，这将大大减少增加点的数目。

问题。首先给出一些定义，然后描述我们感兴趣的问题。

定义1：我们定义一个基于简单边界的模型M为三维空间中三角剖分定义的曲面的集合。曲面之一是闭合的，并且确定了我们研究的区域，而其他曲面有的确定了子区域，有的确定了不连续体（即面具有不闭合的边界）。它们可以沿着三角形边界彼此相交。

定义2：T是M四面体化的点集，当且仅当M的每个顶点是T的四面体顶点。

定义3：T是M的一个约束四面体化，当且仅当每个三角形（边和顶点）是T的四面体表面（边和顶点）。正如上面所写的，不总是有一个与模型相联系的四面体化。

定义4：T是M的一个匹配四面体化，当且仅当M的每个三角形：

●可以是T的四面体的表面；

●或者与T中四面体的表面有一个相关部分，这种情况意味着已增加了Steiner点。

从一个四面体化点集开始，怎样修改使它成为M的一个可接受的四面体网格呢？

初步观察。这个工作从以下观察开始：

我们想对图2.14所示的锥体作四面体化，在这个棱锥的表面，有三角形ABD、ABC。如果我们使用Delaunay三维点集三角剖分算法（例如文献Watson,1981中的方法）产生一个初始解，则由于d（C,D）＜d（A,B）（在一个Delaunay三角剖分中，对于存在4个共面点的情况，生成的两个三角形遵循对角线最短原则），在点A和B之间连接的情况不会出现在棱锥的表面拓扑中，因此也就不会出现在四面体剖分网格中。

图2.14锥体的Delaunay四面体化（Joel Conraud,1995）

如果我们忽略在共面区域这两种三角剖分的不同之处，则生成的两个四面体的边界定义了由六个输入三角形描述的同一个多面体。所以，即使拓扑结构不同，几何形态是相同的。在一个灵活约束的四面体化方法中，上面的四面体网可被视作输入的可接受的近似。

有时，一个仅基于模型节点剖分而成的四面体网格中的三角形，与一个基于边界剖分而成的三角形网格不匹配。在研究这种情况时，则可发现两种另外的方案：

（1）虽然，边界三角形与四面体网格边界不相交，但几何形状不相似。如果图2.15中的点D沿ACB面的垂线向远离E的方向移动，这个锥体被变形为一个凹多面体。按以前的情形生成两个四面体，但这时在输入三角形和四面体形成的多面体表面之间出现一个间隙。在图2.15中的截面上可看出问题。

图2.15（a）改变D的位置（Joel Conraud,1995）

图2.15（b）（E,D,C）被两个四面体共用（Joel Conraud,1995）

（2）边界三角形与一个四面体边界相交。在图2.16的例子中，三角形（C,B,E）与四面体的边（A,D）相交。要么在交点处插入一个点，要么进行网格的局部更新：可用四面体（A,B,C,E）和（D,E,B,C）代替四面体（A,D,C,E）、（A,D,E,B）、（A,D,B,C），产生一个具有8个输入三角形的网格作为四面体的表面，但该表面不再是Delaunay的了。

图2.16（a）由8个三角形确定的多面体（Joel Conraud,1995）

图2.16（b）（A,B,C,D,E,F）的Delaunay四面体化（Joel Conraud,1995）

灵活约束四面体化的目标是仅当上述的第二种情况出现时，去组合测试模型的四面体化的点集并局部修改它。对付第一种情况（不论移不移D点），则用前面的方法，通过局部网格修改或增加点来处理（George等，1991）。

下面，我们给出灵活约束四面体化的如下定义：

定义5：T是一个灵活约束四面体化，当且仅当M的每个三角形：

●可以是T的一个四面体表面；

●或者与T中四面体的表面有一个重合部分；

●或者属于一组n个相连的三角形G_Tgr1，通过它我们可以联系一组n个四面体表面G_Tetra，这组四面体与相互连接的三角形组有相同的表面。

也就是说，灵活约束的四面体化只要求约束面与四面体表面的一部分几何上重合，拓扑可以不同。

下面描述从一个模型的四面体化点集来构建一个灵活约束四面体化的算法。

算法。从一个简单边界的模型M和一个用M中的点通过无约束四面体剖分形成的四面体集T_{Point Set}（M）开始，算法分为四步：

（1）计算无约束四面体剖分T_{Point Set}（M）中与M中的三角形不匹配的三角形集合S_UnMatched

for每个M中的t

if non-matched（t,M,T_{Point Set}）/＊用判断函数判断＊/

then S_UnMatched←S_UnMatched∪｛t}

endif

endfor

这里non-matched（）函数定义如下：

（t是基于三个顶点：v₁,v₂和v₃的三角形）

TetraFaces←φ

If edge（v₁,v₂）不存在于T_{Point Set}中Retur false

else TetraFaces←TetraFaces U tetrafaces（edge（v₁,v₂））

If edge（v₂,v₃）不存在于T_{Point Set}中Return false

else TetraFaces←TetraFaces U tetrafaces（edge（v₂,v₃）

If edge（v₃,v₁）不存在于T_{Point Set}中Return false

else TetraFaces←TetraFaces U tetrafaces（edge（v₃,v₁））

Return#TetraFaces=1

（2）将S_UnMatched分为相连的三角形子集

S_Connex←φ

for S_UnMatched中的每个三角形t

Connex←｛t和M中的邻近t的三角形以及S_UnMatched的成员｝

S_UnMatched←S_UnMatchedConnex/＊把找到的三角形从S_UnMatched中去掉＊/

S_Connex←S_Connex∪｛Connex}/＊把找到的三角形合并到S_Connex中去＊/

endfor

（3）将每个S_Connex的子集分为由边界既存在于M又存在于T_{Point Set}约束的三角形的子集

S_{BoundedConnex}←φ

for S_Connex中的三角形Connex的每一个集合

for Connex的每个三角形t

BoundedConnex←｛t和Connex中邻近t的三角形以及既在M中又在T_{Point Set}没有共享边界的三角形｝

Connex←ConnexBoundedConnex/＊把找到的三角形从Connex中去掉＊/

S_{BoundedConnex}←S_{BoundedConnex}∪｛BoundedConnex｝

/＊把找到的三角形合并到S_{BoundedConnex}中去＊/

endfor

endfor

（4）对于每个子集BoundedConnex，其数据项的个数与基于同一边界的四面体表面集合的数据项个数进行对比

for S_{BoundedConnex}的每一个三角形集合BoundedConnex do

① 计算集合edges（BoundedConnex,M）和edges（BoundedConnex,T_{Point Set}）

edges（BoundedConnex,M）←φ

edges（BoundedConnex,T_{Point Set}）←φ

BoundedConnexVertices←Vertices（BoundedConnex）

for BoundedConnexVertices的每一个顶点vdo

for基于v的T_{Point Set}的每一条边e do

opposite=opposite（e,v）

if opposite是BoundedConnexVertices的成员

then edges（BoundedConnex,T_{Point Set}）←

edges（BoundedConnex,T_{Point Set}）Ue

ifM中在e和opposite之间存在边

then edges（BoundedConnex,M）←

edges（BoundedConnex,M）Ue

endif

endif

endfor

endfor

② 计算共享edges（BoundedConnex,M）的四面体表面集合TetraFaces

TetraFaces←φ

for edges（BoundedConnex,M）的每一个边e do

for基于e的T_{Point Set}的每一个四面体表面tfdo

for i from 1 to 3 do

neighber_edge←edge（tf,i）

if neighber_edge≠e

and member of edges（BoundedConnex,M）

then TetraFaces←TetraFaceUtf

endif

endfor

endfor

endfor

③ 计算在BoundedConnex的边界上的边集合BorderEdges

BorderEdges←φ

for BoundedConnex的每一个三角形t do

for i from 1 to 3 do

edge←edge（t,i）

trgl←neighber（t,i）

if trgl不是BoundedConnex的成员（或不存在）

then BorderEdges←BorderEdgesUedge

endif

endfor

endfor

④ BorderEdges定义一条闭合曲线，TetraFaces的某些元素在其外部，删除它们

for TetraFaces的每一个四面体表面tf do

for i from 1 to 3 do

e←edge（tf,i）

if e不是BorderEdges的成员

then

for基于e的每一个四面体表面tf′ do

iftf′≠tf并且tf′为TetraFaces成员

then TetraFaces

←TetraFaces

{tf′以及TetraFaces中tf′的邻接面

直到遇到BorderEdges的边界）

endif

endfor

endif

endfor

endfor

⑤ 对比＃TetraFaces和＃BoundedConnex

if#TetraFaces=BoundedConnex

then OK

else局部网格修改（或增加Stiener点或同时进行）

endif

endfor

下一节我们将给出上述实现的结果。

实验结果。上面所述算法的一个C+＋版本已被结合到GOCAD软件中，GOCAD软件是一个自然物体的三维建模工具，特别适用于石油工业中地下情况的建模。

从模型点集的一个Delaunay三维三角剖分（在浮点运算中使用Watson的算法。由Schroeder等，1990和George等，1992中所提到了同类型检查方法使得数据结构的有效性在增加点后保持不变）开始，上一节的组合算法使我们可以鉴别在T_{Point Set}和M之间的间隙。此时，还未应用上一节后面所讲的局部网格修改。在后处理过程中，通过在区域内增加一些点来使最终的网格可用。

基本的例子。如果M被简化为定义一个立方体的闭合三角剖分曲面，立方体的每个面分成两个三角形，我们可得到以下结果（图2.17）。

图2.17正方体的四面体化（Joel Conraud,1995）

小窗口中显示的四面体化的网格，元素在其中心周围收缩，以便让用户看到结构内部。这个立方体的三个可见面的三角剖分，在约束面上和四面体边界上是不同的。但两个对象定义的几何形体是一样的。

实例。彩色图版1.1是一个描绘具有超过8000个节点和17000个三角形的地质体的简单三角剖分曲面。如果这个面的网格质量好的话，是不需要对四面体网格进行修改的。

彩图版1.2中，M定义一个闭合的曲面（下面没显示），它和另外三个曲面定义四个地质层位。

彩色图版2.1表示四面体网格的一个切片，我们可以看到三个曲面定义的界面的“伤疤”。

在具有256M的HP9000/800 K200上实现第二、第三个例子需要10min和12min的CPU时间，Delaunay三维三角剖分和M的三角形与四面体表面相比较花费的时间不及总时间的1/4。

讨论。灵活约束的四面体化方法对地下情况建模的目的来说，提供了一个可接受的基于边界模型的近似：在某些情况下，以用尽可能少的元素得到四面体化网格比考虑全部的三角形更重要。如果有人想在表面三角形和四面体T_{Point Set}的边界表面之间得到1对1的关联，则有必要从T_{Point Set}中提取约束和“灵活约束”的元素，从而得到一个新的基于边界的模型M′，它与M几何形状相同，但拓扑不同。第二类型灵活约束改造的T_{Point Set}是一个M′约束的四面体化。

下一步工作将包括比较灵活约束四面体化方法与在T_{Point Set}中与M不匹配的每个三角形增加Steiner点的方法的计算效率比较。

⑷ 求一个用C++写的Delaunay三角剖分间接实现Voronoi图的代码。最好有算法说明谢谢！！ 急用！！

#include<iostream>
#include<cmath>
using namespace std;
#define N 30
typedef struct //定义点的结构体
{
int x,y;
}Point;
class point
{
private:
Point *v;
public:
int distance(Point i,Point j); //计算两点的距离
int w(Point i,Point j,Point k); //计算三条边的长度之和
void minWeightTriangulation(int n,int t[][N],int s[][N]); //用动态规划计算最优值
void print(int s[][N],int i,int j); //输出
};
int point::distance(Point i,Point j)
{
int s=(i.x-j.x)*(i.x-j.x)+(i.y-i.y)*(i.y-i.y);
return sqrt(s);
}
int point::w(Point i,Point j,Point k)
{
return distance(i,j)+distance(j,k)+distance(i,k);
}
void point::minWeightTriangulation(int n,int t[][N],int s[][N]) //用动态规划计算最优值
{
int i=0;
int r=0;
int k=0;
for(i=1;i<=n;i++) t[i][i]=0;
for(r=2;r<=n;r++)
for(i=1;i<=n-r+1;i++)
{
int j=i+r-1;
t[i][j]=t[i+1][j]+w(v[i-1],v[i],v[j]);
s[i][j]=i;
for(k=i+1;k<j;k++)
{
int u=t[i][k]+t[k+1][j]+w(v[i-1],v[k],v[j]);
if(u<t[i][j])
{
t[i][j]=u;
s[i][j]=k;
}
}
}
}
void point::print(int s[][N],int i,int j)
{
if(i==j)
return;
print(s,i,s[i][j]);
print(s,s[i][j]+1,j);
cout<<"三角行：v"<<i-1<<"v"<<s[i][j]<<"v"<<j<<endl;
}
int main()
{
int n,i;
Point v[N]={0,0};
point triangle;
int t[N][N],s[N][N];
cout<<"输入多边形的顶点数：";
cin>>n;
for(i=0;i<n;i++)
{
cout<<"输入第"<<i+1<<"点的坐标：";
cin>>v[i].x>>v[i].y;
}
triangle.minWeightTriangulation(n,t,s);
triangle.print(s,1,n);
return 0;
}