From b452d5a7e1878e737be68d27a3399d38b8faa879 Mon Sep 17 00:00:00 2001 From: Jocelyn MEYRON Date: Tue, 11 Nov 2014 15:51:20 +0100 Subject: [PATCH] use boost::optional for the computation (or not) of an interior point --- .../dual/halfspace_intersection_3.h | 4 +- ...fspace_intersection_with_constructions_3.h | 6 +- .../doc/Convex_hull_3/Convex_hull_3.txt | 6 +- .../halfspace_intersection_3.cpp | 2 +- .../dual/halfspace_intersection_3.h | 91 ++++++++------ ...fspace_intersection_with_constructions_3.h | 111 +++++++++--------- 6 files changed, 120 insertions(+), 100 deletions(-) diff --git a/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h b/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h index 30bb627265d..5b5b73c9fda 100644 --- a/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h +++ b/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h @@ -4,7 +4,7 @@ namespace CGAL { \ingroup PkgConvexHull3Functions \brief computes robustly the intersection of the halfspaces defined by the planes contained in the range [`begin`, `end`) without constructing the dual points. The result is stored in the polyhedron `P`. -`origin` is a point strictly inside the polyhedron. +If `origin` is given then it must be a point strictly inside the polyhedron. If an interior point is not given then it is computed using a linear program and thus is slower than the first approach. This version does not compute the dual points, but instead it uses a traits class for handling predicates for dual points without computing them. \attention Halfspaces are considered as lower halfspaces that is to say if the plane's equation is \f$ a\, x +b\, y +c\, z + d = 0 \f$ then the corresponding halfspace is defined by \f$ a\, x +b\, y +c\, z + d \le 0 \f$ . @@ -22,6 +22,6 @@ This version does not compute the dual points, but instead it uses a traits clas template void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, Polyhedron &P, - typename Polyhedron::Vertex::Point_3 const& origin); + Point_3 origin); } /* namespace CGAL */ diff --git a/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h b/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h index a3f67485dba..f01aa781f7a 100644 --- a/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h +++ b/Convex_hull_3/doc/Convex_hull_3/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h @@ -1,10 +1,10 @@ namespace CGAL { /*! -\ingroup PkgConvexHull3Functions +\Ingroup PkgConvexHull3Functions \brief computes the intersection of the halfspaces defined by the planes contained in the range [`begin`, `end`). The result is stored in the polyhedron `P`. -`origin` is a point strictly inside the polyhedron. +If `origin` is given then it must be a point strictly inside the polyhedron. If an interior point is not given then it is computed using a linear program and thus is slower than the first approach. This version constructs explicitly the dual points using the convex hull algorithm parametrized with the given traits class. \attention Halfspaces are considered as lower halfspaces that is to say if the plane's equation is \f$ a\, x +b\, y +c\, z + d = 0 \f$ then the corresponding halfspace is defined by \f$ a\, x +b\, y +c\, z + d \le 0 \f$ . @@ -24,7 +24,7 @@ template void halfspace_intersection_with_constructions_3(PlaneIterator pbegin, PlaneIterator pend, Polyhedron &P, - typename Polyhedron::Vertex::Point_3 const& origin, + Point_3 origin, const Traits & ch_traits = Default_traits); } /* namespace CGAL */ diff --git a/Convex_hull_3/doc/Convex_hull_3/Convex_hull_3.txt b/Convex_hull_3/doc/Convex_hull_3/Convex_hull_3.txt index 81a7d2c20a9..5fc8ea6684f 100644 --- a/Convex_hull_3/doc/Convex_hull_3/Convex_hull_3.txt +++ b/Convex_hull_3/doc/Convex_hull_3/Convex_hull_3.txt @@ -67,11 +67,15 @@ point set or not. This function is used in postcondition testing for \subsection Convex_hull_3HalfspaceIntersection Halfspace Intersection The functions `halfspace_intersection_3()` and `halfspace_intersection_with_constructions_3()` -uses the convex hull computation and the duality to compute the intersection +uses the convex hull algorithm and the duality to compute the intersection of a list of halfspaces. The first version does not explicitly computes the dual points: the traits class handles this issue. The second one constructs these points and hence is less robust but the computation is faster. +In order to compute, the intersection an interior point is needed. It could be +either given by the user or computed using linear programming. Notice that the +second approach is slower due to the resolution of a linear program. + The following program computes the intersection of halfspaces defined by tangent planes to a sphere: diff --git a/Convex_hull_3/examples/Convex_hull_3/halfspace_intersection_3.cpp b/Convex_hull_3/examples/Convex_hull_3/halfspace_intersection_3.cpp index 5784a3c23b7..7f86d281ec6 100644 --- a/Convex_hull_3/examples/Convex_hull_3/halfspace_intersection_3.cpp +++ b/Convex_hull_3/examples/Convex_hull_3/halfspace_intersection_3.cpp @@ -38,7 +38,7 @@ int main (void) { planes.end(), P, Point(0, 0, 0)); - + return 0; } diff --git a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h index 3383c98b32e..250dfa0d08c 100644 --- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h +++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h @@ -257,37 +257,13 @@ namespace CGAL } // namespace internal } // namespace Convex_hull_3 - // Compute the intersection of halfspaces by computing a point - // inside the polyhedron (using linear programming) - template - void halfspace_intersection_without_origin_3 (PlaneIterator begin, PlaneIterator end, - Polyhedron &P) { - // Types - typedef typename Kernel_traits::Kernel K; - typedef typename Polyhedron::Vertex::Point_3 Point_3; - - // choose exact integral type -#ifdef CGAL_USE_GMP - typedef CGAL::Gmpq ET; -#else - typedef CGAL::MP_Float ET; -#endif - // find a point inside the intersection - typedef Interior_polyhedron_3 Interior_polyhedron; - Interior_polyhedron interior; - bool res = interior.find(begin, end); - CGAL_assertion_msg(res, "halfspace_intersection_without_origin_3: problem when determing a point inside"); - Point_3 origin = interior.inside_point(); - - // compute the intersection - halfspace_intersection_3(begin, end, P, origin); - } - - // Compute the intersection of halfspaces with the origin given by the user + // Compute the intersection of halfspaces. + // If the user gives an origin then the function used it, otherwise, it is + // computed using a linear program. template void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, Polyhedron &P, - typename Polyhedron::Vertex::Point_3 const& origin) { + boost::optional const& origin) { // Checks whether the intersection if a polyhedron CGAL_assertion_msg(Convex_hull_3::internal::is_intersection_dim_3(begin, end), "halfspace_intersection_3: intersection not a polyhedron"); @@ -296,18 +272,61 @@ namespace CGAL typedef Convex_hull_3::Convex_hull_traits_dual_3 Hull_traits_dual_3; typedef Polyhedron_3 Polyhedron_dual_3; typedef Convex_hull_3::internal::Build_primal_polyhedron Builder; + typedef typename Polyhedron::Vertex::Point_3 Point_3; - Hull_traits_dual_3 dual_traits(origin); + if (origin) { + Point_3 p_origin = boost::get(origin); + Hull_traits_dual_3 dual_traits(p_origin); - Polyhedron_dual_3 dual_convex_hull; - CGAL::convex_hull_3(begin, end, dual_convex_hull, dual_traits); - Builder build_primal(dual_convex_hull, origin); - P.delegate(build_primal); + Polyhedron_dual_3 dual_convex_hull; + CGAL::convex_hull_3(begin, end, dual_convex_hull, dual_traits); + Builder build_primal(dual_convex_hull, p_origin); + P.delegate(build_primal); - // Posterior check if the origin is inside the computed polyhedron - Polyhedron Q(P); - CGAL_assertion_msg(!Convex_hull_3::internal::point_inside_convex_polyhedron(Q, origin), "halfspace_intersection_3: origin not in the polyhedron"); + // Posterior check if the origin is inside the computed polyhedron + Polyhedron Q(P); + CGAL_assertion_msg(!Convex_hull_3::internal::point_inside_convex_polyhedron(Q, p_origin), "halfspace_intersection_3: origin not in the polyhedron"); + } else { + // choose exact integral type +#ifdef CGAL_USE_GMP + typedef CGAL::Gmpq ET; +#else + typedef CGAL::MP_Float ET; +#endif + // find a point inside the intersection + typedef Interior_polyhedron_3 Interior_polyhedron; + Interior_polyhedron interior; + bool res = interior.find(begin, end); + CGAL_assertion_msg(res, "halfspace_intersection_without_origin_3: problem when determing a point inside"); + Point_3 origin = interior.inside_point(); + + Hull_traits_dual_3 dual_traits(origin); + + Polyhedron_dual_3 dual_convex_hull; + CGAL::convex_hull_3(begin, end, dual_convex_hull, dual_traits); + Builder build_primal(dual_convex_hull, origin); + P.delegate(build_primal); + + // Posterior check if the origin is inside the computed polyhedron + Polyhedron Q(P); + CGAL_assertion_msg(!Convex_hull_3::internal::point_inside_convex_polyhedron(Q, origin), "halfspace_intersection_3: origin not in the polyhedron"); + } } + + // Compute the intersection of halfspaces by computing a point inside. + template + void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, + Polyhedron &P) { + halfspace_intersection_3(begin, end , P, boost::none); + } + + // Compute the intersection of halfspaces (an interior point is given.) + template + void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, + Polyhedron &P, + typename Polyhedron::Vertex::Point_3 const& origin) { + halfspace_intersection_3(begin, end , P, boost::optional(origin)); + } } // namespace CGAL #endif // CGAL_HALFSPACE_INTERSECTION_3_H diff --git a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h index 34721b17ab1..a35736ac4fd 100644 --- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h +++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h @@ -116,14 +116,32 @@ namespace CGAL void halfspace_intersection_with_constructions_3(PlaneIterator pbegin, PlaneIterator pend, Polyhedron &P, - typename Polyhedron::Vertex::Point_3 const& origin, - const Traits & ch_traits) - { + boost::optional const& origin, + const Traits & ch_traits) { typedef typename Kernel_traits::Kernel K; typedef typename K::Point_3 Point; typedef typename K::Plane_3 Plane; typedef typename CGAL::internal::Build_dual_polyhedron Builder; + Point p_origin; + + if (origin) { + p_origin = boost::get(origin); + } else { + // choose exact integral type +#ifdef CGAL_USE_GMP + typedef CGAL::Gmpq ET; +#else + typedef CGAL::MP_Float ET; +#endif + // find a point inside the intersection + typedef Interior_polyhedron_3 Interior_polyhedron; + Interior_polyhedron interior; + bool res = interior.find(pbegin, pend); + CGAL_assertion_msg(res, "halfspace_intersection_with_constructions_3: problem when determing a point inside"); + p_origin = interior.inside_point(); + } + // construct dual points to apply the convex hull std::vector dual_points; for (PlaneIterator p = pbegin; p != pend; ++p) { @@ -131,54 +149,39 @@ namespace CGAL Plane translated_p(p->a(), p->b(), p->c(), - p->d() + origin.x() * p->a() + origin.y() * p->b() + origin.z() * p->c()); + p->d() + p_origin.x() * p->a() + p_origin.y() * p->b() + p_origin.z() * p->c()); dual_points.push_back(CGAL::ORIGIN + translated_p.orthogonal_vector () / (-translated_p.d())); } Polyhedron ch; CGAL::convex_hull_3(dual_points.begin(), dual_points.end(), ch, ch_traits); - Builder build_dual (ch, origin); + Builder build_dual (ch, p_origin); P.delegate(build_dual); } // Compute the intersection of halfspaces by constructing explicitly // the dual points with the traits class for convex_hull_3 given // as an argument. - // The point inside the polyhedron is computed using linear programming. + // An interior point is given. template - void halfspace_intersection_with_constructions_without_origin_3 (PlaneIterator begin, PlaneIterator end, - Polyhedron &P, - const Traits & ch_traits) { - // Types - typedef typename Kernel_traits::Kernel K; - typedef typename Polyhedron::Vertex::Point_3 Point_3; - - // choose exact integral type -#ifdef CGAL_USE_GMP - typedef CGAL::Gmpq ET; -#else - typedef CGAL::MP_Float ET; -#endif - // find a point inside the intersection - typedef Interior_polyhedron_3 Interior_polyhedron; - Interior_polyhedron interior; - bool res = interior.find(begin, end); - CGAL_assertion_msg(res, "halfspace_intersection_with_constructions_without_origin_3: problem when determing a point inside"); - Point_3 origin = interior.inside_point(); - - // compute the intersection - halfspace_intersection_with_constructions_3(begin, end, P, origin, ch_traits); - } - - // Compute the intersection of halfspaces by constructing explicitly - // the dual points with the default traits class for convex_hull_3 - template void halfspace_intersection_with_constructions_3(PlaneIterator pbegin, PlaneIterator pend, Polyhedron &P, - typename Polyhedron::Vertex::Point_3 const& origin) - { + typename Polyhedron::Vertex::Point_3 const& origin, + const Traits & ch_traits) { + halfspace_intersection_with_constructions_3(pbegin, pend, P, + boost::optional(origin), + ch_traits); + } + + // Compute the intersection of halfspaces by constructing explicitly + // the dual points with the default traits class for convex_hull_3. + template + void halfspace_intersection_with_constructions_3 (PlaneIterator pbegin, + PlaneIterator pend, + Polyhedron &P, + boost::optional const& origin) { typedef typename Kernel_traits::Kernel K; typedef typename K::Point_3 Point_3; typedef typename internal::Convex_hull_3::Default_traits_for_Chull_3::type Traits; @@ -188,30 +191,24 @@ namespace CGAL // Compute the intersection of halfspaces by constructing explicitly // the dual points with the default traits class for convex_hull_3. - // The point inside the polyhedron is computed using linear programming. + // An interior point is given. template - void halfspace_intersection_with_constructions_without_origin_3 (PlaneIterator begin, PlaneIterator end, - Polyhedron &P) { - // Types - typedef typename Kernel_traits::Kernel K; - typedef typename Polyhedron::Vertex::Point_3 Point_3; - typedef typename internal::Convex_hull_3::Default_traits_for_Chull_3::type Traits; + void halfspace_intersection_with_constructions_3 (PlaneIterator pbegin, + PlaneIterator pend, + Polyhedron &P, + typename Polyhedron::Vertex::Point_3 const& origin) { + halfspace_intersection_with_constructions_3(pbegin, pend, P, + boost::optional(origin)); + } - // choose exact integral type -#ifdef CGAL_USE_GMP - typedef CGAL::Gmpq ET; -#else - typedef CGAL::MP_Float ET; -#endif - // find a point inside the intersection - typedef Interior_polyhedron_3 Interior_polyhedron; - Interior_polyhedron interior; - bool res = interior.find(begin, end); - CGAL_assertion_msg(res, "halfspace_intersection_with_constructions_without_origin_3: problem when determing an point inside"); - Point_3 origin = interior.inside_point(); - - // compute the intersection - halfspace_intersection_with_constructions_3(begin, end, P, origin, Traits()); + // Compute the intersection of halfspaces by constructing explicitly + // the dual points with the default traits class for convex_hull_3. + // An interior point is not given. + template + void halfspace_intersection_with_constructions_3 (PlaneIterator pbegin, + PlaneIterator pend, + Polyhedron &P) { + halfspace_intersection_with_constructions_3(pbegin, pend, P, boost::none); } } // namespace CGAL