cdt.cpp

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

   v1=btm_he->getBackVertex();
   v2=btm_he->getFrontVertex();
   this->left_btm_y=yOfLineGivenX(v1->pos.x,v1->pos.y,v2->pos.x,v2->pos.y,lx);
   this->right_btm_y=yOfLineGivenX(v1->pos.x,v1->pos.y,v2->pos.x,v2->pos.y,rx);
}

Mesh::Mesh()
{
   e_head=e_tail=NULL;
   v_num=e_num=ce_num=0;
   maxs.y=maxs.x=HUGE_VAL;
   mins.y=mins.x=-HUGE_VAL;
   v_array=NULL;
   sep_array=NULL;
}

//***************** quick sort functions end  ***********************
// sort vertices by x-coordinate
void Mesh::sortVertices()
{
   // implemented in quick sort
   quicksort(v_array,0,v_num-1);      
   // update order information of vertices
/*not here
   for   (int   vj=0;vj<v_num;vj++)
   {
      v_array[vj]->setIndex(vj);
   }*/

}

EPointer Mesh::addEdge(VPointer v1,VPointer v2,int type)
{
   /*
#ifdef   DEBUG
   // no conflicts
   if   (! (v1->findOutHalfedgeInRing(v2)==NULL   ) )//&&   v2->findOutHalfedgeInRing(v1)==NULL))
   {
      reportError("Duplicate edge!",__LINE__,__FILE__);
   }
#endif
*/
   EPointer   e;
   EListPointer   el;
   
   // virtual edge are not take into account
   if   (type==VIRTUAL)
   {
      e=new Edge(v1,v2,CONSTRAINED);
   }
   else
   {
      e=new Edge(v1,v2,type);

      e_num++;
   
      el=new EdgeList(e);

      if   (type==CONSTRAINED)   ce_num++;

      if   (e_head==NULL)
      {
         e_head=e_tail=el;
      }
      else
      {
         e_tail->next=el;
         e_tail=el;
      }
   }
#ifdef OUTPUT
   printf("Edge added:(%d,%d)\n",v1->getIndex(),v2->getIndex());
#endif
   return e;
}

void   Mesh::deleteEdge(EPointer e)
{
   EListPointer el;
   el=e->iter;

#ifdef   OUTPUT
   printf("Delete edge:(%d,%d)\n",e->he[1]->v->getIndex(),e->he[0]->v->getIndex());
#endif

   if   (e_head==el)
   {
      if   (e_tail==e_head)
         e_tail=NULL;

      e_head=e_head->next;
   }
   else
   {
      EListPointer tmpe;
      tmpe=e_head;
      while   (tmpe->next!=el)
         tmpe=tmpe->next;
      
      if   (e_tail==el)
      {
         e_tail=tmpe;
      }
      
      tmpe->next=el->next;
   }
   
   delete e;
}

/* Read the source file and convert to data structure
   Input: sf is the source file name, the format refers to above
   Output:   0 means no error
*/
int   Mesh::import2OFFFile(char* sf){
   FILE * fp;

   if
   ((fp=fopen(sf,"r"))==NULL)
   {
      reportError("Cannot open this file!",__LINE__,__FILE__);
   }
   else
   {

/*      char   line[LINE_LENGTH];

      fgets(line,LINE_LENGTH,fp);
*/
      omitComment(fp);

      // read OFF, currently only support 'OFF'
      // file name length is maximum 256
      char   keyword[256];
      fscanf(fp,"%s",keyword);

      if   (strcmp((keyword),noff_keyword)!=0 )
      {
         reportError("Not valid format(*.noff)!",__LINE__,__FILE__);
      }

      readToLineEnd(fp);
      omitComment(fp);

      int   _d;
      fscanf(fp,"%d",&_d);
      if   (_d!=2)
      {
         reportError("Not valid dimension!",__LINE__,__FILE__);
      }
      readToLineEnd(fp);
      omitComment(fp);

      int   vnum;
      // read vertices count, faces count and edges count respectively
      // read num
      fscanf(fp,"%d",&vnum);
      if   (vnum<1)
      {
         reportError("Vertex number should above zero.",__LINE__,__FILE__);
      }

      v_array=new VPointer[vnum];

      int   fnum;   // f_num is not checked
      fscanf(fp,"%d",&fnum);

      int cenum;
      fscanf(fp,"%d",&cenum);
      readToLineEnd(fp);
      omitComment(fp);

      double   x;
      double   y;

      //init 
      mins.x=mins.y=HUGE_VAL;
      maxs.x=maxs.y=-HUGE_VAL;

      // read information of vertices
      for (int i=0;i<vnum;i++)
      {
         
         fscanf(fp,"%lf%lf",&x,&y);

         maxs.x=MAX(maxs.x,x);
         maxs.y=MAX(maxs.y,y);
         mins.x=MIN(mins.x,x);
         mins.y=MIN(mins.y,y);

         v_array[i]=new Vertex(x,y,i);
         v_num++;

         readToLineEnd(fp);         
         
         omitComment(fp);
      }

      //vnum should eq v_num now

      //sort the vertex
      sortVertices();

      //get the information of sorting
      int*   reorder=new int[v_num];
      for   (int ni=0;ni<v_num;ni++)
      {
         reorder[v_array[ni]->getIndex()]=ni;
      }

      //update the information of vertices;
      for   (int nj=0;nj<v_num;nj++)
      {
         v_array[nj]->setIndex(nj);
      }
      //omit faces
      int   side_num;
      for   (int   f_i=0;f_i<fnum;f_i++)
      {
         //read side num
         fscanf(fp,"%d",&side_num);
         int vi;
         for   (int v_i=0;v_i<side_num;v_i++)
         {
            fscanf(fp,"%d",&vi);
         }
            readToLineEnd(fp);         
            omitComment(fp);
      }

      //read constrained edge
      for (int j=0;j<cenum;j++)
      {
         int   i1,i2;

         //get first constrained edge
         fscanf(fp,"%d%d",&i1,&i2);
      
         if   (i1<0 || i1>=v_num || i2<0 || i2>=v_num)
               reportError("Incorrect input of edge's endpoint's index. It should not below 0 or large or equal to vertex total count.",__LINE__,__FILE__);      

         addEdge(v_array[reorder[i1]],v_array[reorder[i2]],CONSTRAINED);

         readToLineEnd(fp);
         omitComment(fp);
      }   

      //ce_num shoulde eq cenum now

      fclose(fp);   
      //sort the v_array while let the index in each vertex unchanged
   }

   return(0);
}

