|
shape_doc 0.1
|
|
shape_doc 0.1
|
00001 /***************************************************************************** 00002 * M a t h e m a g i x 00003 ***************************************************************************** 00004 * AlgebraicCurve 00005 * Bernard Mourrain 00006 ***************************************************************************** 00007 * Copyright (C) 2009 INRIA Sophia-Antipolis 00008 ***************************************************************************** 00009 * Comments : 00010 ****************************************************************************/ 00011 # ifndef algebraic_curve_fcts_hpp 00012 # define algebraic_curve_fcts_hpp 00013 00014 # include <realroot/Seq.hpp> 00015 # include <realroot/polynomial_tensor.hpp> 00016 # include <realroot/polynomial_bernstein.hpp> 00017 # include <realroot/solver_uv_continued_fraction.hpp> 00018 # include <shape/point.hpp> 00019 # include <shape/vertex.hpp> 00020 # include <shape/solver_implicit.hpp> 00021 # include <shape/algebraic_curve.hpp> 00022 00023 # define TMPL template<class C,class V> 00024 # define SELF algebraic_curve<C,V> 00025 # define PointSet point_set<C,V> 00026 # define BoundingBox bounding_box<C,V> 00027 //==================================================================== 00028 namespace mmx { 00029 namespace shape { 00030 00031 TMPL Seq<typename topology<C,V>::Point *> 00032 extremal(const SELF& c, const BoundingBox& b) { 00033 typedef typename topology<C,V>::Point Point; 00034 Seq<Point*> res; 00035 solver_implicit<C,V>::extremal(res,c.equation(),b); 00036 return res; 00037 } 00038 00039 TMPL PointSet 00040 plot(const SELF& c, const BoundingBox& bx, int N= 500) { 00041 00042 typedef double Scalar; 00043 typedef typename PointSet::Point Point; 00044 typedef typename SELF::Polynomial Polynomial; 00045 typedef typename Polynomial::Scalar Rational; 00046 typedef solver<Rational, ContFrac<Approximate> > Solver; 00047 00048 Point pt(0,0,0); 00049 PointSet pts; 00050 00051 Rational a0(bx.xmin()), a1(bx.xmax()), b0(bx.ymin()), b1(bx.ymax()); 00052 Rational dx= (a1-a0)/N; 00053 Seq<Polynomial> sy = coefficients(c.equation(),1); 00054 int d = sy.size()-1; 00055 polynomial<Rational, with<MonomialTensor> > P(0,d); 00056 00057 for (int i=0; i<N;i++,a0+=dx) { 00058 pt[0] = as<Scalar>(a0); 00059 for (int u=0;u<=d;u++) P[u] = sy[u](a0); 00060 typename Solver::Solutions sol; Solver::solve(sol,P,b0,b1); 00061 for(unsigned u=0;u<sol.size();u++) { 00062 pt[1] = as<Scalar>((lower(sol[u])+upper(sol[u]))/2); 00063 pts<<pt; 00064 } 00065 } 00066 Rational dy= (b1-b0)/N; 00067 a0= Rational(bx.xmin()); 00068 Seq<Polynomial> sx = coefficients(c.equation(),0); 00069 d = sx.size()-1; 00070 polynomial<Rational, with<MonomialTensor> > Q(0,d); 00071 for (int i=0; i<N;i++,b0+=dy) { 00072 pt[1] = as<Scalar>(b0); 00073 for (int u=0;u<=d;u++) Q[u] = sx[u]((Rational)0,b0); 00074 typename Solver::Solutions sol; Solver::solve(sol,Q,a0,a1); 00075 for(unsigned u=0;u<sol.size();u++) { 00076 pt[0] = as<Scalar>((lower(sol[u])+upper(sol[u]))/2); 00077 pts<<pt; 00078 } 00079 } 00080 00081 return pts; 00082 } 00083 //==================================================================== 00084 } ; // namespace mmx 00085 } ; // namespace mmx 00086 //==================================================================== 00087 # undef TMPL 00088 # undef SELF 00089 # undef PointSet 00090 # undef BoundingBox 00091 # endif // shape_algebraic_curve_fcts_hpp
