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

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * cell2d_voronoi_site2d.hpp
 * 2010-04-05
 * Angelos Mantzaflaris
 *****************************************************************************
 *               Copyright (C) 2010 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/

# ifndef shape_cell2d_voronoi_site2d_hpp
# define shape_cell2d_voronoi_site2d_hpp

# include <realroot/Interval.hpp>
# include <realroot/polynomial_bernstein.hpp>
# include <realroot/Seq.hpp>
# include <shape/point.hpp>
# include <shape/solver_implicit.hpp>
# include <shape/cell2d.hpp>
# include <shape/voronoi_site2d.hpp>

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

# define TMPL template<class C, class V>
# define VoronoiSite2d voronoi_site2d<C,V>
# define SELF cell2d_voronoi_site2d<C,V>
# define Cell2dAlgebraicCurve cell2d_algebraic_curve<C,V>

# undef Cell_t
//====================================================================
namespace mmx {
namespace shape {
//====================================================================
// TMPL struct cell2d_voronoi_site2d;
// //--------------------------------------------------------------------
// TMPL struct with_cell2d_voronoi_site2d { 
//   typedef cell2d_voronoi_site2d<C,V> Cell2dVoronoiSite2d; 
// };
//--------------------------------------------------------------------
struct cell2d_voronoi_site2d_def {};
template<> struct use<cell2d_voronoi_site2d_def> 
//  :public use<cell2d_algebraic_curve_def> 
//   ,public with_voronoi_site2d<K>
//   ,public with_algebraic_set<K>
//   ,public with_cell2d_voronoi_site2d<K>
{
  typedef bounding_box<double,double> Bounding_Box;
  //  typedef typename cell2d_def<K>::Point Vector;
  //  typedef polynomial< Interval<double>, with<Bernstein> > Polynomial;
  typedef polynomial< double, with<Bernstein> > Polynomial;
//    typedef polynomial< mmx::GMP::rational, with<Bernstein> > Polynomial;
  
  typedef solver_implicit<double,double> Solver;
  
//   template<class CELL> static void 
//   subdivisor(CELL & cl, CELL*&, CELL*&, int v, double s);

};

//--------------------------------------------------------------------
// TMPL1
// template<class CELL> void 
// use<cell2d_voronoi_site2d_def,V>::subdivisor(CELL& cl, CELL*& left, CELL*& right, int v, double s)  
// {
//   use<cell2d_algebraic_curve_def,V>::subdivisor(cl,left,right,v,s);
// }

//====================================================================
TMPL struct cell2d_voronoi_site2d : public  cell2d_algebraic_curve<C,V> {
public:
  typedef typename topology2d<C,V>::Point          Point;
  // typedef typename cell2d_voronoi_site2d_def<C,V>::Vertex         Vertex;
  typedef point<C,REF> Vector; 
  //typedef typename use<cell2d_voronoi_site2d_def,V>::Vector         Vector;
  typedef typename topology2d<C,V>::Edge                          Edge;
  typedef cell<C,REF>           Cell;
  //typedef typename cell2d_voronoi_site2d_def<C,V>::AlgebraicSet   AlgebraicSet;
  //typedef typename cell2d_voronoi_site2d_def<C,V>::VoronoiSite2d VoronoiSite2d;
  typedef typename cell2d_voronoi_site2d<C,V>::Polynomial     Polynomial;
  typedef bounding_box<C,REF>  BoundingBox;
  //  typedef typename use<cell2d_voronoi_site2d_def,V>::BoundingBox    BoundingBox;
  //  typedef typename use<cell2d_voronoi_site2d_def,V>::Topology       Topology;
  typedef typename use<cell2d_voronoi_site2d_def,V>::Solver         Solver;
  typedef topology<C,V>     Topology;
  typedef topology2d<C,V>     Topology2d;

  cell2d_voronoi_site2d(const Polynomial &, const BoundingBox&);
  cell2d_voronoi_site2d(const char*, const BoundingBox&);
//  cell2d_voronoi_site2d(const typename AlgebraicSet::Polynomial& , const BoundingBox&);
  cell2d_voronoi_site2d(VoronoiSite2d*, const BoundingBox &);
  cell2d_voronoi_site2d(const SELF& cl)
  : Cell2dAlgebraicCurve(cl) {}

  Polynomial equation(void) { return  this->m_polynomial; }
  
  bool  is_active (void) ;
  bool  is_regular(void) ;
  bool  is_intersected(void) {return (this->nb_intersect()>0);};

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

  virtual bool insert_regular (Topology*);// = 0;
  virtual bool insert_singular(Topology*);// = 0;
  virtual void subdivide(Cell*& left, Cell*& right, int v, double s);
  Vector gradient(const Point& p);

  inline double lower();
  inline double upper();

  void compute_boundary();

private:

//  double upp, low; //upper and lower bound on the value in the domain
};

//====================================================================
  TMPL inline double
  SELF::lower()
  {
    return  double(*std::min_element( this->m_polynomial.begin(), 
                                         this->m_polynomial.end() ) );
  }
//====================================================================
  TMPL inline double
  SELF::upper()
  {
    return double(*std::max_element( this->m_polynomial.begin(),
                                     this->m_polynomial.end() ) );
  }

//====================================================================
  TMPL void 
  SELF::compute_boundary()
  {
    
    //   BoundingBox * b( (BoundingBox*)this );
    //Seq<mmx::GMP::rational> bx;
    //bx<<as<mmx::GMP::rational>( ((BoundingBox*)this)->xmin());
    //bx<<as<mmx::GMP::rational>(this->xmax());
    //bx<<as<mmx::GMP::rational>(this->ymin());
    //bx<<as<mmx::GMP::rational>(this->ymax());
    
    //Solver::edge_point(this->n_intersections, this->m_polynomial, Solver::north_edge, b);
    //Solver::edge_point(this->s_intersections, this->m_polynomial, Solver::south_edge, b);
    //Solver::edge_point(this->w_intersections, this->m_polynomial, Solver::west_edge , b);
    //Solver::edge_point(this->e_intersections, this->m_polynomial, Solver::east_edge , b);
    // Solver::extremal(this->m_singular, this->m_polynomial, b); 

  }

//====================================================================
TMPL inline void 
print(SELF* cv) {
  typedef typename SELF::Point Point;
  std::cout << "point(s) b: "; 
  foreach(Point* p, cv->n_intersections)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  foreach(Point* p, cv->s_intersections)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  foreach(Point* p, cv->w_intersections)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  foreach(Point* p, cv->e_intersections)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  std::cout<<std::endl;
  
  std::cout << "Point(s) s: "; 
  foreach(Point* p, cv->m_singular)  
    std::cout <<" ["<<p->x()<<" "<<p->y()<<"]";
  std::cout<<std::endl;
}
//--------------------------------------------------------------------
TMPL
SELF::cell2d_voronoi_site2d(const Polynomial& pol, const BoundingBox& bx)
  : Cell2dAlgebraicCurve(pol,bx)
{
  //low=upp=NULL;
}

//--------------------------------------------------------------------
TMPL
SELF::cell2d_voronoi_site2d(const char* s, const BoundingBox & b)
  : Cell2dAlgebraicCurve(s,b)
{
  //low=upp=NULL;
}
//--------------------------------------------------------------------  
TMPL
SELF::cell2d_voronoi_site2d(VoronoiSite2d* cv, const BoundingBox & b)
  : Cell2dAlgebraicCurve(cv->distfunc(),b)
{
  //low=upp=NULL;
  //std::cout<<"created voronoi cell" <<std::endl;    

}
//--------------------------------------------------------------------
TMPL bool 
SELF::is_regular() {

    return true;
}
//--------------------------------------------------------------------
TMPL bool 
SELF::is_active() {

    return ( false );

}
//--------------------------------------------------------------------
TMPL bool 
SELF::insert_regular(Topology* t) {

  //std::cout<<"regular voronoi cell" <<std::endl;    
    t->insert((BoundingBox*)this,false);
    return true;

}
//--------------------------------------------------------------------
TMPL bool 
SELF::insert_singular(Topology* t) {

    return true;
  }

//--------------------------------------------------------------------
TMPL void
SELF::subdivide(Cell*& Left, Cell*& Right, int v, double s) {
  cell2d_subdivisor<C,V>::subdivide( *this,(SELF*&)Left,(SELF*&)Right,v,s);
}
//--------------------------------------------------------------------
TMPL  typename SELF::Vector
SELF::gradient(const Point& p) {
  Polynomial 
    dx= diff(this->m_polynomial,0),
    dy= diff(this->m_polynomial,1);
  double 
    u0= (p[0]-this->xmin())/(this->xmax()-this->xmin()),
    u1= (p[1]-this->ymin())/(this->ymax()-this->ymin());
  Vector v;
  v[0] = dx(u0,u1);
  v[1] = dy(u0,u1);
  return v;
}

//--------------------------------------------------------------------
TMPL typename SELF::Point*
SELF::starting_point(int sgn) {//

//std::cout<<"startingPoint.site2d"<<std::endl;

      Seq<Point*> all;
      unsigned a;
      all = this->intersections();
      a   = all.size();

      //std::cout<<this<<" , "<<a<<std::endl;
      //foreach(Point*p,all)
      //  std::cout<<p->x()<<" "<<p->y()<<" "<<p->z()<<std::endl;

      assert(a==2);

      int ev(0);
      int u ;//position of p
      unsigned
          s=this->s_intersections.size() ,// ~0
          e=this->e_intersections.size() ,// ~1
          //n=this->n_intersections.size() ,// ~2
          w=this->w_intersections.size() ;// ~3

      u= ( 0<s   ? 0 :
         ( 0<s+e ? 1 :
         ( 0<a-w ? 2 :
                   3 )));

      int * sz = this->m_polynomial.rep().szs();
      int * st = this->m_polynomial.rep().str();

          switch (u){
          case 0:
              ev= (this->m_polynomial[0] >0 ? 1:-1);
         break;
          case 1:
              ev= (this->m_polynomial[(sz[0]-1)*st[0]] >0 ? 1:-1);      
         break;
          case 2:
              ev= (this->m_polynomial[sz[0]*sz[1]-1]>0 ? 1:-1);
          break;
          case 3:
              ev= (this->m_polynomial[(sz[1]-1)*st[1]] >0 ? 1:-1);
          break;
          }

          if (ev*sgn>0)  
            return all[0]; 
          else                
            return all[1];
}

//--------------------------------------------------------------------
TMPL typename SELF::Point*
SELF::pair( typename SELF::Point* p, int& sgn) {
     // Pair: returns a neighboring point on the component (p,sgn)
     // The candidates are the 2 neighbor intersections of p in *this

      Seq<Point*> all;
      int a;
      all= this->intersections();
      a   = all.size();
      if (a==2)// regular cell with a single branch
         { return (all[0]==p ? all[1]: all[0]);  }

        std::cout<<"...pair Trouble"<< this<<std::endl;
        return (all[0]); 
}

//====================================================================
} ; // namespace shape
} ; // namespace mmx
# undef Solver
# undef SELF
# undef TMPL
# undef Cell2dAlgebraicCurve
# undef STACK
# undef VoronoiSite2d
# endif // shape_cell2d_voronoi_site2d_hpp