int   Mesh::export2OFFFile(char* sf,char* f)
{
   char delimiter=' ';
   char* df="result.noff";
   FILE* fp;   

   if(f!=NULL)
   {
      df=f;   
   }

   if
   ((fp=fopen(df,"w"))==NULL)
   {
      reportError("Cannot create or overwrite this file!",__LINE__,__FILE__);
   }
   else
   {
      // add a mark in the edge: writed
#ifdef   COMMENT_ON
      fprintf(fp,"# generated from source file: %s \n",sf);
#endif
      fprintf(fp,"%s\n",noff_keyword);
      fprintf(fp,"%d\n",2);

      //output vertex count
#ifdef   COMMENT_ON
      fprintf(fp,"%s\n","# vertices num, face num=0, edge num");
#endif
      fprintf(fp,"%d%c%d%c%d\n",v_num,delimiter,0,delimiter,e_num);

      //output vertex information in form of x, y value
      for (int i=0;i<v_num;i++)
      {
#ifdef   COMMENT_ON
         fprintf(fp,"%s %d\n","# vertex index: ",i);
#endif
         fprintf(fp,"%f %f\n",v_array[i]->pos.x,v_array[i]->pos.y);
      }
      

      //output edge information in form of i1,i2 value
      EListPointer el;
      el=e_head;
      while(el!=NULL)
      {
         VPointer v1=el->e->he[0]->getBackVertex();
         VPointer v2=el->e->he[0]->getFrontVertex();
#ifdef   COMMENT_ON

         fprintf(fp,"%s%f,%f%s%f,%f%s\n","# (",v1->pos.x,v1->pos.y,")--(",v2->pos.x,v2->pos.y,")");
#endif
         fprintf(fp,"%d %d\n",v1->getIndex(),v2->getIndex());

         el=el->next;
      }

      fclose(fp);   
   }
   return (0);
}

