cdt.cpp

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

         {
            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;
            }
         }
         //left, then set the right one
         he=tmphe;
      }
      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;
            }
         }
      }
   }
   return   he;
}

// return he:
// if he->v->index is not in (sep_index-1, v->index) when v->=index>sep_index
// return NULL
HePointer Mesh::getHullHalfedgeNextTo(HePointer _he,VPointer v,int l,int sep_index,int r,bool lowside)
{
#ifdef   DEBUG
   assert(_he->getBackVertex()==v);
   assert(v->out_he!=NULL);
#endif
   HePointer he;
   
   int   v_i;
   v_i=v->getIndex();
   int   v_j;

   he=v->out_he;

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

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

   bool   leftv;
   leftv=v_i<sep_index;

   while   (he!=_he)
   {
      //if there is a bug, it will loop forever
      he=he->pre_he;
   }

   //now he=_he;
   if   (leftv)
   {
      if   (lowside)
      {
         he=he->pre_he;
         if   (he->v->getIndex()<v->getIndex()
            ||   he->v->getIndex()>=r
            ||   pointInAngle(v->pos.x,v->pos.y,he->v->pos.x,he->v->pos.y,v->pos.x,v->pos.y-1,_he->v->pos.x,_he->v->pos.y))
         {
            he=NULL;
         }
      }
      else
      {
         he=he->next_he;
         if   (he->v->getIndex()<v->getIndex()
            ||   he->v->getIndex()>=r
            ||   pointInAngle(v->pos.x,v->pos.y,he->v->pos.x,he->v->pos.y,_he->v->pos.x,_he->v->pos.y,v->pos.x,v->pos.y+1))
         {
            he=NULL;
         }
      }
//      v_j=he->v->getIndex();
//      return   (v_j>v_i&&v_j<sep_index)?he:NULL;
      return   he;
   }
   else
   {
      if   (lowside)
      {
         he=he->next_he;
         if   (he->v->getIndex()>v->getIndex()
            ||   he->v->getIndex()<l
            ||   pointInAngle(v->pos.x,v->pos.y,he->v->pos.x,he->v->pos.y,v->pos.x,v->pos.y+1,_he->v->pos.x,_he->v->pos.y))
         {
            he=NULL;
         }
      }
      else
      {
         he=he->pre_he;
         if   (he->v->getIndex()>v->getIndex()
            ||   he->v->getIndex()<l
            ||   pointInAngle(v->pos.x,v->pos.y,he->v->pos.x,he->v->pos.y,_he->v->pos.x,_he->v->pos.y,v->pos.x,v->pos.y-1))
         {
            he=NULL;
         }
      }
//      v_j=he->v->getIndex();
//      return   (v_j>=sep_index&&v_j<v_i)?he:NULL;
      return   he;
   }
}

