diff --git a/AABB_tree/include/CGAL/AABB_tree.h b/AABB_tree/include/CGAL/AABB_tree.h index 1bb4d0b2d172..f3c8fe9ad35c 100644 --- a/AABB_tree/include/CGAL/AABB_tree.h +++ b/AABB_tree/include/CGAL/AABB_tree.h @@ -527,6 +527,36 @@ namespace CGAL { } } + template + void traversal_with_priority(const Query& query, Traversal_traits& traits) const + { + switch(size()) + { + case 0: + break; + case 1: + traits.intersection(query, singleton_data()); + break; + default: // if(size() >= 2) + root_node()->template traversal_with_priority(query, traits, m_primitives.size()); + } + } + + template + void traversal_with_priority_and_group_traversal(const Query& query, Traversal_traits& traits, const std::size_t group_traversal_bound) const + { + switch(size()) + { + case 0: + break; + case 1: + traits.intersection(query, singleton_data()); + break; + default: // if(size() >= 2) + root_node()->template traversal_with_priority_and_group_traversal(m_primitives, query, traits, m_primitives.size(), 0, group_traversal_bound); + } + } + private: typedef AABB_node Node; diff --git a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h index d8edaa91d358..a5c09f42fb46 100644 --- a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h +++ b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h @@ -68,6 +68,20 @@ class AABB_node Traversal_traits& traits, const std::size_t nb_primitives) const; + template + void traversal_with_priority(const Query& query, + Traversal_traits& traits, + const std::size_t nb_primitives) const; + + template + void traversal_with_priority_and_group_traversal(const Primitive_vector& primitives, + const Query& query, + Traversal_traits& traits, + const std::size_t nb_primitives, + std::size_t first_primitive_index, + const std::size_t group_size_bound) const; + + private: typedef AABBTraits AABB_traits; typedef AABB_node Node; @@ -152,6 +166,156 @@ AABB_node::traversal(const Query& query, } } +template +template +void +AABB_node::traversal_with_priority(const Query& query, + Traversal_traits& traits, + const std::size_t nb_primitives) const +{ + // Recursive traversal + switch(nb_primitives) + { + case 2: + traits.intersection(query, left_data()); + if( traits.go_further() ) + { + traits.intersection(query, right_data()); + } + break; + case 3: + traits.intersection(query, left_data()); + if( traits.go_further() && traits.do_intersect(query, right_child()) ) + { + right_child().traversal_with_priority(query, traits, 2); + } + break; + default: + bool ileft, iright; + typename Traversal_traits::Priority pleft, pright; + std::tie(ileft, pleft) = traits.do_intersect_with_priority(query, left_child()); + std::tie(iright, pright) = traits.do_intersect_with_priority(query, right_child()); + CGAL_precondition( (ileft || iright) ? traits.do_intersect(query, *this) : true ); + + if(ileft) + { + if(iright) + { + // Both children have to be inspected. + if(pleft >= pright) + { + // Inspect the left child first, has higher priority. + left_child().traversal_with_priority(query, traits, nb_primitives/2); + if( traits.go_further() ) + right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + } + else + { + // Inspect the right child first, has higher priority. + right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + if( traits.go_further() ) + left_child().traversal_with_priority(query, traits, nb_primitives/2); + } + } + else + { + // Only the left child has to be inspected. + left_child().traversal_with_priority(query, traits, nb_primitives/2); + } + } + else + { + if(iright) + { + // Only the right child has to be inspected. + right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + } + } + } +} + +// TODO: find a better name +template +template +void +AABB_node::traversal_with_priority_and_group_traversal(const Primitive_vector& primitives, + const Query& query, + Traversal_traits& traits, + const std::size_t nb_primitives, + std::size_t first_primitive_index, + const std::size_t group_traversal_bound) const +{ + // Group traversal + CGAL_assertion(group_traversal_bound >= 2); + if ( nb_primitives <= group_traversal_bound ) + { + if ( !traits.do_intersect(query, *this) ) return; + CGAL_assertion(traits.do_intersect(query, *this)); + traits.traverse_group(query, primitives.begin()+first_primitive_index, primitives.begin()+first_primitive_index+nb_primitives); + return; + } + + // Recursive traversal + switch(nb_primitives) + { + case 2: + traits.intersection(query, left_data()); + if( traits.go_further() ) + { + traits.intersection(query, right_data()); + } + break; + case 3: + traits.intersection(query, left_data()); + if( traits.go_further() && traits.do_intersect(query, right_child()) ) + { + right_child().traversal_with_priority_and_group_traversal(primitives, query, traits, 2, first_primitive_index+1, group_traversal_bound); + } + break; + default: + bool ileft, iright; + typename Traversal_traits::Priority pleft, pright; + std::tie(ileft, pleft) = traits.do_intersect_with_priority(query, left_child()); + std::tie(iright, pright) = traits.do_intersect_with_priority(query, right_child()); + CGAL_precondition( (ileft || iright) ? traits.do_intersect(query, *this) : true ); + + if(ileft) + { + if(iright) + { + // Both children have to be inspected. + if(pleft >= pright) + { + // Inspect the left child first, has higher priority. + left_child().traversal_with_priority_and_group_traversal(primitives, query, traits, nb_primitives/2, first_primitive_index, group_traversal_bound); + if( traits.go_further() ) + right_child().traversal_with_priority_and_group_traversal(primitives, query, traits, nb_primitives-nb_primitives/2, first_primitive_index+nb_primitives/2, group_traversal_bound); + } + else + { + // Inspect the right child first, has higher priority. + right_child().traversal_with_priority_and_group_traversal(primitives, query, traits, nb_primitives-nb_primitives/2, first_primitive_index+nb_primitives/2, group_traversal_bound); + if( traits.go_further() ) + left_child().traversal_with_priority_and_group_traversal(primitives, query, traits, nb_primitives/2, first_primitive_index, group_traversal_bound); + } + } + else + { + // Only the left child has to be inspected. + left_child().traversal_with_priority_and_group_traversal(primitives, query, traits, nb_primitives/2, first_primitive_index, group_traversal_bound); + } + } + else + { + if(iright) + { + // Only the right child has to be inspected. + right_child().traversal_with_priority_and_group_traversal(primitives, query, traits, nb_primitives-nb_primitives/2, first_primitive_index+nb_primitives/2, group_traversal_bound); + } + } + } +} + } // end namespace CGAL #endif // CGAL_AABB_NODE_H diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h index a1fe60c5f75b..a69d0d4a965f 100644 --- a/BGL/include/CGAL/boost/graph/parameters_interface.h +++ b/BGL/include/CGAL/boost/graph/parameters_interface.h @@ -121,7 +121,9 @@ CGAL_add_named_parameter(do_not_modify_t, do_not_modify, do_not_modify) CGAL_add_named_parameter(allow_self_intersections_t, allow_self_intersections, allow_self_intersections) CGAL_add_named_parameter(non_manifold_feature_map_t, non_manifold_feature_map, non_manifold_feature_map) CGAL_add_named_parameter(polyhedral_envelope_epsilon_t, polyhedral_envelope_epsilon, polyhedral_envelope_epsilon) +CGAL_add_named_parameter(match_faces_t, match_faces, match_faces) CGAL_add_named_parameter(face_epsilon_map_t, face_epsilon_map, face_epsilon_map) +CGAL_add_named_parameter(use_one_sided_hausdorff_t, use_one_sided_hausdorff, use_one_sided_hausdorff) // List of named parameters that we use in the package 'Surface Mesh Simplification' CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost) diff --git a/Documentation/doc/biblio/geom.bib b/Documentation/doc/biblio/geom.bib index ca708451d71d..4b7a0483d271 100644 --- a/Documentation/doc/biblio/geom.bib +++ b/Documentation/doc/biblio/geom.bib @@ -152050,4 +152050,14 @@ @article{cvl-ew-12 Pages = {215--224}, Year = {2012}, Url = {http://monge.univ-mlv.fr/~colinde/pub/09edgewidth.pdf} + +@inproceedings{tang2009interactive, + title={Interactive Hausdorff distance computation for general polygonal models}, + author={Tang, Min and Lee, Minkyoung and Kim, Young J}, + booktitle={ACM Transactions on Graphics (TOG)}, + volume={28}, + number={3}, + pages={74}, + year={2009}, + organization={ACM} } diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index dd07172543d1..7ac7d589dbaa 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -7,7 +7,17 @@ Release date: December 2021 ### [Polygon Mesh Processing](https://doc.cgal.org/5.4/Manual/packages.html#PkgPolygonMeshProcessing) -- Added the function `CGAL::Polygon_mesh_processing::match_faces()`, which, given two polygon meshes, identifies their common faces as well as as faces present in only either of them. +- Added the function `CGAL::Polygon_mesh_processing::match_faces()`, which, given two polygon meshes, + identifies their common faces as well as faces present in only either of them. + +- Added the functions: `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` that + computes an estimate of the one-sided Hausdorff distance between two triangle meshes which + is bounded by a user-specified error bound; `CGAL::Polygon_mesh_processing::bounded_error_symmetric_Hausdorff_distance()` that computes + an estimate of the symmetric Hausdorff distance bounded by a user-specified error bound; + and `CGAL::Polygon_mesh_processing::is_Hausdorff_distance_larger()` that returns `true` + if the bounded-error Hausdorff distance between two meshes is larger than the user-specified + max distance. + [Release 5.3](https://github.com/CGAL/cgal/releases/tag/v5.3) ----------- diff --git a/Installation/cmake/modules/CGAL_METIS_support.cmake b/Installation/cmake/modules/CGAL_METIS_support.cmake index 563bec6fc7f7..c3fa73714d43 100644 --- a/Installation/cmake/modules/CGAL_METIS_support.cmake +++ b/Installation/cmake/modules/CGAL_METIS_support.cmake @@ -1,6 +1,7 @@ if(METIS_FOUND AND NOT TARGET CGAL::METIS_support) add_library(CGAL::METIS_support INTERFACE IMPORTED) set_target_properties(CGAL::METIS_support PROPERTIES + INTERFACE_COMPILE_DEFINITIONS "CGAL_METIS_ENABLED" INTERFACE_INCLUDE_DIRECTORIES "${METIS_INCLUDE_DIRS}" INTERFACE_LINK_LIBRARIES "${METIS_LIBRARIES}") endif() diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index 9edc40b6e020..e678bf1627fb 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -201,6 +201,9 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - \link measure_grp `CGAL::Polygon_mesh_processing::match_faces()` \endlink \cgalCRPSection{Distance Functions} +- `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` +- `CGAL::Polygon_mesh_processing::bounded_error_symmetric_Hausdorff_distance()` +- `CGAL::Polygon_mesh_processing::is_Hausdorff_distance_larger()` - `CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance()` - `CGAL::Polygon_mesh_processing::approximate_symmetric_Hausdorff_distance()` - `CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set()` diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index bfd466412fba..8757bba98d48 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt @@ -910,10 +910,12 @@ which enables to treat one or several connected components as a face graph. \cgalExample{Polygon_mesh_processing/face_filtered_graph_example.cpp} -\section PMPDistance Approximate Hausdorff Distance +\section PMPDistance Hausdorff Distance This package provides methods to compute (approximate) distances between meshes and point sets. +\subsection ApproxHD Approximate Hausdorff Distance + The function \link approximate_Hausdorff_distance() `approximate_Hausdorff_distance()`\endlink computes an approximation of the Hausdorff distance from a mesh `tm1` to a mesh `tm2`. Given a a sampling of `tm1`, it computes the distance to `tm2` of the farthest sample point to `tm2` \cgalCite{cignoni1998metro}. @@ -938,12 +940,12 @@ computes an approximation of the Hausdorff distance from a mesh to a point set. For each triangle, a lower and upper bound of the Hausdorff distance to the point set are computed. Triangles are refined until the difference between the bounds is lower than a user-defined precision threshold. -\subsection AHDExample Approximate Hausdorff Distance Example +\subsubsection AHDExample Approximate Hausdorff Distance Example In the following example, a mesh is isotropically remeshed and the approximate distance between the input and the output is computed. \cgalExample{Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp} -\subsection PoissonDistanceExample Max Distance Between Point Set and Surface Example +\subsubsection PoissonDistanceExample Max Distance Between Point Set and Surface Example In \ref Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp, a triangulated surface mesh is constructed from a point set using the \link PkgPoissonSurfaceReconstruction3 Poisson reconstruction algorithm \endlink, @@ -952,6 +954,50 @@ with the following code: \snippet Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp PMP_distance_snippet +\subsection BoundedHD Bounded Hausdorff Distance + +The function `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` +computes an estimate of the Hausdorff distance of two triangle meshes which is +bounded by a user-given error bound. Given two meshes `tm1` and `tm2`, it follows +the procedure given by \cgalCite{tang2009interactive}. Namely, a bounded volume hierarchy (BVH) is +built on `tm1` and `tm2` respectively. The BVH on `tm1` is used to iterate over all +triangles in `tm1`. Throughout the traversal, the procedure keeps track of a global +lower and upper bound on the Hausdorff distance respectively. For each triangle +`t` in `tm1`, by traversing the BVH on `tm2`, it is estimated via the global bounds +whether `t` can still contribute to the actual Hausdorff distance. From this +process, a set of candidate triangles is selected. + +The candidate triangles are subsequently subdivided and for each smaller +triangle, the BVH on `tm2` is traversed again. This is repeated until the +triangle is smaller than the user-given error bound, all vertices of the +triangle are projected onto the same triangle in `tm2`, or the triangle's upper +bound is lower than the global lower bound. After creation, the subdivided +triangles are added to the list of candidate triangles. Thereby, all candidate +triangles are processed until a triangle is found in which the Hausdorff +distance is realized or in which it is guaranteed to be realized within the +user-given error bound. + +In the current implementation, the BVH used is an AABB-tree and not the swept sphere +volumes as used in the original implementation. This should explain the runtime difference +observed with the original implementation. + +The function `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` computes +the one-sided Hausdorff distance from `tm1` to `tm2`. This component also provides +the symmetric distance `CGAL::Polygon_mesh_processing::bounded_error_symmetric_Hausdorff_distance()` +and an utility function called `CGAL::Polygon_mesh_processing::is_Hausdorff_distance_larger()` +that returns `true` if the Hausdorff distance between two meshes is larger than the user-defined max distance. + +\subsubsection BHDExample Bounded Hausdorff Distance Example + +In the following examples: (a) the distance of a tetrahedron to a remeshed +version of itself is computed, (b) the distance of two geometries is computed +which is realized strictly in the interior of a triangle of the first geometry, +(c) a perturbation of a user-given mesh is compared to the original user-given +mesh, (d) two user-given meshes are compared, where the second mesh is gradually +moved away from the first one. + +\cgalExample{Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp} + \section PMPDetectFeatures Feature Detection This package provides methods to detect some features of a polygon mesh. diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt index 699614c41d0a..496897e23492 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt @@ -34,4 +34,5 @@ \example Polygon_mesh_processing/locate_example.cpp \example Polygon_mesh_processing/orientation_pipeline_example.cpp \example Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp +\example Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp */ diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index ad9c683900e9..99ebc9a5b6d8 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -37,6 +37,7 @@ include(CGAL_Eigen3_support) # ########################################################## create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp") +create_single_source_cgal_program( "hausdorff_bounded_error_distance_example.cpp") if(TARGET CGAL::Eigen3_support) create_single_source_cgal_program("hole_filling_example.cpp") @@ -122,12 +123,23 @@ if(OpenMesh_FOUND) PRIVATE ${OPENMESH_LIBRARIES}) endif(OpenMesh_FOUND) +find_package(METIS) +include(CGAL_METIS_support) +if(TARGET CGAL::METIS_support) + target_link_libraries(hausdorff_bounded_error_distance_example PUBLIC CGAL::METIS_support) +else() + message(STATUS "Tests, which use the METIS library will not be compiled.") +endif() + find_package(TBB) include(CGAL_TBB_support) if(TARGET CGAL::TBB_support) target_link_libraries(self_intersections_example PUBLIC CGAL::TBB_support) target_link_libraries(hausdorff_distance_remeshing_example PUBLIC CGAL::TBB_support) + target_link_libraries(hausdorff_bounded_error_distance_example + PUBLIC CGAL::TBB_support) + create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp") target_link_libraries(corefinement_parallel_union_meshes PUBLIC CGAL::TBB_support) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp new file mode 100644 index 000000000000..fee57a590a3c --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -0,0 +1,67 @@ +#include +#include +#include + +#include +#include +#include +#include +#include + +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using FT = typename Kernel::FT; +using Point_3 = typename Kernel::Point_3; +using Vector_3 = typename Kernel::Vector_3; + +using TAG = CGAL::Sequential_tag; +using Surface_mesh = CGAL::Surface_mesh; +using Polyhedron = CGAL::Polyhedron_3; +using Affine_transformation_3 = CGAL::Aff_transformation_3; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char** argv) { + + const double error_bound = 1e-4; + const std::string filepath = (argc > 1 ? argv[1] : "data/blobby.off"); + + // We create a tetrahedron, remesh it, and compute the distance. + // The expected distance is error_bound. + std::cout << std::endl << "* remeshing tetrahedron example:" << std::endl; + + Surface_mesh mesh1, mesh2; + CGAL::make_tetrahedron( + Point_3(0, 0, 0), Point_3(2, 0, 0), + Point_3(1, 1, 1), Point_3(1, 0, 2), mesh1); + mesh2 = mesh1; + + using edge_descriptor = typename boost::graph_traits::edge_descriptor; + Surface_mesh::Property_map is_constrained_map = + mesh2.add_property_map("e:is_constrained", true).first; + + const double target_edge_length = 0.05; + PMP::isotropic_remeshing( + mesh2.faces(), target_edge_length, mesh2, + PMP::parameters::edge_is_constrained_map(is_constrained_map)); + + std::cout << "* one-sided bounded-error Hausdorff distance: " << + PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound) << std::endl; + + // We load a mesh, save it in two different containers, and + // translate the second mesh by 1 unit. The expected distance is 1. + std::cout << std::endl << "* moving mesh example:" << std::endl; + + Surface_mesh surface_mesh; + CGAL::IO::read_OFF(filepath, surface_mesh); + + Polyhedron polyhedron; + CGAL::IO::read_OFF(filepath, polyhedron); + + PMP::transform(Affine_transformation_3(CGAL::Translation(), + Vector_3(FT(0), FT(0), FT(1))), polyhedron); + + std::cout << "* symmetric bounded-error Hausdorff distance: " << + PMP::bounded_error_symmetric_Hausdorff_distance(surface_mesh, polyhedron, error_bound) + << std::endl << std::endl; + return EXIT_SUCCESS; +} 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 b945dccf702b..b65c7a414b90 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -8,7 +8,7 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // -// Author(s) : Maxime Gimeno and Sebastien Loriot +// Author(s) : Maxime Gimeno, Sebastien Loriot, Martin Skrodzki, Dmitry Anisimov #ifndef CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H #define CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H @@ -16,7 +16,9 @@ #include #include +#include #include +#include #include #include @@ -26,7 +28,15 @@ #include #include #include +#include #include +#include +#include + +#include +#if defined(CGAL_METIS_ENABLED) +#include +#endif // CGAL_METIS_ENABLED #ifdef CGAL_LINKED_WITH_TBB #include @@ -35,10 +45,12 @@ #endif // CGAL_LINKED_WITH_TBB #include +#include #include #include #include +#include namespace CGAL { namespace Polygon_mesh_processing { @@ -1019,12 +1031,12 @@ double approximate_Hausdorff_distance( CGAL_assertion_code( bool is_triangle = is_triangle_mesh(tm) ); CGAL_assertion_msg (is_triangle, "Mesh is not triangulated. Distance computing impossible."); - #ifdef CGAL_HAUSDORFF_DEBUG - std::cout << "Nb sample points " << sample_points.size() << "\n"; - #endif typedef typename Kernel::Point_3 Point_3; std::vector sample_points (boost::begin(original_sample_points), boost::end(original_sample_points) ); + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "Nb sample points " << sample_points.size() << "\n"; + #endif spatial_sort(sample_points.begin(), sample_points.end()); @@ -1310,8 +1322,1379 @@ double approximate_symmetric_Hausdorff_distance(const TriangleMesh& tm1, tm1, tm2, parameters::all_default(), parameters::all_default()); } +//////////////////////////////////////////////////////////////////////// + +// Use this def in order to get back the parallel version of the one-sided Hausdorff code! +// #define USE_PARALLEL_BEHD + +// Use this def in order to get all DEBUG info related to the bounded-error Hausdorff code! +// #define CGAL_HAUSDORFF_DEBUG + +namespace internal { + +template< class Kernel, + class TriangleMesh1, + class TriangleMesh2, + class VPM1, + class VPM2, + class NamedParameters1, + class NamedParameters2, + class TM1Tree, + class TM2Tree, + class FaceHandle1, + class FaceHandle2 > +std::pair preprocess_bounded_error_Hausdorff_impl( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const bool compare_meshes, + const VPM1& vpm1, + const VPM2& vpm2, + const bool is_one_sided_distance, + const NamedParameters1& np1, + const NamedParameters2& np2, + TM1Tree& tm1_tree, + TM2Tree& tm2_tree, + std::vector& tm1_only, + std::vector& tm2_only) +{ + using FT = typename Kernel::FT; + using Point_3 = typename Kernel::Point_3; + + #ifdef CGAL_HAUSDORFF_DEBUG + using Timer = CGAL::Real_timer; + Timer timer; + timer.start(); + std::cout << "* preprocessing begin ...." << std::endl; + std::cout.precision(17); + #endif + + // Compute the max value that is used as infinity value for the given meshes. + // In our case, it is twice the length of the diagonal of the bbox of two input meshes. + const auto bbox1 = bbox(tm1); + const auto bbox2 = bbox(tm2); + const auto bb = bbox1 + bbox2; + const FT sq_dist = CGAL::squared_distance( + Point_3(bb.xmin(), bb.ymin(), bb.zmin()), + Point_3(bb.xmax(), bb.ymax(), bb.zmax())); + FT infinity_value = CGAL::approximate_sqrt(sq_dist) * FT(2); + CGAL_assertion(infinity_value >= FT(0)); + + // Compare meshes and build trees. + tm1_only.clear(); + tm2_only.clear(); + std::vector< std::pair > common; + + const auto faces1 = faces(tm1); + const auto faces2 = faces(tm2); + + CGAL_precondition(faces1.size() > 0); + CGAL_precondition(faces2.size() > 0); + + // Compare meshes. + bool rebuild = false; + if (compare_meshes) { // exact check + match_faces(tm1, tm2, std::back_inserter(common), + std::back_inserter(tm1_only), std::back_inserter(tm2_only), np1, np2); + + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "- common: " << common.size() << std::endl; + std::cout << "- tm1 only: " << tm1_only.size() << std::endl; + std::cout << "- tm2 only: " << tm2_only.size() << std::endl; + #endif + + if (is_one_sided_distance) { // one-sided distance + + if (tm1_only.size() > 0) { // create TM1 and and full TM2 + tm1_tree.insert(tm1_only.begin(), tm1_only.end(), tm1, vpm1); + tm2_tree.insert(faces2.begin(), faces2.end(), tm2, vpm2); + } else { // do not create trees + CGAL_assertion(tm1_only.size() == 0); + infinity_value = -FT(1); + } + + } else { // symmetric distance + + if (tm1_only.size() == 0 && tm2_only.size() == 0) { // do not create trees + infinity_value = -FT(1); + } else if (common.size() == 0) { // create full TM1 and TM2 + tm1_tree.insert(faces1.begin(), faces1.end(), tm1, vpm1); + tm2_tree.insert(faces2.begin(), faces2.end(), tm2, vpm2); + } else if (tm1_only.size() == 0) { // create TM2 and full TM1 + CGAL_assertion(tm2_only.size() > 0); + CGAL_assertion(tm2_only.size() < faces2.size()); + tm1_tree.insert(faces1.begin(), faces1.end(), tm1, vpm1); + tm2_tree.insert(tm2_only.begin(), tm2_only.end(), tm2, vpm2); + } else if (tm2_only.size() == 0) { // create TM1 and full TM2 + CGAL_assertion(tm1_only.size() > 0); + CGAL_assertion(tm1_only.size() < faces1.size()); + tm1_tree.insert(tm1_only.begin(), tm1_only.end(), tm1, vpm1); + tm2_tree.insert(faces2.begin(), faces2.end(), tm2, vpm2); + } else { // create TM1 and full TM2 and set tag to rebuild them later + CGAL_assertion(tm1_only.size() > 0); + CGAL_assertion(tm1_only.size() < faces1.size()); + tm1_tree.insert(tm1_only.begin(), tm1_only.end(), tm1, vpm1); + tm2_tree.insert(faces2.begin(), faces2.end(), tm2, vpm2); + rebuild = true; + } + } + } else { // create full TM1 and TM2 + tm1_tree.insert(faces1.begin(), faces1.end(), tm1, vpm1); + tm2_tree.insert(faces2.begin(), faces2.end(), tm2, vpm2); + } + + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + std::cout << "* .... end preprocessing" << std::endl; + std::cout << "* preprocessing time (sec.): " << timer.time() << std::endl; + #endif + return std::make_pair(infinity_value, rebuild); +} + +template< class Kernel, + class TriangleMesh1, + class TriangleMesh2, + class VPM1, + class VPM2, + class TM1Tree, + class TM2Tree, + class OutputIterator > +double bounded_error_Hausdorff_impl( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const typename Kernel::FT error_bound, + const VPM1& vpm1, + const VPM2& vpm2, + const typename Kernel::FT infinity_value, + const typename Kernel::FT initial_bound, + const typename Kernel::FT distance_bound, + const TM1Tree& tm1_tree, + const TM2Tree& tm2_tree, + OutputIterator& out) +{ + using FT = typename Kernel::FT; + using Point_3 = typename Kernel::Point_3; + using Triangle_3 = typename Kernel::Triangle_3; + + using TM1_tree = TM1Tree; + using TM2_tree = TM2Tree; + + using TM1_traits = typename TM1_tree::AABB_traits; + using TM2_traits = typename TM2_tree::AABB_traits; + + using TM1_hd_traits = Hausdorff_primitive_traits_tm1; + using TM2_hd_traits = Hausdorff_primitive_traits_tm2; + + using Face_handle_1 = typename boost::graph_traits::face_descriptor; + using Face_handle_2 = typename boost::graph_traits::face_descriptor; + + using Candidate = Candidate_triangle; + + CGAL_precondition(error_bound >= FT(0)); + CGAL_precondition(tm1_tree.size() > 0); + CGAL_precondition(tm2_tree.size() > 0); + + // First, we apply culling. + #ifdef CGAL_HAUSDORFF_DEBUG + using Timer = CGAL::Real_timer; + Timer timer; + timer.start(); + std::cout << "- applying culling" << std::endl; + std::cout.precision(17); + #endif + + // Build traversal traits for tm1_tree. + TM1_hd_traits traversal_traits_tm1( + tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2, + error_bound, infinity_value, initial_bound, distance_bound); + + // Find candidate triangles in TM1, which might realise the Hausdorff bound. + // We build a sorted structure while collecting the candidates. + const Point_3 stub(0, 0, 0); // dummy point given as query since it is not needed + + tm1_tree.traversal_with_priority(stub, traversal_traits_tm1); + auto& candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); + auto global_bounds = traversal_traits_tm1.get_global_bounds(); + + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "- number of candidate triangles: " << candidate_triangles.size() << std::endl; + const FT culling_rate = FT(100) - (FT(candidate_triangles.size()) / FT(tm1_tree.size()) * FT(100)); + std::cout << "- culling rate: " << culling_rate << "%" << std::endl; + #endif + + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + std::cout << "* culling (sec.): " << timer.time() << std::endl; + #endif + + CGAL_assertion(global_bounds.lower >= FT(0)); + CGAL_assertion(global_bounds.upper >= FT(0)); + CGAL_assertion(global_bounds.upper >= global_bounds.lower); + + // If we already reached the user-defined max distance bound, we quit. + if (traversal_traits_tm1.early_quit()) { + CGAL_assertion(distance_bound >= FT(0)); + const double hdist = CGAL::to_double((global_bounds.lower + global_bounds.upper) / FT(2)); + return hdist; + } + CGAL_assertion(!traversal_traits_tm1.early_quit()); + + // Second, we apply subdivision. + #ifdef CGAL_HAUSDORFF_DEBUG + timer.reset(); + timer.start(); + std::cout << "- applying subdivision" << std::endl; + #endif + + // See Section 5.1 in the paper. + const FT squared_error_bound = error_bound * error_bound; + while ( + (global_bounds.upper - global_bounds.lower > error_bound) && + !candidate_triangles.empty()) { + + // Check if we can early quit. + if (distance_bound >= FT(0)) { + const FT hdist = (global_bounds.lower + global_bounds.upper) / FT(2); + const bool early_quit = (hdist >= distance_bound); + if (early_quit) break; + } + + // Get the first triangle and its Hausdorff bounds from the candidate set. + const Candidate triangle_and_bound = candidate_triangles.top(); + // Remove it from the candidate set as it will be processed now. + candidate_triangles.pop(); + + // Only process the triangle if it can contribute to the Hausdorff distance, + // i.e. if its upper bound is higher than the currently known best lower bound + // and the difference between the bounds to be obtained is larger than the + // user given error. + const auto& triangle_bounds = triangle_and_bound.bounds; + + CGAL_assertion(triangle_bounds.lower >= FT(0)); + CGAL_assertion(triangle_bounds.upper >= FT(0)); + CGAL_assertion(triangle_bounds.upper >= triangle_bounds.lower); + + if ( + (triangle_bounds.upper > global_bounds.lower) && + (triangle_bounds.upper - triangle_bounds.lower > error_bound)) { + + // Get the triangle that is to be subdivided and read its vertices. + const Triangle_3& triangle_for_subdivision = triangle_and_bound.triangle; + const Point_3 v0 = triangle_for_subdivision.vertex(0); + const Point_3 v1 = triangle_for_subdivision.vertex(1); + const Point_3 v2 = triangle_for_subdivision.vertex(2); + + // Check second stopping condition: All three vertices of the triangle + // are projected onto the same triangle in TM2. + const auto closest_triangle_v0 = tm2_tree.closest_point_and_primitive(v0); + const auto closest_triangle_v1 = tm2_tree.closest_point_and_primitive(v1); + const auto closest_triangle_v2 = tm2_tree.closest_point_and_primitive(v2); + CGAL_assertion(closest_triangle_v0.second != boost::graph_traits::null_face()); + CGAL_assertion(closest_triangle_v1.second != boost::graph_traits::null_face()); + CGAL_assertion(closest_triangle_v2.second != boost::graph_traits::null_face()); + if ( + (closest_triangle_v0.second == closest_triangle_v1.second) && + (closest_triangle_v1.second == closest_triangle_v2.second)) { + + // The upper bound of this triangle is the actual Hausdorff distance of + // the triangle to the second mesh. Use it as new global lower bound. + // Here, we update the reference to the realizing triangle as this is the best current guess. + global_bounds.lower = triangle_bounds.upper; + global_bounds.lpair.second = triangle_bounds.tm2_uface; + continue; + } + + // Check third stopping condition: All edge lengths of the triangle are + // smaller than the given error bound, we cannot get results beyond this bound. + if ( + CGAL::squared_distance(v0, v1) < squared_error_bound && + CGAL::squared_distance(v0, v2) < squared_error_bound && + CGAL::squared_distance(v1, v2) < squared_error_bound) { + + // The upper bound of this triangle is within error tolerance of + // the actual upper bound, use it. + global_bounds.lower = triangle_bounds.upper; + global_bounds.lpair.second = triangle_bounds.tm2_uface; + continue; + } + + // Subdivide the triangle into four smaller triangles. + const Point_3 v01 = CGAL::midpoint(v0, v1); + const Point_3 v02 = CGAL::midpoint(v0, v2); + const Point_3 v12 = CGAL::midpoint(v1, v2); + const std::array sub_triangles = { + Triangle_3(v0, v01, v02), Triangle_3(v1 , v01, v12), + Triangle_3(v2, v02, v12), Triangle_3(v01, v02, v12) }; + + // Send each of the four triangles to culling on B with the bounds of the parent triangle. + for (std::size_t i = 0; i < 4; ++i) { + + // Call culling on B with the single triangle found. + TM2_hd_traits traversal_traits_tm2( + tm2_tree.traits(), tm2, vpm2, + triangle_bounds, + infinity_value, + infinity_value, + infinity_value); + tm2_tree.traversal_with_priority(sub_triangles[i], traversal_traits_tm2); + + // Update global lower Hausdorff bound according to the obtained local bounds. + const auto local_bounds = traversal_traits_tm2.get_local_bounds(); + + CGAL_assertion(local_bounds.lower >= FT(0)); + CGAL_assertion(local_bounds.upper >= FT(0)); + CGAL_assertion(local_bounds.upper >= local_bounds.lower); + + CGAL_assertion(local_bounds.lpair == local_bounds.default_face_pair()); + CGAL_assertion(local_bounds.upair == local_bounds.default_face_pair()); + + if (local_bounds.lower > global_bounds.lower) { + global_bounds.lower = local_bounds.lower; + global_bounds.lpair.second = local_bounds.tm2_lface; + } + + // Add the subtriangle to the candidate list. + candidate_triangles.push( + Candidate(sub_triangles[i], local_bounds, triangle_and_bound.tm1_face)); + } + + // Update global upper Hausdorff bound after subdivision. + const FT current_max = candidate_triangles.top().bounds.upper; + CGAL_assertion(current_max >= FT(0)); + + if (current_max > global_bounds.lower) { + global_bounds.upper = current_max; + global_bounds.upair.second = candidate_triangles.top().bounds.tm2_uface; + } else { + global_bounds.upper = global_bounds.lower; + global_bounds.upair.second = global_bounds.lpair.second; + } + } + } + + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + std::cout << "* subdivision (sec.): " << timer.time() << std::endl; + #endif + + // Compute linear interpolation between the found lower and upper bounds. + CGAL_assertion(global_bounds.lower >= FT(0)); + CGAL_assertion(global_bounds.upper >= FT(0)); + CGAL_assertion(global_bounds.upper >= global_bounds.lower); + const double hdist = CGAL::to_double((global_bounds.lower + global_bounds.upper) / FT(2)); + + // Get realizing triangles. + CGAL_assertion(global_bounds.lpair.first != boost::graph_traits::null_face()); + CGAL_assertion(global_bounds.lpair.second != boost::graph_traits::null_face()); + CGAL_assertion(global_bounds.upair.first != boost::graph_traits::null_face()); + CGAL_assertion(global_bounds.upair.second != boost::graph_traits::null_face()); + + // Output face pairs, which realize the Hausdorff distance. + *out++ = global_bounds.lpair; + *out++ = global_bounds.upair; + + return hdist; } -} // end of namespace CGAL::Polygon_mesh_processing + +#if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) + +template +struct Triangle_mesh_wrapper { + + const TriangleMesh& tm; const VPM& vpm; + const bool is_tm2; TMTree& tm_tree; + Triangle_mesh_wrapper( + const TriangleMesh& tm, const VPM& vpm, + const bool is_tm2, TMTree& tm_tree) : + tm(tm), vpm(vpm), is_tm2(is_tm2), tm_tree(tm_tree) { } + + void build_tree() { + tm_tree.insert(faces(tm).begin(), faces(tm).end(), tm, vpm); + tm_tree.build(); + if (is_tm2) tm_tree.accelerate_distance_queries(); + else tm_tree.do_not_accelerate_distance_queries(); + } +}; + +template +struct Bounded_error_preprocessing { + + #ifdef CGAL_HAUSDORFF_DEBUG + using Timer = CGAL::Real_timer; + #endif + std::vector& tm_wrappers; + + // Constructor. + Bounded_error_preprocessing( + std::vector& tm_wrappers) : + tm_wrappers(tm_wrappers) { } + + // Split constructor. + Bounded_error_preprocessing( + Bounded_error_preprocessing& s, tbb::split) : + tm_wrappers(s.tm_wrappers) { } + + bool is_tm1_wrapper(const boost::any& operand) const { + return operand.type() == typeid(TM1Wrapper); + } + + bool is_tm2_wrapper(const boost::any& operand) const { + return operand.type() == typeid(TM2Wrapper); + } + + // TODO: make AABB tree build parallel! + void operator()(const tbb::blocked_range& range) { + + #ifdef CGAL_HAUSDORFF_DEBUG + Timer timer; + timer.reset(); + timer.start(); + std::cout.precision(17); + #endif + + for (std::size_t i = range.begin(); i != range.end(); ++i) { + CGAL_assertion(i < tm_wrappers.size()); + auto& tm_wrapper = tm_wrappers[i]; + if (is_tm1_wrapper(tm_wrapper)) { + TM1Wrapper& object = boost::any_cast(tm_wrapper); + object.build_tree(); + } else if (is_tm2_wrapper(tm_wrapper)) { + TM2Wrapper& object = boost::any_cast(tm_wrapper); + object.build_tree(); + } else { + CGAL_assertion_msg(false, "Error: wrong boost any type!"); + } + } + + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + std::cout << "* time operator() preprocessing (sec.): " << timer.time() << std::endl; + #endif + } + + void join(Bounded_error_preprocessing&) { } +}; + +template< class TriangleMesh1, + class TriangleMesh2, + class VPM1, + class VPM2, + class TM1Tree, + class TM2Tree, + class Kernel > +struct Bounded_error_distance_computation { + + using FT = typename Kernel::FT; + #ifdef CGAL_HAUSDORFF_DEBUG + using Timer = CGAL::Real_timer; + #endif + + const std::vector& tm1_parts; const TriangleMesh2& tm2; + const FT error_bound; const VPM1& vpm1; const VPM2& vpm2; + const FT infinity_value; const FT initial_bound; + const std::vector& tm1_trees; const TM2Tree& tm2_tree; + double distance; + + // Constructor. + Bounded_error_distance_computation( + const std::vector& tm1_parts, const TriangleMesh2& tm2, + const FT error_bound, const VPM1& vpm1, const VPM2& vpm2, + const FT infinity_value, const FT initial_bound, + const std::vector& tm1_trees, const TM2Tree& tm2_tree) : + tm1_parts(tm1_parts), tm2(tm2), + error_bound(error_bound), vpm1(vpm1), vpm2(vpm2), + infinity_value(infinity_value), initial_bound(initial_bound), + tm1_trees(tm1_trees), tm2_tree(tm2_tree), distance(-1.0) { + CGAL_assertion(tm1_parts.size() == tm1_trees.size()); + } + + // Split constructor. + Bounded_error_distance_computation( + Bounded_error_distance_computation& s, tbb::split) : + tm1_parts(s.tm1_parts), tm2(s.tm2), + error_bound(s.error_bound), vpm1(s.vpm1), vpm2(s.vpm2), + infinity_value(s.infinity_value), initial_bound(s.initial_bound), + tm1_trees(s.tm1_trees), tm2_tree(s.tm2_tree), distance(-1.0) { + CGAL_assertion(tm1_parts.size() == tm1_trees.size()); + } + + void operator()(const tbb::blocked_range& range) { + + #ifdef CGAL_HAUSDORFF_DEBUG + Timer timer; + timer.reset(); + timer.start(); + std::cout.precision(17); + #endif + + double hdist = -1.0; + auto stub = CGAL::Emptyset_iterator(); + + for (std::size_t i = range.begin(); i != range.end(); ++i) { + CGAL_assertion(i < tm1_parts.size()); + CGAL_assertion(i < tm1_trees.size()); + const auto& tm1 = tm1_parts[i]; + const auto& tm1_tree = tm1_trees[i]; + // TODO: add distance_bound (now it is -FT(1)) in case we use parallel + // for checking if two meshes are close. + const double dist = bounded_error_Hausdorff_impl( + tm1, tm2, error_bound, vpm1, vpm2, + infinity_value, initial_bound, -FT(1), + tm1_tree, tm2_tree, stub); + if (dist > hdist) hdist = dist; + } + if (hdist > distance) distance = hdist; + + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + std::cout << "* time operator() computation (sec.): " << timer.time() << std::endl; + #endif + } + + void join(Bounded_error_distance_computation& rhs) { + distance = (CGAL::max)(rhs.distance, distance); + } +}; + +#endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) + +template< class Concurrency_tag, + class Kernel, + class TriangleMesh1, + class TriangleMesh2, + class VPM1, + class VPM2, + class NamedParameters1, + class NamedParameters2, + class OutputIterator > +double bounded_error_one_sided_Hausdorff_impl( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const typename Kernel::FT error_bound, + const typename Kernel::FT distance_bound, + const bool compare_meshes, + const VPM1& vpm1, + const VPM2& vpm2, + const NamedParameters1& np1, + const NamedParameters2& np2, + OutputIterator& out) +{ + #if !defined(CGAL_LINKED_WITH_TBB) || !defined(CGAL_METIS_ENABLED) + CGAL_static_assertion_msg( + !(boost::is_convertible::value), + "Parallel_tag is enabled but at least TBB or METIS is unavailable."); + #endif + + using FT = typename Kernel::FT; + + using TM1 = TriangleMesh1; + using TM2 = TriangleMesh2; + + using TM1_primitive = AABB_face_graph_triangle_primitive; + using TM2_primitive = AABB_face_graph_triangle_primitive; + + using TM1_traits = AABB_traits; + using TM2_traits = AABB_traits; + + using TM1_tree = AABB_tree; + using TM2_tree = AABB_tree; + + using Face_handle_1 = typename boost::graph_traits::face_descriptor; + using Face_handle_2 = typename boost::graph_traits::face_descriptor; + + // This is parallel version: we split the tm1 into parts, build trees for all parts, and + // run in parallel all BHD computations. The final distance is obtained by taking the max + // between BHDs computed for these parts with respect to tm2. + // This is off by default because the parallel version does not show much of runtime improvement. + // The slowest part is building AABB trees and this is what should be accelerated in the future. + #if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) + using Point_3 = typename Kernel::Point_3; + using TMF = CGAL::Face_filtered_graph; + using TMF_primitive = AABB_face_graph_triangle_primitive; + using TMF_traits = AABB_traits; + using TMF_tree = AABB_tree; + using TM1_wrapper = Triangle_mesh_wrapper; + using TM2_wrapper = Triangle_mesh_wrapper; + + std::vector tm1_parts; + std::vector tm1_trees; + std::vector tm_wrappers; + #endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) + + #ifdef CGAL_HAUSDORFF_DEBUG + using Timer = CGAL::Real_timer; + Timer timer; + std::cout.precision(17); + #endif + + TM1_tree tm1_tree; + TM2_tree tm2_tree; + FT infinity_value = -FT(1); + + #if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) + // TODO: add to NP! + const int nb_cores = 4; + const std::size_t min_nb_faces_to_split = 100; // TODO: increase this number? + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "* num cores: " << nb_cores << std::endl; + #endif + + if ( + boost::is_convertible::value && + nb_cores > 1 && faces(tm1).size() >= min_nb_faces_to_split) { + + // (0) -- Compute infinity value. + #ifdef CGAL_HAUSDORFF_DEBUG + timer.reset(); + timer.start(); + #endif + const auto bbox1 = bbox(tm1); + const auto bbox2 = bbox(tm2); + const auto bb = bbox1 + bbox2; + const FT sq_dist = CGAL::squared_distance( + Point_3(bb.xmin(), bb.ymin(), bb.zmin()), + Point_3(bb.xmax(), bb.ymax(), bb.zmax())); + infinity_value = CGAL::approximate_sqrt(sq_dist) * FT(2); + CGAL_assertion(infinity_value >= FT(0)); + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + const double time0 = timer.time(); + std::cout << "- computing infinity (sec.): " << time0 << std::endl; + #endif + + // (1) -- Create partition of tm1. + #ifdef CGAL_HAUSDORFF_DEBUG + timer.reset(); + timer.start(); + #endif + using Face_property_tag = CGAL::dynamic_face_property_t; + auto face_pid_map = get(Face_property_tag(), tm1); + CGAL::METIS::partition_graph( + tm1, nb_cores, CGAL::parameters:: + face_partition_id_map(face_pid_map)); + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + const double time1 = timer.time(); + std::cout << "- computing partition time (sec.): " << time1 << std::endl; + #endif + + // (2) -- Create a filtered face graph for each part. + #ifdef CGAL_HAUSDORFF_DEBUG + timer.reset(); + timer.start(); + #endif + tm1_parts.reserve(nb_cores); + for (int i = 0; i < nb_cores; ++i) { + tm1_parts.emplace_back(tm1, i, face_pid_map); + // TODO: why is it triggered sometimes? + // CGAL_assertion(tm1_parts.back().is_selection_valid()); + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "- part " << i << " size: " << tm1_parts.back().number_of_faces() << std::endl; + #endif + } + CGAL_assertion(tm1_parts.size() == nb_cores); + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + const double time2 = timer.time(); + std::cout << "- creating graphs time (sec.): " << time2 << std::endl; + #endif + + // (3) -- Preprocess all input data. + #ifdef CGAL_HAUSDORFF_DEBUG + timer.reset(); + timer.start(); + #endif + tm1_trees.resize(tm1_parts.size()); + tm_wrappers.reserve(tm1_parts.size() + 1); + for (std::size_t i = 0; i < tm1_parts.size(); ++i) { + tm_wrappers.push_back(TM1_wrapper(tm1_parts[i], vpm1, false, tm1_trees[i])); + } + tm_wrappers.push_back(TM2_wrapper(tm2, vpm2, true, tm2_tree)); + CGAL_assertion(tm_wrappers.size() == tm1_parts.size() + 1); + Bounded_error_preprocessing bep(tm_wrappers); + tbb::parallel_reduce(tbb::blocked_range(0, tm_wrappers.size()), bep); + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + const double time3 = timer.time(); + std::cout << "- creating trees time (sec.) " << time3 << std::endl; + #endif + + // Final timing. + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "* preprocessing parallel time (sec.) " << + time0 + time1 + time2 + time3 << std::endl; + #endif + + } else // sequential version + #endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) + { + #ifdef CGAL_HAUSDORFF_DEBUG + timer.reset(); + timer.start(); + std::cout << "* preprocessing sequential version " << std::endl; + #endif + bool rebuild = false; + std::vector tm1_only; + std::vector tm2_only; + std::tie(infinity_value, rebuild) = preprocess_bounded_error_Hausdorff_impl( + tm1, tm2, compare_meshes, vpm1, vpm2, true, np1, np2, + tm1_tree, tm2_tree, tm1_only, tm2_only); + CGAL_assertion(!rebuild); + if (infinity_value >= FT(0)) { + tm1_tree.build(); + tm2_tree.build(); + tm1_tree.do_not_accelerate_distance_queries(); + tm2_tree.accelerate_distance_queries(); + } + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + std::cout << "* preprocessing sequential time (sec.) " << timer.time() << std::endl; + #endif + } + + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "* infinity_value: " << infinity_value << std::endl; + #endif + if (infinity_value < FT(0)) { + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "* culling rate: 100%" << std::endl; + #endif + const auto face1 = *(faces(tm1).begin()); + const auto face2 = *(faces(tm2).begin()); + *out++ = std::make_pair(face1, face2); + *out++ = std::make_pair(face1, face2); + return 0.0; // TM1 is part of TM2 so the distance is zero + } + CGAL_assertion(error_bound >= FT(0)); + CGAL_assertion(infinity_value > FT(0)); + const FT initial_bound = error_bound; + std::atomic hdist; + + #ifdef CGAL_HAUSDORFF_DEBUG + timer.reset(); + timer.start(); + #endif + + #if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) + if ( + boost::is_convertible::value && + nb_cores > 1 && faces(tm1).size() >= min_nb_faces_to_split) { + + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "* executing parallel version " << std::endl; + #endif + Bounded_error_distance_computation bedc( + tm1_parts, tm2, error_bound, vpm1, vpm2, + infinity_value, initial_bound, tm1_trees, tm2_tree); + tbb::parallel_reduce(tbb::blocked_range(0, tm1_parts.size()), bedc); + hdist = bedc.distance; + + } else // sequential version + #endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) + { + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "* executing sequential version " << std::endl; + #endif + hdist = bounded_error_Hausdorff_impl( + tm1, tm2, error_bound, vpm1, vpm2, + infinity_value, initial_bound, distance_bound, + tm1_tree, tm2_tree, out); + } + + #ifdef CGAL_HAUSDORFF_DEBUG + timer.stop(); + std::cout << "* computation time (sec.) " << timer.time() << std::endl; + #endif + + CGAL_assertion(hdist >= 0.0); + return hdist; +} + +template< class Concurrency_tag, + class Kernel, + class TriangleMesh1, + class TriangleMesh2, + class VPM1, + class VPM2, + class NamedParameters1, + class NamedParameters2, + class OutputIterator1, + class OutputIterator2 > +double bounded_error_symmetric_Hausdorff_impl( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const typename Kernel::FT error_bound, + const typename Kernel::FT distance_bound, + const bool compare_meshes, + const VPM1& vpm1, + const VPM2& vpm2, + const NamedParameters1& np1, + const NamedParameters2& np2, + OutputIterator1& out1, + OutputIterator2& out2) +{ + #if !defined(CGAL_LINKED_WITH_TBB) || !defined(CGAL_METIS_ENABLED) + CGAL_static_assertion_msg( + !(boost::is_convertible::value), + "Parallel_tag is enabled but at least TBB or METIS is unavailable."); + #endif + + // Optimized version. + // -- We compare meshes only if it is required. + // -- We first build trees and rebuild them only if it is required. + // -- We provide better initial lower bound in the second call to the Hausdorff distance. + using FT = typename Kernel::FT; + + using TM1_primitive = AABB_face_graph_triangle_primitive; + using TM2_primitive = AABB_face_graph_triangle_primitive; + + using TM1_traits = AABB_traits; + using TM2_traits = AABB_traits; + + using TM1_tree = AABB_tree; + using TM2_tree = AABB_tree; + + using Face_handle_1 = typename boost::graph_traits::face_descriptor; + using Face_handle_2 = typename boost::graph_traits::face_descriptor; + + std::vector tm1_only; + std::vector tm2_only; + + // All trees below are built and/or accelerated lazily. + TM1_tree tm1_tree; + TM2_tree tm2_tree; + FT infinity_value = -FT(1); + bool rebuild = false; + std::tie(infinity_value, rebuild) = preprocess_bounded_error_Hausdorff_impl( + tm1, tm2, compare_meshes, vpm1, vpm2, false, np1, np2, + tm1_tree, tm2_tree, tm1_only, tm2_only); + + if (infinity_value < FT(0)) { + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout.precision(17); + std::cout << "* culling rate: 100%" << std::endl; + #endif + const auto face1 = *(faces(tm1).begin()); + const auto face2 = *(faces(tm2).begin()); + *out1++ = std::make_pair(face1, face2); + *out1++ = std::make_pair(face1, face2); + *out2++ = std::make_pair(face2, face1); + *out2++ = std::make_pair(face2, face1); + return 0.0; // TM1 and TM2 are equal so the distance is zero + } + CGAL_assertion(infinity_value > FT(0)); + + // Compute the first one-sided distance. + FT initial_bound = error_bound; + double dista = CGAL::to_double(error_bound); + + if (!compare_meshes || (compare_meshes && tm1_only.size() > 0)) { + dista = bounded_error_Hausdorff_impl( + tm1, tm2, error_bound, vpm1, vpm2, + infinity_value, initial_bound, distance_bound, + tm1_tree, tm2_tree, out1); + } + + // In case this is true, we need to rebuild trees in order to accelerate + // computations for the second call. + if (rebuild) { + CGAL_assertion(compare_meshes); + tm1_tree.clear(); + tm2_tree.clear(); + CGAL_assertion(tm2_only.size() > 0); + CGAL_assertion(tm2_only.size() < faces(tm2).size()); + tm1_tree.insert(faces(tm1).begin(), faces(tm1).end(), tm1, vpm1); + tm2_tree.insert(tm2_only.begin(), tm2_only.end(), tm2, vpm2); + } + + // Compute the second one-sided distance. + initial_bound = static_cast(dista); // TODO: we should better test this optimization! + double distb = CGAL::to_double(error_bound); + + if (!compare_meshes || (compare_meshes && tm2_only.size() > 0)) { + distb = bounded_error_Hausdorff_impl( + tm2, tm1, error_bound, vpm2, vpm1, + infinity_value, initial_bound, distance_bound, + tm2_tree, tm1_tree, out2); + } + + // Return the maximum. + return (CGAL::max)(dista, distb); +} + +template +typename Kernel::FT recursive_hausdorff_subdivision( + const typename Kernel::Point_3& v0, + const typename Kernel::Point_3& v1, + const typename Kernel::Point_3& v2, + const TM2_tree& tm2_tree, + const typename Kernel::FT squared_error_bound) +{ + // If all edge lengths of the triangle are below the error_bound, + // return maximum of the distances of the three points to TM2 (via TM2_tree). + const auto max_squared_edge_length = + (CGAL::max)( + (CGAL::max)( + CGAL::squared_distance(v0, v1), + CGAL::squared_distance(v0, v2)), + CGAL::squared_distance(v1, v2)); + + if (max_squared_edge_length < squared_error_bound) { + return (CGAL::max)( + (CGAL::max)( + CGAL::squared_distance(v0, tm2_tree.closest_point(v0)), + CGAL::squared_distance(v1, tm2_tree.closest_point(v1))), + CGAL::squared_distance(v2, tm2_tree.closest_point(v2))); + } + + // Else subdivide the triangle and proceed recursively. + const auto v01 = midpoint(v0, v1); + const auto v02 = midpoint(v0, v2); + const auto v12 = midpoint(v1, v2); + + return (CGAL::max)( + (CGAL::max)( + recursive_hausdorff_subdivision(v0, v01, v02, tm2_tree, squared_error_bound), + recursive_hausdorff_subdivision(v1, v01, v12, tm2_tree, squared_error_bound)), + (CGAL::max)( + recursive_hausdorff_subdivision(v2 , v02, v12, tm2_tree, squared_error_bound), + recursive_hausdorff_subdivision(v01, v02, v12, tm2_tree, squared_error_bound))); +} + +template< class Concurrency_tag, + class Kernel, + class TriangleMesh1, + class TriangleMesh2, + class VPM1, + class VPM2 > +double bounded_error_Hausdorff_naive_impl( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const typename Kernel::FT error_bound, + const VPM1& vpm1, + const VPM2& vpm2) +{ + using FT = typename Kernel::FT; + using Point_3 = typename Kernel::Point_3; + using Triangle_3 = typename Kernel::Triangle_3; + + using TM2_primitive = AABB_face_graph_triangle_primitive; + using TM2_traits = AABB_traits; + using TM2_tree = AABB_tree; + + using TM1_face_to_triangle_map = Triangle_from_face_descriptor_map; + + // Initially, no lower bound is known. + FT squared_lower_bound = FT(0); + + // Work with squares in the following, only draw sqrt at the very end. + const FT squared_error_bound = error_bound * error_bound; + + // Build an AABB tree on tm2. + TM2_tree tm2_tree(faces(tm2).begin(), faces(tm2).end(), tm2, vpm2); + tm2_tree.build(); + tm2_tree.accelerate_distance_queries(); + + // Build a map to obtain actual triangles from the face descriptors of tm1. + const TM1_face_to_triangle_map face_to_triangle_map(&tm1, vpm1); + + // Iterate over the faces of TM1. + for (const auto& face : faces(tm1)) { + + // Get the vertices of the face and pass them on to a recursive method. + const Triangle_3 triangle = get(face_to_triangle_map, face); + const Point_3 v0 = triangle.vertex(0); + const Point_3 v1 = triangle.vertex(1); + const Point_3 v2 = triangle.vertex(2); + + // Recursively process the current triangle to obtain a lower bound on its Hausdorff distance. + const FT triangle_bound = recursive_hausdorff_subdivision( + v0, v1, v2, tm2_tree, squared_error_bound); + + // Store the largest lower bound. + if (triangle_bound > squared_lower_bound) { + squared_lower_bound = triangle_bound; + } + } + + // Return linear interpolation between found upper and lower bound. + return CGAL::sqrt(CGAL::to_double(squared_lower_bound)); +} + +} // end of namespace internal + +/** + * \ingroup PMP_distance_grp + * returns an estimate on the Hausdorff distance between `tm1` and `tm2` that + * is at most `error_bound` away from the actual Hausdorff distance between + * the two given meshes. + * + * @tparam Concurrency_tag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag` and `Parallel_tag`. + * Currently, the parallel version is not implemented and the + * sequential version is always used whatever tag is chosen! + * + * @tparam TriangleMesh1 a model of the concept `FaceListGraph` + * @tparam TriangleMesh2 a model of the concept `FaceListGraph` + * + * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" + * + * @param tm1 a triangle mesh + * @param tm2 another triangle mesh + * + * @param error_bound a maximum bound by which the Hausdorff distance estimate is + * allowed to deviate from the actual Hausdorff distance. + * + * @param np1 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below + * @param np2 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamNBegin{vertex_point_map} + * \cgalParamDescription{a property map associating points to the vertices of `tm1` and `tm2` (`np1` and `np2`, respectively)} + * \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits::%vertex_descriptor` + * as key type and `%Point_3` as value type} + * \cgalParamDefault{`boost::get(CGAL::vertex_point, tm)`} + * \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` + * must be available in `TriangleMeshX`.} + * \cgalParamNEnd + * \cgalParamNBegin{match_faces} + * \cgalParamDescription{a boolean tag that turns on the preprocessing step that filters out all faces, + * which belong to both meshes and hence do not contribute to the final distance} + * \cgalParamType{Boolean} + * \cgalParamDefault{true} + * \cgalParamExtra{Both `np1` and `np2` must have this tag true in order to activate this preprocessing.} + * \cgalParamNEnd + * \cgalNamedParamsEnd + * + * @return the one-sided Hausdorff distance + */ +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2, + class NamedParameters1, + class NamedParameters2 > +double bounded_error_Hausdorff_distance( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + CGAL_assertion_code( + const bool is_triangle = is_triangle_mesh(tm1) && is_triangle_mesh(tm2)); + CGAL_assertion_msg(is_triangle, + "Both meshes must be triangulated to compute this distance!"); + + using Traits = typename GetGeomTraits::type; + using FT = typename Traits::FT; + + const auto vpm1 = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::vertex_point), + get_const_property_map(vertex_point, tm1)); + const auto vpm2 = parameters::choose_parameter( + parameters::get_parameter(np2, internal_np::vertex_point), + get_const_property_map(vertex_point, tm2)); + + const bool match_faces1 = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::match_faces), true); + const bool match_faces2 = parameters::choose_parameter( + parameters::get_parameter(np2, internal_np::match_faces), true); + const bool match_faces = match_faces1 && match_faces2; + + auto out = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::output_iterator), + CGAL::Emptyset_iterator()); + + CGAL_precondition(error_bound >= 0.0); + const FT error_threshold = static_cast(error_bound); + return internal::bounded_error_one_sided_Hausdorff_impl( + tm1, tm2, error_threshold, -FT(1), match_faces, vpm1, vpm2, np1, np2, out); +} + +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2, + class NamedParameters1 > +double bounded_error_Hausdorff_distance( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound, + const NamedParameters1& np1) +{ + return bounded_error_Hausdorff_distance( + tm1, tm2, error_bound, np1, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2 > +double bounded_error_Hausdorff_distance( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound = 0.0001) +{ + return bounded_error_Hausdorff_distance( + tm1, tm2, error_bound, parameters::all_default()); +} + +/** + * \ingroup PMP_distance_grp + * returns the maximum of `bounded_error_Hausdorff_distance(tm1, tm2, error_bound, np1, np2)` + * and `bounded_error_Hausdorff_distance(tm2, tm1, error_bound, np2, np1)`. + * + * This function optimizes all internal calls to shared data structures in order to + * speed up the computation. + * + * @return the symmetric Hausdorff distance + * @see `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` + */ +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2, + class NamedParameters1, + class NamedParameters2 > +double bounded_error_symmetric_Hausdorff_distance( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + CGAL_assertion_code( + const bool is_triangle = is_triangle_mesh(tm1) && is_triangle_mesh(tm2)); + CGAL_assertion_msg(is_triangle, + "Both meshes must be triangulated to compute this distance!"); + + using Traits = typename GetGeomTraits::type; + using FT = typename Traits::FT; + + const auto vpm1 = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::vertex_point), + get_const_property_map(vertex_point, tm1)); + const auto vpm2 = parameters::choose_parameter( + parameters::get_parameter(np2, internal_np::vertex_point), + get_const_property_map(vertex_point, tm2)); + + const bool match_faces1 = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::match_faces), true); + const bool match_faces2 = parameters::choose_parameter( + parameters::get_parameter(np2, internal_np::match_faces), true); + const bool match_faces = match_faces1 && match_faces2; + + // TODO: should we return a union of these realizing triangles? + auto out1 = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::output_iterator), + CGAL::Emptyset_iterator()); + auto out2 = parameters::choose_parameter( + parameters::get_parameter(np2, internal_np::output_iterator), + CGAL::Emptyset_iterator()); + + CGAL_precondition(error_bound >= 0.0); + const FT error_threshold = static_cast(error_bound); + return internal::bounded_error_symmetric_Hausdorff_impl( + tm1, tm2, error_threshold, -FT(1), match_faces, vpm1, vpm2, np1, np2, out1, out2); +} + +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2, + class NamedParameters1 > +double bounded_error_symmetric_Hausdorff_distance( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound, + const NamedParameters1& np1) +{ + return bounded_error_symmetric_Hausdorff_distance( + tm1, tm2, error_bound, np1, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2> +double bounded_error_symmetric_Hausdorff_distance( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound = 0.0001) +{ + return bounded_error_symmetric_Hausdorff_distance( + tm1, tm2, error_bound, parameters::all_default()); +} + +// TODO: Find better name! +// TODO: Should we use one-sided or symmetric distance here? + +/** + * \ingroup PMP_distance_grp + * returns `true` if the Hausdorff distance between two meshes is larger than + * the user-defined max distance, otherwise it returns `false`. The distance used + * to compute the proximity of the meshes is the bounded-error Hausdorff distance. + * + * Instead of computing the full distance and checking it against the user-provided + * value, this function early quits in case certain criteria show that the meshes + * do not satisfy the provided `distance_bound`. + * + * \cgalNamedParamsBegin + * \cgalParamNBegin{use_one_sided_hausdorff} + * \cgalParamDescription{a boolean tag indicating if the one-sided Hausdorff distance should be used.} + * \cgalParamType{Boolean} + * \cgalParamDefault{`true`} + * \cgalParamExtra{If this tag is set to `false`, the symmetric Hausdorff distance is used.} + * \cgalParamNEnd + * \cgalNamedParamsEnd + * + * @return Boolean `true` or `false` +* @see `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` + */ +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2, + class NamedParameters1, + class NamedParameters2 > +bool is_Hausdorff_distance_larger( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double distance_bound, + const double error_bound, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + CGAL_assertion_code( + const bool is_triangle = is_triangle_mesh(tm1) && is_triangle_mesh(tm2)); + CGAL_assertion_msg(is_triangle, + "Both meshes must be triangulated in order to be compared!"); + + using Traits = typename GetGeomTraits::type; + using FT = typename Traits::FT; + + const auto vpm1 = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::vertex_point), + get_const_property_map(vertex_point, tm1)); + const auto vpm2 = parameters::choose_parameter( + parameters::get_parameter(np2, internal_np::vertex_point), + get_const_property_map(vertex_point, tm2)); + + const bool match_faces1 = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::match_faces), true); + const bool match_faces2 = parameters::choose_parameter( + parameters::get_parameter(np2, internal_np::match_faces), true); + const bool match_faces = match_faces1 && match_faces2; + + const bool use_one_sided = parameters::choose_parameter( + parameters::get_parameter(np1, internal_np::use_one_sided_hausdorff), true); + CGAL_precondition(error_bound >= 0.0); + const FT error_threshold = static_cast(error_bound); + CGAL_precondition(distance_bound >= 0.0); + const FT distance_threshold = static_cast(distance_bound); + auto stub = CGAL::Emptyset_iterator(); + + double hdist = -1.0; + if (use_one_sided) { + hdist = internal::bounded_error_one_sided_Hausdorff_impl( + tm1, tm2, error_threshold, distance_threshold, match_faces, vpm1, vpm2, np1, np2, stub); + } else { + hdist = internal::bounded_error_symmetric_Hausdorff_impl( + tm1, tm2, error_threshold, distance_threshold, match_faces, vpm1, vpm2, np1, np2, stub, stub); + } + CGAL_assertion(hdist >= 0.0); + + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout.precision(17); + std::cout << "- fin distance: " << hdist << std::endl; + std::cout << "- max distance: " << distance_bound << std::endl; + #endif + return hdist > distance_bound; +} + +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2, + class NamedParameters1 > +double is_Hausdorff_distance_larger( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double max_distance, + const double error_bound, + const NamedParameters1& np1) +{ + return is_Hausdorff_distance_larger( + tm1, tm2, max_distance, error_bound, np1, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2 > +double is_Hausdorff_distance_larger( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double max_distance = 1.0, + const double error_bound = 0.0001) +{ + return is_Hausdorff_distance_larger( + tm1, tm2, max_distance, error_bound, parameters::all_default()); +} + +// Implementation of the naive Bounded Error Hausdorff distance. +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2, + class NamedParameters1, + class NamedParameters2 > +double bounded_error_Hausdorff_distance_naive( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + CGAL_assertion_code( + const bool is_triangle = is_triangle_mesh(tm1) && is_triangle_mesh(tm2)); + CGAL_assertion_msg(is_triangle, + "Both meshes must be triangulated to compute this distance!"); + + using Traits = typename GetGeomTraits::type; + using FT = typename Traits::FT; + + using parameters::choose_parameter; + using parameters::get_parameter; + + const auto vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point), + get_const_property_map(vertex_point, tm1)); + const auto vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point), + get_const_property_map(vertex_point, tm2)); + + CGAL_precondition(error_bound >= 0.0); + const FT error_threshold = static_cast(error_bound); + return internal::bounded_error_Hausdorff_naive_impl( + tm1, tm2, error_threshold, vpm1, vpm2); +} + +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2, + class NamedParameters1 > +double bounded_error_Hausdorff_distance_naive( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound, + const NamedParameters1& np1) +{ + return bounded_error_Hausdorff_distance_naive( + tm1, tm2, error_bound, np1, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh1, + class TriangleMesh2 > +double bounded_error_Hausdorff_distance_naive( + const TriangleMesh1& tm1, + const TriangleMesh2& tm2, + const double error_bound = 0.0001) +{ + return bounded_error_Hausdorff_distance_naive( + tm1, tm2, error_bound, parameters::all_default()); +} + +} } // end of namespace CGAL::Polygon_mesh_processing #endif //CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_Hausdorff_distance.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_Hausdorff_distance.h new file mode 100644 index 000000000000..a7f5ddd01eaa --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_Hausdorff_distance.h @@ -0,0 +1,588 @@ +// Copyright (c) 2019 GeometryFactory SARL (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// +// Author(s) : Sebastien Loriot, Martin Skrodzki, Dmitry Anisimov + +#ifndef CGAL_PMP_INTERNAL_AABB_TRAVERSAL_TRAITS_WITH_HAUSDORFF_DISTANCE +#define CGAL_PMP_INTERNAL_AABB_TRAVERSAL_TRAITS_WITH_HAUSDORFF_DISTANCE + +#include + +// STL includes. +#include +#include + +// CGAL includes. +#include +#include +#include + +namespace CGAL { + + // Bounds. + template< class Kernel, + class Face_handle_1, + class Face_handle_2> + struct Bounds { + using FT = typename Kernel::FT; + + FT m_infinity_value = -FT(1); + Bounds(const FT infinity_value) : + m_infinity_value(infinity_value) { } + + FT lower = m_infinity_value; + FT upper = m_infinity_value; + // TODO: update + Face_handle_2 tm2_lface = Face_handle_2(); + Face_handle_2 tm2_uface = Face_handle_2(); + std::pair lpair = default_face_pair(); + std::pair upair = default_face_pair(); + + const std::pair default_face_pair() const { + return std::make_pair(Face_handle_1(), Face_handle_2()); + } + }; + + // Candidate triangle. + template< class Kernel, + class Face_handle_1, + class Face_handle_2> + struct Candidate_triangle { + using FT = typename Kernel::FT; + using Triangle_3 = typename Kernel::Triangle_3; + using Candidate_bounds = Bounds; + + Candidate_triangle( + const Triangle_3& triangle, const Candidate_bounds& bounds, const Face_handle_1& fh) : + triangle(triangle), bounds(bounds), tm1_face(fh) + { } + + Triangle_3 triangle; + Candidate_bounds bounds; + Face_handle_1 tm1_face; + // TODO: no need to use bounds.lower? + bool operator>(const Candidate_triangle& other) const { + CGAL_assertion(bounds.upper >= FT(0)); + CGAL_assertion(other.bounds.upper >= FT(0)); + return bounds.upper < other.bounds.upper; + } + bool operator<(const Candidate_triangle& other) const { + CGAL_assertion(bounds.upper >= FT(0)); + CGAL_assertion(other.bounds.upper >= FT(0)); + return bounds.upper > other.bounds.upper; + } + }; + + // Hausdorff primitive traits on TM2. + template< class AABBTraits, + class Query, + class Kernel, + class TriangleMesh1, + class TriangleMesh2, + class VPM2> + class Hausdorff_primitive_traits_tm2 + { + using FT = typename Kernel::FT; + using Point_3 = typename Kernel::Point_3; + using Vector_3 = typename Kernel::Vector_3; + using Triangle_3 = typename Kernel::Triangle_3; + + using Project_point_3 = typename Kernel::Construct_projected_point_3; + using Face_handle_1 = typename boost::graph_traits::face_descriptor; + using Face_handle_2 = typename boost::graph_traits::face_descriptor; + using Local_bounds = Bounds; + + using TM2_face_to_triangle_map = Triangle_from_face_descriptor_map; + + public: + using Priority = FT; + Hausdorff_primitive_traits_tm2( + const AABBTraits& traits, + const TriangleMesh2& tm2, const VPM2& vpm2, + const Local_bounds& local_bounds, + const FT h_v0_lower_init, + const FT h_v1_lower_init, + const FT h_v2_lower_init) : + m_traits(traits), m_tm2(tm2), m_vpm2(vpm2), + m_face_to_triangle_map(&m_tm2, m_vpm2), + h_local_bounds(local_bounds) { + + // Initialize the global and local bounds with the given values. + h_v0_lower = h_v0_lower_init; + h_v1_lower = h_v1_lower_init; + h_v2_lower = h_v2_lower_init; + } + + // Explore the whole tree, i.e. always enter children if the method + // do_intersect() below determines that it is worthwhile. + bool go_further() const { return true; } + + // Compute the explicit Hausdorff distance to the given primitive. + template + void intersection(const Query& query, const Primitive& primitive) { + + /* Have reached a single triangle, process it. + / Determine the distance according to + / min_{b \in primitive} ( max_{vertex in query} ( d(vertex, b) ) ) + / + / Here, we only have one triangle in B, i.e. tm2. Thus, it suffices to + / compute the distance of the vertices of the query triangles to the + / primitive triangle and use the maximum of the obtained distances. + */ + + // The query object is a triangle from TM1, get its vertices. + const Point_3 v0 = query.vertex(0); + const Point_3 v1 = query.vertex(1); + const Point_3 v2 = query.vertex(2); + + CGAL_assertion(primitive.id() != Face_handle_2()); + const Triangle_3 triangle = get(m_face_to_triangle_map, primitive.id()); + + // Compute distances of the vertices to the primitive triangle in TM2. + const FT v0_dist = CGAL::approximate_sqrt(CGAL::squared_distance(m_project_point(triangle, v0), v0)); + if (v0_dist < h_v0_lower) h_v0_lower = v0_dist; // it is () part of (11) in the paper + + const FT v1_dist = CGAL::approximate_sqrt(CGAL::squared_distance(m_project_point(triangle, v1), v1)); + if (v1_dist < h_v1_lower) h_v1_lower = v1_dist; // it is () part of (11) in the paper + + const FT v2_dist = CGAL::approximate_sqrt(CGAL::squared_distance(m_project_point(triangle, v2), v2)); + if (v2_dist < h_v2_lower) h_v2_lower = v2_dist; // it is () part of (11) in the paper + + // Get the distance as maximizers over all vertices. + const FT distance_lower = (CGAL::max)((CGAL::max)(h_v0_lower, h_v1_lower), h_v2_lower); // it is (11) in the paper + const FT distance_upper = (CGAL::max)((CGAL::max)(v0_dist, v1_dist), v2_dist); // it is () part of (10) in the paper + + CGAL_assertion(distance_lower >= FT(0)); + CGAL_assertion(distance_upper >= FT(0)); + CGAL_assertion(distance_upper >= distance_lower); + + // Since we are at the level of a single triangle in TM2, distance_upper is + // actually the correct Hausdorff distance from the query triangle in + // TM1 to the primitive triangle in TM2. + CGAL_assertion(h_local_bounds.lower >= FT(0)); + if (distance_lower < h_local_bounds.lower) { + h_local_bounds.lower = distance_lower; + h_local_bounds.tm2_lface = primitive.id(); + } + CGAL_assertion(h_local_bounds.upper >= FT(0)); + if (distance_upper < h_local_bounds.upper) { // it is (10) in the paper + h_local_bounds.upper = distance_upper; + h_local_bounds.tm2_uface = primitive.id(); + } + CGAL_assertion(h_local_bounds.upper >= h_local_bounds.lower); + } + + // Determine whether child nodes will still contribute to a smaller + // Hausdorff distance and thus have to be entered. + template + std::pair + do_intersect_with_priority(const Query& query, const Node& node) const { + + // Get the bounding box of the nodes. + const auto bbox = node.bbox(); + + // Get the vertices of the query triangle. + const Point_3 v0 = query.vertex(0); + const Point_3 v1 = query.vertex(1); + const Point_3 v2 = query.vertex(2); + + // Find the axis aligned bbox of the triangle. + const Point_3 tri_min = Point_3( + (CGAL::min)((CGAL::min)(v0.x(), v1.x()), v2.x()), + (CGAL::min)((CGAL::min)(v0.y(), v1.y()), v2.y()), + (CGAL::min)((CGAL::min)(v0.z(), v1.z()), v2.z())); + + const Point_3 tri_max = Point_3( + (CGAL::max)((CGAL::max)(v0.x(), v1.x()), v2.x()), + (CGAL::max)((CGAL::max)(v0.y(), v1.y()), v2.y()), + (CGAL::max)((CGAL::max)(v0.z(), v1.z()), v2.z())); + + // Compute distance of the bounding boxes. + // Distance along the x-axis. + FT dist_x = FT(0); + if (tri_max.x() < bbox.min(0)) { + dist_x = bbox.min(0) - tri_max.x(); + } else if (bbox.max(0) < tri_min.x()) { + dist_x = tri_min.x() - bbox.max(0); + } + + // Distance along the y-axis. + FT dist_y = FT(0); + if (tri_max.y() < bbox.min(1)) { + dist_y = bbox.min(1) - tri_max.y(); + } else if (bbox.max(1) < tri_min.y()) { + dist_y = tri_min.y() - bbox.max(1); + } + + // Distance along the z-axis. + FT dist_z = FT(0); + if (tri_max.z() < bbox.min(2)) { + dist_z = bbox.min(2) - tri_max.z(); + } else if (bbox.max(2) < tri_min.z()) { + dist_z = tri_min.z() - bbox.max(2); + } + + // Lower bound on the distance between the two bounding boxes is given + // as the length of the diagonal of the bounding box between them. + const FT dist = CGAL::approximate_sqrt(Vector_3(dist_x, dist_y, dist_z).squared_length()); + + // See Algorithm 2. + // Check whether investigating the bbox can still lower the Hausdorff + // distance and improve the current global bound. If so, enter the box. + CGAL_assertion(h_local_bounds.lower >= FT(0)); + if (dist <= h_local_bounds.lower) { + return std::make_pair(true , -dist); + } else { + return std::make_pair(false, FT(0)); + } + } + + template + bool do_intersect(const Query& query, const Node& node) const { + return this->do_intersect_with_priority(query, node).first; + } + + // Return the local Hausdorff bounds computed for the passed query triangle. + Local_bounds get_local_bounds() const { + return h_local_bounds; + } + + template + void traverse_group(const Query& query, PrimitiveConstIterator group_begin, PrimitiveConstIterator group_end) { + for (PrimitiveConstIterator it = group_begin; it != group_end; ++it) { + this->intersection(query, *it); + } + } + + private: + // Input data. + const AABBTraits& m_traits; + const TriangleMesh2& m_tm2; + const VPM2& m_vpm2; + const TM2_face_to_triangle_map m_face_to_triangle_map; + + // Local Hausdorff bounds for the query triangle. + Local_bounds h_local_bounds; + FT h_v0_lower, h_v1_lower, h_v2_lower; + Project_point_3 m_project_point; + }; + + // Hausdorff primitive traits on TM1. + template< class AABBTraits, + class Query, + class Kernel, + class TriangleMesh1, + class TriangleMesh2, + class VPM1, + class VPM2> + class Hausdorff_primitive_traits_tm1 + { + using FT = typename Kernel::FT; + using Point_3 = typename Kernel::Point_3; + using Vector_3 = typename Kernel::Vector_3; + using Triangle_3 = typename Kernel::Triangle_3; + + using TM2_primitive = AABB_face_graph_triangle_primitive; + using TM2_traits = AABB_traits; + using TM2_tree = AABB_tree; + using TM2_hd_traits = Hausdorff_primitive_traits_tm2; + + using TM1_face_to_triangle_map = Triangle_from_face_descriptor_map; + + using Face_handle_1 = typename boost::graph_traits::face_descriptor; + using Face_handle_2 = typename boost::graph_traits::face_descriptor; + + using Global_bounds = Bounds; + using Candidate = Candidate_triangle; + using Heap_type = std::priority_queue; + + public: + using Priority = FT; + Hausdorff_primitive_traits_tm1( + const AABBTraits& traits, const TM2_tree& tree, + const TriangleMesh1& tm1, const TriangleMesh2& tm2, + const VPM1& vpm1, const VPM2& vpm2, + const FT error_bound, + const FT infinity_value, + const FT initial_bound, + const FT distance_bound) : + m_traits(traits), + m_tm1(tm1), m_tm2(tm2), + m_vpm1(vpm1), m_vpm2(vpm2), + m_tm2_tree(tree), + m_face_to_triangle_map(&m_tm1, m_vpm1), + m_error_bound(error_bound), + m_infinity_value(infinity_value), + m_initial_bound(initial_bound), + m_distance_bound(distance_bound), + h_global_bounds(m_infinity_value), + m_early_quit(false) { + + CGAL_precondition(m_error_bound >= FT(0)); + CGAL_precondition(m_infinity_value >= FT(0)); + CGAL_precondition(m_initial_bound >= m_error_bound); + + // Initialize the global bounds with 0, they will only grow. + // If we leave zero here, then we are very slow even for big input error bounds! + // Instead, we can use m_error_bound as our initial guess to filter out all pairs, + // which are already within this bound. It makes the code faster for close meshes. + // We also use initial_lower_bound here to accelerate the symmetric distance computation. + h_global_bounds.lower = m_initial_bound; // = FT(0); + h_global_bounds.upper = m_initial_bound; // = FT(0); + } + + // Explore the whole tree, i.e. always enter children if the methods + // do_intersect() below determine that it is worthwhile. + bool go_further() const { + return !m_early_quit; + } + + // Compute the explicit Hausdorff distance to the given primitive. + template + void intersection(const Query&, const Primitive& primitive) { + + if (m_early_quit) return; + + // Set initial tight bounds. + CGAL_assertion(primitive.id() != Face_handle_1()); + std::pair fpair; + const FT max_dist = get_maximum_distance(primitive.id(), fpair); + CGAL_assertion(max_dist >= FT(0)); + CGAL_assertion(fpair.first == primitive.id()); + + Bounds initial_bounds(m_infinity_value); + initial_bounds.lower = max_dist + m_error_bound; + initial_bounds.upper = max_dist + m_error_bound; + initial_bounds.tm2_lface = fpair.second; + initial_bounds.tm2_uface = fpair.second; + + // Call Culling on B with the single triangle found. + TM2_hd_traits traversal_traits_tm2( + m_tm2_tree.traits(), m_tm2, m_vpm2, + initial_bounds, // tighter bounds, in the paper, they start from infinity, see below + // Bounds(m_infinity_value), // starting from infinity + m_infinity_value, + m_infinity_value, + m_infinity_value); + + const Triangle_3 triangle = get(m_face_to_triangle_map, fpair.first); + m_tm2_tree.traversal_with_priority(triangle, traversal_traits_tm2); + + // Update global Hausdorff bounds according to the obtained local bounds. + const auto local_bounds = traversal_traits_tm2.get_local_bounds(); + + CGAL_assertion(local_bounds.lower >= FT(0)); + CGAL_assertion(local_bounds.upper >= FT(0)); + CGAL_assertion(local_bounds.upper >= local_bounds.lower); + CGAL_assertion(local_bounds.lpair == initial_bounds.default_face_pair()); + CGAL_assertion(local_bounds.upair == initial_bounds.default_face_pair()); + + CGAL_assertion(h_global_bounds.lower >= FT(0)); + if (local_bounds.lower > h_global_bounds.lower) { // it is (6) in the paper, see also Algorithm 1 + h_global_bounds.lower = local_bounds.lower; + h_global_bounds.lpair.first = fpair.first; + h_global_bounds.lpair.second = local_bounds.tm2_lface; + } + CGAL_assertion(h_global_bounds.upper >= FT(0)); + if (local_bounds.upper > h_global_bounds.upper) { // it is (6) in the paper, see also Algorithm 1 + h_global_bounds.upper = local_bounds.upper; + h_global_bounds.upair.first = fpair.first; + h_global_bounds.upair.second = local_bounds.tm2_uface; + } + CGAL_assertion(h_global_bounds.upper >= h_global_bounds.lower); + + // Store the triangle given as primitive here as candidate triangle + // together with the local bounds it obtained to send it to subdivision later. + m_candidiate_triangles.push(Candidate(triangle, local_bounds, fpair.first)); + } + + // Determine whether child nodes will still contribute to a larger + // Hausdorff distance and thus have to be entered. + template + std::pair + do_intersect_with_priority(const Query&, const Node& node) { + + // Check if we can stop already here. Since our bounds only grow, in case, we are + // above the user-defined max distance bound, we return. This way, the user can + // early detect that he is behind his thresholds. + if (m_distance_bound >= FT(0) && !m_early_quit) { + + CGAL_assertion(h_global_bounds.lower >= FT(0)); + CGAL_assertion(h_global_bounds.upper >= FT(0)); + CGAL_assertion(h_global_bounds.upper >= h_global_bounds.lower); + + const FT hdist = (h_global_bounds.lower + h_global_bounds.upper) / FT(2); + m_early_quit = (hdist >= m_distance_bound); + // std::cout << "- hdist: " << hdist << std::endl; + // std::cout << "- early quit: " << m_early_quit << std::endl; + } + if (m_early_quit) return std::make_pair(false, FT(0)); + + // Have reached a node, determine whether or not to enter it. + // Get the bounding box of the nodes. + const auto bbox = node.bbox(); + + // Compute its center. + const Point_3 center = Point_3( + (bbox.min(0) + bbox.max(0)) / FT(2), + (bbox.min(1) + bbox.max(1)) / FT(2), + (bbox.min(2) + bbox.max(2)) / FT(2)); + + // Find the point from TM2 closest to the center. + const Point_3 closest = m_tm2_tree.closest_point(center); + + // Compute the difference vector between the bbox center and the closest point in tm2. + Vector_3 difference = Vector_3(closest, center); + + // Shift the vector to be the difference between the farthest corner + // of the bounding box away from the closest point on TM2. + FT diff_x = (bbox.max(0) - bbox.min(0)) / FT(2); + if (difference.x() < 0) diff_x = diff_x * -FT(1); + FT diff_y = (bbox.max(1) - bbox.min(1)) / FT(2); + if (difference.y() < 0) diff_y = diff_y * -FT(1); + FT diff_z = (bbox.max(2) - bbox.min(2)) / FT(2); + if (difference.z() < 0) diff_z = diff_z * -FT(1); + difference = difference + Vector_3(diff_x, diff_y, diff_z); // it is (9) in the paper + + // Compute distance from the farthest corner of the bbox to the closest point in TM2. + const FT dist = CGAL::approximate_sqrt(difference.squared_length()); + + // See Algorithm 1 here. + // If the distance is larger than the global lower bound, enter the node, i.e. return true. + CGAL_assertion(h_global_bounds.lower >= FT(0)); + if (dist > h_global_bounds.lower) { + return std::make_pair(true , +dist); + } else { + return std::make_pair(false, FT(0)); + } + } + + template + bool do_intersect(const Query& query, const Node& node) { + return this->do_intersect_with_priority(query, node).first; + } + + template + void traverse_group(const Query& query, PrimitiveConstIterator group_begin, PrimitiveConstIterator group_end) { + CGAL_assertion_msg(false, "ERROR: we should not call the group traversal on TM1!"); + for (PrimitiveConstIterator it = group_begin; it != group_end; ++it) { + this->intersection(query, *it); + } + } + + bool early_quit() const { + return m_early_quit; + } + + // Return those triangles from TM1, which are candidates for including a + // point realizing the Hausdorff distance. + Heap_type& get_candidate_triangles() { + return m_candidiate_triangles; + } + + // Return the global Hausdorff bounds computed for the passed query triangle. + Global_bounds get_global_bounds() { + + CGAL_assertion(h_global_bounds.lower >= FT(0)); + CGAL_assertion(h_global_bounds.upper >= FT(0)); + CGAL_assertion(h_global_bounds.upper >= h_global_bounds.lower); + + update_global_bounds(); + return h_global_bounds; + } + + // Here, we return the maximum distance from one of the face corners + // to the second mesh. We also return a pair of realizing this distance faces. + FT get_maximum_distance( + const Face_handle_1 tm1_lface, + std::pair& fpair) const { + + const auto triangle = get(m_face_to_triangle_map, tm1_lface); + const Point_3 v0 = triangle.vertex(0); + const Point_3 v1 = triangle.vertex(1); + const Point_3 v2 = triangle.vertex(2); + + const auto pair0 = m_tm2_tree.closest_point_and_primitive(v0); + const auto pair1 = m_tm2_tree.closest_point_and_primitive(v1); + const auto pair2 = m_tm2_tree.closest_point_and_primitive(v2); + + const auto sq_dist0 = std::make_pair( + CGAL::squared_distance(v0, pair0.first), pair0.second); + const auto sq_dist1 = std::make_pair( + CGAL::squared_distance(v1, pair1.first), pair1.second); + const auto sq_dist2 = std::make_pair( + CGAL::squared_distance(v2, pair2.first), pair2.second); + + const auto mdist1 = (sq_dist0.first > sq_dist1.first) ? sq_dist0 : sq_dist1; + const auto mdist2 = (mdist1.first > sq_dist2.first) ? mdist1 : sq_dist2; + + Face_handle_2 tm2_uface = mdist2.second; + fpair = std::make_pair(tm1_lface, tm2_uface); + return CGAL::approximate_sqrt(mdist2.first); + } + + private: + // Input data. + const AABBTraits& m_traits; + const TriangleMesh1& m_tm1; + const TriangleMesh2& m_tm2; + const VPM1& m_vpm1; + const VPM2& m_vpm2; + const TM2_tree& m_tm2_tree; + const TM1_face_to_triangle_map m_face_to_triangle_map; + + // Internal bounds and values. + const FT m_error_bound; + const FT m_infinity_value; + const FT m_initial_bound; + const FT m_distance_bound; + Global_bounds h_global_bounds; + bool m_early_quit; + + // All candidate triangles. + Heap_type m_candidiate_triangles; + + // In case, we did not enter any loop, we set the realizing triangles here. + void update_global_bounds() { + + if (m_candidiate_triangles.size() > 0) { + const auto top = m_candidiate_triangles.top(); + + if (h_global_bounds.lpair.first == Face_handle_1()) + h_global_bounds.lpair.first = top.tm1_face; + if (h_global_bounds.lpair.second == Face_handle_2()) + h_global_bounds.lpair.second = top.bounds.tm2_lface; + + if (h_global_bounds.upair.first == Face_handle_1()) + h_global_bounds.upair.first = top.tm1_face; + if (h_global_bounds.upair.second == Face_handle_2()) + h_global_bounds.upair.second = top.bounds.tm2_uface; + + } else { + + std::pair fpair; + get_maximum_distance(*(faces(m_tm1).begin()), fpair); + CGAL_assertion(fpair.first == *(faces(m_tm1).begin())); + + if (h_global_bounds.lpair.first == Face_handle_1()) + h_global_bounds.lpair.first = fpair.first; + if (h_global_bounds.lpair.second == Face_handle_2()) + h_global_bounds.lpair.second = fpair.second; + + if (h_global_bounds.upair.first == Face_handle_1()) + h_global_bounds.upair.first = fpair.first; + if (h_global_bounds.upair.second == Face_handle_2()) + h_global_bounds.upair.second = fpair.second; + } + } + }; +} + +#endif // CGAL_PMP_INTERNAL_AABB_TRAVERSAL_TRAITS_WITH_HAUSDORFF_DISTANCE diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 4af1dc42dfb4..c3e7dfb367ad 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -54,6 +54,7 @@ else() message(STATUS "Examples that use OpenMesh will not be compiled.") endif() +create_single_source_cgal_program("test_hausdorff_bounded_error_distance.cpp") create_single_source_cgal_program("test_pmp_read_polygon_mesh.cpp") create_single_source_cgal_program("connected_component_polyhedron.cpp") create_single_source_cgal_program("connected_component_surface_mesh.cpp") @@ -109,7 +110,16 @@ create_single_source_cgal_program("triangulate_hole_with_cdt_2_test.cpp") create_single_source_cgal_program("test_pmp_polyhedral_envelope.cpp") # create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp") +find_package(METIS) +include(CGAL_METIS_support) +if(TARGET CGAL::METIS_support) + target_link_libraries(test_hausdorff_bounded_error_distance PUBLIC CGAL::METIS_support) +else() + message(STATUS "Tests, which use the METIS library will not be compiled.") +endif() + if(TARGET CGAL::TBB_support) + target_link_libraries(test_hausdorff_bounded_error_distance PUBLIC CGAL::TBB_support) target_link_libraries(test_pmp_distance PUBLIC CGAL::TBB_support) target_link_libraries(orient_polygon_soup_test PUBLIC CGAL::TBB_support) target_link_libraries(self_intersection_surface_mesh_test diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tetrahedron-remeshed.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tetrahedron-remeshed.off new file mode 100644 index 000000000000..82fe31012a09 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tetrahedron-remeshed.off @@ -0,0 +1,340 @@ +OFF +114 224 0 +0 0 0 +1 0 0 +0 1 0 +0 0 1 +0.5 0.5 0 +0.5 0 0 +0 0 0.5 +0 0.5 0.5 +0.5 0 0.5 +0.5 0.25 0 +0.25 0 0 +0.25 0.25 0 +0 0.25 0.5 +0 0.25 0.25 +0 0 0.25 +0.25 0.25 0.5 +0.5 0.25 0.25 +0.25 0.5 0.25 +0.25 0 0.5 +0.25 0 0.25 +0.5 0 0.25 +0.25 0.75 0 +0 0.5 0 +0.75 0.25 0 +0.75 0 0 +0 0 0.75 +0 0.25 0.75 +0 0.75 0.25 +0.25 0 0.75 +0.75 0 0.25 +0.375 0.125 0 +0.25 0.125 0 +0.375 0.25 0 +0 0.25 0.375 +0 0.125 0.25 +0 0.125 0.375 +0.375 0.25 0.375 +0.375 0.375 0.25 +0.25 0.375 0.375 +0.25 0 0.375 +0.375 0 0.25 +0.375 0 0.375 +0.25 0.5 0 +0.125 0.375 0 +0.125 0.625 0 +0.75 0.125 0 +0.625 0.125 0 +0.625 0.25 0 +0 0.125 0.75 +0 0.25 0.625 +0 0.125 0.625 +0 0.5 0.25 +0 0.625 0.125 +0 0.375 0.125 +0.125 0.625 0.25 +0.25 0.625 0.125 +0.125 0.7500001 0.125 +0.125 0.125 0.75 +0.25 0.125 0.625 +0.125 0.25 0.625 +0.625 0.125 0.25 +0.7500001 0.125 0.125 +0.625 0.25 0.125 +0.625 0 0.25 +0.625 0 0.125 +0.75 0 0.125 +0.125 0 0.75 +0.125 0 0.625 +0.25 0 0.625 +0.125 0 0.25 +0.125 0 0.125 +0.25 0 0.125 +0.5 0.375 0 +0.375 0.375 0 +0.5 0.125 0 +0.375 0 0 +0.125 0 0 +0.125 0.125 0 +0 0.125 0.5 +0 0 0.375 +0 0.375 0.5 +0 0.375 0.375 +0 0.125 0.125 +0 0 0.125 +0.125 0.375 0.5 +0.125 0.5 0.375 +0.375 0.125 0.5 +0.5 0.125 0.375 +0.5 0.375 0.125 +0.375 0.5 0.125 +0.375 0 0.5 +0.5 0 0.375 +0.125 0 0.5 +0.125 0 0.375 +0.375 0 0.125 +0.5 0 0.125 +0.125 0.875 0 +0 0.75 0 +0.375 0.625 0 +0 0.25 0 +0.625 0.375 0 +0.875 0.125 0 +0.875 0 0 +0.625 0 0 +0 0 0.625 +0 0 0.875 +0 0.125 0.875 +0 0.375 0.625 +0 0.625 0.375 +0 0.875 0.125 +0.125 0 0.875 +0.375 0 0.625 +0.625 0 0.375 +0.875 0 0.125 +3 30 31 32 +3 33 34 35 +3 36 37 38 +3 39 40 41 +3 42 43 44 +3 45 46 47 +3 48 49 50 +3 51 52 53 +3 54 55 56 +3 57 58 59 +3 60 61 62 +3 63 64 65 +3 66 67 68 +3 69 70 71 +3 72 32 73 +3 74 75 30 +3 31 76 77 +3 78 35 79 +3 80 81 33 +3 34 82 83 +3 84 38 85 +3 86 87 36 +3 37 88 89 +3 90 41 91 +3 92 93 39 +3 40 94 95 +3 96 44 97 +3 98 73 42 +3 43 77 99 +3 100 47 72 +3 101 102 45 +3 46 103 74 +3 104 50 78 +3 105 106 48 +3 49 107 80 +3 82 53 99 +3 81 108 51 +3 52 109 97 +3 109 56 96 +3 108 85 54 +3 55 89 98 +3 107 59 84 +3 106 110 57 +3 58 111 86 +3 88 62 100 +3 87 112 60 +3 61 113 101 +3 113 65 102 +3 112 91 63 +3 64 95 103 +3 111 68 90 +3 110 105 66 +3 67 104 92 +3 94 71 75 +3 93 79 69 +3 70 83 76 +3 9 30 32 +3 30 10 31 +3 32 31 11 +3 12 33 35 +3 33 13 34 +3 35 34 14 +3 15 36 38 +3 36 16 37 +3 38 37 17 +3 18 39 41 +3 39 19 40 +3 41 40 20 +3 21 42 44 +3 42 11 43 +3 44 43 22 +3 23 45 47 +3 45 24 46 +3 47 46 9 +3 25 48 50 +3 48 26 49 +3 50 49 12 +3 13 51 53 +3 51 27 52 +3 53 52 22 +3 27 54 56 +3 54 17 55 +3 56 55 21 +3 26 57 59 +3 57 28 58 +3 59 58 15 +3 16 60 62 +3 60 29 61 +3 62 61 23 +3 29 63 65 +3 63 20 64 +3 65 64 24 +3 28 66 68 +3 66 25 67 +3 68 67 18 +3 19 69 71 +3 69 14 70 +3 71 70 10 +3 4 72 73 +3 72 9 32 +3 73 32 11 +3 9 74 30 +3 74 5 75 +3 30 75 10 +3 11 31 77 +3 31 10 76 +3 77 76 0 +3 6 78 79 +3 78 12 35 +3 79 35 14 +3 12 80 33 +3 80 7 81 +3 33 81 13 +3 14 34 83 +3 34 13 82 +3 83 82 0 +3 7 84 85 +3 84 15 38 +3 85 38 17 +3 15 86 36 +3 86 8 87 +3 36 87 16 +3 17 37 89 +3 37 16 88 +3 89 88 4 +3 8 90 91 +3 90 18 41 +3 91 41 20 +3 18 92 39 +3 92 6 93 +3 39 93 19 +3 20 40 95 +3 40 19 94 +3 95 94 5 +3 2 96 97 +3 96 21 44 +3 97 44 22 +3 21 98 42 +3 98 4 73 +3 42 73 11 +3 22 43 99 +3 43 11 77 +3 99 77 0 +3 4 100 72 +3 100 23 47 +3 72 47 9 +3 23 101 45 +3 101 1 102 +3 45 102 24 +3 9 46 74 +3 46 24 103 +3 74 103 5 +3 6 104 78 +3 104 25 50 +3 78 50 12 +3 25 105 48 +3 105 3 106 +3 48 106 26 +3 12 49 80 +3 49 26 107 +3 80 107 7 +3 0 82 99 +3 82 13 53 +3 99 53 22 +3 13 81 51 +3 81 7 108 +3 51 108 27 +3 22 52 97 +3 52 27 109 +3 97 109 2 +3 2 109 96 +3 109 27 56 +3 96 56 21 +3 27 108 54 +3 108 7 85 +3 54 85 17 +3 21 55 98 +3 55 17 89 +3 98 89 4 +3 7 107 84 +3 107 26 59 +3 84 59 15 +3 26 106 57 +3 106 3 110 +3 57 110 28 +3 15 58 86 +3 58 28 111 +3 86 111 8 +3 4 88 100 +3 88 16 62 +3 100 62 23 +3 16 87 60 +3 87 8 112 +3 60 112 29 +3 23 61 101 +3 61 29 113 +3 101 113 1 +3 1 113 102 +3 113 29 65 +3 102 65 24 +3 29 112 63 +3 112 8 91 +3 63 91 20 +3 24 64 103 +3 64 20 95 +3 103 95 5 +3 8 111 90 +3 111 28 68 +3 90 68 18 +3 28 110 66 +3 110 3 105 +3 66 105 25 +3 18 67 92 +3 67 25 104 +3 92 104 6 +3 5 94 75 +3 94 19 71 +3 75 71 10 +3 19 93 69 +3 93 6 79 +3 69 79 14 +3 10 70 76 +3 70 14 83 +3 76 83 0 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tetrahedron.off b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tetrahedron.off new file mode 100644 index 000000000000..70848d1c6058 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/data/tetrahedron.off @@ -0,0 +1,10 @@ +OFF +4 4 0 +0 0 0 +1 0 0 +0 1 0 +0 0 1 +3 2 1 0 +3 0 3 2 +3 2 3 1 +3 1 3 0 diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_hausdorff_bounded_error_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_hausdorff_bounded_error_distance.cpp new file mode 100644 index 000000000000..56c1557ca72b --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_hausdorff_bounded_error_distance.cpp @@ -0,0 +1,1257 @@ +// Use it to include parallel computations in the bounded error Hausdorff distance. +// #define USE_PARALLEL_BEHD + +// Use this def in order to get all DEBUG info related to the bounded-error Hausdorff code! +// #define CGAL_HAUSDORFF_DEBUG + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +using SCK = CGAL::Simple_cartesian; +using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel; +using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; + +using Kernel = EPICK; +using FT = typename Kernel::FT; +using Point_3 = typename Kernel::Point_3; +using Vector_3 = typename Kernel::Vector_3; +using Triangle_3 = typename Kernel::Triangle_3; + +using TAG = CGAL::Sequential_tag; +using Surface_mesh = CGAL::Surface_mesh; +using Polyhedron = CGAL::Polyhedron_3; +using Affine_transformation_3 = CGAL::Aff_transformation_3; +using Timer = CGAL::Real_timer; + +using Face_handle = typename boost::graph_traits::face_descriptor; +namespace PMP = CGAL::Polygon_mesh_processing; + +struct Approximate_hd_wrapper { + const double m_num_samples; + std::string name() const { return "approximate"; } + Approximate_hd_wrapper(const double num_samples) : m_num_samples(num_samples) { } + double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + return PMP::approximate_Hausdorff_distance(tm1, tm2, + PMP::parameters::number_of_points_per_area_unit(m_num_samples), + PMP::parameters::number_of_points_per_area_unit(m_num_samples)); + } + double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + return PMP::approximate_symmetric_Hausdorff_distance(tm1, tm2, + PMP::parameters::number_of_points_per_area_unit(m_num_samples), + PMP::parameters::number_of_points_per_area_unit(m_num_samples)); + } +}; + +struct Naive_bounded_error_hd_wrapper { + const double m_error_bound; + std::string name() const { return "naive bounded error"; } + Naive_bounded_error_hd_wrapper(const double error_bound) : m_error_bound(error_bound) { } + double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + return PMP::bounded_error_Hausdorff_distance_naive(tm1, tm2, m_error_bound); + } + double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + const double dista = operator()(tm1, tm2); + const double distb = operator()(tm2, tm1); + return (CGAL::max)(dista, distb); + } +}; + +struct Bounded_error_hd_wrapper { + const double m_error_bound; + std::string name() const { return "bounded error"; } + Bounded_error_hd_wrapper(const double error_bound) : m_error_bound(error_bound) { } + double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + return PMP::bounded_error_Hausdorff_distance(tm1, tm2, m_error_bound); + } + double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + return PMP::bounded_error_symmetric_Hausdorff_distance(tm1, tm2, m_error_bound); + } +}; + +template +void get_mesh(const std::string filepath, PolygonMesh& mesh) { + + mesh.clear(); + std::ifstream input(filepath); + input >> mesh; + std::cout << "* getting mesh with " << faces(mesh).size() << " faces" << std::endl; +} + +template +void get_meshes( + const std::string filepath1, const std::string filepath2, + PolygonMesh1& mesh1, PolygonMesh2& mesh2) { + + get_mesh(filepath1, mesh1); + get_mesh(filepath2, mesh2); +} + +template +void save_mesh(const PolygonMesh& mesh, const std::string filepath) { + + if (!CGAL::IO::write_PLY(filepath + ".ply", mesh)) { + std::cerr << "ERROR: cannot save this file: " << filepath << std::endl; + exit(EXIT_FAILURE); + } +} + +// An easy example of a tetrahedron and its remeshed version. +void remeshing_tetrahedon_example( + const double error_bound, const bool save = false) { + + Timer timer; + Surface_mesh mesh1, mesh2; + std::cout << std::endl << "* (E1) remeshing Tetrahedron example:" << std::endl; + + CGAL::make_tetrahedron( + Point_3(0, 0, 0), Point_3(2, 0, 0), + Point_3(1, 1, 1), Point_3(1, 0, 2), mesh1); + mesh2 = mesh1; + + using edge_descriptor = typename boost::graph_traits::edge_descriptor; + Surface_mesh::Property_map is_constrained_map = + mesh2.add_property_map("e:is_constrained", true).first; + + const double target_edge_length = 0.05; + PMP::isotropic_remeshing( + mesh2.faces(), target_edge_length, mesh2, + PMP::parameters::edge_is_constrained_map(is_constrained_map)); + + if (save) save_mesh(mesh1, "mesh1"); + if (save) save_mesh(mesh2, "mesh2"); + + timer.reset(); + timer.start(); + const double hdist = PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound); + std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; + timer.stop(); + std::cout << "* processing time: " << timer.time() << " s." << std::endl; + assert(hdist == error_bound); +} + +// Example with a point realizing the Hausdorff distance strictly +// lying in the interior of a triangle. +void interior_triangle_example( + const double error_bound, const bool save = false) { + + Timer timer; + Surface_mesh mesh1, mesh2; + std::cout << std::endl << "* (E2) interior Triangle example:" << std::endl; + + mesh1.add_vertex(Point_3(-1, 1, 1)); + mesh1.add_vertex(Point_3( 0, -1, 1)); + mesh1.add_vertex(Point_3( 1, 1, 1)); + mesh1.add_face(mesh1.vertices()); + + auto v0 = mesh2.add_vertex(Point_3(-1.0, 1, 0)); + auto v1 = mesh2.add_vertex(Point_3( 0.0, -1, 0)); + auto v2 = mesh2.add_vertex(Point_3( 1.0, 1, 0)); + auto v3 = mesh2.add_vertex(Point_3( 0.0, 1, -1)); + auto v4 = mesh2.add_vertex(Point_3(-0.5, 0, -1)); + auto v5 = mesh2.add_vertex(Point_3( 0.5, 0, -1)); + mesh2.add_face(v0, v3, v4); + mesh2.add_face(v1, v4, v5); + mesh2.add_face(v2, v5, v3); + + if (save) save_mesh(mesh1, "mesh1"); + if (save) save_mesh(mesh2, "mesh2"); + + timer.reset(); + timer.start(); + const double hdist = PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound); + std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; + timer.stop(); + std::cout << "* processing time: " << timer.time() << " s." << std::endl; + assert(hdist >= 1.0); +} + +// Read a real mesh given by the user, perturb it slightly, and compute the +// Hausdorff distance between the original mesh and its pertubation. +void perturbing_surface_mesh_example( + const std::string filepath, const double error_bound, const bool save = false) { + + Timer timer; + std::cout << std::endl << "* (E3) perturbing Surface Mesh example:" << std::endl; + + Surface_mesh mesh1, mesh2; + get_meshes(filepath, filepath, mesh1, mesh2); + + const double max_size = 0.1; + PMP::random_perturbation( + mesh2.vertices(), mesh2, max_size, CGAL::parameters::do_project(false)); + std::cout << "* perturbing the second mesh" << std::endl; + + if (save) save_mesh(mesh2, "mesh2"); + + timer.reset(); + timer.start(); + const double hdist = PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound); + std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; + timer.stop(); + std::cout << "* processing time: " << timer.time() << " s." << std::endl; + assert(hdist > 0.0); +} + +// Read two meshes and store them in two different face graph containers, +// perturb the second mesh, and compute the Hausdorff distance. +void perturbing_polyhedron_mesh_example( + const std::string filepath, const double error_bound) { + + Timer timer; + std::cout << std::endl << "* (E3) perturbing Polyhedron mesh example:" << std::endl; + + Surface_mesh mesh1; + Polyhedron mesh2; + get_mesh(filepath, mesh1); + get_mesh(filepath, mesh2); + + const double max_size = 0.1; + PMP::random_perturbation( + vertices(mesh2), mesh2, max_size, CGAL::parameters::do_project(false)); + std::cout << "* perturbing the second mesh" << std::endl; + + timer.reset(); + timer.start(); + const double hdist1 = PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound); + const double hdist2 = PMP::bounded_error_Hausdorff_distance(mesh2, mesh1, error_bound); + std::cout << "* bounded-error Hausdorff distance 1->2: " << hdist1 << std::endl; + std::cout << "* bounded-error Hausdorff distance 2->1: " << hdist2 << std::endl; + timer.stop(); + std::cout << "* processing time: " << timer.time() << " s." << std::endl; + assert(hdist1 > 0.0); + assert(hdist2 > 0.0); + assert(hdist2 > hdist1); + + const double hdist3 = PMP::bounded_error_symmetric_Hausdorff_distance(mesh1, mesh2, error_bound); + assert(hdist3 == (CGAL::max)(hdist1, hdist2)); +} + +// Read two meshes given by the user, initially place them at their originally +// given position. Move the second mesh in 300 steps away from the first one. +// Print how the Hausdorff distance changes. +void moving_surface_mesh_example( + const std::string filepath1, const std::string filepath2, + const std::size_t n, const double error_bound, const bool save = false) { + + Timer timer; + std::cout << std::endl << "* (E4) moving Surface Mesh example:" << std::endl; + + Surface_mesh mesh1, mesh2; + get_meshes(filepath1, filepath2, mesh1, mesh2); + + const auto bbox = PMP::bbox(mesh2); + const FT distance = CGAL::approximate_sqrt(CGAL::squared_distance( + Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin()), + Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()))); + + const FT t = FT(1) / FT(100); + if (save) save_mesh(mesh2, "mesh-0"); + + double curr_dist = 0.0; + for (std::size_t i = 0; i < n; ++i) { + PMP::transform(Affine_transformation_3(CGAL::Translation(), + Vector_3(t * distance, t * distance, t * distance)), mesh2); + + timer.reset(); + timer.start(); + const double hdist = PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound); + std::cout << "* position: " << i << std::endl; + std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; + timer.stop(); + std::cout << "* processing time: " << timer.time() << " s." << std::endl; + if (save) save_mesh(mesh2, "mesh-" + std::to_string(i + 1)); + assert(hdist > curr_dist); + curr_dist = hdist; + } +} + +template +void test_0(const FunctionWrapper& functor, const bool save = false) { + + // The same triangle. + // Expected distance is 0. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 0 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + mesh2 = mesh1; + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 0): " << dista << std::endl; + std::cout << "* HInverted distance (expected 0): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 0): " << distc << std::endl; + + assert(dista == 0.0); + assert(distb == 0.0); + assert(distc == naive); +} + +template +void test_1(const FunctionWrapper& functor, const bool save = false) { + + // Two triangles are parallel and 1 unit distance away from each other. + // Expected distance is 1. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 1 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + mesh2.add_vertex(Point_3(0, 0, 1)); + mesh2.add_vertex(Point_3(2, 0, 1)); + mesh2.add_vertex(Point_3(1, 1, 1)); + mesh2.add_face(mesh2.vertices()); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 1): " << dista << std::endl; + std::cout << "* HInverted distance (expected 1): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 1): " << distc << std::endl; + + assert(dista == 1.0); + assert(distb == 1.0); + assert(distc == naive); +} + +template +void test_2(const FunctionWrapper& functor, const bool save = false) { + + // One triangle is orthogonal to the other one and shares a common edge. + // Expected distance is 1. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 2 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + mesh2.add_vertex(Point_3(0, 0, 0)); + mesh2.add_vertex(Point_3(2, 0, 0)); + mesh2.add_vertex(Point_3(1, 0, 1)); + mesh2.add_face(mesh2.vertices()); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 1): " << dista << std::endl; + std::cout << "* HInverted distance (expected 1): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 1): " << distc << std::endl; + + assert(dista == 1.0); + assert(distb == 1.0); + assert(distc == naive); +} + +template +void test_3(const FunctionWrapper& functor, const bool save = false) { + + // One triangle is orthogonal to the other one and shares a common edge + // that is moved 1 unit distance away. + // Expected distances are sqrt(2) and 2. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 3 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + mesh2.add_vertex(Point_3(0, 0, 1)); + mesh2.add_vertex(Point_3(2, 0, 1)); + mesh2.add_vertex(Point_3(1, 0, 2)); + mesh2.add_face(mesh2.vertices()); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected sqrt(2)): " << dista << std::endl; + std::cout << "* HInverted distance (expected 2 ): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 2 ): " << distc << std::endl; + + assert(CGAL::abs(dista - CGAL::sqrt(2.0)) < 1e-5); // error bound is 1e-4 + assert(distb == 2.0); + assert(distc == naive); +} + +template +void test_4(const FunctionWrapper& functor, const bool save = false) { + + // One triangle is orthogonal to the other one and shares a common vertex. + // Expected distance is 1.2247448713915889407. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 4 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + mesh2.add_vertex(Point_3(0, 1, 1)); + mesh2.add_vertex(Point_3(2, 1, 1)); + mesh2.add_vertex(Point_3(1, 1, 0)); + mesh2.add_face(mesh2.vertices()); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 1.22): " << dista << std::endl; + std::cout << "* HInverted distance (expected 1.22): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 1.22): " << distc << std::endl; + + assert(CGAL::abs(dista - 1.224744) < 1e-5); // error bound is 1e-4 + assert(CGAL::abs(distb - 1.224744) < 1e-5); // error bound is 1e-4 + assert(dista == distb); + assert(distc == naive); +} + +template +void test_5(const FunctionWrapper& functor, const bool save = false) { + + // One triangle is orthogonal to the other one and shares a common vertex + // that is moved 1 unit distance away. + // Expected distances are 1.7320508075688771932 and 2.1213203435596423851. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 5 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + mesh2.add_vertex(Point_3(0, 1, 2)); + mesh2.add_vertex(Point_3(2, 1, 2)); + mesh2.add_vertex(Point_3(1, 1, 1)); + mesh2.add_face(mesh2.vertices()); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 1.73): " << dista << std::endl; + std::cout << "* HInverted distance (expected 2.12): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 2.12): " << distc << std::endl; + + assert(CGAL::abs(dista - 1.732050) < 1e-5); // error bound is 1e-4 + assert(CGAL::abs(distb - 2.121320) < 1e-5); // error bound is 1e-4 + assert(distc == naive); +} + +template +void test_6(const FunctionWrapper& functor, const bool save = false) { + + // The first and second mesh have different number of triangles. + // They are parallel and lie at the same plane. The middle triangle is overlapping. + // Expected distances are 0 and 0.70710678118654757274. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 6 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + const auto v0 = mesh2.add_vertex(Point_3(0, 0, 0)); + const auto v1 = mesh2.add_vertex(Point_3(2, 0, 0)); + const auto v2 = mesh2.add_vertex(Point_3(1, 1, 0)); + const auto v3 = mesh2.add_vertex(Point_3(2, 1, 0)); + const auto v4 = mesh2.add_vertex(Point_3(0, 1, 0)); + + mesh2.add_face(v0, v1, v2); + mesh2.add_face(v2, v1, v3); + mesh2.add_face(v0, v2, v4); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 0.0): " << dista << std::endl; + std::cout << "* HInverted distance (expected 0.7): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 0.7): " << distc << std::endl; + + assert(dista == 0.0); + assert(CGAL::abs(distb - 0.707106) < 1e-5); // error bound is 1e-4 + assert(distc == naive); +} + +template +void test_7(const FunctionWrapper& functor, const bool save = false) { + + // One triangle is moved to 0.5 unit distance away from 3 other triangles. + // The first and second meshes are parallel. + // Expected distances are 0.5 and 0.86602540378443859659. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 7 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + const auto v0 = mesh2.add_vertex(Point_3(0, 0, 0.5)); + const auto v1 = mesh2.add_vertex(Point_3(2, 0, 0.5)); + const auto v2 = mesh2.add_vertex(Point_3(1, 1, 0.5)); + const auto v3 = mesh2.add_vertex(Point_3(2, 1, 0.5)); + const auto v4 = mesh2.add_vertex(Point_3(0, 1, 0.5)); + + mesh2.add_face(v0, v1, v2); + mesh2.add_face(v2, v1, v3); + mesh2.add_face(v0, v2, v4); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 0.50): " << dista << std::endl; + std::cout << "* HInverted distance (expected 0.86): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 0.86): " << distc << std::endl; + + assert(dista == 0.5); + assert(CGAL::abs(distb - 0.866025) < 1e-5); // error bound is 1e-4 + assert(distc == naive); +} + +template +void test_8(const FunctionWrapper& functor, const bool save = false) { + + // One mesh has one triangle at zero level, another mesh has two triangles + // where the first one is at level 1 and the second one is at level 2. + // Expected distances are 1 and 2. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 8 ---- " << std::endl; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + if (save) save_mesh(mesh1, "mesh1"); + + const auto v0 = mesh2.add_vertex(Point_3(0, 0, 1)); + const auto v1 = mesh2.add_vertex(Point_3(2, 0, 1)); + const auto v2 = mesh2.add_vertex(Point_3(1, 1, 1)); + const auto v3 = mesh2.add_vertex(Point_3(0, 0, 2)); + const auto v4 = mesh2.add_vertex(Point_3(2, 0, 2)); + const auto v5 = mesh2.add_vertex(Point_3(1, 1, 2)); + + mesh2.add_face(v0, v1, v2); + mesh2.add_face(v3, v4, v5); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 1): " << dista << std::endl; + std::cout << "* HInverted distance (expected 2): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 2): " << distc << std::endl; + + assert(dista == 1.0); + assert(distb == 2.0); + assert(distc == naive); +} + +template +void test_9(const FunctionWrapper& functor, const bool save = false) { + + // Two meshes partially overlap, have 2 triangles in common and each one has + // two its own trianles. All triangles form a Z shape where the height is 1. + // The expected result is 1. + + std::cout.precision(20); + Surface_mesh mesh1, mesh2; + std::cout << " ---- testing 9 ---- " << std::endl; + + auto v0 = mesh1.add_vertex(Point_3(0, 0, 0)); + auto v1 = mesh1.add_vertex(Point_3(1, 0, 0)); + auto v2 = mesh1.add_vertex(Point_3(0, 1, 0)); + auto v3 = mesh1.add_vertex(Point_3(1, 1, 0)); + auto v4 = mesh1.add_vertex(Point_3(1, 0, 1)); + auto v5 = mesh1.add_vertex(Point_3(1, 1, 1)); + mesh1.add_face(v0, v1, v2); + mesh1.add_face(v2, v1, v3); + mesh1.add_face(v1, v4, v3); + mesh1.add_face(v3, v4, v5); + if (save) save_mesh(mesh1, "mesh1"); + + v0 = mesh2.add_vertex(Point_3(2, 0, 1)); + v1 = mesh2.add_vertex(Point_3(1, 0, 0)); + v2 = mesh2.add_vertex(Point_3(2, 1, 1)); + v3 = mesh2.add_vertex(Point_3(1, 1, 0)); + v4 = mesh2.add_vertex(Point_3(1, 0, 1)); + v5 = mesh2.add_vertex(Point_3(1, 1, 1)); + mesh2.add_face(v1, v4, v3); + mesh2.add_face(v3, v4, v5); + mesh2.add_face(v4, v0, v5); + mesh2.add_face(v5, v0, v2); + if (save) save_mesh(mesh2, "mesh2"); + + const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double naive = (CGAL::max)(dista, distb); + const double distc = functor.symmetric(mesh1, mesh2); + + std::cout << "* Hausdorff distance (expected 1): " << dista << std::endl; + std::cout << "* HInverted distance (expected 1): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 1): " << distc << std::endl; + + assert(dista == 1.0); + assert(distb == 1.0); + assert(distc == naive); +} + +template +void test_synthetic_data(const FunctionWrapper& functor) { + + std::cout << std::endl << "-- test synthetic data:" << std::endl << std::endl; + std::cout << "* name -> " << functor.name() << std::endl; + + test_0(functor); // 1 parallel + test_1(functor); + test_2(functor); // 1 edge touching + test_3(functor); + test_4(functor); // 1 vertex touching + test_5(functor); + test_6(functor); // 1 to multiple + test_7(functor); + test_8(functor); + test_9(functor); // overlapping +} + +template< +typename FunctionWrapper1, +typename FunctionWrapper2> +void test_one_versus_another( + const FunctionWrapper1& functor1, + const FunctionWrapper2& functor2) { + + std::cout.precision(20); + const std::string filepath1 = "data/tetrahedron.off"; + const std::string filepath2 = "data/tetrahedron-remeshed.off"; + + std::cout << std::endl << "-- test one versus another (tetrahedron):" << std::endl << std::endl; + std::cout << "* name 1 -> " << functor1.name() << std::endl; + std::cout << "* name 2 -> " << functor2.name() << std::endl; + + Surface_mesh mesh1, mesh2; + get_meshes(filepath1, filepath2, mesh1, mesh2); + + // TEST 0. + // Load and compare. + // The expected distance is 0. + std::cout << " ---- testing 0 ---- " << std::endl; + const double dista0 = functor1(mesh1, mesh2); + const double distb0 = functor2(mesh1, mesh2); + std::cout << "* Hausdorff distance1: " << dista0 << std::endl; + std::cout << "* Hausdorff distance2: " << distb0 << std::endl; + + std::cout << "* traslating by 1 unit ..." << std::endl; + PMP::transform(Affine_transformation_3(CGAL::Translation(), + Vector_3(FT(0), FT(0), FT(1))), mesh2); + + // TEST 1. + // Translate by 1 unit distance and compare. + // The expected distance is 1. + std::cout << " ---- testing 1 ---- " << std::endl; + const double dista1 = functor1(mesh1, mesh2); + const double distb1 = functor2(mesh1, mesh2); + std::cout << "* Hausdorff distance1: " << dista1 << std::endl; + std::cout << "* Hausdorff distance2: " << distb1 << std::endl; + + assert(CGAL::abs(dista0 - distb0) < 1e-3); // error bound is 1e-4 + assert(CGAL::abs(dista1 - distb1) < 1e-3); // error bound is 1e-4 +} + +template< +typename FunctionWrapper1, +typename FunctionWrapper2> +void test_real_meshes( + const std::string filepath1, + const std::string filepath2, + const FunctionWrapper1& functor1, + const FunctionWrapper2& functor2) { + + std::cout.precision(20); + std::cout << std::endl << "-- test real meshes:" << std::endl << std::endl; + std::cout << "* input path 1: " << filepath1 << std::endl; + std::cout << "* input path 2: " << filepath2 << std::endl; + + Surface_mesh mesh1, mesh2; + get_meshes(filepath1, filepath2, mesh1, mesh2); + + std::cout << std::endl; + std::cout << "* name 1 -> " << functor1.name() << std::endl; + std::cout << "* name 2 -> " << functor2.name() << std::endl; + + // Load and compare. + std::cout << std::endl; + std::cout << " ---- testing ---- " << std::endl; + const double dista0 = functor1(mesh1, mesh2); + const double dista1 = functor1(mesh2, mesh1); + const double distb0 = functor2(mesh1, mesh2); + const double distb1 = functor2(mesh2, mesh1); + std::cout << std::endl; + std::cout << "* Hausdorff distance1 f: " << dista0 << std::endl; + std::cout << "* Hausdorff distance1 b: " << dista1 << std::endl; + std::cout << std::endl; + std::cout << "* Hausdorff distance2 f: " << distb0 << std::endl; + std::cout << "* Hausdorff distance2 b: " << distb1 << std::endl; + + assert(CGAL::abs(dista0 - distb0) < 1e-3); // error bound is 1e-4 + assert(CGAL::abs(dista1 - distb1) < 1e-3); // error bound is 1e-4 +} + +template +void test_timings(const std::string filepath, const FunctionWrapper& functor) { + + std::cout.precision(20); + std::cout << std::endl << "-- test timings: " << functor.name() << std::endl << std::endl; + + Timer timer; + Surface_mesh mesh1, mesh2; + get_mesh(filepath, mesh1); + PMP::isotropic_remeshing(faces(mesh1), 0.005, mesh1); + mesh1.collect_garbage(); + mesh2=mesh1; + + + timer.reset(); + timer.start(); + const double dista1 = functor(mesh1, mesh2); + timer.stop(); + double timea = timer.time(); + + timer.reset(); + timer.start(); + const double distb1 = functor(mesh2, mesh1); + timer.stop(); + double timeb = timer.time(); + + timer.reset(); + timer.start(); + const double distc1 = functor.symmetric(mesh1, mesh2); + timer.stop(); + double timeab = timer.time(); + + std::cout << "* time a1 (sec.): " << timea << std::endl; + std::cout << "* time b1 (sec.): " << timeb << std::endl; + std::cout << "* time ab1 naive (sec.): " << timea + timeb << std::endl; + std::cout << "* time ab1 optimized (sec.): " << timeab << std::endl; + + assert(timea > 0.0); + assert(timeb > 0.0); + assert(timeab < timea + timeb); + + PMP::transform(Affine_transformation_3(CGAL::Translation(), + Vector_3(FT(0), FT(0), FT(1))), mesh2); + + timer.reset(); + timer.start(); + const double dista2 = functor(mesh1, mesh2); + timer.stop(); + timea = timer.time(); + + timer.reset(); + timer.start(); + const double distb2 = functor(mesh2, mesh1); + timer.stop(); + timeb = timer.time(); + + timer.reset(); + timer.start(); + const double distc2 = functor.symmetric(mesh1, mesh2); + timer.stop(); + timeab = timer.time(); + + std::cout << "* time a2 (sec.): " << timea << std::endl; + std::cout << "* time b2 (sec.): " << timeb << std::endl; + std::cout << "* time ab2 naive (sec.): " << timea + timeb << std::endl; + std::cout << "* time ab2 optimized (sec.): " << timeab << std::endl; + + assert(timea > 0.0); + assert(timeb > 0.0); + assert(timeab < timea + timeb); + + std::cout << "* dista = " << dista1 << " , " << dista2 << std::endl; + std::cout << "* distb = " << distb1 << " , " << distb2 << std::endl; + std::cout << "* distab = " << distc1 << " , " << distc2 << std::endl; + + assert(dista1 == distb1 && distb1 == distc1); + assert(dista2 == distb2 && distb2 == distc2); +} + +template +void test_bunny( + const FunctionWrapper& functor, const int n = 5, const bool save = false) { + + std::cout.precision(20); + const std::string filepath1 = "data/bunny_16300.off"; // approx 16.3K + const std::string filepath2 = "data/bunny_69400.off"; // approx 69.4K + + std::cout << std::endl << "-- test bunny:" << std::endl << std::endl; + std::cout << "* name -> " << functor.name() << std::endl; + + Surface_mesh mesh1, mesh2; + get_meshes(filepath1, filepath2, mesh1, mesh2); + if (save) save_mesh(mesh1, "mesh1"); + if (save) save_mesh(mesh2, "mesh2"); + + // Get 3 times longest dimension. + const auto bbox = PMP::bbox(mesh2); + const FT dist1 = static_cast(CGAL::abs(bbox.xmax() - bbox.xmin())); + const FT dist2 = static_cast(CGAL::abs(bbox.ymax() - bbox.ymin())); + const FT dist3 = static_cast(CGAL::abs(bbox.zmax() - bbox.zmin())); + const FT dim = FT(3) * (CGAL::max)((CGAL::max)(dist1, dist2), dist3); + + // Get timings. + Timer timer; + std::vector times; + times.reserve(n); + + if (n == 0) { + + const FT distance = dim; + PMP::transform(Affine_transformation_3(CGAL::Translation(), + Vector_3(distance, distance, distance)), mesh2); + if (save) save_mesh(mesh2, "mesh2"); + timer.reset(); + timer.start(); + // const double dista = functor(mesh1, mesh2); + const double distb = functor(mesh2, mesh1); + const double distc = distb; // (CGAL::max)(dista, distb); + timer.stop(); + times.push_back(timer.time()); + std::cout << "* distance / Hausdorff / time (sec.) : " << + distance << " / " << distc << " / " << times.back() << std::endl; + assert(distc > 0.0); + + } else { + + // t is the step where n is the number of steps. + double curr_dist = 1.0; + const FT t = FT(1) / static_cast(n); + for (int k = n; k >= 0; --k) { + auto mesh = mesh2; + const FT distance = k * t * dim; + PMP::transform(Affine_transformation_3(CGAL::Translation(), + Vector_3(distance, distance, distance)), mesh); + if (save) save_mesh(mesh, "mesh2-" + std::to_string(k)); + + timer.reset(); + timer.start(); + const double dista = functor(mesh1, mesh); + const double distb = functor(mesh, mesh1); + const double distc = (CGAL::max)(dista, distb); + timer.stop(); + times.push_back(timer.time()); + std::cout << "* distance / Hausdorff / time (sec.) " << k << " : " << + distance << " / " << distc << " / " << times.back() << std::endl; + assert(distc < curr_dist); + curr_dist = distc; + } + } + + double min_time = +std::numeric_limits::infinity(); + double max_time = -std::numeric_limits::infinity(); + double avg_time = 0.0; + for (const double time : times) { + min_time = (CGAL::min)(min_time, time); + max_time = (CGAL::max)(max_time, time); + avg_time += time; + } + avg_time /= static_cast(times.size()); + + std::cout << std::endl << "* timings (msec.): " << std::endl; + std::cout << "* avg time: " << avg_time * 1000.0 << std::endl; + std::cout << "* min time: " << min_time * 1000.0 << std::endl; + std::cout << "* max time: " << max_time * 1000.0 << std::endl; + + assert(min_time <= max_time); + assert(min_time <= avg_time); + assert(avg_time <= max_time); +} + +Triangle_3 get_triangle(const int index, const Surface_mesh& mesh) { + + const auto mfaces = faces(mesh); + assert(index >= 0); + assert(index < static_cast(mfaces.size())); + const auto face = *(mfaces.begin() + index); + + const auto he = halfedge(face, mesh); + const auto vertices = vertices_around_face(he, mesh); + assert(vertices.size() >= 3); + auto vit = vertices.begin(); + const auto& p0 = mesh.point(*vit); ++vit; + const auto& p1 = mesh.point(*vit); ++vit; + const auto& p2 = mesh.point(*vit); + return Triangle_3(p0, p1, p2); +} + +void compute_realizing_triangles( + const Surface_mesh& mesh1, const Surface_mesh& mesh2, + const double error_bound, const std::string prefix, const bool save = false) { + + std::cout << "* getting realizing triangles: " << std::endl; + std::vector< std::pair > fpairs; + const double hdist = + PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound, + CGAL::parameters::output_iterator(std::back_inserter(fpairs))); + assert(fpairs.size() == 2); // lower bound face pair + upper bound face pair + const int f1 = static_cast(fpairs.back().first); // f1 / f2: upper bound + const int f2 = static_cast(fpairs.back().second); + + std::cout << "* Hausdorff: " << hdist << std::endl; + std::cout << "* f1 / f2: " << f1 << " / " << f2 << std::endl; + + assert(f1 == 0); + assert(f2 == 0 || f2 == 161); + + if (f1 != -1 && save) { + const auto triangle = get_triangle(f1, mesh1); + Surface_mesh sm1; + sm1.add_vertex(triangle[0]); + sm1.add_vertex(triangle[1]); + sm1.add_vertex(triangle[2]); + sm1.add_face(sm1.vertices()); + save_mesh(sm1, prefix + "triangle1"); + } + + if (f2 != -1 && save) { + const auto triangle = get_triangle(f2, mesh2); + Surface_mesh sm2; + sm2.add_vertex(triangle[0]); + sm2.add_vertex(triangle[1]); + sm2.add_vertex(triangle[2]); + sm2.add_face(sm2.vertices()); + save_mesh(sm2, prefix + "triangle2"); + } +} + +void test_realizing_triangles( + const double error_bound, const bool save = false) { + + std::cout.precision(20); + std::cout << std::endl << "-- test realizing triangles:" << std::endl << std::endl; + + // Basic test. + std::cout << " ---- basic test ---- " << std::endl; + Surface_mesh mesh1, mesh2; + + mesh1.add_vertex(Point_3(0, 0, 0)); + mesh1.add_vertex(Point_3(2, 0, 0)); + mesh1.add_vertex(Point_3(1, 1, 0)); + mesh1.add_face(mesh1.vertices()); + + mesh2.add_vertex(Point_3(0, 0, error_bound / 2.0)); + mesh2.add_vertex(Point_3(2, 0, error_bound / 2.0)); + mesh2.add_vertex(Point_3(1, 1, error_bound / 2.0)); + mesh2.add_face(mesh2.vertices()); + + if (save) save_mesh(mesh1, "1mesh1"); + if (save) save_mesh(mesh2, "1mesh2"); + + compute_realizing_triangles(mesh1, mesh2, error_bound, "1", save); + + mesh2.clear(); + mesh2.add_vertex(Point_3(0, 0, 1)); + mesh2.add_vertex(Point_3(2, 0, 1)); + mesh2.add_vertex(Point_3(1, 1, 1)); + mesh2.add_face(mesh2.vertices()); + + if (save) save_mesh(mesh2, "1mesh2"); + compute_realizing_triangles(mesh1, mesh2, error_bound, "1", save); + + // Complex test. + std::cout << std::endl << " ---- complex test ---- " << std::endl; + mesh1.clear(); + mesh2.clear(); + const std::string filepath1 = "data/tetrahedron.off"; + const std::string filepath2 = "data/tetrahedron-remeshed.off"; + + std::array vhs1 = { + mesh1.add_vertex(Point_3(0, 1, 3)), + mesh1.add_vertex(Point_3(1, 1, 3)), + mesh1.add_vertex(Point_3(1, 0, 3)) + }; + mesh1.add_face(vhs1); + + std::array vhs2 = { + mesh2.add_vertex(Point_3(0, 1, 3.5)), + mesh2.add_vertex(Point_3(1, 1, 3.5)), + mesh2.add_vertex(Point_3(1, 0, 3.5)) + }; + mesh2.add_face(vhs2); + + Surface_mesh tmp1,tmp2; + get_meshes(filepath1, filepath2, tmp1, tmp2); + + PMP::transform(Affine_transformation_3( + CGAL::Translation(), Vector_3(0, 0, 10 * error_bound)), tmp2); + + + + mesh1.join(tmp1); + mesh2.join(tmp2); + + if (save) save_mesh(mesh1, "2mesh1"); + if (save) save_mesh(mesh2, "2mesh2"); + + compute_realizing_triangles(mesh1, mesh2, error_bound, "2", save); +} + +#if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) +void test_parallel_version( + const std::string filepath, const double error_bound) { + + std::cout.precision(20); + std::cout << std::endl << "-- test parallel version:" << std::endl << std::endl; + + Timer timer; + Surface_mesh mesh1, mesh2; + get_meshes(filepath, filepath, mesh1, mesh2); + + // One-sided distance. + std::cout << " ---- one-sided distance ---- " << std::endl; + + PMP::transform(Affine_transformation_3(CGAL::Translation(), + Vector_3(FT(0), FT(0), FT(1))), mesh2); + + std::cout << " ---- SEQUENTIAL ---- " << std::endl; + timer.reset(); + timer.start(); + const double dista = PMP::bounded_error_Hausdorff_distance( + mesh1, mesh2, error_bound, + CGAL::parameters::match_faces(false), + CGAL::parameters::match_faces(false)); + timer.stop(); + const double timea = timer.time(); + + std::cout << " ---- PARALLEL ---- " << std::endl; + timer.reset(); + timer.start(); + const double distb = PMP::bounded_error_Hausdorff_distance( + mesh1, mesh2, error_bound, + CGAL::parameters::match_faces(false), + CGAL::parameters::match_faces(false)); + timer.stop(); + const double timeb = timer.time(); + + std::cout << "* time a seq (sec.): " << timea << std::endl; + std::cout << "* time b par (sec.): " << timeb << std::endl; + + std::cout << "* dista seq = " << dista << std::endl; + std::cout << "* distb par = " << distb << std::endl; + + assert(timea > 0.0); + assert(timeb > 0.0); + assert(dista == distb); +} +#endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) + +void test_early_quit(const std::string filepath) { + + std::cout.precision(20); + std::cout << std::endl << "-- test early quit:" << std::endl << std::endl; + + Timer timer; + Surface_mesh mesh1, mesh2; + get_meshes(filepath, filepath, mesh1, mesh2); + + std::cout << std::endl; + std::cout << " ---- distance 0.0 = 0.0 ---- " << std::endl; + timer.reset(); + timer.start(); + assert(!PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 0.0)); + timer.stop(); + const double timea = timer.time(); + + PMP::transform(Affine_transformation_3(CGAL::Translation(), + Vector_3(0.0, 0.0, 0.5)), mesh2); + + std::cout << " ---- distance 0.5 < 1.0 ---- " << std::endl; + timer.reset(); + timer.start(); + assert(!PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 1.0)); + timer.stop(); + const double timeb = timer.time(); + std::cout << " ---- distance 0.5 < 0.6 ---- " << std::endl; + timer.reset(); + timer.start(); + assert(!PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 0.6)); + timer.stop(); + const double timec = timer.time(); + std::cout << " ---- distance 0.5 > 0.4 ---- " << std::endl; + timer.reset(); + timer.start(); + assert(PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 0.4)); + timer.stop(); + const double timed = timer.time(); + std::cout << " ---- distance 0.5 > 0.0 ---- " << std::endl; + timer.reset(); + timer.start(); + assert(PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 0.0)); + timer.stop(); + const double timee = timer.time(); + + std::cout << "* timea 0.0 = 0.0 = " << timea << std::endl; + std::cout << "* timeb 0.5 < 1.0 = " << timeb << std::endl; + std::cout << "* timec 0.5 < 0.6 = " << timec << std::endl; + std::cout << "* timed 0.5 > 0.4 = " << timed << std::endl; + std::cout << "* timee 0.5 > 0.0 = " << timee << std::endl; + + assert(timea > 0.0); + assert(timeb > 0.0); + assert(timec > 0.0); + assert(timed > 0.0); + assert(timee > 0.0); +} + +void run_examples(const double error_bound, const std::string filepath) { + + remeshing_tetrahedon_example(error_bound); + interior_triangle_example(error_bound); + perturbing_surface_mesh_example(filepath, error_bound); + perturbing_polyhedron_mesh_example(filepath, error_bound); + moving_surface_mesh_example(filepath, filepath, 5, error_bound); +} + +int main(int argc, char** argv) { + + // std::string name; + // std::cin >> name; + + const double error_bound = 1e-4; + const double num_samples = 10.; + std::cout << std::endl << "* error bound: " << error_bound << std::endl; + // std::cout << std::endl << "* number of samples: " << num_samples << std::endl; + std::string filepath = (argc > 1 ? argv[1] : "data/blobby.off"); + run_examples(error_bound, filepath); + + // ------------------------------------------------------------------------ // + // Tests. + + // Approximate_hd_wrapper does not work with EPECK! + Approximate_hd_wrapper apprx_hd(num_samples); + Naive_bounded_error_hd_wrapper naive_hd(error_bound); + Bounded_error_hd_wrapper bound_hd(error_bound); + + // --- Testing basic properties. + + // test_synthetic_data(apprx_hd); + // test_synthetic_data(naive_hd); + test_synthetic_data(bound_hd); + + // --- Compare on common meshes. + + // test_one_versus_another(apprx_hd, naive_hd); + // test_one_versus_another(naive_hd, bound_hd); + test_one_versus_another(bound_hd, apprx_hd); + + // --- Compare on real meshes. + + const std::string filepath1 = (argc > 1 ? argv[1] : "data/blobby.off"); + const std::string filepath2 = (argc > 2 ? argv[2] : "data/tetrahedron-remeshed.off"); + + // test_real_meshes(filepath1, filepath2, apprx_hd, naive_hd); + // test_real_meshes(filepath1, filepath2, naive_hd, bound_hd); + test_real_meshes(filepath1, filepath2, bound_hd, apprx_hd); + + // --- Compare timings. + + filepath = (argc > 1 ? argv[1] : "data/blobby.off"); + // test_timings(filepath, apprx_hd); + // test_timings(filepath, naive_hd); + test_timings(filepath, bound_hd); + + // --- Compare with the paper. + + // test_bunny(apprx_hd); + // test_bunny(naive_hd); + // test_bunny(bound_hd, 3); + + // --- Test realizing triangles. + test_realizing_triangles(error_bound); + + #if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) + // --- Test parallelization. + test_parallel_version(filepath, error_bound); + #endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) + + // --- Test early quit. + test_early_quit(filepath); + + // ------------------------------------------------------------------------ // + + std::cout << std::endl; + return EXIT_SUCCESS; +}