Delaunay算法简介及实现

转自：http://blog.csdn.net/chinamming/article/details/16874371

三角剖分：假设V是二维实数域上的有限点集，边e是由点集中的点作为端点构成的封闭线段，E为e的集合，则，点集V的一个三角剖分T=（V,E)是一个平面图G,该平面图满足的条件：

1.除了端点，平面图中的边不包含点集中的任何点。
2.没有相交边。
3.平面图中所有的面都是三角面，且所有三角面的合集是散点集V的凸包。

在实际中运用的最多的三角剖分是Delaunay三角剖分，它是一种特殊的三角剖分。

算法实现：

delaunay.h

#ifndef DELAUNAY_H_INCLUDED 
#define DELAUNAY_H_INCLUDED 
#include <cstdlib> 
#include <iostream> 
#include <cstring> 
#include <string> 
#include <fstream> 
#include <math.h> 
#include <vector>

using namespace std;

typedef struct 
{ 
    double x; 
    double y; 
    double z; 
}Point;//定义点类


typedef vector<Point> PointArray;//定义点类的vector容器

typedef struct 
{ 
    int left; 
    int right; 
    int count;//边的计数，如果计数为0，则删除此边 
}Edge;//定义边类


typedef vector<Edge> EdgeArray;//定义边类的vector容器

typedef struct 
{ 
    int v[3];//三角形的三个顶点 
    Edge s[3];//三角形的三条边 
    double xc;//三角形外接圆圆心的x坐标 
    double yc;//三角形外接圆圆心的y坐标 
    double r;//三角形外接圆的半径 
}Triangle;//定义三角形类


typedef vector<Triangle> TriangleArray;//定义三角形类的vector容器

typedef vector<int> intArray;//定义int类的vector容器

class Delaunay//定义Delaunay类 
{ 
public: 
    Delaunay(Point p1,Point p2,Point p3,Point p4); //Delaunay类的构造函数，创建外边框 
    ~Delaunay();//Delaunay类的析构函数

    bool AddPoint(double xx,double yy,double zz);//向已有剖分图形中加点的函数 
    void Delete_Frame();//删除外边框 
    void Boundary_Recover(int fromPoint,int toPoint);//边界恢复 
    void output();//输出ANSYS命令流文件 
private: 
    void Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3);//计算三角形的外接圆圆心坐标和半径 
    void MakeTriangle(int n1,int n2,int n3);//生成指定顶点的三角形 
    bool inCircle(double xx,double yy,Triangle currentTris);//判断点是否在圆内 
    void DelTriangle(int n,EdgeArray &BoundEdges);//删除指定的三角形

    PointArray m_Pts;//m_Pts用于存储所有点 
    EdgeArray m_Edges;//m_Edges用于存储所有边 
    TriangleArray m_Tris;//m_Tris用于存储所有三角形 
}; 
void GetPoint(double &xx,double &yy,double &zz,string line);//解析从input文件中读取的每一行数据 
#endif // DELAUNAY_H_INCLUDED

delaunay.cpp

#include "delaunay.h"

Delaunay::Delaunay(Point p1,Point p2,Point p3,Point p4) 
{ 
    m_Pts.resize(4); 
    m_Pts[0]=p1; 
    m_Pts[1]=p2; 
    m_Pts[2]=p3; 
    m_Pts[3]=p4;//添加四个外边框点 
    m_Edges.resize(4); 
    Edge l1={0,1,-1}; 
    Edge l2={1,2,-1}; 
    Edge l3={0,3,-1}; 
    Edge l4={2,3,-1}; 
    m_Edges[0]=l1; 
    m_Edges[1]=l2; 
    m_Edges[2]=l3; 
    m_Edges[3]=l4;//添加四个外边框的边 
    MakeTriangle(0,1,2); 
    MakeTriangle(0,2,3);//添加初始的两个三角形 
}

Delaunay::~Delaunay()//清空Delaunay类的数据成员 
{ 
    m_Pts.resize(0); 
    m_Edges.resize(0); 
    m_Tris.resize(0); 
}

void Delaunay::MakeTriangle(int n1,int n2,int n3) 
{ 
    double x_centre,y_centre,radius; 
    Cal_Centre(x_centre,y_centre,radius,n1,n2,n3);//获得顶点为n1,n2,n3的三角形的外接圆圆心坐标和半径 
    Triangle newTriangle={{n1,n2,n3},{{n1,n2,1},{n2,n3,1},{n1,n3,1}},x_centre,y_centre,radius};//生成指定的三角形 
    m_Tris.push_back(newTriangle);//向m_Tris中添加新构造的三角形 
    int EdgeSzie=(int)m_Edges.size();//获得目前的边数 
    int flag; 
    for (int i=0;i<3;i++) 
    { 
        flag=1; 
        for(int j=0;j<EdgeSzie;j++)//通过循环判断新构造的三角形的各边是否已经存在于m_Edges中，如果存在则只增加该边的计数，否则添加新边 
        { 
            if (newTriangle.s[i].left==m_Edges[j].left&&newTriangle.s[i].right==m_Edges[j].right&&

            m_Edges[j].count!=-1) 

            {

                 flag=0;m_Edges[j].count+=1;break;

            } 
            else if(newTriangle.s[i].left==m_Edges[j].left&&newTriangle.s[i].right==m_Edges[j].right&& 

                      m_Edges[j].count==-1)

                   {

                       flag=0;

                       break;

                 } 
        } 
        if (flag==1) m_Edges.push_back(newTriangle.s[i]); 
    } 
}