// In this function, we get the initial candidate pair
// and corresponding halfedge on the hull which we can use to advance
// this procedure request too much exact computation!!!
// see 'diff' below
//# EXACT
void   Mesh::findInitialCandidatePair(WPointer w,
      VPointer &p0,VPointer &q0,
      HePointer &p_he,HePointer &q_he,
      bool   lowside)
{
   int   sign;

   //care!, here I flip the upsign to make this function more general to use
   sign=lowside?   w->downsign:   -(w->upsign);

   int   m=w->lq->r;
   int   l,r;
   l=w->lq->l;
   r=w->rq->r;
   HePointer   he;
   double   diff;
   bool test;

   //in this switch
   // we get the lowset(highest) vertices on either side
   // and corresponding halfedges
   //while there may be situation that more than one vertex is on the same horizontal line
   // then we need to adjust it:
   //for left, we move vertex to right most one and update its halfedge
   //for right, vice versa

   //find the lowest vertex
   switch(sign)
   {
   case   0:
      if   (lowside)
      {
         p0=w->lq->btm_v;
         q0=w->rq->btm_v;
      }
      else
      {
         p0=w->lq->top_v;
         q0=w->rq->top_v;
      }

      //get the halfedge of v between the vertical line through v->index and separate line
      //if lower side, it should be the first counter-clockwisely which is above v, otherwise vice versa
      // also may be horizental
      // possbile is null
      p_he=getHullHalfedge(p0,m,lowside);
      q_he=getHullHalfedge(q0,m,lowside);
      break;
   case   1:
      // q0 can be settled
      if   (lowside)
      {
         //right is higher, q0 is settled
         q0=w->rq->btm_v;
         //w->rq->top_he->v is a vertex in w->lq
         he=w->rq->btm_he;
         p0=he->inv_he->v;
      }
      else
      {
         //right is lower, 
         q0=w->rq->top_v;
         he=w->rq->top_he->inv_he;
         p0=he->inv_he->v;
      }
      q_he=getHullHalfedge(q0,m,lowside);
      //p_he is a halfedge of p0 which lies just above(if lowside=true)he
      //and between p0 and sep[m]
      p_he=getHullHalfedgeNextTo(he,p0,l,m,r,lowside);

      if   (p_he!=NULL&&(p_he->v->getIndex()<m&&p_he->v->getIndex()>p0->getIndex()))
      {
         //whether next v is lower(if lowside=true) than current, if so, advance
         //here although the inverse of '<' is '>=' instead of '>', 
         // it still work because we exclude that case
         diff=p_he->v->pos.y - p0->pos.y;
         if   (diff!=0)
         {
            test=diff<0;
            test=(lowside)?test:!test;
            //if test==true, the very special case occur
            //assume lowside==true, then
            // we know the lowest vertex can be find by repeating comparing y value
            // and it will never go back, it can be only lower and lower to the right side of v
            while   (test)//&&p_he!=NULL)
            {
               //advance
               p0=p_he->v;
               if(lowside)
               {
                  p_he=p_he->inv_he->pre_he;
               }
               else
               {
                  p_he=p_he->inv_he->next_he;
               }

               if   (p_he->v->getIndex()>=m)
               {
               //   p_he=NULL;
                  test=false;
               }
               else
               {
                  //test again
               diff=p_he->v->pos.y - p0->pos.y;
               if   (diff==0)
               {
                  test=false;
               }
               else
               {
                  test=diff<0;
                  test=(lowside)?test:!test;
               }
            }
            }
         }
         /*
         if   (p_he!=NULL   &&(p_he->v->getIndex()>=m||p_he->v->getIndex()<p0->getIndex()) )
         {
            //has found some lowest(if lowside=true) vertex
            //while its hull edge intersect the middle separator line
            p_he=NULL;
         }
         */
      }
      break;
   case   -1:
      // p0 can be settled
      if   (lowside)
      {
         //left is higher
         p0=w->lq->btm_v;
         //w->rq->top_he->v is a vertex in w->lq
         he=w->lq->btm_he->inv_he;
         q0=he->inv_he->v;
      }
      else
      {
         //right is lower
         p0=w->lq->top_v;
         he=w->lq->top_he;
         q0=he->inv_he->v;
      }
      p_he=getHullHalfedge(p0,m,lowside);
      q_he=getHullHalfedgeNextTo(he,q0,l,m,r,lowside);

      if   (q_he!=NULL   && (q_he->v->getIndex()>=m&&q_he->v->getIndex()<q0->getIndex()) )
      {
         diff=q_he->v->pos.y - q0->pos.y;
         if   (diff!=0)
         {
            test=diff<0;
            test=(lowside)?test:!test;
            while   (test)//&&q_he!=NULL)
            {
               //advance
               q0=q_he->v;
               if(lowside)
               {
                  q_he=q_he->inv_he->next_he;
               }
               else
               {
                  q_he=q_he->inv_he->pre_he;
               }

               if   (q_he->v->getIndex()<m)
               {
               //   q_he=NULL;
                  test=false;
               }
               else
               {
                  //test again
                  diff=q_he->v->pos.y - q0->pos.y;
                  if   (diff==0)
                  {
                     test=false;
                  }
                  else
                  {
                     test=diff<0;
                     test=(lowside)?test:!test;
                  }
               }
            }
         }
         /*
         if   (q_he!=NULL   && (q_he->v->getIndex()<m||q_he->v->getIndex()>q0->getIndex()) )
         {
            q_he=NULL;
         }
         */
      }
      break;
   default:   break;
   }

   //if there are vertices which have same lowest y value, 
   //move to the vertex closet to middle separate line
   if   (p_he!=NULL)
   {
      //since we have check in the above
      //here the p_he->v is in the middle of p0 and separater line m
      // we only advance p0 if p_he->v is right to p0, left to m and diff==0
      test=(   p_he->v->getIndex()>p0->getIndex()   
         &&   p_he->v->getIndex()<m
         &&   p_he->v->pos.y == p0->pos.y   );

      while   (test)
      {
         //advance
         p0=p_he->v;
         if(lowside)
         {
            p_he=p_he->inv_he->pre_he;
         }
         else
         {
            p_he=p_he->inv_he->next_he;
         }

         if   (   p_he->v->getIndex()>=m
            ||   p_he->v->getIndex()<p0->getIndex())
         {
         //   p_he=NULL;
            test=false;
         }
         else
         {
            //test again
            test=(p_he->v->pos.y == p0->pos.y);
         }
      }
   }
   //advance q0 in similar way
   if   (q_he!=NULL)
   {
      test=(   q_he->v->getIndex()<q0->getIndex()   
         &&   q_he->v->getIndex()>=m
         &&   q_he->v->pos.y == q0->pos.y   );

      while   (test)
      {
         //advance
         q0=q_he->v;
         if(lowside)
         {
            q_he=q_he->inv_he->next_he;
         }
         else
         {
            q_he=q_he->inv_he->pre_he;
         }

         if   (   q_he->v->getIndex()<m
            ||   q_he->v->getIndex()>q0->getIndex())
         {
         //   q_he=NULL;
            test=false;
         }
         else
         {
            //test again
            test=(q_he->v->pos.y == q0->pos.y);
         }
      }
   }
}

