三角剖分.cpp

平面三角剖分,可以在DOS下显示出三角形顶点的集,数据不是很完善,可以做进一步工作的基础

//版权所有 Rungo(az-zxk) 学号s060269

#include "az_zxk.h"
using namespace std; 

////////////////////////////////////////////////////////////////////////
// CircumCircle() :画外接圆并当点(xp,yp)在由(x1,y1), (x2,y2), (x3,y3)
//                 三点构成的外接圆里时返回true
//      (xc,yc)为外接圆的中心，r为外接圆半径
// Note : 点在边界上时看成在外接圆里
//                                          Interpreter:az-zxk (s060269)
////////////////////////////////////////////////////////////////////////

int CircumCircle(double xp, double yp, double x1, double y1, double x2, 
double y2, double x3, double y3, double &xc, double &yc, double &r)
{
  double m1, m2, mx1, mx2, my1, my2;
  double dx, dy, rsqr, drsqr;

/* 检查符合精度要求的点 */
if(abs(y1 - y2) < EPSILON && abs(y2 - y3) < EPSILON)
  return(false);
if(abs(y2-y1) < EPSILON)//1,2点近似看成在一条直线上
{ 
  m2 = - (x3 - x2) / (y3 - y2);
  mx2 = (x2 + x3) / 2.0;
  my2 = (y2 + y3) / 2.0;
  xc = (x2 + x1) / 2.0;
  yc = m2 * (xc - mx2) + my2;
}
    else if(abs(y3 - y2) < EPSILON)//3,2点近似看成在一条直线上
    {
        m1 = - (x2 - x1) / (y2 - y1);
        mx1 = (x1 + x2) / 2.0;
        my1 = (y1 + y2) / 2.0;
        xc = (x3 + x2) / 2.0;
        yc = m1 * (xc - mx1) + my1;
	}
else
{        m1 = - (x2 - x1) / (y2 - y1); //求斜率
         m2 = - (x3 - x2) / (y3 - y2); //求斜率
         mx1 = (x1 + x2) / 2.0; //边中心
         mx2 = (x2 + x3) / 2.0;
         my1 = (y1 + y2) / 2.0;
         my2 = (y2 + y3) / 2.0;
         xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2); 
         yc = m1 * (xc - mx1) + my1; 
}
  dx = x2 - xc;
  dy = y2 - yc;
  rsqr = dx * dx + dy * dy;
  r = sqrt(rsqr); 
  dx = xp - xc;
  dy = yp - yc;
  drsqr = dx * dx + dy * dy;
  return((drsqr <= rsqr) ? true : false);
}
///////////////////////////////////////////////////////////////////////////////
// Triangulate() :构建三角形其返回植是ntri个三角形网的链表(数组v)
//   这些三角形在数组中按顺时间方向排列
// Note :函数实参已经经过qsort(p,nv,sizeof(XYZ),XYZCompare);按x坐标大小排列
//                                                 Interpreter:az-zxk (s060269)
///////////////////////////////////////////////////////////////////////////////

