include/shape/bcell2d_voronoi_diagram.hpp

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a  g i x
 *****************************************************************************
 * bcell2d_voronoi_diagram.hpp
 * 2010-04-10
 * Angelos Mantzaflaris
 *****************************************************************************
 *               Copyright (C) 2010 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
# ifndef shape_bcell2d_voronoi_diagram
# define shape_bcell2d_voronoi_diagram

# include <shape/topology.hpp>
# include <shape/bcell_list.hpp>
# include <shape/bcell2d_algebraic_curve.hpp>

# include <shape/bcell2d_voronoi_site2d.hpp>
# include <shape/solver_implicit.hpp>

# define EPSILON .000001

# define TMPL template<class C, class V>
# define Shape geometric<V>
# define Cell2dAlgebraicCurve bcell2d_algebraic_curve<C,V>
# define Intersection2dFactory intersection2d_factory<C,V>
# define Solver     solver_implicit<C,V>

# define Cell       bcell<C,V>
# define VSite      bcell2d_voronoi_site2d<C,V>
# define Cell2d     bcell2d<C,V>
# define CellList   bcell2d_list<C,V>
# define SELF bcell2d_voronoi_diagram<C,V>
namespace mmx {
namespace shape {

TMPL class voronoi2d;

template<class C, class V=default_env>
class bcell2d_voronoi_diagram
//: VSite
//: public bounding_box<C,V>
: public bcell2d<C,V>
   {
public:
     typedef typename Cell::BoundingBox      BoundingBox;
     typedef typename bcell2d<C,V>::Point     Point;
     typedef typename mesher2d<C,V>::Edge  Edge;
     typedef topology<C,V>                   Topology;
     //typedef voronoi2d<C,V>               Voronoi2d;
     typedef polynomial< Interval<double>, with<Bernstein> > Polynomial;
     typedef Interval<double> coeff;

  bcell2d_voronoi_diagram(void) ;
  bcell2d_voronoi_diagram(double xmin, double xmax) ;
  bcell2d_voronoi_diagram(double xmin, double xmax, double ymin, double ymax) ;
  bcell2d_voronoi_diagram(double xmin, double xmax, double ymin, double ymax, bool itr) ;
  bcell2d_voronoi_diagram(const BoundingBox& bx);
  virtual ~bcell2d_voronoi_diagram(void) ;

  virtual bool is_regular (void) ;// = 0 ;
  virtual bool is_active  (void) ;// = 0 ;
  virtual bool is_intersected(void) ;// = 0 ;

  virtual bool insert_regular(Topology*);// = 0 ;
  virtual bool process_singular();// = 0 ;

  virtual bool insert_singular(Topology*) {return false;};// not used

  virtual void subdivide     (Cell*& left, Cell*& right, int v, double s);

  virtual void split_position(int& v, double& t); 

  virtual Point* pair(Point * p, int & sgn);
  virtual Point* starting_point( int sgn);

  virtual bool  is_touching (void) {return true; };

  virtual unsigned nb_intersect(void) const;

    //Seq<Cell2d *>  neighbors();
    Cell2d * neighbor(Point * p);
  Seq<Point *> intersections() const;

  virtual Seq<Point *> intersections(int i) const{
      Seq<Point*> l;
      Cell2dAlgebraicCurve* c;
      
      foreach(Cell*m, this->m_objects)
      {
          c = dynamic_cast<Cell2dAlgebraicCurve*>(m);
          l<< c->intersections(i);
      }
      return l;
  }

     void push_back(Cell * cv) 
       {
         m_sites<< unsigned(m_objects.size() );
         m_objects.push_back(cv); 
       }; 

     void push_back(Cell * cv, unsigned k) {
       m_sites << k ; //object cv corresponds to site labeled k

       m_objects.push_back(cv); 
     }; 

     void push_bisector(Cell2dAlgebraicCurve * cv, Seq<unsigned> k) {
       m_sites << k ; 
       m_objects.push_back(cv); 
     }; 


     int  count(void) { return m_objects.size() ; }
     inline Polynomial func(const int i) const{
       return ((VSite*)(this->m_objects[0]))->m_polynomial;
     };

     inline int signature(int & i, int & j) {
       int c(0);
       for ( int u=0; u<i; u++)
         for ( int v=u+1; v<m_sites.size(); v++)
           c++;
       return c + j-i-1; }

     unsigned site(int sgn, unsigned bs=0) { 
       unsigned n=this->count();
       unsigned i,k(bs+1);
       
       for (i=1;i<n;i++)
         if ( k>n-i )
           k-= n-i;
         else
           break;

       //return( sgn>0 ? m_sites[i-1+k]: m_sites[i-1])  ; 
       return ( m_sites[i-1 + (sgn>0?k:0)] )  ; 
     }

     int signof(unsigned st, unsigned bs=0) {

       if (this->m_objects.size()==1 )
         return ( st==m_sites[0]? -1 :1 ); 

       //std::cout<<"Sgnof " <<st<<"(bs="<<bs<<") on "<< *this<<"objs="<<this->m_objects.size()<<std::endl;

       unsigned n=this->count() - 1;
       unsigned i,j(0),k(0);

       for (i=0;i<=bs;i++)
         if ( k<n )
           k++ ;
         else
           j++ ;

       if (st==this->m_sites[j])   return (-1);
       if (st==this->m_sites[k])   return (1);

       std::cout<<"problem at \"signof\" "<<i <<", "<< *this <<"sites= "<<m_sites<<std::endl;
       return 0;
     };

     Seq<unsigned> sites() const { return m_sites;};
     void addsite(unsigned i){m_sites<<i;};

  Seq<Cell *>        m_objects  ;

protected:


  Seq<unsigned>      m_sites;
  bool               m_intersected;
  bool               m_bisector;
  bool               m_treated;

private:

 template<int i> static bool 
 coord(Point*u,Point*v) { return ( (*u)[i]< (*v)[i]); }
  
  static bool comp_up( VSite* a, VSite* b)
    {
      return ( a->upper() < b->upper() );
    }

  static bool intersect( VSite* a, VSite* b)
    {
      return has_sign_variation( a->equation() - b->equation() );
    }


  inline bool over( VSite* a, double bound)
    {
      return ( a->lower() > bound );
    }

     bool is_bisector()
       {
         if ( this->m_objects.size() > 2)
           return false;
//         if ( this->m_bisector==true)
//           return true;
         
         if ( !this->m_bisector  &&  this->m_objects.size() == 2 )
         {// first call: compute f_1-f_2 and intersections
           Polynomial p=  ((VSite*)(this->m_objects[0]))->m_polynomial ;
           p -= ((VSite*)(this->m_objects[1]))->m_polynomial;
           
           Cell2dAlgebraicCurve* cc= 
             new Cell2dAlgebraicCurve( p, (BoundingBox)(*this), false);

           //std::cout<<"bisector: treating "<< *this <<"\n";           


           //0
             if ( this->s_neighbors.size()==0)
             {
               Seq<Point *> ip;
               Solver::edge_point(ip,
               cc->m_polynomial,
               Solver::south_edge, 
               *cc);
               cc->s_intersections << ip;
             }
             else
           foreach( Cell2d *c, this->s_neighbors )
           if  ( ((SELF*)c)->m_bisector ) 
             {
               foreach(Point * p, ((SELF*)c)->intersections(2) ) 
                 if (this->xmin()<p->x() && this->xmax()>p->x())
                   cc->s_intersections << p;
             }
             else
             {
               Seq<Point *> ip;
               Solver::edge_point(ip,
               cc->m_polynomial,
               Solver::south_edge, 
               *cc  );
               foreach(Point * p, ip ) 
                 if (c->xmin()<p->x() && c->xmax()>p->x())
                   cc->s_intersections << p;
               //std::cout<<"ip: "<< ip.size() <<"\n";
             }

           //1
             if (this->e_neighbors.size()==0)
             {
               Seq<Point *> ip;
               Solver::edge_point(ip,
               cc->m_polynomial,
               Solver::east_edge, 
               *cc);
               cc->e_intersections << ip;
             }
             else
             foreach( Cell2d *c, this->e_neighbors )
             if ( ((SELF*)c)->m_bisector ) 
             {
               //Cell2d *c;
               //foreach(Point * p, ((VSite*)( ((SELF*)c)->m_objects[0]))->w_intersections ) 
               foreach(Point * p, ((SELF*)c)->intersections(3) ) 
                 if (this->ymin()<p->y() && this->ymax()>p->y())
                   cc->e_intersections << p;
             }
             else
             {
               Seq<Point *> ip;
               Solver::edge_point(ip,
               cc->m_polynomial,
               Solver::east_edge, 
               *cc  );
               foreach(Point * p, ip ) 
                 if (c->ymin()<p->y() && c->ymax()>p->y())
                   cc->e_intersections << p;
             }

           //2
             if (this->n_neighbors.size()==0 )
             {
               Seq<Point *> ip;
               Solver::edge_point(ip,
               cc->m_polynomial,
               Solver::north_edge, 
               *cc);
               cc->n_intersections << ip;
             }
             else
           foreach( Cell2d *c, this->n_neighbors )
             if ( ((SELF*)c)->m_bisector ) 
             {
               //Cell2d *c;
               //foreach(Point * p, ((VSite*)( ((SELF*)c)->m_objects[0]))->s_intersections ) 
               foreach(Point * p, ((SELF*)c)->intersections(0) ) 
                 if (this->xmin()<p->x() && this->xmax()>p->x())
                   cc->n_intersections << p;
             }
             else
             {
               Seq<Point *> ip;
               Solver::edge_point(ip,
               cc->m_polynomial,
               Solver::north_edge, 
               *cc );
               foreach(Point * p, ip ) 
                 if (c->xmin()<p->x() && c->xmax()>p->x())
                   cc->n_intersections << p;
             }

           //3
             if (this->w_neighbors.size()==0)
             {
               Seq<Point *> ip;
               Solver::edge_point(ip,
               cc->m_polynomial,
               Solver::west_edge, 
               *cc);
               cc->w_intersections << ip;
             }
             else
           foreach( Cell2d *c, this->w_neighbors )
             if ( ((SELF*)c)->m_bisector ) 
             {
               //Cell2d *c;
               //foreach(Point * p, ((VSite*)( ((SELF*)c)->m_objects[0]))->e_intersections ) 
               foreach(Point * p, ((SELF*)c)->intersections(1) ) 
                 if (this->ymin()<p->y() && this->ymax()>p->y())
                   cc->w_intersections << p;

             }
             else
             {
               Seq<Point *> ip;
               Solver::edge_point(ip,
               cc->m_polynomial,
               Solver::west_edge, 
               *cc  );
               foreach(Point * p, ip ) 
                 if (c->ymin()<p->y() && c->ymax()>p->y())
                   cc->w_intersections << p;

               //std::cout<<"west: "<< *ip[0] <<"\n";
             }




           //foreach (Cell* m, this->m_objects) delete m; //segfault
           this->m_objects.clear();
           
           this->m_objects<< (Cell*)(cc);
           //std::cout<<"bisector found "<< this<<std::endl;

           //std::cout<<"Intersections: "<< std::endl;
           //foreach(Point* p, this->intersections() )
           //  std::cout<< *p<<"\n";
           
         }

         this->m_bisector=true;
         return true;

       }
     
public:

bool compute_boundary()
{
  int i(0),j(0), r(m_sites.size()-1);
  foreach ( Cell* obj, this->m_objects)
  {
    Cell2dAlgebraicCurve* cc= dynamic_cast<Cell2dAlgebraicCurve*>(obj);
    if (j<r) // bisector signature (i,j)
      j++;
    else
      { i++; j=i+1;}

    //0
    if ( this->s_neighbors.size()==0)
    {
      Seq<Point *> ip;
        Solver::edge_point(ip,
                           cc->m_polynomial,
                           Solver::south_edge, 
                           *cc);
        cc->s_intersections << ip;
      }
    else
      foreach( Cell2d *c, this->s_neighbors )
        if  ( ((SELF*)c)->m_bisector ) 
        {
          Seq<unsigned> a= ((SELF*)c)->sites();
            if ( a.member( m_sites[i]) && a.member(m_sites[j]) )
            {
              foreach(Point * p, ((SELF*)c)->intersections(2) ) 
                if (this->xmin()<p->x() && this->xmax()>p->x())
                  cc->s_intersections << p;
            }
          }
        else if ( ((SELF*)c)->m_treated )
        {
              int u(0),v(0);
              foreach( Cell *nc, ((SELF*)c)->m_objects )
              {
                if (v<r) // bisector signature (u,v)
                  v++;
                else
                { u++; v=u+1;}
                  if ( ((SELF*)c)->m_sites[u]==m_sites[i] && 
                       ((SELF*)c)->m_sites[v]==m_sites[j] )
              foreach(Point * p, ((SELF*)nc)->intersections(2) ) 
                if (this->xmin()<p->x() && this->xmax()>p->x())
                  cc->s_intersections << p;
          }
       }
        else
        {
          Seq<Point *> ip;
            Solver::edge_point(ip,
                               cc->m_polynomial,
                               Solver::south_edge, 
                               *cc  );
            foreach(Point * p, ip ) 
              if (c->xmin()<p->x() && c->xmax()>p->x())
                cc->s_intersections << p;
            //std::cout<<"ip: "<< ip.size() <<"\n";
          }
        
        //1
        if (this->e_neighbors.size()==0)
        {
          Seq<Point *> ip;
            Solver::edge_point(ip,
                               cc->m_polynomial,
                               Solver::east_edge, 
                               *cc);
               cc->e_intersections << ip;
          }
        else
          foreach( Cell2d *c, this->e_neighbors )
            if ( ((SELF*)c)->m_bisector ) 
            {
              Seq<unsigned> a= ((SELF*)c)->sites();
              if ( a.member( m_sites[i]) && a.member(m_sites[j]) )
              {
                foreach(Point * p, ((SELF*)c)->intersections(3) ) 
                  if (this->ymin()<p->y() && this->ymax()>p->y())
                    cc->e_intersections << p;
              }
            }
            else if ( ((SELF*)c)->m_treated )
            {
              int u(0),v(0);
              foreach( Cell *nc, ((SELF*)c)->m_objects )
              {
                if (v<r) // bisector signature (u,v)
                  v++;
                else
                { u++; v=u+1;}
                  if ( ((SELF*)c)->m_sites[u]==m_sites[i] && 
                       ((SELF*)c)->m_sites[v]==m_sites[j] )
                    foreach(Point * p, ((SELF*)nc)->intersections(3) ) 
                      if (this->ymin()<p->y() && this->ymax()>p->y())
                        cc->e_intersections << p;                    
                    }
              }
            else
            {
              Seq<Point *> ip;
                Solver::edge_point(ip,
                                   cc->m_polynomial,
                                   Solver::east_edge, 
                                   *cc  );
                 foreach(Point * p, ip ) 
                   if (c->ymin()<p->y() && c->ymax()>p->y())
                     cc->e_intersections << p;
              }
            
            //2
             if (this->n_neighbors.size()==0 )
             {
               Seq<Point *> ip;
                 Solver::edge_point(ip,
                                    cc->m_polynomial,
                                    Solver::north_edge, 
                                    *cc);
                 cc->n_intersections << ip;
               }
             else
               foreach( Cell2d *c, this->n_neighbors )
                 if ( ((SELF*)c)->m_bisector ) 
                 {
                   Seq<unsigned> a= ((SELF*)c)->sites();
                     if ( a.member( m_sites[i]) && a.member(m_sites[j]) )
                     {
                   foreach(Point * p, ((SELF*)c)->intersections(0) ) 
                     if (this->xmin()<p->x() && this->xmax()>p->x())
                       cc->n_intersections << p;
                     }
                   }
            else if ( ((SELF*)c)->m_treated )
            {
              int u(0),v(0);
              foreach( Cell *nc, ((SELF*)c)->m_objects )
              {
                if (v<r) // bisector signature (u,v)
                  v++;
                else
                { u++; v=u+1;}
                  if ( ((SELF*)c)->m_sites[u]==m_sites[i] && 
                       ((SELF*)c)->m_sites[v]==m_sites[j] )
                   foreach(Point * p, ((SELF*)nc)->intersections(0) ) 
                      if (this->xmin()<p->x() && this->xmax()>p->x())
                        cc->n_intersections << p;                    
                }
              }
                 else
                 {
                   Seq<Point *> ip;
                     Solver::edge_point(ip,
                                        cc->m_polynomial,
                                        Solver::north_edge, 
                                        *cc );
                     foreach(Point * p, ip ) 
                       if (c->xmin()<p->x() && c->xmax()>p->x())
                         cc->n_intersections << p;
                   }
                 
                 //3
                 if (this->w_neighbors.size()==0)
                 {
                   Seq<Point *> ip;
                 Solver::edge_point(ip,
                                    cc->m_polynomial,
                                    Solver::west_edge, 
                                    *cc);
                 cc->w_intersections << ip;
               }
                 else
               foreach( Cell2d *c, this->w_neighbors )
                 if ( ((SELF*)c)->m_bisector ) 
                 {
                   Seq<unsigned> a= ((SELF*)c)->sites();
                     if ( a.member( m_sites[i]) && a.member(m_sites[j]) )
                     {
                   //Cell2d *c;
                   //foreach(Point * p, ((VSite*)( ((SELF*)c)->m_objects[0]))->e_intersections ) 
                   foreach(Point * p, ((SELF*)c)->intersections(1) ) 
                     if (this->ymin()<p->y() && this->ymax()>p->y())
                       cc->w_intersections << p;
                     }
                   }
            else if ( ((SELF*)c)->m_treated )
            {
              int u(0),v(0);
              foreach( Cell *nc, ((SELF*)c)->m_objects )
              {
                if (v<r) // bisector signature (u,v)
                  v++;
                else
                { u++; v=u+1;}
                  if ( ((SELF*)c)->m_sites[u]==m_sites[i] && 
                       ((SELF*)c)->m_sites[v]==m_sites[j] )
                   foreach(Point * p, ((SELF*)nc)->intersections(1) ) 
                     if (this->ymin()<p->y() && this->ymax()>p->y())
                       cc->w_intersections << p;                    
                    }
              }
                 else
                 {
                   Seq<Point *> ip;
                     Solver::edge_point(ip,
                                        cc->m_polynomial,
                                        Solver::west_edge, 
                                        *cc  );
                     foreach(Point * p, ip ) 
                       if (c->ymin()<p->y() && c->ymax()>p->y())
                         cc->w_intersections << p;
                     
                     //std::cout<<"west: "<< *ip[0] <<"\n";
                   }
                 
        }
  this->m_treated=true;
  return true;
}//compute_boundary
               
} ;


//--------------------------------------------------------------------

/*
TMPL
class intersection2d_factory {
public:
  typedef bounding_box<C,V>      BoundingBox;
  typedef typename bcell2d<C,V>::Point Point;
  static Intersection2dFactory * instance(void) {

        if(!m_instance)
            m_instance = new Intersection2dFactory ;

        return m_instance ;
    }

  void compute(Seq<Point*>& res, Shape * c1, Shape * c2, const BoundingBox& bx) {
    
    if(Cell2dAlgebraicCurve * ca1 = dynamic_cast<Cell2dAlgebraicCurve *>(c1)) {
      if(Cell2dAlgebraicCurve * ca2 = dynamic_cast<Cell2dAlgebraicCurve *>(c2)) {

        typedef polynomial< double, with<Bernstein> > MultivariateDenseBernstein;
        MultivariateDenseBernstein p1, p2;
        let::assign(p1.rep(), ca1->m_polynomial.rep() );
        let::assign(p2.rep(), ca2->m_polynomial.rep() );
        Solver::intersection(res, p1, p2, bx);
      }
    }
  }

private:
    intersection2d_factory(void) {} ;

private:
    static Intersection2dFactory * m_instance;
} ;

TMPL Intersection2dFactory*  Intersection2dFactory::m_instance = NULL;

*/

TMPL
SELF::bcell2d_voronoi_diagram(void): m_intersected(false), m_treated(false) {}

TMPL
SELF::bcell2d_voronoi_diagram(double xmin, double xmax, double ymin, double ymax) :  bcell2d<C,V>(xmin, xmax, ymin, ymax), m_intersected(false), m_bisector(false), m_treated(false) {}

TMPL
SELF::bcell2d_voronoi_diagram(double xmin, double xmax, double ymin, double ymax, bool itr) : bcell2d<C,V>(xmin, xmax, ymin, ymax), m_intersected(itr), m_bisector(false), m_treated(false) {}

TMPL
SELF::bcell2d_voronoi_diagram(const BoundingBox& bx): bcell2d<C,V>(bx), m_intersected(false), m_bisector(false), m_treated(false)  {};

TMPL
SELF::~bcell2d_voronoi_diagram(void) {
    foreach (Cell* m, this->m_objects) delete m;
}
  
TMPL bool 
SELF::is_regular()
{
  //std::cout<<"bisector? "<< this->is_bisector()<<std::endl;
  if ( this->is_bisector() ) 
  {
    return (this->m_objects[0])->is_regular();
  }
  else
  {
    return false;
  }
}

TMPL bool 
SELF::is_intersected() {
  
  if(this->m_objects.size() >1 && !m_intersected) {
//      std::cout<<"Intersecting inside box "<< this <<std::endl;
      for(unsigned i=0; i<this->m_objects.size();i++)
          for(unsigned j=i+1; j<this->m_objects.size(); j++)
            Intersection2dFactory::instance()->compute(this->m_singular, (Shape*)this->m_objects[i], (Shape*)this->m_objects[j], (BoundingBox)*this);
      m_intersected = true;      
  }
  
  if (this->m_singular.size() > 0) return true;
  return false;
}


TMPL bool 
SELF::insert_regular(Topology* s) {

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

  
  //std::cout<<"VD, regular: "<< *this<<", #="<<l.size() <<std::endl;
  // foreach( Point* e, l)
  //    std::cout<<*e <<", at "<< e <<std::endl;

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

  if ( l.size()==2 )
  {
    s->insert( l[0] );
    s->insert( l[1] );
    s->insert( new Edge(l[0],l[1]) );
    return true;
  }

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

    s->insert( l[0] );
    foreach( Cell* c, this->m_objects )
      if ( ((Cell2d*)c)->nb_intersect()==1 )
      {
        sz= ((VSite*)c)->m_polynomial.rep().szs();
        st= ((VSite*)c)->m_polynomial.rep().str();
        if ( ((VSite*)c)->m_polynomial[0]<EPSILON )
        {
          q= new Point(this->xmin(),this->ymin(),0); 
          std::cout<< "1.add ("<< *q <<")->("<< *l[0] << ") in "<< *this<<std::endl;
          s->insert(q);
          s->insert( new Edge(l[0],q) );
          ((VSite*)c)->n_intersections<< q;

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

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


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

        } else if ( ((VSite*)c)->m_polynomial[sz[0]*sz[1]-1]<EPSILON )
        {
          q= new Point(this->xmin(),this->ymax(),0);
          std::cout<< "3.add ("<< *q <<")->("<< *l[0] << ") in "<< *this<<std::endl;
          s->insert(q);
          s->insert( new Edge(l[0],q) );
          ((VSite*)c)->w_intersections<< q;
          return true;
        } else if ( ((VSite*)c)->m_polynomial[(sz[1]-1)*st[1]]<EPSILON )
        {
          q= new Point(this->xmax(),this->ymin(),0);
          std::cout<< "4.add ("<< *q <<")->("<< *l[0] << ") in "<< *this<<std::endl;
          s->insert(q);
          s->insert( new Edge(l[0],q) );
          ((VSite*)c)->e_intersections<< q;
          return true;
        }  
      }
  }

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

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

