00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 # ifndef shape_rational_curve_hpp
00012 # define shape_rational_curve_hpp
00013
00014 # include <realroot/GMP.hpp>
00015 # include <realroot/ring_monomial_tensor.hpp>
00016 # include <shape/graphic.hpp>
00017 # include <shape/point.hpp>
00018 # include <shape/parametric_curve.hpp>
00019
00020 # define TMPL_DEF template<class C, class V=default_env>
00021 # define TMPL template<class C, class V>
00022 # define TMPL1 template<class V>
00023 # define REF REF_OF(V)
00024 # define SELF rational_curve<C,V>
00025 # define GRAPHIC graphic<C,V>
00026 # define AXEL viewer<axel,V>
00027
00028 namespace mmx {
00029 namespace shape {
00030
00031 struct rational_curve_def {};
00032
00037 template<class C> struct use<rational_curve_def, C> {
00038 typedef C Range;
00039 typedef polynomial< C, with<MonomialTensor> > Polynomial;
00040 };
00041
00044 TMPL_DEF
00045 class rational_curve : public parametric_curve<C,V> {
00046 public:
00047 typedef bounding_box<C,V> BoundingBox;
00048 typedef parametric_curve<C,V> ParametricCurve;
00049 typedef typename ParametricCurve::Scalar Scalar;
00050 typedef typename ParametricCurve::Point Point;
00051 typedef typename use<rational_curve_def,C,V>::Polynomial Polynomial;
00052
00053 rational_curve(void) :ParametricCurve(), m_tmin(0), m_tmax(1) {}
00054
00055 rational_curve(const BoundingBox & box):
00056 ParametricCurve(box), m_tmin(0), m_tmax(1) {}
00057
00058 rational_curve(const Polynomial& X, const Polynomial& Y, const Polynomial& Z, const Polynomial& W=1)
00059 : ParametricCurve(), m_tmin(0), m_tmax(1), m_w(W)
00060 {
00061 m_p<<X; m_p<<Y; m_p<<Z;
00062 }
00063
00064 rational_curve(double m, double M, const Polynomial& X, const Polynomial& Y, const Polynomial& Z, const Polynomial& W=1)
00065 : ParametricCurve(), m_tmin(m), m_tmax(M), m_w(W)
00066 {
00067 m_p<<X; m_p<<Y; m_p<<Z;
00068 }
00069
00070 rational_curve(const SELF& c)
00071 : ParametricCurve(), m_tmin(c.tmin()), m_tmax(c.tmax()), m_w(c.m_w)
00072 {
00073 m_p=c.m_p;
00074 }
00075
00076 int dimension() const { return m_p.size();}
00077
00078 ~rational_curve(void) {};
00079
00080 double tmin(void) const { return this->m_tmin; }
00081 double tmax(void) const { return this->m_tmax; }
00082 void set_range(double tmin, double tmax);
00083
00084 Point* eval (const double& t) const;
00085 Point* operator() (double t);
00086 void eval (double t, double * x, double * y, double * z) const ;
00087 void eval (Point&, Scalar ) const;
00088
00089 void subdivide (ParametricCurve*& a, ParametricCurve*& b, double t ) const ;
00090
00091 const Polynomial& denominator() const { return m_w; }
00092 const Polynomial& numerator(int i) const { return m_p[i]; }
00093
00094 void set_denominator(const Polynomial& d) { m_w=d; }
00095 void set_numerator (const Polynomial& p, int i) { if (m_p.size() > i) m_p[i]=p; else m_p.push_back(p); }
00096
00097 Seq<Point *> critical_points(void) ;
00098 Seq<Point *> extremal_points(void) ;
00099 Seq<Point *> singular_points(void) ;
00100
00101 private:
00102 double m_tmin ;
00103 double m_tmax ;
00104
00105 Polynomial m_w;
00106 Seq<Polynomial> m_p;
00107 };
00108
00109 TMPL typename SELF::Point*
00110 SELF::eval(const double& t) const {
00111 Scalar v=as<Scalar>(t);
00112 double w = as<double>(m_w(v));
00113 Point* p = new Point(m_p[0](v)/w, m_p[1](v)/w, m_p[2](v)/w);
00114 return p ;
00115 }
00116
00117 TMPL void
00118 SELF::eval(Point& v, Scalar p) const {
00119 Scalar w=m_w(p);
00120 v.setx(m_p[0](p)/w);
00121 v.sety(m_p[1](p)/w);
00122 v.setz(m_p[2](p)/w);
00123 }
00124
00125 TMPL typename SELF::Point*
00126 SELF::operator() (double t) { return eval(t); }
00127
00128
00129 TMPL void
00130 SELF::set_range(double tmin, double tmax) {
00131 m_tmin=tmin;
00132 m_tmax=tmax;
00133 }
00134
00135 TMPL Seq<typename SELF::Point *>
00136 SELF::critical_points(void) {
00137 Seq<Point *> l ;
00138
00139 return l ;
00140 }
00141
00142 TMPL Seq<typename SELF::Point *>
00143 SELF::extremal_points(void) {
00144 Seq<Point *> l ;
00145
00146 return l ;
00147 }
00148
00149 TMPL Seq<typename SELF::Point *>
00150 SELF::singular_points(void)
00151 {
00152 Seq<Point *> l ;
00153
00154 return l ;
00155 }
00156 TMPL void
00157 SELF::subdivide(ParametricCurve*& a, ParametricCurve*& b, double t ) const {
00158 SELF
00159 * ar = new SELF(*this),
00160 * br = new SELF(*this);
00161 ar->set_range(tmin(),t);
00162 br->set_range(t,tmax());
00163
00164 a = dynamic_cast<ParametricCurve*>(ar);
00165 b = dynamic_cast<ParametricCurve*>(br);
00166
00167 }
00168
00169 TMPL struct viewer; struct axel;
00170 TMPL AXEL&
00171 operator<<(AXEL& os, const SELF& c) {
00172 os<<"\n <curve type=\"rational\" color=\""<<(int)(255*os.color.r)<<" "<<(int)(255*os.color.g)<<" "<<(int)(255*os.color.b)<<"\">\n";
00173 os<<" <domain>"<< c.tmin()<<" "<<c.tmax()<<"</domain>\n";
00174 for(int i=0; i<c.dimension();i++){
00175 os<<" <polynomial>";
00176 print(os,c.numerator(i),variables("t"));
00177 os<<"</polynomial>\n";
00178 }
00179 os<<" <polynomial>";
00180 print(os,c.denominator(),variables("t"));
00181 os<<"</polynomial>\n";
00182 os<<" </curve>\n";
00183 return os;
00184 }
00185
00186 TMPL GRAPHIC*
00187 as_graphic(const SELF& c, unsigned N= 100) {
00188 typedef typename curve_pl<C,V>::Point Point;
00189 typedef typename curve_pl<C,V>::Edge Edge;
00190 typedef C Scalar;
00191
00192
00193 GRAPHIC* r = new GRAPHIC(N+1, 2*N,0);
00194
00195 Scalar t=c.tmin(),dt=(c.tmax()-c.tmin())/N;
00196 for(unsigned i=0;i<=N;i++) {
00197 Scalar w=c.denominator()(t);
00198 for(unsigned j=0;j<3;j++)
00199 r->vertices[3*i+j]=c.numerator(j)(t)/w;
00200 t+=dt;
00201 }
00202
00203 for(unsigned i=0;i<N;i++) {
00204 r->edges[2*i]=i;
00205 r->edges[2*i+1]=i+1;
00206 }
00207 return r;
00208 }
00209
00210
00211 } ;
00212
00213 } ;
00214
00215 # undef TMPL
00216 # undef TMPL1
00217 # undef TMPL_DEF
00218 # undef REF
00219 # undef Scalar
00220 # undef Point
00221 # undef SELF
00222 # undef ParametricCurve
00223 # undef AXEL
00224 # undef GRAPHIC
00225 # endif // shape_rational_curve_hpp