Delaunay三角剖分算法（二维）

// DelaunayT.cpp: implementation of the CDelaunayT class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "DelaunayT.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CDelaunayT::CDelaunayT()
{
 nPoints = 0;
 nEdges = 0;
 nTriangles = 0;
}

CDelaunayT::~CDelaunayT()
{

}

/*************************************
功能：对给定新点，首先在TriPoint[]数组中查找一下看是否已经存在，若已存在则
      不做任何改动，若不存在，首先将该新点加入TriPoint[]数组中，
参数：
返回参数：
*************************************/
void CDelaunayT::Triangulation(POINT newPoint)
{
 Vertex addPoint;
 addPoint.x = newPoint.x; 
 addPoint.y = newPoint.y;

    // 判断是否新点已经存在了，若不存在，就加入点集并用插值法
 // 进行Delaunay三角化，存在就不作处理
 if(!TheSamePoint(addPoint))   
 {
  // 将新点插入点集
  nOutLines = 0;
  TriPoint[nPoints].x = newPoint.x;
  TriPoint[nPoints].y = newPoint.y;
  TriPoint[nPoints].index = nPoints;
  nPoints ++;
  
  // 顶点数为2时，创建第一条边
  if(nPoints == 2)
   CreateNewEdge(TriPoint[0], TriPoint[1]);
  
  // 顶点数为3时，首先判断新点是否在第一条边上，若是在就将
  // 第一条边一分为二，若是不在就创造第一个三角形
  else if(nPoints == 3)
  {
   if(TheSameLine(TriEdge[0], TriPoint[2]) == 0)
       CutoTwoLine(TriEdge[0], TriPoint[2]);
   else
       CreateNewTriangle(TriEdge[0], TriPoint[2]);    
  }

  // 当顶点数大于3时
  else if(nPoints >3)
  {
   int i;
   // 首先对所有三角形都遍历一次
   int TriangleNum = nTriangles;
   for(i=0; i<nTriangles; i++)   
   {
    // 判断新加入点是否在被测三角形的外接圆里  
    while((InCircumcircle(Trigon[i], TriPoint[nPoints -1])) && (nTriangles > 0) && i<nTriangles)                                                   
    {  
     // 如果在外接圆里，首先将三角形的三条边的三角形索引值设为-2
                    AmendEdge(Trigon[i], -2);
     
     InOutLines(TriEdge[Trigon[i].BorderA], TriPoint[nPoints -1]);             
     InOutLines(TriEdge[Trigon[i].BorderB], TriPoint[nPoints -1]);
     InOutLines(TriEdge[Trigon[i].BorderC], TriPoint[nPoints -1]);
                   
     //删除这个三角形
     DeleteTriangle(i);    
    }
   }
           
   PutInOrder(OutLine, nOutLines);

   //删除OutLine[]中的边
   bool b;
   for(i=0; i<nOutLines; i++)
   { 
    DeleteEdge(i);
    b = (TriEdge[OutLine[i]].LeftTriangle == -2
     && TriEdge[OutLine[i]].RightTriangle == -2)
     || (TheSameLine(TriEdge[OutLine[i]], TriPoint[nPoints -1]) >= 0
     && TriEdge[OutLine[i]].RightTriangle == -2)
     || (TheSameLine(TriEdge[OutLine[i]], TriPoint[nPoints -1]) <= 0
     && TriEdge[OutLine[i]].LeftTriangle == -2);
    while(b && nOutLines > 1 && (OutLine[i] < nEdges -1))
    {
     nOutLines -- ;
     DeleteEdge(i);
     
     b = (TriEdge[OutLine[i]].LeftTriangle == -2
      && TriEdge[OutLine[i]].RightTriangle == -2)
      || (TheSameLine(TriEdge[OutLine[i]], TriPoint[nPoints -1]) >= 0
      && TriEdge[OutLine[i]].RightTriangle == -2)
      || (TheSameLine(TriEdge[OutLine[i]], TriPoint[nPoints -1]) <= 0
      && TriEdge[OutLine[i]].LeftTriangle == -2);
    }
   }
    

   // 创造新的三角形及其三条边

      int EdgeNum = nEdges;
   for(i=0; i<EdgeNum; i++)
   {
    float side;
    side = TheSameLine(TriEdge[i], TriPoint[nPoints-1]);
       if((side > 0) && (TriEdge[i].LeftTriangle == -1 || TriEdge[i].LeftTriangle == -2)
         || (side < 0) && (TriEdge[i].RightTriangle == -1 || TriEdge[i].RightTriangle == -2))
     CreateNewTriangle(TriEdge[i], TriPoint[nPoints-1]);
   }
  }
 }
}

