include/shape/voronoi_site2d.hpp

Go to the documentation of this file.

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

# ifndef shape_voronoi_site2d_hpp
# define shape_voronoi_site2d_hpp

# include <realroot/GMP.hpp>
# include <shape/point.hpp>
# include <realroot/polynomial.hpp>
# include <realroot/ring_sparse.hpp>

#include <realroot/solver_uv_sleeve.hpp>
# include <shape/algebraic_set.hpp>
# include <shape/viewer_axel.hpp>

# define TMPL template<class C, class V>
# define SELF voronoi_site2d<C,V>
# define REF REF_OF(V)
# define Point point<double>
# define Viewer viewer<V>
//====================================================================
namespace mmx {
namespace shape {
//====================================================================

template<class C, class V=default_env>
class voronoi_site2d 
{

public:

  typedef typename algebraic_set<C,V>::Polynomial Polynomial;
  //typedef polynomial<C, with<Sparse, DegRevLex> > Polynomial;

  typedef bounding_box<C,REF>  BoundingBox;

  voronoi_site2d(void) {} ;
  voronoi_site2d(const Point & p, const int& k= 2);
  voronoi_site2d(const Point & p, double a, double b,double c);
  voronoi_site2d(const double & w, const Point & p);
  voronoi_site2d(const Point & p, const char * s);
  voronoi_site2d(const Point & p, const Polynomial &) ;

  voronoi_site2d(const Point & p, const double & r) ;

  ~voronoi_site2d(void) {};

  Polynomial distfunc() const  { return this->func ; }

  Point      coordinates() const { return this->coords ; }
  double      x() const { return this->coords.x() ; }
  double      y() const { return this->coords.y() ; } 
  double      z() const { return this->coords.z() ; }

  double      upper_bound(BoundingBox& bx) const; 
  Polynomial  offset( const C& z){ return (this->func - Polynomial(z,0,0) ); }

//Polynomial lp(const int& n, const Point& p );

Polynomial lp(const int& n, const Point& p )
{
//  typedef typename SELF::Polynomial Polynomial;

  Polynomial a= Polynomial(1,0,0);

  C c= C(1);
  if( p.x()!=0.0)
    for(int i=0; i<=n; i++ ) {
      a+= Polynomial( c*binomial<C>(n,i),n-i,0);
      c*= -p.x(); }
  else a+= Polynomial(1,n,0);
  c=1.0;
  if( p.y()!=0.0)
    for(int i=0; i<=n; i++ ) {
      a+= Polynomial( c*binomial<double>(n,i),n-i,1);
      c*= -p.y(); }
  else a+= Polynomial(1,n,1);

  a-=C(1);
  return a;
}


private:

  Point coords;
  Polynomial func;

} ;


// Point and custom distance field
TMPL
SELF::voronoi_site2d(const Point & p, const char * s)
 { 
   // typedef typename use<numeric_def,V>::Rational Coefficients;
  polynomial<double, with<Sparse,DegRevLex> >t(s, variables("x y"));
  Polynomial e; let::assign(func,t);
  // func= Polynomial(s, variables("x y"));

  //std::cout<< func; 
  coords= p; 
}

// Anisotropic site 
// with quad. form given  by matrix M*tr(M) where
// M = [ sqrt(a)  0     ] 
//     [   c    sqrt(b) ]
// i.e. symmetric and positive definite, distance= tr(p)*M*tr(M)*p
TMPL
SELF::voronoi_site2d(const Point & p, double a,double b,double c)
{
  //std::cout<<"seed "<< a<<" "<< b <<" "<< c <<std::endl;

  func= Polynomial(a,2,0)+Polynomial(b+c*c,2,1) + 

    ( c!=0? Polynomial(2*c*std::sqrt(a),1,0)*Polynomial(1,1,1):0) + 
    Polynomial(-2*std::sqrt(a)*c*p.x()-2*p.y()*c*c-2*p.y()*b ,1,1) -
    (p.x()!=0? Polynomial(2*a*p.x(), 1 , 0): 0) - 
    ( c!=0? Polynomial(2*std::sqrt(a)*c*p.y(), 1, 0 ): 0 ) + 
    Polynomial( p.y()*p.y()*b+a*p.x()*p.x()+p.y()*p.y()*c*c+
                2*std::sqrt(a)*c*p.y()*p.x(),0,0)  ;


  std::cout<<"site "<< func <<std::endl;
  coords= p; 

}

// Classic point-site
TMPL
SELF::voronoi_site2d(const Point & p, const int& k)
{ 
  func= this->lp(k , p );
  coords= p;

}

// Point-site with weight (power or Laguerre diagram)
TMPL
SELF::voronoi_site2d(const double & w, const Point & p)
{ 
  func= Polynomial(1,2,0)+Polynomial(1,2,1);

  func += Polynomial(-2*p.x(),1,0) // (coeff,degree,variable)
    +  Polynomial(-2*p.y(),1,1)
    +  Polynomial(p.x()*p.x()+p.y()*p.y(),0,0)
    -  Polynomial( w*w ,0,0);//adding weight 

  coords= p;
  std::cout<<"Point: "<< p<<", w: "<< w <<std::endl;
}

// Point and costum distance field
TMPL
SELF::voronoi_site2d(const Point& p, const Polynomial& f) 
    {
      func  = f;
      coords= p; 
    }

// VD of disks, (additively weighted or Apollonius diagram)
TMPL
SELF::voronoi_site2d(const Point & p, const double & r)
{
  func= Polynomial(1,2,0)+Polynomial(1,2,1) - Polynomial(1,2,2);
  
  func += Polynomial(-2*p.x(),1,0) // (coeff,degree,variable)
    +  Polynomial(-2*p.y(),1,1)
    +  Polynomial(p.x()*p.x()+p.y()*p.y() -r*r ,0,0)
    -  Polynomial( 2*r ,1,2);
  
      coords= p; 

      //BoundingBox bx(0,1,0,1);
      //std::cout<<"upp: "<< this->upper_bound(bx) <<std::endl;


      //std::cout<< func <<"\n";
  }

TMPL double
SELF::upper_bound(BoundingBox& bx) const
{
  //double ub;

  //AkritasBound<C> lb;


    Interval<double> ev(0.3434,0.3434);
    Seq< Interval<double> > bb;
    bb<< Interval<double>(bx.xmin(), bx.xmax() ) ;    
    bb<< Interval<double>(bx.ymin(), bx.ymax() ) ;    
    polynomial<Interval<double>,with<Sparse, DegRevLex> > pol(1.4,1,1);
    polynomial<double, with<Sparse, DegRevLex> > upol;
    let::assign(pol,func);
    //std::cout<<"  ev= "<< ev<<"\n";
    //std::cout<<"  bb= "<< bb<<"\n";

  if ( func.nbvar()==2 )
  {
    sparse::eval_at( ev , pol.rep() ,bb );
    //tensor::eval( ev , func.rep(),bb , 2 );


    //std::cout<<"pol= "<< pol<<"\n";


    return ev.upper();
  }
  else  // ( func.nbvar()==3 )
  {



    pol= sparse::subs(pol, 0, bb[0] ) ;
    std::cout<<"sub= " << pol <<"\n";
    pol= sparse::subs(pol, 1, bb[1] ) ;
    std::cout<<"sub= " << pol <<"\n";

    
    typedef typename mmx::polynomial<double, with<Sparse, DegRevLex> >::Monomial  mon;

        for ( typename mmx::polynomial<Interval<double>,with<Sparse, DegRevLex> >::const_iterator it =pol.begin();   it != pol.end(); ++it )
                  upol+= mon(it->coeff().upper(), it->rep() ) ;

    
    std::cout<<"upol=" << upol <<"\n";

    //std::cout<<"upol=" << mon() <<"\n";

    //Cauchy<double> bound;
    //solve(f,Sleeve<Approximate>())

    return 1;// bound.upper_bound(upol);
  }
}

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


# undef TMPL 
# undef SELF
# undef REF
# undef Point
# undef Viewer
# endif // shape_voronoi_site2d_hpp

// 
// voronoi_site2d.hpp ends here

---