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

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * TopologyCurve
 * 2008-05-16
 * Julien Wintz
 *****************************************************************************
 *               Copyright (C) 2006 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
# ifndef shape_semialgebraic2d_hpp
# define shape_semialgebraic2d_hpp

# include <shape/topology2d.hpp>


# include <stack>
# define  STACK std::stack<Cell2d*>

# define TMPL template<class C, class V>
# define Shape geometric<V>
# define AlgebraicCurve algebraic_curve<C,V>
# define SemiAlgebraicCurve semialgebraic_curve<C,V>
# define ParametricCurve parametric_curve<C,V>
# define Cell2dAlgebraicCurve cell2d_algebraic_curve<C,V>
# define Cell2dParametricCurve cell2d_parametric_curve<C,V>
# define Cell2dSemiAlgebraicCurve cell2d_semialgebraic_curve<C,V>
# define Cell2dInter cell2d_intersection<C,V>
# define SELF semialgebraic2d<C,V>
# undef Cell
//====================================================================
namespace mmx {
namespace shape {



template<class C, class V=default_env>
class semialgebraic2d: public topology2d<C,V> {

  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:
  
  semialgebraic2d(const BoundingBox& bx): topology2d<C,V>(bx){ };
  semialgebraic2d(Curve * curve) ;
  ~semialgebraic2d(void) {   delete this->m_tree ;  };

  void run(void)  ;
  
  //    virtual void insert(Point*);
  //    virtual void insert(Edge *);
  //    virtual void insert(Point * p1, Point * p2);

  virtual void insert_singular(Point*);
  
protected:
  bool is_regular (Cell * cell) ;
  bool insert_regular(Cell * cell) ;
  bool singularity(Cell * cell) ;
  bool   subdivide(Cell * cell, Node * node) ;
  
  
private:
  Kdtree         *  m_tree ;
  std::list<Node *> m_nodes ;
  
};
  
//--------------------------------------------------------------------

TMPL void 
SELF::run() {
//
     std::cout<<"Running subdivision" <<std::endl;

    topology2d<C,V>::run();

     std::cout<<"Computing Semi-algebraic set." <<std::endl;

     Seq<Cell2d*> nlist;
     this->b_leaves.dfs(nlist);// remove empty border cells
     Cell2d* pr;
     pr= nlist[nlist.size()-1];
     //if (0)
     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);
   assert( nlist.size()>1);

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

   // Start exploration
   foreach(Cell2d* cl, nlist)
     if ( cl->is_active()    &&
          //cl->is_regular()   &&
          cl->nb_intersect()==2 
          && !cl->is_border() )
       stack.push(cl);  

  sgn=1;
//if (0)
while ( !stack.empty() )
   {
   s= stack.top();
   aux=s->m_gnode->aux();

   //-- only (+) faces
   //thus if cell is visited once then it is removed
   if (aux!=0) 
   {
     stack.pop();
     continue;
   }

   // Recovering face (s,sgn)
   f= new Face();

   // Get starting point (CCW)
   //std::cout<<"Getting ("<<(sgn==1 ? "+": "-")
   //  <<") face starting at "<< *s<<  std::endl;
   p= s->starting_point(sgn);

   // Walking on CCW face
   p= s->pair(p,sgn);
   f->insert(p);
   b=s;

//int c(0); 
   do  {
//if (++c>320) {while(!stack.empty()) stack.pop(); exit(0);break;}
       //std::cout<<"Next"<<b <<std::endl;
       //if ( !b->is_corner())
       //std::cout<<"aux="<<b->m_gnode->aux()<<std::endl;

      if( this->m_leaves.member(b) ) {
      aux=b->m_gnode->aux();
      b->m_gnode->aux(aux + (sgn==1 ? 2:-1) ); }

      t= b->neighbor(p);

      if ( t==NULL)
      { // border cell reached

        //std::cout<<"Border cell "<< *b <<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) );
        
        b=this->b_leaves.next(b);

          //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
            //std::cout<<"Entered at "<< *b <<std::endl;
            p= b->starting_point(sgn);
            //std::cout<<"first point "<< *p <<std::endl;
            f->insert(p);
            p= b->pair(p,sgn);
            //std::cout<<"pair  point "<< *p <<std::endl;
            f->insert(p);
          }
      }
      else
      { // next normal cell
          b=t;
          if ( b->m_singular.size()==1 )
              f->insert(b->m_singular[0]);
          p= b->pair(p,sgn);
          f->insert(p);
      }

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

  //std::cout<<"Face Added" << std::endl;
  this->m_faces<< f;

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

   }// End exploration

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

}

TMPL bool 
SELF::is_regular(Cell * cl) {
  return cl->is_regular();
}

TMPL bool 
SELF::insert_regular(Cell * cl) {
  return cl->insert_regular(this);
}

TMPL bool 
SELF::singularity(Cell * cl) {
  //  this->insert((BoundingBox*)cl);
  return cl->insert_singular(this) ;
}

TMPL void 
SELF::insert_singular(Point * p) {
  this->m_specials<<p;
}

TMPL bool 
SELF::subdivide(Cell * cl, Node * node) {
  int v=0;
  
  Cell* left, * right;
  v = cl->subdivide(left, right) ;
  
  node->m_left  = new Node(node, left,  Node::LEFT,  v); 
  m_nodes << node->m_left ;
  node->m_right = new Node(node, right, Node::RIGHT, v); 
  m_nodes << node->m_right;

  return true ;
}


//====================================================================
} // namespace shape
} // namespace mmx
//====================================================================
# undef TMPL
# undef AlgebraicCurve
# undef SemiAlgebraicCurve
# undef Cell
# undef Cell2dAlgebraicCurve
# undef Cell2dParametricCurve
# undef Cell2dSemiAlgebraicCurve
# undef Cell2dInter
# undef Cell2dFactory
# undef ParametricCurve
# undef Curve
# undef SELF 
# undef Shape
# undef STACK
# endif