[Small Feature] Add bounded error Hausdorff distance

AABB_tree/include/CGAL/AABB_tree.h

@@ -527,6 +527,36 @@ namespace CGAL {
       }
     }
 
+    template <class Query, class Traversal_traits>
+    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<Traversal_traits,Query>(query, traits, m_primitives.size());
+      }
+    }
+
+    template <class Query, class Traversal_traits>
+    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<AABBTraits> Node;
 

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<class Traversal_traits, class Query>
+  void traversal_with_priority(const Query& query,
+                               Traversal_traits& traits,
+                               const std::size_t nb_primitives) const;
+
+  template<class Primitive_vector, class Traversal_traits, class Query>
+  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<AABB_traits> Node;
@@ -152,6 +166,156 @@ AABB_node<Tr>::traversal(const Query& query,
   }
 }
 
+template<typename Tr>
+template<class Traversal_traits, class Query>
+void
+AABB_node<Tr>::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<typename Tr>
+template<class Primitive_vector, class Traversal_traits, class Query>
+void
+AABB_node<Tr>::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

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)

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}
 }

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

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

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

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.

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
 */

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)