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

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a  g i x
 *****************************************************************************
 * Cell
 * 2008-03-20
 * Julien Wintz & Bernard Mourrain
 *****************************************************************************
 *               Copyright (C) 2008 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/

# ifndef shape_cell2d_List_hpp
# define shape_cell2d_List_hpp

# include <shape/topology.hpp>
# include <shape/cell_list.hpp>
# include <shape/cell2d_algebraic_curve.hpp>
# include <shape/solver_implicit.hpp>

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

# define Cell       cell<C,V>
# define Cell2d     cell2d<C,V>
# define CellList   cell_list<C,V>
# define Cell2dList cell2d_list<C,V>
namespace mmx {
namespace shape {


TMPL
class cell2d_list: public cell2d<C,V>   {
public:
  typedef typename Cell::BoundingBox  BoundingBox;
  typedef typename cell2d<C,V>::Point   Point;
  typedef topology<C,V>                 Topology;

  cell2d_list(void) ;
  cell2d_list(double xmin, double xmax) ;
  cell2d_list(double xmin, double xmax, double ymin, double ymax) ;
  cell2d_list(double xmin, double xmax, double ymin, double ymax, bool itr) ;
  cell2d_list(const BoundingBox& bx);
  virtual ~cell2d_list(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 insert_singular(Topology*);// = 0 ;

  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_objects.push_back(cv); }; 
  int  count(void) { return m_objects.size() ; }

protected:
  Seq<Cell *>        m_objects  ;

  //Seq<Point*>        m_singulars;
  bool               m_intersected;

private:

