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

Merge pull request #7370 from sloriot/PMP-soup_self_intersections 

Add self-intersection test and report functions for triangle soups
This commit is contained in:
Laurent Rineau 2023-04-21 11:20:52 +02:00
parent 91f0834dc5 2393258370
commit ad30839e6b
6 changed files with 728 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

115
BGL/include/CGAL/boost/graph/named_params_helper.h
View File

@ -110,22 +110,75 @@ public:
typedef typename boost::property_traits<PMap>::value_type type;
};
template <typename PolygonMesh,
typename VPM_from_NP>
struct GetVertexPointMap_impl
{
typedef VPM_from_NP type;
typedef VPM_from_NP const_type;
template<class NamedParameters>
static const_type
get_const_map(const NamedParameters& np, const PolygonMesh&)
{
return parameters::get_parameter(np, internal_np::vertex_point);
}
template<class NamedParameters>
static type
get_map(const NamedParameters& np, PolygonMesh&)
{
return parameters::get_parameter(np, internal_np::vertex_point);
}
};
template <typename PolygonMesh>
struct GetVertexPointMap_impl<PolygonMesh, internal_np::Param_not_found>
{
typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::const_type const_type;
typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::type type;
template<class NamedParameters>
static const_type
get_const_map(const NamedParameters& /* np */, const PolygonMesh& pm)
{
return get_const_property_map(boost::vertex_point, pm);
}
template<class NamedParameters>
static type
get_map(const NamedParameters& /* np */, PolygonMesh& pm)
{
return get_property_map(boost::vertex_point, pm);
}
};
template <typename PolygonMesh,
typename NamedParameters = parameters::Default_named_parameters>
class GetVertexPointMap
{
typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::const_type
DefaultVPMap_const;
typedef typename property_map_selector<PolygonMesh, boost::vertex_point_t>::type
DefaultVPMap;
typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_point_t,
NamedParameters,
internal_np::Param_not_found>::type VPM_from_NP;
typedef GetVertexPointMap_impl<PolygonMesh, VPM_from_NP> Impl;
public:
typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_point_t,
NamedParameters,
DefaultVPMap>::type type;
typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_point_t,
NamedParameters,
DefaultVPMap_const>::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<typename PolygonMesh, typename NamedParameters>
@ -138,10 +191,15 @@ public:
typedef typename CGAL::Kernel_traits<Point>::Kernel Kernel;
};
template <typename PolygonMesh,
typename NamedParametersGT = parameters::Default_named_parameters,
typename NamedParametersVPM = NamedParametersGT>
class GetGeomTraits
template<typename PolygonMesh, class GT, class NamedParametersVPM>
struct GetGeomTraits_impl
{
typedef GT type;
};
template<typename PolygonMesh, class NamedParametersVPM>
struct GetGeomTraits_impl<PolygonMesh, internal_np::Param_not_found, NamedParametersVPM>
{
typedef typename CGAL::graph_has_property<PolygonMesh, boost::vertex_point_t>::type Has_internal_pmap;
@ -154,12 +212,20 @@ class GetGeomTraits
typedef typename boost::mpl::if_c<Has_internal_pmap::value ||
!std::is_same<internal_np::Param_not_found, NP_vpm>::value,
typename GetK<PolygonMesh, NamedParametersVPM>::Kernel,
Fake_GT>::type DefaultKernel;
Fake_GT>::type type;
};
public:
template <typename PolygonMesh,
typename NamedParametersGT = parameters::Default_named_parameters,
typename NamedParametersVPM = NamedParametersGT>
struct GetGeomTraits
{
typedef typename internal_np::Lookup_named_param_def<internal_np::geom_traits_t,
NamedParametersGT,
DefaultKernel>::type type;
internal_np::Param_not_found>::type GT_from_NP;
typedef typename GetGeomTraits_impl<PolygonMesh,
GT_from_NP,
NamedParametersVPM>::type type;
};
// Define the following structs:
@ -278,6 +344,21 @@ public:
typedef typename CGAL::Identity_property_map<const Dummy_point> const_type;
};
template <typename PointRange, typename NamedParameters>
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<PointRange, NamedParameters>::type
>::value_type
>::type
> ::type type;
};
template <class PointRange, class NamedParameters, class PointMap = Default, class NormalMap = Default>
struct Point_set_processing_3_np_helper
{

2
Installation/CHANGES.md
View File

@ -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.

2
Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
View File

@ -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()`

460
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h
View File

@ -43,6 +43,7 @@
#endif
#include <boost/iterator/function_output_iterator.hpp>
#include <boost/range/irange.hpp>
#include <exception>
#include <sstream>
@ -54,6 +55,127 @@ namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
template <class TM>
struct Triangle_mesh_and_triangle_soup_wrapper
{
typedef typename boost::graph_traits<TM>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<TM>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TM>::halfedge_descriptor halfedge_descriptor; // private
static void get_face_vertices(face_descriptor fd, std::array<vertex_descriptor,3>& vh, const TM& tm)
{
CGAL_assertion(boost::graph_traits<TM>::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<vertex_descriptor, 4>& vh, const TM& tm)
{
CGAL_assertion(boost::graph_traits<TM>::null_face() != f);
CGAL_assertion(boost::graph_traits<TM>::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 <class PointRange, class TriangleRange>
struct Triangle_mesh_and_triangle_soup_wrapper<
std::pair<const PointRange&,
const TriangleRange&>>
{
typedef std::size_t face_descriptor;
typedef std::size_t vertex_descriptor;
typedef std::pair<const PointRange&, const TriangleRange& > Soup;
static void get_face_vertices(face_descriptor fd, std::array<vertex_descriptor,3>& 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<vertex_descriptor, 4>& 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<typename TriangleRange::const_iterator>::value_type& t : soup.second)
if (t.size()!=3)
return false;
return true;
}
};
template<typename Output_iterator>
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 <class GT, class TM, class VPM>
bool do_faces_intersect(typename boost::graph_traits<TM>::halfedge_descriptor h,
typename boost::graph_traits<TM>::halfedge_descriptor g,
bool do_faces_intersect(typename Triangle_mesh_and_triangle_soup_wrapper<TM>::face_descriptor fh,
typename Triangle_mesh_and_triangle_soup_wrapper<TM>::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<TM>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TM>::halfedge_descriptor halfedge_descriptor;
typedef Triangle_mesh_and_triangle_soup_wrapper<TM> 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<vertex_descriptor, 3> 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<vertex_descriptor, 4> 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 <class Box, class TM, class VPM, class GT,
class OutputIterator>
struct Strict_intersect_faces // "strict" as in "not sharing a subface"
{
typedef typename boost::graph_traits<TM>::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<GT>(h, g, m_tmesh, m_vpmap, m_construct_segment, m_construct_triangle, m_do_intersect))
if(do_faces_intersect<GT>(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<TM> 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<TM>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TM>::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<double, 3, face_descriptor, Box_policy> Box;
@ -245,9 +354,9 @@ self_intersections_impl(const FaceRange& face_range,
typedef typename GetGeomTraits<TM, NamedParameters>::type GT;
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename GetVertexPointMap<TM, NamedParameters>::const_type VPM;
VPM vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(boost::vertex_point, tmesh));
typedef GetVertexPointMap<TM, NamedParameters> 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<NamedParameters, internal_np::maximum_number_t>::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<vertex_descriptor, 3> vh;
Wrapper::get_face_vertices(f, vh, tmesh);
typename boost::property_traits<VPM>::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<unsigned int> 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<ConcurrencyTag>(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. <br>
* 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<ConcurrencyTag>(faces(tmesh), tmesh, np);
}
#ifndef DOXYGEN_RUNNING
template <class PointRange, class VPM>
struct Property_map_for_soup
{
typedef std::size_t key_type;
typedef typename boost::property_traits<VPM>::value_type value_type;
//typedef typename boost::property_traits<VPM>::category category;
typedef boost::readable_property_map_tag category;
typedef typename boost::property_traits<VPM>::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<PointRange, VPM>& 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<std::size_t,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 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.<br>
* 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 <class ConcurrencyTag = Sequential_tag,
class PointRange,
class TriangleRange,
class TriangleIdPairOutputIterator,
class CGAL_NP_TEMPLATE_PARAMETERS>
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<PointRange, CGAL_NP_CLASS>::const_type Point_map_base;
Point_map_base pm_base = choose_parameter<Point_map_base>(get_parameter(np, internal_np::point_map));
typedef Property_map_for_soup<PointRange, Point_map_base> Point_map;
typedef typename GetPolygonSoupGeomTraits<PointRange, CGAL_NP_CLASS>::type GT;
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
const bool do_limit = !(is_default_parameter<CGAL_NP_CLASS, internal_np::maximum_number_t>::value);
if (do_limit)
{
return self_intersections<ConcurrencyTag>(boost::irange<std::size_t>(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<ConcurrencyTag>(boost::irange<std::size_t>(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 <class ConcurrencyTag = Sequential_tag,
class PointRange,
class TriangleRange,
class CGAL_NP_TEMPLATE_PARAMETERS>
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<PointRange, CGAL_NP_CLASS>::const_type Point_map_base;
Point_map_base pm_base = choose_parameter<Point_map_base>(get_parameter(np, internal_np::point_map));
typedef Property_map_for_soup<PointRange, Point_map_base> Point_map;
typedef typename GetPolygonSoupGeomTraits<PointRange, CGAL_NP_CLASS>::type GT;
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
internal::self_intersections_impl<ConcurrencyTag>(boost::irange<std::size_t>(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

1
Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
View File

@ -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")

234
Polygon_mesh_processing/test/Polygon_mesh_processing/self_intersection_triangle_soup_test.cpp Normal file
View File

@ -0,0 +1,234 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <CGAL/tags.h>
#include <CGAL/Timer.h>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <sstream>
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 <typename K>
int test_self_intersections(const std::string filename,
const bool expected)
{
std::vector<typename K::Point_3> points;
std::vector< std::array<std::size_t, 3> > 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<std::pair<std::size_t, std::size_t> > intersected_tris;
PMP::triangle_soup_self_intersections<CGAL::Parallel_if_available_tag>(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<CGAL::Parallel_if_available_tag>(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 <typename K>
int test_limited_self_intersections(const std::string& filename)
{
std::vector<typename K::Point_3> points;
std::vector< std::array<std::size_t, 3> > 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<std::pair<std::size_t, std::size_t> > intersected_tris;
#ifdef CGAL_LINKED_WITH_TBB
PMP::triangle_soup_self_intersections()<CGAL::Parallel_if_available_tag>(
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."<<std::endl;
return 1;
}
intersected_tris.clear();
timer.reset();
#endif
PMP::triangle_soup_self_intersections<CGAL::Sequential_tag>(
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"<<std::endl;
return 1;
}
return 0;
}
int main(int argc, char** argv)
{
// If file(s) are provided, the associated expected result must also be provided.
// Note that this expected value is a Boolean that is passed in command line
// with either 'true' or 'false' (and not integers), that is for example:
// > 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<EPICK>(filename, expected);
std::cout << "First test (EPECK):" << std::endl;
r += test_self_intersections<EPECK>(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<EPICK>(filename, expected);
std::cout << "Second test (EPECK):" << std::endl;
r += test_self_intersections<EPECK>(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<EPICK>(filename, expected);
std::cout << "Third test (EPECK):" << std::endl;
r += test_self_intersections<EPECK>(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<EPICK>(filename, expected);
std::cout << "Fourth test (EPECK):" << std::endl;
r += test_self_intersections<EPECK>(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<EPICK>(filename);
std::cout << "Test with maximum_number (EPECK):" << std::endl;
r += test_limited_self_intersections<EPECK>(filename);
// extra hand written tests
{
// shared edge with same point: no self-intersection
typedef EPICK::Point_3 Point_3;
std::vector<EPICK::Point_3> 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<std::size_t, 3> > 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<EPICK::Point_3> 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<std::size_t, 3> > 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<EPICK::Point_3> 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<std::size_t, 3> > 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<EPICK::Point_3> 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<std::size_t, 3> > 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<EPICK::Point_3> 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<std::size_t, 3> > 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<EPICK::Point_3> 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<std::size_t, 3> > triangles={{0,1,2},{0,1,3},{0,1,4},{0,1,5}};
assert(PMP::does_triangle_soup_self_intersect(points, triangles));
}
return r;
}