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 4135165f892..4d88721cc77 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 @@ -9,7 +9,7 @@ If `origin` is given then it must be a point strictly inside the polyhedron. If This version does not construct the dual points explicitely but uses a special traits class for the function `CGAL::convex_hull_3()` to handle predicates on dual points without constructing 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$ . -\attention The value type of `PlaneIterator` and the point type of `origin` must come from the same \cgal %Kernel. +\attention The point type of `origin` and the point type of the vertices of `Polyhedron` must come from the same \cgal %Kernel. \pre if provided, `origin` is inside the intersection of halfspaces defined by the range `[begin, end)`. \pre The computed intersection must be a bounded convex polyhedron. @@ -23,6 +23,6 @@ This version does not construct the dual points explicitely but uses a special t template void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, Polyhedron &P, - boost::optional origin = boost::none); + boost::optional::value_type>::Kernel::Point_3> > origin = boost::none); } /* 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 a7e645ac64d..bd0c5ff4535 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 @@ -61,7 +61,7 @@ The following program reads points from an input file and computes their convex hull. We assume that the points are not all identical and not all collinear, thus we directly use a polyhedron as output. In the example you see that the convex hull function can -write in any model of the concept `MutableFaceListGraph`. +write in any model of the concept `MutableFaceGraph`. \cgalExample{Convex_hull_3/quickhull_3.cpp} @@ -93,8 +93,8 @@ second approach is slower due to the resolution of a linear program. \subsubsection Convex_hull_3HalfspaceInntersectionExample Example -The following program computes the intersection of halfspaces defined by -tangent planes to a sphere: +The following example computes the intersection of halfspaces defined by +tangent planes to a sphere. \cgalExample{Convex_hull_3/halfspace_intersection_3.cpp} @@ -125,6 +125,12 @@ not all of them are vertices of the hull. \subsection Convex_hull_3Example_2 Example +The following example shows how to compute a convex hull with a triangulation. +The vertices incident to the infinite vertex are on the convex hull. + +The function `star_to_face_graph()` can be used to obtain a polyhedral surface +that is model of the concept `MutableFaceGraph`, e.g. `Polyhedron_3` and `Surface_mesh`. + \cgalExample{Convex_hull_3/dynamic_hull_3.cpp} \section Convex_hull_3Performance Performance diff --git a/Convex_hull_3/examples/Convex_hull_3/dynamic_hull_3.cpp b/Convex_hull_3/examples/Convex_hull_3/dynamic_hull_3.cpp index c08c0a683de..478e88f6c02 100644 --- a/Convex_hull_3/examples/Convex_hull_3/dynamic_hull_3.cpp +++ b/Convex_hull_3/examples/Convex_hull_3/dynamic_hull_3.cpp @@ -1,7 +1,7 @@ #include #include #include -#include +#include #include #include @@ -11,7 +11,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef K::Point_3 Point_3; typedef CGAL::Delaunay_triangulation_3 Delaunay; typedef Delaunay::Vertex_handle Vertex_handle; -typedef CGAL::Polyhedron_3 Polyhedron_3; +typedef CGAL::Surface_mesh Surface_mesh; int main() { @@ -39,7 +39,7 @@ int main() //copy the convex hull of points into a polyhedron and use it //to get the number of points on the convex hull - Polyhedron_3 chull; + Surface_mesh chull; CGAL::star_to_face_graph(T, T.infinite_vertex(), chull); std::cout << "After removal of 25 points, there are " 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 100f5d68ea6..c66820f5dac 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 @@ -1,13 +1,15 @@ #include #include #include +#include #include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef K::Plane_3 Plane; typedef K::Point_3 Point; -typedef CGAL::Polyhedron_3 Polyhedron_3; +typedef CGAL::Surface_mesh Surface_mesh; + // compute the tangent plane of a point template @@ -31,15 +33,16 @@ int main (void) { } // define polyhedron to hold the intersection - Polyhedron_3 P; + Surface_mesh chull; // compute the intersection // if no point inside the intersection is provided, one // will be automatically found using linear programming CGAL::halfspace_intersection_3(planes.begin(), planes.end(), - P, - boost::make_optional(Point(0, 0, 0)) ); + chull ); + + std::cout << "The convex hull contains " << num_vertices(chull) << " vertices" << std::endl; 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 49b9af65206..490b335db94 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 @@ -23,18 +23,18 @@ #define CGAL_HALFSPACE_INTERSECTION_3_H #include -#include #include #include #include #include #include - +#include // For interior_polyhedron_3 #include #include - - #include +#include +#include +#include namespace CGAL { @@ -44,130 +44,123 @@ namespace CGAL { // Build the primal polyhedron associated to a dual polyhedron // We also need the `origin` which represents a point inside the primal polyhedron - template - class Build_primal_polyhedron : - public CGAL::Modifier_base { - typedef typename Polyhedron::HalfedgeDS HDS; - const Polyhedron_dual & _dual; + template + void + build_primal_polyhedron(Polyhedron& primal, + const Polyhedron_dual & _dual, + const Point_3& origin) + { + typedef typename Kernel_traits::Kernel Kernel; + typedef typename K::RT RT; - // Origin - typedef typename K::Point_3 Primal_point_3; - Primal_point_3 origin; + // Typedefs for dual + typedef typename Polyhedron_dual::Facet Facet; + typedef typename Polyhedron_dual::Facet_const_handle + Facet_const_handle; + typedef typename Polyhedron_dual::Facet_const_iterator + Facet_const_iterator; + typedef typename Polyhedron_dual::Vertex_const_iterator + Vertex_const_iterator; - public: - Build_primal_polyhedron (const Polyhedron_dual & dual, - Primal_point_3 o = - Primal_point_3(CGAL::ORIGIN)) : _dual (dual), origin(o) - {} + // typedefs for primal + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - void operator () (HDS &hds) - { - typedef typename K::RT RT; - typedef typename K::Point_3 Point_3; + // Typedefs for intersection + typedef typename K::Plane_3 Plane_3; + typedef typename K::Line_3 Line_3; + typedef boost::optional< boost::variant< Point_3, + Line_3, + Plane_3 > > result_inter; - // Typedefs for dual - typedef typename Polyhedron_dual::Facet Facet; - typedef typename Polyhedron_dual::Facet_const_handle - Facet_const_handle; - typedef typename Polyhedron_dual::Facet_const_iterator - Facet_const_iterator; - typedef typename Polyhedron_dual::Vertex_const_iterator - Vertex_const_iterator; + + boost::unordered_map primal_vertices; + size_t n = 0; - // Typedefs for primal - typename CGAL::Polyhedron_incremental_builder_3 B(hds, true); + // First, computing the primal vertices + for (Facet_const_iterator it = _dual.facets_begin(); + it != _dual.facets_end(); ++it, ++n) { + typename Facet::Halfedge_const_handle h = it->halfedge(); + // Build the dual plane corresponding to the current facet + Plane_3 p1 = h->vertex()->point(); + Plane_3 p2 = h->next()->vertex()->point(); + Plane_3 p3 = h->next()->next()->vertex()->point(); - // Typedefs for intersection - typedef typename K::Plane_3 Plane_3; - typedef typename K::Line_3 Line_3; - typedef boost::optional< boost::variant< Point_3, - Line_3, - Plane_3 > > result_inter; + RT dp1 = p1.d() + origin.x() * p1.a() + + origin.y() * p1.b() + origin.z() * p1.c(); + RT dp2 = p2.d() + origin.x() * p2.a() + + origin.y() * p2.b() + origin.z() * p2.c(); + RT dp3 = p3.d() + origin.x() * p3.a() + + origin.y() * p3.b() + origin.z() * p3.c(); - B.begin_surface(_dual.size_of_facets(), - _dual.size_of_vertices(), - _dual.size_of_halfedges()); + Plane_3 pp1(p1.a(), p1.b(), p1.c(), dp1); + Plane_3 pp2(p2.a(), p2.b(), p2.c(), dp2); + Plane_3 pp3(p3.a(), p3.b(), p3.c(), dp3); - std::map primal_vertices; - size_t n = 0; + // Compute the intersection + result_inter result = CGAL::intersection(pp1, pp2, pp3); + CGAL_assertion_msg(bool(result), + "halfspace_intersection_3: no intersection"); + CGAL_assertion_msg(boost::get(& *result) != NULL, + "halfspace_intersection_3: intersection is not a point"); - // First, computing the primal vertices - for (Facet_const_iterator it = _dual.facets_begin(); - it != _dual.facets_end(); ++it, ++n) { - typename Facet::Halfedge_const_handle h = it->halfedge(); - // Build the dual plane corresponding to the current facet - Plane_3 p1 = h->vertex()->point(); - Plane_3 p2 = h->next()->vertex()->point(); - Plane_3 p3 = h->next()->next()->vertex()->point(); + const Point_3* pp = boost::get(& *result); - RT dp1 = p1.d() + origin.x() * p1.a() - + origin.y() * p1.b() + origin.z() * p1.c(); - RT dp2 = p2.d() + origin.x() * p2.a() - + origin.y() * p2.b() + origin.z() * p2.c(); - RT dp3 = p3.d() + origin.x() * p3.a() - + origin.y() * p3.b() + origin.z() * p3.c(); + // Primal vertex associated to the current dual plane + Point_3 ppp(origin.x() + pp->x(), + origin.y() + pp->y(), + origin.z() + pp->z()); - Plane_3 pp1(p1.a(), p1.b(), p1.c(), dp1); - Plane_3 pp2(p2.a(), p2.b(), p2.c(), dp2); - Plane_3 pp3(p3.a(), p3.b(), p3.c(), dp3); + + primal_vertices[it] = add_vertex(ppp, primal); + } - // Compute the intersection - result_inter result = CGAL::intersection(pp1, pp2, pp3); - CGAL_assertion_msg(bool(result), - "halfspace_intersection_3: no intersection"); - CGAL_assertion_msg(boost::get(& *result) != NULL, - "halfspace_intersection_3: intersection is not a point"); + // Then, add facets to the primal polyhedron + // To do this, for each dual vertex, we circulate around this vertex + // and we add an edge between each facet we encounter - const Point_3* pp = boost::get(& *result); + for (Vertex_const_iterator it = _dual.vertices_begin(); + it != _dual.vertices_end(); ++it) { + std::vector vertices; + typename Polyhedron_dual::Halfedge_around_vertex_const_circulator + h0 = it->vertex_begin(), hf = h0; + + //B.begin_facet(); + do { + vertices.push_back(primal_vertices[hf->facet()]); + } while (--hf != h0); + Euler::add_face(vertices,primal); + //B.end_facet(); + + } - // Primal vertex associated to the current dual plane - Point_3 ppp(origin.x() + pp->x(), - origin.y() + pp->y(), - origin.z() + pp->z()); - B.add_vertex(ppp); - primal_vertices[it] = n; - } - - // Then, add facets to the primal polyhedron - // To do this, for each dual vertex, we circulate around this vertex - // and we add an edge between each facet we encounter - for (Vertex_const_iterator it = _dual.vertices_begin(); - it != _dual.vertices_end(); ++it) { - typename Polyhedron_dual::Halfedge_around_vertex_const_circulator - h0 = it->vertex_begin(), hf = h0; - - B.begin_facet(); - do { - B.add_vertex_to_facet(primal_vertices[hf->facet()]); - } while (--hf != h0); - B.end_facet(); - } - - B.end_surface(); - } - }; + } // Test if a point is inside a convex polyhedron - template + template bool point_inside_convex_polyhedron (const Polyhedron &P, - typename Polyhedron::Traits::Point_3 const& p) { - // Compute the equations of the facets of the polyhedron - typedef typename Polyhedron::Facet_const_iterator Facet_iterator; - for(Facet_iterator fit=P.facets_begin(), fit_end=P.facets_end(); - fit!=fit_end; ++fit) + Point const& p) { + // Compute the equations of the facets of the polyhedron + typedef boost::graph_traits::face_descriptor face_descriptor; + typedef boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typename boost::property_map::type vpmap = get(CGAL::vertex_point, P); + + BOOST_FOREACH(face_descriptor fd, faces(P)) { - typename Polyhedron::Halfedge_const_handle h = fit->halfedge(); - typename Polyhedron::Traits::Point_3 const& p1=h->vertex()->point(); - typename Polyhedron::Traits::Point_3 const& p2=h->next()->vertex()->point(); - h=h->next()->next(); - typename Polyhedron::Traits::Point_3 p3 = h->vertex()->point(); - while( h!=fit->halfedge() && collinear(p1,p2,p3) ) + halfedge_descriptor h = halfedge(fd,P), done(h); + Point_3 const& p1 = get(vpmap, target(h,P)); + h = next(h,p); + Point_3 const& p2 = get(vpmap, target(h,P)); + h = next(h,p); + Point_3 const& p3 = get(vpmap, target(h,P)); + + while( h!=done && collinear(p1,p2,p3) ) { - h=h->next(); - p3 = h->vertex()->point(); + h = next(h,p); + p3 = get(vpmap, target(h,P)); } - if (h==fit->halfedge()) continue; //degenerate facet, skip it + if (h==done) continue; //degenerate facet, skip it if ( orientation (p1, p2, p3, p) != CGAL::NEGATIVE ) return false; } @@ -246,15 +239,14 @@ namespace CGAL template void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, Polyhedron &P, - boost::optional origin = boost::none) { + boost::optional::value_type>::Kernel::Point_3> origin = boost::none) { // Checks whether the intersection is a polyhedron CGAL_assertion_msg(Convex_hull_3::internal::is_intersection_dim_3(begin, end), "halfspace_intersection_3: intersection not a polyhedron"); // Types - typedef typename Kernel_traits::Kernel K; + typedef typename Kernel_traits::value_type>::Kernel K; 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; // if a point inside is not provided find one using linear programming if (!origin) { @@ -278,8 +270,7 @@ namespace CGAL 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); + Convex_hull_3::internal::build_primal_polyhedron(P, dual_convex_hull, *origin); // Posterior check if the origin is inside the computed polyhedron // The check is done only if the number type is not float or double because in that @@ -296,7 +287,7 @@ namespace CGAL template void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end, Polyhedron &P, - typename Polyhedron::Vertex::Point_3 const& origin) { + typename Kernel_traits::value_type>::Kernel::Point_3 const& origin) { halfspace_intersection_3(begin, end , P, boost::make_optional(origin)); } #endif diff --git a/Convex_hull_3/include/CGAL/convex_hull_3.h b/Convex_hull_3/include/CGAL/convex_hull_3.h index a8118bac224..b3f7e3f8fa1 100644 --- a/Convex_hull_3/include/CGAL/convex_hull_3.h +++ b/Convex_hull_3/include/CGAL/convex_hull_3.h @@ -38,7 +38,6 @@ #include #include #include -#include #include #include #include @@ -51,6 +50,8 @@ #include #include +#include + #ifndef CGAL_CH_NO_POSTCONDITIONS #include #endif // CGAL_CH_NO_POSTCONDITIONS @@ -327,7 +328,7 @@ find_visible_set(TDS_2& tds, const typename Traits::Point_3& point, typename TDS_2::Face_handle start, std::list& visible, - std::map& outside, + boost::unordered_map& outside, const Traits& traits) { typedef typename Traits::Plane_3 Plane_3; @@ -473,7 +474,7 @@ ch_quickhull_3_scan(TDS_2& tds, typedef typename Traits::Point_3 Point_3; typedef std::list Outside_set; typedef typename std::list::iterator Outside_set_iterator; - typedef std::map Border_edges; + typedef boost::unordered_map Border_edges; std::list visible_set; typename std::list::iterator vis_set_it;