C# Delaunay三角剖分.docx
《C# Delaunay三角剖分.docx》由会员分享，可在线阅读，更多相关《C# Delaunay三角剖分.docx（19页珍藏版）》请在冰豆网上搜索。
C#Delaunay三角剖分
Delaunay三角剖分
在实际中运用的最多的三角剖分是Delaunay三角剖分。
首先，我们来了解一下Delaunay边。
Delaunay边的定义为：
假设E中的一条边e（其端点为a,b），若e满足条件：
存在一个圆经过a,b两点，圆内不含点集中任何其他的点，这一特性又称空圆特性，则称之为Delaunay边：
Delaunay三角剖分的定义为：
如果点集的一个三角剖分只包含Delaunay边，那么该三角剖分称为Delaunay三角剖分。
要满足Delaunay三角剖分的定义，必须符合下面两个重要的准则：
1）空圆特性：
Delaunay三角网是唯一的，在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在；
2）最大化最小角特性：
在散点集可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小角最大。
从这个意义上讲，Delaunay三角网是“最接近于规则化的”的三角网。
具体来说是指在两个相邻的三角形构成凸四边形的对角线，在相互交换后，六个内角的最小角不再增大。
经典的Delaunay剖分算法主要有两类[1]：
1）增量算法：
又称为Delaunay空洞算法或加点法，其思路为从一个三角形开始，每次增加一个点，保证每一步得到的当前三角形是局部优化的三角形。
2）局部变换法：
又称为换边或换面法，其思路为构造非优化的三角网，然后对两个共边三角形形成的凸四边形迭代换边优化。
迄今为止关于Delaunay剖分已经出现了很多算法，主要有分治算法、逐步插入法、三角网生长法等。
其中三角网生长算法由于效率较低,目前较少采用;分治算法最为高效,但算法相对比较复杂;逐点插入法实现简单,但它的时间复杂度差[2]。
特别是近些年，随着计算机水平的不断提升，又出现了各种各样的改进算法。
本节将主要根据逐步插入法的原理，通过对给予的数据高程点进行Delaunay三角剖分。
其基本步骤为：
1）获取点集坐标数组；2）获取点集外围边界；3）根据边界及内部点生成三角网。
现新建一个项目，并命名为“Delaunay三角剖分”，同时添加相关的引用及定义“启动CAD（）”函数。
1、获取点集坐标数组
添加一个按钮，设置其Name和Text属性都为“获取点集坐标数组”，为其Click事件添加代码如下：
privatevoid获取点集坐标数组_Click（objectsender,EventArgse）
{
Microsoft.VisualBasic.Interaction.AppActivate（AcadApp.Caption）;
AcadSelectionSetmySelectionSet;
mySelectionSet=AcadDoc.SelectionSets.Add（"NewSelectionSet001"）;
Int16[]FilterType=newInt16[1];
object[]FilterData=newobject[1];
FilterType[0]=0;
FilterData[0]="POINT";
mySelectionSet.SelectOnScreen（FilterType,FilterData）;
double[]arrPoints=newdouble[3*mySelectionSet.Count];
intcount=0;
foreach（AcadObjectacadObjinmySelectionSet）
{
if（acadObj.ObjectName=="AcDbPoint"）
{
count++;
double[]PointCoord;
PointCoord=（Double[]）（（（AcadPoint）acadObj）.Coordinates）;
arrPoints[3*count-3]=PointCoord[0];
arrPoints[3*count-2]=PointCoord[1];
arrPoints[3*count-1]=PointCoord[2];
}
}
MessageBox.Show（"共选择点的个数为：
"+count.ToString（））;
AcadDoc.SelectionSets.Item（"NewSelectionSet001"）.Delete（）;
}
其中，arrPoints数组保存各点的坐标值。
2、获取点集外围边界
根据前面获取的外围点集来获取外围边界，其主要思路为：
1）搜寻X坐标最小的点，在此称为startPoint，从该点开始搜寻下一点，本例中以逆时针方向来搜寻边界；2）现定义一个虚拟点（该点并不存在）位置为startPoint正北方向，该点与startPoint形成一个向量Vector1，下一点（判断是否为边界点的点，或称选取的点）与startPoint也形成一个向量Vector2，则判断该选取的点是外围边界点所要满足的条件为向量Vector2与向量Vector1之间的夹角最大；3）按同样的方式，新形成的向量Vector1为新选取的外围边界点与上一个外围边界点形成的向量，向量Vector2为选取的点与新选取的外围边界点形成的向量，判断该选取的点是外围边界点所要满足的条件仍然为向量Vector2与向量Vector1之间的夹角最大；4）按此思路循环进行，直到找到一个边界点与startPoint相同。
为了数据查询、调用方便以及数据结构的管理，在调用任何一个点、边或三角形时只需要通过索引即可调取。
故在全局变量中定义如下类用于管理边和三角形：
publicclassEdge
{
publicintStart;//边的起点
publicintEnd;//边的终点
publicintLeftTri=-1;//左边三角形索引
publicintRightTri=-1;//右边三角形索引
}
publicclassTri
{
publicintNodeA;//第一个节点的索引
publicintNodeB;//第二个节点的索引
publicintNodeC;//第三个节点的索引
publicintAdjTriA=-1;//第一个邻接三角形索引
publicintAdjTriB=-1;//第二个邻接三角形索引
publicintAdjTriC=-1;//第三个邻接三角形索引
}
此外，在全局变量中定义边和三角形数组（由于在此用到了ArrayList，所以需要添加引用“usingSystem.Collections;”），其代码如下：
privateArrayListarrayEdges=newArrayList（）;
privateArrayListarrayTris=newArrayList（）;
下面从startPoint开始通过寻找最大夹角来获取外围边界点，其代码如下：
//获取点集外围边界
inti,startIndex=0,tempIndex,lastIndex,pointCount=arrPoints.Length/3;
for（i=1;i{
if（arrPoints[3*i]{
startIndex=i;
}
}
Edgeedge=newEdge（）;
edge.Start=startIndex;
lastIndex=startIndex-1;
tempIndex=startIndex;
double[]vector1=newdouble[2],vector2=newdouble[2];
vector1[0]=0;vector1[1]=100;
doublevector1Length,vector2Length,angleTemp,angleMax,lengthMin;
angleMax=0;
while（lastIndex!
=startIndex）
{
vector1Length=Math.Sqrt（vector1[0]*vector1[0]+vector1[1]*vector1[1]）;
lengthMin=300;
for（i=0;i{
if（i!
=edge.Start）
{
vector2[0]=arrPoints[3*i]-arrPoints[3*tempIndex];
vector2[1]=arrPoints[3*i+1]-arrPoints[3*tempIndex+1];
vector2Length=Math.Sqrt（vector2[0]*vector2[0]+vector2[1]*vector2[1]）;
angleTemp=Math.Acos（（vector1[0]*vector2[0]+vector1[1]*vector2[1]）/（vector1Length*vector2Length））;
if（angleTemp>angleMax）
{
angleMax=angleTemp;
edge.End=i;
lengthMin=vector2Length;
}
elseif（angleTemp==angleMax&&vector2Length{
edge.End=i;
lengthMin=vector2Length;
}
}
}
arrayEdges.Add（edge）;
lastIndex=edge.End;
edge=newEdge（）;
edge.Start=lastIndex;
angleMax=0;
vector1[0]=arrPoints[3*tempIndex]-arrPoints[3*lastIndex];
vector1[1]=arrPoints[3*tempIndex+1]-arrPoints[3*lastIndex+1];
tempIndex=lastIndex;
}
其中，Vector1和Vector2分别表示两个向量，根据向量的数量积公式即可求出两向量间的夹角，该公式如下：
为了显示该边界线，还可以绘制出该边界线，其代码如下：
//绘制边界
Edge[]edges=newEdge[arrayEdges.Count];
arrayEdges.CopyTo（edges）;
for（i=0;i{
double[]p1=newdouble[3],p2=newdouble[3];
p1[0]=arrPoints[3*edges[i].Start];
p1[1]=arrPoints[3*edges[i].Start+1];
p1[2]=arrPoints[3*edges[i].Start+2];
p2[0]=arrPoints[3*edges[i].End];
p2[1]=arrPoints[3*edges[i].End+1];
p2[2]=arrPoints[3*edges[i].End+2];
AcadDoc.ModelSpace.AddLine（p1,p2）;
}
运行程序，其显示结果如下图所示：
3、根据边界及内部点生成三角网
生成三角网主要是根据边界及内部点，先从边界开始搜寻满足Delaunay准则的点形成三角形。
由于前面已约定边界搜索方式按逆时针方向，所以新形成的三角形都在边的左侧。
在新形成的三角形中又包含新的边，然后新的边又可以继续搜索新的点来形成三角形，直到最后所有点都被搜索完为止。
现取一条边来分析，如下图所示，start—end是已经存在的一条边。
首先要对选取的点P1进行判断，判断其在向量Vector1的左侧还是右侧。
由于本例中已约定按逆时针方向来搜寻，即每次均在边的左侧生成新的三角网，所以如果点在向量Vector1的右侧，则不考虑该点。
若点在向量Vector1的左侧，则计算点P1与边的start点和end点形成的两个向量Vector2和Vector3，同时计算Vector1和Vector2以及Vector2和Vector3之间的夹角。
生成TIN的实现代码如下：
//生成TIN
intnewIndex,edgeCount=arrayEdges.Count;
doubleangleV1V2Temp,angleV2V3Temp,angleV1V2Max,angleV2V3Max,vector3Length;
double[]vector3=newdouble[2];
boolisTriExist;
for（i=0;i{
edge=newEdge（）;
edge=（Edge）arrayEdges[i];//取出一条边
//判断三角形是否存在（若本边的左三角已存在，则计算右三角）？
if（edge.LeftTri==-1）
{
newIndex=-1;//新选取点的索引
angleV1V2Max=0;angleV2V3Max=0;
vector1[0]=arrPoints[3*edge.End]-arrPoints[3*edge.Start];
vector1[1]=arrPoints[3*edge.End+1]-arrPoints[3*edge.Start+1];
for（intj=0;j{
if（j!
=edge.Start&&j!
=edge.End）//排除边的端点
{
vector2[0]=arrPoints[3*j]-arrPoints[3*edge.Start];
vector2[1]=arrPoints[3*j+1]-arrPoints[3*edge.Start+1];
if（vector1[0]*vector2[1]-vector1[1]*vector2[0]>0）//如果点j在向量vector1的左侧
{
vector3[0]=arrPoints[3*j]-arrPoints[3*edge.End];
vector3[1]=arrPoints[3*j+1]-arrPoints[3*edge.End+1];
vector1Length=Math.Sqrt（vector1[0]*vector1[0]+vector1[1]*vector1[1]）;
vector2Length=Math.Sqrt（vector2[0]*vector2[0]+vector2[1]*vector2[1]）;
vector3Length=Math.Sqrt（vector3[0]*vector3[0]+vector3[1]*vector3[1]）;
angleV1V2Temp=Math.Acos（（vector1[0]*vector2[0]+vector1[1]*vector2[1]）
/（vector1Length*vector2Length））;
angleV2V3Temp=Math.Acos（（vector2[0]*vector3[0]+vector2[1]*vector3[1]）
/（vector2Length*vector3Length））;
if（angleV2V3Temp>angleV2V3Max）
{
angleV2V3Max=angleV2V3Temp;
angleV1V2Max=angleV1V2Temp;
newIndex=j;
}
elseif（angleV2V3Temp==angleV2V3Max&&angleV1V2Max>=angleV1V2Temp）
{
angleV1V2Max=angleV1V2Temp;
newIndex=j;
}
}
}
}
if（newIndex!
=-1）//若找到了这么一个满足要求的点就记录三角形
{
Tritri=newTri（）;
tri.NodeA=edge.Start;
tri.NodeB=edge.End;
tri.NodeC=newIndex;
edge.LeftTri=arrayTris.Count;//设置边的左侧三角形索引
isTriExist=false;
//记录边1
for（intk=0;k{
EdgetempEdge=（Edge）arrayEdges[k];
if（tempEdge.Start==edge.Start&&tempEdge.End==newIndex）
{
tempEdge.RightTri=arrayTris.Count;
tri.AdjTriB=tempEdge.LeftTri;
isTriExist=true;
break;
}
elseif（tempEdge.Start==newIndex&&tempEdge.End==edge.Start）
{
tempEdge.LeftTri=arrayTris.Count;
tri.AdjTriB=tempEdge.RightTri;
isTriExist=true;
break;
}
}
if（isTriExist==false）//若不存在这条边就新建一条边
{
EdgenewEdge=newEdge（）;
newEdge.Start=newIndex;
newEdge.End=edge.Start;
newEdge.LeftTri=arrayTris.Count;
arrayEdges.Add（newEdge）;
}
isTriExist=false;
//记录边2
for（intk=0;k{
EdgetempEdge=（Edge）arrayEdges[k];
if（tempEdge.Start==newIndex&&tempEdge.End==edge.End）
{
tempEdge.RightTri=arrayTris.Count;
tri.AdjTriA=tempEdge.LeftTri;
isTriExist=true;
break;
}
elseif（tempEdge.Start==edge.End&&tempEdge.End==newIndex）
{
tempEdge.LeftTri=arrayTris.Count;
tri.AdjTriA=tempEdge.RightTri;
isTriExist=true;
break;
}
}
if（isTriExist==false）//若不存在这条边就新建一条边
{
EdgenewEdge=newEdge（）;
newEdge.Start=edge.End;
newEdge.End=newIndex;
newEdge.LeftTri=arrayTris.Count;
arrayEdges.Add（newEdge）;
}
tri.AdjTriC=edge.RightTri;//如果edge的右三角形不存在，则左三角也不存在
arrayTris.Add（tri）;
}
}
elseif（edge.RightTri==-1）
{
newIndex=-1;//新选取点的索引
angleV1V2Max=0;angleV2V3Max=0;
vector1[0]=arrPoints[3*edge.End]-arrPoints[3*edge.Start];
vector1[1]=arrPoints[3*edge.End+1]-arrPoints[3*edge.Start+1];
for（intj=0;j{
if（j!
=edge.Start&&j!
=edge.End）//排除边的端点
{
vector2[0]=arrPoints[3*j]-arrPoints[3*edge.Start];
vector2[1]=arrPoints[3*j+1]-arrPoints[3*edge.Start+1];
if（vector1[0]*vector2[1]-vector1[1]*vector2[0]<0）
{
vector3[0]=arrPoints[3*j]-arrPoints[3*edge.End];
vector3[1]=arrPoints[3*j+1]-arrPoints[3*edge.End+1];
vector1Length=Math.Sqrt（vector1[0]*vector1[0]+vector1[1]*vector1[1]）;
vector2Length=Math.Sqrt（vector2[0]*vector2[0]+vector2[1]*vector2[1]）;
vector3Length=Math.Sqrt（vector3[0]*vector3[0]+vector3[1]*vector3[1]）;
angleV1V2Temp=Math.Acos（（vector1[0]*vector2[0]+vector1[1]*vector2[1]）
/（vector1Length*vector2Length））;
angleV2V3Temp=Math.Acos（（vector2[0]*vector3[0]+vector2[1]*
升级会员 