Update traits classes to new spec

This commit is contained in:
Mael Rouxel-Labbé 2021-01-29 12:22:16 +01:00
parent dfd59a5504
commit e6dba91060
3 changed files with 283 additions and 161 deletions

View File

@ -17,6 +17,7 @@
#include <CGAL/license/Triangulation_on_sphere_2.h>
#include <CGAL/Algebraic_kernel_for_spheres_2_3.h>
#include <CGAL/Has_conversion.h>
#include <CGAL/Spherical_kernel_3.h>
#include <CGAL/enum.h>
@ -24,76 +25,129 @@
namespace CGAL {
namespace internal {
template <typename Traits>
template <typename LK>
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 <typename Traits>
template <typename LK>
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 <typename Traits>
template <typename LK>
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 <typename LK, typename SK>
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<LK, SK>::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 <typename LK_,
@ -112,6 +166,7 @@ public:
typedef typename LK::Point_3 Point_on_sphere_2;
typedef typename LK::Point_3 Point_3;
typedef typename LK::Segment_3 Segment_3;
typedef typename LK::Triangle_3 Triangle_3;
typedef typename SK::Circular_arc_3 Arc_on_sphere_2;
@ -120,18 +175,28 @@ public:
typedef typename LK::Compare_xyz_3 Compare_on_sphere_2;
typedef internal::Equal_on_sphere_2<LK> Equal_on_sphere_2;
typedef internal::Orientation_on_sphere_2<LK> 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<LK, SK> 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:

View File

@ -20,144 +20,190 @@
#include <CGAL/number_utils_classes.h>
namespace CGAL {
namespace internal {
template < typename K>
template <typename LK>
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 <typename Traits, typename Predicate_>
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 <typename K>
class Projection_sphere_traits_3
: public Delaunay_triangulation_sphere_traits_2<K>
{
protected:
typedef Delaunay_triangulation_sphere_traits_2<K> Base;
typedef Projection_sphere_traits_3<K> Self;
public:
typedef typename K::FT FT;
typedef Point_with_scale<K> Point_on_sphere_2;
typedef typename K::Point_3 Point_3;
typedef Functor_projection_adaptor<Self, typename Base::Construct_circumcenter_on_sphere_2> Construct_circumcenter_on_sphere_2;
typedef Functor_projection_adaptor<Self, typename Base::Equal_on_sphere_2> Equal_on_sphere_2;
typedef Functor_projection_adaptor<Self, typename Base::Inside_cone_2> Inside_cone_2;
typedef Functor_projection_adaptor<Self, typename Base::Orientation_on_sphere_2> Orientation_on_sphere_2;
typedef Functor_projection_adaptor<Self, typename Base::Power_test_2> 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<const Base_point&>(*this) != center);
// @todo use LK
Vector_3 v = static_cast<const Base_point&>(*this) - center;
v = (radius / CGAL::approximate_sqrt(v * v)) * v;
return center + v;
}
};
template <typename P3, typename PoS2>
struct Construct_point_with_scale
: public std::unary_function<P3, PoS2>
{
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 <typename LK, typename Functor_>
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<LK> 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 <typename LK>
class Projection_sphere_traits_3
: public Delaunay_triangulation_sphere_traits_2<LK>
{
protected:
typedef Delaunay_triangulation_sphere_traits_2<LK> Base;
typedef Projection_sphere_traits_3<LK> Self;
public:
typedef typename LK::FT FT;
typedef typename LK::Point_3 Point_3;
typedef internal::Point_with_scale<LK> Point_on_sphere_2;
// predicates
typedef internal::Functor_projection_adaptor<LK, typename Base::Collinear_are_strictly_ordered_on_great_circle_2> Collinear_are_strictly_ordered_on_great_circle_2;
typedef internal::Functor_projection_adaptor<LK, typename Base::Compare_on_sphere_2> Compare_on_sphere_2;
typedef internal::Functor_projection_adaptor<LK, typename Base::Equal_on_sphere_2> Equal_on_sphere_2;
typedef internal::Functor_projection_adaptor<LK, typename Base::Orientation_on_sphere_2> Orientation_on_sphere_2;
typedef internal::Functor_projection_adaptor<LK, typename Base::Side_of_oriented_circle_on_sphere_2> Side_of_oriented_circle_on_sphere_2;
// constructions
typedef internal::Construct_point_with_scale<Point_3, Point_on_sphere_2> Construct_point_on_sphere_2;
typedef internal::Functor_projection_adaptor<LK, typename Base::Construct_point_3> Construct_point_3;
typedef internal::Functor_projection_adaptor<LK, typename Base::Construct_segment_3> Construct_segment_3;
typedef internal::Functor_projection_adaptor<LK, typename Base::Construct_triangle_3> Construct_triangle_3;
typedef internal::Functor_projection_adaptor<LK, typename Base::Construct_arc_on_sphere_2> Construct_arc_on_sphere_2;
typedef internal::Functor_projection_adaptor<LK, typename Base::Construct_circumcenter_on_sphere_2> 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<Point_3, Point_on_sphere_2>
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

View File

@ -452,7 +452,7 @@ Orientation
Triangulation_on_sphere_2<Gt, Tds>::
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 <typename Gt, typename Tds>