tested group traversal, its sequential version is slower due to more candidate triangles

AABB_tree/include/CGAL/AABB_tree.h

@@ -543,7 +543,7 @@ public:
     }
 
     template <class Query, class Traversal_traits>
-    void traversal_with_priority_and_group_traversal(const Query& query, Traversal_traits& traits, double group_traversal_bound) const
+    void traversal_with_priority_and_group_traversal(const Query& query, Traversal_traits& traits, const std::size_t group_traversal_bound) const
     {
       switch(size())
       {
@@ -553,8 +553,7 @@ public:
         traits.intersection(query, singleton_data());
         break;
       default: // if(size() >= 2)
-        if ( traits.do_intersect(query, *root_node()) )
+        root_node()->template traversal_with_priority_and_group_traversal(m_primitives, query, traits, m_primitives.size(), 0, group_traversal_bound);
-          root_node()->template traversal_with_priority_and_group_traversal(m_primitives, query, traits, m_primitives.size(), 0, group_traversal_bound);
       }
     }
 

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

@@ -234,8 +234,7 @@ AABB_node<Tr>::traversal_with_priority(const Query& query,
   }
 }
 
-
+// TODO: find a better name
-//TODO: find a better name
 template<typename Tr>
 template<class Primitive_vector, class Traversal_traits, class Query>
 void
@@ -244,10 +243,14 @@ AABB_node<Tr>::traversal_with_priority_and_group_traversal(const Primitive_vecto
                                                            Traversal_traits& traits,
                                                            const std::size_t nb_primitives,
                                                            std::size_t first_primitive_index,
-                                                           const std::size_t group_size_bound) const
+                                                           const std::size_t group_traversal_bound) const
 {
-  if (nb_primitives <= group_size_bound)
+  // 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;
   }
@@ -266,7 +269,7 @@ AABB_node<Tr>::traversal_with_priority_and_group_traversal(const Primitive_vecto
     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_size_bound);
+      right_child().traversal_with_priority_and_group_traversal(primitives, query, traits, 2, first_primitive_index+1, group_traversal_bound);
     }
     break;
   default:
@@ -284,22 +287,22 @@ AABB_node<Tr>::traversal_with_priority_and_group_traversal(const Primitive_vecto
         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_size_bound);
+          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_size_bound);
+            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_size_bound);
+          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_size_bound);
+            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_size_bound);
+        left_child().traversal_with_priority_and_group_traversal(primitives, query, traits, nb_primitives/2, first_primitive_index, group_traversal_bound);
       }
     }
     else
@@ -307,7 +310,7 @@ AABB_node<Tr>::traversal_with_priority_and_group_traversal(const Primitive_vecto
       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_size_bound);
+        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);
       }
     }
   }

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h