    std::cout<< "nb_in= "<<  this->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;
}
  
TMPL bool 
SELF::process_singular() {

//  std::cout<<"VD, Inserting singular"<<*this<<std::endl;  
//  ((voronoi2d<C,V>*)s)->m_singular_cells<< this; // does not work with forward declaration of voronoi2d
 
//  foreach( Point* q, this->intersections() )
//    std::cout<< *q<<"   adr  "<<q<<std::endl; 


     //CellList* l= new CellList( (BoundingBox)(*this) );
     Seq<Cell2dAlgebraicCurve*> l;
     Cell2dAlgebraicCurve* cc;

     //std::cout<<"compute arrangement of bisectors"<<std::endl;
     //More than one bisector in the bcell: Compute arrangement
     Polynomial p;
     for ( unsigned i=0; i< this->m_objects.size(); i++ )
       for ( unsigned j=i+1; j< this->m_objects.size(); j++ )
       {
         p  = ((VSite*)(this->m_objects[i]))->m_polynomial ;
         p -= ((VSite*)(this->m_objects[j]))->m_polynomial ;
         
         cc= new Cell2dAlgebraicCurve( p,  (BoundingBox)(*this), false  );
         //if ( cc->is_active() ) //commented:Put all curves, even inactive!
           l<< cc;
       }

     this->m_objects.clear();
     this->m_objects<< l;
     this->is_intersected();
//     foreach(Cell * m, this->m_objects) m->insert_singular(s);

     return true;

     this->m_bisector=true;// .. 

     //foreach(Cell * m, l ) m->insert_singular(s);

     //if (this->m_singular.size()>0 ) 
     //  std::cout<<"VORONOI VERTEX COMPUTED in"<< *this <<std::endl;
   
}
  
TMPL void
SELF::split_position(int& v, double& s) {
  double sx = (this->xmax()-this->xmin());
  double sy = (this->ymax()-this->ymin());
  if(sx<sy) {
    v=1;
    s=(this->ymax()+this->ymin())/2;
  } else {
    v=0;
    s=(this->xmax()+this->xmin())/2;
  }
}

TMPL void
SELF::subdivide(Cell *& left, Cell *& right, int  v, double c) {

  //std::cout<<"Subdividing "<< this << "sites="<<this->m_sites <<std::endl;

  typedef SELF Cell_t;

  if(v==1) {
    left =(Cell*)new Cell_t(this->xmin(), this->xmax(), this->ymin(), c, m_intersected) ;
    right=(Cell*)new Cell_t(this->xmin(), this->xmax(), c, this->ymax(), m_intersected) ;
    
    foreach(Point * p, this->m_singular) {
      if(p->y() <=  c) 
        ((Cell_t*) left)->m_singular << p ;
      else
        ((Cell_t*)right)->m_singular << p ;
    }

    /*  Update neighbors  */
    this->connect1( (Cell_t*)left, (Cell_t*)right);
    ((Cell_t*)left)->join1((Cell_t*)right);

  } else {//v==0
    left = (Cell*)new Cell_t(this->xmin(), c, this->ymin(), this->ymax(), m_intersected) ;
    right= (Cell*)new Cell_t(c, this->xmax(), this->ymin(), this->ymax(), m_intersected) ;

    foreach(Point * p, this->m_singular) {
      if(p->x() <= c ) 
        ((Cell_t*)left)->m_singular << p ;
      else
        ((Cell_t*)right)->m_singular << p ;
    }

    /*  Update neighbors  */
    this->connect0((Cell_t*)left, (Cell_t*)right);
    ((Cell_t*)left)->join0((Cell_t*)right);

  }

  /* disconnect parent */
  this->disconnect( );

  if (!this->is_bisector() ) 
  {

  Seq<VSite*> ll, rr;
  Cell * cv_left, * cv_right;

  foreach(Cell* cv, this->m_objects) {
    cv->subdivide( cv_left, cv_right);
      ll<< (VSite*)cv_left;
      rr<< (VSite*)cv_right;
  }

  /* Filtering sites that are "far away" */
  //1. find minimum upper bound 
  //2. remove all bcells whose lower bound is bigger than that
  //std::cout<<"Filtering Cells "<< cv_left<<" , "<<cv_right<<std::endl;

  /* Left bcell */
  double mm;
  unsigned cnt;
  mm=  ( (VSite*)ll.min(this->comp_up) )->upper();
  
  cnt=0;
  foreach(VSite* vs, ll)
  {
    //std::cout<<"check "<< m_sites[cnt]<<std::endl;
    if ( !this->over(vs,mm) )
    {
      ((Cell_t*)left)->push_back( vs, m_sites[cnt] );
      //std::cout<<"added "<< m_sites[cnt]<<std::endl;
    }
    cnt++;
    //else
    //{
    //delete vs;
    //ll.erase( ll.search(vs) ) ;
    //}
  }

  /* Right bcell */
  mm=  ( (VSite*)rr.min(this->comp_up) )->upper();
  
  cnt=0;
  foreach(VSite* vs, rr)
  {
    if ( !this->over(vs,mm) )
    {
      ((Cell_t*)right)->push_back( vs, m_sites[cnt] );
      //((Cell_t*)right)->m_objects<< (Cell*)vs;
    }
    //else
    //{
    //delete vs;
    //ll.erase( ll.search(vs) ) ;
    //}
    cnt++;
  }
  }
  else
  {//bisector box

    Cell * cv_left, * cv_right;
    //std::cout<<"bisector subdiv " <<std::endl;
    this->m_objects[0]->subdivide( cv_left, cv_right);
    ((Cell_t*)left)->push_bisector( (Cell2dAlgebraicCurve*)cv_left,  m_sites );
    ((Cell_t*)right)->push_bisector((Cell2dAlgebraicCurve*)cv_right, m_sites );

    ((Cell_t*)left)->m_bisector=true;
    ((Cell_t*)right)->m_bisector=true;
    
  }

//    std::cout<<"ok " <<std::endl;  
}

//----------

TMPL  bool 
SELF::is_active()  {
  if ( this->is_bisector() ) 
  {
    return (this->m_objects[0])->is_active();
  }
  else
  {
    return true;
  }

}

TMPL unsigned 
SELF::nb_intersect(void) const {
    unsigned sum(0);
    Cell2d* c;
    foreach (Cell* m, m_objects)
    { 
        c = dynamic_cast<Cell2d*>(m);
        sum+= c->nb_intersect();
    }
    return sum;
  }


TMPL typename SELF::Point* 
SELF::pair(Point * p, int & sgn)
{

//  std::cout<<(sgn>0?"+": "-")<<" pair of "<< *this<<"("<< p->x()<<","<<p->y()<<")"<<std::endl;

  if ( this->is_intersected() )  
  {
    //std::cout<<"Reached Vertex bcell "<< *this << " ("<<sgn<<")"<<std::endl;
    //std::cout<<"sites are "<< this->m_sites <<std::endl;
    //foreach( Point* q, this->intersections() )
    //std::cout<< *q<<std::endl; 
  
    //Find bisector containing p
    
    unsigned st, cnt(0);
    foreach (Cell* mm, this->m_objects)
    { 
      Cell2d* m = dynamic_cast<Cell2d*>(mm);

      if ( m->intersections().member(p) )
//      foreach( Point* q, m->intersections() )
//        if ( abs(p->x()-q->x())<EPSILON && abs(p->y()-q->y())<EPSILON )
        { 
          // p found in bisector m.
          //std::cout<< "Current bisector: "<< cnt <<std::endl;
          // deduce site from sgn
          st=this->site(sgn,cnt);
          //std::cout<< "Site is: "<< st <<std::endl;
          // jump to a different bisector containing st
          unsigned cnt2(0);
          foreach (Cell* cc, this->m_objects)
          { 
            Cell2d* c = dynamic_cast<Cell2d*>(cc);
            if ( c!=m                   &&
               (st==this->site(1,cnt2)  ||
                st==this->site(-1,cnt2) ))
            {
              sgn= this->signof(st,cnt2);
//              std::cout<< "New bisector: "<< cnt2 <<std::endl;
//              std::cout<< "New sign: "<< sgn <<std::endl;
//              std::cout<< "start="<< *(c->starting_point(sgn))<<std::endl;
//              std::cout<< "end  ="<< *(c->pair(c->starting_point(sgn), sgn )) <<std::endl;

              return c->pair(c->starting_point(sgn), sgn );
            } 
            cnt2++;
          }
        }
      cnt++;
    }

  } else {
    
    Cell2d* c;
    foreach (Cell* m, m_objects)
    { 
      c = dynamic_cast<Cell2d*>(m);
      if ( c->intersections().member(p) )
//    foreach( Point* q, c->intersections() )
//      if ( abs(p->x()-q->x())<EPSILON && abs(p->y()-q->y())<EPSILON )
      { 
        //return c->pair(q,sgn);
        return c->pair(p,sgn);
      }
    }
  }
  std::cout<<"... Cell list pair trouble"<<std::endl;
  return NULL;
} 

TMPL typename SELF::Point* 
SELF::starting_point( int sgn)
{
//  std::cout<<"starting point(voronoi diagram)"<<std::endl;

   Cell2d* c;
   foreach (Cell* m, m_objects)
    { 
        if ( m->is_active() )
        {
            c = dynamic_cast<Cell2d*>(m);           
            return ((VSite*)c)->starting_point(sgn);
        }
    }
    return NULL; 
}

TMPL Cell2d*
SELF::neighbor( Point* p)
{
  foreach( Cell2d *c, this->s_neighbors     )
    if ( c->intersections(2).member(p) )
//    foreach( Point* q, c->intersections(2) )
//    if ( abs(p->x()-q->x())<EPSILON && abs(p->y()-q->y())<EPSILON )
      return c;

  foreach( Cell2d *c, this->e_neighbors     )
    if ( c->intersections(3).member(p) )
//    foreach( Point* q, c->intersections(3) )
//    if ( abs(p->x()-q->x())<EPSILON && abs(p->y()-q->y())<EPSILON )
      return c;

  foreach( Cell2d *c, this->n_neighbors     )
    if ( c->intersections(0).member(p) )
//    foreach( Point* q, c->intersections(0) )
//    if ( abs(p->x()-q->x())<EPSILON && abs(p->y()-q->y())<EPSILON )
      return c;

  foreach( Cell2d *c, this->w_neighbors     )
    if ( c->intersections(1).member(p) )
//    foreach( Point* q, c->intersections(1) )
//    if ( abs(p->x()-q->x())<EPSILON && abs(p->y()-q->y())<EPSILON )
      return c;
  
/*

  std::cout<<"Point "<<"("<<*p<<") not found in neighbs of "<< *this<<"( #neibs="<<this->neighbors().size() <<")"<<std::endl;

      foreach( Cell2d* c, this->neighbors() )
      {
        std::cout<< *c <<" ints:"<< ((SELF*)c)->intersections().size()  <<std::endl;
        foreach(Point* q, c->intersections() )
          std::cout<<"Point "<<q<<" ("<<*q<<")"<<std::endl;
      }
*/

    return NULL;
}

TMPL Seq<typename mmx::shape::bcell2d<C,V>::Point*>
SELF::intersections() const{
    Seq<Point *> s,e,n,w,r;
    //std::cout<<"intersections, "<<*this<<std::endl;
    Cell2d* cl;
    foreach (Cell* m, m_objects)
    {
        cl = dynamic_cast<Cell2d*>(m);           
        s<< cl->s_intersections;
        e<< cl->e_intersections;
        n<< cl->n_intersections;
        w<< cl->w_intersections;
    }
    s.sort(this->coord<0>);
    e.sort(this->coord<1>);
    n.sort(this->coord<0>);
    w.sort(this->coord<1>);

    r<<s;
    r<<e;
    r<<n.reversed();
    r<<w.reversed();


    return ( r ); 
}


} ; // namespace shape
} ; // namespace mmx

# undef TMPL
# undef Solver
# undef Cell
# undef Cell2d
# undef CellList
# undef SELF
# undef AlgebraicCurve
# undef Intersection2dFactory
# undef Cell2dAlgebraicCurve
# undef Shape
# undef VSite

# endif 

---