Gsoc2019 pmph dist martinskrodzki

AABB_tree/include/CGAL/AABB_tree.h

@@ -520,6 +520,22 @@ 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());
+			}
+		}
+
 	private:
 		typedef AABB_node<AABBTraits> Node;
 

AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h

@@ -86,6 +86,11 @@ 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;
+
 private:
   typedef AABBTraits AABB_traits;
   typedef AABB_node<AABB_traits> Node;
@@ -201,6 +206,67 @@ 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());
+
+    if (ileft)
+    {
+      if (iright)
+      {
+        // Both children have to be inspected
+        if(pleft >= pright)
+        {
+          // Inspect left first, has higher priority
+          left_child().traversal_with_priority(query, traits, nb_primitives/2);
+          if( traits.go_further() ) //&& traits.do_intersect(query, right_child()) ) // TODO shall we call again do_intersect? (Benchmarks show it slows down the Hausdorff Distance)
+            right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2);
+        }
+        else
+        {
+          // Inspect right first, has higher priority
+          right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2);
+          if( traits.go_further() ) //&& traits.do_intersect(query, left_child()) ) // TODO shall we call again do_intersect? (Benchmarks show it slows down the Hausdorff Distance)
+            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);
+  }
+}
+
 } // end namespace CGAL
 
 #endif // CGAL_AABB_NODE_H

Documentation/doc/biblio/geom.bib

@@ -152016,3 +152016,14 @@ @inproceedings{boissonnat2009Delaunay
   HAL_ID = {inria-00412437},
   HAL_VERSION = {v1},
 }
+
+@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}
+}

Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt

@@ -176,6 +176,7 @@ and provides a list of the parameters that are used in this package.
 - \link measure_grp `CGAL::Polygon_mesh_processing::centroid()` \endlink
 
 \cgalCRPSection{Distance Functions}
