Home
Fractals
Tutorials
Books
Archive
My blog
My LinkedIn Profile

BOOKS i'm reading

Napoleon Hill Keys to Success: The 17 Principles of Personal Achievement, Napoleon Hill, ISBN: 978-0452272811
The 4-Hour Workweek: Escape 9-5, Live Anywhere, and Join the New Rich (Expanded and Updated), Timothy Ferriss, ISBN: 978-0307465351
The Fountainhead, Ayn Rand, ISBN: 0452273331

/*
 * Module ID: hedge.cpp
 * Titre    : Declaration de la classe halfEdge.
 * Utilite  : Classe utilisee pour la construction du diagramme de Voronoi
 *
 * Auteur   : Olivier Langlois (olivier@olivierlanglois.net)
 * Date     : 24 Fevrier 1998
 */

#include "debug.h"
#include "collect/idlist.h"
#include "geo/point.h"
#include "geo/segment.h"
#include "geo/line.h"
#include "geo/circle.h"
#include "collect/colvar.h"
#include "geo/dcel.h"
#include "geo/hedge.h"
#include "prime.h"
#include "collect/contx.h"
#ifdef __GNUG__
#include "umath.h"
#else
#include <math.h>
#endif

GVectorimplement(halfEdgeP)

GVector(halfEdgeP) halfEdge::table_;
halfEdge *halfEdge::leftEnd  = NULL;
halfEdge *halfEdge::rightEnd = NULL;

halfEdge::halfEdge( edge *e, char dir )
{
  Vertex   = NULL;
  Previous = NULL;
  Next     = NULL;
  ve = e;
  direction = dir;

  V1isSetted = BOOL_FALSE;
  V2isSetted = BOOL_FALSE;

  if( ve )
  {
     ve->addReference();
     point s1 = ve->getF1()->getLocation();
     point s2 = ve->getF2()->getLocation();

     segment s(s1,s2);
     c = s1.xcoord()*s.dx() + s1.ycoord()*s.dy() +
          (s.dx()*s.dx() + s.dy()*s.dy())*0.5;
     float adx = s.dx()>0?s.dx():-s.dx();
     float ady = s.dy()>0?s.dy():-s.dy();
     if( adx>ady )
        c /= s.dx();
     else
        c/= s.dy();
  }
}

halfEdge::~halfEdge()
{
  if( Previous && Next )
  {
     Previous->Next = Next;
     Next->Previous = Previous;
  }

  if( ve && !ve->removeReference() /*&&
      (V1isSetted == BOOL_TRUE || V2isSetted == BOOL_TRUE)*/)
  {
     if( V1isSetted == BOOL_FALSE && direction == 0)
     {
        vertex *v1 = ve->getV1();
        vertex *v2 = ve->getV2();

        // Si les 2 vertices sont du meme cote, une correction est nessecaire.
        if( v1->xcoord() >  v2->xcoord() ||
            (v1->xcoord() == v2->xcoord() && v1->ycoord() < v2->ycoord()) )
          ve->setV1( new vertex(v1->reflect( *v2 )) );
     }
     else if( V2isSetted == BOOL_FALSE && direction == 1)
     {
        vertex *v1 = ve->getV1();
        vertex *v2 = ve->getV2();

        // Si les 2 vertices sont du meme cote, une correction est nessecaire.
        if( v2->xcoord() <  v1->xcoord() ||
            (v2->xcoord() == v1->xcoord() && v2->ycoord() > v1->ycoord()) )
          ve->setV2( new vertex(v2->reflect( *v1 )) );
     }

     ve->insert();
  }
  removeYStar();

  // Enleve les reference de cette arete de la hash table.
  if( references() )
  {
     halfEdge *cur;
     for( int i = 0; i < table_.length(); i++ )
        if( (cur = table_(i)) == this )
        {
          if( i == 0 )
             table_(0) = leftEnd;
          else if( i == table_.length() - 1 )
             table_(i) = rightEnd;
          else table_(i) = NULL;
          if( !removeReference() ) break;
        }
  }

}

/******************************************************************************
 *
 * Nom       : init
 *
 ****************************************************************************/