/*************************************
功能：检测点ap是否在以p1,p2,p3为顶点的三角形的外接圆里，
      在则返回1，否则返回0；
参数：POINT型：p1, p2, p3, ap
返回参数：0/1 
*************************************/
bool CDelaunayT::InCircumcircle(Triangle T, Vertex ap)
{
 Vertex Center, p1, p2, p3;
 float r, d;

 p1 = TriPoint[T.NodeA];
 p2 = TriPoint[T.NodeB];
 p3 = TriPoint[T.NodeC];
 
 float cx, cy;
 float a1, b1, c1, a2, b2, c2;
 
 a1 = (float)(p2.y - p1.y);    b1 = (float)(p1.x - p2.x);
 c1 = (float)((p2.y + p1.y)*(p2.y - p1.y) + (p2.x + p1.x)*(p2.x - p1.x))/(float)2.0;
 
 a2 = (float)(p3.y - p1.y);  b2 = (float)(p1.x - p3.x);
 c2 = (float)((p3.y + p1.y)*(p3.y - p1.y) + (p3.x + p1.x)*(p3.x - p1.x))/(float)2.0;
 
 cx = (a1*c2 - a2*c1)/(b1*a2 - b2*a1);
 cy = (b1*c2 - b2*c1)/(b1*a2 - b2*a1);
   
 Center.x = (int)cx;
 Center.y = (int)cy;
 
 r = Distance(Center, p1);
 d = Distance(Center, ap);

 if(d<=r)
  return 1;
 else return 0; 
}

/********* 求新是否与顶点集中的点重合 ************/
bool CDelaunayT::TheSamePoint(Vertex ap)
{
 int i;
 bool Y = 0;
 for(i=0; i<nPoints; i++)
 {
  if((TriPoint[i].x == ap.x) && (TriPoint[i].y == ap.y))
  {
   Y = 1;
   break;
  }
 }
 return Y;
}

/*求新点是否在被测边上, 返回正数表示点在边的左边，负数表示在边的右边，0表示在边上*/
float CDelaunayT::TheSameLine(Edge e, Vertex ap)
{
 int u1, v1, u2, v2;

 u1 =  TriPoint[e.End].x - TriPoint[e.Start].x;
 v1 =  TriPoint[e.End].y - TriPoint[e.Start].y;
 u2 =  ap.x - TriPoint[e.Start].x;
 v2 =  ap.y - TriPoint[e.Start].y;

    return float(u1 * v2 - v1 * u2);
}

/********* 求两点之间的距离 *********/
float CDelaunayT::Distance(Vertex p1, Vertex p2)
{
 return (float)(sqrt(float((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y))));
}

/***********  创建一条新边  ***********/
void CDelaunayT::CreateNewEdge(Vertex p1, Vertex p2)
{
 if((p2.y > p1.y) || (p1.y == p2.y && p2.x > p1.x))
 {
  TriEdge[nEdges].Start = p1.index;
  TriEdge[nEdges].End = p2.index;
 }
 else
 {
  TriEdge[nEdges].Start = p2.index;
  TriEdge[nEdges].End = p1.index;
 }

 TriEdge[nEdges].LeftTriangle = -1;
 TriEdge[nEdges].RightTriangle = -1;
 TriEdge[nEdges].index = nEdges;
 nEdges ++;
}

/******* 新点在被测边上，将被测边一分为二 *******/
void CDelaunayT::CutoTwoLine(Edge e, Vertex ap)
{
 float v, u, s, maxDis;
 s = Distance(TriPoint[e.Start], TriPoint[e.End]);
 maxDis = s;
 u = Distance(TriPoint[e.Start], ap);
 if(u>maxDis)
  maxDis = u;
 v = Distance(TriPoint[e.End], ap);
 if(v>maxDis)
  maxDis = v;
 if(maxDis == s)
 {  
  DeleteEdge(e.index);
  CreateNewEdge(TriPoint[e.Start], ap);
  CreateNewEdge(TriPoint[e.End], ap);  
 }
 else
 {
  if(u<v)
      CreateNewEdge(TriPoint[e.Start], ap);
        else
         CreateNewEdge(TriPoint[e.End], ap);
 }
}

/*************************************
功能：若被测边不在OutLine[]中，就将其加入其中，若是在，
      则删除被测边，并将最后一条边补足被测边位置。
参数：被测边的索引值
返回参数：无，改变保存在OutLine[]或TriEdge[]数组中。
*************************************/
void CDelaunayT::InOutLines(Edge e, Vertex ap)
{
 int i;
 bool Y=0;
 float side;
 side = TheSameLine(e, ap);
 if((e.LeftTriangle == -2 && e.RightTriangle == -2)
  || (TheSameLine(e, ap) >= 0 && e.RightTriangle == -2)
  || (TheSameLine(e, ap) <= 0 && e.LeftTriangle == -2))
 {
  for(i=0; i<nOutLines; i++)
  {
   if(e.index == OutLine[i])
   {
    Y = 1;
    break;
   }
  }
  if(!Y)
  {
   OutLine[nOutLines] = e.index;
   nOutLines ++;
  }
 }
}

