cdt.cpp

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

         if   (p_he!=NULL)
         {
            //(p_he!=NULL)
            //candidate search given one
            tmphe=p_he->next_he;
            v_i=tmphe->v->getIndex();
            if   (v_i>p0->getIndex()   ||   v_i<l)
               tmphe=NULL;

            while(tmphe!=NULL
               &&   (vertexToVertexVertex(tmphe->v,p0,q0)==1))
            {
               p0=tmphe->v;
               tmphe=tmphe->inv_he->next_he;

               v_i=tmphe->v->getIndex();
               if   (v_i>p0->getIndex()   ||   v_i<l)
                  tmphe=NULL;
            }
            p_he=tmphe;
         }
         if   (q_he!=NULL)
         {
            int   v_i=q_he->v->getIndex();

            if   (v_i<m   ||v_i>=r)
               q_he=NULL;
         }
         // p_he is in wedge, must advance
         while   (q_he!=NULL   
            &&   (vertexToVertexVertex(q_he->v,p0,q0)!=-1))
         {
            q0=q_he->v;
            q_he=q_he->inv_he->next_he;
            int   v_i=q_he->v->getIndex();

            if   (v_i<m   ||v_i>q0->getIndex())
               q_he=NULL;
            
            if   (p_he!=NULL)
            {
               //advance ,same to above
               //candidate search given one
//!!!changed. fix bug 1204, it is not same to above, actually, here we don't need to initialize q_he like above
//               tmphe=p_he->next_he;
//               v_i=tmphe->v->getIndex();
//               if   (v_i>p0->getIndex()   ||   v_i<l)
//                  tmphe=NULL;
//replaced with
//since we have checked tmphe->v->index above, we don't need to do it here
               tmphe=p_he;

               while(tmphe!=NULL
                  &&   (vertexToVertexVertex(tmphe->v,p0,q0)==1))
               {
                  p0=tmphe->v;
                  tmphe=tmphe->inv_he->next_he;

                  v_i=tmphe->v->getIndex();
                  if   (v_i>p0->getIndex()   ||   v_i<l)
                     tmphe=NULL;
               }
               p_he=tmphe;
            }
         }
         //find it,test whether exists
         tmphe=p0->findOutHalfedgeInRing(q0);
         return   (tmphe==NULL)?addEdge(p0,q0):tmphe->e;
      }
   }
   else
   {
      if   (!posSlope)
      {
         if   (q_he!=NULL)
         {
            //(q_he!=NULL)
            //candidate search given one
            tmphe=q_he->next_he;
            v_i=tmphe->v->getIndex();
            if   (v_i<q0->getIndex()   ||   v_i>=r)
               tmphe=NULL;

            while(tmphe!=NULL
               &&   (vertexToVertexVertex(tmphe->v,p0,q0)==-1))
            {
               q0=tmphe->v;
               tmphe=tmphe->inv_he->next_he;
               v_i=tmphe->v->getIndex();
               if   (v_i<q0->getIndex()   ||   v_i>=r)
                  tmphe=NULL;
            }
            q_he=tmphe;
         }
         if   (p_he!=NULL)
         {
            int   v_i=p_he->v->getIndex();

            if   (v_i>=m   || v_i<l)
               p_he=NULL;
         }
         // p_he is in wedge, must advance
         while   (p_he!=NULL   
            &&   (vertexToVertexVertex(p_he->v,p0,q0)!=1))
         {
            p0=p_he->v;
            p_he=p_he->inv_he->next_he;
            int   v_i=p_he->v->getIndex();

            if   (v_i>=m   || v_i<p0->getIndex())
               p_he=NULL;
            
            if   (q_he!=NULL)
            {
               //advance ,same to above
               //candidate search given one
//!!!changed. fix bug 1204, it is not same to above, actually, here we don't need to initialize q_he like above
//               tmphe=q_he->pre_he;
//               v_i=tmphe->v->getIndex();
//               if   (v_i<q0->getIndex()   ||   v_i>=r)
//                  tmphe=NULL;
//replaced with
//since we have checked tmphe->v->index above, we don't need to do it here
               tmphe=q_he;

               while(tmphe!=NULL
                  &&   (vertexToVertexVertex(tmphe->v,p0,q0)==-1))
               {
                  q0=tmphe->v;
                  tmphe=tmphe->inv_he->next_he;
                  v_i=tmphe->v->getIndex();
                  if   (v_i<q0->getIndex()   ||   v_i>=r)
                     tmphe=NULL;
               }
               q_he=tmphe;
            }
         }
         //find it,test whether exists
         tmphe=p0->findOutHalfedgeInRing(q0);
         return   (tmphe==NULL)?addEdge(p0,q0):tmphe->e;
      }
      else
      {
         //(p_he!=NULL)
         //actually the following just change the above counter-part 
         // all next_he to pre_he, all 1 to -1

         if   (p_he!=NULL)
         {
            //candidate search given one
            tmphe=p_he->pre_he;
            v_i=tmphe->v->getIndex();
            if   (v_i>p0->getIndex()   ||   v_i<l)
               tmphe=NULL;

            while(tmphe!=NULL
               &&   (vertexToVertexVertex(tmphe->v,p0,q0)==-1))
            {
               p0=tmphe->v;
               tmphe=tmphe->inv_he->pre_he;

               v_i=tmphe->v->getIndex();
               if   (v_i>p0->getIndex()   ||   v_i<l)
                  tmphe=NULL;
            }
            p_he=tmphe;
         }
         if   (q_he!=NULL)
         {
            int   v_i=q_he->v->getIndex();
            if   (v_i<m   ||v_i>=r)
               q_he=NULL;
         }
         // p_he is in wedge, must advance
         while   (q_he!=NULL   
            &&   (vertexToVertexVertex(q_he->v,p0,q0)!=1))
         {
            q0=q_he->v;
            q_he=q_he->inv_he->pre_he;

            int   v_i=q_he->v->getIndex();
            if   (v_i<m   ||v_i>q0->getIndex())
               q_he=NULL;
            
            if   (p_he!=NULL)
            {
               //advance ,same to above
               //candidate search given one
//!!!changed. fix bug 1204, it is not same to above, actually, here we don't need to initialize q_he like above
//               tmphe=p_he->next_he;
//               v_i=tmphe->v->getIndex();
//               if   (v_i>p0->getIndex()   ||   v_i<l)
//                  tmphe=NULL;
//replaced with
//since we have checked tmphe->v->index above, we don't need to do it here
               tmphe=p_he;

               while(tmphe!=NULL
                  &&   (vertexToVertexVertex(tmphe->v,p0,q0)==-1))
               {
                  p0=tmphe->v;
                  tmphe=tmphe->inv_he->pre_he;

                  v_i=tmphe->v->getIndex();
                  if   (v_i>p0->getIndex()   ||   v_i<l)
                     tmphe=NULL;
               }
               p_he=tmphe;
            }
         }
         //find it,test whether exists
         tmphe=p0->findOutHalfedgeInRing(q0);
         return   (tmphe==NULL)?addEdge(p0,q0):tmphe->e;
      }
   }
}

