diff --git a/BGL/include/CGAL/boost/graph/named_params_helper.h b/BGL/include/CGAL/boost/graph/named_params_helper.h index 1d68f2a6f87..1e0efd6faa6 100644 --- a/BGL/include/CGAL/boost/graph/named_params_helper.h +++ b/BGL/include/CGAL/boost/graph/named_params_helper.h @@ -110,22 +110,75 @@ public: typedef typename boost::property_traits::value_type type; }; + +template +struct GetVertexPointMap_impl +{ + typedef VPM_from_NP type; + typedef VPM_from_NP const_type; + + template + static const_type + get_const_map(const NamedParameters& np, const PolygonMesh&) + { + return parameters::get_parameter(np, internal_np::vertex_point); + } + + template + static type + get_map(const NamedParameters& np, PolygonMesh&) + { + return parameters::get_parameter(np, internal_np::vertex_point); + } +}; + +template +struct GetVertexPointMap_impl +{ + typedef typename property_map_selector::const_type const_type; + typedef typename property_map_selector::type type; + + template + static const_type + get_const_map(const NamedParameters& /* np */, const PolygonMesh& pm) + { + return get_const_property_map(boost::vertex_point, pm); + } + + template + static type + get_map(const NamedParameters& /* np */, PolygonMesh& pm) + { + return get_property_map(boost::vertex_point, pm); + } +}; + template class GetVertexPointMap { - typedef typename property_map_selector::const_type - DefaultVPMap_const; - typedef typename property_map_selector::type - DefaultVPMap; + typedef typename internal_np::Lookup_named_param_def::type VPM_from_NP; + + typedef GetVertexPointMap_impl Impl; public: - typedef typename internal_np::Lookup_named_param_def::type type; - typedef typename internal_np::Lookup_named_param_def::type const_type; + typedef typename Impl::type type; + typedef typename Impl::const_type const_type; + + static const_type + get_const_map(const NamedParameters& np, const PolygonMesh& pm) + { + return Impl::get_const_map(np, pm); + } + + static type + get_map(const NamedParameters& np, PolygonMesh& pm) + { + return Impl::get_map(np, pm); + } }; template @@ -138,10 +191,15 @@ public: typedef typename CGAL::Kernel_traits::Kernel Kernel; }; -template -class GetGeomTraits + +template +struct GetGeomTraits_impl +{ + typedef GT type; +}; + +template +struct GetGeomTraits_impl { typedef typename CGAL::graph_has_property::type Has_internal_pmap; @@ -154,12 +212,20 @@ class GetGeomTraits typedef typename boost::mpl::if_c::value, typename GetK::Kernel, - Fake_GT>::type DefaultKernel; + Fake_GT>::type type; +}; -public: +template +struct GetGeomTraits +{ typedef typename internal_np::Lookup_named_param_def::type type; + internal_np::Param_not_found>::type GT_from_NP; + typedef typename GetGeomTraits_impl::type type; }; // Define the following structs: @@ -278,6 +344,21 @@ public: typedef typename CGAL::Identity_property_map const_type; }; +template +struct GetPolygonSoupGeomTraits +{ + typedef typename internal_np::Lookup_named_param_def < + internal_np::geom_traits_t, + NamedParameters, + typename CGAL::Kernel_traits< + typename boost::property_traits< + typename GetPointMap::type + >::value_type + >::type + > ::type type; +}; + + template struct Point_set_processing_3_np_helper { diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index f52e8537932..5d3bbdd0cf8 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -45,6 +45,8 @@ CGAL tetrahedral Delaunay refinement algorithm. - Added a named parameter to `CGAL::Polygon_mesh_processing::smooth_shape()` to disable scaling to compensate volume loss. +- Added the functions `CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect()` and `CGAL::Polygon_mesh_processing::triangle_soup_self_intersections()` to identify and report self-intersections in a triangle soup, similarly to existing functions on triangle meshes. + ### [3D Simplicial Mesh Data Structure](https://doc.cgal.org/5.6/Manual/packages.html#PkgSMDS3) (new package) - This new package wraps all the existing code that deals with a `MeshComplex_3InTriangulation_3` to describe 3D simplicial meshes, and makes the data structure independent from the tetrahedral mesh generation package. diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index 07a2d21a57b..c347d4a7ae3 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -152,6 +152,8 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. \cgalCRPSection{Intersection Functions} - `CGAL::Polygon_mesh_processing::does_self_intersect()` - `CGAL::Polygon_mesh_processing::self_intersections()` +- `CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect()` +- `CGAL::Polygon_mesh_processing::triangle_soup_self_intersections()` - \link PMP_intersection_grp `CGAL::Polygon_mesh_processing::do_intersect()` \endlink - `CGAL::Polygon_mesh_processing::intersecting_meshes()` diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h index cc07ed99f65..3524b46c775 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h @@ -43,6 +43,7 @@ #endif #include +#include #include #include @@ -54,6 +55,127 @@ namespace CGAL { namespace Polygon_mesh_processing { namespace internal { +template +struct Triangle_mesh_and_triangle_soup_wrapper +{ + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; // private + + static void get_face_vertices(face_descriptor fd, std::array& vh, const TM& tm) + { + CGAL_assertion(boost::graph_traits::null_face() != fd); + halfedge_descriptor h = halfedge(fd, tm); + vh[0]=source(h, tm); + vh[1]=target(h, tm); + vh[2]=target(next(h, tm), tm); + } + + static bool faces_have_a_shared_edge(face_descriptor f, face_descriptor g, std::array& vh, const TM& tm) + { + CGAL_assertion(boost::graph_traits::null_face() != f); + CGAL_assertion(boost::graph_traits::null_face() != g); + halfedge_descriptor h=halfedge(f, tm); + for(unsigned int i=0; i<3; ++i) + { + halfedge_descriptor opp_h = opposite(h, tm); + if(face(opp_h, tm) == g) + { + vh[0]=source(h, tm); + vh[1]=target(h, tm); + vh[2]=target(next(h, tm), tm); + vh[3]=target(next(opp_h, tm), tm); + return true; + } + h = next(h, tm); + } + return false; + } + + static bool is_pure_triangle(const TM& tm) + { + return is_triangle_mesh(tm); + } +}; + +template +struct Triangle_mesh_and_triangle_soup_wrapper< + std::pair> +{ + typedef std::size_t face_descriptor; + typedef std::size_t vertex_descriptor; + + typedef std::pair Soup; + + static void get_face_vertices(face_descriptor fd, std::array& vh, const Soup& soup) + { + const auto& face = soup.second[fd]; + vh[0]=face[0]; + vh[1]=face[1]; + vh[2]=face[2]; + } + + static bool faces_have_a_shared_edge(face_descriptor fd, face_descriptor gd, std::array& vh, const Soup& soup) + { + const auto& f = soup.second[fd]; + const auto& g = soup.second[gd]; + + for(unsigned int i=0; i<2; ++i) // no need to check f[2] if neither f[0] nor f[1] are shared + { + for(unsigned int j=0; j<3; ++j) + { + if (f[i]==g[j]) + { + vh[0]=f[i]; + vh[1]=f[i+1]; + vh[2]=f[(i+2)%3]; + + if (vh[1]==g[(j+1)%3]) + { + vh[3]=g[(j+2)%3]; + return true; + } + if (vh[1]==g[(j+2)%3]) + { + vh[3]=g[(j+1)%3]; + return true; + } + + if (i==0) + { + vh[1]=f[i]; + vh[2]=f[(i+1)%3]; + vh[0]=f[(i+2)%3]; + if (vh[0]==g[(j+1)%3]) + { + vh[3]=g[(j+2)%3]; + return true; + } + if (vh[0]==g[(j+2)%3]) + { + vh[3]=g[(j+1)%3]; + return true; + } + } + + return false; + } + } + } + + return false; + } + + static bool is_pure_triangle(const Soup& soup) + { + for (const typename std::iterator_traits::value_type& t : soup.second) + if (t.size()!=3) + return false; + return true; + } +}; + template struct Throw_at_count_reached_functor { @@ -82,59 +204,49 @@ struct Throw_at_count_reached_functor { // Checks for 'real' intersections, i.e. not simply a shared vertex or edge template -bool do_faces_intersect(typename boost::graph_traits::halfedge_descriptor h, - typename boost::graph_traits::halfedge_descriptor g, +bool do_faces_intersect(typename Triangle_mesh_and_triangle_soup_wrapper::face_descriptor fh, + typename Triangle_mesh_and_triangle_soup_wrapper::face_descriptor fg, const TM& tmesh, const VPM vpmap, const typename GT::Construct_segment_3& construct_segment, const typename GT::Construct_triangle_3& construct_triangle, const typename GT::Do_intersect_3& do_intersect) { - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef Triangle_mesh_and_triangle_soup_wrapper Wrapper; + typedef typename Wrapper::vertex_descriptor vertex_descriptor; - typedef typename GT::Segment_3 Segment; - typedef typename GT::Triangle_3 Triangle; + typedef typename GT::Segment_3 Segment; + typedef typename GT::Triangle_3 Triangle; - CGAL_assertion(!is_border(h, tmesh)); - CGAL_assertion(!is_border(g, tmesh)); - vertex_descriptor hv[3], gv[3]; - hv[0] = target(h, tmesh); - hv[1] = target(next(h, tmesh), tmesh); - hv[2] = source(h, tmesh); - - gv[0] = target(g, tmesh); - gv[1] = target(next(g, tmesh), tmesh); - gv[2] = source(g, tmesh); + std::array hv, gv; + Wrapper::get_face_vertices(fh, hv, tmesh); + Wrapper::get_face_vertices(fg, gv, tmesh); // check for shared edge - for(unsigned int i=0; i<3; ++i) + std::array verts; + if (Wrapper::faces_have_a_shared_edge(fh, fg, verts, tmesh)) { - halfedge_descriptor opp_h = opposite(h, tmesh); - if(face(opp_h, tmesh) == face(g, tmesh)) - { - // there is an intersection if the four points are coplanar and the triangles overlap - if(CGAL::coplanar(get(vpmap, hv[i]), - get(vpmap, hv[(i+1)%3]), - get(vpmap, hv[(i+2)%3]), - get(vpmap, target(next(opp_h, tmesh), tmesh))) && - CGAL::coplanar_orientation(get(vpmap, hv[(i+2)%3]), - get(vpmap, hv[i]), - get(vpmap, hv[(i+1)%3]), - get(vpmap, target(next(opp_h, tmesh), tmesh))) - == CGAL::POSITIVE) - { - return true; - } - else - { - // there is a shared edge but no intersection - return false; - } - } + if (verts[2]==verts[3]) return false; // only for a soup of triangles - h = next(h, tmesh); + // there is an intersection if the four points are coplanar and the triangles overlap + if(CGAL::coplanar(get(vpmap, verts[0]), + get(vpmap, verts[1]), + get(vpmap, verts[2]), + get(vpmap, verts[3])) && + CGAL::coplanar_orientation(get(vpmap, verts[0]), + get(vpmap, verts[1]), + get(vpmap, verts[2]), + get(vpmap, verts[3])) + == CGAL::POSITIVE) + { + return true; + } + else + { + // there is a shared edge but no intersection + return false; + } } // check for shared vertex --> maybe intersection, maybe not @@ -188,8 +300,6 @@ template struct Strict_intersect_faces // "strict" as in "not sharing a subface" { - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - mutable OutputIterator m_iterator; const TM& m_tmesh; const VPM m_vpmap; @@ -209,10 +319,7 @@ struct Strict_intersect_faces // "strict" as in "not sharing a subface" void operator()(const Box* b, const Box* c) const { - const halfedge_descriptor h = halfedge(b->info(), m_tmesh); - const halfedge_descriptor g = halfedge(c->info(), m_tmesh); - - if(do_faces_intersect(h, g, m_tmesh, m_vpmap, m_construct_segment, m_construct_triangle, m_do_intersect)) + if(do_faces_intersect(b->info(), c->info(), m_tmesh, m_vpmap, m_construct_segment, m_construct_triangle, m_do_intersect)) *m_iterator++ = std::make_pair(b->info(), c->info()); } }; @@ -229,15 +336,17 @@ self_intersections_impl(const FaceRange& face_range, const bool throw_on_SI, const NamedParameters& np) { - CGAL_precondition(CGAL::is_triangle_mesh(tmesh)); + typedef TriangleMesh TM; + typedef Triangle_mesh_and_triangle_soup_wrapper Wrapper; + + CGAL_precondition(Wrapper::is_pure_triangle(tmesh)); using CGAL::parameters::choose_parameter; using CGAL::parameters::get_parameter; using CGAL::parameters::is_default_parameter; - typedef TriangleMesh TM; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename Wrapper::face_descriptor face_descriptor; + typedef typename Wrapper::vertex_descriptor vertex_descriptor; typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy; typedef CGAL::Box_intersection_d::Box_with_info_d Box; @@ -245,9 +354,9 @@ self_intersections_impl(const FaceRange& face_range, typedef typename GetGeomTraits::type GT; GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); - typedef typename GetVertexPointMap::const_type VPM; - VPM vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), - get_const_property_map(boost::vertex_point, tmesh)); + typedef GetVertexPointMap VPM_helper; + typedef typename VPM_helper::const_type VPM; + VPM vpmap = VPM_helper::get_const_map(np, tmesh); const bool do_limit = !(is_default_parameter::value); const unsigned int maximum_number = choose_parameter(get_parameter(np, internal_np::maximum_number), 0); @@ -268,11 +377,13 @@ self_intersections_impl(const FaceRange& face_range, // This loop is very cheap, so there is hardly anything to gain from parallelizing it for(face_descriptor f : face_range) { - halfedge_descriptor h = halfedge(f, tmesh); + std::array vh; + Wrapper::get_face_vertices(f, vh, tmesh); + typename boost::property_traits::reference - p = get(vpmap, target(h,tmesh)), - q = get(vpmap, target(next(h, tmesh), tmesh)), - r = get(vpmap, target(prev(h, tmesh), tmesh)); + p = get(vpmap, vh[0]), + q = get(vpmap, vh[1]), + r = get(vpmap, vh[2]); // tiny fixme: if f is degenerate, we might still have a real intersection between f // and another face f', but right now we are not creating a box for f and thus not returning those @@ -342,7 +453,7 @@ self_intersections_impl(const FaceRange& face_range, std::atomic atomic_counter(counter); Throw_functor throwing_count_functor(atomic_counter, maximum_number, std::back_inserter(face_pairs)); Throwing_after_count_output_iterator count_filter(throwing_count_functor); - Filtered_intersecting_faces_filter limited_callback(tmesh, vpmap, gt, count_filter); + Filtered_intersecting_faces_filter limited_callback(tmesh, vpmap, gt, count_filter); CGAL::box_self_intersection_d(box_ptr.begin(), box_ptr.end(), limited_callback, cutoff); } catch(const CGAL::internal::Throw_at_output_exception&) @@ -408,9 +519,9 @@ self_intersections_impl(const FaceRange& face_range, * * collects intersections between a subset of faces of a triangulated surface mesh. * Two faces are said to intersect if the corresponding triangles intersect - * and the intersection is not an edge nor a vertex incident to both faces. + * and the intersection is neither an edge nor a vertex incident to both faces. * - * This function depends on the package \ref PkgBoxIntersectionD + * This function depends on the package \ref PkgBoxIntersectionD. * * @pre `CGAL::is_triangle_mesh(tmesh)` * @@ -475,9 +586,9 @@ self_intersections(const FaceRange& face_range, * * collects intersections between all the faces of a triangulated surface mesh. * Two faces are said to intersect if the corresponding triangles intersect - * and the intersection is not an edge nor a vertex incident to both faces. + * and the intersection is neither an edge nor a vertex incident to both faces. * - * This function depends on the package \ref PkgBoxIntersectionD + * This function depends on the package \ref PkgBoxIntersectionD. * * @pre `CGAL::is_triangle_mesh(tmesh)` * @@ -490,9 +601,9 @@ self_intersections(const FaceRange& face_range, * * @param tmesh the triangulated surface mesh to be checked * @param out output iterator to be filled with all pairs of non-adjacent faces that intersect. - In case `tmesh` contains some degenerate faces, for each degenerate face `f` a pair `(f,f)` - will be put in `out` before any other self intersection between non-degenerate faces. - These are the only pairs where degenerate faces will be reported. + * In case `tmesh` contains some degenerate faces, for each degenerate face `f` a pair `(f,f)` + * will be put in `out` before any other self intersection between non-degenerate faces.
+ * Note that these are the only pairs where degenerate faces will be reported. * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below * * \cgalNamedParamsBegin @@ -541,7 +652,7 @@ self_intersections(const TriangleMesh& tmesh, * * \brief tests if a set of faces of a triangulated surface mesh self-intersects. * - * This function depends on the package \ref PkgBoxIntersectionD + * This function depends on the package \ref PkgBoxIntersectionD. * * @pre `CGAL::is_triangle_mesh(tmesh)` * @@ -584,8 +695,6 @@ bool does_self_intersect(const FaceRange& face_range, const TriangleMesh& tmesh, const NamedParameters& np = parameters::default_values()) { - CGAL_precondition(CGAL::is_triangle_mesh(tmesh)); - try { CGAL::Emptyset_iterator unused_out; @@ -614,7 +723,7 @@ bool does_self_intersect(const FaceRange& face_range, * * \brief tests if a triangulated surface mesh self-intersects. * - * This function depends on the package \ref PkgBoxIntersectionD + * This function depends on the package \ref PkgBoxIntersectionD. * * @pre `CGAL::is_triangle_mesh(tmesh)` * @@ -656,6 +765,219 @@ bool does_self_intersect(const TriangleMesh& tmesh, return does_self_intersect(faces(tmesh), tmesh, np); } + +#ifndef DOXYGEN_RUNNING +template +struct Property_map_for_soup +{ + typedef std::size_t key_type; + typedef typename boost::property_traits::value_type value_type; + //typedef typename boost::property_traits::category category; + typedef boost::readable_property_map_tag category; + typedef typename boost::property_traits::reference reference; + + const PointRange& points; + VPM vpm; + + Property_map_for_soup(const PointRange& points, VPM vpm) + : points(points) + , vpm(vpm) + {} + + inline friend + reference get(const Property_map_for_soup& map, key_type k) + { + return get(map.vpm, map.points[k]); + } +}; +#endif + +/** + * \ingroup PMP_intersection_grp + * + * collects intersections between all the triangles in a triangle soup. + * + * Two triangles of the soup are said to intersect if the corresponding geometric triangles intersect + * and the intersection is neither an edge nor a vertex of both triangles + * (with the same point ids, ignoring the orientation for an edge). + * + * This function depends on the package \ref PkgBoxIntersectionD. + * + * @tparam ConcurrencyTag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`. + * @tparam PointRange a model of the concept `RandomAccessContainer` + * whose value type is the point type + * @tparam TriangleRange a model of the concept `RandomAccessContainer` whose + * value type is a model of the concept `RandomAccessContainer` whose value type is `std::size_t` + * @tparam TriangleIdPairOutputIterator a model of `OutputIterator` holding objects of type + * `std::pair` + * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" + * + * @param points points of the soup of triangles + * @param triangles each element in the range describes a triangle using the indices of the points in `points` + * @param out output iterator to be filled with all pairs of ids of triangles intersecting (the id of a triangle is its position in `triangles`). + * In case the triangle soup contains some degenerate triangles, for each degenerate triangle `t` with id `i` a pair `(i,i)` + * will be put in `out` before any other self intersection between non-degenerate triangles.
+ * Note that these are the only pairs where degenerate triangles will be reported. + * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamNBegin{maximum_number} + * \cgalParamDescription{the maximum number of self intersections that will be detected and returned by the function.} + * \cgalParamType{unsigned int} + * \cgalParamDefault{No limit.} + * \cgalParamExtra{In parallel mode, the number of returned self-intersections is at least `maximum_number` + * (and not exactly that number) as no strong synchronization is put on threads for performance reasons.} + * \cgalParamNEnd + * + * \cgalParamNBegin{point_map} + * \cgalParamDescription{a property map associating points to the elements of the range `points`} + * \cgalParamType{a model of `ReadablePropertyMap` whose value type is a point type from a \cgal `Kernel`.} + * \cgalParamDefault{`CGAL::Identity_property_map`} + * \cgalParamNEnd + * + * \cgalParamNBegin{geom_traits} + * \cgalParamDescription{an instance of a geometric traits class} + * \cgalParamType{a class model of `PMPSelfIntersectionTraits`} + * \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} + * \cgalParamExtra{The geometric traits class must be compatible with the point type of the point map.} + * \cgalParamNEnd + * \cgalNamedParamsEnd + * + * @return `out` + * + * @sa `does_triangle_soup_self_intersect()` + * @sa `self_intersections()` + * @sa `does_self_intersect()` + */ +template +TriangleIdPairOutputIterator +triangle_soup_self_intersections(const PointRange& points, + const TriangleRange& triangles, + TriangleIdPairOutputIterator out, + const CGAL_NP_CLASS& np = parameters::default_values()) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + using parameters::is_default_parameter; + + typedef typename CGAL::GetPointMap::const_type Point_map_base; + Point_map_base pm_base = choose_parameter(get_parameter(np, internal_np::point_map)); + typedef Property_map_for_soup Point_map; + typedef typename GetPolygonSoupGeomTraits::type GT; + GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); + + const bool do_limit = !(is_default_parameter::value); + if (do_limit) + { + return self_intersections(boost::irange(0, triangles.size()), + std::make_pair(std::cref(points), std::cref(triangles)), + out, + parameters::vertex_point_map(Point_map(points,pm_base)). + geom_traits(gt). + maximum_number(choose_parameter(get_parameter(np, internal_np::maximum_number), 0))); + } + + return self_intersections(boost::irange(0, triangles.size()), + std::make_pair(std::cref(points), std::cref(triangles)), + out, + parameters::vertex_point_map(Point_map(points,pm_base)). + geom_traits(gt)); +} + +/** + * \ingroup PMP_intersection_grp + * + * \brief tests if a triangle soup self-intersects. + * + * A triangle soup self-intersects if at least two triangles of the soup intersect. + * Two triangles of the soup are said to intersect if the corresponding geometric triangles intersect + * and the intersection is neither an edge nor a vertex of both triangles + * (with the same point ids, ignoring the orientation for an edge). + * + * This function depends on the package \ref PkgBoxIntersectionD. + * + * @tparam ConcurrencyTag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`. + * @tparam PointRange a model of the concept `RandomAccessContainer` + * whose value type is the point type + * @tparam TriangleRange a model of the concept `RandomAccessContainer` whose + * value type is a model of the concept `RandomAccessContainer` whose value type is `std::size_t` + * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" + * + * @param points points of the soup of triangles + * @param triangles each element in the range describes a triangle using the indices of the points in `points` + * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamNBegin{point_map} + * \cgalParamDescription{a property map associating points to the elements of the range `points`} + * \cgalParamType{a model of `ReadablePropertyMap` whose value type is a point type from a \cgal `Kernel`.} + * \cgalParamDefault{`CGAL::Identity_property_map`} + * \cgalParamNEnd + * + * \cgalParamNBegin{geom_traits} + * \cgalParamDescription{an instance of a geometric traits class} + * \cgalParamType{a class model of `PMPSelfIntersectionTraits`} + * \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} + * \cgalParamExtra{The geometric traits class must be compatible with the point type of the point map.} + * \cgalParamNEnd + * \cgalNamedParamsEnd + * + * @return `true` if the triangle soup self-intersects, and `false` otherwise. + * + * @sa `triangle_soup_self_intersections()` + * @sa `self_intersections()` + * @sa `does_self_intersect()` + */ +template +bool does_triangle_soup_self_intersect(const PointRange& points, + const TriangleRange& triangles, + const CGAL_NP_CLASS& np = parameters::default_values()) +{ + try + { + using parameters::choose_parameter; + using parameters::get_parameter; + + CGAL::Emptyset_iterator unused_out; + typedef typename CGAL::GetPointMap::const_type Point_map_base; + Point_map_base pm_base = choose_parameter(get_parameter(np, internal_np::point_map)); + typedef Property_map_for_soup Point_map; + typedef typename GetPolygonSoupGeomTraits::type GT; + GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); + + internal::self_intersections_impl(boost::irange(0, triangles.size()), + std::make_pair(std::cref(points), std::cref(triangles)), + unused_out, true /*throw*/, + parameters::vertex_point_map(Point_map(points,pm_base)) + .geom_traits(gt)); + } + catch (const CGAL::internal::Throw_at_output_exception&) + { + return true; + } + #if defined(CGAL_LINKED_WITH_TBB) && TBB_USE_CAPTURED_EXCEPTION + catch (const tbb::captured_exception& e) + { + const char* ti1 = e.name(); + const char* ti2 = typeid(const CGAL::internal::Throw_at_output_exception&).name(); + const std::string tn1(ti1); + const std::string tn2(ti2); + if (tn1 == tn2) return true; + else throw; + } + #endif + return false; +} + }// namespace Polygon_mesh_processing }// namespace CGAL diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 9848f9fdf9d..422c72854c2 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -21,6 +21,7 @@ create_single_source_cgal_program("point_inside_surface_mesh_test.cpp") create_single_source_cgal_program("polygon_mesh_slicer_test.cpp") create_single_source_cgal_program("self_intersection_polyhedron_test.cpp") create_single_source_cgal_program("self_intersection_surface_mesh_test.cpp") +create_single_source_cgal_program("self_intersection_triangle_soup_test.cpp") create_single_source_cgal_program("pmp_do_intersect_test.cpp") create_single_source_cgal_program("test_is_polygon_soup_a_polygon_mesh.cpp") create_single_source_cgal_program("test_stitching.cpp") diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/self_intersection_triangle_soup_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/self_intersection_triangle_soup_test.cpp new file mode 100644 index 00000000000..e067e13c129 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/self_intersection_triangle_soup_test.cpp @@ -0,0 +1,234 @@ +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel EPICK; +typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK; + +namespace PMP = ::CGAL::Polygon_mesh_processing; +namespace CP = ::CGAL::parameters; + +template +int test_self_intersections(const std::string filename, + const bool expected) +{ + std::vector points; + std::vector< std::array > triangles; + + bool read_ok = CGAL::IO::read_polygon_soup(filename, points, triangles); + + if ( !read_ok ) { + std::cerr << "Error: cannot read file: " << filename << std::endl; + return 1; + } + + std::cout << "Reading file: " << filename << std::endl; + + CGAL::Timer timer; + timer.start(); + + std::vector > intersected_tris; + + PMP::triangle_soup_self_intersections(points, triangles, std::back_inserter(intersected_tris)); + bool intersecting_1 = !intersected_tris.empty(); + + std::cout << "self_intersections test took " << timer.time() << " sec." << std::endl; + std::cout << intersected_tris.size() << " pairs of triangles are intersecting." << std::endl; + + timer.reset(); + + bool intersecting_2 = + PMP::does_triangle_soup_self_intersect(points, triangles); + + std::cout << "does_self_intersect test took " << timer.time() << " sec." << std::endl; + std::cout << (intersecting_2 ? "There is a self-intersection." : + "There are no self-intersections.") << std::endl; + + assert(intersecting_1 == intersecting_2); + assert(intersecting_1 == expected); + + std::cout << filename << " passed the tests." << std::endl << std::endl; + + return 0; +} + +template +int test_limited_self_intersections(const std::string& filename) +{ + std::vector points; + std::vector< std::array > triangles; + + bool read_ok = CGAL::IO::read_polygon_soup(filename, points, triangles); + + if ( !read_ok ) { + std::cerr << "Error: cannot read file: " << filename << std::endl; + return 1; + } + + CGAL::Timer timer; + timer.start(); + + std::vector > intersected_tris; + +#ifdef CGAL_LINKED_WITH_TBB + PMP::triangle_soup_self_intersections()( + points, triangles, + std::back_inserter(intersected_tris), CGAL::parameters::maximum_number(40)); + std::cout << "self_intersections test for 40 SI took " << timer.time() << " sec." << std::endl; + std::cout << "Found " << intersected_tris.size() << " SIs." << std::endl; + if(intersected_tris.size() < 40) + { + std::cerr<<"Not enough intersections found in parallel."<( + points, triangles, + std::back_inserter(intersected_tris), CGAL::parameters::maximum_number(40)); + std::cout << "self_intersections test for 40 SI took " << timer.time() << " sec." << std::endl; + timer.reset(); + if(intersected_tris.size() != 40) + { + std::cerr<<"Too many intersections found in sequential"< self_intersection_surface_mesh_test data/U.off false + + // First test ---------------------------------------------------------------- + bool expected = false; + std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off"); + if(argc > 1) { + assert(argc > 2); + std::stringstream ss(argv[2]); + ss >> std::boolalpha >> expected; + assert(!ss.fail()); // make sure that argv[2] is either 'true' or 'false' + } + + std::cout << "First test (EPICK):" << std::endl; + int r = test_self_intersections(filename, expected); + + std::cout << "First test (EPECK):" << std::endl; + r += test_self_intersections(filename, expected); + + // Second test --------------------------------------------------------------- + expected = true; + filename = (argc > 3) ? argv[3] : CGAL::data_file_path("meshes/mannequin-devil.off"); + if(argc > 3) { + assert(argc > 4); + std::stringstream ss(argv[4]); + ss >> std::boolalpha >> expected; + assert(!ss.fail()); + } + + std::cout << "Second test (EPICK):" << std::endl; + r += test_self_intersections(filename, expected); + + std::cout << "Second test (EPECK):" << std::endl; + r += test_self_intersections(filename, expected); + + // Third test ---------------------------------------------------------------- + expected = true; + filename = (argc > 5) ? argv[5] : "data/overlapping_triangles.off"; + if(argc > 5) { + assert(argc > 6); + std::stringstream ss(argv[6]); + ss >> std::boolalpha >> expected; + assert(!ss.fail()); + } + + std::cout << "Third test (EPICK):" << std::endl; + r += test_self_intersections(filename, expected); + + std::cout << "Third test (EPECK):" << std::endl; + r += test_self_intersections(filename, expected); + + // Fourth test ---------------------------------------------------------------- + expected = true; + filename = (argc > 7) ? argv[7] : "data_degeneracies/degtri_single.off"; + if(argc > 7) { + assert(argc > 8); + std::stringstream ss(argv[8]); + ss >> std::boolalpha >> expected; + assert(!ss.fail()); + } + + std::cout << "Fourth test (EPICK):" << std::endl; + r += test_self_intersections(filename, expected); + + std::cout << "Fourth test (EPECK):" << std::endl; + r += test_self_intersections(filename, expected); + + filename = (argc > 9) ? argv[9] : CGAL::data_file_path("meshes/mannequin-devil.off"); + + std::cout << "Test with maximum_number (EPICK):" << std::endl; + r += test_limited_self_intersections(filename); + + std::cout << "Test with maximum_number (EPECK):" << std::endl; + r += test_limited_self_intersections(filename); + + // extra hand written tests + { + // shared edge with same point: no self-intersection + typedef EPICK::Point_3 Point_3; + std::vector points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,-1,0)}; + std::vector< std::array > triangles={{0,1,2}, {0,1,3}}; + assert(!PMP::does_triangle_soup_self_intersect(points, triangles)); + } + { + // shared edge with duplicated points: self-intersection + typedef EPICK::Point_3 Point_3; + std::vector points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0),Point_3(0,0,0), Point_3(1,0,0), Point_3(0,-1,0)}; + std::vector< std::array > triangles={{0,1,2}, {3,4,5}}; + assert(PMP::does_triangle_soup_self_intersect(points, triangles)); + } + { + // shared vertex with same point: no self-intersection + typedef EPICK::Point_3 Point_3; + std::vector points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,2,0), Point_3(1,2,0)}; + std::vector< std::array > triangles={{0,1,2}, {3,4,2}}; + assert(!PMP::does_triangle_soup_self_intersect(points, triangles)); + } + { + // shared vertex with duplicated points: self-intersection + typedef EPICK::Point_3 Point_3; + std::vector points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,2,0), Point_3(1,2,0), Point_3(0,1,0)}; + std::vector< std::array > triangles={{0,1,2}, {3,4,5}}; + assert(PMP::does_triangle_soup_self_intersect(points, triangles)); + } + { + // 4 triangles around a shared edge: no self-intersection + typedef EPICK::Point_3 Point_3; + std::vector points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,-1,0), Point_3(0,0,1), Point_3(0,0,-1)}; + std::vector< std::array > triangles={{0,1,2},{0,1,3},{0,1,4},{0,1,5}}; + assert(!PMP::does_triangle_soup_self_intersect(points, triangles)); + } + { + // 4 triangles around a shared edge but two triangles intersecting: self-intersection + typedef EPICK::Point_3 Point_3; + std::vector points = {Point_3(0,0,0), Point_3(1,0,0), Point_3(0,1,0), Point_3(0,0.5,0), Point_3(0,0,1), Point_3(0,0,-1)}; + std::vector< std::array > triangles={{0,1,2},{0,1,3},{0,1,4},{0,1,5}}; + assert(PMP::does_triangle_soup_self_intersect(points, triangles)); + } + + return r; +}