include/shape/voronoi2dimpl.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * Voronoi2d
 * 2011-04-05 
* Angelos Mantzaflaris
 *****************************************************************************
 *               Copyright (C) 2010 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
# ifndef shape_voronoi2dimpl_hpp
# define shape_voronoi2dimpl_hpp

# include <shape/topology.hpp>
# include <shape/voronoi_site2d.hpp>
# include <shape/bcell2d_voronoi_impl2d.hpp>
# include <shape/bcell2d_voronoi_impl_diagram.hpp>
# include <shape/tpl3d.hpp>

# include <shape/semialgebraic_domain2d_axel.hpp>
# include <stack>

# define  STACK std::stack<Cell2d *>

# define TMPL template<class C, class V>
# define Shape geometric<V>
# define VoronoiCurve voronoi_bisector<C,V>
# define VoronoiSite2d voronoi_site2d<C,V >
# define VSite voronoi_site2d<C,V >
# define Cell2dVoronoiImpl2d bcell2d_voronoi_impl2d<C,V>
# define SELF voronoi2dimpl<C,V>
# define Viewer viewer<V>
# define VCell bcell2d_voronoi_impl_diagram<C,V>
# undef Cell
//====================================================================
namespace mmx { namespace shape {
//====================================================================

template<class C, class V=default_env>
class voronoi2dimpl
{

public:

  typedef topology<C,REF_OF(V)>           Topology;

  typedef typename Topology::Point        Point;  
  typedef typename Topology::Edge         Edge;
  typedef typename Topology::Face         Face;
  typedef typename Topology::BoundingBox BoundingBox;
  typedef typename Topology::Cell         Cell;
  typedef typename mesher2d<C,V>::Cell2d       Cell2d;
  typedef typename mesher2d<C,V>::Curve        Curve;
  
  typedef node<Cell*>               Node;
  typedef kdtree<Cell*>             Kdtree;

public:
  
  voronoi2dimpl(const BoundingBox& bx): m_maxprec(0.1), m_minprec(0.01), m_bbx(bx) { };

  voronoi2dimpl(Curve * curve) ;
  ~voronoi2dimpl(void) {};//   delete this->m_tree ;  };

  void run(void)  ;
  void push_back(VoronoiSite2d *v);

  unsigned const size() { return m_sites.size(); };

  void set_precision (double);
  void set_smoothness(double);

protected:
  bool is_regular (Cell * bcell) 
    { 
      return bcell->is_regular();
      //return mesher2d<C,V>::is_regular(bcell); 
    };