void halfEdge::init( size_t N )
{
#ifdef __GNUG__
  table_.resize(setTable((unsigned)(2*usqrt(N))));
#else
  table_.resize(setTable(2*sqrt(N)));
#endif
  for( int i = 1; i < table_.length() - 1; i++ )
     table_(i) = NULL;

  leftEnd  = new halfEdge();
  rightEnd = new halfEdge();
  leftEnd->Previous  = leftEnd;
  leftEnd->Next      = rightEnd;
  rightEnd->Previous = leftEnd;
  rightEnd->Next     = rightEnd;

  leftEnd->addReference();
  rightEnd->addReference();

  table_(0)                 = leftEnd;
  table_(table_.length()-1) = rightEnd;
}

/******************************************************************************
 *
 * Nom       : setTable
 *
 ****************************************************************************/
size_t halfEdge::setTable(size_t N)
{
  size_t bucketNum, i = 0;

  do
  {
     bucketNum = primeTab[i];
  }
  while( primeTab[i++] < N );

  D( cout << "bucketNum   : " << bucketNum << endl; )
  return bucketNum;
}

/******************************************************************************
 *
 * Nom       : terminate
 *
 ****************************************************************************/
void halfEdge::terminate( void )
{
  halfEdge *cur;
  while( (cur = leftEnd->Next) != rightEnd )
     delete cur;
  delete leftEnd;
  delete rightEnd;
}

void halfEdge::setYStar( vertex *v, const point &p )
{
  Vertex = v;
  v->addReference();
  ystar = circle( *(point *)v, p );
}

void halfEdge::removeYStar( void )
{
  if( Vertex )
  {
     if( !Vertex->removeReference() )
        delete Vertex;
     Vertex = NULL;
  }
}

face *halfEdge::rightSite()
{
  if(direction)
     return ve->getF1();
  else return ve->getF2();
}

face *halfEdge::leftSite()
{
  if(direction)
     return ve->getF2();
  else return ve->getF1();
}

void halfEdge::setV1( vertex *v )
{
  if( V1isSetted == BOOL_TRUE )
     throw GeoEx(GX_VERTEX);

  // Switch les vertices si necessaire
  vertex *v2 = ve->getV2();
  if( *v2 == *v && V2isSetted == BOOL_FALSE )
  {
     v2 = ve->getV1();
     ve->setV2(v2);
  }
  ve->setV1(v);
  V1isSetted = BOOL_TRUE;
}

void halfEdge::setV2( vertex *v )
{
  if( V2isSetted == BOOL_TRUE )
     throw GeoEx(GX_VERTEX);

  // Switch les vertices si necessaire
  vertex *v1 = ve->getV1();
  if( *v1 == *v && V1isSetted == BOOL_FALSE )
  {
     v1 = ve->getV2();
     ve->setV1(v1);
  }
  ve->setV2(v);
  V2isSetted = BOOL_TRUE;
}

int halfEdge::operator==(const halfEdge &h) const
{
  return point::cmp_yx( *Vertex, *h.Vertex ) == 0;
/*
  return (Vertex->ycoord()+ystar.radius() ==
             h.Vertex->ycoord()+h.ystar.radius() &&
             Vertex->xcoord() == h.Vertex->xcoord());
*/
}

int halfEdge::operator!=(const halfEdge &h) const
{
  return (!operator==(h));
}

int halfEdge::operator <(const halfEdge &he) const
{
  // Note : Plus ystar est petit, plus la priorite est grande.
//  return point::cmp_yx( *Vertex, *h.Vertex ) > 0;
  halfEdge &h = const_cast< halfEdge & >(he);
  return (Vertex->ycoord()+const_cast<circle &>(ystar).radius() >
             h.Vertex->ycoord()+h.ystar.radius() ||
             (Vertex->ycoord()+const_cast<circle &>(ystar).radius() ==
              h.Vertex->ycoord()+h.ystar.radius() &&
              Vertex->xcoord() > h.Vertex->xcoord()));

}

int halfEdge::operator >(const halfEdge &he) const
{
  // Note : Plus ystar est petit, plus la priorite est grande.
//  return point::cmp_yx( *Vertex, *h.Vertex ) < 0;
  halfEdge &h = const_cast< halfEdge & >(he);
  return (Vertex->ycoord()+const_cast<circle &>(ystar).radius() <
             h.Vertex->ycoord()+h.ystar.radius() ||
             (Vertex->ycoord()+const_cast<circle &>(ystar).radius() ==
              h.Vertex->ycoord()+h.ystar.radius() &&
              Vertex->xcoord() < h.Vertex->xcoord()));

}