@@ -1366,6 +1366,7 @@ bounded_error_Hausdorff_impl(
   std::cout.precision(20);
   std::vector<Face_handle> tm1_only, tm2_only;
   std::vector< std::pair<Face_handle, Face_handle> > common;
+  const std::size_t group_traversal_bound = 1; // TODO: Should we put it in the NP?
 
   const auto faces1 = faces(tm1);
   const auto faces2 = faces(tm2);
@@ -1408,26 +1409,35 @@ bounded_error_Hausdorff_impl(
 
   // Build traversal traits for tm1_tree.
   TM1_hd_traits traversal_traits_tm1(
-    tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2, error_bound);
+    tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2, error_bound, group_traversal_bound);
 
   // We skip computing internal Kd tree because we do not compute distance queries.
   tm1_tree.do_not_accelerate_distance_queries();
 
-  // TODO: Initialize the distances on all the vertices first and store those.
   if (false) { // inexact check
     FT max_dist = -FT(1);
-    std::pair<Face_handle, Face_handle> fpair;
+
-    for (const auto& face1 : faces1) {
+    // super slow version with all distances
-      const FT dist = traversal_traits_tm1.get_maximum_distance(face1, fpair);
+    if (false) {
-      CGAL_assertion(fpair.first == face1);
+      std::pair<Face_handle, Face_handle> stub;
-      max_dist = (CGAL::max)(max_dist, dist);
+      for (const auto& face1 : faces1) {
+        const FT dist = traversal_traits_tm1.get_maximum_distance(face1, stub);
+        max_dist = (CGAL::max)(max_dist, dist);
+      }
+    }
+
+    // fast version with bboxes only
+    if (true) {
+      CGAL_assertion_msg(false, "TODO: ADD BBOX COMPARISON!");
     }
 
     if (max_dist <= error_bound) {
       timer.stop();
       // std::cout << "* culling rate: 100%" << std::endl;
       // std::cout << "* preprocessing time (sec.): " << timer.time() << std::endl;
+      std::pair<Face_handle, Face_handle> fpair;
       traversal_traits_tm1.get_maximum_distance(*(faces1.begin()), fpair);
+      CGAL_assertion(fpair.first == *(faces1.begin()));
       return std::make_pair(error_bound, fpair);
     }
   }
@@ -1439,8 +1449,11 @@ bounded_error_Hausdorff_impl(
   timer.reset();
   timer.start();
   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);
+  if (group_traversal_bound > 1) {
-  //tm1_tree.traversal_with_priority_and_group_traversal(stub, traversal_traits_tm1, 10); //TODO: experiment with group traversal
+    tm1_tree.traversal_with_priority_and_group_traversal(stub, traversal_traits_tm1, group_traversal_bound); // with group traversal
+  } else {
+    tm1_tree.traversal_with_priority(stub, traversal_traits_tm1);
+  }
   timer.stop();
   // std::cout << "* TM1 traversal (sec.): " << timer.time() << std::endl;
 
@@ -1449,7 +1462,7 @@ bounded_error_Hausdorff_impl(
   auto global_bounds = traversal_traits_tm1.get_global_bounds();
 
   // std::cout << "* number of all triangles for TM1: " << tm1_tree.size() << std::endl;
-  // std::cout << "* number of candidate triangles before: " << candidate_triangles.size() << std::endl;
+  std::cout << "* number of candidate triangles before: " << candidate_triangles.size() << std::endl;
 
   // std::cout << "* culling rate: " <<
   //   FT(100) - (FT(candidate_triangles.size()) / FT(tm1_tree.size()) * FT(100)) << "%" << std::endl;
@@ -1533,8 +1546,12 @@ bounded_error_Hausdorff_impl(
           infinity_value<FT>(),
           infinity_value<FT>(),
           infinity_value<FT>());
-        tm2_tree.traversal_with_priority(sub_triangles[i], traversal_traits_tm2);
+
-        //tm2_tree.traversal_with_priority_and_group_traversal(sub_triangles[i], traversal_traits_tm2, 10); //TODO: experiment with group traversal
+        if (group_traversal_bound > 1) {
+          tm2_tree.traversal_with_priority_and_group_traversal(sub_triangles[i], traversal_traits_tm2, group_traversal_bound); // with group traversal
+        } else {
+          tm2_tree.traversal_with_priority(sub_triangles[i], traversal_traits_tm2); // without group traversal
+        }
 
         // Update global lower Hausdorff bound according to the obtained local bounds.
         const auto local_bounds = traversal_traits_tm2.get_local_bounds();

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_Hausdorff_distance.h

@@ -235,11 +235,11 @@ namespace CGAL {
       return h_local_bounds;
     }
 
-    template <class PrimitiveConstIterator>
+    template<class PrimitiveConstIterator>
-    void traverse_group(const Query& query, PrimitiveConstIterator group_begin, PrimitiveConstIterator group_end)
+    void traverse_group(const Query& query, PrimitiveConstIterator group_begin, PrimitiveConstIterator group_end) {
-    {
+      for (PrimitiveConstIterator it = group_begin; it != group_end; ++it) {
-      for (PrimitiveConstIterator it=group_begin; it!=group_end; ++it)
         this->intersection(query, *it);
+      }
     }
 
   private:
@@ -287,13 +287,15 @@ namespace CGAL {
       const AABBTraits& traits, const TM2_tree& tree,
       const TriangleMesh& tm1, const TriangleMesh& tm2,
       const VPM1& vpm1, const VPM2& vpm2,
-      const FT error_bound) :
+      const FT error_bound,
+      const std::size_t group_traversal_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_error_bound(error_bound),
+    m_group_traversal_bound(group_traversal_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!
@@ -332,7 +334,12 @@ namespace CGAL {
         infinity_value<FT>());
 
       const Triangle_3 triangle = get(m_face_to_triangle_map, fpair.first);
-      m_tm2_tree.traversal_with_priority(triangle, traversal_traits_tm2);
+
+      if (m_group_traversal_bound > 1) {
+        m_tm2_tree.traversal_with_priority_and_group_traversal(triangle, traversal_traits_tm2, m_group_traversal_bound); // with group traversal
+      } else {
+        m_tm2_tree.traversal_with_priority(triangle, traversal_traits_tm2); // without group traversal
+      }
 
       // Update global Hausdorff bounds according to the obtained local bounds.
       const auto local_bounds = traversal_traits_tm2.get_local_bounds();
@@ -404,11 +411,11 @@ namespace CGAL {
       return this->do_intersect_with_priority(query, node).first;
     }
 
-    template <class PrimitiveConstIterator>
+    template<class PrimitiveConstIterator>
-    void traverse_group(const Query& query, PrimitiveConstIterator group_begin, PrimitiveConstIterator group_end)
+    void traverse_group(const Query& query, PrimitiveConstIterator group_begin, PrimitiveConstIterator group_end) {
-    {
+      for (PrimitiveConstIterator it = group_begin; it != group_end; ++it) {
-      for (PrimitiveConstIterator it=group_begin; it!=group_end; ++it)
         this->intersection(query, *it);
+      }
     }
 
     // Return those triangles from TM1, which are candidates for including a
@@ -465,6 +472,7 @@ namespace CGAL {
 
     // Global Hausdorff bounds for the query triangle.
     const FT m_error_bound;
+    const std::size_t m_group_traversal_bound;
     Global_bounds h_global_bounds;
 
     // All candidate triangles.

Polygon_mesh_processing/test/Polygon_mesh_processing/test_hausdorff_bounded_error_distance.cpp

@@ -732,14 +732,14 @@ int main(int argc, char** argv) {
 
   // test_synthetic_data(apprx_hd);
   // test_synthetic_data(naive_hd);
-  test_synthetic_data(bound_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);
+  // test_one_versus_another(bound_hd, apprx_hd);
 
 
   // --- Compare on real meshes.
@@ -749,7 +749,7 @@ int main(int argc, char** argv) {
 
   // 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);
+  // test_real_meshes(filepath1, filepath2, bound_hd, apprx_hd);
 
 
   // --- Compare timings.
@@ -757,7 +757,7 @@ int main(int argc, char** argv) {
   filepath = (argc > 1 ? argv[1] : "data/blobby-remeshed.off");
   // test_timings(filepath, apprx_hd);
   // test_timings(filepath, naive_hd);
-  test_timings(filepath, bound_hd);
+  // test_timings(filepath, bound_hd);
 
 
   // --- Compare with the paper.
@@ -768,7 +768,7 @@ int main(int argc, char** argv) {
 
 
   // --- Test realizing triangles.
-  test_realizing_triangles(error_bound);
+  // test_realizing_triangles(error_bound);
 
 
   // ------------------------------------------------------------------------ //