+- `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()`
 - `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

@@ -729,10 +729,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 (Approximate) 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}.
@@ -757,12 +759,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,
@@ -771,6 +773,40 @@ 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, an AABB tree is
+built on tm1 and tm2 respectively. The tree 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 tree 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 tree 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.
+
+\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

@@ -23,4 +23,5 @@
 \example Polygon_mesh_processing/detect_features_example.cpp
 \example Polygon_mesh_processing/manifoldness_repair_example.cpp
 \example Polygon_mesh_processing/repair_polygon_soup_example.cpp
+\example Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp
 */

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -20,7 +20,7 @@ find_package( CGAL QUIET COMPONENTS  )
 if ( NOT CGAL_FOUND )
 
   message(STATUS "This project requires the CGAL library, and will not be compiled.")
-  return()  
+  return()
 
 endif()
 
@@ -31,7 +31,7 @@ if ( NOT Boost_FOUND )
 
   message(STATUS "This project requires the Boost library, and will not be compiled.")
 
-  return()  
+  return()
 
 endif()
 
@@ -55,8 +55,10 @@ find_package(Eigen3 3.2.0) #(requires 3.2.0 or greater)
 
 find_package( TBB )
 create_single_source_cgal_program( "hausdorff_distance_remeshing_example.cpp")
+create_single_source_cgal_program( "hausdorff_bounded_error_distance_example.cpp")
 if( TBB_FOUND )
   CGAL_target_use_TBB(hausdorff_distance_remeshing_example)
+  CGAL_target_use_TBB(hausdorff_bounded_error_distance_example)
 else()
   message( STATUS "NOTICE: Intel TBB was not found. hausdorff_distance_example will use sequential code." )
 endif()
@@ -118,6 +120,3 @@ target_link_libraries( stitch_borders_example_OM PRIVATE ${OPENMESH_LIBRARIES} )
 create_single_source_cgal_program( "triangulate_faces_example_OM.cpp")
 target_link_libraries( triangulate_faces_example_OM PRIVATE ${OPENMESH_LIBRARIES} )
 endif(OpenMesh_FOUND)
-
-
-

Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp

@@ -0,0 +1,151 @@
+#include <CGAL/Aff_transformation_3.h>
+#include <CGAL/aff_transformation_tags.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+#include <CGAL/Polygon_mesh_processing/remesh.h>
+#include <CGAL/Polygon_mesh_processing/random_perturbation.h>
+#include <CGAL/Polygon_mesh_processing/transform.h>
+#include <CGAL/Real_timer.h>
+
+#if defined(CGAL_LINKED_WITH_TBB)
+#define TAG CGAL::Parallel_tag
+#else
+#define TAG CGAL::Sequential_tag
+#endif
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::Point_3 Point_3;
+typedef K::Triangle_3 Triangle_3;
+typedef K::Vector_3 Vector_3;
+typedef CGAL::Surface_mesh<Point_3> Mesh;
+typedef Mesh::Vertex_index Vertex_index;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char** argv)
+{
+  // Objects to hold the meshes
+  Mesh tm1, tm2;
+
+  // Timer to measure the runtime of h-distance Distance_computation
+  CGAL::Real_timer time;
+
+  // Set an error bound
+  double error_bound = 0.0001;
+
+// ----------------------------------------------------------------------------
+
+  // Easy Example of a tetrahedron and a remeshed version of itself
+
+  // Create the Tetrahedron
+  CGAL::make_tetrahedron(Point_3(.0,.0,.0),
+                         Point_3(2,.0,.0),
+                         Point_3(1,1,1),
+                         Point_3(1,.0,2),
+                         tm1);
+  // Copy it and remesh it
+  tm2=tm1;
+  PMP::isotropic_remeshing(tm2.faces(),.05, tm2);
+  // Compute the Hausdorff distance
+  time.reset();
+  time.start();
+  std::cout << "Approximated Hausdorff distance: "
+            << PMP::bounded_error_Hausdorff_distance
+                  <TAG>(tm1, tm2, error_bound)
+            << std::endl;
+  time.stop();
+  std::cout << "Processing took " << time.time() << "s." << std::endl;
+
+// ----------------------------------------------------------------------------
+
+  // Example with point realizing the Hausdorff distance strictly lying in the
+
+  // interior of a triangle
+  tm1 = Mesh();
+  tm1.add_vertex( Point_3(-1.,1.,1.) );
+  tm1.add_vertex( Point_3(0.,-1.,1.) );
+  tm1.add_vertex( Point_3(1.,1.,1.) );
+  tm1.add_face( tm1.vertices() );
+
+  tm2 = Mesh();
+  Vertex_index w0 = tm2.add_vertex( Point_3(-1.,1.,0.) );
+  Vertex_index w1 = tm2.add_vertex( Point_3(0.,-1.,0.) );
+  Vertex_index w2 = tm2.add_vertex( Point_3(1.,1.,0.) );
+  Vertex_index w3 = tm2.add_vertex( Point_3(0.,1.,-1.) );
+  Vertex_index w4 = tm2.add_vertex( Point_3(-0.5,0.,-1.) );
+  Vertex_index w5 = tm2.add_vertex( Point_3(0.5,0.,-1.) );
+  tm2.add_face( w0, w3, w4 );
+  tm2.add_face( w1, w4, w5 );
+  tm2.add_face( w2, w5, w3 );
+
+  // Compute the Hausdorff distance
+  time.reset();
+  time.start();
+  std::cout << "Approximated Hausdorff distance: "
+            << PMP::bounded_error_Hausdorff_distance
+                  <TAG>(tm1, tm2, error_bound)
+            << std::endl;
+  time.stop();
+  std::cout << "Processing took " << time.time() << "s." << std::endl;
+
+// ----------------------------------------------------------------------------
+
+  // Read a real meshes given by the user, perturb it slightly and compute the
+  // Hausdorff distance between the original mesh and its pertubation
+
+  std::ifstream input( argv[1] );
+  input >> tm1;
+  std::cout << "Read a mesh with " << tm1.number_of_faces() << " triangles." << std::endl;
+
+  // Copy the mesh and perturb it slightly
+  tm2 = tm1;
+  bool do_project = false;
+  PMP::random_perturbation( tm2.vertices(), tm2, 0.1, CGAL::parameters::do_project(do_project));
+  std::cout << "Perturbed the input mesh, now computing the Hausdorff distance." << std::endl;
+
+  // Compute the Hausdorff distance
+  time.reset();
+  time.start();
+  std::cout << "Approximated Hausdorff distance: "
+            << PMP::bounded_error_Hausdorff_distance
+                  <TAG>(tm1, tm2, error_bound)
+            << std::endl;
+  time.stop();
+  std::cout << "Processing took " << time.time() << "s." << std::endl;
+
+// ----------------------------------------------------------------------------
+
+  // 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.
+
+  std::ifstream input1( argv[1] );
+  input1 >> tm1;
+  std::cout << "Read a mesh with " << tm1.number_of_faces() << " triangles." << std::endl;
+
+  std::ifstream input2( argv[2] );
+  input2 >> tm2;
+  std::cout << "Read a mesh with " << tm2.number_of_faces() << " triangles." << std::endl;
+
+  CGAL::Bbox_3 bb = PMP::bbox(tm2);
+  double dist = CGAL::approximate_sqrt( Vector_3(bb.xmax() - bb.xmin(), bb.ymax() - bb.ymin(), bb.zmax() - bb.zmin()).squared_length() );
+
+  for (int i=0; i<300; i++) {
+    PMP::transform(
+       CGAL::Aff_transformation_3<K> ( CGAL::Translation(), Vector_3( 0.01*dist, 0.01*dist, 0.01*dist ) ),
+       tm2
+     );
+    time.reset();
+    time.start();
+    std::cout << "Position: " << i << std::endl;
+    std::cout << "Approximated Hausdorff distance: "
+              << PMP::bounded_error_Hausdorff_distance
+                    <TAG>(tm1, tm2, error_bound)
+              << std::endl;
+    time.stop();
+    std::cout << "Processing took " << time.time() << "s." << std::endl;
+  }
+}