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: chains.cpp
 * Titre    : Definition des classes de chaines.
 * Utilite  : Structure de donnee pour la location de point
 *            comme decrit dans 'Computational Geometry' de
 *            Franco P. Preparata & Michael Ian Shamos (P. 48)
 *
 * Note     : Les chaines monotones utilisees sur un diagramme de Voronoi
 *            permettent de trouver le plus proche voisin d'un point x en
 *            O(nlogn) comparativement a O(n^2) avec une methode brutale.
 *
 * Auteur   : Olivier Langlois <olivier@olivierlanglois.net>
 * Date     : 24 Mars 1998
 *
 * Revision :
 *
 * 001        27-Mar-1998 - Olivier Langlois
 *            Les fstreams ne fonctionnent pas avec le compilateur GNU 2.7.2.1
 *            lorsque compile avec l'option -frtti.
 */

#ifdef DEBUG
#ifndef __GNUC__
#include "diagwin.h"
#include <typeinfo.h>
#else
#include <std/typeinfo.h>
#include "diagdos.h"
#endif
#endif
#include "collect/colvar.h"
#include "collect/dict.h"
#include "../tools/collect/dynastck.cpp"
#include "../tools/collect/cirqueue.cpp"
#include "../tools/collect/bintree.cpp"
#include "../tools/collect/rbtree.cpp"
#include "../tools/collect/pq.cpp"
#include "geo/point.h"
#include "geo/dcel.h"
#include "geo/chains.h"
#ifdef __GNUG__
#include "umath.h"
#else
#include <math.h>
#endif
#include "mathc.h"

DEFINE_COLLECTABLE( ChainMesh, __CHAINMESH )
DEFINE_COLLECTABLE( Chain    , __CHAIN     )
DEFINE_COLLECTABLE( ChainsSet, __CHAINSSET )

DCEL *Chain::edgesPtr = NULL;

#ifdef __GNUG__
void *restorePtr;
void ChainMesh::saveGuts( FILE *os )
#else
void ChainMesh::saveGuts( ofstream &os )
#endif
{
#ifdef __GNUG__
  fwrite( &myAtom, sizeof(ClassID), 1, os );
  edgeIdx.saveGuts(os);
  fwrite( &myAtom, sizeof(ClassID), 1, os );
#else
  os.write((char *)&myAtom, sizeof(ClassID));
  edgeIdx.saveGuts(os);
  os.write((char *)&myAtom, sizeof(ClassID));
#endif
}
#ifdef __GNUG__
void ChainMesh::restoreGuts( FILE *is )
#else
void ChainMesh::restoreGuts( ifstream &is )
#endif
{
  ClassID id;

  // Verifie la validite du fichier
#ifdef __GNUG__
  fread( &id, sizeof(ClassID), 1, is );
#else
  is.read( (char *)&id, sizeof(ClassID) );
#endif
  if( id != myAtom ) throw int(0);

  edgeIdx.restoreGuts(is);

  // Verifie de nouveau l'intregrite du fichier
#ifdef __GNUG__
  fread( &id, sizeof(ClassID), 1, is );
#else
  is.read( (char *)&id, sizeof(ClassID) );
#endif
  if( id != myAtom ) throw int(0);
}

int ChainMesh::operator==(const ChainMesh &b) const
{
  vertex *v1  = edgesPtr->getEdge(&edgeIdx)->getV1();
  vertex *v2  = edgesPtr->getEdge(&edgeIdx)->getV2();
  vertex *bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV1();
  vertex *bv2 = b.edgesPtr->getEdge(&b.edgeIdx)->getV2();

  return ((*v1 == *bv1 && *v2 == *bv2) ||
             (*v1 == *bv2 && *v2 == *bv1)  );
}

int ChainMesh::operator<(const ChainMesh &b) const
{
  vertex *v1  = edgesPtr->getEdge(&edgeIdx)->getV1();
  vertex *bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV1();

  if( v1->isEqual(bv1) == BOOL_TRUE )
     bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV2();

  return( point::cmp_yx(*v1,*bv1)<0);
}

