cdt.cpp

这是个生成约束delaunay triangulations源程序。对学习约束delaunay triangulations很有帮助！

// file name: cdt.cpp
//
//      main program to compute a
//      constrained delaunay triangultaion (CDT).
//
// author: Zhihua Wang ([e-mail])
// date: Dec 12, 2002

/*************************************************************


    Purpose:
      This program implements a divide-and-conquer algorithm
      to compute the Constrained Delaunay Triangulation (CDT) for a given 
     planar straight line graph.
      For detail, please refer to the accompanying paper "cdt.ps"
     about this algorithm.

    Usage:
          % ./cdt input_file [output_file]

      where the format of the input_file is described
      in the accompanying doc/format.txt in this directory.
      Basically, the inputFile format is a subset of .noff format.

      If the output_file is not given, the default output file will
   be called "result.noff".  Output file format is the same
   as the input file format.


   More Details:
   The input file comprises a sequence of n vertices and 
      m line-segments in the plane.

        The output file describes a triangulation represented 
      by the same sequence of n vertices, and a sequence
      of M line-segments (M >= m).  All the input line-
      segments will be included in the output.

    Algorithm:

    (1) Initialization:
        1) Read and store input as a half-edge data structure;
            the necessary topological information will be computed
       (e.g., the clockwise order of all the edges incident
      upon a given vertex).

           1.1) Read the vertices and save them in an array.
           1.2) Sort vertices by X-coordinate and construct a 
      map into the original array of vertices.
           1.3) Read the edges and construct the 
                 half-edges using the map above.

        2) Conceptually, add n+1 vertical lines which separate each vertex.  
           (We assume that no two vertices have same X-value.)

        3) Add two horizontal line-segments which lie (respectively) above
          and below all the vertices; the line-segments
     is sufficiently long so that it crosses all the vertical lines.

    (2) Division:
           Create a list of quads, where each quad is the bounded region 
           between two consecutive vertical lines and between the
      two horizontal line-segments.


    (3) Conquer:
   We pair up the quads so that each pair are adjacent to each
   other.  For each pair of quads, we want to merge them 
        and construct the CDT of the merged quad.

    NOTE:
      The function with a comment line "//#EXACT " before it 
      means it requires exact computation.

*************************************************************/

#include "cdt.h"
#include <math.h>
#include <string.h>
#include <GL/glut.h>

//*****************************************************************
//extension of IO for Debuging
//*****************************************************************
//print information about slab
void   printSlab(SPointer s_head)
{

      printf("%s\n","# information of slab and quads");
//      cout<<"# information of slab and quads"<<endl;
      SPointer   tmps;
      tmps=s_head;
      while   (tmps!=NULL)
      {
//         cout<<""# slab ("<<tmps->l<<','<<tmps->r<<")"<<endl;
         printf("%s%d,%d%s\n","# slab (",tmps->l,tmps->r,")");

         QPointer   tmpq;
         int   cnt=0;
         tmpq=tmps->q_head;
         while   (tmpq!=NULL)
         {
            cnt++;
            printf("%s%d%s%d\n","# quad (",cnt,") highest vertex's index is ",tmpq->top_v->getIndex());
            printf("%s%d%s%d\n","# quad (",cnt,") lowest vertex's index is ",tmpq->btm_v->getIndex());
//            printf("%s%d%s%f\n","# quad (",cnt,") top y =",tmpq->topy);
//            printf("%s%d%s%f\n","# quad (",cnt,") bottom y =",tmpq->btmy);
            printf("%s%d%s%d%s%d\n","# quad (",cnt,") top half edge is ",tmpq->top_he->getBackVertex()->getIndex(),"-->",tmpq->top_he->getFrontVertex()->getIndex());
            printf("%s%d%s%d%s%d\n","# quad (",cnt,") bottom half edge is ",tmpq->btm_he->getBackVertex()->getIndex(),"-->",tmpq->btm_he->getFrontVertex()->getIndex());
/*            cout<<"# quad ("<<cnt<<") highest vertex's index is "<<tmpq->hv->index<<endl;
            cout<<"# quad ("<<cnt<<") lowest vertex's index is "<<tmpq->lv->index<<endl;
//            printf("%s%d%s%f\n","# quad (",cnt,") top y =",tmpq->topy);
//            printf("%s%d%s%f\n","# quad (",cnt,") bottom y =",tmpq->btmy);
            printf("%s%d%s%d%s%d\n","# quad (",cnt,") top half edge is ",tmpq->top->getBackVertex()->index,"-->",tmpq->top->fv->index);
            cout<<"# quad ("<<cnt<<") top half edge is "<<tmpq->top->getBackVertex()->index<<"-->"<<tmpq->top->fv->index<<endl;
            printf("%s%d%s%d%s%d\n","# quad (",cnt,") bottom half edge is ",tmpq->btm->getBackVertex()->index,"-->",tmpq->btm->fv->index);
            cout<<"# quad ("<<cnt<<") bottom half edge is "<<tmpq->btm->getBackVertex()->index<<"-->"<<tmpq->btm->fv->index<<endl;
*/
         // forward one link
            tmpq=tmpq->next;
         }
         tmps=tmps->next;   
      }   
}

