From f0bf17c1a0257541e684afb2e55622c5d65ebdf5 Mon Sep 17 00:00:00 2001 From: Nico Kruithof Date: Wed, 3 May 2006 13:28:33 +0000 Subject: [PATCH] More changes to the interface --- .../Skin_surface_3/skin_surface_simple.C | 2 - .../CGAL/Skin_surface_quadratic_surface_3.h | 432 +++++++++--------- .../Triangulated_mixed_complex_observer_3.h | 44 +- .../test/Skin_surface_3/surface_test.C | 5 +- 4 files changed, 240 insertions(+), 243 deletions(-) diff --git a/Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.C b/Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.C index 144f96472a6..5a7e43adc73 100644 --- a/Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.C +++ b/Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.C @@ -1,5 +1,3 @@ -// NGHK: remove later -#define CGAL_PROFILE // examples/Skin_surface_3/skin_surface_simple.C #include #include diff --git a/Skin_surface_3/include/CGAL/Skin_surface_quadratic_surface_3.h b/Skin_surface_3/include/CGAL/Skin_surface_quadratic_surface_3.h index a749d042188..490499be5cb 100644 --- a/Skin_surface_3/include/CGAL/Skin_surface_quadratic_surface_3.h +++ b/Skin_surface_3/include/CGAL/Skin_surface_quadratic_surface_3.h @@ -24,13 +24,6 @@ CGAL_BEGIN_NAMESPACE -///////////////////////////////////////////////// -// // -// Virtual base class of a quadratic surface // -// This implementation is specific for // -// skin surfaces // -// // -///////////////////////////////////////////////// template < class Kernel > class Skin_surface_quadratic_surface_3 { public: @@ -41,9 +34,32 @@ public: typedef typename K::RT RT; typedef Weighted_point Weighted_point; - Skin_surface_quadratic_surface_3() { + Skin_surface_quadratic_surface_3(RT Qinput[], Point p, RT c) + : p(p), c(c) + { + for (int i=0; i<6; i++) Q[i] = Qinput[i]; } - virtual ~Skin_surface_quadratic_surface_3() {}; + + template + RT value(Input_point const &x) const { + typedef Cartesian_converter Converter; + + Vector v = Converter()(x) - p; + + return + v.x()*(Q[0]*v.x()) + + v.y()*(Q[1]*v.x()+Q[2]*v.y()) + + v.z()*(Q[3]*v.x()+Q[4]*v.y()+Q[5]*v.z()) + + c; + } + + + Vector gradient(Point const &x) { + std::cout << "NGHK: NOT YET IMPLEMENTED" << std::endl; + // NGHK: TODO: + return (x-p); + } + /// Construct the intersection point with the segment (p0,p1) template Input_point to_surface(const Input_point &p0, const Input_point &p1) { @@ -69,247 +85,215 @@ public: } return mid; }; -// virtual Point to_surface(Point const &p, Vector const &v) = 0; - inline Point to_surface(const Segment &s) { - return to_surface(s.source(), s.target()); - } - -// // Gradient descent -// virtual Point to_surface(Point const &p0) = 0; - - // compute the function value of p - template - RT value(const Input_point &p) const { - typedef Cartesian_converter Converter; - return _value(Converter()(p)); - } - virtual RT _value(const Point &p) const = 0; - // Compute the gradient in p - virtual Vector gradient(const Point &p) = 0; - - // return the dimension of the delaunay simplex: - virtual int dimension() const = 0; - - // return a continuous density function on the skin surface: - virtual RT sq_density(const Point &p) = 0; - -}; - -///////////////////////////////////////////////// -// // -// Sphere // -// // -// orient*|x-p|^2/s - P == 0 // -// orient*|x-p|^2/s - P == 0 // -// // -// p: center of the sphere // -// P: squared radius of the sphere // -// orient: orientation of the sphere // -// // -///////////////////////////////////////////////// -template < class Kernel > -class Skin_surface_sphere_3 : public Skin_surface_quadratic_surface_3 { -public: - typedef Skin_surface_quadratic_surface_3 Parent; - typedef typename Parent::Point Point; - typedef typename Parent::Vector Vector; - typedef typename Parent::RT RT; - typedef typename Parent::Weighted_point Weighted_point; - - Skin_surface_sphere_3(Weighted_point wp, RT s, int orient) - : Skin_surface_quadratic_surface_3(), - wp(wp), s(s), orient(orient) { - assert((orient == -1) || (orient == 1)); - } - - RT _value(Point const &x) const { - return orient * (squared_distance(x, wp)/s - wp.weight()); - } - - Point to_surface(const Point &p, const Vector &v) { - Vector pc = p-wp; - RT sq_d = v*v; - RT top = -pc*v/sq_d; - - RT extr_val = squared_distance(p+top*v, wp) - s*wp.weight(); - - if (extr_val > 0) { - // std::cerr << "im. intersection[" << dimension() << "] " - // << extr_val << "\n"; - return p; - } - RT d = sqrt(-extr_val/sq_d); - - // t should be in [0,1] - RT t, t1; - t = top + d; t1 = top - d; - if (t1*t1 < t*t) t = t1; - - CGAL_assertion(abs(value(p + t*v)) < 0.001); - return p + t*v; - } - Point to_surface(Point const &p) { - Vector pc = p-wp; - return wp + sqrt(s*wp.weight()/(pc*pc))*pc; - } - - - Vector gradient(Point const &p) { - return orient*(p-wp); - } - - int dimension() const { - if (orient == 1) return 0; else return 3; - } - - RT sq_density(Point const &p) { - assert(wp.weight() > 0); - return s*wp.weight(); - } private: - Weighted_point wp; - RT s; - int orient; + RT Q[6]; + Point p; + RT c; }; -/* ************************************************* - * Hyperboloid - * - * Point p. - * Let y = (p-wp)*t/|t| - * x = |p-wp|^2 - y^2 - * - * orient*((1-s) x^2 - s y^2 + s(1-s) W = 0 - * orient*((1-s)|p-wp|^2 - y^2 + s(1-s)W = 0 - * - ***************************************************/ +// template < class Kernel > +// class Skin_surface_quadratic_surface_3 { +// public: +// typedef Kernel K; +// typedef typename K::Point_3 Point; +// typedef typename K::Vector_3 Vector; +// typedef typename K::Segment_3 Segment; +// typedef typename K::RT RT; +// typedef Weighted_point Weighted_point; -template < class Kernel > -class Skin_surface_hyperboloid_3 : public Skin_surface_quadratic_surface_3 { -public: - typedef Skin_surface_quadratic_surface_3 Parent; - typedef typename Parent::Point Point; - typedef typename Parent::Vector Vector; - typedef typename Parent::RT RT; - typedef typename Parent::Weighted_point Weighted_point; +// Skin_surface_quadratic_surface_3() { +// } +// virtual ~Skin_surface_quadratic_surface_3() {}; +// /// Construct the intersection point with the segment (p0,p1) +// template +// Input_point to_surface(const Input_point &p0, const Input_point &p1) { +// // NGHK: We use the kernel trick again: DOCUMENT!!!! +// typedef Cartesian_converter Converter; +// Converter conv; - Skin_surface_hyperboloid_3(Weighted_point wp, Vector t, RT s, int orient) - : Skin_surface_quadratic_surface_3(), wp(wp), t(t), s(s), orient(orient) { - assert((orient == -1) || (orient == 1)); - sq_t = t*t; - } +// Input_point pp0 = p0, pp1 = p1, mid; +// double sq_d = to_double(squared_distance(pp0,pp1)); + +// if (value(conv(pp1)) < value(conv(pp0))) { +// std::swap(pp0, pp1); +// } + +// while (sq_d > 1.e-10) { +// mid = midpoint(pp0,pp1); +// if (value(conv(mid)) > 0) { +// pp1 = mid; +// } else { +// pp0 = mid; +// } +// sq_d /= 4; +// } +// return mid; +// }; +// // virtual Point to_surface(Point const &p, Vector const &v) = 0; +// inline Point to_surface(const Segment &s) { +// return to_surface(s.source(), s.target()); +// } + +// // // Gradient descent +// // virtual Point to_surface(Point const &p0) = 0; + +// // compute the function value of p +// template +// RT value(const Input_point &p) const { +// typedef Cartesian_converter Converter; +// return _value(Converter()(p)); +// } +// virtual RT _value(const Point &p) const = 0; +// // Compute the gradient in p +// virtual Vector gradient(const Point &p) = 0; + +// // // NGHK: REMOVE: +// virtual int dimension() const = 0; + +// // // return a continuous density function on the skin surface: +// // virtual RT sq_density(const Point &p) = 0; + +// }; + + +// /* ************************************************* +// * Hyperboloid +// * +// * Point p. +// * Let y = (p-wp)*t/|t| +// * x = |p-wp|^2 - y^2 +// * +// * orient*((1-s) x^2 - s y^2 + s(1-s) W = 0 +// * orient*((1-s)|p-wp|^2 - y^2 + s(1-s)W = 0 +// * +// ***************************************************/ + +// template < class Kernel > +// class Skin_surface_hyperboloid_3 : public Skin_surface_quadratic_surface_3 { +// public: +// typedef Skin_surface_quadratic_surface_3 Parent; +// typedef typename Parent::Point Point; +// typedef typename Parent::Vector Vector; +// typedef typename Parent::RT RT; +// typedef typename Parent::Weighted_point Weighted_point; + +// Skin_surface_hyperboloid_3(Weighted_point wp, Vector t, RT s, int orient) +// : Skin_surface_quadratic_surface_3(), wp(wp), t(t), s(s), orient(orient) { +// assert((orient == -1) || (orient == 1)); +// sq_t = t*t; +// } - RT _value(Point const &x) const { - Vector dir = x-wp; - RT tmp = dir*t; - tmp = tmp*tmp/sq_t; +// RT _value(Point const &x) const { +// Vector dir = x-wp; +// RT tmp = dir*t; +// tmp = tmp*tmp/sq_t; - return orient * (dir*dir/s - tmp/(s*(1-s))) + wp.weight(); - } -// Point to_surface(Point const &p0, Point const &p1) { -// assert(value(p0) * value(p1) <= 0); +// return orient * (dir*dir/s - tmp/(s*(1-s))) + wp.weight(); +// } +// // Point to_surface(Point const &p0, Point const &p1) { +// // assert(value(p0) * value(p1) <= 0); -// Vector p0c = p0-wp; -// Vector p1p0 = p1-p0; -// RT sq_d = p1p0*p1p0; +// // Vector p0c = p0-wp; +// // Vector p1p0 = p1-p0; +// // RT sq_d = p1p0*p1p0; -// RT p0_sym = p0c*t; +// // RT p0_sym = p0c*t; +// // RT p1_sym = (p1-wp)*t; +// // RT d_sym = p0_sym-p1_sym; + +// // RT den = ((1-s)*sq_d - d_sym*d_sym/sq_t); +// // RT top = -((1-s)*(p0c*p1p0) + p0_sym*d_sym/sq_t) / den; +// // RT extr_val = +// // orient*((1-s)*p0c*p0c-p0_sym*p0_sym/sq_t-s*(1-s)*wp.weight()); + +// // { +// // Point extr = p0+top*p1p0; +// // Vector dir = extr-wp; +// // RT tmp = dir*t; tmp *= tmp/sq_t; +// // extr_val = orient*(-tmp + (1-s)*dir*dir + s*(1-s)*wp.weight()); +// // } +// // RT d = sqrt(-orient*extr_val/den); + +// // RT t, t1; +// // t = top + d; t1 = top - d; +// // if ((2*t1-1)*(2*t1-1) < (2*t-1)*(2*t-1)) t = t1; + +// // if (t < 0) { +// // std::cerr << "Hl[" << dimension() <<"] " << t << "\n"; +// // return p0; +// // } else if (t > 1) { +// // std::cerr << "Hh[" << dimension() <<"] " << t << "\n"; +// // return p1; +// // } else { +// // assert (std::abs(value(p0 + t*p1p0)) < 0.001); +// // return p0 + t*p1p0; +// // } +// // } +// Point to_surface(Point const &p, Vector const &v){ +// Vector pc = p-wp; +// Point p1 = p + v; +// RT sq_d = v*v; + +// RT p_sym = pc*t; // RT p1_sym = (p1-wp)*t; -// RT d_sym = p0_sym-p1_sym; +// RT d_sym = p_sym-p1_sym; // RT den = ((1-s)*sq_d - d_sym*d_sym/sq_t); -// RT top = -((1-s)*(p0c*p1p0) + p0_sym*d_sym/sq_t) / den; +// RT top = -((1-s)*(pc*v) + p_sym*d_sym/sq_t) / den; // RT extr_val = -// orient*((1-s)*p0c*p0c-p0_sym*p0_sym/sq_t-s*(1-s)*wp.weight()); +// orient*((1-s)*pc*pc-p_sym*p_sym/sq_t-s*(1-s)*wp.weight()); // { -// Point extr = p0+top*p1p0; +// Point extr = p+top*v; // Vector dir = extr-wp; // RT tmp = dir*t; tmp *= tmp/sq_t; // extr_val = orient*(-tmp + (1-s)*dir*dir + s*(1-s)*wp.weight()); // } -// RT d = sqrt(-orient*extr_val/den); +// RT d = -orient*extr_val/den; +// if (d<0) return p; +// d = sqrt(d); // RT t, t1; // t = top + d; t1 = top - d; -// if ((2*t1-1)*(2*t1-1) < (2*t-1)*(2*t-1)) t = t1; +// if (t1*t1 < t*t) t = t1; -// if (t < 0) { -// std::cerr << "Hl[" << dimension() <<"] " << t << "\n"; -// return p0; -// } else if (t > 1) { -// std::cerr << "Hh[" << dimension() <<"] " << t << "\n"; -// return p1; -// } else { -// assert (std::abs(value(p0 + t*p1p0)) < 0.001); -// return p0 + t*p1p0; -// } +// assert (std::abs(value(p + t*v)) < 0.001); +// return p + t*v; // } - Point to_surface(Point const &p, Vector const &v){ - Vector pc = p-wp; - Point p1 = p + v; - RT sq_d = v*v; +// Point to_surface(Point const &p0) { +// return to_surface(p0, gradient(p0)); +// } +// Vector gradient(Point const &p) { +// // -s x + (1-s) y +// Vector v = p - wp; +// Vector vt = (v*t)/(t*t)*t; +// return orient*((1-s)*v - vt); +// } + +// int dimension() const { +// if (orient == 1) return 1; else return 2; +// } + +// RT sq_density(Point const &p) { +// Vector v = p - wp; +// RT vt = v*t; +// RT scale = 1-s; +// if (s < RT(.5)) { +// scale = (scale-s)/(scale*scale); +// } else { +// scale = (s-scale)/(s*s); +// } +// assert(scale >= 0); - RT p_sym = pc*t; - RT p1_sym = (p1-wp)*t; - RT d_sym = p_sym-p1_sym; - - RT den = ((1-s)*sq_d - d_sym*d_sym/sq_t); - RT top = -((1-s)*(pc*v) + p_sym*d_sym/sq_t) / den; - RT extr_val = - orient*((1-s)*pc*pc-p_sym*p_sym/sq_t-s*(1-s)*wp.weight()); - - { - Point extr = p+top*v; - Vector dir = extr-wp; - RT tmp = dir*t; tmp *= tmp/sq_t; - extr_val = orient*(-tmp + (1-s)*dir*dir + s*(1-s)*wp.weight()); - } - RT d = -orient*extr_val/den; - if (d<0) return p; - d = sqrt(d); - - RT t, t1; - t = top + d; t1 = top - d; - if (t1*t1 < t*t) t = t1; - - assert (std::abs(value(p + t*v)) < 0.001); - return p + t*v; - } - Point to_surface(Point const &p0) { - return to_surface(p0, gradient(p0)); - } - Vector gradient(Point const &p) { - // -s x + (1-s) y - Vector v = p - wp; - Vector vt = (v*t)/(t*t)*t; - return orient*((1-s)*v - vt); - } - - int dimension() const { - if (orient == 1) return 1; else return 2; - } - - RT sq_density(Point const &p) { - Vector v = p - wp; - RT vt = v*t; - RT scale = 1-s; - if (s < RT(.5)) { - scale = (scale-s)/(scale*scale); - } else { - scale = (s-scale)/(s*s); - } - assert(scale >= 0); - - return squared_distance(wp, p) - scale*vt*vt/(t*t); - } -private: - Weighted_point wp; - Vector t; - RT s, sq_t; - int orient; +// return squared_distance(wp, p) - scale*vt*vt/(t*t); +// } +// private: +// Weighted_point wp; +// Vector t; +// RT s, sq_t; +// int orient; -}; +// }; CGAL_END_NAMESPACE diff --git a/Skin_surface_3/include/CGAL/Triangulated_mixed_complex_observer_3.h b/Skin_surface_3/include/CGAL/Triangulated_mixed_complex_observer_3.h index 98dd09836d6..c14799b36c9 100644 --- a/Skin_surface_3/include/CGAL/Triangulated_mixed_complex_observer_3.h +++ b/Skin_surface_3/include/CGAL/Triangulated_mixed_complex_observer_3.h @@ -67,9 +67,9 @@ public: // typedef typename Triangulated_mixed_complex_traits::RT TMC_RT; // typedef Weighted_point TMC_Weighted_point; -// typedef Skin_surface_quadratic_surface_3 QuadrSurface; - typedef Skin_surface_sphere_3 Sphere_surface; - typedef Skin_surface_hyperboloid_3 Hyperboloid_surface; + typedef Skin_surface_quadratic_surface_3 Quadratic_surface; +// typedef Skin_surface_sphere_3 Sphere_surface; +// typedef Skin_surface_hyperboloid_3 Hyperboloid_surface; typedef Weighted_converter_3< Cartesian_converter < typename Regular_traits::Bare_point::R, @@ -111,8 +111,9 @@ public: { vh = s; create_sphere(r2s_converter(vh->point().point()), - r2s_converter(vh->point().weight()), - shrink); + -r2s_converter(vh->point().weight()), + shrink, + 1); break; } case 1: @@ -170,7 +171,8 @@ public: r2s_converter(ch->vertex(1)->point()), r2s_converter(ch->vertex(2)->point()), r2s_converter(ch->vertex(3)->point())), - shrink-1); + 1-shrink, + -1); } } } @@ -183,22 +185,34 @@ public: void create_sphere(const Surface_point &c, const Surface_RT &w, - const Surface_RT &s) { - if (s >= 0) { - surf = new Sphere_surface(Surface_weighted_point(c,w), s, 1); - } else { - surf = new Sphere_surface(Surface_weighted_point(c,w), -s, -1); - } + const Surface_RT &s, + const int orient) { + Q[1] = Q[3] = Q[4] = 0; + Q[0] = Q[2] = Q[5] = orient*(1-s); + + surf = new Quadratic_surface(Q, c, s*(1-s)*w); } void create_hyperboloid(const Surface_point &c, const Surface_RT &w, - const Surface_vector &v, + const Surface_vector &t, const Surface_RT &s, - int orient) { - surf = new Hyperboloid_surface(Surface_weighted_point(c,w), v, s, orient); + const int orient) { + Surface_RT den = t*t; + + Q[0] = orient*(- t.x()*t.x()/den + (1-s)); + + Q[1] = orient*(-2*t.y()*t.x()/den); + Q[2] = orient*(- t.y()*t.y()/den + (1-s)); + + Q[3] = orient*(-2*t.z()*t.x()/den); + Q[4] = orient*(-2*t.z()*t.y()/den); + Q[5] = orient*(- t.z()*t.z()/den + (1-s)); + + surf = new Quadratic_surface(Q, c, s*(1-s)*w); } + Surface_RT Q[6]; R2S_converter r2s_converter; }; diff --git a/Skin_surface_3/test/Skin_surface_3/surface_test.C b/Skin_surface_3/test/Skin_surface_3/surface_test.C index 34d6d4bfaac..24728c4eea9 100644 --- a/Skin_surface_3/test/Skin_surface_3/surface_test.C +++ b/Skin_surface_3/test/Skin_surface_3/surface_test.C @@ -53,7 +53,7 @@ typedef Triangulated_mixed_complex::Finite_vertices_iterator #include int main(int argc, char *argv[]) { - double shrink = .85; + double shrink = .5; std::list l; while (argc>1) { @@ -71,6 +71,7 @@ int main(int argc, char *argv[]) { for (Tmc_Finite_vertices_iterator vit = tmc.finite_vertices_begin(); vit != tmc.finite_vertices_end(); vit++) { + std::cout << std::endl; if (tmc.is_infinite(vit->cell())) { std::cerr << "ERROR: is_infinite (main)" << std::endl; } @@ -82,7 +83,7 @@ int main(int argc, char *argv[]) { cell != cells.end(); cell++) { if (!tmc.is_infinite(*cell)) { Quadratic_surface::RT val2 = (*cell)->surf->value(vit->point()); - CGAL_assertion(val == val2); + CGAL_assertion(val == val2); } } }