/Users/mourrain/Devel/mmx/shape/include/shape/voronoi2d.hpp

Go to the documentation of this file.

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

# include <shape/topology2d.hpp>
# include <shape/voronoi_site2d.hpp>
# include <shape/cell2d_voronoi_site2d.hpp>
# include <shape/cell2d_voronoi_diagram.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 Cell2dVoronoiSite2d cell2d_voronoi_site2d<C,V>
# define SELF voronoi2d<C,V>
# define Viewer viewer<V>
//# define Cell  cell2d_voronoi_diagram<C,V>
# define VCell cell2d_voronoi_diagram<C,V>
# undef Cell
//====================================================================
namespace mmx {
namespace shape {

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

public:

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

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

  voronoi2d(Curve * curve) ;
  ~voronoi2d(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 * cell) 
    { 
      return cell->is_regular();
      //return topology2d<C,V>::is_regular(cell); 
    };

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

  Seq<Point*> l;
    l= ((VCell*)cell)->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, "<< *cell<<std::endl;

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

    foreach( Cell* c, ((VCell*)cell)->m_objects )
      if ( ((Cell2d*)c)->nb_intersect()==1 )
      {
        sz= ((Cell2dVoronoiSite2d*)c)->m_polynomial.rep().szs();
        st= ((Cell2dVoronoiSite2d*)c)->m_polynomial.rep().str();
        if ( ((Cell2dVoronoiSite2d*)c)->m_polynomial[0]<EPSILON )
        {
          q= new Point(cell->xmin(),cell->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);

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

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

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

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


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

        } else if ( ((Cell2dVoronoiSite2d*)c)->m_polynomial[sz[0]*sz[1]-1]<EPSILON )
        {
          q= new Point(cell->xmin(),cell->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);

          ((Cell2dVoronoiSite2d*)c)->w_intersections<< q;
          return true;
        } else if ( ((Cell2dVoronoiSite2d*)c)->m_polynomial[(sz[1]-1)*st[1]]<EPSILON )
        {
          q= new Point(cell->xmax(),cell->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);
          ((Cell2dVoronoiSite2d*)c)->e_intersections<< q;
          return true;
        }  
      }
  }

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

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

//    std::cout<< "nb_in= "<<  ((VCell*)cell)->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 * cell) 
//     { 
//       std::cout<< "Voronoi Singularity"<<std::endl;
// return topology2d<C,V>::singularity(cell); };

  bool subdivide(Cell * cell) 
    { return topology2d<C,V>::subdivide(cell,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
     //topology2d<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);
  //std::cout<<"voronoi diagram"<<std::endl;
  foreach( VSite * vcurve, this->m_objects) 
  {
    //VoronoiSite2d * vcurve = dynamic_cast<VoronoiSite2d *>(o);
    l->push_back( new Cell2dVoronoiSite2d(vcurve,this->m_bbx) );
  }

  std::cout<<"RUNing "<< l <<std::endl;

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

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

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

//    Node * node = this->m_nodes.front() ;
    Cell   * cl = cell_stack.top();
    cell_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;
            ((cell<C,V> *)cl)->subdivide(left, right) ;
          cell_stack.push(left);
          cell_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) ;
          cell_stack.push(left);
          cell_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 */
  Seq<Cell2d*> nlist;

  std::cout<<"border " << std::endl;
  // Border graph
  this->b_leaves.dfs(nlist);

  foreach(Cell2d* cl2, nlist)
  {
    if (cl2->s_neighbors.size()==0)
      foreach(Cell2d* nb, cl2->e_neighbors )
        if (nlist.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 (nlist.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 (nlist.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 (nlist.member(nb)&& nb->w_neighbors.size()==0)
                 this->b_leaves.push_edge( cl2,nb ) ;
  }
  nlist.clear();






     // Treating singular cells
     foreach( VCell* q, m_singular_cells )
     { q->compute_boundary();
       //m_singular_cells.erase(q);
     }

     //Computing Voronoi Cells:
     this->b_leaves.dfs(nlist);// remove empty border cells
     std::cout<<"Border size= "<< nlist.size() <<std::endl;

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

     foreach(Cell2d* cl2, nlist)
      {
          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 Ok." <<std::endl;

      // Leaf graph
      nlist.clear();
      this->m_leaves.dfs(nlist);
      foreach(Cell2d* cl2, nlist)
      {
          //Cell * cl= dynamic_cast<Cell*>(cl);
          foreach(Cell2d* nb, cl2->neighbors() )
              this->m_leaves.push_edge( cl2,nb ) ;
      }
      //remove interior cells
      if (0)//done at subdivision time
      foreach(Cell2d* cl, nlist) {
        if (!cl->is_active() )//|| !cl->is_touching() ) 
        {  
          this->m_leaves.delete_vertex(cl);
        }
      }
      
     
 std::cout<<"Leaf graph Ok." <<std::endl;



 //Remove inactive cells OK ..
 // AND misleading edges 
 //Recover connected components .. ( for all regular cells 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 aux;
   int sgn(1);
   unsigned site;//the site whose cell is computed

   assert( nlist.size()>1);

   Cell2d *s= NULL,
          *t= NULL,
          *b= NULL;
   STACK stack;

   // Start exploration
   foreach(Cell2d* cl, nlist)
   {
     if ( cl->is_active()        &&
          (!this->b_leaves.member(cl))   &&
          ((VCell*)cl)->count() ==1 
        )
     { cl->m_gnode->aux(0);
       stack.push(cl);  }
   }

   std::cout<<"#stack: "<<nlist.size()<<std::endl;

//   if (0) //Don't compute 2Dcells
while ( !stack.empty() )
   {
   s= stack.top();
   aux=s->m_gnode->aux();

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


   // Recovering site 
   f= new Face();

   // Get starting point (CCW)
   std::cout<<"Getting voronoi cell of "<<site <<"("<<(sgn==1 ? "+": "-")
            <<"), starting "<< *s<<  std::endl;
   p= s->starting_point(sgn);
//   std::cout<<"starting at "<< *p <<std::endl;
   // Walking on CCW face
   p= s->pair(p,sgn);
//   std::cout<<"p= "<< *p <<std::endl;
   f->insert(p);
   b=s;
   
//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) ) {
      aux=b->m_gnode->aux();
      b->m_gnode->aux(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 cell 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 cell
        //std::cout<<"next normal cell "<<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!=s); 
   //(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()== 6)  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 Cell2dFactory
# undef Curve
# undef SELF 
# undef VoronoiSite2d
# undef VoronoiCurve
# undef Viewer
# undef STACK
# undef Shape
# endif