/*
 * Module ID: voronoi.cpp
 * Titre    : Construction d'un diagrame de Voronoi.
 * Utilite  : Utilise l'algorithme de Steven Fortune des laboratoires
 *            Bell AT&T comme decrit dans 'Computational Geometry in C' de
 *            J. O'Rourke. (P. 179)
 *
 * Limite   : Cet algorithme n'est pas le plus approprié lorsque l'on
 *            desire ajouter ou supprimer des sites de l'ensemble de sites
 *            de depart.
 *            De plus, avec certaines conditions (inconnues) le diagramme
 *            resultant de cet algorithme est errone.
 *
 * Auteur   : Olivier Langlois (olivier@olivierlanglois.net)
 * Date     : 07 Mars 1998
 */

#ifdef __BORLANDC__
// Indique au compilateur que les instances de template utilisees dans ce
// module seront definies dans un autre module.
#pragma option -Jgx
#endif

#ifdef _MSC_VER
#define __WIN32__ _WIN32
#endif

#ifdef __GNUG__
#include <std/typeinfo.h>
#endif
#include "collect/idlist.h"
#include "geo/point.h"
#include "geo/segment.h"
#include "geo/circle.h"
#include "../tools/collect/pq.cpp"
#include "collect/colvar.h"
#include "collect/dict.h"
#include "geo/dcel.h"
#include "collect/vector.h"
#include "geo/hedge.h"
#include "geo/voronoi.h"
#ifdef __GNUG__
#include "umath.h"
#else
#include <math.h>
#endif
#include "debug.h"
#ifdef __WIN32__
#include "diagwin.h"
#else
#include "diagdos.h"
#endif

Boolean yStarIsSmallest( circle &star, IDList &sl)
{
  IDListIterator it(sl);
  point *p;
  while( ( p = (point *)++it) &&
              star.inside(*p) == BOOL_FALSE)
     if( p->ycoord() > star.center().ycoord()+star.radius()  ||
         (p->ycoord() == star.center().ycoord()+star.radius() &&
          p->xcoord() >  star.center().xcoord()) ) return BOOL_TRUE;
  return BOOL_FALSE;
}

/******************************************************************************
 *
 * Nom       : Voronoi
 *
 * Utilite   : Contruire un diagrame de Voronoi.
 *
 * Reference : voronoi.h
 *
 ****************************************************************************/
