Delaunay三角剖分—三角网生长算法

Delaunay三角网是俄国数学家B.Delaunay于1934年发现的。关于Delaunay三角网构建的研究有许多，但由于本课题具有数据量大的特征，不宜直接沿用已有构建方法，笔者针对本课题数据特征，研究获得了适应本课题，速度较快的构建方法。Delaunay三角网有一个特性，每个三角网形成的外接圆都不包含其他参考点。利用这一个性质，我们可以直接构成Delaunay三角网：
一、建立第一个三角形
1、判断用来建立TIN的总脚点数，小于3则报错退出。
2、第一点的选择：
链表的第一节点，命名为Pt1;
3、第二点的选择：
A.非Pt1点； B.距Pt1最近
命名为Pt2
4、第三点的选择
A.非Pt1，Pt2点；
B.与Pt1，Pt2点组成的三角形的外接圆内无其他节点；
C.与Pt1，Pt2组成的三角形中的角Pt1Pt3Pt2最大。
命名为Pt3
5、生成三边，加入边表
6、生成第一个三角形，组建三角形表
二、扩展TIN
1、从边表头取一边，要求：该边flag标志为假（只在一个三角形中）
2、从点链表中搜索一点，要求：
A、与边中的Pixel3在边的异侧；
B、该点与边组成的三角形的外接圆内无其他点
C、满足上面两条件的点中角Pt1Pt3Pt2最大的点为Pt3。
3、判断新生成的边，如果边表中没有，则加入边表尾，设定标志flag为假，如果有，则设定该边flag为真。
4、将生成的三角形假如三角形表
5、设定选中的边的标志flag为真
6、转至1，直至边表边的标志flag全部为真。

注：如果需要进行限制三角剖分，则可利用重心法取出不要的三角形，必要时，可对边界进行限制，不让生成的边与边界相交

数据结构：

struct Pixel    //脚点数据
{
double x,y,z,g;
bool flag;
};
struct List //数据链表
{
Pixel *pixel;
List *next;
};
struct Line //三角形边
{
Pixel *pixel1;  //三角形边一端点
Pixel *pixel2;  //三角形边另一端点
Pixel *pixel3;  //三角形边所对顶点
bool flag;
};
struct Linelist //三角形边表
{
Line *line;
Linelist *next;
};
struct Triangle //三角形表
{
Line *line1;
Line *line2;
Line *line3;
Triangle *next;
};