 template<int i> static bool 
 coord(Point*u,Point*v) { return ( (*u)[i]< (*v)[i]); }
  
} ;

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

TMPL
class intersection2d_factory {
public:
  typedef bounding_box<C,V>      BoundingBox;
  typedef typename cell2d<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
Cell2dList::cell2d_list(void): m_intersected(false) {}

TMPL
Cell2dList::cell2d_list(double xmin, double xmax, double ymin, double ymax) :  cell2d<C,V>(xmin, xmax, ymin, ymax), m_intersected(false) {}

TMPL
Cell2dList::cell2d_list(double xmin, double xmax, double ymin, double ymax, bool itr) : cell2d<C,V>(xmin, xmax, ymin, ymax), m_intersected(itr) {}

TMPL
Cell2dList::cell2d_list(const BoundingBox& bx): cell2d<C,V>(bx), m_intersected(false)  {};

TMPL
Cell2dList::~cell2d_list(void) {
    foreach (Cell* m, this->m_objects) delete m;
}
  
TMPL bool 
Cell2dList::is_regular(){
  foreach (Cell* m, this->m_objects) 
    if(!m-> is_regular()) return false;

  if ( is_intersected() )
      return ( (this->m_singular.size()<2) &&
               (this->nb_intersect()<5) );
  else
      return (this->nb_intersect()<3) ;
}

TMPL bool 
Cell2dList::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 
Cell2dList::insert_regular(Topology* s) {
  foreach(Cell*  m, this->m_objects)   m->insert_regular(s);
  if (this->m_singular.size() > 0) 
    foreach(Point* p, this->m_singular) s->insert_singular(p);
  return true;
}
  
TMPL bool 
Cell2dList::insert_singular(Topology* s) {
  //s->singular(this);
  foreach(Cell * m, this->m_objects) m->insert_singular(s);
  return true;
}
TMPL void
Cell2dList::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
Cell2dList::subdivide(Cell *& left, Cell *& right, int  v, double c) {

  typedef Cell2dList 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( );

  Cell * cv_left, * cv_right;
  foreach(Cell* cv, this->m_objects) {
      cv->subdivide( cv_left, cv_right);
    if(cv_left->is_active())
        ((Cell_t*)left)->m_objects<<cv_left;
    if(cv_right->is_active())
        ((Cell_t*)right)->m_objects<<cv_right;
  }
}

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

TMPL  bool 
Cell2dList::is_active()  {
    foreach (Cell* m, m_objects) 
    if(m->is_active()) return true;
  return false;
}

TMPL unsigned 
Cell2dList::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 Seq<Cell2d *> 
// Cell2dList:: neighbors()
//    {
//     Seq<Cell2d *> r;
//     Cell2d* c;
//     foreach (Cell* m, m_objects)
//     { 
//         c = dynamic_cast<Cell2d*>(m);
//         //std::cout<< c <<" has "<< c->neighbors()<<std::endl;
//         r<< c->neighbors();
//     }
//     //std::cout<< " Total="<< r.size()<<std::endl;
//   return ( r ); 
//    }


TMPL typename Cell2dList::Point* 
Cell2dList::pair(Point * p, int &  sgn) // (!) sgn is updated inside 
{
    if ( this->is_intersected() ) {
        //
      //std::cout<<"Reached intersection cell "<< *this<<std::endl;
      // std::cout<<"p= "<< *p <<std::endl;
      //foreach (Point*q, this->intersections() ) 
      //std::cout<<"q= "<< *q <<std::endl;

        Seq<Point*> l0,l1,l2,l3;
        Cell2dAlgebraicCurve* c, *v(NULL);
        double ev(0); int i, u(0);

        foreach (Cell* m, m_objects)
        {
            c = dynamic_cast<Cell2dAlgebraicCurve*>(m);
            int * sz = c->m_polynomial.rep().szs();
            int * st = c->m_polynomial.rep().str();
      
            l0<< c->s_intersections;
            if (c->s_intersections.member(p))
            {u=0; ev= (c->m_polynomial[0] >0 ? 1:-1); v=c;}
            l1<< c->e_intersections;
            if (c->e_intersections.member(p))
            {u=1; ev= (c->m_polynomial[(sz[0]-1)*st[0]] >0 ? 1:-1); v=c;}
            l2<< c->n_intersections;
            if (c->n_intersections.member(p))
            {u=2;ev=(c->m_polynomial[sz[0]*sz[1]-1]>0 ? 1:-1); v=c;}
            l3<< c->w_intersections;
            if (c->w_intersections.member(p))
             {u=3;ev=(c->m_polynomial[(sz[1]-1)*st[1]] >0 ? 1:-1); v=c;}
        }
//        std::cout<<"IN: u="<<u<<", ev="<<ev<<std::endl;
        l0.sort(this->coord<0>);
        l1.sort(this->coord<1>);
        l2.sort(this->coord<0>);
        l3.sort(this->coord<1>);
        l0<<l1; l0<<l2.reversed(); l0<<l3.reversed();
        
        foreach (Point*q, v->intersections(u) ) 
         {
             if (q==p) break;
             ev*=-1;
         }
        

        i=l0.search(p);
        int a=l0.size();

        Point *q;
        if (ev*sgn>0)
            q=l0[(i!=0 ?i-1:a-1)];
        else
            q= l0[(i!=a-1 ?i+1:0)];

        //std::cout<<"res="<< *q <<std::endl;

        int ev2(444);
        foreach (Cell* m, m_objects)
        {
            c = dynamic_cast<Cell2dAlgebraicCurve*>(m);
            int * sz = c->m_polynomial.rep().szs();
            int * st = c->m_polynomial.rep().str();

            if (c->s_intersections.member(q))
            {u=0; ev2= (c->m_polynomial[0] >0 ? 1:-1); v=c;}
            if (c->e_intersections.member(q))
            {u=1; ev2= (c->m_polynomial[(sz[0]-1)*st[0]] >0 ? 1:-1); v=c;}
            if (c->n_intersections.member(q))
            {u=2;ev2=(c->m_polynomial[sz[0]*sz[1]-1]>0 ? 1:-1); v=c;}
            if (c->w_intersections.member(q))
             {u=3;ev2=(c->m_polynomial[(sz[1]-1)*st[1]] >0 ? 1:-1); v=c;}
        }
        foreach (Point*w, v->intersections(u) ) 
         {
             if (q==w) break;
             ev2*=-1;
         }

//        std::cout<<"Entered from branch "<<sgn<<std::endl;

        if (ev*sgn>0)
            sgn=-ev2;
        else
            sgn=ev2;

//        std::cout<<"exiting to "<<sgn<<std::endl;
//        std::cout<<"OUT: u="<<u<<", ev2="<<ev2<<std::endl;

        return q;

    } else {

        Cell2d* c;
        foreach (Cell* m, m_objects)
        { 
            c = dynamic_cast<Cell2d*>(m);
            if ( c->intersections().member(p) )
                return c->pair(p,sgn);
        }
    }
    std::cout<<"... Cell list pair trouble"<<std::endl;
    return NULL;
} 

TMPL typename Cell2dList::Point* 
Cell2dList::starting_point( int sgn)
{
   Cell2d* c;
   foreach (Cell* m, m_objects)
    { 
        if ( m->is_active() )
        {
            c = dynamic_cast<Cell2d*>(m);           
            return c->starting_point(sgn);
        }
    }
    return NULL; 
}

TMPL Cell2d*
Cell2dList::neighbor( Point* p)
{
    //Cell2d* cl;
      foreach( Cell2d *c, this->s_neighbors     )
        if ( c->intersections(2).member(p) )
          { 
              return c;
          }
      foreach( Cell2d *c, this->e_neighbors     )
        if ( c->intersections(3).member(p) )
          { 
              return c;
          }
      foreach( Cell2d *c, this->n_neighbors     )
        if ( c->intersections(0).member(p) )
          { 
              return c;
          }
      foreach( Cell2d *c, this->w_neighbors     )
        if ( c->intersections(1).member(p) )
          { 
              return c;
          }

      //std::cout<<"Point ("<<p->x()<<","<<p->y()<<") not found on neighbors of "<< this<<"("<<this->neighbors().size() <<")"<<std::endl; 
    return NULL;
}

TMPL Seq<typename mmx::shape::cell2d<C,V>::Point*>
Cell2dList::intersections() const{
    Seq<Point *> s,e,n,w,r;

    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 Cell2dList
# undef AlgebraicCurve
# undef Intersection2dFactory
# undef Shape
# endif // shape_cell2d_List_hpp