/**** 删除三角形，并将最后的三角形补足被删除三角形位置 ****/
void CDelaunayT::DeleteTriangle(int TrigleIndex)
{
 if((nTriangles == 1) || (TrigleIndex == nTriangles - 1))
  nTriangles --;
 else
 {
  Trigon[TrigleIndex].BorderA = Trigon[nTriangles - 1].BorderA;
  Trigon[TrigleIndex].BorderB = Trigon[nTriangles - 1].BorderB;
  Trigon[TrigleIndex].BorderC = Trigon[nTriangles - 1].BorderC;
  Trigon[TrigleIndex].NodeA = Trigon[nTriangles - 1].NodeA;
  Trigon[TrigleIndex].NodeB = Trigon[nTriangles - 1].NodeB;
  Trigon[TrigleIndex].NodeC = Trigon[nTriangles - 1].NodeC;
  Trigon[TrigleIndex].index = TrigleIndex;

  AmendEdge(Trigon[nTriangles - 1], TrigleIndex); 
  nTriangles --;
 }
 
}

/****** 形成三角形, 并保存三角形以及其新增的两条边 ******/
void CDelaunayT::CreateNewTriangle(Edge &e, Vertex ap)
{
 int e1_index, e2_index;
 e1_index = IsEdgeExist(TriPoint[e.Start], ap);
 e2_index = IsEdgeExist(TriPoint[e.End], ap);
 

 if((e1_index == -1) || (e2_index == -1))
 {
  if(e1_index == -1)
  {
   CreateNewEdge(TriPoint[e.Start], ap);
   e1_index = nEdges-1;
  }
  if(e2_index == -1)
  {
   CreateNewEdge(TriPoint[e.End], ap);
   e2_index = nEdges-1;
  }
 }

 Trigon[nTriangles].NodeA = ap.index;
 Trigon[nTriangles].NodeB = e.Start;
 Trigon[nTriangles].NodeC = e.End;
 Trigon[nTriangles].BorderA = e.index;
 Trigon[nTriangles].BorderB = e2_index;
 Trigon[nTriangles].BorderC = e1_index;
 Trigon[nTriangles].index = nTriangles;
 nTriangles++;
 
 //修正e1, e2两条边的左、右三角形索引值
 AmendEdge(Trigon[nTriangles-1], nTriangles-1);
}

/** 判断以顶点p1, p2为端点的边是否已经存在，Y=1 表示已经存在了**/
int CDelaunayT::IsEdgeExist(Vertex p1, Vertex p2)
{
 int i;
 int pos = -1;
 int xp1, yp1, xp2, yp2;
 for(i=0; i<nEdges; i++)
 {
  xp1 = TriPoint[TriEdge[i].Start].x;
  yp1 = TriPoint[TriEdge[i].Start].y;
  xp2 = TriPoint[TriEdge[i].End].x;
  yp2 = TriPoint[TriEdge[i].End].y;
  if((xp1 == p1.x && yp1 == p1.y && xp2 == p2.x && yp2 == p2.y)
   || (xp1 == p2.x && yp1 == p2.y && xp2 == p1.x && yp2 == p1.y))
  {
   pos = i;
   break;
  }
 }
 return pos;
}

/***** 当删除三角形时，改变三角形三条边的三角形索引值 *****/
void CDelaunayT::AmendEdge(Triangle T, int value)
{
    float As, Bs, Cs;
 As = TheSameLine(TriEdge[T.BorderA], TriPoint[T.NodeA]);
 Bs = TheSameLine(TriEdge[T.BorderB], TriPoint[T.NodeB]);
 Cs = TheSameLine(TriEdge[T.BorderC], TriPoint[T.NodeC]);
 
 if(As > 0)                                                                                               
  TriEdge[T.BorderA].LeftTriangle = value;       
 else TriEdge[T.BorderA].RightTriangle = value;
 
 if(Bs > 0)
  TriEdge[T.BorderB].LeftTriangle = value;
 else TriEdge[T.BorderB].RightTriangle = value;
 
 if(Cs > 0)
  TriEdge[T.BorderC].LeftTriangle = value;
 else TriEdge[T.BorderC].RightTriangle = value;
}

