songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Cleaned up; eliminated local kernels

This commit is contained in:
Efi Fogel 2024-01-29 21:55:39 +02:00
parent 1c7cdad3f0
commit 7381e39e1f
4 changed files with 1090 additions and 1188 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

482
Envelope_3/include/CGAL/Env_plane_traits_3.h
View File

@ -7,7 +7,8 @@
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Baruch Zukerman <baruchzu@post.tau.ac.il>
// Author(s) : Baruch Zukerman <baruchzu@post.tau.ac.il>
// Efi Fogel <efifogel@gmail.com>
#ifndef CGAL_ENV_PLANE_TRAITS_3_H
#define CGAL_ENV_PLANE_TRAITS_3_H
@ -27,169 +28,149 @@
namespace CGAL {
template <class Kernel_>
class Env_plane_traits_3 : public Arr_linear_traits_2<Kernel_>
{
template <typename Kernel_>
class Env_plane_traits_3 : public Arr_linear_traits_2<Kernel_> {
public:
typedef Kernel_ Kernel;
typedef typename Kernel::FT FT;
typedef Arr_linear_traits_2<Kernel> Base;
typedef Env_plane_traits_3<Kernel> Self;
using Kernel = Kernel_;
using FT = typename Kernel::FT;
using Base = Arr_linear_traits_2<Kernel>;
using Self = Env_plane_traits_3<Kernel>;
typedef typename Base::Multiplicity Multiplicity;
typedef typename Base::Point_2 Point_2;
typedef typename Base::Curve_2 Curve_2;
typedef typename Base::X_monotone_curve_2 X_monotone_curve_2;
typedef typename Kernel::Plane_3 Plane_3;
typedef typename Kernel::Vector_2 Vector_2;
typedef typename Kernel::Vector_3 Vector_3;
typedef typename Kernel::Segment_2 Segment_2;
typedef typename Kernel::Ray_2 Ray_2;
typedef typename Kernel::Line_2 Line_2;
typedef typename Kernel::Line_3 Line_3;
typedef std::pair<Curve_2, Multiplicity> Intersection_curve;
using Multiplicity = typename Base::Multiplicity;
using Point_2 = typename Base::Point_2;
using Curve_2 = typename Base::Curve_2;
using X_monotone_curve_2 = typename Base::X_monotone_curve_2;
using Plane_3 = typename Kernel::Plane_3;
using Vector_2 = typename Kernel::Vector_2;
using Vector_3 = typename Kernel::Vector_3;
using Segment_2 = typename Kernel::Segment_2;
using Ray_2 = typename Kernel::Ray_2;
using Line_2 = typename Kernel::Line_2;
using Line_3 = typename Kernel::Line_3;
using Intersection_curve = std::pair<Curve_2, Multiplicity>;
typedef typename Base::Left_side_category Left_side_category;
typedef typename Base::Bottom_side_category Bottom_side_category;
typedef typename Base::Top_side_category Top_side_category;
typedef typename Base::Right_side_category Right_side_category;
using Left_side_category = typename Base::Left_side_category;
using Bottom_side_category = typename Base::Bottom_side_category;
using Top_side_category = typename Base::Top_side_category;
using Right_side_category = typename Base::Right_side_category;
class Is_vertical_3
{
//!
class Is_vertical_3 {
public:
bool operator()(const Plane_3& h) const
{
return CGAL::is_zero(h.c());
}
bool operator()(const Plane_3& h) const { return CGAL::is_zero(h.c()); }
};
Is_vertical_3 is_vertical_3_object() const
{
return Is_vertical_3();
}
//! Obtain an Is_vertical_3 functor object.
Is_vertical_3 is_vertical_3_object() const { return Is_vertical_3(); }
class _Env_plane
{
//
class _Env_plane {
protected:
Plane_3 m_plane;
Line_2 m_line;
bool m_is_all_plane; // true -> all plane, false -> halfplane
bool m_is_vert;
Plane_3 m_plane;
Line_2 m_line;
bool m_is_all_plane; // true -> all plane, false -> halfplane
bool m_is_vert;
public:
_Env_plane()
{}
_Env_plane() {}
_Env_plane(const Plane_3& h) : m_plane(h),
m_is_all_plane(true)
{
_Env_plane(const Plane_3& h) : m_plane(h), m_is_all_plane(true) {
Self s;
m_is_vert = s.is_vertical_3_object()(h);
}
_Env_plane(const Plane_3& h, const Line_2& l) : m_plane(h),
m_line(l),
m_is_all_plane(false),
m_is_vert(false)
{
_Env_plane(const Plane_3& h, const Line_2& l) :
m_plane(h), m_line(l), m_is_all_plane(false), m_is_vert(false) {
CGAL_precondition_code(Self s);
CGAL_precondition(!s.is_vertical_3_object()(h));
CGAL_precondition(! s.is_vertical_3_object()(h));
}
bool is_vertical() const
{
return m_is_vert;
}
bool is_vertical() const { return m_is_vert; }
const Plane_3& plane() const
{
return m_plane;
}
const Plane_3& plane() const { return m_plane; }
operator Plane_3() const { return m_plane; }
operator Plane_3 () const
{
return (m_plane);
}
const Line_2& line() const
{
CGAL_assertion(!m_is_all_plane);
const Line_2& line() const {
CGAL_assertion(! m_is_all_plane);
return m_line;
}
bool is_all_plane() const
{
return m_is_all_plane;
}
bool is_all_plane() const { return m_is_all_plane; }
};
typedef _Env_plane Xy_monotone_surface_3;
typedef _Env_plane Surface_3;
using Xy_monotone_surface_3 = _Env_plane;
using Surface_3 = _Env_plane;
class Make_xy_monotone_3
{
/*! Subdivide a given surface into \f$xy\f$-monotone parts.
*/
class Make_xy_monotone_3 {
public:
template <class OutputIterator>
OutputIterator operator()(const Surface_3& s,
bool /* is_lower */,
OutputIterator o) const
{
template <typename OutputIterator>
OutputIterator operator()(const Surface_3& s, bool /* is_lower */,
OutputIterator o) const {
*o++ = s;
return o;
}
};
//! Obtain a Make_xy_monotone_3 functor object.
Make_xy_monotone_3 make_xy_monotone_3_object() const
{
return Make_xy_monotone_3();
}
{ return Make_xy_monotone_3(); }
/*! Determine the relative \f$z\f$-order of two given \f$xy\f$-monotone
* surfaces at the \f$xy\f$-coordinates of a given point or \f$x\f$-monotone
* curve.
*/
class Compare_z_at_xy_3 {
protected:
using Traits_3 = Env_plane_traits_3<Kernel>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Compare_z_at_xy_3(const Traits_3& traits) : m_traits(traits) {}
friend class Env_plane_traits_3<Kernel>;
class Compare_z_at_xy_3
{
public:
//
Comparison_result operator()(const Point_2& p,
const Xy_monotone_surface_3& h1,
const Xy_monotone_surface_3& h2) const
{
const Xy_monotone_surface_3& h2) const {
const Plane_3& plane1 = h1.plane();
const Plane_3& plane2 = h2.plane();
Sign sign_of_c1c2 = CGAL::sign(plane1.c() * plane2.c());
Sign sign_of_expr =
CGAL::sign ((p.x()*plane1.a() + p.y()*plane1.b() +
plane1.d())*plane2.c() -
(p.x()*plane2.a() + p.y()*plane2.b() +
plane2.d())*plane1.c());
CGAL::sign((p.x()*plane1.a() + p.y()*plane1.b() +
plane1.d())*plane2.c() -
(p.x()*plane2.a() + p.y()*plane2.b() +
plane2.d())*plane1.c());
int i = -1 * static_cast<int>(sign_of_c1c2) *
static_cast<int>(sign_of_expr);
return static_cast<Comparison_result>(i);
}
//
Comparison_result operator()(const X_monotone_curve_2& cv,
const Xy_monotone_surface_3& h1,
const Xy_monotone_surface_3& h2) const
{
Kernel k;
Point_2 p;
if(cv.is_segment())
p = k.construct_midpoint_2_object()(cv.left(), cv.right());
else
if(cv.is_ray())
p = k.construct_point_on_2_object()(cv.ray(), 1);
else
{
CGAL_assertion(cv.is_line());
p = k.construct_point_on_2_object()(cv.line(), 1);
}
const Xy_monotone_surface_3& h2) const {
const Kernel& kernel = m_traits;
Point_2 p = (cv.is_segment()) ?
kernel.construct_midpoint_2_object()(cv.left(), cv.right()) :
((cv.is_ray()) ?
kernel.construct_point_on_2_object()(cv.ray(), 1) :
kernel.construct_point_on_2_object()(cv.line(), 1));
return this->operator()(p, h1, h2);
}
//
Comparison_result operator()(const Xy_monotone_surface_3& h1,
const Xy_monotone_surface_3& h2) const
{
const Xy_monotone_surface_3& h2) const {
CGAL_assertion(h1.is_all_plane() && h2.is_all_plane());
const Plane_3& p1 = h1.plane();
@ -201,35 +182,63 @@ public:
}
};
//! Obtain a Compare_z_at_xy_3 functor object.
Compare_z_at_xy_3 compare_z_at_xy_3_object() const
{
return Compare_z_at_xy_3();
}
{ return Compare_z_at_xy_3(*this); }
/*! Determine the relative \f$z\f$-order of the two given \f$xy\f$-monotone
* surfaces immediately above their projected intersection curve (a planar
* point \f$p\f$ is above an \f$x\f$-monotone curve \f$c\f$ if it is in the
* \f$x\f$-range of \f$c\f$, and lies to its left when the curve is traversed
* from its \f$xy\f$-lexicographically smaller endpoint to its larger
* endpoint).
*/
class Compare_z_at_xy_above_3 {
protected:
using Traits_3 = Env_plane_traits_3<Kernel>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Compare_z_at_xy_above_3(const Traits_3& traits) : m_traits(traits) {}
friend class Env_plane_traits_3<Kernel>;
class Compare_z_at_xy_above_3
{
public:
/*! Determine the relative \f$z\f$-order of the two given \f$xy\f$-monotone
* surfaces immediately above their projected intersection curve, which is
* also given.
*
* \param cv the intersection curve.
* \param h1 the first surface.
* \param h2 the second surface.
* \pre `h1` and `h2` are defined "above" `cv`, and their relative
* \f$z\f$-order is the same for some small enough neighborhood of points
* above `cv`.
*/
Comparison_result operator()(const X_monotone_curve_2& cv,
const Xy_monotone_surface_3& h1,
const Xy_monotone_surface_3& h2) const
{
const Xy_monotone_surface_3& h2) const {
const Plane_3& plane1 = h1.plane();
const Plane_3& plane2 = h2.plane();
const FT& a1 = plane1.a(),
b1 = plane1.b(),
c1 = plane1.c();
const FT& a1 = plane1.a();
const FT& b1 = plane1.b();
const FT& c1 = plane1.c();
const FT& a2 = plane2.a(),
b2 = plane2.b(),
c2 = plane2.c();
const FT& a2 = plane2.a();
const FT& b2 = plane2.b();
const FT& c2 = plane2.c();
// our line is a3*x + b3*y + c3 = 0
// it is assumed that the planes intersect over this line
const Line_2& line = cv.supp_line();
const FT& a3 = line.a(),
b3 = line.b(),
c3 = line.c();
const FT& a3 = line.a();
const FT& b3 = line.b();
const FT& c3 = line.c();
// if the line was parallel to the y-axis (i.e x = const),
// then it was enough to compare dz/dx of both planes
@ -266,62 +275,97 @@ public:
// are transformed to (v1,w1) and (v2,w2), so we need that w2 > w1
// (otherwise the result should be multiplied by -1)
Kernel k;
Point_2 p1 (k.construct_point_on_2_object()(line, 0));
Point_2 p2 (k.construct_point_on_2_object()(line, 1));
const Kernel& kernel = m_traits;
Point_2 p1 = kernel.construct_point_on_2_object()(line, 0);
Point_2 p2 = kernel.construct_point_on_2_object()(line, 1);
if(k.compare_xy_2_object()(p1, p2) == LARGER)
std::swap(p1, p2);
if (kernel.compare_xy_2_object()(p1, p2) == LARGER) std::swap(p1, p2);
CGAL_assertion(k.compare_xy_2_object()(p1, p2) == SMALLER);
CGAL_assertion(kernel.compare_xy_2_object()(p1, p2) == SMALLER);
const FT& x1 = p1.x(),
y1 = p1.y(),
x2 = p2.x(),
y2 = p2.y();
const FT& x1 = p1.x();
const FT& y1 = p1.y();
const FT& x2 = p2.x();
const FT& y2 = p2.y();
Sign s2 = CGAL_NTS sign(-b3*x1+a3*y1-(-b3*x2+a3*y2));
return s1 * s2;
}
};
//! Obtain a Compare_z_at_xy_above_3 functor object.
Compare_z_at_xy_above_3 compare_z_at_xy_above_3_object() const
{
return Compare_z_at_xy_above_3();
}
{ return Compare_z_at_xy_above_3(*this); }
/*! Determine the relative \f$z\f$-order of the two given \f$xy\f$-monotone
* surfaces immediately below their projected intersection curve (a planar
* point \f$p\f$ is below an \f$x\f$-monotone curve \f$c\f$ if it is in the
* \f$x\f$-range of \f$c\f$, and lies to its left when the curve is traversed
* from its \f$xy\f$-lexicographically smaller endpoint to its larger
* endpoint).
*/
class Compare_z_at_xy_below_3 {
protected:
using Traits_3 = Env_plane_traits_3<Kernel>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Compare_z_at_xy_below_3(const Traits_3& traits) : m_traits(traits) {}
friend class Env_plane_traits_3<Kernel>;
class Compare_z_at_xy_below_3
{
public:
/*! Determine the relative \f$z\f$-order of the two given \f$xy\f$-monotone
* surfaces immediately below their projected intersection curve, which is
* also given.
*
* \param cv the intersection curve.
* \param h1 the first surface.
* \param h2 the second surface.
* \pre `h1` and `h2` are defined "above" `cv`, and their relative
* \f$z\f$-order is the same for some small enough neighborhood of points
* below `cv`.
*/
Comparison_result operator()(const X_monotone_curve_2& cv,
const Xy_monotone_surface_3& h1,
const Xy_monotone_surface_3& h2) const
{
Compare_z_at_xy_above_3 cmp_above;
const Xy_monotone_surface_3& h2) const {
auto cmp_above = m_traits.compare_z_at_xy_above_3_object();
return CGAL::opposite(cmp_above(cv, h1, h2));
}
};
//! Obtain a Compare_z_at_xy_below_3 functor object.
Compare_z_at_xy_below_3 compare_z_at_xy_below_3_object() const
{
return Compare_z_at_xy_below_3();
}
{ return Compare_z_at_xy_below_3(*this); }
/*! Compute all planar \f$x\f$-monotone curves and possibly isolated planar
* points that form the projection of the boundary of the given
* \f$xy\f$-monotone surface s onto the \f$xy\f$-plane.
*/
class Construct_projected_boundary_2 {
protected:
using Traits_3 = Env_plane_traits_3<Kernel>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Construct_projected_boundary_2(const Traits_3& traits) : m_traits(traits) {}
friend class Env_plane_traits_3<Kernel>;
class Construct_projected_boundary_2
{
public:
template <class OutputIterator>
template <typename OutputIterator>
OutputIterator operator()(const Xy_monotone_surface_3& s,
OutputIterator o) const
{
if(s.is_all_plane())
{
if(!s.is_vertical())
return o;
OutputIterator o) const {
if (s.is_all_plane()) {
if (! s.is_vertical()) return o;
const Plane_3& h = s.plane();
Line_2 proj_line(h.a(), h.b(), h.d());
@ -331,10 +375,10 @@ public:
}
// s is half-plane
Kernel k;
const Point_2& p1 = k.construct_point_on_2_object()(s.line(), 0);
const Point_2& p2 = k.construct_point_on_2_object()(s.line(), 1);
Comparison_result res = k.compare_xy_2_object()(p1, p2);
const Kernel& kernel = m_traits;
const Point_2& p1 = kernel.construct_point_on_2_object()(s.line(), 0);
const Point_2& p2 = kernel.construct_point_on_2_object()(s.line(), 1);
Comparison_result res = kernel.compare_xy_2_object()(p1, p2);
Oriented_side side =
(res == SMALLER) ? ON_POSITIVE_SIDE : ON_NEGATIVE_SIDE;
@ -343,118 +387,106 @@ public:
}
};
//
Construct_projected_boundary_2
construct_projected_boundary_2_object() const
{
return Construct_projected_boundary_2();
}
{ return Construct_projected_boundary_2(*this); }
/*! compute the projection of the intersections of the \f$xy\f$-monotone
* surfaces onto the \f$xy\f$-plane,
*/
class Construct_projected_intersections_2 {
protected:
using Traits_3 = Env_plane_traits_3<Kernel>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Construct_projected_intersections_2(const Traits_3& traits) :
m_traits(traits)
{}
friend class Env_plane_traits_3<Kernel>;
class Construct_projected_intersections_2
{
public:
template <class OutputIterator>
template <typename OutputIterator>
OutputIterator operator()(const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2,
OutputIterator o) const
{
Kernel k;
OutputIterator o) const {
const Kernel& kernel = m_traits;
const Plane_3& h1 = s1.plane();
const Plane_3& h2 = s2.plane();
if(s1.is_vertical() && s2.is_vertical())
{
if (s1.is_vertical() && s2.is_vertical()) {
Line_2 l1(h1.a(), h1.b(), h1.d());
Line_2 l2(h2.a(), h2.b(), h2.d());
auto obj = k.intersect_2_object()(l1, l2);
auto obj = kernel.intersect_2_object()(l1, l2);
if(const Point_2* p = std::get_if<Point_2>(&(*obj)))
*o++ = *p;
if (const auto* p = std::get_if<Point_2>(&(*obj))) *o++ = *p;
// otherwise, the vertical planes are parallel or overlap, so we return
// nothing.
return o;
}
if(s1.is_all_plane() && s2.is_all_plane())
{
auto obj = k.intersect_3_object()(h1, h2);
if (s1.is_all_plane() && s2.is_all_plane()) {
auto obj = kernel.intersect_3_object()(h1, h2);
CGAL_assertion(obj != std::nullopt);
if(const Line_3* l = std::get_if<Line_3>(&(*obj)))
*o++ = Intersection_curve(project_xy(*l, k), 1);
if (const auto* l = std::get_if<Line_3>(&(*obj)))
*o++ = Intersection_curve(project_xy(*l, kernel), 1);
return o;
}
if(s1.is_all_plane() && !s2.is_all_plane())
{
auto obj = plane_half_plane_proj_intersection(h1,
h2,
s2.line(),
k);
if(obj ==std::nullopt)
return o;
if(const Line_2* line = std::get_if<Line_2>(&(*obj)))
{
if (s1.is_all_plane() && ! s2.is_all_plane()) {
auto obj = plane_half_plane_proj_intersection(h1, h2, s2.line(), kernel);
if (obj == std::nullopt) return o;
if (const auto* line = std::get_if<Line_2>(&(*obj))) {
*o++ = Intersection_curve(*line, 1);
return o;
}
if(const Ray_2* ray = std::get_if<Ray_2>(&(*obj)))
{
if (const auto* ray = std::get_if<Ray_2>(&(*obj))) {
*o++ = Intersection_curve(*ray, 1);
return o;
}
return o;
}
if(!s2.is_all_plane() && s2.is_all_plane())
{
auto obj = plane_half_plane_proj_intersection(h2,
h1,
s1.line(),
k);
if(obj == std::nullopt)
return o;
if(const Line_2* line = std::get_if<Line_2>(&(*obj)))
{
if (! s2.is_all_plane() && s2.is_all_plane()) {
auto obj = plane_half_plane_proj_intersection(h2, h1, s1.line(), kernel);
if (obj == std::nullopt) return o;
if (const auto* line = std::get_if<Line_2>(&(*obj))) {
*o++ = Intersection_curve(*line, 1);
return o;
}
if(const Ray_2* ray = std::get_if<Ray_2>(&(*obj)))
{
if (const auto* ray = std::get_if<Ray_2>(&(*obj))) {
*o++ = Intersection_curve(*ray, 1);
return o;
}
return o;
}
CGAL_assertion(!s2.is_all_plane() && !s2.is_all_plane());
auto obj =
half_plane_half_plane_proj_intersection(h1, s1.line(), h2, s2.line(), k);
CGAL_assertion(! s2.is_all_plane() && ! s2.is_all_plane());
auto obj = half_plane_half_plane_proj_intersection(h1, s1.line(),
h2, s2.line(), kernel);
if(obj ==std::nullopt )
return o;
if (obj == std::nullopt) return o;
if(const Line_2* line = std::get_if<Line_2>(&(*obj)))
{
if (const auto* line = std::get_if<Line_2>(&(*obj))) {
*o++ = Intersection_curve(*line, 1);
return o;
}
if(const Ray_2* ray = std::get_if<Ray_2>(&(*obj)))
{
if (const auto* ray = std::get_if<Ray_2>(&(*obj))) {
*o++ = Intersection_curve(*ray, 1);
return o;
}
if(const Segment_2* seg = std::get_if<Segment_2>(&(*obj)))
{
if (const auto* seg = std::get_if<Segment_2>(&(*obj))) {
*o++ = Intersection_curve(*seg, 1);
return o;
}
if(const Point_2* p = std::get_if<Point_2>(&(*obj)))
{
if (const auto* p = std::get_if<Point_2>(&(*obj))) {
*o++ = *p;
return o;
}
@ -462,12 +494,10 @@ public:
}
};
//! Obtain a Construct_projected_intersections_2 functor object.
Construct_projected_intersections_2
construct_projected_intersections_2_object() const
{
return Construct_projected_intersections_2();
}
{ return Construct_projected_intersections_2(*this); }
};
} //namespace CGAL

493
Envelope_3/include/CGAL/Env_sphere_traits_3.h
View File

@ -7,10 +7,10 @@
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Michal Meyerovitch <gorgymic@post.tau.ac.il>
// Baruch Zukerman <baruchzu@post.tau.ac.il>
// Ron Wein <wein@post.tau.ac.il>
// Efi Fogel <fogel@post.tau.ac.il>
// Author(s) : Michal Meyerovitch <gorgymic@post.tau.ac.il>
// Baruch Zukerman <baruchzu@post.tau.ac.il>
// Ron Wein <wein@post.tau.ac.il>
// Efi Fogel <efifogel@gmail.com>
#ifndef CGAL_ENV_SPHERE_TRAITS_3_H
#define CGAL_ENV_SPHERE_TRAITS_3_H
@ -32,96 +32,115 @@ namespace CGAL {
template <typename ConicTraits_2>
class Env_sphere_traits_3 : public ConicTraits_2 {
public:
typedef ConicTraits_2 Traits_2;
typedef Env_sphere_traits_3<Traits_2> Self;
using Traits_2 = ConicTraits_2;
typedef typename Traits_2::Point_2 Point_2;
typedef typename Traits_2::Curve_2 Curve_2;
typedef typename Traits_2::X_monotone_curve_2 X_monotone_curve_2;
typedef typename Traits_2::Multiplicity Multiplicity;
using Point_2 = typename Traits_2::Point_2;
using Curve_2 = typename Traits_2::Curve_2;
using X_monotone_curve_2 = typename Traits_2::X_monotone_curve_2;
using Multiplicity = typename Traits_2::Multiplicity;
typedef typename Traits_2::Rat_kernel Rat_kernel;
typedef typename Traits_2::Alg_kernel Alg_kernel;
typedef typename Traits_2::Nt_traits Nt_traits;
using Rat_kernel = typename Traits_2::Rat_kernel;
using Alg_kernel = typename Traits_2::Alg_kernel;
using Nt_traits = typename Traits_2::Nt_traits;
typedef typename Rat_kernel::FT Rational;
typedef typename Rat_kernel::Point_2 Rat_point_2;
typedef typename Rat_kernel::Segment_2 Rat_segment_2;
typedef typename Rat_kernel::Line_2 Rat_line_2;
typedef typename Rat_kernel::Circle_2 Rat_circle_2;
typedef typename Rat_kernel::Point_3 Rat_point_3;
using Rational = typename Rat_kernel::FT;
using Rat_point_2 = typename Rat_kernel::Point_2;
using Rat_segment_2 = typename Rat_kernel::Segment_2;
using Rat_line_2 = typename Rat_kernel::Line_2;
using Rat_circle_2 = typename Rat_kernel::Circle_2;
using Rat_point_3 = typename Rat_kernel::Point_3;
typedef typename Alg_kernel::FT Algebraic;
typedef typename Alg_kernel::Point_2 Alg_point_2;
typedef typename Alg_kernel::Circle_2 Alg_circle_2;
using Algebraic = typename Alg_kernel::FT;
using Alg_point_2 = typename Alg_kernel::Point_2;
using Alg_circle_2 = typename Alg_kernel::Circle_2;
typedef typename Rat_kernel::Sphere_3 Surface_3;
using Surface_3 = typename Rat_kernel::Sphere_3;
// here we refer to the lower part of the sphere only
typedef Surface_3 Xy_monotone_surface_3;
using Xy_monotone_surface_3 = Surface_3;
protected:
typedef std::pair<X_monotone_curve_2, Multiplicity> Intersection_curve;
using Intersection_curve = std::pair<X_monotone_curve_2, Multiplicity>;
public:
/*! Subdivide a given surface into \f$xy\f$-monotone parts.
*/
class Make_xy_monotone_3 {
protected:
const Self& parent;
using Traits_3 = Env_sphere_traits_3<Traits_2>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Make_xy_monotone_3(const Traits_3& traits) : m_traits(traits) {}
friend class Env_sphere_traits_3<Traits_2>;
public:
Make_xy_monotone_3(const Self * p) : parent(*p) {}
// create xy-monotone surfaces from a general surface
// return a past-the-end iterator
template <class OutputIterator>
OutputIterator operator()(const Surface_3& s,
bool is_lower,
OutputIterator o) const
{
template <typename OutputIterator>
OutputIterator operator()(const Surface_3& s, bool is_lower,
OutputIterator o) const {
// our half sphere is of same type as our full sphere since we always
// need only the lower/upper part of each sphere
parent.m_is_lower = is_lower;
m_traits.m_is_lower = is_lower;
*o++ = s;
return o;
}
};
/*! Get a Make_xy_monotone_3 functor object. */
Make_xy_monotone_3
make_xy_monotone_3_object() const { return Make_xy_monotone_3(this); }
/*! Obtain a Make_xy_monotone_3 functor object. */
Make_xy_monotone_3 make_xy_monotone_3_object() const
{ return Make_xy_monotone_3(*this); }
/*! Compute all planar \f$x\f$-monotone curves and possibly isolated planar
* points that form the projection of the boundary of the given
* \f$xy\f$-monotone surface s onto the \f$xy\f$-plane.
*/
class Construct_projected_boundary_2 {
protected:
const Self& parent;
using Traits_3 = Env_sphere_traits_3<Traits_2>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Construct_projected_boundary_2(const Traits_3& traits) : m_traits(traits) {}
friend class Env_sphere_traits_3<Traits_2>;
public:
Construct_projected_boundary_2(const Self* p) : parent(*p) {}
// insert into the OutputIterator all the (2d) curves of the boundary of
// the vertical projection of the surface on the xy-plane
// the OutputIterator value type is X_monotone_curve_2
template <class OutputIterator>
template <typename OutputIterator>
OutputIterator
operator()(const Xy_monotone_surface_3& s, OutputIterator o) const {
const auto gt_2 = parent.geometry_traits_2();
const Traits_2& gt2 = m_traits;
// the projected boundary in a circle, with a projected center,
// and same radius
Rat_point_2 proj_center = parent.project(s.center());
Rat_point_2 proj_center = m_traits.project(s.center());
Rat_circle_2 circ(proj_center, s.squared_radius());
Curve_2 curve = gt_2->construct_curve_2_object()(circ);
typedef std::variant<X_monotone_curve_2, Point_2> Variant;
Curve_2 curve = gt2.construct_curve_2_object()(circ);
using Variant = std::variant<X_monotone_curve_2, Point_2>;
Variant objs[2];
CGAL_assertion_code(Variant* p = )
parent.make_x_monotone_2_object()(curve, objs);
m_traits.make_x_monotone_2_object()(curve, objs);
CGAL_assertion(p == objs + 2);
const X_monotone_curve_2* cv1 = std::get_if<X_monotone_curve_2>(&(objs[0]));
const X_monotone_curve_2* cv2 = std::get_if<X_monotone_curve_2>(&(objs[1]));
const auto* cv1 = std::get_if<X_monotone_curve_2>(&(objs[0]));
const auto* cv2 = std::get_if<X_monotone_curve_2>(&(objs[1]));
CGAL_assertion(cv1!=nullptr);
CGAL_assertion(cv2!=nullptr);
CGAL_assertion(cv1 != nullptr);
CGAL_assertion(cv2 != nullptr);
if (cv1->is_lower()) {
CGAL_assertion(cv2->is_upper());
@ -138,39 +157,52 @@ public:
}
};
/*! Get a Construct_projected_boundary_2 functor object. */
/*! Obtain a Construct_projected_boundary_2 functor object. */
Construct_projected_boundary_2
construct_projected_boundary_2_object() const
{ return Construct_projected_boundary_2(this); }
{ return Construct_projected_boundary_2(*this); }
/*! compute the projection of the intersections of the \f$xy\f$-monotone
* surfaces onto the \f$xy\f$-plane,
*/
class Construct_projected_intersections_2 {
protected:
const Self& parent;
using Traits_3 = Env_sphere_traits_3<Traits_2>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Construct_projected_intersections_2(const Traits_3& traits) :
m_traits(traits)
{}
friend class Env_sphere_traits_3<Traits_2>;
public:
Construct_projected_intersections_2(const Self * p) : parent(*p) {}
// insert into OutputIterator all the (2d) projections on the xy plane of
// the intersection objects between the 2 surfaces
// the data type of OutputIterator is Object
template <typename OutputIterator>
OutputIterator
operator()(const Xy_monotone_surface_3& s1, const Xy_monotone_surface_3& s2,
OutputIterator o) const {
const auto gt_2 = parent.geometry_traits_2();
auto ctr_cv = gt_2->construct_curve_2_object();
auto nt_traits = gt_2->nt_traits();
OutputIterator operator()(const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2,
OutputIterator o) const {
const Traits_2& gt2 = m_traits;
auto ctr_cv = gt2.construct_curve_2_object();
auto nt_traits = gt2.nt_traits();
Rat_point_3 p1 = s1.center();
Rat_point_3 p2 = s2.center();
const Rational a1 = p1.x();
const Rational b1 = p1.y();
const Rational c1 = p1.z();
const Rational a2 = p2.x();
const Rational b2 = p2.y();
const Rational c2 = p2.z();
const Rational sqr_r1 = s1.squared_radius();
const Rational sqr_r2 = s2.squared_radius();
const Rat_point_3& p1 = s1.center();
const Rat_point_3& p2 = s2.center();
const Rational& a1 = p1.x();
const Rational& b1 = p1.y();
const Rational& c1 = p1.z();
const Rational& a2 = p2.x();
const Rational& b2 = p2.y();
const Rational& c2 = p2.z();
const Rational& sqr_r1 = s1.squared_radius();
const Rational& sqr_r2 = s2.squared_radius();
// // the spheres intersect iff d(p1, p2) <= (r1+r2)
// // squaring this twice, we get the condition
@ -225,7 +257,7 @@ public:
Rational B = -8*b1*sqr_a_diff;
Rational C = 4*sqr_a_diff*(sqr_a1+sqr_b1-sqr_r1) + m*m - 4*m*a1*a_diff;
Algebraic ys[2];
Algebraic ys[2];
Algebraic* ys_end = nt_traits->solve_quadratic_equation(A, B, C, ys);
std::ptrdiff_t n_ys = ys_end - ys;
@ -254,7 +286,7 @@ public:
// 2(a1-a2)x - m = 0
Curve_2 res = ctr_cv(0, 0, 0, 2*a_diff, 0, -m, COLLINEAR, end1, end2);
parent.add_curve_to_output(res, o);
m_traits.add_curve_to_output(res, o);
//*o++ = Intersection_curve(res, TRANSVERSAL);
}
else {
@ -328,7 +360,7 @@ public:
// 2(a1-a2)x + 2(b1-b2)y - m = 0
Curve_2 res =
ctr_cv(0,0,0, 2*a_diff, 2*b_diff, -m, COLLINEAR, end1, end2);
parent.add_curve_to_output(res, o);
m_traits.add_curve_to_output(res, o);
//*o++ = Intersection_curve(res, TRANSVERSAL);
}
}
@ -401,14 +433,14 @@ public:
// if the full spheres do not intersect, the equation we get has no
// real solution, so we should check it:
bool ellipse_is_point = false;
if (! parent.is_valid_conic_equation(R, S, T, U, V, W, ellipse_is_point))
if (! m_traits.is_valid_conic_equation(R, S, T, U, V, W, ellipse_is_point))
return o;
// we need only a part of the ellipse (as stated in (**)) so we
// construct the cutting line, which is:
// equation (*) <= min(c1,c2) -- for lower envelope
// equation (*) >= max(c1,c2) -- for upper envelope
Rational z_plane = (parent.m_is_lower) ?
Rational z_plane = (m_traits.m_is_lower) ?
((c1 < c2) ? c1 : c2) : ((c1 > c2) ? c1 : c2);
@ -423,7 +455,7 @@ public:
// for upper envelope, we should multiply the line equation by -1
int envelope_coef = 1;
if (! parent.m_is_lower) envelope_coef = -1;
if (! m_traits.m_is_lower) envelope_coef = -1;
Sign sign_c_diff = CGAL_NTS sign(c_diff);
Rational la = envelope_coef*2*a_diff*sign_c_diff;
@ -454,13 +486,12 @@ public:
// intersection - in which case lc >= 0
// or there is no intersection at all between the 2 half spheres -
// in which case lc < 0
if (CGAL_NTS compare(a_diff, zero) == EQUAL &&
CGAL_NTS compare(b_diff, zero) == EQUAL)
{
if ((CGAL_NTS compare(a_diff, zero) == EQUAL) &&
(CGAL_NTS compare(b_diff, zero) == EQUAL)) {
Sign sign_lc = CGAL_NTS sign(lc);
if (sign_lc != NEGATIVE) {
Curve_2 res = ctr_cv(R, S, T, U, V, W);
parent.add_curve_to_output(res, o);
m_traits.add_curve_to_output(res, o);
//*o++ = Intersection_curve(res, TRANSVERSAL);
}
return o;
@ -510,7 +541,7 @@ public:
Curve_2 inter_cv = ctr_cv(R, S, T, U, V, W);
CGAL_precondition_code(x_mid_n_y_points = )
gt_2->points_at_x(inter_cv, x_mid_point, x_mid_y_points);
gt2.points_at_x(inter_cv, x_mid_point, x_mid_y_points);
CGAL_precondition(x_mid_n_y_points > 0);
@ -566,7 +597,7 @@ public:
Curve_2 inter_cv = ctr_cv(R, S, T, U, V, W);
CGAL_precondition_code(int y_mid_n_x_points =)
gt_2->points_at_y(inter_cv, y_mid_point, y_mid_x_points);
gt2.points_at_y(inter_cv, y_mid_point, y_mid_x_points);
CGAL_precondition(y_mid_n_x_points > 0);
Algebraic x1 = y_mid_x_points[0].x(), x2 = y_mid_x_points[1].x();
@ -600,7 +631,7 @@ public:
Alg_point_2 vtan_ps[2];
CGAL_assertion_code(std::size_t n_vtan_ps =)
gt_2->vertical_tangency_points(inter_cv, vtan_ps);
gt2.vertical_tangency_points(inter_cv, vtan_ps);
CGAL_assertion(n_vtan_ps == 2);
@ -609,7 +640,7 @@ public:
Sign lval_sign = CGAL_NTS sign(lval);
if (lval_sign == POSITIVE) {
// the full ellipse is in the positive side
parent.add_curve_to_output(inter_cv, o);
m_traits.add_curve_to_output(inter_cv, o);
//*o++ = Intersection_curve(inter_cv, TRANSVERSAL);
return o;
}
@ -630,7 +661,7 @@ public:
lval_sign = CGAL_NTS sign(lval);
CGAL_assertion(lval_sign != ZERO);
if (lval_sign == POSITIVE) parent.add_curve_to_output(inter_cv, o);
if (lval_sign == POSITIVE) m_traits.add_curve_to_output(inter_cv, o);
//*o++ = Intersection_curve(inter_cv, TRANSVERSAL);
else *o++ = Point_2(source);
@ -644,14 +675,14 @@ public:
// If the mid-point forms a left-turn with the source and the target
// points, the orientation is positive (going counterclockwise).
// Otherwise, it is negative (going clockwise).
auto alg_kernel = gt_2->alg_kernel();
auto alg_kernel = gt2.alg_kernel();
auto orient_f = alg_kernel->orientation_2_object();
Orientation orient = (orient_f(source, pmid, target) == LEFT_TURN) ?
CGAL::COUNTERCLOCKWISE : CGAL::CLOCKWISE;
Curve_2 res = ctr_cv(R, S, T, U, V, W, orient, source, target);
CGAL_assertion(res.is_valid());
parent.add_curve_to_output(res, o);
m_traits.add_curve_to_output(res, o);
//*o++ = Intersection_curve(res, TRANSVERSAL);
}
@ -659,18 +690,30 @@ public:
}
};
/*! Get a Construct_projected_intersections_2 functor object. */
/*! Obtain a Construct_projected_intersections_2 functor object. */
Construct_projected_intersections_2
construct_projected_intersections_2_object() const
{ return Construct_projected_intersections_2(this); }
{ return Construct_projected_intersections_2(*this); }
/*! Determine the relative \f$z\f$-order of two given \f$xy\f$-monotone
* surfaces at the \f$xy\f$-coordinates of a given point or \f$x\f$-monotone
* curve.
*/
class Compare_z_at_xy_3 {
protected:
const Self& parent;
using Traits_3 = Env_sphere_traits_3<Traits_2>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Compare_z_at_xy_3(const Traits_3& traits) : m_traits(traits) {}
friend class Env_sphere_traits_3<Traits_2>;
public:
Compare_z_at_xy_3(const Self* p) : parent(*p) {}
// check which of the surfaces is closer to the envelope at the xy
// coordinates of point (i.e. lower if computing the lower envelope, or
// upper if computing the upper envelope)
@ -697,9 +740,8 @@ public:
const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2) const {
// we compute a middle point on cv and use the previous function
Point_2 mid = parent.construct_middle_point(cv);
Comparison_result res = parent.compare_z_at_xy_3_object()(mid, s1, s2);
return res;
Point_2 mid = m_traits.construct_middle_point(cv);
return m_traits.compare_z_at_xy_3_object()(mid, s1, s2);
}
protected:
@ -710,12 +752,12 @@ public:
const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2) const {
// find the z coordinates of surface 1 over p
Algebraic z1 = parent.compute_envelope_z_in_point(p, s1);
Algebraic z1 = m_traits.compute_envelope_z_in_point(p, s1);
// find the z coordinates of surface 2 over p
Algebraic z2 = parent.compute_envelope_z_in_point(p, s2);
Algebraic z2 = m_traits.compute_envelope_z_in_point(p, s2);
Sign res = CGAL_NTS sign(z1 - z2);
if (parent.m_is_lower) return res;
if (m_traits.m_is_lower) return res;
else return -res;
}
@ -735,19 +777,19 @@ public:
Comparison_result
compare_in_point_second_method(const Point_2& p,
const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2) const
{
Rat_point_3 p1 = s1.center();
Rat_point_3 p2 = s2.center();
const Rational a1 = p1.x();
const Rational b1 = p1.y();
const Rational c1 = p1.z();
const Rational a2 = p2.x();
const Rational b2 = p2.y();
const Rational c2 = p2.z();
const Rational sqr_r1 = s1.squared_radius();
const Rational sqr_r2 = s2.squared_radius();
const Algebraic x1 = p.x(), y1 = p.y();
const Xy_monotone_surface_3& s2) const {
const Rat_point_3& p1 = s1.center();
const Rat_point_3& p2 = s2.center();
const Rational& a1 = p1.x();
const Rational& b1 = p1.y();
const Rational& c1 = p1.z();
const Rational& a2 = p2.x();
const Rational& b2 = p2.y();
const Rational& c2 = p2.z();
const Rational& sqr_r1 = s1.squared_radius();
const Rational& sqr_r2 = s2.squared_radius();
const Algebraic& x1 = p.x();
const Algebraic& y1 = p.y();
Rational c_diff = c1 - c2;
Algebraic x_diff1 = x1 - a1, y_diff1 = y1 - b1;
@ -773,17 +815,32 @@ public:
}
};
/*! Get a Compare_z_at_xy_3 functor object. */
/*! Obtain a Compare_z_at_xy_3 functor object. */
Compare_z_at_xy_3
compare_z_at_xy_3_object() const { return Compare_z_at_xy_3(this); }
compare_z_at_xy_3_object() const { return Compare_z_at_xy_3(*this); }
/*! Determine the relative \f$z\f$-order of the two given \f$xy\f$-monotone
* surfaces immediately above their projected intersection curve (a planar
* point \f$p\f$ is above an \f$x\f$-monotone curve \f$c\f$ if it is in the
* \f$x\f$-range of \f$c\f$, and lies to its left when the curve is traversed
* from its \f$xy\f$-lexicographically smaller endpoint to its larger
* endpoint).
*/
class Compare_z_at_xy_above_3 {
protected:
const Self& parent;
using Traits_3 = Env_sphere_traits_3<Traits_2>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Compare_z_at_xy_above_3(const Traits_3& traits) : m_traits(traits) {}
friend class Env_sphere_traits_3<Traits_2>;
public:
Compare_z_at_xy_above_3(const Self* p) : parent(*p) {}
// check which of the surfaces is closer to the envelope on the points above
// the curve cv (i.e. lower if computing the lower envelope, or upper if
// computing the upper envelope)
@ -795,35 +852,46 @@ public:
operator()(const X_monotone_curve_2& cv,
const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2) const
{ return parent.compare_on_side(cv, s1, s2, false); }
{ return m_traits.compare_on_side(cv, s1, s2, false); }
};
/*! Get a Compare_z_at_xy_above_3 functor object. */
/*! Obtain a Compare_z_at_xy_above_3 functor object. */
Compare_z_at_xy_above_3
compare_z_at_xy_above_3_object() const
{ return Compare_z_at_xy_above_3(this); }
{ return Compare_z_at_xy_above_3(*this); }
/*! Determine the relative \f$z\f$-order of the two given \f$xy\f$-monotone
* surfaces immediately below their projected intersection curve (a planar
* point \f$p\f$ is below an \f$x\f$-monotone curve \f$c\f$ if it is in the
* \f$x\f$-range of \f$c\f$, and lies to its left when the curve is traversed
* from its \f$xy\f$-lexicographically smaller endpoint to its larger
* endpoint).
*/
class Compare_z_at_xy_below_3 {
protected:
const Self& parent;
using Traits_3 = Env_sphere_traits_3<Traits_2>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Compare_z_at_xy_below_3(const Traits_3& traits) : m_traits(traits) {}
friend class Env_sphere_traits_3<Traits_2>;
public:
Compare_z_at_xy_below_3(const Self* p) : parent(*p) {}
Comparison_result
operator()(const X_monotone_curve_2& cv,
const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2) const
{
Comparison_result res = parent.compare_on_side(cv, s1, s2, true);
return res;
}
Comparison_result operator()(const X_monotone_curve_2& cv,
const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2) const
{ return m_traits.compare_on_side(cv, s1, s2, true); }
};
/*! Get a Compare_z_at_xy_below_3 functor object. */
/*! Obtain a Compare_z_at_xy_below_3 functor object. */
Compare_z_at_xy_below_3
compare_z_at_xy_below_3_object() const
{ return Compare_z_at_xy_below_3(this); }
{ return Compare_z_at_xy_below_3(*this); }
/***************************************************************************/
@ -832,31 +900,39 @@ public:
// checks if point is in the xy-range of surf
class Is_defined_over {
protected:
const Self& parent;
using Traits_3 = Env_sphere_traits_3<Traits_2>;
//! The traits (in case it has state).
const Traits_3& m_traits;
/*! Constructor
* \param traits the traits
*/
Is_defined_over(const Traits_3& traits) : m_traits(traits) {}
friend class Env_sphere_traits_3<Traits_2>;
public:
Is_defined_over(const Self* p) : parent(*p) {}
// checks if point is in the xy-range of surf
bool operator()(const Point_2& p, const Xy_monotone_surface_3& s) const {
// project the surface on the plane
Rat_point_2 proj_center = parent.project(s.center());
Rat_point_2 proj_center = m_traits.project(s.center());
Rat_circle_2 boundary(proj_center, s.squared_radius());
const auto gt_2 = parent.geometry_traits_2();
auto nt_traits = gt_2->nt_traits();
const Traits_2& gt2 = m_traits;
auto nt_traits = gt2.nt_traits();
Alg_point_2 aproj_center(proj_center.x(), proj_center.y());
Alg_circle_2 aboundary(aproj_center, nt_traits->convert(s.squared_radius()));
// check if the projected point is inside the projected boundary
auto alg_kernel = gt_2->alg_kernel();
auto alg_kernel = gt2.alg_kernel();
return (! alg_kernel->has_on_unbounded_side_2_object()(aboundary, p));
}
};
/*! Get a Is_defined_over functor object. */
/*! Obtain a Is_defined_over functor object. */
Is_defined_over is_defined_over_object() const
{ return Is_defined_over(this); }
{ return Is_defined_over(*this); }
/***************************************************************************/
@ -867,8 +943,7 @@ public:
Comparison_result compare_on_side(const X_monotone_curve_2& cv,
const Xy_monotone_surface_3& s1,
const Xy_monotone_surface_3& s2,
bool compare_on_right) const
{
bool compare_on_right) const {
// cv(x,y) : r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0
// let p be the leftmost endpoint of cv, p=(x0, y0)
// the tangence of cv at p is a line. on the infinitesimal region
@ -881,7 +956,8 @@ public:
// since we assume that such point represents better what is going
// on all internal curve points
Point_2 cv_point = construct_middle_point(cv);
Algebraic x0 = cv_point.x(), y0 = cv_point.y();
const Algebraic& x0 = cv_point.x();
const Algebraic& y0 = cv_point.y();
// d(cv)/dx : 2r*x + 2s*y*dy/dx + t*y + t*x*dy/dx +u + v*dy/dx = 0
// in point p=(x0,y0) we get
@ -895,8 +971,12 @@ public:
// n != 0: y - y0 = y'(x-x0) ==> -y'x + y + (y'x0 - y0) = 0
// and in general we have:
// -m*x + n*y + (m*x0 -n*y0) = 0 (with integer coordinates)
const Rational r = cv.r(), s = cv.s(), t = cv.t(),
u = cv.u(), v = cv.v(), w = cv.w();
const Rational& r = cv.r();
const Rational& s = cv.s();
const Rational& t = cv.t();
const Rational& u = cv.u();
const Rational& v = cv.v();
const Rational& w = cv.w();
Algebraic m = -1 * (2*r*x0 + t*y0 + u);
Algebraic n = 2*s*y0 + t*x0 + v;
// line coefficients: A3, B3, C3
@ -931,14 +1011,14 @@ public:
// Di = -(x0-ai)x0 - (y0-bi)y0 - (z0-ci)z0
//
// and we solve the problem as for triangles
Rat_point_3 p1 = s1.center();
Rat_point_3 p2 = s2.center();
const Rational a1 = p1.x();
const Rational b1 = p1.y();
const Rational c1 = p1.z();
const Rational a2 = p2.x();
const Rational b2 = p2.y();
const Rational c2 = p2.z();
const Rat_point_3& p1 = s1.center();
const Rat_point_3& p2 = s2.center();
const Rational& a1 = p1.x();
const Rational& b1 = p1.y();
const Rational& c1 = p1.z();
const Rational& a2 = p2.x();
const Rational& b2 = p2.y();
const Rational& c2 = p2.z();
Algebraic A1 = x0 - a1, B1 = y0 - b1, C1 = z0 - c1;
Algebraic A2 = x0 - a2, B2 = y0 - b2, C2 = z0 - c2;
if (C1 != 0 && C2 != 0) {
@ -971,6 +1051,7 @@ public:
return EQUAL;
}
//
Rat_point_2 project(const Rat_point_3& p) const
{ return Rat_point_2(p.x(), p.y()); }
@ -979,15 +1060,15 @@ public:
// precondition: s is defined at p
Algebraic compute_envelope_z_in_point(const Point_2& p,
const Xy_monotone_surface_3& s) const {
Algebraic res;
// the point coordinates
const Algebraic x1 = p.x(), y1 = p.y();
// the surface equations
Rat_point_3 center = s.center();
const Rational a = center.x(), b = center.y(), c = center.z();
const Rational sqr_r = s.squared_radius();
const Rat_point_3& center = s.center();
const Rational& a = center.x();
const Rational& b = center.y();
const Rational& c = center.z();
const Rational& sqr_r = s.squared_radius();
// we substitute x1 and y1 in the equation of s
// (x-a)^2 + (y-b)^2 + (z-c)^2 = r^2
@ -999,11 +1080,12 @@ public:
Algebraic B = -2*c;
Algebraic C = x_diff*x_diff + y_diff*y_diff + c*c - sqr_r;
Algebraic zs[2];
Algebraic zs[2];
Algebraic* zs_end;
std::ptrdiff_t n_zs;
auto nt_traits = m_geometry_traits_2->nt_traits();
const Traits_2& gt2 = *this;
auto nt_traits = gt2.nt_traits();
zs_end = nt_traits->solve_quadratic_equation(A, B, C, zs);
n_zs = zs_end - zs;
@ -1013,6 +1095,7 @@ public:
if (n_zs == 1) return zs[0];
CGAL_assertion(n_zs == 2);
Algebraic res;
Comparison_result comp = CGAL_NTS compare(zs[0], zs[1]);
if (m_is_lower) res = ((comp == SMALLER) ? zs[0] : zs[1]);
else res = ((comp == LARGER) ? zs[0] : zs[1]);
@ -1021,8 +1104,10 @@ public:
// construct the point in the middle of cv
Point_2 construct_middle_point(const X_monotone_curve_2& cv) const {
const Traits_2& gt2 = *this;
// get the x-value of the middle point
auto alg_kernel = m_geometry_traits_2->alg_kernel();
auto alg_kernel = gt2.alg_kernel();
Alg_point_2 mid_x =
alg_kernel->construct_midpoint_2_object()(cv.source(), cv.target());
@ -1030,13 +1115,14 @@ public:
// maybe we want it there?
// if (cv.is_segment()) return mid_x;
if (cv.is_vertical()) return Point_2(mid_x);
return Point_2(m_geometry_traits_2->point_at_x(cv, mid_x));
return Point_2(gt2.point_at_x(cv, mid_x));
}
// for the test
Point_2 construct_middle_point(const Point_2& p1, const Point_2& p2) const {
auto alg_kernel = m_geometry_traits_2->alg_kernel();
const Traits_2& gt2 = *this;
auto alg_kernel = gt2.alg_kernel();
return Point_2(alg_kernel->construct_midpoint_2_object()(p1, p2));
}
@ -1044,7 +1130,7 @@ public:
// r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0
// has real solutions
// is_point is set to true if the equation represents just one point
template <class NT>
template <typename NT>
bool is_valid_conic_equation(const NT& r, const NT& s, const NT& t,
const NT& u, const NT& v, const NT& w,
bool& is_point) const {
@ -1103,30 +1189,31 @@ public:
// for the test:
Point_2 vertical_ray_shoot_2(const Point_2& pt, const X_monotone_curve_2& cv)
const
{
const {
const Traits_2& gt2 = *this;
if (cv.is_vertical()) {
auto alg_kernel = m_geometry_traits_2->alg_kernel();
if (! alg_kernel->less_y_2_object()(cv.left(), pt))
return cv.left();
else {
CGAL_assertion(alg_kernel->less_y_2_object()(cv.right(), pt));
return cv.right();
}
auto alg_kernel = gt2.alg_kernel();
if (! alg_kernel->less_y_2_object()(cv.left(), pt)) return cv.left();
CGAL_assertion(alg_kernel->less_y_2_object()(cv.right(), pt));
return cv.right();
}
else return m_geometry_traits_2->point_at_x(cv, pt);
return gt2.point_at_x(cv, pt);
}
//
template <typename OutputIterator>
OutputIterator add_curve_to_output(const Curve_2& c, OutputIterator oi) const {
std::variant<X_monotone_curve_2, Point_2> objs[2];
std::variant<X_monotone_curve_2, Point_2>* p_obj = this->make_x_monotone_2_object()(c, objs);
for(std::variant<X_monotone_curve_2, Point_2>* o = objs; o != p_obj; ++o) {
if(const X_monotone_curve_2* cv = std::get_if<X_monotone_curve_2>(o))
std::variant<X_monotone_curve_2, Point_2>* p_obj =
this->make_x_monotone_2_object()(c, objs);
for (std::variant<X_monotone_curve_2, Point_2>* o = objs; o != p_obj; ++o) {
if (const auto* cv = std::get_if<X_monotone_curve_2>(o))
*oi++ = Intersection_curve(*cv, 1);
else {
const Point_2* pt = std::get_if<Point_2>(o);
const auto* pt = std::get_if<Point_2>(o);
CGAL_assertion(pt!=nullptr);
*oi++ = *pt;
}
@ -1134,55 +1221,39 @@ public:
return oi;
}
typedef std::shared_ptr<Traits_2> Shared_geometry_traits_2;
/*! Default constructor. */
Env_sphere_traits_3() :
m_is_lower(true),
m_geometry_traits_2(new Traits_2)
{}
/*! Constructor from a conic 2D geometry traits. */
Env_sphere_traits_3(Shared_geometry_traits_2 geometry_traits_2) :
m_is_lower(true),
m_geometry_traits_2(geometry_traits_2)
{}
/*! Obtain the undelying conic 2D geometry traits.
*/
const Shared_geometry_traits_2 geometry_traits_2() const
{ return m_geometry_traits_2; }
Env_sphere_traits_3() : m_is_lower(true) {}
protected:
mutable bool m_is_lower;
private:
//! The conic geometry-traits.
const Shared_geometry_traits_2 m_geometry_traits_2;
};
/*!
* Compare two spheres: first compare their center points in an
* xyz-lexicographic order, then by their radii.
/*! Compare two spheres: first compare their center points in an
* xyz-lexicographic order, then by their squared radii.
*
* \todo This should be deprecated and instead the traits should provide this
* functionality.
*/
template <typename Kernel>
bool operator<(const CGAL::Sphere_3<Kernel>& a, const CGAL::Sphere_3<Kernel>& b)
{
Kernel k;
Comparison_result res = k.compare_xyz_3_object()(a.center(), b.center());
if (res == EQUAL) res = CGAL::compare (a.squared_radius(), b.squared_radius());
if (res == EQUAL) res = CGAL::compare(a.squared_radius(), b.squared_radius());
return (res == SMALLER);
}
/*!
* Compare two spheres for equality.
/*! Compare two spheres for equality.
*
* \todo This should be deprecated and instead the traits should provide this
* functionality.
*/
template <typename Kernel>
bool operator== (const typename Kernel::Sphere_3& a,
const typename Kernel::Sphere_3& b) {
bool operator==(const typename Kernel::Sphere_3& a,
const typename Kernel::Sphere_3& b) {
Kernel k;
if (! k.equal_3_object() (a.center(), b.center())) return (false);
return (CGAL::compare (a.squared_radius(), b.squared_radius()) == EQUAL);
if (! k.equal_3_object()(a.center(), b.center())) return (false);
return (CGAL::compare(a.squared_radius(), b.squared_radius()) == EQUAL);
}
} //namespace CGAL

86
Envelope_3/include/CGAL/Env_surface_data_traits_3.h
View File

@ -8,7 +8,9 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Ron Wein <wein@post.tau.ac.il>
// Author(s) : Ron Wein <wein@post.tau.ac.il>
// Efi Fogel <efifogel@gmail.com>
#ifndef CGAL_ENV_SURFACE_DATA_TRAITS_3_H
#define CGAL_ENV_SURFACE_DATA_TRAITS_3_H
@ -19,9 +21,10 @@
* Definition of the env_surface_data_traits_3<> class template.
*/
#include <CGAL/Arr_geometry_traits/Curve_data_aux.h>
#include <list>
#include <CGAL/Arr_geometry_traits/Curve_data_aux.h>
namespace CGAL {
/*! \class
@ -31,95 +34,80 @@ namespace CGAL {
* It can attach data objects to Surface_3 and to Xy_monotone_surface_3 objects
* (possibly of two different types).
*/
template <class Traits_, class XyMonotoneSurfaceData_,
class SurfaceData_ = XyMonotoneSurfaceData_,
class Convert_ = _Default_convert_func<SurfaceData_,
XyMonotoneSurfaceData_> >
class Env_surface_data_traits_3 : public Traits_
{
template <typename Traits_, typename XyMonotoneSurfaceData,
typename SurfaceData = XyMonotoneSurfaceData,
typename Convert_ = _Default_convert_func<SurfaceData,
XyMonotoneSurfaceData>>
class Env_surface_data_traits_3 : public Traits_ {
public:
using Base_traits_3 = Traits_;
using Xy_monotone_surface_data = XyMonotoneSurfaceData;
using Surface_data = SurfaceData;
using Convert = Convert_;
typedef Traits_ Base_traits_3;
typedef XyMonotoneSurfaceData_ Xy_monotone_surface_data;
typedef SurfaceData_ Surface_data;
typedef Convert_ Convert;
typedef typename Base_traits_3::Surface_3 Base_surface_3;
typedef typename Base_traits_3::Xy_monotone_surface_3
Base_xy_monotone_surface_3;
using Base_surface_3 = typename Base_traits_3::Surface_3;
using Base_xy_monotone_surface_3 =
typename Base_traits_3::Xy_monotone_surface_3;
// Representation of a surface with an additional data field:
typedef _Curve_data_ex<Base_surface_3, Surface_data> Surface_3;
using Surface_3 = _Curve_data_ex<Base_surface_3, Surface_data>;
// Representation of an xy-monotone surface with an additional data field:
typedef _Curve_data_ex<Base_xy_monotone_surface_3,
Xy_monotone_surface_data> Xy_monotone_surface_3;
using Xy_monotone_surface_3 =
_Curve_data_ex<Base_xy_monotone_surface_3, Xy_monotone_surface_data>;
public:
/// \name Construction.
//@{
/*! Default constructor. */
Env_surface_data_traits_3 ()
{}
Env_surface_data_traits_3() {}
/*! Constructor from a base-traits class. */
Env_surface_data_traits_3 (const Base_traits_3 & traits) :
Base_traits_3 (traits)
Env_surface_data_traits_3(const Base_traits_3& traits) :
Base_traits_3(traits)
{}
//@}
/// \name Overridden functors.
//@{
class Make_xy_monotone_3
{
class Make_xy_monotone_3 {
private:
const Base_traits_3 * base;
const Base_traits_3* base;
public:
/*! Constructor. */
Make_xy_monotone_3 (const Base_traits_3 * _base) : base (_base)
{}
Make_xy_monotone_3(const Base_traits_3* _base) : base (_base) {}
/*!
* Subdivide the given surface into xy-monotone surfaces and insert them
/*! Subdivide the given surface into xy-monotone surfaces and insert them
* to the given output iterator.
* \param S The surface.
* \param oi The output iterator,
* whose value-type is Xy_monotone_surface_2.
* \return The past-the-end iterator.
*/
template<class OutputIterator>
OutputIterator operator() (const Surface_3& S, bool is_lower,
OutputIterator oi) const
{
template <typename OutputIterator>
OutputIterator operator()(const Surface_3& S, bool is_lower,
OutputIterator oi) const {
// Make the original surface xy-monotone.
std::list<Base_xy_monotone_surface_3> xy_surfaces;
typename std::list<Base_xy_monotone_surface_3>::iterator xys_it;
std::list<Base_xy_monotone_surface_3> xy_surfaces;
typename std::list<Base_xy_monotone_surface_3>::iterator xys_it;
base->make_xy_monotone_3_object()
(S, is_lower, std::back_inserter (xy_surfaces));
base->make_xy_monotone_3_object()(S, is_lower,
std::back_inserter(xy_surfaces));
// Attach the data to each of the resulting xy-monotone surfaces.
for (xys_it = xy_surfaces.begin(); xys_it != xy_surfaces.end(); ++xys_it)
{
*oi = Xy_monotone_surface_3 (*xys_it,
Convert() (S.data()));
++oi;
}
*oi++ = Xy_monotone_surface_3(*xys_it, Convert() (S.data()));
return (oi);
return oi;
}
};
/*! Get a Make_xy_monotone_3 functor object. */
Make_xy_monotone_3 make_xy_monotone_3_object() const
{
return Make_xy_monotone_3 (this);
}
{ return Make_xy_monotone_3(this); }
//@}

1217
Envelope_3/include/CGAL/Env_triangle_traits_3.h
View File

File diff suppressed because it is too large Load Diff