EPointer   Mesh::findBridge(WPointer w,bool lowside)
{
   VPointer   p0,q0;
   HePointer   p_he,q_he;
   //get the initial start candidate pair
   // p0 is the lowest vertex in window which lies left to separate line
   findInitialCandidatePair(w,p0,q0,p_he,q_he,lowside);

   return   advancePairToBridge(w,p0,q0,p_he,q_he,lowside);
}

EPointer   Mesh::findNextBaseEdge(WPointer w,EPointer e)
{
   VPointer   p0,p;
   VPointer   q0,q;
   HePointer   p_he,q_he;

   int   m=w->lq->r;
   int   l=w->lq->l;
   int   r=w->rq->r;

   p0=e->he[1]->v;
   q0=e->he[0]->v;
   //set p_he as (p0->p) which is the first halfedge above e
   p_he=e->he[0]->pre_he;
   q_he=e->he[1]->next_he;

   p=p_he->v;
   q=q_he->v;

   VPointer   tmpv;
   HePointer   tmphe;
   //if p0 has only one edge(=e), q0 must have more than one edge 
   //because otherwise low_e==high_e 
   //so we just create new base above it using the safe q
   if   (p_he==e->he[0])
   {
      //find safe q
      tmphe=q_he->inv_he->next_he;
      tmpv=tmphe->v;
      while   (   q_he->e->type!=CONSTRAINED
         &&      (vertexToCircle(tmpv,q,p0,q0)==1)   )
      {
   #ifdef   DEBUG
         assert(q_he->v!=p0); //can't rotate back
   #endif
         deleteEdge(q_he->e);

         q_he=tmphe;
         q=tmpv;
         tmphe=q_he->inv_he->next_he;
         tmpv=tmphe->v;
      }
      return   addEdge(p0,q);
   }

   if   (q_he==e->he[1])
   {
      //find safe p
      tmphe=p_he->inv_he->pre_he;
      tmpv=tmphe->v;
      while   (   p_he->e->type!=CONSTRAINED
         &&      (vertexToCircle(tmpv,p,p0,q0)==1)   )
      {
   #ifdef   DEBUG
         assert(p_he->v!=q0); //can't rotate back
   #endif
         deleteEdge(p_he->e);

         //forward
         p_he=tmphe;
         p=tmpv;
         tmphe=p_he->inv_he->pre_he;
         tmpv=tmphe->v;
      }
      return   addEdge(p,q0);
   }

   // since it would not be the current base edge
   //if p_he->e intersect the middle separate line
   // it must be a constrained base edge
//   if   (p_he->v->getIndex()>=m)   return   (p_he->e);
//   if   (q_he->v->getIndex()<m)      return   (q_he->e);

   int valid;
   //find safe p
   tmphe=p_he->pre_he;
   //tmphe must above p_he
   tmpv=tmphe->v;
   
   //!!!changed 1204
//   if   (   tmpv->getIndex()>=r
   if   (   tmpv->getIndex()>=m
      ||   tmpv->getIndex()<l)
   {
      valid=-1;
   }
   else
   {
      valid=vertexToVertexVertex(tmpv,q0,p0);
   }
#ifdef   DEBUG
   if   (valid==0)   reportError("Collinear degeneracy find",__LINE__,__FILE__);
#endif

   while   (   p_he->e->type!=CONSTRAINED
      &&      valid==1
      &&      (vertexToCircle(tmpv,p,p0,q0)==1)   )
   {
#ifdef   DEBUG
      assert(p_he->v!=q0); //can't rotate back
#endif
      deleteEdge( p_he->e);

      //forward
      p_he=tmphe;
      p=tmpv;
      tmphe=p_he->pre_he;
      tmpv=tmphe->v;
      if   (   tmpv->getIndex()>=r
         ||   tmpv->getIndex()<l)
      {
         valid=-1;
      }
      else
      {
         valid=vertexToVertexVertex(tmpv,q0,p0);
      }
#ifdef   DEBUG
      if   (valid==0)   reportError("Collinear degeneracy find",__LINE__,__FILE__);
#endif
   }

   //find safe q
   tmphe=q_he->next_he;
   tmpv=tmphe->v;

   //!!!changed 1204
   if   (   tmpv->getIndex()>=r
//      ||   tmpv->getIndex()<l)
      ||   tmpv->getIndex()<m)
   {
      valid=-1;
   }
   else
   {
      valid=vertexToVertexVertex(tmpv,q0,p0);
   }
#ifdef   DEBUG
   if   (valid==0)   reportError("Collinear degeneracy find",__LINE__,__FILE__);
#endif
   while   (   q_he->e->type!=CONSTRAINED
      &&      valid==1
      &&      (vertexToCircle(tmpv,q,p0,q0)==1)   )
   {
#ifdef   DEBUG
      assert(q_he->v!=p0); //can't rotate back
#endif
      deleteEdge( q_he->e);

      q_he=tmphe;
      q=tmpv;
      tmphe=q_he->next_he;
      tmpv=tmphe->v;
      if   (   tmpv->getIndex()>=r
         ||   tmpv->getIndex()<l)
      {
         valid=-1;
      }
      else
      {
         valid=vertexToVertexVertex(tmpv,q0,p0);
      }
#ifdef   DEBUG
      if   (valid==0)   reportError("Collinear degeneracy find",__LINE__,__FILE__);
#endif
   }


   //test again to speed up
   // since it would not be the current base edge
   //if p_he->e intersect the middle separate line
   // it must be a constrained base edge
   // the above statement is wrong
//   if   (p_he->v->getIndex()>=m)   return   (p_he->e);
//   if   (q_he->v->getIndex()<m)      return   (q_he->e);

   if   (p==q)
   {
      //if p_he->e intersect the middle separate line
      // it must be a constrained base edge
      return   (p->getIndex()>=m?p_he->e:q_he->e);
   }

   //otherwise compare p,q
   if   (vertexToCircle(p,q,p0,q0)==1)
   {
      //p in O(q,p0,q0), select p
      return   addEdge(p,q0);
   }
   else
   {
      return   addEdge(p0,q);
   }
}