EPointer   Mesh::advancePairToBridge(WPointer w,VPointer p0,VPointer q0,HePointer p_he,HePointer q_he,bool lowside)
{
   double   dy;
   dy=(q0->pos.y - p0->pos.y);
   HePointer he;

   if   (q_he==NULL   &&   p_he==NULL)
   {
      he=p0->findOutHalfedgeInRing(q0);
      return   (he==NULL)?addEdge(p0,q0):he->e;
   }

   if   (dy==0   )
   {
      he=p0->findOutHalfedgeInRing(q0);
      return   (he==NULL)?addEdge(p0,q0):he->e;
   }

   //else dy<>0
   int   l,m,r;
   l=w->lq->l;
   m=w->rq->l;
   r=w->rq->r;

   bool   posSlope;
   posSlope=dy>0;

   HePointer   tmphe;
   int   v_i;
   if   (lowside)
   {
      if   (posSlope)
      {
         if   (q_he!=NULL)
         {
            //(q_he!=NULL)
            //candidate search given one
            tmphe=q_he->pre_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->pre_he;
               v_i=tmphe->v->getIndex();
               if   (v_i<q0->getIndex()   ||   v_i>=r)
                  tmphe=NULL;
            }
            q_he=tmphe;
         }
         if   (p_he!=NULL)
         {
         // p_he is in wedge, must advance
            int   v_i=p_he->v->getIndex();

            if   (v_i>=m   || v_i<l)
               p_he=NULL;
         }
         while   (p_he!=NULL   
            &&   (vertexToVertexVertex(p_he->v,p0,q0)!=-1))
         {
            p0=p_he->v;
            p_he=p_he->inv_he->pre_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->pre_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
      {