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

Go to the documentation of this file.

/*****************************************************************************
 * M a t h e m a g i x
 *****************************************************************************
 * RationalCurve
 * 2010-04-03
 *****************************************************************************
 *               Copyright (C) 2006 INRIA Sophia-Antipolis
 *****************************************************************************
 * Comments :
 ****************************************************************************/
# ifndef shape_curve_rational_hpp
# define shape_curve_rational_hpp

# include <realroot/GMP.hpp>
# include <realroot/polynomial_bernstein.hpp>
# include <shape/graphic.hpp>
# include <shape/point.hpp>
# include <shape/parametric_curve.hpp>

# define TMPL_DEF template<class C, class V=default_env>
# define TMPL template<class C, class V>
# define TMPL1 template<class V>
# define REF     REF_OF(V)
# define SELF curve_rational<C,V>
# define GRAPHIC graphic<C,V>
# define AXEL viewer<axel,V>
//====================================================================
namespace mmx {
namespace shape {

struct curve_rational_def {};

template<> struct use<curve_rational_def> {
  typedef double Range;
  typedef polynomial< double, with<Bernstein> > Polynomial;
};

TMPL_DEF
class curve_rational : public parametric_curve<C,V> {
public:
  typedef bounding_box<C,V>   BoundingBox;
  typedef parametric_curve<C,V> ParametricCurve;
  typedef typename ParametricCurve::Scalar Scalar;
  typedef typename ParametricCurve::Point Point;
  typedef typename use<curve_rational_def,V>::Polynomial Polynomial;

  curve_rational(void) :ParametricCurve(), m_tmin(0), m_tmax(1) {}

  curve_rational(const BoundingBox & box): 
    ParametricCurve(box), m_tmin(0), m_tmax(1) {}
  
  curve_rational(const Polynomial& X, const Polynomial& Y, const Polynomial& Z, const Polynomial& W=1)
    : ParametricCurve(), m_tmin(0), m_tmax(1), m_w(W) 
  {
    m_p<<X; m_p<<Y; m_p<<Z;
  }

  curve_rational(double m, double M, const Polynomial& X, const Polynomial& Y, const Polynomial& Z, const Polynomial& W=1)
    : ParametricCurve(), m_tmin(m), m_tmax(M), m_w(W) 
  {
    m_p<<X; m_p<<Y; m_p<<Z;
  }

  curve_rational(const SELF& c)
    : ParametricCurve(), m_tmin(c.tmin()), m_tmax(c.tmax()),  m_w(c.m_w) 
  {
    m_p=c.m_p; 
  }

  int dimension() const { return m_p.size();}

 ~curve_rational(void) {};

  double tmin(void) const { return this->m_tmin; }
  double tmax(void) const { return this->m_tmax; }
  void   set_range(double tmin, double tmax);

  Point* eval (const double& t) const;
  Point* operator() (double t);
  void   eval   (double t, double * x, double * y, double * z)  const ;
  void   eval   (Point&,  Scalar ) const;

  void   subdivide (ParametricCurve*& a, ParametricCurve*& b, double t ) const ;

  const Polynomial& denominator()    const { return m_w; }
  const Polynomial& numerator(int i) const { return m_p[i]; }

  Seq<Point *> critical_points(void) ;
  Seq<Point *> extremal_points(void) ;
  Seq<Point *> singular_points(void) ;
  
private:
  double m_tmin ;
  double m_tmax ;
  
  Polynomial m_w;
  Seq<Polynomial> m_p;
};

TMPL typename SELF::Point* 
SELF::eval(const double& t) const {
  Scalar v=as<Scalar>(t);
  double w = as<double>(m_w(v));
  Point* p = new Point(m_p[0](v)/w, m_p[1](v)/w, m_p[2](v)/w);
  return p ;
}

TMPL void 
SELF::eval(Point& v, Scalar p) const {
  Scalar w=m_w(p);
  v.setx(m_p[0](p)/w);
  v.sety(m_p[1](p)/w);
  v.setz(m_p[2](p)/w);
}

TMPL typename SELF::Point* 
SELF::operator() (double t) { return eval(t); }


TMPL void 
SELF::set_range(double tmin, double tmax) {
  m_tmin=tmin;                                           
  m_tmax=tmax;                                           
}

  TMPL Seq<typename SELF::Point *> 
SELF::critical_points(void) {
  Seq<Point *> l ;
  // ...
  return l ;
}

TMPL Seq<typename SELF::Point *> 
SELF::extremal_points(void) {
   Seq<Point *> l ;
   // ...
   return l ;
}

TMPL Seq<typename SELF::Point *> 
SELF::singular_points(void)
{
   Seq<Point *> l ;
   // ...
   return l ;
}
TMPL void
SELF::subdivide(ParametricCurve*& a, ParametricCurve*& b, double t ) const {
  SELF
    * ar = new SELF(*this),
    * br = new SELF(*this);
  ar->set_range(tmin(),t);
  br->set_range(t,tmax());

  a = dynamic_cast<ParametricCurve*>(ar);
  b = dynamic_cast<ParametricCurve*>(br);

}

TMPL struct viewer; struct axel;
TMPL AXEL&
operator<<(AXEL& os, const SELF& c) {
  os<<"\n <curve type=\"rational\"  color=\""<<(int)(255*os.color.r)<<" "<<(int)(255*os.color.g)<<" "<<(int)(255*os.color.b)<<"\">\n";
  os<<"   <domain>"<< c.tmin()<<" "<<c.tmax()<<"</domain>\n";
  for(int i=0; i<c.dimension();i++){
    os<<"   <polynomial>";
    print(os,c.numerator(i),variables("t"));
    os<<"</polynomial>\n";
  }
  os<<"   <polynomial>";
  print(os,c.denominator(),variables("t"));
  os<<"</polynomial>\n";
  os<<" </curve>\n";
  return os;
}
//--------------------------------------------------------------------
TMPL GRAPHIC*
as_graphic(const SELF& c, unsigned N= 100) {
  typedef typename curve_pl<C,V>::Point Point;
  typedef typename curve_pl<C,V>::Edge  Edge;
  typedef C                    Scalar;


  GRAPHIC* r = new GRAPHIC( GRAPHIC::E_LINE, N+1, 2*N);
  
  Scalar t=c.tmin(),dt=(c.tmax()-c.tmin())/N;
  for(unsigned i=0;i<=N;i++) {
    Scalar w=c.denominator()(t);
    for(unsigned j=0;j<3;j++) 
      r->vertices[3*i+j]=c.numerator(j)(t)/w;
    t+=dt;
  }

  for(unsigned i=0;i<N;i++) {
    r->indices[2*i]=i;
    r->indices[2*i+1]=i+1;
  }
  return r;
}

//====================================================================
} ; // namespace shape
//====================================================================
} ; // namespace mmx
//====================================================================
# undef TMPL
# undef TMPL1
# undef TMPL_DEF
# undef REF
# undef Scalar
# undef Point
# undef SELF
# undef ParametricCurve
# undef AXEL
# undef GRAPHIC
# endif // shape_curve_rational_hpp