From 2577fd912e4f3cb00dd62c79b45197047e0ee68c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 18 Jul 2016 15:02:06 +0200 Subject: [PATCH] WIP to add Hausdorff distance to a mesh --- .../PackageDescription.txt | 10 + .../CGAL/Polygon_mesh_processing/distance.h | 319 +++++++++++++++--- .../Polygon_mesh_processing/CMakeLists.txt | 9 + .../test_pmp_distance.cpp | 62 ++++ 4 files changed, 359 insertions(+), 41 deletions(-) create mode 100644 Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index dd4bbd9fd2d..b17816a7e2b 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -40,6 +40,12 @@ /// Functions to orient polygon soups and to stitch geometrically identical boundaries. /// \ingroup PkgPolygonMeshProcessing +/// \defgroup PMP_distance_grp Combinatorial Repairing +/// Functions to compute the distance between meshes. +/// \ingroup PkgPolygonMeshProcessing + + + /*! \addtogroup PkgPolygonMeshProcessing @@ -126,6 +132,10 @@ and provides a list of the parameters that are used in this package. - \link measure_grp `CGAL::Polygon_mesh_processing::edge_length()` \endlink - \link measure_grp `CGAL::Polygon_mesh_processing::face_border_length()` \endlink +## Distance Functions ## +- `CGAL::Polygon_mesh_processing::approximated_Hausdorff_distance()` +- `CGAL::Polygon_mesh_processing::approximated_symmetric_Hausdorff_distance()` + ## Miscellaneous ## - `CGAL::Polygon_mesh_slicer` - `CGAL::Side_of_triangle_mesh` diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h index ba90c8e4031..0490a784929 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -28,8 +28,18 @@ #include #include #include +#include +#include + + +#include + +#ifdef CGAL_LINKED_WITH_TBB +#include +#include +#include +#endif // CGAL_LINKED_WITH_TBB -/// \cond SKIP_IN_MANUAL namespace CGAL{ namespace Polygon_mesh_processing { @@ -37,17 +47,17 @@ namespace internal{ template OutputIterator -triangle_grid_sampling(const typename Kernel::Point_3& p0, - const typename Kernel::Point_3& p1, - const typename Kernel::Point_3& p2, - double distance, - OutputIterator out) +triangle_grid_sampling( const typename Kernel::Point_3& p0, + const typename Kernel::Point_3& p1, + const typename Kernel::Point_3& p2, + double distance, + OutputIterator out) { const double d_p0p1 = CGAL::sqrt( CGAL::squared_distance(p0, p1) ); const double d_p0p2 = CGAL::sqrt( CGAL::squared_distance(p0, p2) ); const double n = (std::max)(std::ceil( d_p0p1 / distance ), - std::ceil( d_p0p2 / distance )); + std::ceil( d_p0p2 / distance )); for (double i=1; i +#ifdef CGAL_LINKED_WITH_TBB +template +struct Distance_computation{ + const AABB_tree& tree; + const std::vector& sample_points; + cpp11::atomic* distance; + + Distance_computation(const AABB_tree& tree, const std::vector& sample_points, cpp11::atomic* d) + : tree(tree) + , sample_points(sample_points) + , distance(d) + {} + + void + operator()(const tbb::blocked_range& range) const + { + double hdist = 0; + Point_3 hint = sample_points.front(); + + for( std::size_t i = range.begin(); i != range.end(); ++i) + { + hint = tree.closest_point(sample_points[i], hint); + double d = CGAL::sqrt( squared_distance(hint,sample_points[i]) ); + if (d>hdist) hdist=d; + } + + if (hdist > distance->load(cpp11::memory_order_acquire)) + distance->store(hdist, cpp11::memory_order_release); + } +}; +#endif + +/// \todo update me to work with a triangulated surface mesh +template +double approximated_Hausdorff_distance( + std::vector& sample_points, + const TriangleRange& triangles) +{ + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "Nb sample points " << sample_points.size() << "\n"; + #endif + spatial_sort(sample_points.begin(), sample_points.end()); + + /// \todo shall we use Simple_cartesian for the distance computation? + /// check if this can be problematic to have non-exact predicates. + typedef typename TriangleRange::const_iterator Triangle_iterator; + typedef AABB_triangle_primitive Primitive; + typedef AABB_traits Traits; + typedef AABB_tree< Traits > Tree; + + Tree tree(triangles.begin(), triangles.end()); + tree.accelerate_distance_queries(); + tree.build(); + +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) + { + cpp11::atomic distance(0); + Distance_computation f(tree, sample_points, &distance); + tbb::parallel_for(tbb::blocked_range(0, sample_points.size()), f); + return distance; + } + else +#endif + { + double hdist = 0; + typename Traits::Point_3 hint = sample_points.front(); + BOOST_FOREACH(const typename Kernel::Point_3& pt, sample_points) + { + hint = tree.closest_point(pt, hint); + double d = CGAL::sqrt( squared_distance(hint,pt) ); + if (d>hdist) hdist=d; + } + return hdist; + } +} + +/// \todo remove me (except if you find an easy way avoid the ambiguity with TriangleMesh overload +template double approximated_Hausdorff_distance( const TriangleRange1& triangles_1, const TriangleRange2& triangles_2, @@ -124,32 +215,15 @@ double approximated_Hausdorff_distance( ) { std::vector sample_points; + sample_triangles( triangles_1, targeted_precision, std::back_inserter(sample_points) ); - // std::ofstream out("/tmp/samples.xyz"); - // std::copy(sample_points.begin(), sample_points.end(), std::ostream_iterator(out,"\n")); - - typedef typename TriangleRange2::const_iterator Triangle_iterator; - typedef AABB_triangle_primitive Primitive; - typedef AABB_traits Traits; - typedef AABB_tree< Traits > Tree; - - Tree tree(triangles_2.begin(), triangles_2.end()); - tree.accelerate_distance_queries(); - - double hdist = 0; - BOOST_FOREACH(const typename Kernel::Point_3& pt, sample_points) - { - double d = CGAL::sqrt( tree.squared_distance(pt) ); - // if ( d > 1e-1 ) std::cout << pt << "\n"; - if (d>hdist) hdist=d; - } - - return hdist; + return approximated_Hausdorff_distance(sample_points, triangles_2); } -template +/// \todo remove me (except if you find an easy way avoid the ambiguity with TriangleMesh overload +template double approximated_symmetric_Hausdorff_distance( const TriangleRange1& triangles_1, const TriangleRange2& triangles_2, @@ -157,26 +231,189 @@ double approximated_symmetric_Hausdorff_distance( ) { return (std::max)( - approximated_Hausdorff_distance(triangles_1, triangles_2, targeted_precision), - approximated_Hausdorff_distance(triangles_2, triangles_1, targeted_precision) + approximated_Hausdorff_distance(triangles_1, triangles_2, targeted_precision), + approximated_Hausdorff_distance(triangles_2, triangles_1, targeted_precision) ); } -///\todo add approximated_Hausdorff_distance(FaceGraph, FaceGraph) -///\todo add approximated_Hausdorff_distance(TriangleRange, FaceGraph) -///\todo add approximated_Hausdorff_distance(FaceGraph, TriangleRange) -///\todo add approximated_Hausdorff_distance(PointRange, FaceGraph) -///\todo add approximated_Hausdorff_distance(PointRange, TriangleRange) -///\todo add barycentric_triangle_sampling(Triangle, PointPerAreaUnit) -///\todo add random_triangle_sampling(Triangle, PointPerAreaUnit) -/// \todo maybe we should use a sampling using epec to avoid being too far from the sampled triangle - +/// \todo add a function `sample_triangle_mesh()` (or better name) that provide different sampling methods: +/// - random uniform ( points/area unit) +/// - grid (warning mininum 1 pt / triangle) (parameter = distance between points?) +/// - monte carlo? (random points in each triangle proportional to the face area with 1 pt mini per triangle) +/// The goal is to test different strategies and put the better one in `approximated_Hausdorff_distance()` +/// for particular cases one can still use a specific sampling method together with `max_distance_to_triangle_mesh()` + +/// \todo add a plugin in the demo to display the distance betweem 2 meshes as a texture (if not complicated) + +// documented functions + +/** + * \ingroup PMP_distance_grp + * computes the approximated Hausdorff distance of `tm1` from `tm2` by + * by generating a uniform random point sampling on `tm1`, and by then + * returning the distance of the furthest point from `tm2`. + * + * A parallel version is provided and requires the executable to be + * linked against the Intel TBB library. + * To control the number of threads used, the user may use the `tbb::task_scheduler_init` class. + * See the TBB documentation + * for more details. + * + * @tparam Concurrency_tag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag` + * and `Parallel_tag`. + * @tparam TriangleMesh a model of the concept `FaceListGraph` that has an internal property map + * for `CGAL::vertex_point_t` + * @tparam NamedParameters1 a sequence of \ref namedparameters for `tm1` + * @tparam NamedParameters2 a sequence of \ref namedparameters for `tm2` + * + * @param tm1 the triangle mesh that will be sampled + * @param tm2 the triangle mesh to compute the distance to + * @param precision TODO: either the number of points per squared area unit or something else depending on the result of the testing + * @param np1 optional sequence of \ref namedparameters for `tm1` among the ones listed below + * @param np2 optional sequence of \ref namedparameters for `tm2` among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh` \cgalParamEnd + * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd + * \cgalNamedParamsEnd + * \todo When using TBB, runtime is slower. The chunk size should probably be tuned. see PMP/test/test_pmp_distance.cpp + */ +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1, + class NamedParameters2> +double approximated_Hausdorff_distance( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double precision, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + typedef typename GetGeomTraits::type Geom_traits; + + /// \todo implement the missing function + return approximated_Hausdorff_distance( + tm1, tm2, + precision, + choose_const_pmap(get_param(np1, boost::vertex_point), + tm1, + vertex_point), + choose_const_pmap(get_param(np2, boost::vertex_point), + tm2, + vertex_point) + + ); +} + +/** + * \ingroup PMP_distance_grp + * computes the approximated symmetric Hausdorff distance between `tm1` and `tm2`. + * It returns the maximum of `approximated_Hausdorff_distance(tm1,tm2)` and + * `approximated_Hausdorff_distance(tm1,tm2)`. + * + * \copydetails CGAL::Polygon_mesh_processing::approximated_Hausdorff_distance() + */ +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1, + class NamedParameters2> +double approximated_symmetric_Hausdorff_distance( + const TriangleMesh& tm1, + const TriangleMesh& tm2, + double precision, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + return (std::max)( + approximated_Hausdorff_distance(tm1,tm2,precision,np1,np2), + approximated_Hausdorff_distance(tm2,tm1,precision,np2,np1) + ); +} + +/// \todo document and implement me +template< class Concurrency_tag, + class TriangleMesh, + class PointRange, + class NamedParameters> +double max_distance_to_triangle_mesh(const PointRange& points, + const TriangleMesh& tm, + const NamedParameters& np) +{ + return 0; +} + +/// \todo document and implement me +/// see with @sgiraudot for the implementation +/// that should be put in a seperate header file +/// with copyright INRIA +/// Maybe find a better name too +template< class Concurrency_tag, + class TriangleMesh, + class PointRange, + class NamedParameters> +double max_distance_to_point_set(const TriangleMesh& tm, + const PointRange& points, + const NamedParameters& np) +{ + return 0; +} + +// convenience functions with default parameters + +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1> +double approximated_Hausdorff_distance( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double precision, + const NamedParameters1& np1) +{ + return approximated_Hausdorff_distance( + tm1, tm2, precision, np1, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1> +double approximated_Hausdorff_distance( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double precision) +{ + return approximated_Hausdorff_distance( + tm1, tm2, precision, + parameters::all_default(), + parameters::all_default()); +} +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1> +double approximated_symmetric_Hausdorff_distance(const TriangleMesh& tm1, + const TriangleMesh& tm2, + double precision, + const NamedParameters1& np1) +{ + return approximated_symmetric_Hausdorff_distance( + tm1, tm2, precision, np1, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1> +double approximated_symmetric_Hausdorff_distance(const TriangleMesh& tm1, + const TriangleMesh& tm2, + double precision) +{ + return approximated_symmetric_Hausdorff_distance( + tm1, tm2, precision, + parameters::all_default(), + parameters::all_default()); +} } } // end of namespace CGAL::Polygon_mesh_processing -/// \endcond #endif //CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 7eb4eaa9848..69c5bb193ee 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -51,6 +51,14 @@ include_directories( BEFORE ../../include ) find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) +find_package( TBB ) +if( TBB_FOUND ) + include( ${TBB_USE_FILE} ) + list( APPEND CGAL_3RD_PARTY_LIBRARIES ${TBB_LIBRARIES} ) +else() + message( STATUS "NOTICE: Intel TBB was not found. test_pmp_distance will use sequential code." ) +endif() + if (EIGEN3_FOUND) # Executables that require Eigen 3.1 include( ${EIGEN3_USE_FILE} ) @@ -77,6 +85,7 @@ if (EIGEN3_FOUND) create_single_source_cgal_program("measures_test.cpp") create_single_source_cgal_program("triangulate_faces_test.cpp") create_single_source_cgal_program("test_pmp_remove_border_edge.cpp") + create_single_source_cgal_program("test_pmp_distance.cpp") if(NOT (${EIGEN3_VERSION} VERSION_LESS 3.2.0)) create_single_source_cgal_program("fairing_test.cpp") diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp new file mode 100644 index 00000000000..a35c6f3f656 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp @@ -0,0 +1,62 @@ +#include +#include +#include +#include +#include + + +#include + +#include +#include + +// typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Simple_cartesian K; +typedef CGAL::Surface_mesh Mesh; + + +int main(int argc, char** argv) +{ + Mesh m1,m2; + + std::ifstream input(argv[1]); + input >> m1; + input.close(); + + input.open(argv[2]); + input >> m2; + + std::cout << "First mesh has " << num_faces(m1) << " faces\n"; + std::cout << "Second mesh has " << num_faces(m2) << " faces\n"; + + std::vector t1; + t1.reserve(num_faces(m1)); + CGAL::Triangle_from_face_descriptor_map map1(&m1); + BOOST_FOREACH(Mesh::Face_index f, faces(m1)) + t1.push_back(get(map1,f)); + + std::vector t2; + t2.reserve(num_faces(m2)); + CGAL::Triangle_from_face_descriptor_map map2(&m2); + BOOST_FOREACH(Mesh::Face_index f, faces(m2)) + t2.push_back(get(map2,f)); + + + CGAL::Timer time; + time.start(); + std::cout << "Distance between meshes (parallel)" + << CGAL::Polygon_mesh_processing::approximated_Hausdorff_distance(t1,t2,0.001) + << "\n"; + time.stop(); + std::cout << "done in " << time.time() << "s.\n"; + + time.reset(); + time.start(); + std::cout << "Distance between meshes (sequential)" + << CGAL::Polygon_mesh_processing::approximated_Hausdorff_distance(t1,t2,0.001) + << "\n"; + time.stop(); + std::cout << "done in " << time.time() << "s.\n"; +} + +