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事件添加代码如下：

private void 获取点集坐标数组_Click(object sender, EventArgs e)

AcadSelectionSet mySelectionSet;

mySelectionSet = AcadDoc.SelectionSets.Add("NewSelectionSet001");

Int16[] FilterType = new Int16[1];

object[] FilterData = new object[1];

FilterType[0] = 0;

FilterData[0] = "POINT";

mySelectionSet.SelectOnScreen(FilterType, FilterData);

double[] arrPoints = new double[3 * mySelectionSet.Count];

int count = 0;

foreach (AcadObject acadObj in mySelectionSet)

{

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相同。

为了数据查询、调用方便以及数据结构的管理，在调用任何一个点、边或三角形时只需要通过索引即可调取。故在全局变量中定义如下类用于管理边和三角形：

public class Edge

{

public int Start;//边的起点

public int End;//边的终点

public int LeftTri = -1;//左边三角形索引

public int RightTri = -1;//右边三角形索引

}

public class Tri

{

public int NodeA;//第一个节点的索引

public int NodeB;//第二个节点的索引

public int NodeC;//第三个节点的索引

public int AdjTriA = -1;//第一个邻接三角形索引

public int AdjTriB = -1;//第二个邻接三角形索引

public int AdjTriC = -1;//第三个邻接三角形索引

}

此外，在全局变量中定义边和三角形数组（由于在此用到了ArrayList，所以需要添加引用“using System.Collections;”），其代码如下：

private ArrayList arrayEdges = new ArrayList();

private ArrayList arrayTris = new ArrayList();

下面从startPoint开始通过寻找最大夹角来获取外围边界点，其代码如下：

//获取点集外围边界

int i, startIndex = 0, tempIndex, lastIndex, pointCount = arrPoints.Length / 3;

for (i = 1; i < pointCount; i++) //寻找X值最小的点号

{

if (arrPoints[3 * i] < arrPoints[3 * startIndex])

{

startIndex = i;

}

}

Edge edge = new Edge();

edge.Start = startIndex;

lastIndex = startIndex - 1;

tempIndex = startIndex;

double[] vector1 = new double[2], vector2 = new double[2];

vector1[0] = 0; vector1[1] = 100;

double vector1Length, 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 < pointCount; 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;

}

else if (angleTemp == angleMax && vector2Length < lengthMin)

{

edge.End = i;

lengthMin = vector2Length;

}

}

}

arrayEdges.Add(edge);

lastIndex=edge.End;

edge = new Edge();

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 = new Edge[arrayEdges.Count];

arrayEdges.CopyTo(edges);

for(i=0;i<arrayEdges.Count;i++)

{

double[] p1=new double[3],p2=new double[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

int newIndex, edgeCount = arrayEdges.Count;

double angleV1V2Temp, angleV2V3Temp, angleV1V2Max, angleV2V3Max, vector3Length;

double[] vector3 = new double[2];

bool isTriExist;

for (i = 0; i < arrayEdges.Count; i++)

{

edge=new Edge();

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 (int j = 0; j < pointCount; 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;

}

else if (angleV2V3Temp == angleV2V3Max && angleV1V2Max >= angleV1V2Temp)

{

angleV1V2Max = angleV1V2Temp;

newIndex = j;

}

}

}

}

if (newIndex != -1)//若找到了这么一个满足要求的点就记录三角形

