/Users/mourrain/Devel/mmx/shape/include/shape/cell2d_intersection.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_intersection_hpp
# define shape_cell2d_intersection_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 Cell2dInter cell2d_intersection<C,V>
namespace mmx {
namespace shape {

TMPL

class cell2d_intersection: public cell2d_list<C,V>   {
public:

  typedef typename Cell::BoundingBox  BoundingBox;
  typedef typename cell2d<C,V>::Point   Point;
  typedef topology<C,V>                 Topology;

  cell2d_intersection(void) ;
  cell2d_intersection(double xmin, double xmax) ;
  cell2d_intersection(double xmin, double xmax, double ymin, double ymax) ;
  cell2d_intersection(double xmin, double xmax, double ymin, double ymax, bool itr) ;
  cell2d_intersection(const BoundingBox& bx);
  virtual ~cell2d_intersection(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);

  bool is_touching(void);

  virtual unsigned nb_intersect(void) const;

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

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

} ;

TMPL
Cell2dInter::cell2d_intersection(void) {}

TMPL
Cell2dInter::cell2d_intersection(double xmin, double xmax, double ymin, double ymax) :  cell2d_list<C,V>(xmin, xmax, ymin, ymax) {}

TMPL
Cell2dInter::cell2d_intersection(double xmin, double xmax, double ymin, double ymax, bool itr) : cell2d_list<C,V>(xmin, xmax, ymin, ymax) {}

TMPL
Cell2dInter::cell2d_intersection(const BoundingBox& bx): cell2d_list<C,V>(bx) {};

TMPL
Cell2dInter::~cell2d_intersection(void) {
    foreach (Cell* m, this->m_objects) delete m;
}


TMPL bool 
Cell2dInter::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 
Cell2dInter::is_intersected() {
  
  if(this->m_objects.size() >1 && !this->m_intersected) {
    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);
    this->m_intersected = true;
  }
  
  if (this->m_singular.size() > 0) return true;
  return false;
}


TMPL bool 
Cell2dInter::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 
Cell2dInter::insert_singular(Topology* s) {
  //s->singular(this);
//  std::cout<<"Ins Singular"<<std::endl;

  foreach(Cell * m, this->m_objects) m->insert_singular(s);
  return true;
}
TMPL void
Cell2dInter::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
Cell2dInter::subdivide(Cell *& left, Cell *& right, int  v, double c) {

  typedef Cell2dInter Cell_t;

  if(v==1) {
    left =(Cell*)new Cell_t(this->xmin(), this->xmax(), this->ymin(), c, this->m_intersected) ;
    right=(Cell*)new Cell_t(this->xmin(), this->xmax(), c, this->ymax(), this->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 {
    left = (Cell*)new Cell_t(this->xmin(), c, this->ymin(), this->ymax(), this->m_intersected) ;
    right= (Cell*)new Cell_t(c, this->xmax(), this->ymin(), this->ymax(), this->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( );

  //Note: for a semi-alg. set, we keep all the 
  //equations in the cell, even if inactive.

  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 
Cell2dInter::is_active() {

  //return ( cell2d_list<C,V>::is_active() );

  if ( !cell2d_list<C,V>::is_active() )
    return false;

  if( this->is_intersected() ) return true;
  if( ! this->is_regular() ) return true;

  if( this->is_touching() )//is touching the set?
    return true;
  else
    return false;
}


TMPL  bool 
Cell2dInter::is_touching() {

  int c0(0),c1(0),c2(0),c3(0), cnt=this->count();//cnt: number of eqs.

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

    if (c->m_polynomial[0] >0) c0++;
    if (c->m_polynomial[(sz[0]-1)*st[0]] >0) c1++;
    if (c->m_polynomial[sz[0]*sz[1]-1]>0) c2++;
    if (c->m_polynomial[(sz[1]-1)*st[1]] >0) c3++;
  }
  return ( (c0==cnt) || (c1==cnt) || (c2==cnt) || (c3==cnt) );
}
 

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



TMPL typename Cell2dInter::Point* 
Cell2dInter::pair(Point * p, int sgn)
{
    if ( this->is_intersected() ) {
        //
//         std::cout<<"Reached intersection cell "<<this<<std::endl;
//         std::cout<<"p= "<< p->x()<<","<<p->y()<<std::endl;
//         foreach (Point*q, this->intersections() ) 
//             std::cout<<"q= "<< q->x()<<","<<q->y()<<std::endl;

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

        foreach (Cell* m, this->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->x()<<","<<q->y()<<std::endl;

        int ev2(444);
        foreach (Cell* m, this->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, this->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 Cell2dInter::Point* 
Cell2dInter::starting_point( int sgn)
{
   Cell2d* c;
   foreach (Cell* m, this->m_objects)
    { 
        if ( m->is_active() )
        {
            c = dynamic_cast<Cell2d*>(m);           
            return c->starting_point(sgn);
        }
    }
    return NULL; 
}

TMPL Cell2d*
Cell2dInter::neighbor( Point* p)
{

      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*>
Cell2dInter::intersections() const{
    Seq<Point *> s,e,n,w,r;

    Cell2d* cl;
    foreach (Cell* m, this->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 Cell2dInter
# undef AlgebraicCurve
# undef Intersection2dFactory
# undef Cell2dAlgebraicCurve
# undef Shape
# endif // shape_cell2d_intersection_hpp