/******* 删除一条边 ********/
void CDelaunayT::DeleteEdge(int Num)
{
 int EdgeIndex = OutLine[Num];
 int LeftT, RightT;
 if((nEdges > 1) && (EdgeIndex < nEdges -1)) //将最后一条边的索引值改为即将删除的边的索引值
 {
  TriEdge[EdgeIndex].Start = TriEdge[nEdges-1].Start;
  TriEdge[EdgeIndex].End = TriEdge[nEdges-1].End;
  TriEdge[EdgeIndex].LeftTriangle = TriEdge[nEdges-1].LeftTriangle;
  TriEdge[EdgeIndex].RightTriangle = TriEdge[nEdges-1].RightTriangle;
  TriEdge[EdgeIndex].index = EdgeIndex;
  
  LeftT = TriEdge[EdgeIndex].LeftTriangle;
  RightT = TriEdge[EdgeIndex].RightTriangle;
  
  if(LeftT>=0) //修改左边三角形信息
  {
   if(Trigon[LeftT].BorderA == nEdges - 1)
    Trigon[LeftT].BorderA = EdgeIndex;
   else if(Trigon[LeftT].BorderB == nEdges - 1)
    Trigon[LeftT].BorderB = EdgeIndex;
   else Trigon[LeftT].BorderC = EdgeIndex;
  }
  
  if(RightT>=0) //修改右边三角形信息
  {
   if(Trigon[RightT].BorderA == nEdges - 1)
    Trigon[RightT].BorderA = EdgeIndex;
   else if(Trigon[RightT].BorderB == nEdges - 1)
    Trigon[RightT].BorderB = EdgeIndex;
   else Trigon[RightT].BorderC = EdgeIndex;
  }
 }

 nEdges--;
}

void CDelaunayT::PutInOrder(int *Group, int groupNum)
{
 int i, j;
 int tmp;
 for(i=0; i<groupNum - 1; i++)
  for(j=i+1; j<groupNum; j++)
   if(Group[i]>Group[j])
   {
    tmp = Group[i];
    Group[i] = Group[j];
    Group[j] = tmp;
   }
}

//DelaunayT.h: interface for the CDelaunayT class.
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_DelaunayT_H__AF5BEF13_CC96_4084_BDFC_E5F40B629356__INCLUDED_)
#define AFX_DelaunayT_H__AF5BEF13_CC96_4084_BDFC_E5F40B629356__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
#define MAXPOINTS  100
#define MAXEDGES  350
#define MAXTRIANGLE  250
#define MAXOUTLINE 25

typedef struct Vertex
{
 int x, y;// 离散点的坐标
 int index;//点的索引
};

typedef struct Triangle
{
 int NodeA; // 三角形的顶点A的坐标索引
 int NodeB; // 三角形的顶点B的坐标索引
 int NodeC; // 三角形的顶点C的坐标索引
 int BorderA; // 三角形的顶点A的对边的索引
 int BorderB; // 三角形的顶点B的对边的索引
 int BorderC; // 三角形的顶点C的对边的索引
 int index; // 三角形的索引
};

typedef struct Edge
{
 int Start; // 边的起点的索引
 int End; // 边的终点的索引
 int LeftTriangle; // 边的左三角形索引
 int RightTriangle; // 边的右三角形索引
 int index; // 边的索引
};

class CDelaunayT 
{
public:
 CDelaunayT();
 virtual ~CDelaunayT();

public:
 int nPoints;  // 顶点数
 int nEdges;  // 边数
 int nTriangles;  // 三角形数
 int nOutLines;  // 需修改的边数
 int OutLine[MAXOUTLINE];  // 保存需修改的边
 Vertex TriPoint[MAXPOINTS];  // 保存所有顶点的点结构数组
 Edge TriEdge[MAXEDGES];  // 保存所有边的边结构数组
 Triangle Trigon[MAXTRIANGLE];  // 保存所有三角形的三角形结构数组
   
public:
 void PutInOrder(int *Group, int groupNum);
 // 中心函数
 void Triangulation(POINT newPoint);  // 三角化

 // 操作函数
 bool TheSamePoint(Vertex ap);  // 判断一个点是否已经存在
    int IsEdgeExist(Vertex p1, Vertex p2);  // 判断一条边是否已经存在
    void CutoTwoLine(Edge e, Vertex ap);  // 将一条边一分为二
 void CreateNewEdge(Vertex p1, Vertex p2);  // 创造一条新边
 void CreateNewTriangle(Edge &e, Vertex ap );  // 创造一个新三角形
 void InOutLines(Edge e, Vertex ap);  // 查找某条边是否已在OutLine[]数组里
 void DeleteEdge(int Num);   // 删除一条边
 void DeleteTriangle(int TrigleIndex);  // 删除一个三角形
 void AmendEdge(Triangle T, int value);   // 修复一条边
 // 数学函数
 float Distance(Vertex p1, Vertex p2);  // 求两点之间的距离
 float TheSameLine(Edge e, Vertex ap);  // 判断一个点与一条直线的关系
 bool InCircumcircle(Triangle T, Vertex ap);  // 判断给定点是否在给定三角形内

};

#endif // !defined(AFX_DelaunayT_H__AF5BEF13_CC96_4084_BDFC_E5F40B629356__INCLUDED_)