// Copyright (c) 2005 Rijksuniversiteit Groningen (Netherlands) // All rights reserved. // // This file is part of CGAL (www.cgal.org). // You can redistribute it and/or modify it under the terms of the GNU // General Public License as published by the Free Software Foundation, // either version 3 of the License, or (at your option) any later version. // // Licensees holding a valid commercial license may use this file in // accordance with the commercial license agreement provided with the software. // // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // // $URL$ // $Id$ // // // Author(s) : Nico Kruithof #ifndef CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_H #define CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_H #include namespace CGAL { template < class SkinSurfaceQuadraticSurfaceTraits_3 > class Skin_surface_quadratic_surface_3 { public: typedef SkinSurfaceQuadraticSurfaceTraits_3 K; typedef typename K::FT FT; typedef typename K::Point_3 Point; typedef typename K::Vector_3 Vector; typedef typename K::Segment_3 Segment; typedef CGAL::Weighted_point Weighted_point; Skin_surface_quadratic_surface_3() : dim(-1), p(0,0,0), c(0) { for (int i=0; i<6; i++) Q[i] = 0; } Skin_surface_quadratic_surface_3(FT Qinput[], Point p, FT c, int d) : dim(10+d), p(p), c(c) { for (int i=0; i<6; i++) Q[i] = Qinput[i]; } int dim; Skin_surface_quadratic_surface_3(Weighted_point wp0, FT s) : dim(0), p(wp0.point()), c(-s*(1-s)*wp0.weight()) { CGAL_PROFILER(std::string("Constructor : ") + std::string(CGAL_PRETTY_FUNCTION)); Q[1] = Q[3] = Q[4] = 0; Q[0] = Q[2] = Q[5] = (1-s); } Skin_surface_quadratic_surface_3(Weighted_point wp0, Weighted_point wp1, FT s) : dim(1) { CGAL_PROFILER(std::string("Constructor : ") + std::string(CGAL_PRETTY_FUNCTION)); Regular_triangulation_euclidean_traits_3 reg_traits; p = reg_traits.construct_weighted_circumcenter_3_object()(wp0,wp1); c = s*(1-s)*reg_traits.compute_squared_radius_smallest_orthogonal_sphere_3_object()(wp0,wp1); Vector t = wp0-wp1; FT den = t*t; Q[0] = (- t.x()*t.x()/den + (1-s)); Q[1] = (-2*t.y()*t.x()/den); Q[2] = (- t.y()*t.y()/den + (1-s)); Q[3] = (-2*t.z()*t.x()/den); Q[4] = (-2*t.z()*t.y()/den); Q[5] = (- t.z()*t.z()/den + (1-s)); } Skin_surface_quadratic_surface_3(Weighted_point wp0, Weighted_point wp1, Weighted_point wp2, FT s) : dim(2) { CGAL_PROFILER(std::string("Constructor : ") + std::string(CGAL_PRETTY_FUNCTION)); Regular_triangulation_euclidean_traits_3 reg_traits; p = reg_traits.construct_weighted_circumcenter_3_object()(wp0,wp1,wp2); c = s*(1-s)*reg_traits.compute_squared_radius_smallest_orthogonal_sphere_3_object()(wp0,wp1,wp2); Vector t = K().construct_orthogonal_vector_3_object()(wp0,wp1,wp2); FT den = t*t; Q[0] = -(- t.x()*t.x()/den + s); Q[1] = -(-2*t.y()*t.x()/den); Q[2] = -(- t.y()*t.y()/den + s); Q[3] = -(-2*t.z()*t.x()/den); Q[4] = -(-2*t.z()*t.y()/den); Q[5] = -(- t.z()*t.z()/den + s); } Skin_surface_quadratic_surface_3(Weighted_point wp0, Weighted_point wp1, Weighted_point wp2, Weighted_point wp3, FT s) : dim(3) { CGAL_PROFILER(std::string("Constructor : ") + std::string(CGAL_PRETTY_FUNCTION)); Regular_triangulation_euclidean_traits_3 reg_traits; p = reg_traits.construct_weighted_circumcenter_3_object()(wp0,wp1,wp2,wp3); c = s*(1-s)*reg_traits.compute_squared_radius_smallest_orthogonal_sphere_3_object()(wp0,wp1,wp2,wp3); Q[1] = Q[3] = Q[4] = 0; Q[0] = Q[2] = Q[5] = -s; } template FT value(Input_point const &x) const { typedef Cartesian_converter Converter; FT vx = Converter()(x.x()) - p.x(); FT vy = Converter()(x.y()) - p.y(); FT vz = Converter()(x.z()) - p.z(); return vx*(Q[0]*vx) + vy*(Q[1]*vx+Q[2]*vy) + vz*(Q[3]*vx+Q[4]*vy+Q[5]*vz) + c; } template Sign sign(Input_point const &x) const { return CGAL_NTS sign(value(x)); } template typename Input_point::R::Vector_3 gradient(Input_point const &x) { // NGHK: We use the kernel trick again: DOCUMENT!!!! typedef Cartesian_converter Converter; typedef Cartesian_converter Inv_converter; Converter conv; Inv_converter inv_conv; Point xp = conv(x); return inv_conv(compute_gradient(xp)); } /// 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; 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 midpoint(pp0,pp1); } private: //template Vector compute_gradient(Point const &x) { FT vx = x.x() - p.x(); FT vy = x.y() - p.y(); FT vz = x.z() - p.z(); return Vector(2*Q[0]*vx + Q[1]*vy + Q[3]*vz, Q[1]*vx + 2*Q[2]*vy + Q[4]*vz, Q[3]*vx + Q[4]*vy + 2*Q[5]*vz); } FT Q[6]; Point p; FT c; }; } //namespace CGAL #endif // CGAL_SKIN_SURFACE_QUADRATIC_SURFACE_H