//*****************************************************************
//extension of geomtry primitive operations 
//*****************************************************************
/*
  The side of v to he
*/
int   vertexToEdge(VPointer v,HePointer he)
{
   Point   q=v->pos;
   Point   p1=he->getBackVertex()->pos;
   Point   p2=he->getFrontVertex()->pos;

   return pointToLine(q.x,q.y,p1.x,p1.y,p2.x,p2.y);
}

/*
  The side of v to edge (v1->v2)
*/
int   vertexToVertexVertex(VPointer v,VPointer v1,VPointer v2)
{
   Point   q=v->pos;
   Point   p1=v1->pos;
   Point   p2=v2->pos;

   return pointToLine(q.x,q.y,p1.x,p1.y,p2.x,p2.y);
}

/*
  The side of v to circumcircle of v1,v2,v3
*/
int   vertexToCircle(VPointer v,VPointer v1,VPointer v2,VPointer v3)
{
   Point   q=v->pos;
   Point   p1=v1->pos;
   Point   p2=v2->pos;
   Point   p3=v3->pos;
#ifdef   DETECT_COLLINEAR_DEGENERACY
   if   (pointToLine(p1.x,p1.y,p2.x,p2.y,p3.x,p3.y)==0)
      reportError("Three collinear points can no define a circumcircle.",__LINE__,__FILE__);
#endif
   return   pointToCircle(q.x,q.y,p1.x,p1.y,p2.x,p2.y,p3.x,p3.y);
}

/*
  test on whether edge he lies in angle from he1 to he2 
*/
int    edgeInAngle(HePointer he,HePointer he1,HePointer he2)
{
/*
#ifdef   DEBUG
   if   (!   (he->getBackVertex()==he1->getBackVertex() && he->getBackVertex()==he2->getBackVertex()) )
   {
      reportError("Wrong use of call to function edgeInAngle().",__LINE__,__FILE__);
   }
#endif
   */
//   return   pointInAngle(he->getBackVertex()->x,he->getBackVertex()->y,he->fv->x,he->fv->y,he1->fv->x,he1->fv->y,he2->fv->x,he2->fv->y);
   Point   o=he1->getBackVertex()->pos;
   Point   q=he->getFrontVertex()->pos;
   Point   p1=he1->getFrontVertex()->pos;
   Point   p2=he2->getFrontVertex()->pos;

   return   pointInAngle(o.x,o.y,q.x,q.y,p1.x,p1.y,p2.x,p2.y);
}

/* 
  test whether vertex v lies in the angle from (vo->v1) to (vo->v2)
*/
int   vertexInAngle(VPointer vo,VPointer v,VPointer v1,VPointer v2)
{
   Point   o=vo->pos;
   Point   q=v->pos;
   Point   p1=v1->pos;
   Point   p2=v2->pos;

   return   pointInAngle(o.x,o.y,q.x,q.y,p1.x,p1.y,p2.x,p2.y);
}