int ChainMesh::operator>(const ChainMesh &b) const
{
  vertex *v1  = edgesPtr->getEdge(&edgeIdx)->getV1();
  vertex *bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV1();

  if( v1->isEqual(bv1) == BOOL_TRUE )
     bv1 = b.edgesPtr->getEdge(&b.edgeIdx)->getV2();

  return( point::cmp_yx(*v1,*bv1)>0);
}

/******************************************************************************
 *
 * Nom       : isInMesh
 *
 * Utilite   : Verifie si le point est a la meme hauteur (ycoord) que le
 *             maillon.
 *
 * Reference : chains.h
 *
 ****************************************************************************/
int ChainMesh::isInMesh( const point &p )
{
  vertex *v1 = edgesPtr->getEdge(&edgeIdx)->getV1();
  vertex *v2 = edgesPtr->getEdge(&edgeIdx)->getV2();

  double resultV1 = point::cmp_yx(p, *v1);
  double resultV2 = point::cmp_yx(p, *v2);

  if( resultV1 != resultV2 ) return 0;
  else if( resultV1 < 0 && resultV2 < 0 ) return -1;
  else return 1;
}

/******************************************************************************
 *
 * Nom       : findMesh
 *
 * Utilite   : Trouve le maillon qui est a la meme hauteur que p.
 *
 * Reference : chains.h
 *
 ****************************************************************************/
CollectableULong *Chain::findMesh( const point &p )
{
  TreeNode< RBData<ChainMesh> > *n = Root;

  while(n != Sentinel)
  {
     int result = n->Data->RealData->isInMesh(p);

     if (result < 0)
        n = n->Less;
     else if ( result > 0 )
        n = n->More;
     else
        break;
  }

  if( n != Sentinel ) return n->Data->RealData->getEdgeIdx();
  else
  {
     return NULL;
  }
}
#ifdef __GNUG__
void Chain::saveGuts( FILE *os )
#else
void Chain::saveGuts( ofstream &os )
#endif
{
#ifdef __GNUG__
  fwrite( &myAtom, sizeof(ClassID), 1, os );
  fwrite( &chainNum, sizeof(unsigned long), 1, os );
  RedBlackTree<ChainMesh>::saveGuts(os);
  fwrite( &myAtom, sizeof(ClassID), 1, os );
#else
  os.write((char *)&myAtom, sizeof(ClassID));
  os.write((char *)&chainNum, sizeof(unsigned long));
  RedBlackTree<ChainMesh>::saveGuts(os);
  os.write((char *)&myAtom, sizeof(ClassID));
#endif
}
#ifdef __GNUG__
void Chain::restoreGuts( FILE *is )
#else
void Chain::restoreGuts( ifstream &is )
#endif
{
  ClassID id;

  // Verifie la validite du fichier
#ifdef __GNUG__
  fread( &id, sizeof(ClassID), 1, is );
#else
  is.read( (char *)&id, sizeof(ClassID) );
#endif
  if( id != myAtom ) throw int(0);
#ifdef __GNUG__
  fread( &chainNum, sizeof(unsigned long), 1, is );
#else
  is.read((char *)&chainNum, sizeof(unsigned long));
#endif
  RedBlackTree<ChainMesh>::restoreGuts(is);

  // Verifie de nouveau l'intregrite du fichier
#ifdef __GNUG__
  fread( &id, sizeof(ClassID), 1, is );
#else
  is.read( (char *)&id, sizeof(ClassID) );
#endif
  if( id != myAtom ) throw int(0);
}

void Chain::saveChainMesh( const ChainMesh &item, void *param )
{
#ifdef __GNUG__
  FILE *f = (FILE *)param;
#else
  ofstream &f = *(ofstream *)(ostream *)param;
#endif
  const_cast< ChainMesh & >(item).saveGuts( f );
}