{

Tri tri = new Tri();

tri.NodeA = edge.Start;

tri.NodeB = edge.End;

tri.NodeC = newIndex;

edge.LeftTri = arrayTris.Count;//设置边的左侧三角形索引

isTriExist = false;

//记录边1

for (int k = 0; k < arrayEdges.Count; k++)

{

Edge tempEdge = (Edge)arrayEdges[k];

if (tempEdge.Start == edge.Start && tempEdge.End == newIndex)

{

tempEdge.RightTri = arrayTris.Count;

tri.AdjTriB = tempEdge.LeftTri;

isTriExist = true;

break;

}

else if (tempEdge.Start == newIndex && tempEdge.End == edge.Start)

{

tempEdge.LeftTri = arrayTris.Count;

tri.AdjTriB = tempEdge.RightTri;

isTriExist = true;

break;

}

}

if (isTriExist == false)//若不存在这条边就新建一条边

{

Edge newEdge = new Edge();

newEdge.Start = newIndex;

newEdge.End = edge.Start;

newEdge.LeftTri = arrayTris.Count;

arrayEdges.Add(newEdge);

}

isTriExist = false;

//记录边2

for (int k = 0; k < arrayEdges.Count; k++)

{

Edge tempEdge = (Edge)arrayEdges[k];

if (tempEdge.Start == newIndex && tempEdge.End == edge.End)

{

tempEdge.RightTri = arrayTris.Count;

tri.AdjTriA = tempEdge.LeftTri;

isTriExist = true;

break;

}

else if (tempEdge.Start == edge.End && tempEdge.End == newIndex)

{

tempEdge.LeftTri = arrayTris.Count;

tri.AdjTriA = tempEdge.RightTri;

isTriExist = true;

break;

}

}

if (isTriExist == false)//若不存在这条边就新建一条边

{

Edge newEdge = new Edge();

newEdge.Start = edge.End ;

newEdge.End = newIndex;

newEdge.LeftTri = arrayTris.Count;

arrayEdges.Add(newEdge);

}

tri.AdjTriC = edge.RightTri;//如果edge的右三角形不存在，则左三角也不存在

arrayTris.Add(tri);

}

}

else if (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 (int j = 0; j < pointCount; 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] * vector3[1])

/ (vector2Length * vector3Length));

if (angleV2V3Temp > angleV2V3Max)

{

angleV1V2Max = angleV1V2Temp;

angleV2V3Max = angleV2V3Temp;

newIndex = j;

}

else if (angleV2V3Temp == angleV2V3Max && angleV1V2Max <= angleV1V2Temp)

{

angleV1V2Max = angleV1V2Temp;

newIndex = j;

}

}

}

}

if (newIndex != -1)//若找到了这么一个满足要求的点就记录三角形

{

Tri tri = new Tri();

tri.NodeA = edge.Start;

tri.NodeB = edge.End;

tri.NodeC = newIndex;

edge.RightTri = arrayTris.Count;//设置边的左侧三角形索引

isTriExist = false;

//记录边1

for (int k = 0; k < arrayEdges.Count; k++)

{

Edge tempEdge = (Edge)arrayEdges[k];

if (tempEdge.Start == newIndex && tempEdge.End ==edge.Start )

{

tempEdge.RightTri = arrayTris.Count;

tri.AdjTriB = tempEdge.LeftTri;

isTriExist = true;

break;

}

else if (tempEdge.Start == edge.Start && tempEdge.End == newIndex)

{

tempEdge.LeftTri = arrayTris.Count;

tri.AdjTriB = tempEdge.RightTri;

isTriExist = true;

break;

}

}

if (isTriExist == false)//若不存在这条边就新建一条边

{

Edge newEdge = new Edge();

newEdge.Start = edge.Start;

newEdge.End = newIndex;

newEdge.LeftTri = arrayTris.Count;

arrayEdges.Add(newEdge);

}

isTriExist = false;

//记录边2

for (int k = 0; k < arrayEdges.Count; k++)

{

Edge tempEdge = (Edge)arrayEdges[k];

if (tempEdge.Start == edge.End && tempEdge.End == newIndex)

{

tempEdge.RightTri = arrayTris.Count;

tri.AdjTriA = tempEdge.LeftTri;

isTriExist = true;

break;

}

else if (tempEdge.Start == newIndex && tempEdge.End == edge.End)

{

tempEdge.LeftTri = arrayTris.Count;

tri.AdjTriA = tempEdge.RightTri;

isTriExist = true;

break;

}

}

if (isTriExist == false)//若不存在这条边就新建一条边

{

Edge newEdge = new Edge();

newEdge.Start = newIndex;

newEdge.End = edge.End;

newEdge.LeftTri = arrayTris.Count;

arrayEdges.Add(newEdge);

}

tri.AdjTriC = edge.LeftTri;//如果edge的右三角形不存在，则左三角也不存在

arrayTris.Add(tri);

}

}

}

其中，主要通过边的左侧三角形索引和右侧三角形索引来判断是否需要通过该边来创建三角形。若该边一侧的三角形索引为-1则表示这一侧没有三角形，然后通过搜寻新的符合Delaunay准则的点形成新的三角形，同时判断新形成的两条边本身是否存在。若该边已经存在，则设置该边的左侧或右侧三角形索引为新创建的三角形索引；若该边不存在，则创建一条新的边，同时设置改边的左侧或右侧三角形索引。整个程序中形成三角形的过程为从边界线向内扩展的方式，如下示意图所示：