From fe21d983456fb9bc81abe6e0bb3dd6e67541e15f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 31 Mar 2022 12:48:18 +0200 Subject: [PATCH 01/14] Improve Hausdorff distance debug code + more assertions --- .../CGAL/Polygon_mesh_processing/distance.h | 47 ++++++++++++++----- ...traversal_traits_with_Hausdorff_distance.h | 47 +++++++++---------- 2 files changed, 58 insertions(+), 36 deletions(-) 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 1e1735ddbbb..1d0e17f7b63 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1445,6 +1445,14 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, using Point_3 = typename Kernel::Point_3; using Triangle_3 = typename Kernel::Triangle_3; +#ifdef CGAL_HAUSDORFF_DEBUG + std::cout << " -- Bounded Hausdorff --" << std::endl; + std::cout << "error bound: " << sq_error_bound << " (" << approximate_sqrt(sq_error_bound) << ")" << std::endl; + std::cout << "initial bound: " << sq_initial_bound << " (" << approximate_sqrt(sq_initial_bound) << ")" << std::endl; + std::cout << "distance bound: " << sq_distance_bound << " (" << approximate_sqrt(sq_distance_bound) << ")" << std::endl; + std::cout << "inf val: " << infinity_value << " (" << approximate_sqrt(infinity_value) << ")" << std::endl; +#endif + using TM1_hd_traits = Hausdorff_primitive_traits_tm1; using TM2_hd_traits = Hausdorff_primitive_traits_tm2; @@ -1453,6 +1461,7 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, using Candidate = Candidate_triangle; + CGAL_precondition(sq_initial_bound >= sq_error_bound); CGAL_precondition(sq_distance_bound != FT(0)); // value is -1 if unused CGAL_precondition(sq_error_bound >= FT(0)); CGAL_precondition(tm1_tree.size() > 0); @@ -1511,8 +1520,8 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, // Second, we apply subdivision. #ifdef CGAL_HAUSDORFF_DEBUG timer.reset(); - timer.start(); std::cout << "- applying subdivision" << std::endl; + timer.start(); std::size_t explored_candidates_count = 0; #endif @@ -1523,14 +1532,17 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, std::cout << "===" << std::endl; std::cout << candidate_triangles.size() << " candidates" << std::endl; std::cout << "- infinity_value: " << infinity_value << std::endl; - std::cout << "- initial_bound: " << sq_initial_bound << std::endl; - std::cout << "- distance_bound: " << sq_distance_bound << std::endl; + std::cout << "- sq_error_bound: " << sq_error_bound << std::endl; + std::cout << "- sq_initial_bound: " << sq_initial_bound << std::endl; + std::cout << "- sq_distance_bound: " << sq_distance_bound << std::endl; std::cout << "- global_bounds.lower: " << global_bounds.lower << std::endl; std::cout << "- global_bounds.upper: " << global_bounds.upper << std::endl; std::cout << "- diff = " << (global_bounds.upper - global_bounds.lower) << ", below bound? " << ((global_bounds.upper - global_bounds.lower) <= sq_error_bound) << std::endl; #endif + CGAL_assertion(global_bounds.upper >= global_bounds.lower); + if((global_bounds.upper - global_bounds.lower <= sq_error_bound)) break; @@ -1563,6 +1575,8 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, std::cout << triangle_and_bounds.triangle.vertex(2) << std::endl; std::cout << "triangle_bounds.lower: " << triangle_bounds.lower << std::endl; std::cout << "triangle_bounds.upper: " << triangle_bounds.upper << std::endl; + std::cout << "- diff = " << (triangle_bounds.upper - triangle_bounds.lower) << ", below bound? " + << ((triangle_bounds.upper - triangle_bounds.lower) <= sq_error_bound) << std::endl; #endif CGAL_assertion(triangle_bounds.lower >= FT(0)); @@ -1574,6 +1588,10 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, if(triangle_bounds.upper < global_bounds.lower) continue; + // The check we want is |d1 - d2| < error_bound, but that would require square roots. + // It's cheaper to require |d1^2 - d2^2| < error_bound^2, which is a sufficient condition: + // without loss of generality, assume that d1 > d2, (d1 - d2)^2 = d1^2 + d2^2 + 2d1d2, + // and d1 and d2 are positive thus if d1^2 - d2^2 < error_bound^2, then (d1 - d2) < error_bound if((triangle_bounds.upper - triangle_bounds.lower) <= sq_error_bound) { #ifdef CGAL_HAUSDORFF_DEBUG_PP @@ -1592,15 +1610,16 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, const Point_3& v1 = triangle_for_subdivision.vertex(1); const Point_3& v2 = triangle_for_subdivision.vertex(2); - // Third stopping condition: all edge lengths of the triangle are smaller than the given error bound - // @todo can we even be here considering it implies the triangle bounds are within the error bound, - // and thus we would have already continue'd? + // Stopping condition: all edge lengths of the triangle are smaller than the given error bound. if(CGAL::squared_distance(v0, v1) < sq_error_bound && CGAL::squared_distance(v0, v2) < sq_error_bound && CGAL::squared_distance(v1, v2) < sq_error_bound) { #ifdef CGAL_HAUSDORFF_DEBUG_PP std::cout << "Third stopping condition, small triangle" << std::endl; + std::cout << CGAL::squared_distance(v0, v1) << " vs " << sq_error_bound << std::endl; + std::cout << CGAL::squared_distance(v0, v2) << " vs " << sq_error_bound << std::endl; + std::cout << CGAL::squared_distance(v1, v2) << " vs " << sq_error_bound << std::endl; #endif // By definition, lower_bound(t1, TM2) is smaller than h(t1, TM2). @@ -1672,12 +1691,15 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, // Update global lower Hausdorff bound according to the obtained local bounds. const auto& sub_triangle_bounds = traversal_traits_tm2.get_local_bounds(); + #ifdef CGAL_HAUSDORFF_DEBUG_PP std::cout << "Subdivided triangle bounds: " << sub_triangle_bounds.lower << " " << sub_triangle_bounds.upper << std::endl; #endif CGAL_assertion(sub_triangle_bounds.lower >= FT(0)); CGAL_assertion(sub_triangle_bounds.upper >= sub_triangle_bounds.lower); + CGAL_assertion(sub_triangle_bounds.tm2_lface != boost::graph_traits::null_face()); + CGAL_assertion(sub_triangle_bounds.tm2_uface != boost::graph_traits::null_face()); // The global lower bound is the max of the per-face lower bounds if(sub_triangle_bounds.lower > global_bounds.lower) @@ -1700,18 +1722,16 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, const Candidate& top_candidate = candidate_triangles.top(); const FT current_upmost = top_candidate.bounds.upper; #ifdef CGAL_HAUSDORFF_DEBUG_PP - std::cout << "current candidates count: " << candidate_triangles.size() << std::endl; std::cout << "current upper bound = " << current_upmost << std::endl; + std::cout << "global_bounds.lower = " << global_bounds.lower << std::endl; #endif + CGAL_assertion(is_positive(current_upmost)); - // Below can happen if the subtriangle returned something like [l;u], - // with l and u both below the global error bound. The lowest bound will not have been updated - // since it has been initialized with the error bound. if(current_upmost < global_bounds.lower) { #ifdef CGAL_HAUSDORFF_DEBUG_PP - std::cout << "upmost is below lowest, end." << std::endl; + std::cout << "Top of the queue is lower than the lowest!" << std::endl; #endif global_bounds.upper = global_bounds.lower; // not really needed since lower is returned but doesn't hurt @@ -1725,6 +1745,8 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, break; } + CGAL_assertion(current_upmost >= global_bounds.lower); + global_bounds.upper = current_upmost; global_bounds.upair.first = top_candidate.tm1_face; global_bounds.upair.second = top_candidate.bounds.tm2_uface; @@ -1741,6 +1763,7 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, timer.stop(); std::cout << "* subdivision (sec.): " << timer.time() << std::endl; std::cout << "Explored " << explored_candidates_count << " candidates" << std::endl; + std::cout << "Final global bounds: " << global_bounds.lower << " " << global_bounds.upper << std::endl; #endif // Get realizing triangles. @@ -2196,7 +2219,7 @@ bounded_error_squared_one_sided_Hausdorff_distance_impl(const TriangleMesh1& tm1 #endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) { #ifdef CGAL_HAUSDORFF_DEBUG - std::cout << "* executing sequential version " << std::endl; + std::cout << "* executing sequential version" << std::endl; #endif sq_hdist = bounded_error_squared_Hausdorff_distance_impl( tm1, tm2, vpm1, vpm2, tm1_tree, tm2_tree, 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 index 08e0e6e99a9..d60f4523d0a 100644 --- 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 @@ -25,12 +25,6 @@ #include #include -#ifdef CGAL_HAUSDORFF_DEBUG_PP - #ifndef CGAL_HAUSDORFF_DEBUG - #define CGAL_HAUSDORFF_DEBUG - #endif -#endif - namespace CGAL { namespace Polygon_mesh_processing { namespace internal { @@ -188,11 +182,14 @@ public: if(m_early_exit) return; -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL std::cout << "Intersection with TM2's " << primitive.id() << std::endl; std::cout << "Initial local bounds " << m_local_bounds.lower << " " << m_local_bounds.upper << std::endl; #endif + CGAL_assertion(m_local_bounds.lower >= FT(0)); + CGAL_assertion(m_local_bounds.upper >= FT(0)); + /* Have reached a single triangle, process it. / Determine the upper distance according to / min_{b \in primitive} ( max_{vertex in query} ( d(vertex, b) ) ) @@ -210,7 +207,7 @@ public: CGAL_assertion(primitive.id() != Face_handle_2()); const Triangle_3 triangle = get(m_face_to_triangle_map, primitive.id()); -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL std::cout << "Geometry: " << triangle << std::endl; #endif @@ -238,7 +235,7 @@ public: // h_lower_i(query, TM2) := max_{v in query} min_{1<=j<=i} d(v, primitive_j) const FT distance_lower = (CGAL::max)((CGAL::max)(m_v0_lower, m_v1_lower), m_v2_lower); -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL std::cout << "Distance from vertices of t1 to t2: " << v0_dist << " " << v1_dist << " " << v2_dist << std::endl; #endif @@ -247,10 +244,9 @@ public: // With each new TM2 face, the min value m_v{k}_lower can become smaller, // and thus also the value max_{v in query} min_{1<=j<=i} d(v, primitive_j) - CGAL_assertion(m_local_bounds.lower >= FT(0)); if(distance_lower < m_local_bounds.lower) { -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL std::cout << "new best lower (" << distance_lower << ") with TM2 face: " << triangle << std::endl; #endif m_local_bounds.lower = distance_lower; @@ -259,20 +255,21 @@ public: // This is the 'min_{1<=j<=i}' part in: // h_upper_i(query, TM2) = min_{1<=j<=i} max_{v in query} (v, primitive_j), Equation (10) - CGAL_assertion(m_local_bounds.upper >= FT(0)); if(distance_upper < m_local_bounds.upper) { -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL std::cout << "new best upper (" << distance_upper << ") with TM2 face: " << triangle << std::endl; #endif m_local_bounds.upper = distance_upper; m_local_bounds.tm2_uface = primitive.id(); } -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL std::cout << "Distance from vertices of t1 to t2: " << v0_dist << " " << v1_dist << " " << v2_dist << std::endl; std::cout << "Current local bounds " << m_local_bounds.lower << " " << m_local_bounds.upper << std::endl; #endif + + CGAL_assertion(m_local_bounds.lower >= FT(0)); CGAL_assertion(m_local_bounds.lower <= m_local_bounds.upper); // #define CGAL_PMP_HDIST_NO_CULLING_DURING_TRAVERSAL @@ -281,8 +278,8 @@ public: // whereas the rhs can only go up with every additional TM1 face if(m_local_bounds.upper < m_global_bounds.lower) // Section 4.1, first ยง { -#ifdef CGAL_HAUSDORFF_DEBUG_PP - std::cout << "Quitting early (TM2 traversal): " << m_global_bounds.lower << std::endl; +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL + std::cout << "Quitting early (TM2 traversal), global lower " << m_global_bounds.lower << " greater than local upper " << m_local_bounds.upper << std::endl; #endif m_early_exit = true; } @@ -298,7 +295,7 @@ public: if(m_early_exit) return std::make_pair(false, FT(0)); -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL std::cout << "Do_intersect TM2 node with bbox: " << node.bbox() << std::endl; #endif @@ -333,8 +330,8 @@ public: // between the query and the TM2 primitives that are children of this node. // If this lower bound is greater than the current upper bound for this query, // then none of these primitives will reduce the Hausdorff distance between the query and TM2. -#ifdef CGAL_HAUSDORFF_DEBUG_PP - std::cout << "Culling TM1? dist vs local bound upper " << sq_dist << " " << m_local_bounds.upper << std::endl; +#ifdef CGAL_HAUSDORFF_DEBUG_TM2_TRAVERSAL + std::cout << "Culling TM2? dist vs local bound upper " << sq_dist << " " << m_local_bounds.upper << std::endl; #endif CGAL_assertion(m_local_bounds.upper >= FT(0)); if(sq_dist > m_local_bounds.upper) @@ -544,7 +541,7 @@ public: if(m_early_exit) return; -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM1_TRAVERSAL std::cout << "Intersection with TM1's " << primitive.id() << std::endl; std::cout << "Initial global bounds " << m_global_bounds.lower << " " << m_global_bounds.upper << std::endl; #endif @@ -554,7 +551,7 @@ public: const Face_handle_1 tm1_face = primitive.id(); const Triangle_3 triangle = get(m_face_to_triangle_map, tm1_face); -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM1_TRAVERSAL std::cout << "Geometry: " << triangle << std::endl; #endif @@ -566,12 +563,14 @@ public: // Post traversal, we have computed h_lower(query, TM2) and h_upper(query, TM2) const auto& local_bounds = traversal_traits_tm2.get_local_bounds(); -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM1_TRAVERSAL std::cout << "Bounds for TM1 primitive: " << local_bounds.lower << " " << local_bounds.upper << std::endl; #endif CGAL_assertion(local_bounds.lower >= FT(0)); CGAL_assertion(local_bounds.upper >= local_bounds.lower); + CGAL_assertion(local_bounds.tm2_lface != boost::graph_traits::null_face()); + CGAL_assertion(local_bounds.tm2_uface != boost::graph_traits::null_face()); // Update global Hausdorff bounds according to the obtained local bounds. // h_lower(TM1, TM2) = max_{query in TM1} h_lower(query, TM2) @@ -617,7 +616,7 @@ public: if(m_early_exit) return std::make_pair(false, FT(0)); -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM1_TRAVERSAL std::cout << "Do_intersect TM1 node with bbox " << node.bbox() << std::endl; #endif @@ -640,7 +639,7 @@ public: // If the upper bound is smaller than the current global lower bound, // it is pointless to visit this node (and its children) because a larger distance // has been found somewhere else. -#ifdef CGAL_HAUSDORFF_DEBUG_PP +#ifdef CGAL_HAUSDORFF_DEBUG_TM1_TRAVERSAL std::cout << "Culling TM1? dist & global lower bound: " << sq_dist << " " << m_global_bounds.lower << std::endl; #endif CGAL_assertion(m_global_bounds.lower >= FT(0)); From bd29b976f5050f60b0c70e1a7f368e31095e67b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 31 Mar 2022 12:48:44 +0200 Subject: [PATCH 02/14] Use tighter initialization during computation of subdivided faces' local bounds --- .../CGAL/Polygon_mesh_processing/distance.h | 24 +++++++++---------- 1 file changed, 11 insertions(+), 13 deletions(-) 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 1d0e17f7b63..4594f76c72f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1670,22 +1670,20 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, // (We also have that error_bound is a lower bound.) const Bbox_3 sub_t1_bbox = sub_triangles[i].bbox(); - // Because: + // The lower bound is: // h_lower(t1, TM2) := max_{v in t1} min_{t2 in TM2} d(v, t2) - // adding more vertices can only increase the max (and the lower bound). - // + // The upper bound is: // h_upper(t1, TM2) := min_{t2 in TM2} max_{v in t1} d(v, t2) - // If t2m is the face of TM2 realizing the minimum of max_{v in t1} d(v, t2), - // then subdividing t1 can only decrease this upper bound: let v' be a new vertex v' - // of a triangle subdividing t1 is such that - // min_{t2 in TM2} d(v', t2) > h_upper(t1, TM2) - // Because the maximum of the distance over t1 is necessarily reached at a vertex, - // there must exist a vertex v1 in t1 such that d(v1, t2) > d(v', t2), which contradicts - // the definition of v'. - // Thus, subdivision can only decrease the upper bound. - // @FIXME using 'triangle_bounds' actually create bugs?... - Local_bounds bounds(infinity_value); + // The value max_{p in t1} d(p, t2) is realized at a vertex of t1. + // Thus, when splitting t1 into four subtriangles, the distance at the three new vertices + // is smaller than max_{v in t1} d(v, t2) + // Thus, subdivision can only decrease the min, and the upper bound. + Local_bounds bounds(triangle_bounds.upper); + + // Ensure 'uface' is initialized in case the upper bound is not changed by the subdivision + bounds.tm2_uface = triangle_bounds.tm2_uface; + TM2_hd_traits traversal_traits_tm2(sub_t1_bbox, tm2, vpm2, bounds, global_bounds, infinity_value); tm2_tree.traversal_with_priority(sub_triangles[i], traversal_traits_tm2); From f46d4d1fae85c1abc6883d6d08f7099b8b47b636 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 31 Mar 2022 12:49:31 +0200 Subject: [PATCH 03/14] Do not add triangles that cannot realize the distance to the priority queue --- .../include/CGAL/Polygon_mesh_processing/distance.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) 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 4594f76c72f..f9bd7e82425 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1712,8 +1712,9 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, // which can go down, so it is only recomputed once splitting is finished, // using the top value of the PQ - // Add the subtriangle to the candidate list. - candidate_triangles.emplace(sub_triangles[i], sub_triangle_bounds, triangle_and_bounds.tm1_face); + // Add the subtriangle to the candidate list if it is meaningful + if(sub_triangle_bounds.upper >= global_bounds.lower) + candidate_triangles.emplace(sub_triangles[i], sub_triangle_bounds, triangle_and_bounds.tm1_face); } // Update global upper Hausdorff bound after subdivision. From 14a9abcca6e48bb3ec54bca7260446f07c83bf31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 31 Mar 2022 12:50:05 +0200 Subject: [PATCH 04/14] Fix assertion: global lbound is init. w/ the initial bound, not the error bound --- .../include/CGAL/Polygon_mesh_processing/distance.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) 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 f9bd7e82425..a9ad7c689ed 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1733,14 +1733,15 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, std::cout << "Top of the queue is lower than the lowest!" << std::endl; #endif + // This can happen if the actual Hausdorff distance is smaller than sq_initial_bound, + // which can happen either if the error bound is really large, or if the initial bound + // is non-zero (symmetric Hausdorff distance). In this case, we return 'sq_initial_bound' + CGAL_assertion(global_bounds.lower == sq_initial_bound); + global_bounds.upper = global_bounds.lower; // not really needed since lower is returned but doesn't hurt global_bounds.upair.first = global_bounds.lpair.first; global_bounds.upair.second = global_bounds.lpair.second; - // Current upmost being equal to the lower is fine, but if it's strictly below, it must - // be because we crossed the error bound, or there is some issue... - CGAL_assertion(global_bounds.lower == sq_error_bound); - break; } From 18c5997fe9fa58b24ee8e28d8f0e056195cead2e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 31 Mar 2022 12:52:38 +0200 Subject: [PATCH 05/14] Re-enable Hausdorff checks in Surface_mesh_simplification --- .../edge_collapse_garland_heckbert_variations.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Surface_mesh_simplification/test/Surface_mesh_simplification/edge_collapse_garland_heckbert_variations.cpp b/Surface_mesh_simplification/test/Surface_mesh_simplification/edge_collapse_garland_heckbert_variations.cpp index 52fc0248fda..535254708bb 100644 --- a/Surface_mesh_simplification/test/Surface_mesh_simplification/edge_collapse_garland_heckbert_variations.cpp +++ b/Surface_mesh_simplification/test/Surface_mesh_simplification/edge_collapse_garland_heckbert_variations.cpp @@ -280,9 +280,9 @@ void run(const std::pair& input) time_all_policies(input.first, out); -// std::array range {{ 0.15 }}; + std::array range {{ 0.15 }}; // std::array range {{ 0.7, 0.6, 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15 }}; -// hausdorff_errors(input.first, out, range.begin(), range.end()); + hausdorff_errors(input.first, out, range.begin(), range.end()); gather_face_aspect_ratio(input.first, out); } From 6fcc5f64a5553402700712faf11c4e258cca06a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 4 Apr 2022 12:16:21 +0200 Subject: [PATCH 06/14] Do not .top() an empty queue --- .../include/CGAL/Polygon_mesh_processing/distance.h | 4 ++++ 1 file changed, 4 insertions(+) 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 a9ad7c689ed..e09a7299f9f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1717,6 +1717,10 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, candidate_triangles.emplace(sub_triangles[i], sub_triangle_bounds, triangle_and_bounds.tm1_face); } + // In case all subdividing triangles of the last queue entry are useless + if(candidate_triangles.empty()) + break; + // Update global upper Hausdorff bound after subdivision. const Candidate& top_candidate = candidate_triangles.top(); const FT current_upmost = top_candidate.bounds.upper; From 30e0a5d0212b21c893cdfcdc3eccb6e46cf8eee5 Mon Sep 17 00:00:00 2001 From: Mael Date: Tue, 5 Apr 2022 16:43:29 +0200 Subject: [PATCH 07/14] Restore consistency between bounds and candidates queue --- .../include/CGAL/Polygon_mesh_processing/distance.h | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) 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 e09a7299f9f..552587f0f29 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1712,14 +1712,10 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, // which can go down, so it is only recomputed once splitting is finished, // using the top value of the PQ - // Add the subtriangle to the candidate list if it is meaningful - if(sub_triangle_bounds.upper >= global_bounds.lower) - candidate_triangles.emplace(sub_triangles[i], sub_triangle_bounds, triangle_and_bounds.tm1_face); + candidate_triangles.emplace(sub_triangles[i], sub_triangle_bounds, triangle_and_bounds.tm1_face); } - // In case all subdividing triangles of the last queue entry are useless - if(candidate_triangles.empty()) - break; + CGAL_assertion(!candidate_triangles.empty()); // Update global upper Hausdorff bound after subdivision. const Candidate& top_candidate = candidate_triangles.top(); From aa5fd2e0ce2ff0bda9af1f6aefcdeee8c2b4a8e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 7 Apr 2022 11:40:30 +0200 Subject: [PATCH 08/14] Switch back to non-squared values for comparisons --- .../CGAL/Polygon_mesh_processing/distance.h | 87 +++++++++++-------- ...traversal_traits_with_Hausdorff_distance.h | 7 +- 2 files changed, 52 insertions(+), 42 deletions(-) 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 552587f0f29..cc61ab329cb 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1435,7 +1435,7 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, const VPM2 vpm2, const TM1Tree& tm1_tree, const TM2Tree& tm2_tree, - const typename Kernel::FT sq_error_bound, + const typename Kernel::FT error_bound, const typename Kernel::FT sq_initial_bound, const typename Kernel::FT sq_distance_bound, const typename Kernel::FT infinity_value, @@ -1445,9 +1445,11 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, using Point_3 = typename Kernel::Point_3; using Triangle_3 = typename Kernel::Triangle_3; + const FT sq_error_bound = square(FT(error_bound)); + #ifdef CGAL_HAUSDORFF_DEBUG std::cout << " -- Bounded Hausdorff --" << std::endl; - std::cout << "error bound: " << sq_error_bound << " (" << approximate_sqrt(sq_error_bound) << ")" << std::endl; + std::cout << "error bound: " << error_bound << " (square: " << sq_error_bound << ")" << std::endl; std::cout << "initial bound: " << sq_initial_bound << " (" << approximate_sqrt(sq_initial_bound) << ")" << std::endl; std::cout << "distance bound: " << sq_distance_bound << " (" << approximate_sqrt(sq_distance_bound) << ")" << std::endl; std::cout << "inf val: " << infinity_value << " (" << approximate_sqrt(infinity_value) << ")" << std::endl; @@ -1478,7 +1480,7 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, // Build traversal traits for tm1_tree. TM1_hd_traits traversal_traits_tm1(tm2_tree, tm1, tm2, vpm1, vpm2, - sq_error_bound, infinity_value, sq_initial_bound, sq_distance_bound); + infinity_value, sq_initial_bound, sq_distance_bound); // Find candidate triangles in TM1, which might realise the Hausdorff bound. // We build a sorted structure while collecting the candidates. @@ -1532,18 +1534,22 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, std::cout << "===" << std::endl; std::cout << candidate_triangles.size() << " candidates" << std::endl; std::cout << "- infinity_value: " << infinity_value << std::endl; - std::cout << "- sq_error_bound: " << sq_error_bound << std::endl; + std::cout << "- error_bound: " << error_bound << std::endl; std::cout << "- sq_initial_bound: " << sq_initial_bound << std::endl; std::cout << "- sq_distance_bound: " << sq_distance_bound << std::endl; std::cout << "- global_bounds.lower: " << global_bounds.lower << std::endl; std::cout << "- global_bounds.upper: " << global_bounds.upper << std::endl; - std::cout << "- diff = " << (global_bounds.upper - global_bounds.lower) << ", below bound? " - << ((global_bounds.upper - global_bounds.lower) <= sq_error_bound) << std::endl; + std::cout << "- diff = " << CGAL::approximate_sqrt(global_bounds.upper) - + CGAL::approximate_sqrt(global_bounds.lower) << ", below bound? " + << ((CGAL::approximate_sqrt(global_bounds.upper) - + CGAL::approximate_sqrt(global_bounds.lower)) <= error_bound) << std::endl; #endif + CGAL_assertion(global_bounds.lower >= FT(0)); CGAL_assertion(global_bounds.upper >= global_bounds.lower); - if((global_bounds.upper - global_bounds.lower <= sq_error_bound)) + // @todo could cache those sqrts + if(CGAL::approximate_sqrt(global_bounds.upper) - CGAL::approximate_sqrt(global_bounds.lower) <= error_bound) break; // Check if we can early quit. @@ -1575,8 +1581,10 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, std::cout << triangle_and_bounds.triangle.vertex(2) << std::endl; std::cout << "triangle_bounds.lower: " << triangle_bounds.lower << std::endl; std::cout << "triangle_bounds.upper: " << triangle_bounds.upper << std::endl; - std::cout << "- diff = " << (triangle_bounds.upper - triangle_bounds.lower) << ", below bound? " - << ((triangle_bounds.upper - triangle_bounds.lower) <= sq_error_bound) << std::endl; + std::cout << "- diff = " << CGAL::approximate_sqrt(triangle_bounds.upper) - + CGAL::approximate_sqrt(triangle_bounds.lower) << ", below bound? " + << ((CGAL::approximate_sqrt(triangle_bounds.upper) - + CGAL::approximate_sqrt(triangle_bounds.lower)) <= error_bound) << std::endl; #endif CGAL_assertion(triangle_bounds.lower >= FT(0)); @@ -1586,13 +1594,14 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, // Might have been a good candidate when added to the queue, but rendered useless by later insertions if(triangle_bounds.upper < global_bounds.lower) + { +#ifdef CGAL_HAUSDORFF_DEBUG_PP + std::cout << "Upper bound is lower than global.lower" << std::endl; +#endif continue; + } - // The check we want is |d1 - d2| < error_bound, but that would require square roots. - // It's cheaper to require |d1^2 - d2^2| < error_bound^2, which is a sufficient condition: - // without loss of generality, assume that d1 > d2, (d1 - d2)^2 = d1^2 + d2^2 + 2d1d2, - // and d1 and d2 are positive thus if d1^2 - d2^2 < error_bound^2, then (d1 - d2) < error_bound - if((triangle_bounds.upper - triangle_bounds.lower) <= sq_error_bound) + if((CGAL::approximate_sqrt(triangle_bounds.upper) - CGAL::approximate_sqrt(triangle_bounds.lower)) <= error_bound) { #ifdef CGAL_HAUSDORFF_DEBUG_PP std::cout << "Candidate triangle bounds are tight enough: " << triangle_bounds.lower << " " << triangle_bounds.upper << std::endl; @@ -1764,8 +1773,16 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, std::cout << "* subdivision (sec.): " << timer.time() << std::endl; std::cout << "Explored " << explored_candidates_count << " candidates" << std::endl; std::cout << "Final global bounds: " << global_bounds.lower << " " << global_bounds.upper << std::endl; + std::cout << "Final global bounds (sqrt): " << CGAL::approximate_sqrt(global_bounds.lower) << " " + << CGAL::approximate_sqrt(global_bounds.upper) << std::endl; + std::cout << "Difference: " << CGAL::approximate_sqrt(global_bounds.upper) - + CGAL::approximate_sqrt(global_bounds.lower) << std::endl; #endif + CGAL_assertion(global_bounds.lower >= FT(0)); + CGAL_assertion(global_bounds.upper >= global_bounds.lower); + CGAL_assertion(CGAL::approximate_sqrt(global_bounds.upper) - CGAL::approximate_sqrt(global_bounds.lower) <= error_bound); + // 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()); @@ -1880,7 +1897,7 @@ struct Bounded_error_squared_distance_computation const std::vector& tm1_parts; const TriangleMesh2& tm2; - const FT sq_error_bound; + const double error_bound; const VPM1 vpm1; const VPM2 vpm2; const FT infinity_value; const FT sq_initial_bound; @@ -1891,14 +1908,14 @@ struct Bounded_error_squared_distance_computation // Constructor. Bounded_error_squared_distance_computation(const std::vector& tm1_parts, const TriangleMesh2& tm2, - const FT sq_error_bound, + const double error_bound, const VPM1 vpm1, const VPM2 vpm2, const FT infinity_value, const FT sq_initial_bound, const std::vector& tm1_trees, const TM2Tree& tm2_tree) : tm1_parts(tm1_parts), tm2(tm2), - sq_error_bound(sq_error_bound), + error_bound(error_bound), vpm1(vpm1), vpm2(vpm2), infinity_value(infinity_value), sq_initial_bound(sq_initial_bound), tm1_trees(tm1_trees), tm2_tree(tm2_tree), @@ -1910,7 +1927,7 @@ struct Bounded_error_squared_distance_computation // Split constructor. Bounded_error_squared_distance_computation(Bounded_error_squared_distance_computation& s, tbb::split) : tm1_parts(s.tm1_parts), tm2(s.tm2), - sq_error_bound(s.sq_error_bound), + error_bound(s.error_bound), vpm1(s.vpm1), vpm2(s.vpm2), infinity_value(s.infinity_value), sq_initial_bound(s.sq_initial_bound), tm1_trees(s.tm1_trees), tm2_tree(s.tm2_tree), @@ -1942,7 +1959,7 @@ struct Bounded_error_squared_distance_computation // for checking if two meshes are close. const FT sqd = bounded_error_squared_Hausdorff_distance_impl( tm1, tm2, vpm1, vpm2, tm1_tree, tm2_tree, - sq_error_bound, sq_initial_bound, FT(-1) /*sq_distance_bound*/, infinity_value, + error_bound, sq_initial_bound, FT(-1) /*sq_distance_bound*/, infinity_value, stub); if(sqd > sq_dist) sq_dist = sqd; @@ -1977,7 +1994,7 @@ template FT(0)); - CGAL_assertion(sq_error_bound >= FT(0)); + CGAL_assertion(error_bound >= 0.); - const FT sq_initial_bound = sq_error_bound; + const FT sq_initial_bound = square(FT(error_bound)); FT sq_hdist = FT(-1); #ifdef CGAL_HAUSDORFF_DEBUG @@ -2209,7 +2226,7 @@ bounded_error_squared_one_sided_Hausdorff_distance_impl(const TriangleMesh1& tm1 using Comp = Bounded_error_squared_distance_computation; - Comp bedc(tm1_parts, tm2, sq_error_bound, vpm1, vpm2, + Comp bedc(tm1_parts, tm2, error_bound, vpm1, vpm2, infinity_value, sq_initial_bound, tm1_trees, tm2_tree); tbb::parallel_reduce(tbb::blocked_range(0, tm1_parts.size()), bedc); @@ -2223,7 +2240,7 @@ bounded_error_squared_one_sided_Hausdorff_distance_impl(const TriangleMesh1& tm1 #endif sq_hdist = bounded_error_squared_Hausdorff_distance_impl( tm1, tm2, vpm1, vpm2, tm1_tree, tm2_tree, - sq_error_bound, sq_initial_bound, sq_distance_bound, infinity_value, out); + error_bound, sq_initial_bound, sq_distance_bound, infinity_value, out); } #ifdef CGAL_HAUSDORFF_DEBUG @@ -2251,7 +2268,7 @@ template tm1_only; std::vector tm2_only; + const FT sq_error_bound = square(FT(error_bound)); FT infinity_value = FT(-1); // All trees below are built and/or accelerated lazily. @@ -2312,6 +2330,7 @@ bounded_error_squared_symmetric_Hausdorff_distance_impl(const TriangleMesh1& tm1 return 0.; // TM1 and TM2 are equal so the distance is zero } + CGAL_assertion(is_positive(infinity_value)); // Compute the first one-sided distance. @@ -2322,7 +2341,7 @@ bounded_error_squared_symmetric_Hausdorff_distance_impl(const TriangleMesh1& tm1 { sq_dista = bounded_error_squared_Hausdorff_distance_impl( tm1, tm2, vpm1, vpm2, tm1_tree, tm2_tree, - sq_error_bound, sq_initial_bound, sq_distance_bound, infinity_value, out1); + error_bound, sq_initial_bound, sq_distance_bound, infinity_value, out1); } // In case this is true, we need to rebuild trees in order to accelerate @@ -2346,7 +2365,7 @@ bounded_error_squared_symmetric_Hausdorff_distance_impl(const TriangleMesh1& tm1 { sq_distb = bounded_error_squared_Hausdorff_distance_impl( tm2, tm1, vpm2, vpm1, tm2_tree, tm1_tree, - sq_error_bound, sq_initial_bound, sq_distance_bound, infinity_value, out2); + error_bound, sq_initial_bound, sq_distance_bound, infinity_value, out2); } return (CGAL::max)(sq_dista, sq_distb); @@ -2519,10 +2538,9 @@ double bounded_error_Hausdorff_distance(const TriangleMesh1& tm1, CGAL::Emptyset_iterator()); CGAL_precondition(error_bound >= 0.); - const FT sq_error_bound = square(FT(error_bound)); const FT sq_hdist = internal::bounded_error_squared_one_sided_Hausdorff_distance_impl( - tm1, tm2, sq_error_bound, FT(-1) /*distance threshold*/, match_faces, vpm1, vpm2, np1, np2, out); + tm1, tm2, error_bound, FT(-1) /*distance threshold*/, match_faces, vpm1, vpm2, np1, np2, out); return to_double(approximate_sqrt(sq_hdist)); } @@ -2573,10 +2591,9 @@ double bounded_error_symmetric_Hausdorff_distance(const TriangleMesh1& tm1, CGAL::Emptyset_iterator()); CGAL_precondition(error_bound >= 0.); - const FT sq_error_bound = square(FT(error_bound)); const FT sq_hdist = internal::bounded_error_squared_symmetric_Hausdorff_distance_impl( - tm1, tm2, sq_error_bound, FT(-1) /*distance_threshold*/, match_faces, vpm1, vpm2, np1, np2, out1, out2); + tm1, tm2, error_bound, FT(-1) /*distance_threshold*/, match_faces, vpm1, vpm2, np1, np2, out1, out2); return to_double(approximate_sqrt(sq_hdist)); } @@ -2637,7 +2654,6 @@ bool is_Hausdorff_distance_larger(const TriangleMesh1& tm1, const bool use_one_sided = choose_parameter(get_parameter(np1, internal_np::use_one_sided_hausdorff), true); CGAL_precondition(error_bound >= 0.); - const FT sq_error_bound = square(FT(error_bound)); CGAL_precondition(distance_bound > 0.); const FT sq_distance_bound = square(FT(distance_bound)); @@ -2647,12 +2663,12 @@ bool is_Hausdorff_distance_larger(const TriangleMesh1& tm1, if(use_one_sided) { sq_hdist = internal::bounded_error_squared_one_sided_Hausdorff_distance_impl( - tm1, tm2, sq_error_bound, sq_distance_bound, match_faces, vpm1, vpm2, np1, np2, stub); + tm1, tm2, error_bound, sq_distance_bound, match_faces, vpm1, vpm2, np1, np2, stub); } else { sq_hdist = internal::bounded_error_squared_symmetric_Hausdorff_distance_impl( - tm1, tm2, sq_error_bound, sq_distance_bound, match_faces, vpm1, vpm2, np1, np2, stub, stub); + tm1, tm2, error_bound, sq_distance_bound, match_faces, vpm1, vpm2, np1, np2, stub, stub); } #ifdef CGAL_HAUSDORFF_DEBUG @@ -2690,10 +2706,9 @@ double bounded_error_Hausdorff_distance_naive(const TriangleMesh1& tm1, get_const_property_map(vertex_point, tm2)); CGAL_precondition(error_bound >= 0.); - const FT sq_error_bound = square(FT(error_bound)); const FT sq_hdist = internal::bounded_error_squared_Hausdorff_distance_naive_impl( - tm1, tm2, sq_error_bound, vpm1, vpm2); + tm1, tm2, error_bound, vpm1, vpm2); return to_double(approximate_sqrt(sq_hdist)); } 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 index d60f4523d0a..ccf77246ee6 100644 --- 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 @@ -405,7 +405,6 @@ private: const TM1_face_to_triangle_map m_face_to_triangle_map; // Internal bounds and values. - const FT m_sq_error_bound; const FT m_sq_initial_bound; const FT m_sq_distance_bound; const FT m_infinity_value; @@ -421,7 +420,6 @@ public: const TriangleMesh2& tm2, const VPM1 vpm1, const VPM2 vpm2, - const FT sq_error_bound, const FT infinity_value, const FT sq_initial_bound, const FT sq_distance_bound) @@ -429,20 +427,17 @@ public: m_vpm1(vpm1), m_vpm2(vpm2), m_tm2_tree(tree), m_face_to_triangle_map(&m_tm1, m_vpm1), - m_sq_error_bound(sq_error_bound), m_sq_initial_bound(sq_initial_bound), m_sq_distance_bound(sq_distance_bound), m_infinity_value(infinity_value), m_global_bounds(m_infinity_value), m_early_exit(false) { - CGAL_precondition(m_sq_error_bound >= FT(0)); CGAL_precondition(m_infinity_value >= FT(0)); - CGAL_precondition(m_sq_initial_bound >= m_sq_error_bound); // Bounds grow with every face of TM1 (Equation (6)). // If we initialize to zero here, then we are very slow even for big input error bounds! - // Instead, we can use m_sq_error_bound as our initial guess to filter out all pairs + // Instead, we can use the 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. m_global_bounds.lower = m_sq_initial_bound; m_global_bounds.upper = m_sq_initial_bound; From c38758db66890424f4a3d871ae63a0c1eb31b19a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 7 Apr 2022 11:40:56 +0200 Subject: [PATCH 09/14] Replace custom stop criterion with paper's --- .../CGAL/Polygon_mesh_processing/distance.h | 40 ++++++++----------- 1 file changed, 16 insertions(+), 24 deletions(-) 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 cc61ab329cb..ec41be9e0b2 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1619,31 +1619,23 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, const Point_3& v1 = triangle_for_subdivision.vertex(1); const Point_3& v2 = triangle_for_subdivision.vertex(2); - // Stopping condition: all edge lengths of the triangle are smaller than the given error bound. - if(CGAL::squared_distance(v0, v1) < sq_error_bound && - CGAL::squared_distance(v0, v2) < sq_error_bound && - CGAL::squared_distance(v1, v2) < sq_error_bound) - { -#ifdef CGAL_HAUSDORFF_DEBUG_PP - std::cout << "Third stopping condition, small triangle" << std::endl; - std::cout << CGAL::squared_distance(v0, v1) << " vs " << sq_error_bound << std::endl; - std::cout << CGAL::squared_distance(v0, v2) << " vs " << sq_error_bound << std::endl; - std::cout << CGAL::squared_distance(v1, v2) << " vs " << sq_error_bound << std::endl; -#endif + // 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; - // By definition, lower_bound(t1, TM2) is smaller than h(t1, TM2). - // Let `v` is the vertex of t1 and `p_l` the point on TM2 realizing the lower bound, i.e., - // d(v, p_l) = max_{v in t1} min_{t2 in TM2) d(v, t2). - // Let `p` be the pair of points resp. on t1 and TM2 realizing the Hausdorff distance, i.e., - // h(t1, TM2) = d(p.first, p.second), - // Since we are here: - // d(p.first, v) < error_bound (1) - // From the lower bound - // d(v, p_l) <= d(p.first, p.second) - // Because d(p.first, p.second) is the min distance from p.first to TM2, - // d(p.first, p.second) <= d(p.first, p_l) - // And by triangular inequality - // d(p.first, p.second) <= d(v, p_l) + d(p.first, v) <= d(v, p_l) + error_bound = lower_bound + epsilon continue; } From 699c0aae9cb7b119b528c60fdce9ba0b734fae5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 7 Apr 2022 11:41:15 +0200 Subject: [PATCH 10/14] Do not pollute the queue with meaningless triangles --- .../include/CGAL/Polygon_mesh_processing/distance.h | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) 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 ec41be9e0b2..a0c02c296e1 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1700,6 +1700,10 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, CGAL_assertion(sub_triangle_bounds.tm2_lface != boost::graph_traits::null_face()); CGAL_assertion(sub_triangle_bounds.tm2_uface != boost::graph_traits::null_face()); + // Add the subtriangle to the candidate list if it is meaningful + if(sub_triangle_bounds.upper < global_bounds.lower) + continue; + // The global lower bound is the max of the per-face lower bounds if(sub_triangle_bounds.lower > global_bounds.lower) { @@ -1716,7 +1720,9 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, candidate_triangles.emplace(sub_triangles[i], sub_triangle_bounds, triangle_and_bounds.tm1_face); } - CGAL_assertion(!candidate_triangles.empty()); + // In case all subdividing triangles of the last queue entry are useless + if(candidate_triangles.empty()) + break; // Update global upper Hausdorff bound after subdivision. const Candidate& top_candidate = candidate_triangles.top(); From ac8755df07147d274a1a46f6d62ee000148da7a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 7 Apr 2022 11:41:38 +0200 Subject: [PATCH 11/14] Remove assertion that can sometimes fail due to numerical errors --- .../include/CGAL/Polygon_mesh_processing/distance.h | 5 ----- 1 file changed, 5 deletions(-) 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 a0c02c296e1..3a51a926118 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1740,11 +1740,6 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, std::cout << "Top of the queue is lower than the lowest!" << std::endl; #endif - // This can happen if the actual Hausdorff distance is smaller than sq_initial_bound, - // which can happen either if the error bound is really large, or if the initial bound - // is non-zero (symmetric Hausdorff distance). In this case, we return 'sq_initial_bound' - CGAL_assertion(global_bounds.lower == sq_initial_bound); - global_bounds.upper = global_bounds.lower; // not really needed since lower is returned but doesn't hurt global_bounds.upair.first = global_bounds.lpair.first; global_bounds.upair.second = global_bounds.lpair.second; From 898382be0f20ce5caf69e5b7113d794af734a19c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 Apr 2022 13:27:52 +0200 Subject: [PATCH 12/14] Fix unused warning --- .../include/CGAL/Polygon_mesh_processing/distance.h | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) 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 3a51a926118..88dc0c830e9 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1445,11 +1445,9 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, using Point_3 = typename Kernel::Point_3; using Triangle_3 = typename Kernel::Triangle_3; - const FT sq_error_bound = square(FT(error_bound)); - #ifdef CGAL_HAUSDORFF_DEBUG std::cout << " -- Bounded Hausdorff --" << std::endl; - std::cout << "error bound: " << error_bound << " (square: " << sq_error_bound << ")" << std::endl; + std::cout << "error bound: " << error_bound << std::endl; std::cout << "initial bound: " << sq_initial_bound << " (" << approximate_sqrt(sq_initial_bound) << ")" << std::endl; std::cout << "distance bound: " << sq_distance_bound << " (" << approximate_sqrt(sq_distance_bound) << ")" << std::endl; std::cout << "inf val: " << infinity_value << " (" << approximate_sqrt(infinity_value) << ")" << std::endl; @@ -1463,9 +1461,8 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, using Candidate = Candidate_triangle; - CGAL_precondition(sq_initial_bound >= sq_error_bound); + CGAL_precondition(sq_initial_bound >= square(FT(error_bound))); CGAL_precondition(sq_distance_bound != FT(0)); // value is -1 if unused - CGAL_precondition(sq_error_bound >= FT(0)); CGAL_precondition(tm1_tree.size() > 0); CGAL_precondition(tm2_tree.size() > 0); From 720c52f00acab4a8eee294f2bf4d23eff1bfffbc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 Apr 2022 13:28:04 +0200 Subject: [PATCH 13/14] Do not miss global_bounds.upper updates by not pushing subdivide faces into PQ --- .../include/CGAL/Polygon_mesh_processing/distance.h | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) 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 88dc0c830e9..1c23b587093 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1697,10 +1697,6 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, CGAL_assertion(sub_triangle_bounds.tm2_lface != boost::graph_traits::null_face()); CGAL_assertion(sub_triangle_bounds.tm2_uface != boost::graph_traits::null_face()); - // Add the subtriangle to the candidate list if it is meaningful - if(sub_triangle_bounds.upper < global_bounds.lower) - continue; - // The global lower bound is the max of the per-face lower bounds if(sub_triangle_bounds.lower > global_bounds.lower) { @@ -1717,16 +1713,13 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, candidate_triangles.emplace(sub_triangles[i], sub_triangle_bounds, triangle_and_bounds.tm1_face); } - // In case all subdividing triangles of the last queue entry are useless - if(candidate_triangles.empty()) - break; - // Update global upper Hausdorff bound after subdivision. const Candidate& top_candidate = candidate_triangles.top(); const FT current_upmost = top_candidate.bounds.upper; #ifdef CGAL_HAUSDORFF_DEBUG_PP - std::cout << "current upper bound = " << current_upmost << std::endl; std::cout << "global_bounds.lower = " << global_bounds.lower << std::endl; + std::cout << "global_bounds.upper = " << global_bounds.upper << std::endl; + std::cout << "current upper bound = " << current_upmost << std::endl; #endif CGAL_assertion(is_positive(current_upmost)); From cb9cf1d0e3410d680e79cebad13bcbf52c98193a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 8 Apr 2022 13:29:40 +0200 Subject: [PATCH 14/14] Add some debug info --- .../include/CGAL/Polygon_mesh_processing/distance.h | 4 ++++ 1 file changed, 4 insertions(+) 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 1c23b587093..d851e07cc93 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1627,6 +1627,10 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1, if((closest_triangle_v0.second == closest_triangle_v1.second) && (closest_triangle_v1.second == closest_triangle_v2.second)) { +#ifdef CGAL_HAUSDORFF_DEBUG_PP + std::cout << "Projects onto the same TM2 face" << std::endl; +#endif + // 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.