Delaunay三角剖分算法附C++实现

下面介绍Delaunay三角剖分算法：

一. 生成凸包
生成凸包的算法在我的另一个博文有详细介绍

二. 凸包切分
 在凸包链表中每次寻找一个由相邻两条凸包边组成的三角形，在该三角形的内部和边界上都不包含凸包上的任何其它点。将这个点去掉后得到新的凸包链表。重复这个过程，直到凸包链表中只剩三个离散点为止。将凸包链表中的最后三个离散点构成一个三角形，结束凸包三角剖分过程，这一过程只对凸包中的点进行处理。

List<Triangle> DivideHull(List<Point> pts)
{
List<Point> hpts;
for(int i = 0; i< pts.size() ; i++)
hpts.push_back(pts[i]);
List<Triangle> tins; //保存得到的三角形
while(hpts.size() >2) //一直到最后剩下离散的三个点
{
int tag = 0;
float minangle = 180; //每次构成相邻的边，优先找角度最小的三角形

for(int i = index; i< hpts.size() ; i++)
{
float tri_angle = 180.0;

if(i == 0){
tri_angle = Gemetry::angle3D(hpts.last(),hpts[i],hpts[i+1]);

}
else if(i == hpts.size()-1){
tri_angle = Gemetry::angle3D(hpts[i-1],hpts[i],hpts[0]);

}
else{
tri_angle = Gemetry::angle3D(hpts[i-1],hpts[i],hpts[i+1]);

}
if(tri_angle < minangle)
{
tag = i;
minangle = tri_angle;
}

}
int tagb = tag-1;
int tage = tag+1;
if(tag == 0)
tagb =hpts.size()-1;
if(tag == hpts.size()-1)
tage = 0;
tins.push_back(Triangle(hpts[tagb],hpts[tag],hpts[tage]));
hpts.removeAt(tag);

}

return tins;
}
三. 离散点内插
 算法流程：

1、选择一个尚未构成三角形的离散点

2、在已经生成的三角形中，找出该离散点的三角形，有两种情况：在三角形内部和在三角形边上。

3、如果离散点在三角形的内部，则将该三角形以及三角形的边删除，然后将三个顶点以及离散点分别连接，形成三个新的三角形。如果离散点在三角形的边上，记录点所在的边E，根据拓扑关系，找出该边的左右相邻三角形T1，T2，添加四条新边和四个新三角形NT，删除T1，T2以及边E，基本上操作是相同的。

对于新生成的三角形，不能直接加入到三角形数组里，需要挨个对其边进行空外接圆检测。具体做法为：对于新生成的三角形的边E，找出该边相邻的两个三角形，判断该边一侧的对角的顶点是否位于另外一个三角形的外接圆的里面。如果是，则将边E删除，再将两个对角连接起来，形成两个新的三角形T3，T4。

对于新生成三角形T3，T4，同样需要进行空外接圆检测，一直迭代下去，直到所有新生成的三角形都通过空外接圆检测为止。

4、重复1、2、3，直到所有非凸壳离散点都插入完为止。

List<Triangle> getDelaunay(List<Triangle> hulltins, List<Point> pts)
{

for(int i = 0; i<pts.size(); i++)
{
List<Triangle> delTin; //保存要删除的三角形
for(int j =0; j< hulltins.size(); j++)
{

if(hulltins[j].isInTriangle(pts[i]) == true) //该点在三角形内部
{
delTin.push_back(hulltins[j]);
}
if(hulltins[j].isOnTriangle(pts[i]) == true){ //点在三角形边上
delTin.push_back(hulltins[j]);

}

List<Line> borderLines;//保存离散点的相邻边
List<Triangle> newTri;//保存新得到的三角形，之后用于检测空圆
if(delTin.size() == 1) //在三角形内部
{

Line l1 = delTin[0].l1;
Line l2 = delTin[0].l2;
Line l3 = delTin[0].l3;
borderLines.push_back(l1);
borderLines.push_back(l2);
borderLines.push_back(l3);

for(int j = 0; j< borderLines.size();j++)
{
newTri.push_back(Triangle(borderLines[j].p1,borderLines[j].p2,pts[i]));
}
hulltins.removeOne(delTin[0]);
}else if(delTin.size() == 2)//在三角形边上，则要找到四条相临边，去掉重复的边
{

Line l[3];
l[0] = delTin[0].l1;
l[1] = delTin[0].l2;
l[2] = delTin[0].l3;
int index = 0;//保存重复边的下标
//四条相临边，去掉重复的边
for( int m =0;m< 3; m++)
if(delTin[1].containsLine(l[m]) ==0)
{
borderLines.push_back(l[m]);
}else
{
index = delTin[1].containsLine(l[m])-1;
}
for( int m =0;m< 3; m++)
{
if(m!=index)
borderLines.push_back(delTin[1].l[m]);
}

for(int j = 0; j< borderLines.size();j++)
{
newTri.push_back(Triangle(borderLines[j].p1,borderLines[j].p2,pts[i]));
}
hulltins.removeOne(delTin[0]);
hulltins.removeOne(delTin[1]);
}
// 对新得到的三角形进行检测空圆
delTin.clear();//之后保存需要删除的新生成的三角形
for( int s = 0; s< newTri.size(); s++){ //每一个三角形
for(int j = 0; j< 3 ; j++){ //每一条边
Line line = newTri[s].l[j];
for( int m = 0; m< hulltins.size() ; m++)//在三角形数组里搜索包含该边的三角形
{
{
if(hulltins[m].containsLine(line) )//如果三角形包含该边
{
Circle tinCircle = Circle::genTriCircle(hulltins[m]);
if(tinCircle.isInCircle(vec3(pts[i]))) //如果当前离散点在三角形的外接圆上
{
delTin.push_back(newTri[s]); //后面要删除该三角形

int x = hulltins[m].containsLine(line)-1;
//对于这个三角形的三个边，去掉重复边还有两个边，分别于当前离散点构成一个新的三角形
for( int k = 0; k<3;k ++)
{
if(x!=k){

newTri.push_back(Triangle(hulltins[m].l[k].p1,hulltins[m].l[k].p2,pts[i])); //对于
}
}
hulltins.removeAt(m);//去除该三角形

}

}
}
}
}

}
//最终把新的三角形加入到三角形数组里
for( int m =0; m< newTri.size() ; m++)
{
hulltins.push_back(newTri[m]);
}
for( int m = 0; m< delTin.size() ; m++)
{
hulltins.removeOne(delTin[m]);
}

}

return hulltins;
}