void Delaunay::Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3) 
{ 
    double x1,x2,x3,y1,y2,y3; 
    x1=m_Pts[n1].x; 
    y1=m_Pts[n1].y; 
    x2=m_Pts[n2].x; 
    y2=m_Pts[n2].y; 
    x3=m_Pts[n3].x; 
    y3=m_Pts[n3].y; 
    x_centre=((y2-y1)*(y3*y3-y1*y1+x3*x3-x1*x1)-(y3-y1)*(y2*y2-y1*y1+x2*x2-x1*x1))/(2*(x3-x1)*(y2-y1)-2*((x2-x1)*(y3-y1)));//计算外接圆圆心的x坐标 
    y_centre=((x2-x1)*(x3*x3-x1*x1+y3*y3-y1*y1)-(x3-x1)*(x2*x2-x1*x1+y2*y2-y1*y1))/(2*(y3-y1)*(x2-x1)-2*((y2-y1)*(x3-x1)));//计算外接圆圆心的y坐标 
    radius= sqrt((x1 - x_centre)*(x1 - x_centre) + (y1 - y_centre)*(y1 - y_centre));//计算外接圆的半径 
}

bool Delaunay::AddPoint(double xx,double yy,double zz) 
{ 
    EdgeArray BoundEdges;//BoundEdges用于存储在删除三角形后留下的边框，用于构造新的三角形 
    Point newPoint={xx,yy,zz}; 
    m_Pts.push_back(newPoint);//向m_Pts中添加新点 
    intArray badTriangle;//badTriangle用于存储不符合空圆规则的三角形的索引号 
    int TriSize=(int)m_Tris.size();//获得目前的三角形数 
    for (int i=0;i<TriSize;i++)//通过循环找到所有不符合空圆规则的三角形，并将其索引号存在badTriangle中 
    { 
        if (inCircle(xx,yy,m_Tris[i])==true) badTriangle.push_back(i); 
    } 
    for (int i=0;i<(int)badTriangle.size();i++)//通过循环删除所有不符合空圆规则的三角形，同时保留边框 
    { 
        DelTriangle(badTriangle[i],BoundEdges); 
        for (int j=i+1;j<(int)badTriangle.size();j++) badTriangle[j]-=1; 
    } 
    int PtSize=(int)m_Pts.size();//获得目前的点数 
    for (int i=0;i<(int)BoundEdges.size();i++)//生成新的三角形 
    { 
        if (PtSize-1<BoundEdges[i].left) MakeTriangle(PtSize-1,BoundEdges[i].left,BoundEdges[i].right); 
        else if (PtSize-1>BoundEdges[i].left && PtSize-1<BoundEdges[i].right) MakeTriangle(BoundEdges[i].left,PtSize-1,BoundEdges[i].right); 
        else MakeTriangle(BoundEdges[i].left,BoundEdges[i].right,PtSize-1); 
    } 
    return true; 
}

bool Delaunay::inCircle(double xx,double yy,Triangle currentTris)//判断点是否在三角形的外接圆内 
{ 
    double dis=sqrt((currentTris.xc-xx)*(currentTris.xc-xx)+(currentTris.yc-yy)*(currentTris.yc-yy)); 
    if (dis>currentTris.r) return false; 
    else return true; 
}

void Delaunay::DelTriangle(int n,EdgeArray &BoundEdges) 
{ 
    for (int i=0;i<3;i++) 
    { 
        for (int j=0;j<(int)m_Edges.size();j++) 
        { 
            if (m_Edges[j].left==m_Tris[n].s[i].left&&m_Edges[j].right==m_Tris[n].s[i].right) 
            { 
                if (m_Edges[j].count==2)//若要删除三角形的一边的计数为2，则将其计数减1，并将其压入BoundEdges容器中 
                { 
                    m_Edges[j].count=1; 
                    BoundEdges.push_back(m_Edges[j]); 
                } 
                else if (m_Edges[j].count==-1) BoundEdges.push_back(m_Edges[j]);//如果是外边框，则直接压入BoundEdges容器中 
                else if (m_Edges[j].count==1)//如果删除三角形的一边的计数为1，则删除该边，同时查看BoundEdges中是否有此边，若有，则删除 
                { 
                    for (int k=0;k<(int)BoundEdges.size();k++) 
                    { 
                        if (BoundEdges[k].left==m_Edges[j].left&&BoundEdges[k].right==m_Edges[j].right) 
                        { 
                            BoundEdges.erase(BoundEdges.begin()+k); 
                            break; 
                        } 
                    } 
                    m_Edges.erase(m_Edges.begin()+j); 
                    j--; 
                } 
                break; 
            } 
        } 
    } 
    m_Tris.erase(m_Tris.begin()+n);//删除该三角形 
}

