mirror of https://github.com/CGAL/cgal
Improve Hausdorff distance debug code + more assertions
This commit is contained in:
Mael Rouxel-Labbé
parent
ae28ae8a31
commit
fe21d98345
|
|
@ -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<Point_3, Kernel, TriangleMesh1, TriangleMesh2, VPM1, VPM2>;
|
||||
using TM2_hd_traits = Hausdorff_primitive_traits_tm2<Triangle_3, Kernel, TriangleMesh1, TriangleMesh2, VPM2>;
|
||||
|
||||
|
|
@ -1453,6 +1461,7 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1,
|
|||
|
||||
using Candidate = Candidate_triangle<Kernel, Face_handle_1, Face_handle_2>;
|
||||
|
||||
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<TriangleMesh2>::null_face());
|
||||
CGAL_assertion(sub_triangle_bounds.tm2_uface != boost::graph_traits<TriangleMesh2>::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<Kernel>(
|
||||
tm1, tm2, vpm1, vpm2, tm1_tree, tm2_tree,
|
||||
|
|
|
|||
|
|
@ -25,12 +25,6 @@
|
|||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
#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<TriangleMesh2>::null_face());
|
||||
CGAL_assertion(local_bounds.tm2_uface != boost::graph_traits<TriangleMesh2>::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));
|
||||
|
|
|
|||
Loading…
Reference in New Issue