int halfEdge::operator<=(const halfEdge &h) const
{
  return ( operator<(h) || operator==(h) );
}

int halfEdge::operator>=(const halfEdge &h) const
{
  return ( operator>(h) || operator==(h) );
}

/******************************************************************************
 *
 * Nom       : insert
 *
 * Utilite   : Insere une nouvelle arete.
 *
 ****************************************************************************/
void halfEdge::insert( halfEdge *lb )
{
    Previous  = lb;
    Next      = lb->Next;
    (lb->Next)->Previous = this;
    lb->Next  = this;
}

/******************************************************************************
 *
 * Nom       : leftbnd
 *
 * Utilite   : Retrouve l'arete dont le point est le plus proche.
 *
 ****************************************************************************/
halfEdge *halfEdge::leftbnd(const point &p)
{
    halfEdge *e;
    halfEdge *closestEdge;
    halfEdge *tmpE;
    unsigned buck, i;

    D( e = leftEnd; )
    D( do {e = e->Next;if(e->ve) cout << e->ve->getSegment() << (e->direction?"right":"left") << e->ve->getF1()->getLocation() << e->ve->getF2()->getLocation() << endl;} while( e !=rightEnd ); )

    // Utilise la hash table pour etre proche de la halfEdge
    buck = (unsigned)((((float)p.hash())/UINT_MAX) * table_.length());
    if(buck>=table_.length()) buck = table_.length() - 1;
    e = table_(buck);

    if( e == NULL )
    {
         for(i=1; 1; i++)
         {
            if ((e=table_(buck-i)) !=  NULL) break;
            if ((e=table_(buck+i)) !=  NULL) break;
         }
    }

    // recherche lineairement la liste
    if( e==leftEnd  || (e != rightEnd && e->atRight(p) == BOOL_TRUE) )
    {
      closestEdge = e;
      do
      {
         e = e->Next;
      }
      while( e != rightEnd && e->atRight(p) == BOOL_TRUE );
      closestEdge = e->Previous;

      // Regarde a gauche
      while(
                closestEdge->Previous != leftEnd &&
                closestEdge->Previous->atRight(p) == BOOL_TRUE &&
                closestEdge->Previous->rightSite()->getLocation().distance(p) <
                closestEdge->rightSite()->getLocation().distance(p) )
         closestEdge = closestEdge->Previous;

      // Regarde a droite
      if( e != rightEnd )
         do
         {
            e = e->Next;
         }
         while( e != rightEnd && e->atRight(p) == BOOL_FALSE );

      if( e != rightEnd &&
            e->atRight(p) == BOOL_TRUE &&
            e->rightSite()->getLocation().distance(p) <
            closestEdge->rightSite()->getLocation().distance(p) )
         closestEdge = e;
/*
      if( e != rightEnd )
      {
         while(
                  e != rightEnd &&
                  e->atRight(p) == BOOL_TRUE &&
                  e->rightSite()->getLocation().distance(p) <
                  closestEdge->rightSite()->getLocation().distance(p) )
         {
            closestEdge = e;
            e = e->Next;
         }
      }
*/
    }
    else
    {
      tmpE = e;
      do{ e = e->Previous; } while(e != leftEnd && e->atRight(p) == BOOL_FALSE);

      // Regarde a gauche
      closestEdge = e;
      while(
                closestEdge->Previous != leftEnd &&
                closestEdge->Previous->atRight(p) == BOOL_TRUE &&
                closestEdge->Previous->rightSite()->getLocation().distance(p) <
                closestEdge->rightSite()->getLocation().distance(p) )
         closestEdge = closestEdge->Previous;

      // Regarde a droite
      if( tmpE != rightEnd )
         do
         {
            tmpE = tmpE->Next;
         }
         while( tmpE != rightEnd && tmpE->atRight(p) == BOOL_FALSE );

      if( tmpE != rightEnd && closestEdge != leftEnd )
      {
         if( tmpE != rightEnd &&
                     tmpE->atRight(p) == BOOL_TRUE &&
                     tmpE->rightSite()->getLocation().distance(p) <
                     closestEdge->rightSite()->getLocation().distance(p) )
            closestEdge = tmpE;
/*
         while(
                  tmpE != rightEnd &&
                  tmpE->atRight(p) == BOOL_TRUE &&
                  tmpE->rightSite()->getLocation().distance(p) <
                  closestEdge->rightSite()->getLocation().distance(p) )
         {
            closestEdge = tmpE;
            tmpE = tmpE->Next;
         }
*/
      }
    }

    // Mise a jour de la hash table et du nombre de references
    if( (tmpE = table_(buck)) != closestEdge )
    {
      if( tmpE !=  NULL)
         tmpE->removeReference();

      table_(buck) = closestEdge;
      closestEdge->addReference();
    }
    return closestEdge;
}