void Chain::restoreChainMesh( const ChainMesh &item, void *param )
{
#ifdef __GNUG__
  FILE *f = (FILE *)param;
#else
  ifstream &f = *(ifstream *)(istream *)param;
#endif
  const_cast<ChainMesh &>(item).restoreGuts( f );
  const_cast<ChainMesh &>(item).edgesPtr = edgesPtr;
}

  ChainsSet::ChainsSet()
  : BinaryTree<Chain>(&ChainsSet::saveChain,&ChainsSet::restoreChain),
     edgesPtr(NULL) {}
  ChainsSet::ChainsSet( DCEL *list )
  : BinaryTree<Chain>(&ChainsSet::saveChain,&ChainsSet::restoreChain),
     edgesPtr(list) {init();}

// Classe utilisee par ChainsSet::init
class weight : public Collectable
{
  public:
  weight() : w(1), min(0), max(0) {}
  size_t w, min, max;
};

void ChainsSet::init()
{
  Dictionary weightTable;
  IDList     sortedVertices;
  CollectableULong ymin = edgesPtr->Vertices().getyminIdx();
  CollectableULong ymax = edgesPtr->Vertices().getymaxIdx();

  // Initialisation de weightTable
  DictIterator it(edgesPtr->Edges());
  while( ++it == BOOL_TRUE )
     weightTable.Insert( it.getKey()->copy(), new weight );

  DictIterator VerticesIt(edgesPtr->Vertices());

  while( ++VerticesIt == BOOL_TRUE )
  {
     // Insere toutes les vertices dans la liste sauf ymin & ymax
     if( VerticesIt.getKey()->isEqual(&ymin) == BOOL_FALSE &&
          VerticesIt.getKey()->isEqual(&ymax) == BOOL_FALSE  )

     // Ici, on suppose que les vertices ne se trouve pas dans une autre liste
     sortedVertices.insert((point *)VerticesIt.getValue());
  }
  sortedVertices.sort(sortSites);

  unsigned long Win, Wout;

  /*
   * d1 : Arete sortante la plus a gauche
   */
  edge *d1;
  segment sd1;

  // Premiere passe
  segment curSeg;
  IDListIterator svit(sortedVertices);
  vertex *curVertex;
  edge   *curEdge;
  CollectableULong *idx;
  D( while( (curVertex = (vertex *)++svit) ) cout << *curVertex << endl; )
  D( svit.reset(); )
  while( (curVertex = (vertex *)++svit) )
  {
     IDListIterator eit(*curVertex->getEdgeIdxList());
     Wout = 0;
     Win  = 0;
     d1 = NULL;
     while ( (idx = (CollectableULong *)++eit) )
     {
        curEdge = edgesPtr->getEdge(idx);
        vertex *secondVertex = curEdge->getV1();
        if( secondVertex->isEqual(curVertex) == BOOL_TRUE )
          curSeg = segment(*curVertex,*curEdge->getV2());
        else
          curSeg = segment(*curVertex,*secondVertex);

        D( cout << "Angle :" << curSeg.angle() << endl; )

        if( curSeg.angle() >= -0.00001 && curSeg.angle() < M_PI-0.00001 )
        {
          // Arete sortante
          if( !d1 || sd1.angle() < curSeg.angle() )
          {
             d1  = curEdge;
             sd1 = curSeg;
          }
          Wout += ((weight *)
                     weightTable.getValue(idx))->w;
        }
        else
        {
          // Arete rentrante
          Win  += ((weight *)
                     weightTable.getValue(idx))->w;
        }
     } // WHILE(edge)
     // A la place d'une exception, une arete de correction pourrait etre
     // cree.
     if( !Win || ! Wout ) throw GeoEx(GX_REGPSLG);
     if( Win > Wout )
     {
        ((weight *)
        weightTable.getValue(d1->getKey()))->w = Win-Wout+1;
     }
  }   // WHILE(vertices)

  // A la fin de la premiere passe, on peut connaitre le nombre de chaines.
  PQ<edgeInfo> edgePQ;
  edgeInfo     Sentinel;

  // Une clef valide ne peut pas etre negative.
  Sentinel.Key = -100;
#ifdef __BORLANDC__
  PQ<edgeInfo>::setSentinel(&Sentinel);
#else
  edgePQ.setSentinel(&Sentinel);
#endif
  edgeInfo    *curEdgeInfo;

  size_t numChains = 0;
  int curChain     = 0;
  curVertex = edgesPtr->getVertex(&ymax);
  IDListIterator ymaxEdgeit
                      (*curVertex->getEdgeIdxList());
  while( (idx = (CollectableULong *)++ymaxEdgeit) )
  {
     numChains += ((weight *)weightTable.getValue(idx))->w;

     curEdge = edgesPtr->getEdge(idx);
     vertex *secondVertex = curEdge->getV1();
     if( secondVertex->isEqual(curVertex) == BOOL_TRUE )
        curSeg = segment(*curVertex,*curEdge->getV2());
     else
        curSeg = segment(*curVertex,*secondVertex);

     double Angle = curSeg.angle();
     curEdgeInfo = new edgeInfo;
     // Si l'angle est PI, il est change pour -PI
     // Cela est necessaire pour l'ordre des elements de edgePQ
     curEdgeInfo->Key = ((Angle > 0)?-Angle:Angle);
     curEdgeInfo->e = curEdge;
     edgePQ.insert(curEdgeInfo);
  }
  D( cout << "numChains :" << numChains << endl; )
  for( unsigned long i = 0; i < numChains; i++ )
     Insert(new Chain(i));

  balance();

  while( (curEdgeInfo = edgePQ.removeFirst()) != &Sentinel )
  {
     weight *w = (weight *)weightTable.getValue(curEdgeInfo->e->getKey());
     w->min = curChain;
     w->max = curChain + w->w - 1;
     curChain += w->w;
     delete curEdgeInfo;
  }

  // Deuxieme passe

  while( (curVertex = (vertex *)--svit) )
  {
     IDListIterator eit(*curVertex->getEdgeIdxList());
     Wout     = 0;
     Win      = 0;
     curChain = numChains;

     while ( (idx = (CollectableULong *)++eit) )
     {
        curEdge = edgesPtr->getEdge(idx);
        vertex *secondVertex = curEdge->getV1();
        if( secondVertex->isEqual(curVertex) == BOOL_TRUE )
          curSeg = segment(*curVertex,*curEdge->getV2());
        else
          curSeg = segment(*curVertex,*secondVertex);

        D( cout << "Angle :" << curSeg.angle() << endl; )

        // La constante M_PI est plus precise que le nombre retourne par
        // segment::angle() lorsque la reponse est PI c'est pourquoi le test
        // == n'est pas utilise
        if( curSeg.angle() >= -0.00001 && curSeg.angle() < M_PI-0.00001 )
        {
          // Arete sortante
          weight *w = (weight *)
                          weightTable.getValue(idx);
          Wout += w->w;
          insert( new ChainMesh(curEdge->getKey(), edgesPtr), w->min, w->max );
          if( w->min < curChain )
             curChain = w->min;
        }
        else
        {
          // Arete rentrante
          double Angle = curSeg.angle();
          curEdgeInfo = new edgeInfo;
          // Si l'angle est PI, il est change pour -PI
          // Cela est necessaire pour l'ordre des elements de edgePQ
          curEdgeInfo->Key = ((Angle > 0)?-Angle:Angle);
          curEdgeInfo->e = curEdge;
          edgePQ.insert(curEdgeInfo);

          Win  += ((weight *)
                     weightTable.getValue(idx))->w;
        }
     } // WHILE(edge)

     if( Win < Wout )
     {
        weight *w = (weight *)weightTable.getValue(edgePQ.first()->e->getKey());
        D( cout << "Win " << Win << " Wout " << Wout << " Wd2 " << w->w << endl; )
        w->w = Wout-Win+w->w;
     }

     // Assigne les min/max aretes entrante
     while( (curEdgeInfo = edgePQ.removeFirst()) != &Sentinel )
     {
#ifdef __WIN32__
        D( if( curChain < 0 || curChain > numChains) DiagOutWin("Erreur").DisplayMsg("curChain Corrompu", DIAG_WARNING); )
#endif
        D( cout << "Segment :" << curEdgeInfo->e->getSegment() << endl; )
        weight *w = (weight *)weightTable.getValue(curEdgeInfo->e->getKey());
        w->min = curChain;
        w->max = curChain + w->w - 1;
        curChain += w->w;
        delete curEdgeInfo;
     }
  }   // WHILE(vertices)

  // Insere les aretes connectees a ymin.
  curVertex = edgesPtr->getVertex(&ymin);
  IDListIterator eit(*curVertex->getEdgeIdxList());
  while ( (idx = (CollectableULong *)++eit) )
  {
     weight *w = (weight *)
                     weightTable.getValue(idx);
     insert( new ChainMesh(idx, edgesPtr), w->min, w->max );
  }

  // Vide la liste de ses vertices pour eviter qu'ils ne soient effaces
  sortedVertices.clearWithoutDestroy();
}