以下是程序中关于建网的部分：
// DelaunayTIN.cpp: implementation of the CDelaunayTIN class.
//
//
//功能：   用给定的数据链表数据，组建Delaunay不规则三角网
//输入参数：数据链表list;区域范围(XMin,YMin),(XMax,YMax)
//输出参数：不规则三角网首三角形地址

 Triangle * CSimRegular::CreateDelaunayTIN(List *list)
{
    //组建第一个三角形
     List *node;
     Pixel *pt1,*pt2,*pt3;
    bool flag;
    Triangle *TIN;
    pt1=list->pixel;
    pt2=list->next->pixel;
    node=list->next->next;
    while(node!=NULL)
    {
        if(
         (pt1->x-node->pixel->x)*(pt1->x-node->pixel->x)+(pt1->y-node->pixel->y)*(pt1->y-node->pixel->y)
       <(pt1->x-pt2->x)*(pt1->x-pt2->x)+(pt1->y-pt2->y)*(pt1->y-pt2->y)         
            )
        {
            pt2=node->pixel;
        }
        node=node->next;
    }
    node=list->next;
    pt3=NULL;
    while(node!=NULL)
    {
        if(node->pixel==pt1 || node->pixel==pt2)
        {
            node=node->next;
            continue;
        }

        if(pt3==NULL)
        {
            pt3=node->pixel;
        }
        else
        {
            float dist11=sqrt((pt1->x-node->pixel->x)*(pt1->x-node->pixel->x)+(pt1->y-node->pixel->y)*(pt1->y-node->pixel->y));
            float dist12=sqrt((pt2->x-node->pixel->x)*(pt2->x-node->pixel->x)+(pt2->y-node->pixel->y)*(pt2->y-node->pixel->y));
            float dist12_3=sqrt((pt1->x-pt2->x)*(pt1->x-pt2->x)+(pt1->y-pt2->y)*(pt1->y-pt2->y));
            float dist21=sqrt((pt1->x-pt3->x)*(pt1->x-pt3->x)+(pt1->y-pt3->y)*(pt1->y-pt3->y));
            float dist22=sqrt((pt3->x-pt2->x)*(pt3->x-pt2->x)+(pt3->y-pt2->y)*(pt3->y-pt2->y));
            if((pow(dist11,2)+pow(dist12,2)-pow(dist12_3,2))/(2*dist11*dist12)
                <(pow(dist21,2)+pow(dist22,2)-pow(dist12_3,2))/(2*dist21*dist22)) //余弦判断角度

            {
                pt3=node->pixel;
            }
        }
        node=node->next;
    }
    //LineList
    Linelist *linehead,*linenode,*linelast;
    Line *ln1,*ln2,*ln3;
    linenode=new Linelist;
    linenode->line=new Line;
    linenode->line->pixel1=pt1;
    linenode->line->pixel2=pt2;
    linenode->line->pixel3=pt3;
    linenode->line->flag=false;
    linenode->next=NULL;
    linehead=linelast=linenode;
    ln1=linenode->line;
    linenode=new Linelist;
    linenode->line=new Line;
    linenode->line->pixel1=pt2;
    linenode->line->pixel2=pt3;
    linenode->line->pixel3=pt1;
    linenode->line->flag=false;
    linenode->next=NULL;
    linelast->next=linenode;
    linelast=linenode;
    ln2=linenode->line;
    linenode=new Linelist;
    linenode->line=new Line;
    linenode->line->pixel1=pt3;
    linenode->line->pixel2=pt1;
    linenode->line->pixel3=pt2;
    linenode->line->flag=false;
    linenode->next=NULL;
    linelast->next=linenode;
    linelast=linenode;
    ln3=linenode->line;
    //first Triangle
    Triangle *tglhead,*tglnode,*tgllast;
    tglnode=new Triangle;
    tglnode->line1=ln1;
    tglnode->line2=ln2;
    tglnode->line3=ln3;
    tglnode->next=NULL;
    tglhead=tgllast=tglnode;


    //expend tin;
    Linelist *linetmp,*linetemp;
    List *pixeltmp;
    float x1,y1,x2,y2,x3,y3;
    linetmp=linehead;
    while(linetmp!=NULL)
    {
        if(linetmp->line->flag==true)
        {
            linetmp=linetmp->next;
            continue;
        }
        ln1=linetmp->line;
        pt1=linetmp->line->pixel1;
        pt2=linetmp->line->pixel2;
        x1=linetmp->line->pixel1->x;
        y1=linetmp->line->pixel1->y;
        x2=linetmp->line->pixel2->x;
        y2=linetmp->line->pixel2->y;
        x3=linetmp->line->pixel3->x;
        y3=linetmp->line->pixel3->y;


        pixeltmp=list;
        pt3=NULL;
        while(pixeltmp!=NULL)
        {
            if(pixeltmp->pixel==pt1 || pixeltmp->pixel==pt2)
            {
                pixeltmp=pixeltmp->next;
                continue;
            }
            if(((y2-y1)*pixeltmp->pixel->x+(x1-x2)*pixeltmp->pixel->y+(x2*y1-x1*y2))
            *((y2-y1)*x3+(x1-x2)*y3+(x2*y1-x1*y2))>=0)
            {
                pixeltmp=pixeltmp->next;
                continue;
            }


            if(pt3==NULL)pt3=pixeltmp->pixel;
            else
            {
                float dist11=sqrt((pt1->x-pixeltmp->pixel->x)*(pt1->x-pixeltmp->pixel->x)+(pt1->y-pixeltmp->pixel->y)*(pt1->y-pixeltmp->pixel->y));
                float dist12=sqrt((pt2->x-pixeltmp->pixel->x)*(pt2->x-pixeltmp->pixel->x)+(pt2->y-pixeltmp->pixel->y)*(pt2->y-pixeltmp->pixel->y));
                float dist12_3=sqrt((pt1->x-pt2->x)*(pt1->x-pt2->x)+(pt1->y-pt2->y)*(pt1->y-pt2->y));
                float dist21=sqrt((pt1->x-pt3->x)*(pt1->x-pt3->x)+(pt1->y-pt3->y)*(pt1->y-pt3->y));
                float dist22=sqrt((pt3->x-pt2->x)*(pt3->x-pt2->x)+(pt3->y-pt2->y)*(pt3->y-pt2->y));
                if((pow(dist11,2)+pow(dist12,2)-pow(dist12_3,2))/(2*dist11*dist12)
                    <(pow(dist21,2)+pow(dist22,2)-pow(dist12_3,2))/(2*dist21*dist22)) //余弦判断角度
                {
                    pt3=pixeltmp->pixel;
                }
            }
            pixeltmp=pixeltmp->next;
        }
        if(pt3!=NULL)
        {
            linetemp=linehead;
            flag=false;
            while(linetemp!=NULL)
            {
                if((pt1==linetemp->line->pixel1 && pt3==linetemp->line->pixel2)
                    || (pt3==linetemp->line->pixel1 && pt1==linetemp->line->pixel2))
                {
                    linetemp->line->flag=true;
                    flag=true;
                    ln2=linetemp->line;
                    break;
                }
                linetemp=linetemp->next;
            }
            if(!flag)
            {
                linenode=new Linelist;
                linenode->line=new Line;
                linenode->line->pixel1=pt3;
                linenode->line->pixel2=pt1;
                linenode->line->pixel3=pt2;
                linenode->line->flag=false;
                linenode->next=NULL;
                linelast->next=linenode;
                linelast=linenode;
                ln2=linenode->line;
            }
            linetemp=linehead;
            flag=false;
            while(linetemp!=NULL)
            {
                if((pt2==linetemp->line->pixel1 && pt3==linetemp->line->pixel2)
                    || (pt3==linetemp->line->pixel1 && pt2==linetemp->line->pixel2))
                {
                    linetemp->line->flag=true;
                    flag=true;
                    ln3=linetemp->line;
                    break;
                }
                linetemp=linetemp->next;
            }
            if(!flag)
            {
                linenode=new Linelist;
                linenode->line=new Line;
                linenode->line->pixel1=pt2;
                linenode->line->pixel2=pt3;
                linenode->line->pixel3=pt1;
                linenode->line->flag=false;
                linenode->next=NULL;
                linelast->next=linenode;
                linelast=linenode;
                ln3=linenode->line;
            }
            tglnode=new Triangle;
            tglnode->line1=ln1;
            tglnode->line2=ln2;
            tglnode->line3=ln3;
            tglnode->next=NULL;
            tgllast->next=tglnode;
            tgllast=tglnode;
        }


        linetmp->line->flag=true;
        linetmp=linetmp->next;
    }
    TIN=tglhead;
       return TIN;