/*
// get the halfedge from v1 to v2
HPointer getHalfedge(VPointer v1,VPointer v2)
{
   HPointer tmpe;
   tmpe=v1->out;
   if   (tmpe==NULL)
   {
      return   NULL;
   }
   else
   {
      if   (tmpe->fv ==v2)   return   tmpe;
      tmpe=tmpe->clk;
      while   (tmpe->fv !=v2 &&   tmpe!=v1->out)
      {
         tmpe=tmpe->clk;
      }
      if   (tmpe==v1->out)
      {   return NULL;}
      else
      {   return   tmpe;}
   }
}
*/
//*****************************************************************
//implementation of primitives
//*****************************************************************
Vertex::Vertex(double _x,double _y,int   _i)
:index(_i)
{
   pos.x=_x;
   pos.y=_y;
   out_he=NULL;
   valence=0;
}
/*
Vertex::Vertex(const Point &p,int _i)
:pos(p),index(_i)
{
   out_he=NULL;
   valence=0;
}
*/
Vertex::~Vertex()
{
   /*
   while   (out_he!=NULL)
   {
      delete   out_he->e;
   }
   */
}

Halfedge::~Halfedge()
{
   //getBackVertex()->deleteFromOutHalfedgeList(this);
}


// utility to maintain the halfedge ring of vertex
// assume he is not added yet
void Vertex::addToOutHalfedgeList(HePointer he)
{
   if   (out_he==NULL)
   {
      he->next_he=he->pre_he=he;
      out_he=he;
      return;
   }

   if   (out_he->next_he==out_he)
   {
      out_he->next_he=he;
      out_he->pre_he=he;
      he->next_he=out_he;
      he->pre_he=out_he;
      return;
   }

   // insert in order of ring_circulation_order
   HePointer n_he;

   n_he=out_he->next_he;
   int flag;
   flag=edgeInAngle(he,n_he,n_he->pre_he);

   // he is in the angle span from n_he->pre_he to n_he
   while   (flag==0)
      //&& n_he!=out_he)   //we can always find it
   {
      //forward
      n_he=n_he->next_he;
      flag=edgeInAngle(he,n_he,n_he->pre_he);
   }

   if (flag!=1)
   {
      reportError("Collinear Degeneracy!",__LINE__,__FILE__);
   }
   //else we find inserting place, in fact we can always find it
   //INSERT_INBETWEEN(he,n_he->pre_he,n_he);
   he->next_he=n_he;
   he->pre_he=n_he->pre_he;
   n_he->pre_he->next_he=he;
   n_he->pre_he=he;
}

// not an interface for user
// assume he's back vertex is this v, v->out_he is not NULL
void Vertex::deleteFromOutHalfedgeList(HePointer he)
{
/*
#ifdef DEBUG
   assert(he);
   assert(he->getBackVertex()==this);   // this means this v has at least one out halfedge, so out_he!=NULL
#endif//DEBUG
*/
   if ( out_he == out_he->next_he )
   {out_he = NULL;}
    else if ( he == out_he )
      out_he = out_he->next_he;      // move out_he to next place
   
   //delete he from the list
   he->next_he->pre_he = he->pre_he;
   he->pre_he->next_he = he->next_he;
//   delete he;
}

// find the halfedge which point from this v to v2
HePointer Vertex::findOutHalfedgeInRing(VPointer v2)
{
#ifdef DEBUG
   assert(v2);
#endif//DEBUG
   if   (out_he==NULL)   return   NULL;
   HePointer   he=out_he->pre_he;
   while   (he!=out_he   &&   he->v!=v2)
      he=he->pre_he;
   if   (he->v==v2)   return   he;
   return   NULL;
}

// find the halfedge which point from v2 to this v 
HePointer Vertex::findInHalfedgeInRing(VPointer v2)
{
   HePointer he=findOutHalfedgeInRing(v2);
   if   (he!=NULL)   return   he->inv_he;
   return   NULL;
}