int Triangulate(int nv, XYZ pxyz[], ITRIANGLE v[], int &ntri)
{
  int *complete = NULL;
  IEDGE *edges = NULL; 
  IEDGE *p_EdgeTemp;
  int nedge = 0;
  int trimax, emax = 200;
  int status = 0;
  int inside;
  int i, j, k;
  double xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r;
  double xmin, xmax, ymin, ymax, xmid, ymid;
  double dx, dy, dmax; 

/* 给三角形的数组分配空间 */
  trimax = 4 * nv;
  complete = new int[trimax];
/* 给边界的数组分配空间 */
  edges = new IEDGE[emax];
/* 将点数组分成小集*/
  xmin = pxyz[0].x;
  ymin = pxyz[0].y;
  xmax = xmin;
  ymax = ymin;
  for(i = 1; i < nv; i++)
  {
    if (pxyz[i].x < xmin) xmin = pxyz[i].x;
    if (pxyz[i].x > xmax) xmax = pxyz[i].x;
    if (pxyz[i].y < ymin) ymin = pxyz[i].y;
    if (pxyz[i].y > ymax) ymax = pxyz[i].y;
  }
  dx = xmax - xmin;
  dy = ymax - ymin;
  dmax = (dx > dy) ? dx : dy;
  xmid = (xmax + xmin) / 2.0;
  ymid = (ymax + ymin) / 2.0;
/*
   构建第一个三角形，将其添加到定点数组的最后，
   这样它就成为三角形数组的第一个三角形
*/
  pxyz[nv+0].x = xmid - 20 * dmax;
  pxyz[nv+0].y = ymid - dmax;
  pxyz[nv+1].x = xmid;
  pxyz[nv+1].y = ymid + 20 * dmax;
  pxyz[nv+2].x = xmid + 20 * dmax;
  pxyz[nv+2].y = ymid - dmax;
  v[0].p1 = nv;
  v[0].p2 = nv+1;
  v[0].p3 = nv+2;//第一个三角形的三个顶点
  complete[0] = false;
  ntri = 1;
/*保证每个点在现有三角网里只算一次*/
 for(i = 0; i < nv; i++)
 {
    xp = pxyz[i].x;
    yp = pxyz[i].y;//结点p的坐标
    nedge = 0;
/*建立边界“缓冲器”，当结点p在外接圆内，则将这个三角形的三边加入边界“缓冲器”，
并且将这个三角形从数组切掉 */
  for(j = 0; j < ntri; j++)
  {
    if(complete[j])
      continue;
    x1 = pxyz[v[j].p1].x;
    y1 = pxyz[v[j].p1].y;
    x2 = pxyz[v[j].p2].x;
    y2 = pxyz[v[j].p2].y;
    x3 = pxyz[v[j].p3].x;
    y3 = pxyz[v[j].p3].y;
    inside = CircumCircle(xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r);
    if (xc + r < xp)// if (xc + r + EPSILON < xp)
      complete[j] = true;
    if(inside)/* 检查保证未超过边界数组的最大容量 */
	{
      if(nedge + 3 >= emax)
	  {
        emax += 100;
        p_EdgeTemp = new IEDGE[emax];
        for (int i = 0; i < nv; i++)
		{
          p_EdgeTemp[i] = edges[i];   
        }
        delete []edges;
        edges = p_EdgeTemp;
      }
      edges[nedge+0].p1 = v[j].p1;
      edges[nedge+0].p2 = v[j].p2;
      edges[nedge+1].p1 = v[j].p2;
      edges[nedge+1].p2 = v[j].p3;
      edges[nedge+2].p1 = v[j].p3;
      edges[nedge+2].p2 = v[j].p1;
      nedge += 3;
      v[j] = v[ntri-1];
      complete[j] = complete[ntri-1];
      ntri--;
      j--;
	}
  }
/* 如果三角形都是按逆时针放入数组的，则内边界元素被相反的放置*/
  for(j = 0; j < nedge - 1; j++)
  {
    for(k = j + 1; k < nedge; k++)
	{
      if((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1))
	  {
        edges[j].p1 = -1;
        edges[j].p2 = -1;
        edges[k].p1 = -1;
        edges[k].p2 = -1;
      }
      if((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2))
	  {
        edges[j].p1 = -1;
        edges[j].p2 = -1;
        edges[k].p1 = -1;
        edges[k].p2 = -1;
      }
    }
  }
/*为当前点建立新的三角形，所有的边按顺时间方向排列*/
  for(j = 0; j < nedge; j++) 
  {
    if(edges[j].p1 < 0 || edges[j].p2 < 0)
      continue;
    v[ntri].p1 = edges[j].p1;
    v[ntri].p2 = edges[j].p2;
    v[ntri].p3 = i;
    complete[ntri] = false;
    ntri++;
  }
 }

 for(i = 0; i < ntri; i++) 
  {
    if(v[i].p1 >= nv || v[i].p2 >= nv || v[i].p3 >= nv)
	{
      v[i] = v[ntri-1];
      ntri--;
      i--;
    }
  }
  delete[] edges;
  delete[] complete;
  return 0;
} 