void Delaunay::output()//向“output.log"文件中写入ANSYS命令流 
{ 
    ofstream outfile("output.log"); 
    if (!outfile) 
    { 
        cout<<"Unable to output nodes!"; 
        exit(1); 
    } 
    outfile<<"/PREP7"<<endl; 
    for (int i=0;i<(int)m_Pts.size();i++) 
    { 
        outfile<<"K,"<<i+1<<","<<m_Pts[i].x<<","<<m_Pts[i].y<<","<<m_Pts[i].z<<endl; 
    } 
    for (int i=0;i<(int)m_Edges.size();i++) 
    { 
        outfile<<"L,"<<m_Edges[i].left+1<<","<<m_Edges[i].right+1<<endl; 
    } 
    outfile.close(); 
} 
void GetPoint(double &xx,double &yy,double &zz,string line)//从字符串line中解析出点的x,y,z坐标 
{ 
    int flag=0; 
    string tmp=""; 
    char *cstr; 
    for (int i=(int)line.find(',')+1;i<(int)line.size();i++) 
    { 
        if (line[i]==',') 
        { 
            cstr=new char[tmp.size()+1]; 
            strcpy(cstr,tmp.c_str()); 
            if (flag==0) {xx=atof(cstr);tmp.resize(0);flag++;} 
            else if (flag==1) {yy=atof(cstr);tmp.resize(0);flag++;} 
            continue; 
        } 
        tmp=tmp+line[i]; 
    } 
    if (fabs(xx)<1.0e-6) xx=0.0; 
    if (fabs(yy)<1.0e-6) yy=0.0; 
    if (fabs(zz)<1.0e-6) zz=0.0; 
} 
void Delaunay::Delete_Frame()//删除外边框 
{ 
    EdgeArray BoundEdges; 
    for (int i=0;i<4;i++) m_Pts.erase(m_Pts.begin()); 
    for (int i=0;i<(int)m_Tris.size();i++) 
    { 
        if (m_Tris[i].v[0]==0||m_Tris[i].v[0]==1||m_Tris[i].v[0]==2||m_Tris[i].v[0]==3) 
        { 
            DelTriangle(i,BoundEdges); 
            BoundEdges.resize(0); 
            i--; 
        } 
        else 
        { 
            for (int j=0;j<3;j++) 
            { 
                m_Tris[i].v[j]-=4; 
                m_Tris[i].s[j].left-=4; 
                m_Tris[i].s[j].right-=4; 
            } 
        } 
    } 
    for (int i=0;i<4;i++) m_Edges.erase(m_Edges.begin()); 
    for (int i=0;i<(int)m_Edges.size();i++) 
    { 
        m_Edges[i].left-=4; 
        m_Edges[i].right-=4; 
    } 
}

void Delaunay::Boundary_Recover(int fromPoint,int toPoint)//恢复由指定点组成的边界 
{ 
    EdgeArray BoundEdges; 
    for (int i=0;i<(int)m_Tris.size();i++) 
    { 
        if (m_Tris[i].v[0]>=(fromPoint-1)&&m_Tris[i].v[2]<=(toPoint-1)) 
        { 
            DelTriangle(i,BoundEdges); 
            BoundEdges.resize(0); 
            i--; 
        } 
    } 
}

main.cpp

#include "delaunay.h" 
int main() 
{ 
    ifstream infile("input.txt");//打开"input.txt"文件 
    if (!infile)//判断文件是否正常打开 
    { 
        cout<<"Unable to input nodes!"; 
        exit(1); 
    } 
    string line; 
    PointArray p; 
    double xx,yy,zz; 
    int nodeSize; 
    for (int i=0;i<4;i++)//读入4外边框点 
    { 
        getline(infile,line); 
        GetPoint(xx,yy,zz,line); 
        Point tmp={xx,yy,zz}; 
        p.push_back(tmp); 
    } 
    Delaunay MyMesh(p[0],p[1],p[2],p[3]);//实例化Delaunay类 
    getline(infile,line);//读入节点数，用于后面循环 
    char *cstr; 
    cstr=new char[line.size()+1]; 
    strcpy(cstr,line.c_str()); 
    nodeSize=atoi(cstr); 
    for (int i=0;i<nodeSize;i++)//读入每个节点的坐标 
    { 
        getline(infile,line); 
        GetPoint(xx,yy,zz,line); 
        MyMesh.AddPoint(xx,yy,zz); 
    } 
    infile.close(); 
    MyMesh.Delete_Frame();//删除外边框 
    MyMesh.Boundary_Recover(203,466); 
    MyMesh.Boundary_Recover(467,487); 
    MyMesh.Boundary_Recover(488,511); 
    MyMesh.Boundary_Recover(512,537);//以上都是恢复指定边界 
    MyMesh.output();//将相应ANSYS命令流输出 
    return 0; 
}