// not an interface to user, only called by new Edge()
// add halfedge and update all the topology information
// assume v->v2 is not added
HePointer   Vertex::addHalfedge(VPointer v2,EPointer e)
{
   
   HePointer   he=new   Halfedge(v2,e);

   Vector   d;
   d.x=v2->pos.x-this->pos.x;
   d.y=v2->pos.y-this->pos.y;

   he->d=d;

   valence++;
   //this will update next_he and pre_he of he
   addToOutHalfedgeList(he);

   //inv_he will be updated by new Edge()
   return he;
}

Edge::Edge(VPointer v1,VPointer v2,int _type)
:type(_type)
{
   iter=NULL;
   if   (v1->getIndex()<v2->getIndex())
   {
      // add halfedge
      he[0]=v1->addHalfedge(v2,this);
      he[1]=v2->addHalfedge(v1,this);
   }
   else
   {
      he[0]=v2->addHalfedge(v1,this);
      he[1]=v1->addHalfedge(v2,this);
   }
   he[0]->inv_he=he[1];
   he[1]->inv_he=he[0];
}

//delete related halfedges
Edge::~Edge()
{
   he[0]->v->deleteFromOutHalfedgeList(he[1]);
   he[1]->v->deleteFromOutHalfedgeList(he[0]);

   delete he[0];
   delete he[1];

}

Slab::Slab(int li,int ri)
:l(li),r(ri)
{
   q_head=q_tail=NULL;
   next=NULL;
}

//assume q_head is not null
void Slab::pushBackQuad(QPointer q)
{
/*#ifdef   DEBUG
   assert(q_tail!=NULL);
#endif
*/   
   if   (q_tail==NULL)
   {
      q_head=q_tail=q;
   }
   else
   {
      q_tail->next=q;
      q_tail=q;
   }
}

Quad::Quad(VPointer hv,VPointer lv,HePointer t,HePointer b,int   _l,int _r)
:top_v(hv),btm_v(lv),top_he(t),btm_he(b),l(_l),r(_r)
{
   left_top_y=left_btm_y=right_top_y=right_btm_y=0;
   next=NULL;
}

Quad::Quad()
{
   top_v=btm_v=NULL;
   top_he=btm_he=NULL;
   l=r=0;
   left_top_y=left_btm_y=right_top_y=right_btm_y=0;
   next=NULL;
}

void Quad::updateLongEdge(EPointer e)
{
   double ey;   // y value of intersection point of e and vertical line through x 

   VPointer v1=e->he[1]->getFrontVertex();
   VPointer v2=e->he[0]->getFrontVertex();

   // y value of the intersection of e and the vertical line through (x,y)
   ey=yOfLineGivenX(v1->pos.x,v1->pos.y,v2->pos.x,v2->pos.y,top_v->pos.x);


   if   (ey>top_v->pos.y   &&   ey<top_y)
   {
      top_he=e->he[1];
      top_y=ey;
   }
   else
   if   (ey<btm_v->pos.y   &&   ey>btm_y)
   {
      btm_he=e->he[0];
      btm_y=ey;
   }
   else
   if   (ey==top_y   ||   ey==btm_y)
   {
      reportError("Intersection detected!",__LINE__,__FILE__);
   }
#ifdef   DEBUG
   else
   if   (ey==top_v->pos.y   ||   ey==btm_v->pos.y)
   {
      reportWarningString("Collinear degeneracy but does not matter the program!");
   }
#endif
}

void Quad::computeCornorY(double lx,double rx)
{
   VPointer v1=top_he->getFrontVertex();
   VPointer v2=top_he->getBackVertex();

   // y value of the intersection of e and the vertical line through (x,y)
   this->left_top_y=yOfLineGivenX(v1->pos.x,v1->pos.y,v2->pos.x,v2->pos.y,lx);
   this->right_top_y=yOfLineGivenX(v1->pos.x,v1->pos.y,v2->pos.x,v2->pos.y,rx);