/******************************************************************************
 *
 * Nom       : insert
 *
 * Utilite   : Insere une maille.
 *
 * Reference : chains.h
 *
 ****************************************************************************/
void ChainsSet::insert( ChainMesh *m, unsigned long min, unsigned long  max)
{
  TreeNode<Chain> *n = Root;
  D( cout << "min :" << min << " max :" << max << endl; )

  while(n != Sentinel &&
          (n->Data->getChainNum() < min  || n->Data->getChainNum() > max) )
  {
     D( cout << "chainNum :" << n->Data->getChainNum() << endl; )
     if (max < n->Data->getChainNum())
        n = n->Less;
     else
        n = n->More;
  }

  if( n != Sentinel ) n->Data->Insert(m);
  else
  {
     D( cout << "Chaine non trouvee" << endl; )
     delete m;
  }
}

/******************************************************************************
 *
 * Nom       : location
 *
 * Utilite   : Trouve la location de p.
 *
 * Reference : chains.h
 *
 ****************************************************************************/
face *ChainsSet::location( const point &p )
{
  TreeNode<Chain> *n = Root;
  edge *e = NULL;
  vertex *v1;
  vertex *v2;
  Boolean atRight;
  face   *loc;

  while(1)
  {
     CollectableULong *result = n->Data->findMesh(p);
     if( result )
     {
        e = edgesPtr->getEdge(result);
        v1 = e->getV1();
        v2 = e->getV2();

        if( point::cmp_yx( *v1, *v2 ) > 0 )
        {
          vertex *tmp = v1;
          v1 = v2;
          v2 = tmp;
        }

        if( right_turn( *v1, *v2, p ) == BOOL_TRUE )
          atRight = BOOL_TRUE;
        // Le point est contenu dans le segment, Le point n'est dans aucune
        // face mais les 2 faces les plus proches sont e.F1 et e.F2 et puisque
        // dans la methode DCEL::setBox la face outside est toujour F1 et dans
        // la majorite des applications elle moins utile que les autres faces,
        // F2 est retournee dans cette situation.
        else if( segment( *v1, *v2 ).contains(p) == BOOL_TRUE )
          return e->getF2();
        else
          atRight = BOOL_FALSE;

        if ( atRight == BOOL_FALSE && n->Less != Sentinel )
          n = n->Less;
        else if ( atRight == BOOL_TRUE && n->More != Sentinel )
          n = n->More;
        else
        {
          face *f1 = e->getF1();
          face *f2 = e->getF2();
          if( f1->getLocation().distance(p) < f2->getLocation().distance(p) )
             loc = f1;
          else
             loc = f2;
          break;
        }
     } // IF(result)
     // Sinon continue a chercher l'arbre primaire
     else if( e )
     {
        if ( atRight == BOOL_FALSE && n->Less != Sentinel )
          n = n->Less;
        else if ( atRight == BOOL_TRUE && n->More != Sentinel )
          n = n->More;
        else
        {
          face *f1 = e->getF1();
          face *f2 = e->getF2();
          if( f1->getLocation().distance(p) < f2->getLocation().distance(p) )
             loc = f1;
          else
             loc = f2;
          break;
        }
     }
     // Pour eviter un loop infini mais cela ne devrait pas arriver.
     else return NULL;
  }

  return loc;
}
#ifdef __GNUG__
void ChainsSet::saveGuts( FILE *os )
#else
void ChainsSet::saveGuts( ofstream &os )
#endif
{
#ifdef __GNUG__
  fwrite( &myAtom, sizeof(ClassID), 1, os );
  edgesPtr->saveGuts(os);
  BinaryTree<Chain>::saveGuts(os);
  fwrite( &myAtom, sizeof(ClassID), 1, os );
#else
  os.write((char *)&myAtom, sizeof(ClassID));
  edgesPtr->saveGuts(os);
  BinaryTree<Chain>::saveGuts(os);
  os.write((char *)&myAtom, sizeof(ClassID));
#endif
}
#ifdef __GNUG__
void ChainsSet::restoreGuts( FILE *is )
#else
void ChainsSet::restoreGuts( ifstream &is )
#endif
{
  ClassID id;

  // Verifie la validite du fichier
#ifdef __GNUG__
  fread( &id, sizeof(ClassID), 1, is );
#else
  is.read( (char *)&id, sizeof(ClassID) );
#endif
  if( id != myAtom ) throw TreeEx(BTX_FILE);

  delete edgesPtr;
  edgesPtr = new DCEL;
  edgesPtr->restoreGuts(is);
  Chain::edgesPtr = edgesPtr;
  BinaryTree<Chain>::restoreGuts(is);

  // Verifie de nouveau l'intregrite du fichier
#ifdef __GNUG__
  fread( &id, sizeof(ClassID), 1, is );
#else
  is.read( (char *)&id, sizeof(ClassID) );
#endif
  if( id != myAtom ) throw TreeEx(BTX_FILE);
}

