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