DCEL *Voronoi( IDList &sl )
{
  face *newsite, *bot, *top, *bottomsite, *temp;
  vertex *v, *p;
  point *tmpPoint;
  circle newintstar;
  halfEdge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
  edge *e;
  line  l;
  char  direction;
  float xmax = 0;

  // Prepare un 'plane sweep' en triant les sites
  sl.sort(sortSites);

  IDListIterator it(sl);

  while( tmpPoint = (point *)++it )
  {
     D( cout << *tmpPoint << endl; )
     if( tmpPoint->xcoord() < point::xmin ) point::xmin = tmpPoint->xcoord();
     if( tmpPoint->xcoord() > xmax        ) xmax        = tmpPoint->xcoord();
  }
  if( xmax - point::xmin )
     point::deltax = (xmax - point::xmin);

  DCEL *edgeList = new DCEL;
  halfEdge::init( sl.entries() );

  PQ<halfEdge> pQueue;
  // Assigne un Sentinel a la PQ.
  lbnd = new halfEdge;
  lbnd->setYStar( new vertex(point( -HUGE_VAL, -HUGE_VAL )),
                             point( -HUGE_VAL, -HUGE_VAL ));
  pQueue.setSentinel( lbnd );

  tmpPoint = (point *)sl.get();
  if(tmpPoint)
  {
     bottomsite = new face(*tmpPoint);
     delete tmpPoint;
  }
  else bottomsite = NULL;

  tmpPoint = (point *)sl.get();
  if(tmpPoint)
  {
     newsite    = new face(*tmpPoint);
     delete tmpPoint;
  }
  else newsite = NULL;

try{
  while(1)
  {
     if(pQueue.isEmpty() == BOOL_FALSE)
     {
        newintstar = pQueue.first()->getYStar();
     }

     if( newsite != NULL
          && (pQueue.isEmpty() == BOOL_TRUE
           /*
            * Chaque vertex du diagramme est le centre d'un cercle dont
            * ses 3 sites sont cocirculaires.
            * newintstar sert a determiner si le nouveau site est a
            * l'interieur du rayon du cercle des 3 derniers sites traites
            */
          || newintstar.inside(newsite->getLocation()) == BOOL_TRUE
          || yStarIsSmallest( newintstar, sl) == BOOL_FALSE ))

     {
        /* Le nouveau site est le plus petit */
        D( if( newsite ) cout << " newsite :" << newsite->getLocation() << endl; )
        lbnd = halfEdge::leftbnd(newsite->getLocation());
/*
        if( sl.first() &&
             lbnd->ve   &&
             lbnd->rightSite()->getLocation().distance(*(point *)sl.first()) <
             lbnd->rightSite()->getLocation().distance(newsite->getLocation()) &&
             lbnd->next() != halfEdge::rightEnd &&
             lbnd->next()->atRight(newsite->getLocation()) == BOOL_TRUE )
          lbnd = lbnd->next();
*/
        rbnd = lbnd->next();
//      D( if( lbnd == halfEdge::leftEnd || lbnd == halfEdge::rightEnd ) DiagOutWin("Voronoi").DisplayMsg("Mauvais lbnd"); )

        if( !lbnd->ve ) bot = bottomsite;
        else bot = lbnd->rightSite();

        l = p_bisector(bot->getLocation(), newsite->getLocation());
        D( cout << "Bissectrice :" << l << endl;                           )
        D( if( lbnd->ve ) cout << "lbnd :" << lbnd->ve->getSegment() << ((lbnd->direction)?"right":"left") << lbnd->ve->getF1()->getLocation() << lbnd->ve->getF2()->getLocation() << endl; )
        D( if( rbnd->ve ) cout << "rbnd :" << rbnd->ve->getSegment() << ((rbnd->direction)?"right":"left") << rbnd->ve->getF1()->getLocation() << rbnd->ve->getF2()->getLocation() << endl; )

        e = new edge( edgeList, bot, newsite, l );
        bisector = new halfEdge( e, 0 );
        bisector->insert( lbnd );
        if ((p = lbnd->intersection(bisector)) != NULL)
        {
          pQueue.remove(lbnd);
          lbnd->removeYStar();
          lbnd->setYStar(p,newsite->getLocation());
          pQueue.insert(lbnd);
        }
        lbnd = bisector;
        bisector = new halfEdge( e, 1 );
        bisector->insert(lbnd);
        if ((p = bisector->intersection(rbnd)) != NULL)
        {
          bisector->setYStar(p,newsite->getLocation());
          pQueue.insert(bisector);
        }
        tmpPoint = (point *)sl.get();
        if(tmpPoint)
        {
          newsite    = new face(*tmpPoint);
          delete tmpPoint;
        }
        else
        {
          newsite = NULL;
        }
     }
     else if ( pQueue.isEmpty() == BOOL_FALSE )
     /* l'intersection est la plus petite */
     {
        lbnd  = pQueue.removeFirst();
        llbnd = lbnd->previous();
        rbnd  = lbnd->next();
        rrbnd = rbnd->next();
        bot   = lbnd->leftSite();
        top   = rbnd->rightSite();
        v     = lbnd->Vertex;

        D( cout << "Vertex :" << *(point *)v << endl; )
        D( if( lbnd->ve ) cout << "lbnd :" << lbnd->ve->getSegment() << ((lbnd->direction)?"right":"left") << lbnd->ve->getF1()->getLocation() << lbnd->ve->getF2()->getLocation() << endl; )
        D( if( rbnd->ve ) cout << "rbnd :" << rbnd->ve->getSegment() << ((rbnd->direction)?"right":"left") << rbnd->ve->getF1()->getLocation() << rbnd->ve->getF2()->getLocation() << endl; )
        D( if( llbnd->ve) cout << "llbnd:" << lbnd->ve->getSegment() << ((llbnd->direction)?"right":"left") << llbnd->ve->getF1()->getLocation() << llbnd->ve->getF2()->getLocation() << endl; )
        D( if( rrbnd->ve) cout << "rrbnd:" << rbnd->ve->getSegment() << ((rrbnd->direction)?"right":"left") << rrbnd->ve->getF1()->getLocation() << rrbnd->ve->getF2()->getLocation() << endl; )
      D( if( pQueue.first() ) cout << "Next intersection: " << pQueue.first()->ve->getSegment() << endl; )

        D( if(line(lbnd->ve->getSegment()).contains(*v) == BOOL_FALSE) cout << "Mauvais vertex" << endl; )
        if( lbnd->direction == 0 )
        // 1
          lbnd->setV1(v);
        else
          lbnd->setV2(v);


        D( if(line(rbnd->ve->getSegment()).contains(*v) == BOOL_FALSE) cout << "Mauvais vertex" << endl; )
        if( rbnd->direction == 0 )
        // 2
          rbnd->setV1(v);
        else
          rbnd->setV2(v);

        D( cout << "Site 1:" << bot->getLocation() << " Site 2:" << top->getLocation() << endl; )

        if( bot->getLocation().ycoord() > top->getLocation().ycoord() )
        { temp = bot; bot = top; top = temp; direction = 1; }
        else direction = 0;

        l = p_bisector(bot->getLocation(), top->getLocation());
        e = new edge( edgeList, bot, top, l );
        bisector = new halfEdge(e, direction);
        D( cout << "Bissectrice :" << l << endl; )
        /*
         * Chaque vertex du diagramme de Voronoi possede exactement 3
         * aretes. ( Computational Geometry, [PREPARATA,SHAMOS] P.206 )
         */
        // 3
        if( direction == 1 )
          bisector->setV1(v);
        else
          bisector->setV2(v);

        D( if(l.contains(*v) == BOOL_FALSE) cout << "Mauvais vertex" << endl; )
        D( if(lbnd == halfEdge::leftEnd) cout << "lbnd == halfEdge::leftEnd" << endl; )
        D( if(rbnd == halfEdge::rightEnd) cout << "rbnd == halfEdge::rightEnd" << endl; )

        // Insere le edge dans le DCEL
        delete lbnd;
        pQueue.remove(rbnd);
        delete rbnd;

        bisector->insert(llbnd);

        if ((p = llbnd->intersection(bisector)) != NULL)
        {
          pQueue.remove(llbnd);
          llbnd->removeYStar();
          llbnd->setYStar(p,bot->getLocation());
          pQueue.insert(llbnd);
        }
        if ((p = bisector->intersection(rrbnd)) != NULL)
        {
          bisector->setYStar(p,bot->getLocation());
          pQueue.insert(bisector);
        }
     }
     else break;
  }
} catch( GeoEx &ex )
{
#ifdef __WIN32__
  ex.Explain( DiagOutWin("Erreur") );
#else
  DiagOutDOS out(cerr);
  ex.Explain( out );
#endif
}
  // Les halfEdges restant sont les edges non connectees du diagramme.
/*
  for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd))
  {
     e = lbnd -> ELedge;
     out_ep(e);
  }
*/
  pQueue.clear();
  delete pQueue.getSentinel();
  halfEdge::terminate();

  // Sanity Check
  /*
   * Un diagramme de Voronoi de N sites possede au plus 2N - 5 vertices
   * et 3N - 6 aretes.
   * ( Computational Geometry, [PREPARATA,SHAMOS] P.211 )
   */
  D( size_t N = edgeList->Faces().entries(); )
  D( cout << "Sites :" << N << " Vertices :" << edgeList->Vertices().entries() << " Aretes :" << edgeList->Edges().entries() << endl; )

  // Ne tient pas compte des vertices non attaches au diagrame
  // Cela peut parfois donner une fausse alarme.
  D( if( edgeList->Vertices().entries() > 2*N - 5 ) cout << "Voronoi: N = " << N << " numVertices = " << edgeList->Vertices().entries() << endl; )
  D( if( edgeList->Edges().entries() > 3*N - 6 ) cout << "Voronoi: N = " << N << " numEdges = " << edgeList->Edges().entries() << endl; )

  return edgeList;
}

#ifdef __GNUG__
template class PQ<halfEdge>;
#endif