00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 # ifndef shape_rational_surface_hpp
00013 # define shape_rational_surface_hppx
00014
00015 # include <realroot/ring_sparse.hpp>
00016 # include <realroot/ring_bernstein_tensor.hpp>
00017 # include <shape/algebraic_set.hpp>
00018 # include <shape/parametric_surface.hpp>
00019 # include <math.h>
00020 # define TMPL_DEF template<class C, int N=3, class V= default_env>
00021 # define TMPL template<class C, int N, class V>
00022 # define ParametricSurface parametric_surface<C,V>
00023 # define SELF rational_surface<C,N,V>
00024
00025 namespace mmx { namespace shape {
00026
00027 TMPL_DEF
00028 class rational_surface : public ParametricSurface {
00029
00030 public:
00031 typedef C Scalar;
00032 typedef typename use<point_def,C,V>::Point Point;
00033 typedef typename point_set<C>::PointIterator PointIterator;
00034
00035 typedef typename algebraic_set<C,V>::Polynomial Polynomial;
00036
00037 rational_surface(void): m_w(1) {};
00038
00039
00040 rational_surface(const Polynomial& X, const Polynomial& Y, const Polynomial& Z, const Polynomial& W=1):
00041 ParametricSurface() {
00042 set_range(0,1,0,1);
00043 m_w=W;
00044 m_p[0]=X; m_p[1]=Y; m_p[2]=Z;
00045 }
00046
00047 Point* eval(Scalar u, Scalar v) const ;
00048 void eval(Point&, Scalar u, Scalar v) const ;
00049 void eval(Scalar&, Scalar&, Scalar&, Scalar u, Scalar v) const ;
00050
00051 void normal(Point&, Scalar u, Scalar v) const ;
00052 void normal(Scalar&, Scalar&, Scalar&, Scalar u, Scalar v) const ;
00053
00054 Scalar umin() const;
00055 Scalar umax() const;
00056 Scalar vmin() const;
00057 Scalar vmax() const;
00058
00059 Point* operator() (const Scalar& u, const Scalar& v) const ;
00060
00061 void set_range(Scalar umin, Scalar umax, Scalar vmin, Scalar vmax) ;
00062 void add_to_domain(const Point& p) ;
00063 Seq<Point> domain(void) const;
00064 Point domain(int i) const {return m_domain[i];}
00065
00066 Polynomial denominator() const { return m_w; }
00067 Polynomial numerator(int i) const { return m_p[i]; }
00068
00069
00070 void set_denominator(const Polynomial& d) { m_w=d; }
00071 void set_numerator (const Polynomial& p, int i) { m_p[i]=p; }
00072
00073 private:
00074 Seq<Point> m_domain;
00075 Polynomial m_w, m_p[N];
00076 } ;
00077
00078
00079 TMPL typename SELF::Point*
00080 SELF::eval(Scalar u, Scalar v) const {
00081 Scalar w = m_w(u,v);
00082 Point* p = new Point(m_p[0](u,v)/w, m_p[1](u,v)/w, m_p[2](u,v)/w);
00083 return p ;
00084 }
00085
00086 TMPL void
00087 SELF::eval(Point& p, Scalar u, Scalar v) const {
00088 Scalar w=m_w(u,v);
00089 if (w != 0) {
00090 p.setx(m_p[0](u,v)/w);
00091 p.sety(m_p[1](u,v)/w);
00092 p.setz(m_p[2](u,v)/w);
00093 } else {
00094 std::cerr << "Division by 0 in rational_surface::eval(Point& p, Scalar u, Scalar v) const" << std::endl;
00095 }
00096 }
00097
00098 TMPL void
00099 SELF::normal(Point& p, Scalar u, Scalar v) const {
00100 normal(p[0],p[1],p[2],u,v);
00101
00102 }
00103 TMPL void
00104 SELF::normal(Scalar& x, Scalar& y, Scalar &z, Scalar u, Scalar v) const {
00105 Polynomial dw0 = diff(m_w,0) , dw1 = diff(m_w,1),
00106 dx0 = diff(m_p[0],0),dx1 = diff(m_p[0],1),
00107 dy0 = diff(m_p[1],0), dy1 = diff(m_p[1],1),
00108 dz0 = diff(m_p[2],0), dz1 = diff(m_p[2],1);
00109
00110 Scalar w=m_w(u,v);
00111 if (w != 0) {
00112
00113 Scalar x0 = (dx0(u,v)* w - m_p[0](u,v)*dw0(u,v))/(w*w) ;
00114 Scalar x1 = (dx1(u,v)* w - m_p[0](u,v)*dw1(u,v))/(w*w);
00115 Scalar y0 = (dy0(u,v)* w - m_p[1](u,v)*dw0(u,v))/(w*w) ;
00116 Scalar y1 = (dy1(u,v)* w - m_p[1](u,v)*dw1(u,v))/(w*w);
00117 Scalar z0 = (dz0(u,v)* w - m_p[2](u,v)*dw0(u,v))/(w*w) ;
00118 Scalar z1 = (dx1(u,v)* w - m_p[2](u,v)*dw1(u,v))/(w*w);
00119
00120 x= (y0 * z1) - (z0 * y1);
00121 y= -(x0 * z1) + (z0 * x1);
00122 z= (x0 * y1) - (y0 * x1);
00123 } else {
00124 std::cerr << "Division by 0 in rational_surface::eval(Scalar& x, Scalar& y, Scalar &z, Scalar u, Scalar v const" << std::endl;
00125 }
00126
00127
00128
00129
00130 Scalar a = sqrtf(x*x+y*y+z*z);
00131 if (a != 0) {
00132 x /= a;
00133 y /= a;
00134 z /= a;
00135
00136
00137 } else{
00138 std::cerr << "Division by 0 in rational_surface::eval(Scalar& x, Scalar& y, Scalar &z, Scalar u, Scalar v const" << std::endl;
00139 }
00140
00141 }
00142 TMPL void
00143 SELF::eval(Scalar& x, Scalar& y, Scalar &z, Scalar u, Scalar v) const {
00144 Scalar w=m_w(u,v);
00145 if (w != 0) {
00146 x=m_p[0](u,v)/w;
00147 y=m_p[1](u,v)/w;
00148 z=m_p[2](u,v)/w;
00149 } else {
00150 std::cerr << "Division by 0 in rational_surface::eval(Scalar& x, Scalar& y, Scalar &z, Scalar u, Scalar v const" << std::endl;
00151 }
00152 }
00153
00154 TMPL C
00155 SELF::umin(void) const {
00156 Scalar r=m_domain[0].x();
00157 for(unsigned i=1; i<m_domain.size();i++)
00158 r=std::min(r,m_domain[i].x());
00159 return r;
00160 }
00161
00162 TMPL C
00163 SELF::umax(void) const {
00164 Scalar r=m_domain[0].x();
00165 for(unsigned i=1; i<m_domain.size();i++)
00166 r=std::max(r,m_domain[i].x());
00167 return r;
00168 }
00169
00170 TMPL C
00171 SELF::vmin(void) const {
00172 Scalar r=m_domain[0].y();
00173 for(unsigned i=1; i<m_domain.size();i++)
00174 r=std::min(r,m_domain[i].y());
00175 return r;
00176 }
00177
00178 TMPL C
00179 SELF::vmax(void) const {
00180 Scalar r=m_domain[0].y();
00181 for(unsigned i=1; i<m_domain.size();i++)
00182 r=std::max(r,m_domain[i].y());
00183 return r;
00184 }
00185
00186 TMPL typename SELF::Point*
00187 SELF::operator() (const Scalar& u, const Scalar& v) const {
00188 return this->eval(u,v);
00189 }
00190
00191 TMPL void
00192 SELF::set_range(Scalar umin, Scalar umax, Scalar vmin, Scalar vmax) {
00193 m_domain.resize(0);
00194 m_domain << Point(umin,vmin) << Point(umax,vmin)
00195 << Point(umin,vmax) << Point(umax,vmin);
00196 }
00197
00198 TMPL void
00199 SELF::add_to_domain(const Point& p) {
00200 this->m_domain <<p;
00201 }
00202
00203
00204 TMPL Seq<typename SELF::Point>
00205 SELF::domain() const {
00206 return m_domain;
00207 }
00208
00209 template<class MESH, class C, int N, class V>
00210 void as_graphic(MESH* mesh, const SELF& surf, unsigned m = 50, unsigned n = 50) {
00211
00212 typedef typename SELF::Scalar Scalar;
00213 typedef typename SELF::Point Point;
00214 typedef typename SELF::Polynomial Polynomial;
00215 typedef typename MESH::Point MPoint;
00216 std::cout<<"domain: "<<surf.domain()<<std::endl;
00217 std::cout<<"numer: "<<surf.numerator(0)<<" "<<surf.numerator(1)<<" "<<surf.numerator(2)<<std::endl;
00218 std::cout<<"denom: "<<surf.denominator()<<std::endl;
00219 if (surf.domain().size() != 0 && surf.domain().size() != 3 && surf.domain().size() != 4) return;
00220
00221 int c=0;
00222
00223
00224 Point A1, dA1,
00225 U, B0, B1, dB0, P(0,0,0);
00226 if (surf.domain().size() == 0) {
00227 if (m == 0 || n == 0) {
00228 std::cerr << "Error in as_graphic: division by 0";
00229 return;
00230 }
00231 double umin = -1.0+10e-6, umax = -umin;
00232
00233 double ds = (umax - umin) / (double)m;
00234 double dt = (umax - umin) / (double)n;
00235
00236 for (double s = umin; s <= umax; s += ds) {
00237 for (double t = umin; t <= umax; t += dt) {
00238 surf.eval(P,s/(1-s*s), t/(1-t*t));
00239
00240 mesh->push_back_vertex(P.x(),P.y(),P.z());
00241 if (s > umin) {
00242 if (t > umin) {
00243 mesh->push_back_face(c-n-1, c-1, c);
00244 }
00245 if (t < umax) {
00246 mesh->push_back_face(c-n, c-n-1, c);
00247 }
00248 }
00249 c++;
00250 }
00251 }
00252 } else {
00253 Point A0 = surf.domain(0), dA0 = (surf.domain(1)-A0)/m;
00254
00255 if (surf.domain().size()==3) {
00256 A1 = surf.domain(2);
00257 dA1 = dA0;
00258 B0=A0; B1=A1;
00259 for (int i=0;i<=(int)m;i++) {
00260 U=B0;
00261 Point dB0= (B1-B0)/m;
00262 for (int j=0;j<=(int)m-i;j++) {
00263 surf.eval(P,U.x(),U.y());
00264
00265
00266 mesh->push_back_vertex(P.x(),P.y(),P.z());
00267 if (i> 0) {
00268 if (j>0) mesh->push_back_face(c-m+i-2, c-1, c);
00269
00270 mesh->push_back_face(c-m+i-1, c-m+i-2, c);
00271 }
00272 c++;
00273 U=U+dB0;
00274 }
00275 B0 = B0 + dA0;
00276 B1 = B1 + dA0;
00277 }
00278 } else if (surf.domain().size()==4) {
00279 A1 = surf.domain(3);
00280 dA1 = (surf.domain(2)-A1)/m;
00281 B0=A0; B1=A1;
00282 for (unsigned i=0;i<=m;i++) {
00283 U=B0;
00284 Point dB0= (B1-B0)/n;
00285 for (unsigned j=0;j<=n;j++) {
00286 surf.eval(P,U.x(),U.y());
00287
00288 mesh->push_back_vertex(P.x(),P.y(),P.z());
00289 if (i> 0) {
00290 if (j>0) mesh->push_back_face(c-n-1, c-1, c);
00291 if (j<n) mesh->push_back_face(c-n, c-n-1, c);
00292 }
00293 c++;
00294 U=U+dB0;
00295 }
00296 B0 = B0 + dA0;
00297 B1 = B1 + dA1;
00298 }
00299 }
00300 }
00301 std::cout<< mesh->normal_count()<< std::endl;
00302 }
00303
00304 }
00305
00306 }
00307
00308 # undef TMPL
00309 # undef TMPL1
00310 # undef TMPL_DEF
00311 # undef Scalar
00312 # undef Point
00313 # undef SELF
00314 # undef ParametricSurface
00315 # undef AXEL
00316 # endif // rational_surface_hpp