void   Mesh::display()
{
   //display

   VPointer   v1,v2;
   glBegin(GL_LINES);
   //glLineWidth(2.0);

   EListPointer   el;
   el=e_head;
   while   (el!=NULL)
   {
      if   (el->e->type!=VIRTUAL)
      {
         if   (el->e->type==CONSTRAINED)
         {
            glColor3f(1.0,.0,.0);
         }
         else
         {
            glColor3f(1.0,1.0,1.0);
         }
         v1=el->e->he[1]->v;
         v2=el->e->he[0]->v;

         glVertex2d(v1->pos.x,v1->pos.y);
         glVertex2d(v2->pos.x,v2->pos.y);
      }

      el=el->next;
   }
   glEnd();

   glBegin(GL_POINTS);
   glPointSize(2.0);
   glColor3f(.0,1.0,.0);
   for   (int i=0;i<v_num;i++)
      glVertex2d(v_array[i]->pos.x,v_array[i]->pos.y);
   glEnd();
}

void   Mesh::buildCDTOfWindow(WPointer w)
{
   EPointer low_e;
   EPointer high_e;
   bool   lowside=true;

   low_e=findBridge(w,lowside);
   high_e=findBridge(w,!lowside);

   //from low_e, repeatly add base edge until meet high_e
   while   (low_e!=high_e)
   {
      low_e=findNextBaseEdge(w,low_e);
   }
}

void   Mesh::mergeSlabPair(SPointer   ls)
{
   SPointer   rs;
   SPointer   ms;

   QPointer   tmpq;
   QPointer   newq;

   int      m,l,r;   // middle separator line's index

   rs=ls->next;

#ifdef DEBUG
   assert(ls->r==rs->l);