  bool insert_regular(Cell * bcell) 
    { 
      //return bcell->insert_regular(this);

  Seq<Point*> l;
    l= ((VCell*)bcell)->intersections();

  int * sz;
  int * st;
  Point *q;

  if ( l.size()==2 )
  {
    this->m_graph.push_vertex(l[0]);
    this->m_graph.push_vertex(l[1]);
    this->m_graph.push_uedge(l[0],l[1]);

    return true;
  }

  if ( l.size()==4 ) // two dublicated
  {
    this->m_graph.push_vertex(l[0]);
    this->m_graph.push_vertex(l[1]);
    this->m_graph.push_uedge(l[0],l[1]);
    return true;
  }
  
  if ( l.size()==1)
  {
    std::cout<< "SIZE ONE, "<< *bcell<<std::endl;

    this->m_graph.push_vertex(l[0]);

    foreach( Cell3d* c, ((VCell*)bcell)->m_objects )
      if ( ((Cell2d*)c)->nb_intersect()==1 )
      {
        sz= ((Cell2dVoronoiImpl2d*)c)->m_polynomial.rep().szs();
        st= ((Cell2dVoronoiImpl2d*)c)->m_polynomial.rep().str();
        if ( ((Cell2dVoronoiImpl2d*)c)->m_polynomial[0]<EPSILON )
        {
          q= new Point(bcell->xmin(),bcell->ymin(),0); 
//          std::cout<< "1.add ("<< *q <<")->("<< *l[0] << ") in "<< *this<<std::endl;

          this->m_graph.push_vertex(q);
          this->m_graph.push_uedge(l[0],q);

          ((Cell2dVoronoiImpl2d*)c)->n_intersections<< q;

          foreach( Cell2d *nb, ((Cell2d*)bcell)->s_neighbors )
            foreach( Cell3d* cc, ((VCell*)nb)->m_objects )
            if ( cc->is_active() )
            {
              ((Cell2dVoronoiImpl2d*)cc)->n_intersections<< q;
              std::cout<<"This Intersections: "<< ((VCell*)bcell)->intersections().size() << std::endl;
              std::cout<<"Neib Intersections: "<< nb->intersections().size() << std::endl;
              return true;
            }
          foreach( Cell2d *nb, ((Cell2d*)bcell)->w_neighbors )
            foreach( Cell3d* cc, ((VCell*)nb)->m_objects )
            if ( cc->is_active() )
            {
              ((Cell2dVoronoiImpl2d*)cc)->e_intersections<< q;
              std::cout<<"This Intersections: "<< ((VCell*)bcell)->intersections().size() << std::endl;
              std::cout<<"Neib Intersections: "<< nb->intersections().size() << std::endl;
              return true;
            }

        } else if ( ((Cell2dVoronoiImpl2d*)c)->m_polynomial[(sz[0]-1)*st[0]]<EPSILON )
        {
          q= new Point(bcell->xmax(),bcell->ymax(),0);
          std::cout<< "2.add ("<< *q <<")->("<< *l[0] << ") in "<< *bcell<<std::endl;

          this->m_graph.push_vertex(q);
          this->m_graph.push_uedge(l[0],q);
          ((Cell2dVoronoiImpl2d*)c)->s_intersections<< q;


          foreach( Cell2d *nb, ((VCell*)bcell)->e_neighbors )
            foreach( Cell3d* cc, ((VCell*)nb)->m_objects )
            if ( cc->is_active() )
            {
              ((Cell2dVoronoiImpl2d*)cc)->w_intersections<< q;
              std::cout<<"This Intersections: "<< ((VCell*)bcell)->intersections().size() << std::endl;
              std::cout<<"Neib Intersections: "<< nb->intersections().size() << std::endl;
              return true;
            }
          foreach( Cell2d *nb, ((VCell*)bcell)->n_neighbors )
            foreach( Cell3d* cc, ((VCell*)nb)->m_objects )
            if ( cc->is_active() )
            {
              ((Cell2dVoronoiImpl2d*)cc)->s_intersections<< q;
              std::cout<<"This Intersections: "<< ((VCell*)bcell)->intersections().size() << std::endl;
              std::cout<<"Neib Intersections: "<< nb->intersections().size() << std::endl;
              return true;
            }

        } else if ( ((Cell2dVoronoiImpl2d*)c)->m_polynomial[sz[0]*sz[1]-1]<EPSILON )
        {
          q= new Point(bcell->xmin(),bcell->ymax(),0);
//          std::cout<< "3.add ("<< *q <<")->("<< *l[0] << ") in "<< *this<<std::endl;
          this->m_graph.push_vertex(q);
          this->m_graph.push_uedge(l[0],q);

          ((Cell2dVoronoiImpl2d*)c)->w_intersections<< q;
          return true;
        } else if ( ((Cell2dVoronoiImpl2d*)c)->m_polynomial[(sz[1]-1)*st[1]]<EPSILON )
        {
          q= new Point(bcell->xmax(),bcell->ymin(),0);
//          std::cout<< "4.add ("<< *q <<")->("<< *l[0] << ") in "<< *this<<std::endl;
          this->m_graph.push_vertex(q);
          this->m_graph.push_uedge(l[0],q);
          ((Cell2dVoronoiImpl2d*)c)->e_intersections<< q;
          return true;
        }  
      }
  }

  if ( l.size()==0)
  {
    std::cout<< "SIZE ZERO, "<< *bcell<<std::endl;
    return true;

    foreach( Cell3d* c, ((VCell*)bcell)->m_objects )
      {
        sz= ((Cell2dVoronoiImpl2d*)c)->m_polynomial.rep().szs();
        st= ((Cell2dVoronoiImpl2d*)c)->m_polynomial.rep().str();
        if ( ((Cell2dVoronoiImpl2d*)c)->m_polynomial[0]<EPSILON )
        {
          q= new Point(bcell->xmin(),bcell->ymin(),0); 
          std::cout<< "1.add ("<< *q <<") in "<< *bcell<<std::endl;
          this->m_graph.push_vertex(q);
          ((Cell2dVoronoiImpl2d*)c)->n_intersections<< q;
          return true;
        } else if ( ((Cell2dVoronoiImpl2d*)c)->m_polynomial[(sz[0]-1)*st[0]]<EPSILON )
        {
          q= new Point(bcell->xmax(),bcell->ymax(),0);
          std::cout<< "1.add ("<< *q <<") in "<< *bcell<<std::endl;
          this->m_graph.push_vertex(q);
          ((Cell2dVoronoiImpl2d*)c)->s_intersections<< q;
          return true;
        } else if ( ((Cell2dVoronoiImpl2d*)c)->m_polynomial[sz[0]*sz[1]-1]<EPSILON )
        {
          q= new Point(bcell->xmin(),bcell->ymax(),0);
          std::cout<< "1.add ("<< *q <<") in "<< *bcell<<std::endl;
          this->m_graph.push_uedge(l[0],q);
          ((Cell2dVoronoiImpl2d*)c)->w_intersections<< q;
          return true;
        } else if ( ((Cell2dVoronoiImpl2d*)c)->m_polynomial[(sz[1]-1)*st[1]]<EPSILON )
        {
          q= new Point(bcell->xmax(),bcell->ymin(),0);
          std::cout<< "1.add ("<< *q <<") in "<< *bcell<<std::endl;
          this->m_graph.push_vertex(q);
          this->m_graph.push_uedge(l[0],q);
          ((Cell2dVoronoiImpl2d*)c)->e_intersections<< q;
          return true;
        }  
      }
  }

//    std::cout<< "nb_in= "<<  ((VCell*)bcell)->nb_intersect() <<std::endl;
    //std::cout<< "box  = "<< *this <<std::endl;
        //foreach(Point* q, this->intersections() ) 
        //std::cout<< " "<< *q <<std::endl;
//    print((Cell2dAlgebraicCurve*)m_objects[0]);
    return true;

    };


//   bool singularity(Cell * bcell) 
//     { 
//       std::cout<< "Voronoi Singularity"<<std::endl;
// return mesher2d<C,V>::singularity(bcell); };