/* return 0 means no error*/
int Mesh::initSlabs()
{
   // the index itself shows it's position
//   sortVertices();

#define   OFFSET 1

   //index is -2,-1
   vlt=new Vertex(mins.x-3*OFFSET,maxs.y+OFFSET,-2);
   vlb=new Vertex(mins.x-2*OFFSET,mins.y-OFFSET,-1);

   //index is v_num+1,v_num
   vrt=new Vertex(maxs.x+3*OFFSET,maxs.y+OFFSET,v_num+1);
   vrb=new Vertex(maxs.x+2*OFFSET,mins.y-OFFSET,v_num);

   // 2 long horizontal lines
   sky_e=addEdge(vlt,vrt,VIRTUAL);
   eth_e=addEdge(vlb,vrb,VIRTUAL);

   //init separate lines
   sep_array=new double[v_num+1];

   sep_array[0]=mins.x-OFFSET;
   sep_array[v_num]=maxs.x+OFFSET;
#ifdef DEBUG
      printf("%s%d%s%f\n","# add separator line ",0," x:=",sep_array[0]);
      printf("%s%d%s%f\n","# add separator line ",v_num," x:=",sep_array[v_num]);
#endif

   double x1;
   double x2;
   x1=v_array[0]->pos.x;

   for (int v_i=1;v_i<v_num;v_i++)
   {   
      x2=v_array[v_i]->pos.x;
      sep_array[v_i]=(x1+x2)/2;
      x1=x2;
#ifdef DEBUG
      printf("%s%d%s%f\n","# add separator line ",v_i," x:=",sep_array[v_i]);
#endif
   }

   //init slabs
   SPointer*slab_array=new SPointer[v_num];
   // create an array to accelerate the initialization of quad and slab
   QPointer tmpq;
   VPointer tmpv;

   //other slabs
   for (int v_j=0;v_j<v_num;v_j++)
   {
      slab_array[v_j]=new Slab(v_j,v_j+1);
      //initially, each quad has sky and earth as its top and bottom halfedge
      tmpv=v_array[v_j];
      tmpq=new Quad(tmpv,tmpv,sky_e->he[1],eth_e->he[0],v_j,v_j+1);
      tmpq->top_y=maxs.y+OFFSET;
      tmpq->btm_y=mins.y-OFFSET;

      slab_array[v_j]->q_head=slab_array[v_j]->q_tail=tmpq;
   }

   for (int k=0;k<v_num-1;k++)
   {
      slab_array[k]->next=slab_array[k+1];
   }

   // init slab list's head
   s_head=slab_array[0];

   //   init quad for each slab by comparing each constrained edges indices to the indices of separator lines
   EListPointer   tmpe;
   tmpe=e_head;
   while(tmpe!=NULL)
   {
      int   l,r;
      l=tmpe->e->he[1]->v->getIndex();
      r=tmpe->e->he[0]->v->getIndex();

      for(int m=l+1;m<r;m++)
      {
         // update slab m's initial quad with constrained edge tmpe->e
         slab_array[m]->q_head->updateLongEdge(tmpe->e);
      }
      tmpe=tmpe->next;
   }

   // update quad's cornor's y value
   for (int s_i=0;s_i<v_num;s_i++)
      slab_array[s_i]->q_head->computeCornorY(sep_array[s_i],sep_array[s_i+1]);

   //release slab_array's space, not the slab's space
   delete []slab_array;
   return 0;
}

//get the hull edge here is different from Halfedge::getNext()
// here, we want to find one halfedge he for v such that
//if v->index <sep_index
// he->v is (for lowside=true,below otherwise) 
// the first halfedge just right to the lower vertical ray start from v
// if no such halfedge exist, return null
//# EXACT
HePointer   Mesh::getHullHalfedge(VPointer v,int sep_index,bool lowside)
{
   HePointer   he;
   HePointer   tmphe;

   Vector   in;

   int   v_i;
   v_i=v->getIndex();
   int   v_j;

   he=v->out_he;
   if   (he==NULL)   return   NULL;

   v_j=he->v->getIndex();

   bool   leftv;
   leftv=v_i<sep_index;

   //if just on the same horizontal
   if   (he->d.x==0)
   {
      return   he;
   }

   bool inangle;
   //else
   //here we used the property that no two vertices has the same x value
   // so the loop below won't run forever
   if   (lowside)
   {
      in.x=Y.x;
      in.y=-Y.y;
      if   (leftv)
      {
         tmphe=he->next_he;
         if   (tmphe!=he)   //test if there is only one halfedge
         {
            inangle=pointInAngle(0,0,in.x,in.y,tmphe->d.x,tmphe->d.y,he->d.x,he->d.y) == 1;
            while   (!inangle)
            {
               he=tmphe;
               tmphe=tmphe->next_he;
               inangle=pointInAngle(0,0,in.x,in.y,tmphe->d.x,tmphe->d.y,he->d.x,he->d.y) == 1;
            }
         }
      }
      else
      {
         tmphe=he->next_he;
         if   (tmphe!=he)   //test if there is only one halfedge
         {
            inangle=pointInAngle(0,0,in.x,in.y,tmphe->d.x,tmphe->d.y,he->d.x,he->d.y) == 1;
            while   (!inangle)
            {
               he=tmphe;
               tmphe=tmphe->next_he;
               inangle=pointInAngle(0,0,in.x,in.y,tmphe->d.x,tmphe->d.y,he->d.x,he->d.y) == 1;
            }
            //right, then set the left one
            he=tmphe;
         }
      }
   }
   else
   {
      in.x=Y.x;
      in.y=Y.y;
      if   (leftv)
      {
         tmphe=he->next_he;
         if   (tmphe!=he)   //test if there is only one halfedge