void ChainsSet::saveChain( const Chain &item, void *param )
{
#ifdef __GNUG__
  FILE *f = (FILE *)param;
#else
  ofstream &f = *(ofstream *)(ostream *)param;
#endif
  const_cast< Chain & >(item).saveGuts( f );
}

void ChainsSet::restoreChain( const Chain &item, void *param )
{
#ifdef __GNUG__
  FILE *f = (FILE *)param;
#else
  ifstream &f = *(ifstream *)(istream *)param;
#endif
  const_cast< Chain & >(item).restoreGuts( f );
}

#ifdef __GNUG__
template class PQ<edgeInfo>;
template class BinaryTree<Chain>;
template class BinaryTree<RBData<ChainMesh> >;
template class RBData<ChainMesh>;
template class RedBlackTree<ChainMesh>;
template class Queue<TreeNode<Chain> *>;
template class cirQueue<TreeNode<Chain> *>;
template class Queue<TreeNode<RBData<ChainMesh> > *>;
template class cirQueue<TreeNode<RBData<ChainMesh> > *>;
template class TreeNode< Chain >;
template class TreeNode< RBData<ChainMesh> >;
typedef TreeNode<Chain>* nodeChainPtr;
typedef TreeNode<RBData<ChainMesh> >* nodeChainMeshPtr;
template void saveNodeToFile( const nodeChainPtr &,void * );
template void saveNodeToFile( const nodeChainMeshPtr &,void * );
template void saveBTNodeToFile(
          const Chain &item,
          void *param
          );
template void saveBTNodeToFile(
          const RBData<ChainMesh> &item,
          void *param
          );
template void saveRBTDataToFile( const ChainMesh &, void * );
template void saveRBTNodeToFile( const RBData<ChainMesh> &item, void *param );
template void restoreRBTNodeFromFile( const RBData<ChainMesh> &item,
                      void *param );
template void ProcessData( const RBData<ChainMesh> &item, void *param );

#endif

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