/*void randomize()
{
  srand((time_t) time(NULL));  
}

int random(int n)
{
  return rand()%n; 
}
*/
int XYZCompare(const void *v1, const void *v2)
{
  XYZ *p1, *p2;
    
  p1 = (XYZ*)v1;
  p2 = (XYZ*)v2;
  if(p1->x < p2->x)
    return(-1);
  else if(p1->x > p2->x)
         return(1);
       else
         return(0);
}

void outputtriangle(int &nv, XYZ p[], ITRIANGLE v[], int &ntri)
{
  int X, Y, i = 0;
  int max = 10;
  double x, y;
  for(i = 0; i < ntri; i++)
  {
    cout<<(int)p[v[i].p1].x<<" "<< (int)p[v[i].p1].y<<" "<< (int)p[v[i].p2].x
    <<" "<< (int)p[v[i].p2].y<<"\n";
    cout<<(int)p[v[i].p2].x<<" "<< (int)p[v[i].p2].y<<" "<< (int)p[v[i].p3].x
    <<" "<< (int)p[v[i].p3].y<<"\n";
    cout<<(int)p[v[i].p3].x<<" "<< (int)p[v[i].p3].y<<" "<< (int)p[v[i].p1].x
    <<" "<< (int)p[v[i].p1].y<<"\n";
	cout<<"\n";
  }
}

int main()
{
  ITRIANGLE *v = NULL;
  int max = 10 , n = 1;
  XYZ *p = new XYZ[max]; 
  XYZ *p_Temp = NULL;   
  int nv = 0;
  int X, Y;
  int i;
  int ntri = 0;
  double x, y;
  bool b_Ok = false;
  //randomize();
  nv = 0;
  p = new XYZ[max];
  cout<<"请输入点坐标:(测试共"<<n_MaxPoints<<"点)"<<"\n";//测试点获得
  while (nv != n_MaxPoints)
  {
    do
	{ 
	  b_Ok = true;
      cout <<"第"<<n<<"点坐标"<<"\n";	  
	  cout <<"x=";
	  cin >> x;
	  //x = (double)random(500);//randomize()和random(int n)只是附加测试用的
	  cout <<"y=";
      cin >> y;
	  //y = (double)random(500);
	  n++;
	  for(int n_Cpt = 0; n_Cpt <= nv; n_Cpt++)
	  {
        if((x == p[n_Cpt].x) && (y == p[n_Cpt].y)) b_Ok = false;
      }// 避免点重复
	}while(!b_Ok);
    if (nv >= max)
	{
      max = max * 2;            
      p_Temp = new XYZ[max]; 
      for (int i = 0; i < nv; i++) 
	  {
        p_Temp[i] = p[i];  
      }
      delete []p;  
      p = p_Temp; 
	}   
    p[nv].x = x * 1.0;
    p[nv].y = y * 1.0;
    nv++;
  }
  p_Temp = new XYZ[nv + 3]; 
  for (i = 0; i < nv; i++) 
  {
    p_Temp[i] = p[i];      
  }
  delete []p;           
  p = p_Temp;
  v = new ITRIANGLE[3 * nv]; 
  qsort(p, nv, sizeof(XYZ), XYZCompare); //将点数组按x大小排列
  Triangulate(nv, p, v, ntri);//构建三角网
  outputtriangle(nv, p, v, ntri);  //输出
  delete []p;//空间释放
  delete []v;
  p = NULL;
  v = NULL;
  system("pause"); 
  return 0;
}