  bool subdivide(Cell * bcell) 
    { return mesher2d<C,V>::subdivide(bcell,NULL); };

private:

public:

  std::list<VoronoiSite2d *> m_sites ;

  std::list<VCell *> m_singular_cells;  

  Graph<Point*>       m_graph;  // topological graph
  
  Graph<Cell2d*>      m_leaves; // active leaves graph
  Graph<Cell2d*>      b_leaves; // border leaves CCW loop

  double            m_maxprec ;
  double            m_minprec ;
  
  BoundingBox       m_bbx;
  Seq<VSite*>       m_objects;

  Seq<Face *>         m_faces;
};
  
//--------------------------------------------------------------------

TMPL void 
SELF::push_back(VoronoiSite2d *v)
{
  this->m_objects.push_back(v); 
  this->m_sites.push_back(v);

  //std::cout<<"added "<< v->coordinates() << std::endl;
}

TMPL void SELF::set_smoothness(double eps) {
  m_maxprec = eps;
}

TMPL void SELF::set_precision(double eps) {
  m_minprec = eps;
}

TMPL void 
SELF::run() {
     // assume this->m_sites.size()>1
     //mesher2d<C,V>::run();
     //std::cout<<"Topology run OK." <<std::endl;


  std::cout<<"Running, objs="<<this->m_objects.size() << std::endl;

  //Cell* cl = Cell2dFactory::instance()->create(this->m_objects,this->m_bbx);
  VCell * l = new VCell(this->m_bbx);
  BoundingBox bx= this->m_bbx;
  foreach( VSite * vcurve, this->m_objects) 
  {
    //VoronoiSite2d * vcurve = dynamic_cast<VoronoiSite2d *>(o);
    this->m_bbx.set_zmin(0);
    bx.set_zmax(1);

    l->push_back( new Cell2dVoronoiImpl2d(vcurve,this->m_bbx) );
  }

//   double maxsz = this->m_maxprec* this->size();
//   double minsz = this->m_minprec* this->size();
  double maxsz = this->m_maxprec;
  double minsz = this->m_minprec;

  
  std::stack<Cell *> bcell_stack;
  bcell_stack.push(l);

  while(!bcell_stack.empty()  ) {   //  

//    Node * node = this->m_nodes.front() ;
    Cell   * cl = bcell_stack.top();
    bcell_stack.pop();
    Cell2d * cl2= dynamic_cast<Cell2d*>(cl);

    //std::cout<<"node "<<cl <<std::endl;

    if(cl->is_active()) {

      //std::cout<<"active.. "<<std::endl;

      if(cl->size() > maxsz) 
        { 
          //this->subdivide(cl, node) ;
          Cell * left, * right;
            ((bcell<C,V> *)cl)->subdivide(left, right) ;
          bcell_stack.push(left);
          bcell_stack.push(right);
        }
      else if(this->is_regular(cl) )// && cl->is_intersected())
      { 
        //std::cout<<"ins regular "<<std::endl;
        this->insert_regular(cl) ;
        //std::cout<<"ok "<<std::endl;
        cl2->m_gnode= 
          this->m_leaves.push_vertex(cl2);
        if (  cl2->is_border() )
          this->b_leaves.push_vertex(cl2);       
      }
      else if(cl->size() > minsz)
      { 
        //std::cout<<"subdivide "<<std::endl;
          Cell* left, * right;
          cl->subdivide(left, right) ;
          bcell_stack.push(left);
          bcell_stack.push(right);
      }
      else
      { 
        //std::cout<<"singularity "<<std::endl;
        //this->singularity(cl) ;
//        std::cout<<"VD, Inserting singular"<<*cl2<<std::endl;  
        this->m_singular_cells<< (VCell*)cl2; // does not work with forwa
        ( (VCell*)cl2 )->process_singular() ;

        cl2->m_gnode= 
          this->m_leaves.push_vertex(cl2);
        if (  cl2->is_border() )
          this->b_leaves.push_vertex(cl2);       
      }
        }
    else    //inactive leaf -- end if
    {
        if (  cl2->is_border() )
            this->b_leaves.push_vertex(cl2);       
    }
  }

  /* add edges */

  // Border graph
  //this->b_leaves.dfs(nlist);
  GraphDfsIter<Cell2d*> dfs(b_leaves);
  Cell2d* cl2= NULL;

  //foreach(Cell2d* cl2, nlist)
  for (dfs.first(); !dfs.isDone(); dfs.next() )
  {
    cl2= dfs.currentItem()->get_data();
    //std::cout<< *cl2 <<"\n";

    if (cl2->s_neighbors.size()==0)
      foreach(Cell2d* nb, cl2->e_neighbors )
        if ( this->b_leaves.member(nb) && nb->s_neighbors.size()==0)
          this->b_leaves.push_edge( cl2,nb ) ;
    if (cl2->e_neighbors.size()==0)
      foreach(Cell2d* nb, cl2->n_neighbors )
        if (this->b_leaves.member(nb) && nb->e_neighbors.size()==0)
          this->b_leaves.push_edge( cl2,nb ) ;
    if (cl2->n_neighbors.size()==0)
      foreach(Cell2d* nb, cl2->w_neighbors )
        if (this->b_leaves.member(nb)&& nb->n_neighbors.size()==0)
          this->b_leaves.push_edge( cl2,nb ) ;
    if (cl2->w_neighbors.size()==0)
             foreach(Cell2d* nb, cl2->s_neighbors )
               if (this->b_leaves.member(nb)&& nb->w_neighbors.size()==0)
                 this->b_leaves.push_edge( cl2,nb ) ;
  }
  //nlist.clear();


     // Treating singular bcells
     foreach( VCell* q, m_singular_cells )
     { 

       std::cout<<"NEED boundary \n";
       //q->compute_boundary();

     }

     //Computing Voronoi Cells:
     Seq<Cell2d*> nlist;
     this->b_leaves.dfs(nlist);// remove empty border bcells

     Cell2d* pr;
     pr= nlist[nlist.size()-1];

     foreach(Cell2d* cl2, nlist)
       //for (dfs.first(); !dfs.isDone(); dfs.next() )
      {
        // cl2= dfs.currentItem()->get_data();
        // std::cout<< *cl2 <<"\n";

        if ( cl2->is_corner() )
          pr=cl2;
        else {
          if ( (cl2->s_neighbors.size()==0 && 
                    cl2->intersections(0).size()==0 ) ||
                   (cl2->e_neighbors.size()==0 && 
                    cl2->intersections(1).size()==0 ) ||
                   (cl2->n_neighbors.size()==0 && 
                    cl2->intersections(2).size()==0 ) ||
                   (cl2->w_neighbors.size()==0 && 
                    cl2->intersections(3).size()==0 )  )
              {
                  this->b_leaves.push_edge(pr, this->b_leaves.next(cl2) ) ;
                  this->b_leaves.delete_vertex( cl2 ) ;
              }
              else
                  pr=cl2;
          }
      }

     std::cout<<"Border size= "<< nlist.size() <<std::endl;

      // Leaf graph
      //nlist.clear();
      //this->m_leaves.dfs(nlist);

      dfs= GraphDfsIter<Cell2d*>(m_leaves);

      for (dfs.first(); !dfs.isDone(); dfs.next() )
        //foreach(Cell2d* cl2, nlist)
      {
        cl2= dfs.currentItem()->get_data();
          //Cell * cl= dynamic_cast<Cell*>(cl);
          foreach(Cell2d* nb, cl2->neighbors() )
              this->m_leaves.push_edge( cl2,nb ) ;
      }
      
      std::cout<<"Leaf graph Ok." <<std::endl;



 //Remove inactive bcells OK ..
 // AND misleading edges 
 //Recover connected components .. ( for all regular bcells and sign +,-)
 // FOR ALL CC's:
 //1. check if CC is actually SCC .. 
 //2. walk on boundry and output face

   //get boundary

//   nlist.clear();
//   this->m_leaves.dfs(nlist);

   Point *p= NULL;
   Face * f= NULL;
   int sgn(1);
   unsigned site;//the site whose bcell is computed

   assert( nlist.size()>1);

   Cell2d// *s= NULL,
          *t= NULL,
          *b= NULL;


   // Start exploration
   //if (0) //Don't compute 2Dbcells
 for (dfs.first(); !dfs.isDone(); dfs.next() )
   {
     cl2= dfs.currentItem()->get_data();

     if (!( cl2->is_active()        &&
          (!this->b_leaves.member(cl2))   &&
          ((VCell*)cl2)->count() ==1 
           ))
     {  continue;  }
     

   //Every bisector box is visited twice
   switch ( cl2->m_gnode->aux() )
   {
   case 0:
     sgn=1 ;
   case -1 :
     sgn=1 ;
     break ;
   case  2 :
     sgn=-1;
   break   ;
   default :
     continue;
   }
   site=((VCell*)cl2)->site(sgn);//computing Voronoi Cell of "site"


   // Recovering site 
   f= new Face();

   // Get starting point (CCW)
   std::cout<<"Getting voronoi bcell of "<<site <<"("<<(sgn==1 ? "+": "-")
            <<"), starting "<< *cl2<<  std::endl;
   p= cl2->starting_point(sgn);
//   std::cout<<"starting at "<< *p <<std::endl;
   // Walking on CCW face
   p= cl2->pair(p,sgn);
//   std::cout<<"p= "<< *p <<std::endl;
   f->insert(p);
   b=cl2;
   
//int c(0); 
   do  {
//if (++c>137) {while(!stack.empty()) stack.pop(); exit(0);break;}
     //std::cout<<"Next"<< *b <<" sites= "<< ((VCell*)b)->sites()<<std::endl;
     //std::cout<<"."<<std::endl;
     //if ( !b->is_corner())


      if( this->m_leaves.member(b) ) {
      b->m_gnode->aux( b->m_gnode->aux() + (sgn==1 ? 2:-1) );
      //std::cout<<"aux="<<b->m_gnode->aux()<<std::endl; 
      }
      //else
      //  std::cout<<"Passed: "<< (*b) <<"\n";


      t= b->neighbor(p);

      //if ( t!=NULL)
      //  std::cout<<*t<<" is neighbor of "<<*b<<" for point "<<*p <<std::endl;

        // std::cout<<"Ints,b "<<b<<": "<<b->nb_intersect() <<std::endl;
        // foreach( Point* q, b->intersections() )
        //   std::cout<<"b: "<<": "<<*q  <<std::endl;
        // foreach( Cell2d* c, b->neighbors() )
        // {
        //   VCell* cc= dynamic_cast<VCell*>( c );
        //   std::cout<<"   Ints,t "<<"("<<cc->count() <<",sites="<<cc->sites()<<") "<<": "<<c->nb_intersect()<<", "<< *c <<", adr "<< cc <<std::endl;
        // foreach( Point* q, cc->intersections() )
        //   std::cout<<"   t: "<<": "<<*q  <<std::endl;

        //   foreach(Cell * m, cc->m_objects)        
        //     std::cout<<"   t-obj:"<< m <<std::endl;


        // }

      // if (b->m_singular.size() >0) 
      //  std::cout<<"Voronoi Vertex! "<< *b <<", " <<std::endl;

      //if (0)
      if ( t==NULL)
      { // border bcell reached

//        std::cout<<"Border Cell "<<*b <<", #="<< ((VCell*)b)->nb_intersect() <<std::endl;
        //foreach( Point* q, b->intersections() )
        //    std::cout<< *q<<"   adr  "<<q<<std::endl; 

     //check meeting corner
          if          (b->s_neighbors.size()==0 &&
                       b->e_neighbors.size()==0 && b->intersections(1).size()==0)
           f->insert(new Point(b->xmax(),b->ymin(),0.0) );
          if    (b->e_neighbors.size()==0 &&
                       b->n_neighbors.size()==0 && b->intersections(2).size()==0)
           f->insert(new Point(b->xmax(),b->ymax(),0.0) );
          if   (b->n_neighbors.size()==0 &&
                       b->w_neighbors.size()==0 && b->intersections(3).size()==0)
           f->insert(new Point(b->xmin(),b->ymax(),0.0) );
          if   (b->w_neighbors.size()==0 &&
                       b->s_neighbors.size()==0 && b->intersections(0).size()==0)
             f->insert(new Point(b->xmin(),b->ymin(),0.0) );

//          std::cout<<"Moving on border. ("<< this->m_faces.size()<<")"<< std::endl;
          b=this->b_leaves.next(b);
//          std::cout<<"Moved to "<<*b <<" sites="<<((VCell*)b)->sites()<<std::endl;

          //check leaving corner
          if          (b->s_neighbors.size()==0 &&
                       b->e_neighbors.size()==0 && b->intersections(0).size()==0 )
          { f->insert(new Point(b->xmax(),b->ymin(),0.0) );
          }
          else if   (b->e_neighbors.size()==0 &&
                       b->n_neighbors.size()==0 && b->intersections(1).size()==0)
          { f->insert(new Point(b->xmax(),b->ymax(),0.0) );
          }

          else if   (b->n_neighbors.size()==0 &&
                       b->w_neighbors.size()==0 && b->intersections(2).size()==0)
          { f->insert(new Point(b->xmin(),b->ymax(),0.0) );
          }
          else if   (b->w_neighbors.size()==0 &&
                       b->s_neighbors.size()==0 && b->intersections(3).size()==0)
          {   f->insert(new Point(b->xmin(),b->ymin(),0.0) );
          }//end check corner

          if (  this->m_leaves.member(b) ) 
          {   //entering point. Changing bisector
            sgn=((VCell*)b)->signof(site);
            p= b->starting_point(sgn);
            f->insert(p);
            //std::cout<<"Entering from "<<*b<<" ,sgn="<<sgn<<std::endl;
            //std::cout<<"insert "<<*p<<std::endl;
            p= b->pair(p,sgn);
            f->insert(p);
            //std::cout<<"insert "<<*p<<std::endl;
            //std::cout<<"To : "<< *(b->neighbor(p)) <<std::endl;
          }
      }
      else
      { // next normal bcell
        //std::cout<<"next normal bcell "<<std::endl;
          b=t;
          if ( b->m_singular.size() > 0 )
              f->insert(b->m_singular[0]);
          p= b->pair(p,sgn);
          f->insert(p);
          //std::cout<<"p="<<*p<<std::endl;
      }

   } while (b!=cl2); 
   //(b->m_gnode->aux()==1 || b->m_gnode->aux()==-1 || b->m_gnode->aux()==2) );

   std::cout<<"Face ADDED ("<< this->m_faces.size()+1<<")" << std::endl;
   f->set_index(site);
   this->m_faces<< f;

//  if (this->m_faces.size()== 14)  break;

   }// End exploration

  std::cout<<" # faces= "<< this->m_faces.size() << std::endl;


}

TMPL Viewer&
operator<<(Viewer& out, const SELF& tp) {

  foreach(Shape* vs, tp.m_objects)
  {
    out<<" <point color=\"255 127 0\"> "<<
      ((VoronoiSite2d *)vs)->x() << " "<<
      ((VoronoiSite2d *)vs)->y() << " "<<
      ((VoronoiSite2d *)vs)->z() << "</point>\n";
  }

  return out;
}

//====================================================================
} // namespace shape
} // namespace mmx
//====================================================================
# undef TMPL
# undef AlgebraicCurve
# undef VoronoiSite
# undef VCell
# undef Cell
# undef Cell2dVoronoiBisector
# undef Cell2dVoronoiImpl2d
# undef Cell2dFactory
# undef Curve
# undef SELF 
# undef VoronoiSite2d
# undef VoronoiCurve
# undef Viewer
# undef STACK
# undef Shape
# endif

---