/******************************************************************************
 *
 * Nom       : intersection
 *
 ****************************************************************************/
vertex *halfEdge::intersection( halfEdge *e2 )
{
  // Si un des 2 edges n'est pas initialise. (Ex.: Sentinel)
  if(ve == NULL || e2->ve == NULL)
     return NULL;

  // Si les 2 halfEdges sont complementaire
  if(ve->getF2Idx() == e2->ve->getF2Idx() ) return NULL;

  point p;
  if( ve->intersection( *e2->ve, p ) == BOOL_FALSE ) return NULL;

  halfEdge *he = (halfEdge *)
  (point::cmp_yx(ve->getF2()->getLocation(),e2->ve->getF2()->getLocation())< 0)?
  this:e2;
  D( cout << "Site :" << he->ve->getF2()->getLocation() << endl; )

  // Verifie si l'intersection est a droite du site
  int right_of_site = p.xcoord() > he->ve->getF2()->getLocation().xcoord();



  // Verifie que l'intersection se trouve du bon cote
  if ((right_of_site &&  he->direction == 0) ||
        (!right_of_site && he->direction == 1)) return NULL;

  return new vertex(p);
}

Boolean halfEdge::atRight( const point &p ) const
{
  Boolean result, fast;
  int right_of_site;
  point topsite    = ve->getF2()->getLocation();

  right_of_site = p.xcoord() > topsite.xcoord();
  if( right_of_site && direction == 0) return BOOL_TRUE;
  if(!right_of_site && direction == 1) return BOOL_FALSE;

  point bottomsite = ve->getF1()->getLocation();
  segment s1( bottomsite, topsite );

  if( s1.is_vertical() || s1.slope() >= 1 )
  {
        double yl, t1, t2, t3;
        yl = c - (s1.dx()/s1.dy())*p.xcoord();
        t1 = p.ycoord() - yl;
        t2 = p.xcoord() - topsite.xcoord();
        t3 = yl - topsite.ycoord();
        result = (t1*t1 > t2*t2 + t3*t3)?BOOL_TRUE:BOOL_FALSE;
  }
  else
  {
        double dxp = p.xcoord() - topsite.xcoord();
        double dyp = p.ycoord() - topsite.ycoord();
        double m = s1.slope();
        fast = BOOL_FALSE;

      // Si le point est a gauche && la pente est negative
      // || le point est a droite && la pente est positive || nulle
        if( (!right_of_site && m<0.0) || (right_of_site && m >= 0.0) )
        {
          result = (dyp >= m*dxp)?BOOL_TRUE:BOOL_FALSE;
          fast   = result;
        }
        else
        {
          result = (p.xcoord() + p.ycoord()*m > c)?BOOL_TRUE:BOOL_FALSE;
          if(m<0.0) result = (result == BOOL_TRUE)?BOOL_FALSE:BOOL_TRUE;
          if( result == BOOL_FALSE ) fast = BOOL_TRUE;
        }
        if( fast == BOOL_FALSE )
        {
          double dxs = topsite.xcoord() - bottomsite.xcoord();
          result  = (m * (dxp*dxp - dyp*dyp) <
                         dxs*dyp*(1.0+2.0*dxp/dxs + m*m))?BOOL_TRUE:BOOL_FALSE;
          if(m<0.0) (result == BOOL_TRUE)?BOOL_FALSE:BOOL_TRUE;
        }
  }

  if( direction == 0 ) return result;
  else return (result==BOOL_TRUE)?BOOL_FALSE:BOOL_TRUE;
}

Home :: Fractals :: Tutorials :: Books :: Archive :: My blog :: My LinkedIn Profile :: Contact