diff --git a/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_sphere_traits_2.h b/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_sphere_traits_2.h index 77cecae2da1..ba7b8f142f1 100644 --- a/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_sphere_traits_2.h +++ b/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_sphere_traits_2.h @@ -17,6 +17,7 @@ #include #include +#include #include #include @@ -24,76 +25,129 @@ namespace CGAL { namespace internal { -template +template class Orientation_on_sphere_2 { public: - typedef typename Traits::Point_3 Point_3; + typedef typename LK::Point_3 Point_3; + typedef typename LK::Point_3 Point_on_sphere_2; typedef Comparison_result result_type; - Orientation_on_sphere_2(const Point_3& center, const Traits& traits) - : _center(center), _traits(traits) + Orientation_on_sphere_2(const Point_3& center, const LK& lk) + : _center(center), _lk(lk) { } - Comparison_result operator()(const Point_3& p, const Point_3& q, const Point_3& r) const - { return _traits.orientation_3_object()(_center, p, q, r); } + Comparison_result operator()(const Point_on_sphere_2& p, const Point_on_sphere_2& q, const Point_on_sphere_2& r) const + { return _lk.orientation_3_object()(_center, p, q, r); } protected: const Point_3& _center; - const Traits& _traits; + const LK& _lk; }; -template +template class Equal_on_sphere_2 { public: - typedef typename Traits::Point_3 Point_3; + typedef typename LK::Point_3 Point_3; + typedef typename LK::Point_3 Point_on_sphere_2; typedef bool result_type; - Equal_on_sphere_2(const Point_3& center, const Traits& traits) - : _center(center), _traits(traits) + Equal_on_sphere_2(const Point_3& center, const LK& lk) + : _center(center), _lk(lk) { } - bool operator()(const Point_3& p, const Point_3 q) const + bool operator()(const Point_on_sphere_2& p, const Point_on_sphere_2& q) const { - return _traits.collinear_3_object()(_center, p, q) && - !_traits.collinear_are_ordered_along_line_3_object()(p, _center, q); // @fixme correct? strictly? + return _lk.collinear_3_object()(_center, p, q) && + !_lk.collinear_are_strictly_ordered_along_line_3_object()(p, _center, q); } protected: const Point_3& _center; - const Traits& _traits; + const LK& _lk; }; -template +template class Inside_cone_2 { public: - typedef typename Traits::Point_3 Point_3; + typedef typename LK::Point_3 Point_3; + typedef typename LK::Point_3 Point_on_sphere_2; typedef bool result_type; - Inside_cone_2(const Point_3& center, const Traits& traits) - : _center(center), _traits(traits) + Inside_cone_2(const Point_3& center, const LK& lk) + : _center(center), _lk(lk) { } - bool operator()(const Point_3& p, const Point_3& q, const Point_3& r) const + bool operator()(const Point_on_sphere_2& p, const Point_on_sphere_2& q, const Point_on_sphere_2& r) const { - if(collinear(_center, p, r) || collinear(_center, q, r) || orientation(_center, p, q, r) != COLLINEAR) + // @todo first two checks might be unnecessary because we should always have r != p|q (using equal_on_sphere) + if(collinear(_center, p, r) || // @todo use LK + collinear(_center, q, r) || + orientation(_center, p, q, r) != COLLINEAR) + { return false; + } if(collinear(_center, p, q)) return true; - // @fixme (?) what's going on here? - return _traits.coplanar_orientation_3_object()(_center, p, q, r) == - (POSITIVE == _traits.coplanar_orientation_3_object()(_center, q, p, r)); + const Orientation op = _lk.coplanar_orientation_3_object()(_center, p, q, r); + const Orientation oq = _lk.coplanar_orientation_3_object()(_center, q, p, r); // @todo add an early exit + CGAL_assertion(op != COLLINEAR && oq != COLLINEAR); + + return (op == POSITIVE) && (oq == POSITIVE); } protected: const Point_3& _center; - const Traits& _traits; + const LK& _lk; }; +template +class Construct_arc_on_sphere_2 +{ +public: + typedef typename LK::FT FT; + typedef typename LK::Point_3 Point_3; + typedef typename LK::Point_3 Point_on_sphere_2; + typedef typename SK::Circular_arc_3 result_type; + + Construct_arc_on_sphere_2(const Point_3& center, const FT radius, const LK& lk, const SK& sk) + : _center(center), _radius(radius), _lk(lk), _sk(sk) + { } + + result_type operator()(const Point_on_sphere_2& p, const Point_on_sphere_2& q) const + { + // the orientation of the plane (and thus of the circle) does not really matter + // since the construction of the circular arc uses the positive normal of the circle's + // supporting plane... + typedef typename CGAL::internal::Converter_selector::type Converter; + Converter cv; + + typename SK::Point_3 sc = cv(_center); // @tmp cache that, I guess + typename SK::Point_3 sp = cv(p); + typename SK::Point_3 sq = cv(q); + + typename SK::Plane_3 pl = _sk.construct_plane_3_object()(sp, sc, sq); + typename SK::Circle_3 c = _sk.construct_circle_3_object()(sc, square(_radius), pl); + + typename SK::Circular_arc_point_3 cp = _sk.construct_circular_arc_point_3_object()(sp); + typename SK::Circular_arc_point_3 cq = _sk.construct_circular_arc_point_3_object()(sq); + + // @fixme ensure in arc_dual and arc_segment that 'cp' and 'cq' are in the correct order + return _sk.construct_circular_arc_3_object()(c, cp, cq); + } + +protected: + const Point_3& _center; + const FT _radius; + const LK& _lk; + const SK& _sk; +}; + + } // namespace internal template Equal_on_sphere_2; typedef internal::Orientation_on_sphere_2 Orientation_on_sphere_2; - typedef typename LK::Orientation_3 Side_of_oriented_circle_2; // @todo 'on_sphere'? + typedef typename LK::Orientation_3 Side_of_oriented_circle_on_sphere_2; // constructions + typedef typename LK::Construct_point_3 Construct_point_on_sphere_2; typedef typename LK::Construct_point_3 Construct_point_3; typedef typename LK::Construct_segment_3 Construct_segment_3; typedef typename LK::Construct_triangle_3 Construct_triangle_3; - typedef typename SK::Construct_circular_arc_3 Construct_arc_on_sphere_2; - // @fixme should project on the sphere - typedef typename LK::Construct_circumcenter_3 Construct_circumcenter_on_sphere_2; + // For the hilbert sort, not actually needed when the traits derives from LK @todo + typedef typename LK::Compute_x_3 Compute_x_3; + typedef typename LK::Compute_y_3 Compute_y_3; + typedef typename LK::Compute_z_3 Compute_z_3; + Compute_x_3 compute_x_3_object() const { return Compute_x_3(); } + Compute_y_3 compute_y_3_object() const { return Compute_y_3(); } + Compute_z_3 compute_z_3_object() const { return Compute_z_3(); } - typedef typename LK::Compute_squared_distance_3 Compute_squared_distance_3; + typedef internal::Construct_arc_on_sphere_2 Construct_arc_on_sphere_2; + + // @todo needs to offer dual_on_sphere that projects on the sphere + typedef void Construct_circumcenter_on_sphere_2; + + typedef typename LK::Construct_circumcenter_3 Construct_circumcenter_3; Delaunay_triangulation_sphere_traits_2(const Point_3& center = CGAL::ORIGIN, const FT radius = 1, @@ -145,9 +210,15 @@ public: private: void initialize_bounds() { - const FT minDist = _radius * std::pow(2, -23); // @fixme + // @fixme + // - check the correctness of the implementation + // - if LK can represent algebraic coordinates, there is no need for that + const FT minDist = _radius * std::pow(2, -23); const FT minRadius = _radius * (1 - std::pow(2, -50)); const FT maxRadius = _radius * (1 + std::pow(2, -50)); + std::cout << "minDist = " << minDist << std::endl; + std::cout << "minRadius = " << minRadius << std::endl; + std::cout << "maxRadius = " << maxRadius << std::endl; _minDistSquared = CGAL::square(minDist); _minRadiusSquared = CGAL::square(minRadius); _maxRadiusSquared = CGAL::square(maxRadius); @@ -170,11 +241,15 @@ public: orientation_on_sphere_2_object() const { return Orientation_on_sphere_2(_center, LK()); } - Side_of_oriented_circle_2 - side_of_oriented_circle_2_object() const + Side_of_oriented_circle_on_sphere_2 + side_of_oriented_circle_on_sphere_2_object() const { return typename LK::Orientation_3(); } // @tmp public: + Construct_point_on_sphere_2 + construct_point_on_sphere_2_object() const + { return Construct_point_on_sphere_2(); } // @tmp + Construct_point_3 construct_point_3_object() const { return typename LK::Construct_point_3(); } // @tmp @@ -189,15 +264,15 @@ public: Construct_arc_on_sphere_2 construct_arc_on_sphere_2_object() const - { return typename SK::Construct_circular_arc_3(); } // @tmp + { return Construct_arc_on_sphere_2(_center, _radius, LK(), SK()); } // @tmp Construct_circumcenter_on_sphere_2 construct_circumcenter_on_sphere_2_object() const { return typename LK::Construct_circumcenter_3(); } // @tmp - Compute_squared_distance_3 - compute_squared_distance_3_object() const - { return typename LK::Compute_squared_distance_3(); } // @tmp + Construct_circumcenter_3 + construct_circumcenter_3_object() const + { return typename LK::Construct_circumcenter_3(); } // @tmp public: void set_center(const Point_3& center) { _center = center; } @@ -206,15 +281,16 @@ public: FT radius() const { return _radius; } public: - bool is_on_sphere(const Point_3& p) const + bool is_on_sphere(const Point_on_sphere_2& p) const { - const FT sq_dist = compute_squared_distance_3_object()(p, _center); - return (_minRadiusSquared < sq_dist && sq_dist < _maxRadiusSquared); + const FT sq_dist = CGAL::squared_distance(p, _center); // @todo use LK + std::cout << _minRadiusSquared << " " << sq_dist << " " << _maxRadiusSquared << std::endl; + return (_minRadiusSquared <= sq_dist && sq_dist < _maxRadiusSquared); } - bool are_points_too_close(const Point_3& p, const Point_3& q) const + bool are_points_too_close(const Point_on_sphere_2& p, const Point_on_sphere_2& q) const { - return (compute_squared_distance_3_object()(p, q) <= _minDistSquared); + return (CGAL::squared_distance(p, q) <= _minDistSquared); } protected: diff --git a/Triangulation_on_sphere_2/include/CGAL/Projection_sphere_traits_3.h b/Triangulation_on_sphere_2/include/CGAL/Projection_sphere_traits_3.h index cb1fc267521..d45876cc769 100644 --- a/Triangulation_on_sphere_2/include/CGAL/Projection_sphere_traits_3.h +++ b/Triangulation_on_sphere_2/include/CGAL/Projection_sphere_traits_3.h @@ -20,144 +20,190 @@ #include namespace CGAL { +namespace internal { -template < typename K> +template class Point_with_scale - : public K::Point_3 + : public LK::Point_3 { public: - typedef typename K::FT FT; - typedef typename K::Point_3 Base_point; + typedef typename LK::FT FT; + typedef typename LK::Point_3 Base_point; + typedef typename LK::Vector_3 Vector_3; +public: Point_with_scale() : Base_point() { } - Point_with_scale(const Base_point& p, - const typename K::Point_3& sphere_center) + Point_with_scale(const Base_point& p) : Base_point(p) - { - compute_scale(p-(sphere_center-ORIGIN)); - } - - FT scale() const { return _scale; } - -private: - void compute_scale(const FT x, const FT y, const FT z) - { - const FT tmp = x*x + y*y + z*z; - if(tmp == 0) - _scale = FT(0); - else - _scale = FT(1) / CGAL::approximate_sqrt(tmp); - } - - void compute_scale(const Base_point& p) - { - return compute_scale(p.x(), p.y(), p.z()); - } - - FT _scale; -}; - -//adaptor for calling the Predicate_ with the points projected on the sphere -template -class Functor_projection_adaptor - : public Predicate_ -{ - typedef Predicate_ Predicate; - typedef Predicate_ Base; - - typedef typename Traits::FT FT; - typedef typename Traits::Point_on_sphere_2 Point; - typedef typename Traits::Point_3 Point_3; - -public: - typedef typename Predicate::result_type result_type; - - Functor_projection_adaptor(const Predicate& p, const FT r) : Base(p), _radius(r) { } - - result_type operator()(const Point& p0, const Point& p1) - { return Base::operator()(project(p0), project(p1)); } - - result_type operator()(const Point& p0, const Point& p1, const Point& p2) - { return Base::operator()(project(p0), project(p1), project(p2)); } - - result_type operator ()(const Point& p0, const Point& p1, const Point& p2, const Point& p3) - { return Base::operator()(project(p0), project(p1), project(p2), project(p3)); } - - result_type operator()(const Point& p0, const Point& p1, const Point& p2, const Point& p3, const Point& p4) - { return Base::operator()(project(p0), project(p1), project(p2), project(p3), project(p4)); } - -private: - Base_point project(const Point& p) - { - const FT scale = _radius * p.scale(); - return Base_point(scale*p.x(), scale*p.y(), scale*p.z()); - } - - const FT _radius; -}; - -template -class Projection_sphere_traits_3 - : public Delaunay_triangulation_sphere_traits_2 -{ -protected: - typedef Delaunay_triangulation_sphere_traits_2 Base; - typedef Projection_sphere_traits_3 Self; - -public: - typedef typename K::FT FT; - typedef Point_with_scale Point_on_sphere_2; - typedef typename K::Point_3 Point_3; - - typedef Functor_projection_adaptor Construct_circumcenter_on_sphere_2; - typedef Functor_projection_adaptor Equal_on_sphere_2; - typedef Functor_projection_adaptor Inside_cone_2; - typedef Functor_projection_adaptor Orientation_on_sphere_2; - typedef Functor_projection_adaptor Power_test_2; - - Projection_sphere_traits_3(const Point_3& sphere = CGAL::ORIGIN, - const FT radius = 1, - const K& k = K()) - : Base(sphere, radius, k) { } public: - Construct_circumcenter_on_sphere_2 - construct_circumcenter_on_sphere_2_object() const - { return Construct_circumcenter_on_sphere_2(Base::construct_circumcenter_on_sphere_2(), Base::_radius); } + Base_point project(const Base_point& center, + const FT radius) const + { + CGAL_precondition(static_cast(*this) != center); + + // @todo use LK + Vector_3 v = static_cast(*this) - center; + v = (radius / CGAL::approximate_sqrt(v * v)) * v; + return center + v; + } +}; + +template +struct Construct_point_with_scale + : public std::unary_function +{ + PoS2 operator()(const P3& pt) const { return PoS2(pt); } +}; + +// Adaptor for calling the functors with the points projected on the sphere +// +// @todo filter this +template +class Functor_projection_adaptor + : public Functor_ +{ + typedef Functor_ Functor; + typedef Functor_ Base; + + typedef typename LK::FT FT; + typedef typename LK::Point_3 Point_3; + typedef internal::Point_with_scale Point; + +public: + Functor_projection_adaptor(const Functor& f, const Point_3& center, const FT radius) + : Base(f), _c(center), _r(radius) + { } + +public: + using Base::operator(); + + decltype(auto) operator()(const Point& p0, const Point& p1) + { return Base::operator()(p0.project(_c, _r), p1.project(_c, _r)); } + + decltype(auto) operator()(const Point& p0, const Point& p1, const Point& p2) + { return Base::operator()(p0.project(_c, _r), p1.project(_c, _r), p2.project(_c, _r)); } + + decltype(auto) operator ()(const Point& p0, const Point& p1, const Point& p2, const Point& p3) + { return Base::operator()(p0.project(_c, _r), p1.project(_c, _r), p2.project(_c, _r), p3.project(_c, _r)); } + + decltype(auto) operator()(const Point& p0, const Point& p1, const Point& p2, const Point& p3, const Point& p4) + { return Base::operator()(p0.project(_c, _r), p1.project(_c, _r), p2.project(_c, _r), p3.project(_c, _r), p4.project(_c, _r)); } + +private: + const Point_3& _c; + const FT _r; +}; + +} // namespace internal + +template +class Projection_sphere_traits_3 + : public Delaunay_triangulation_sphere_traits_2 +{ +protected: + typedef Delaunay_triangulation_sphere_traits_2 Base; + typedef Projection_sphere_traits_3 Self; + +public: + typedef typename LK::FT FT; + typedef typename LK::Point_3 Point_3; + typedef internal::Point_with_scale Point_on_sphere_2; + + // predicates + typedef internal::Functor_projection_adaptor Collinear_are_strictly_ordered_on_great_circle_2; + typedef internal::Functor_projection_adaptor Compare_on_sphere_2; + typedef internal::Functor_projection_adaptor Equal_on_sphere_2; + typedef internal::Functor_projection_adaptor Orientation_on_sphere_2; + typedef internal::Functor_projection_adaptor Side_of_oriented_circle_on_sphere_2; + + // constructions + typedef internal::Construct_point_with_scale Construct_point_on_sphere_2; + typedef internal::Functor_projection_adaptor Construct_point_3; + typedef internal::Functor_projection_adaptor Construct_segment_3; + typedef internal::Functor_projection_adaptor Construct_triangle_3; + + typedef internal::Functor_projection_adaptor Construct_arc_on_sphere_2; + typedef internal::Functor_projection_adaptor Construct_circumcenter_on_sphere_2; + +public: + Projection_sphere_traits_3(const Point_3& sphere = CGAL::ORIGIN, + const FT radius = 1, + const LK& lk = LK()) + : Base(sphere, radius, lk) + { } + + // predicates +public: + Collinear_are_strictly_ordered_on_great_circle_2 + collinear_are_strictly_ordered_on_great_circle_2_object() const + { return Collinear_are_strictly_ordered_on_great_circle_2( + Base::collinear_are_strictly_ordered_on_great_circle_2_object(), + Base::center(), Base::radius()); } + + Compare_on_sphere_2 + compare_on_sphere_2_object() const + { return Compare_on_sphere_2(Base::compare_on_sphere_2_object(), + Base::center(), Base::radius()); } Equal_on_sphere_2 equal_on_sphere_2_object() const - { return Equal_on_sphere_2(Base::equal_on_sphere_2(), Base::_radius); } - - Inside_cone_2 - inside_cone_2_object() const - { return Inside_cone_2(Base::inside_cone_2_object(), Base::_radius); } + { return Equal_on_sphere_2(Base::equal_on_sphere_2_object(), + Base::center(), Base::radius()); } Orientation_on_sphere_2 - orientation_on_sphere_2() const - { return Orientation_2(Base::orientation_on_sphere_2_object(), Base::_radius); } + orientation_on_sphere_2_object() const + { return Orientation_on_sphere_2(Base::orientation_on_sphere_2_object(), + Base::center(), Base::radius()); } - Power_test_2 - power_test_2_object() const - { return Power_test_2(Base::power_test_2_object(), Base::_radius); } + Side_of_oriented_circle_on_sphere_2 + side_of_oriented_circle_on_sphere_2_object() const + { return Side_of_oriented_circle_on_sphere_2(Base::side_of_oriented_circle_on_sphere_2_object(), + Base::center(), Base::radius()); } + + // constructions +public: + Construct_point_on_sphere_2 + construct_point_on_sphere_2_object() const + { return Construct_point_on_sphere_2(); } + + Construct_point_3 + construct_point_3_object() const + { return Construct_point_3(Base::construct_point_3_object(), + Base::center(), Base::radius()); } + + Construct_segment_3 + construct_segment_3_object() const + { return Construct_segment_3(Base::construct_segment_3_object(), + Base::center(), Base::radius()); } + + Construct_triangle_3 + construct_triangle_3_object() const + { return Construct_triangle_3(Base::construct_triangle_3_object(), + Base::center(), Base::radius()); } + + Construct_arc_on_sphere_2 + construct_arc_on_sphere_2_object() const + { return Construct_arc_on_sphere_2(Base::construct_arc_on_sphere_2_object(), + Base::center(), Base::radius()); } + + Construct_circumcenter_on_sphere_2 + construct_circumcenter_on_sphere_2_object() const + { return Construct_circumcenter_on_sphere_2(Base::construct_circumcenter_on_sphere_2_object(), + Base::center(), Base::radius()); } public: - // @fixme rename that? - struct Construct_projected_point_3 - : public std::unary_function + bool is_on_sphere(const Point_on_sphere_2& /*p*/) const { - Construct_projected_point_3(const Point_3& sc) : sphere_center(sc) { } + return true; + } - Point_on_sphere_2 operator()(const Point_3& pt) const { return Point_on_sphere_2(pt, sphere_center); } - - private: - const Point_3& sphere_center; - }; - - Construct_projected_point_3 - construct_projected_point_3_object() const - { return Construct_projected_point_3(_center); } + bool are_points_too_close(const Point_on_sphere_2& p, const Point_on_sphere_2& q) const + { + return Base::are_points_too_close(p.project(Base::center(), Base::radius()), + q.project(Base::center(), Base::radius())); + } }; } // namespace CGAL diff --git a/Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h b/Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h index 130b16851a2..63083269c39 100644 --- a/Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h +++ b/Triangulation_on_sphere_2/include/CGAL/Triangulation_on_sphere_2.h @@ -452,7 +452,7 @@ Orientation Triangulation_on_sphere_2:: orientation(const Point&p, const Point& q, const Point& r, const Point& s) const { - return geom_traits().side_of_oriented_circle_2_object()(p, q, r, s); + return geom_traits().side_of_oriented_circle_on_sphere_2_object()(p, q, r, s); } template