From 82e8e4bb46630b9581c39c7b98885df4411027ce Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Tue, 21 May 2019 11:56:22 +0200 Subject: [PATCH 01/48] Initial version of bounded_error_hausdorff method. --- .../Polygon_mesh_processing/CMakeLists.txt | 9 +- ...usdorff_bounded_error_distance_example.cpp | 33 +++ .../CGAL/Polygon_mesh_processing/distance.h | 201 +++++++++++++++++- 3 files changed, 237 insertions(+), 6 deletions(-) create mode 100644 Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 2f853f6cbe6..4053c3a4ced 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -20,7 +20,7 @@ find_package( CGAL QUIET COMPONENTS ) if ( NOT CGAL_FOUND ) message(STATUS "This project requires the CGAL library, and will not be compiled.") - return() + return() endif() @@ -31,7 +31,7 @@ if ( NOT Boost_FOUND ) message(STATUS "This project requires the Boost library, and will not be compiled.") - return() + return() endif() @@ -55,8 +55,10 @@ find_package(Eigen3 3.2.0) #(requires 3.2.0 or greater) find_package( TBB ) create_single_source_cgal_program( "hausdorff_distance_remeshing_example.cpp") +create_single_source_cgal_program( "hausdorff_bounded_error_distance_example.cpp") if( TBB_FOUND ) CGAL_target_use_TBB(hausdorff_distance_remeshing_example) + CGAL_target_use_TBB(hausdorff_bounded_error_distance_example) else() message( STATUS "NOTICE: Intel TBB was not found. hausdorff_distance_example will use sequential code." ) endif() @@ -118,6 +120,3 @@ target_link_libraries( stitch_borders_example_OM PRIVATE ${OPENMESH_LIBRARIES} ) create_single_source_cgal_program( "triangulate_faces_example_OM.cpp") target_link_libraries( triangulate_faces_example_OM PRIVATE ${OPENMESH_LIBRARIES} ) endif(OpenMesh_FOUND) - - - diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp new file mode 100644 index 00000000000..dd32a2b5e4e --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -0,0 +1,33 @@ +#include +#include +#include +#include + +#if defined(CGAL_LINKED_WITH_TBB) +#define TAG CGAL::Parallel_tag +#else +#define TAG CGAL::Sequential_tag +#endif + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef K::Point_3 Point; +typedef CGAL::Surface_mesh Mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main() +{ + Mesh tm1, tm2; + CGAL::make_tetrahedron(Point(.0,.0,.0), + Point(2,.0,.0), + Point(1,1,1), + Point(1,.0,2), + tm1); + tm2=tm1; + CGAL::Polygon_mesh_processing::isotropic_remeshing(tm2.faces(),.05, tm2); + + std::cout << "Approximated Hausdorff distance: " + << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance + (tm1, tm2, 0.001) + << std::endl; +} 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 944fe523fb3..1ef31ac4f8c 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -809,8 +810,206 @@ double approximate_symmetric_Hausdorff_distance(const TriangleMesh& tm1, tm1, tm2, parameters::all_default(), parameters::all_default()); } +//////////////////////////////////////////////////////////////////////// + +namespace internal { +/* +#if defined(CGAL_LINKED_WITH_TBB) +template +struct Distance_computation{ + const AABB_tree& tree; + const std::vector& sample_points; + Point_3 initial_hint; + tbb::atomic* distance; + + Distance_computation( + const AABB_tree& tree, + const Point_3& p, + const std::vector& sample_points, + tbb::atomic* d) + : tree(tree) + , sample_points(sample_points) + , initial_hint(p) + , distance(d) + {} + + void + operator()(const tbb::blocked_range& range) const + { + Point_3 hint = initial_hint; + double hdist = 0; + for( std::size_t i = range.begin(); i != range.end(); ++i) + { + hint = tree.closest_point(sample_points[i], hint); + typename Kernel_traits::Kernel::Compute_squared_distance_3 squared_distance; + double d = to_double(CGAL::approximate_sqrt( squared_distance(hint,sample_points[i]) )); + if (d>hdist) hdist=d; + } + + // update max value stored in distance + double current_value = *distance; + while( current_value < hdist ) + { + current_value = distance->compare_and_swap(hdist, current_value); + } + } +}; +#endif +*/ +template +double bounded_error_Hausdorff_impl( + const TriangleMesh& tm1, + const TriangleMesh& tm2, + const typename Kernel::FT& error_bound, + VPM1 vpm1, + VPM2 vpm2) +{ + CGAL_assertion_code( bool is_triangle = is_triangle_mesh(tm1) && is_triangle_mesh(tm2) ); + CGAL_assertion_msg (is_triangle, + "One of the meshes is not triangulated. Distance computing impossible."); + + typedef typename Kernel::Point_3 Point_3; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; + + std::vector tm1_vertices; + tm1_vertices.reserve(num_vertices(tm1)); + tm1_vertices.insert(tm1_vertices.end(),vertices(tm1).begin(),vertices(tm1).end()); + + // Sort vertices along a Hilbert curve + spatial_sort( tm1_vertices.begin(), + tm1_vertices.end(), + Search_traits_3(vpm1) ); + + typedef AABB_face_graph_triangle_primitive TM2_primitive; + typedef AABB_tree< AABB_traits > TM2_tree; + + // Build an AABB tree on the second mesh + TM2_tree tm2_tree( faces(tm2).begin(), faces(tm2).end(), tm2, vpm2 ); + tm2_tree.build(); + tm2_tree.accelerate_distance_queries(); + Point_3 hint = get(vpm2, *vertices(tm2).begin()); + +#if !defined(CGAL_LINKED_WITH_TBB) + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + // if (boost::is_convertible::value) + // { + // tbb::atomic distance; + // distance=0; + // Distance_computation f(tm2_tree + // , hint, sample_points, &distance); + // tbb::parallel_for(tbb::blocked_range(0, sample_points.size()), f); + // return distance; + // } + // else +#endif + { + double hdist = 0; + for(const typename Kernel::Point_3& pt : sample_points) + { + hint = tree.closest_point(pt, hint); + typename Kernel::Compute_squared_distance_3 squared_distance; + typename Kernel::FT dist = squared_distance(hint,pt); + double d = to_double(CGAL::approximate_sqrt(dist)); + if(d>hdist) + hdist=d; + } + return hdist; + } + return 0.; } -} // end of namespace CGAL::Polygon_mesh_processing + +} //end of namespace internal + +/** + * \ingroup PMP_distance_grp + * computes the approximate Hausdorff distance from `tm1` to `tm2` by returning + * the distance of the farthest point from `tm2` amongst a sampling of `tm1` + * generated with the function `sample_triangle_mesh()` with + * `tm1` and `np1` as parameter. + * + * A parallel version is provided and requires the executable to be + * linked against the Intel TBB library. + * To control the number of threads used, the user may use the `tbb::task_scheduler_init` class. + * See the TBB documentation + * for more details. + * + * @tparam Concurrency_tag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag` + * and `Parallel_tag`. + * @tparam TriangleMesh a model of the concept `FaceListGraph` + * @tparam NamedParameters1 a sequence of \ref pmp_namedparameters "Named Parameters" for `tm1` + * @tparam NamedParameters2 a sequence of \ref pmp_namedparameters "Named Parameters" for `tm2` + * + * @param tm1 the triangle mesh that will be sampled + * @param tm2 the triangle mesh to compute the distance to + * @param np1 optional sequence of \ref pmp_namedparameters "Named Parameters" for `tm1` passed to `sample_triangle_mesh()`. + * + * @param np2 optional sequence of \ref pmp_namedparameters "Named Parameters" for `tm2` among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm2` + * If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` must be available in `TriangleMesh` + * and in all places where `vertex_point_map` is used. + * \cgalParamEnd + * \cgalNamedParamsEnd + * The function `CGAL::parameters::all_default()` can be used to indicate to use the default values for + * `np1` and specify custom values for `np2` + */ +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1, + class NamedParameters2> +double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double error_bound, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + typedef typename GetGeomTraits::type Geom_traits; + + typedef typename GetVertexPointMap::const_type Vpm1; + typedef typename GetVertexPointMap::const_type Vpm2; + + using boost::choose_param; + using boost::get_param; + + Vpm1 vpm1 = choose_param(get_param(np1, internal_np::vertex_point), + get_const_property_map(vertex_point, tm1)); + Vpm2 vpm2 = choose_param(get_param(np2, internal_np::vertex_point), + get_const_property_map(vertex_point, tm2)); + + return internal::bounded_error_Hausdorff_impl(tm1, tm2, error_bound, vpm1, vpm2); +} + +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1> +double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double error_bound, + const NamedParameters1& np1) +{ + return bounded_error_Hausdorff_distance(tm1, tm2, error_bound, np1, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh> +double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double error_bound) +{ + return bounded_error_Hausdorff_distance(tm1, tm2, error_bound, parameters::all_default() ); +} + +} } // end of namespace CGAL::Polygon_mesh_processing #endif //CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H From 04a05c0db1ca947481a6964802cab22bf8ec490f Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Tue, 21 May 2019 14:00:43 +0200 Subject: [PATCH 02/48] Implemented distance computation from mesh_1 vertices to mesh_2 triangles. --- .../CGAL/Polygon_mesh_processing/distance.h | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 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 1ef31ac4f8c..fe37c15ff3c 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -874,7 +874,10 @@ double bounded_error_Hausdorff_impl( typedef typename Kernel::Point_3 Point_3; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; + typedef CGAL::dynamic_vertex_property_t> Vertex_property_tag; + typedef typename boost::property_map::const_type Vertex_closest_triangle_map; std::vector tm1_vertices; tm1_vertices.reserve(num_vertices(tm1)); @@ -892,7 +895,9 @@ double bounded_error_Hausdorff_impl( TM2_tree tm2_tree( faces(tm2).begin(), faces(tm2).end(), tm2, vpm2 ); tm2_tree.build(); tm2_tree.accelerate_distance_queries(); - Point_3 hint = get(vpm2, *vertices(tm2).begin()); + std::pair hint = tm2_tree.any_reference_point_and_id(); + + Vertex_closest_triangle_map vctm = get(Vertex_property_tag(), tm1); #if !defined(CGAL_LINKED_WITH_TBB) CGAL_static_assertion_msg (!(boost::is_convertible::value), @@ -910,17 +915,17 @@ double bounded_error_Hausdorff_impl( // else #endif { - double hdist = 0; - for(const typename Kernel::Point_3& pt : sample_points) + for(vertex_descriptor vd : tm1_vertices) { - hint = tree.closest_point(pt, hint); + typename boost::property_traits::reference pt = get(vpm1, vd); + hint = tm2_tree.closest_point_and_primitive(pt, hint); typename Kernel::Compute_squared_distance_3 squared_distance; - typename Kernel::FT dist = squared_distance(hint,pt); + typename Kernel::FT dist = squared_distance(hint.first, pt); double d = to_double(CGAL::approximate_sqrt(dist)); - if(d>hdist) - hdist=d; + + + put(vctm, vd, std::make_pair(d, hint.second)); } - return hdist; } return 0.; } From 951e4a7f44cf635b630b39cf70fa1e6f1c111066 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Tue, 21 May 2019 17:11:18 +0200 Subject: [PATCH 03/48] Implement first rough approximation. --- .../CGAL/Polygon_mesh_processing/distance.h | 78 ++++++++++++++++--- 1 file changed, 67 insertions(+), 11 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 fe37c15ff3c..ef9998b13ea 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -49,6 +49,8 @@ #include +#include + namespace CGAL{ namespace Polygon_mesh_processing { namespace internal{ @@ -864,7 +866,7 @@ template TM2_primitive; + typedef AABB_tree< AABB_traits > TM2_tree; typedef typename Kernel::Point_3 Point_3; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typedef std::pair Hausdorff_bounds; typedef CGAL::Spatial_sort_traits_adapter_3 Search_traits_3; typedef CGAL::dynamic_vertex_property_t> Vertex_property_tag; + typedef CGAL::dynamic_face_property_t Face_property_tag; + typedef typename boost::property_map::const_type Vertex_closest_triangle_map; + typedef typename boost::property_map::const_type Triangle_hausdorff_bounds; + + typename Kernel::Compute_squared_distance_3 squared_distance; + typename Kernel::Construct_projected_point_3 project_point; + typename Kernel::FT dist; std::vector tm1_vertices; tm1_vertices.reserve(num_vertices(tm1)); @@ -888,9 +903,6 @@ double bounded_error_Hausdorff_impl( tm1_vertices.end(), Search_traits_3(vpm1) ); - typedef AABB_face_graph_triangle_primitive TM2_primitive; - typedef AABB_tree< AABB_traits > TM2_tree; - // Build an AABB tree on the second mesh TM2_tree tm2_tree( faces(tm2).begin(), faces(tm2).end(), tm2, vpm2 ); tm2_tree.build(); @@ -898,11 +910,13 @@ double bounded_error_Hausdorff_impl( std::pair hint = tm2_tree.any_reference_point_and_id(); Vertex_closest_triangle_map vctm = get(Vertex_property_tag(), tm1); + Triangle_hausdorff_bounds thb = get(Face_property_tag(), tm1); #if !defined(CGAL_LINKED_WITH_TBB) CGAL_static_assertion_msg (!(boost::is_convertible::value), "Parallel_tag is enabled but TBB is unavailable."); #else + // TODO implement parallelized version of the below here. // if (boost::is_convertible::value) // { // tbb::atomic distance; @@ -915,19 +929,61 @@ double bounded_error_Hausdorff_impl( // else #endif { + // For each vertex in the first mesh, find the closest triangle in the + // second mesh, store it and also store the distance to this triangle + // in a dynamic vertex property for(vertex_descriptor vd : tm1_vertices) { typename boost::property_traits::reference pt = get(vpm1, vd); hint = tm2_tree.closest_point_and_primitive(pt, hint); - typename Kernel::Compute_squared_distance_3 squared_distance; - typename Kernel::FT dist = squared_distance(hint.first, pt); - double d = to_double(CGAL::approximate_sqrt(dist)); - - + dist = squared_distance(hint.first, pt); + double d = to_double(dist); put(vctm, vd, std::make_pair(d, hint.second)); } + + Triangle_from_face_descriptor_map face_to_triangle_map(&tm2, vpm2); + double h_lower = 0.; + double h_upper = std::numeric_limits::infinity(); + + for(face_descriptor fd : faces(tm1)) + { + double h_triangle_lower = 0.; + double h_triangle_upper = std::numeric_limits::infinity(); + halfedge_descriptor hd = halfedge(fd, tm1); + + std::array face_vertices = {source(hd,tm1), target(hd,tm1), target(next(hd, tm1),tm1)}; + + std::array,3> vertex_properties = { + get(vctm, face_vertices[0]), + get(vctm, face_vertices[1]), + get(vctm, face_vertices[2])}; + + std::array triangles_in_B = { + get(face_to_triangle_map, vertex_properties[0].second), + get(face_to_triangle_map, vertex_properties[1].second), + get(face_to_triangle_map, vertex_properties[2].second)}; + + for(int i=0; i<3; ++i) + { + // Iterate over the vertices by i + h_triangle_lower = (std::max)(h_triangle_lower, vertex_properties[i].first); + // Iterate over the triangles by i + double face_distance_1 = vertex_properties[i].second==vertex_properties[(i+1)%3].second + ? vertex_properties[(i+1)%3].first + : squared_distance(project_point(triangles_in_B[i], get(vpm1, face_vertices[(i+1)%3])), get(vpm1, face_vertices[(i+1)%3])); + double face_distance_2 = vertex_properties[i].second==vertex_properties[(i+2)%3].second + ? vertex_properties[(i+2)%3].first + : squared_distance(project_point(triangles_in_B[i], get(vpm1, face_vertices[(i+2)%3])), get(vpm1, face_vertices[(i+2)%3])); + + h_triangle_upper = (std::min)((std::max)((std::max)(face_distance_1, face_distance_2), vertex_properties[i].first), h_triangle_upper); + } + put(thb, fd, Hausdorff_bounds(h_triangle_lower, h_triangle_upper)); + h_lower = (std::max)(h_lower, h_triangle_lower); + h_upper = (std::max)(h_upper, h_triangle_upper); + } + + return (h_lower+h_upper)/2.; } - return 0.; } } //end of namespace internal @@ -991,7 +1047,7 @@ double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, Vpm2 vpm2 = choose_param(get_param(np2, internal_np::vertex_point), get_const_property_map(vertex_point, tm2)); - return internal::bounded_error_Hausdorff_impl(tm1, tm2, error_bound, vpm1, vpm2); + return internal::bounded_error_Hausdorff_impl(tm1, tm2, error_bound*error_bound, vpm1, vpm2); } template< class Concurrency_tag, From 9df0d3fc8029142890f2d17a7888dbd32944a187 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Tue, 21 May 2019 18:30:15 +0200 Subject: [PATCH 04/48] Fix initialization of global upper bound, include TODOs for further steps. --- .../CGAL/Polygon_mesh_processing/distance.h | 26 +++++++++++++++++-- 1 file changed, 24 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 ef9998b13ea..36cce31a096 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -943,8 +943,11 @@ double bounded_error_Hausdorff_impl( Triangle_from_face_descriptor_map face_to_triangle_map(&tm2, vpm2); double h_lower = 0.; - double h_upper = std::numeric_limits::infinity(); + double h_upper = 0.; + // For each triangle in the first mesh, initialize its local upper and + // lower bound and store these in a dynamic face property for furture + // reference for(face_descriptor fd : faces(tm1)) { double h_triangle_lower = 0.; @@ -982,7 +985,26 @@ double bounded_error_Hausdorff_impl( h_upper = (std::max)(h_upper, h_triangle_upper); } - return (h_lower+h_upper)/2.; + // Initialize an array of candidate triangles in A to be procesed in the + // following + std::vector candidate_triangles; + + for(face_descriptor fd : faces(tm1)) + { + Hausdorff_bounds triangle_bounds = get(thb, fd); + if (triangle_bounds.second > h_lower) { + + // TODO culling on B + + candidate_triangles.push_back(fd); + } + } + + // TODO Iterate over candidate_triangles and kill those which cannot contribute anymore + + // TODO Send the remaining triangles to the Subdivision + + return (CGAL::approximate_sqrt(h_lower)+CGAL::approximate_sqrt(h_upper))/2.; } } From c28f0076fa758910398dffdc952e932354d0d2c4 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Thu, 6 Jun 2019 09:51:16 +0200 Subject: [PATCH 05/48] Add detailed comments on the code written so far and merge face processing loop into previous loop. --- .../CGAL/Polygon_mesh_processing/distance.h | 66 +++++++++++++------ 1 file changed, 46 insertions(+), 20 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 36cce31a096..22a6be53523 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -894,6 +894,7 @@ double bounded_error_Hausdorff_impl( typename Kernel::Construct_projected_point_3 project_point; typename Kernel::FT dist; + // Store all vertices of tm1 in a vector std::vector tm1_vertices; tm1_vertices.reserve(num_vertices(tm1)); tm1_vertices.insert(tm1_vertices.end(),vertices(tm1).begin(),vertices(tm1).end()); @@ -903,13 +904,16 @@ double bounded_error_Hausdorff_impl( tm1_vertices.end(), Search_traits_3(vpm1) ); - // Build an AABB tree on the second mesh + // Build an AABB tree on tm2 TM2_tree tm2_tree( faces(tm2).begin(), faces(tm2).end(), tm2, vpm2 ); tm2_tree.build(); tm2_tree.accelerate_distance_queries(); std::pair hint = tm2_tree.any_reference_point_and_id(); + // For each vertex in tm1, store the distance to the closest triangle of tm2 Vertex_closest_triangle_map vctm = get(Vertex_property_tag(), tm1); + // For each triangle in tm1, sotre its respective local lower and upper bound + // on the Hausdorff measure Triangle_hausdorff_bounds thb = get(Face_property_tag(), tm1); #if !defined(CGAL_LINKED_WITH_TBB) @@ -929,38 +933,55 @@ double bounded_error_Hausdorff_impl( // else #endif { - // For each vertex in the first mesh, find the closest triangle in the - // second mesh, store it and also store the distance to this triangle - // in a dynamic vertex property + /* + / For each vertex in the first mesh, find the closest triangle in the + / second mesh, store it and also store the distance to this triangle + / in a dynamic vertex property + */ for(vertex_descriptor vd : tm1_vertices) { + // Get the point represented by the vertex typename boost::property_traits::reference pt = get(vpm1, vd); + // Use the AABB tree to find the closest point and face in tm2 hint = tm2_tree.closest_point_and_primitive(pt, hint); + // Compute the distance of the point to the closest point in tm2 dist = squared_distance(hint.first, pt); double d = to_double(dist); + // Store the distance and the closest triangle in the corresponding map put(vctm, vd, std::make_pair(d, hint.second)); } + // Maps the faces of tm2 to actual triangles Triangle_from_face_descriptor_map face_to_triangle_map(&tm2, vpm2); + // Initialize global bounds on the Hausdorff measure double h_lower = 0.; double h_upper = 0.; + // Initialize an array of candidate triangles in A to be procesed in the + // following + std::vector candidate_triangles; - // For each triangle in the first mesh, initialize its local upper and - // lower bound and store these in a dynamic face property for furture - // reference + /* + / For each triangle in the first mesh, initialize its local upper and + / lower bound and store these in a dynamic face property for furture + / reference + */ for(face_descriptor fd : faces(tm1)) { + // Initialize the local bounds for the current face fd double h_triangle_lower = 0.; double h_triangle_upper = std::numeric_limits::infinity(); - halfedge_descriptor hd = halfedge(fd, tm1); + // Create a halfedge descriptor for the current face and store the vertices + halfedge_descriptor hd = halfedge(fd, tm1); std::array face_vertices = {source(hd,tm1), target(hd,tm1), target(next(hd, tm1),tm1)}; + // Get the distance and closest triangle in tm2 for each vertex of fd std::array,3> vertex_properties = { get(vctm, face_vertices[0]), get(vctm, face_vertices[1]), get(vctm, face_vertices[2])}; + // Convert the closest faces of tm2 to triangles std::array triangles_in_B = { get(face_to_triangle_map, vertex_properties[0].second), get(face_to_triangle_map, vertex_properties[1].second), @@ -968,9 +989,13 @@ double bounded_error_Hausdorff_impl( for(int i=0; i<3; ++i) { - // Iterate over the vertices by i + // Iterate over the vertices by i, the distance to the closest point in + // tm2 computed above is a lower bound for the local triangle h_triangle_lower = (std::max)(h_triangle_lower, vertex_properties[i].first); - // Iterate over the triangles by i + + // Iterate over the triangles by i, if the triangles are the same, we do + // not need to compute the distance, only compute it if the triangles + // differ double face_distance_1 = vertex_properties[i].second==vertex_properties[(i+1)%3].second ? vertex_properties[(i+1)%3].first : squared_distance(project_point(triangles_in_B[i], get(vpm1, face_vertices[(i+1)%3])), get(vpm1, face_vertices[(i+1)%3])); @@ -978,21 +1003,22 @@ double bounded_error_Hausdorff_impl( ? vertex_properties[(i+2)%3].first : squared_distance(project_point(triangles_in_B[i], get(vpm1, face_vertices[(i+2)%3])), get(vpm1, face_vertices[(i+2)%3])); - h_triangle_upper = (std::min)((std::max)((std::max)(face_distance_1, face_distance_2), vertex_properties[i].first), h_triangle_upper); + // Update the local lower bound of the triangle + h_triangle_upper = (std::min)( + (std::max)( + (std::max)(face_distance_1, face_distance_2), + vertex_properties[i].first), + h_triangle_upper); } + + // Store the computed lower and upper bound in a dynamic face property put(thb, fd, Hausdorff_bounds(h_triangle_lower, h_triangle_upper)); h_lower = (std::max)(h_lower, h_triangle_lower); h_upper = (std::max)(h_upper, h_triangle_upper); - } - // Initialize an array of candidate triangles in A to be procesed in the - // following - std::vector candidate_triangles; - - for(face_descriptor fd : faces(tm1)) - { - Hausdorff_bounds triangle_bounds = get(thb, fd); - if (triangle_bounds.second > h_lower) { + // Only process the triangle further if it can still contribute to a + // Hausdorff distance + if (h_triangle_upper > h_lower) { // TODO culling on B From a67fa389f505447c5564dd67a89f2cec7357965a Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sun, 9 Jun 2019 17:17:54 +0200 Subject: [PATCH 06/48] Add traversal trait for Hausdorff-based AABB traversal on the first triangle mesh. --- .../CGAL/Polygon_mesh_processing/distance.h | 20 ++++- ...traversal_traits_with_Hausdorff_distance.h | 80 +++++++++++++++++++ 2 files changed, 98 insertions(+), 2 deletions(-) create mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_Hausdorff_distance.h 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 22a6be53523..85d4ff65c29 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -17,14 +17,13 @@ // SPDX-License-Identifier: GPL-3.0+ // // -// Author(s) : Maxime Gimeno and Sebastien Loriot +// Author(s) : Maxime Gimeno, Sebastien Loriot, Martin Skrodzki #ifndef CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H #define CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H #include - #include #include #include @@ -41,6 +40,8 @@ #include #include +#include + #ifdef CGAL_LINKED_WITH_TBB #include #include @@ -874,8 +875,12 @@ double bounded_error_Hausdorff_impl( CGAL_assertion_msg (is_triangle, "One of the meshes is not triangulated. Distance computing impossible."); + typedef AABB_face_graph_triangle_primitive TM1_primitive; typedef AABB_face_graph_triangle_primitive TM2_primitive; + typedef AABB_tree< AABB_traits > TM1_tree; typedef AABB_tree< AABB_traits > TM2_tree; + typedef typename AABB_tree< AABB_traits >::AABB_traits Tree_traits; + typedef typename Kernel::Point_3 Point_3; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -904,6 +909,15 @@ double bounded_error_Hausdorff_impl( tm1_vertices.end(), Search_traits_3(vpm1) ); + // Build an AABB tree on tm1 + TM1_tree tm1_tree( faces(tm1).begin(), faces(tm1).end(), tm1, vpm1 ); + tm1_tree.build(); + tm1_tree.accelerate_distance_queries(); + + // Build traversal traits for tm1_tree + Hausdorff_primitive_traits traversal_traits( tm1_tree.traits() ); + tm1_tree.traversal( Point_3(0,0,0), traversal_traits ); + // Build an AABB tree on tm2 TM2_tree tm2_tree( faces(tm2).begin(), faces(tm2).end(), tm2, vpm2 ); tm2_tree.build(); @@ -1026,6 +1040,8 @@ double bounded_error_Hausdorff_impl( } } + + // TODO Iterate over candidate_triangles and kill those which cannot contribute anymore // TODO Send the remaining triangles to the Subdivision 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 new file mode 100644 index 00000000000..58bf1d85805 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_Hausdorff_distance.h @@ -0,0 +1,80 @@ +// Copyright (c) 2019 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0+ +// +// +// Author(s) : Sebastien Loriot, Martin Skrodzki +// + +#ifndef CGAL_PMP_INTERNAL_AABB_TRAVERSAL_TRAITS_WITH_HAUSDORFF_DISTANCE +#define CGAL_PMP_INTERNAL_AABB_TRAVERSAL_TRAITS_WITH_HAUSDORFF_DISTANCE + +#include + +namespace CGAL { + + /** + * @class Hausdorff_primitive_traits + */ + template + class Hausdorff_primitive_traits + { + typedef typename AABBTraits::Primitive Primitive; + typedef ::CGAL::AABB_node Node; + + public: + Hausdorff_primitive_traits(const AABBTraits& traits) + : m_traits(traits) {} + + // Explore the whole tree, i.e. always enter children if the methods + // do_intersect() below determines that it is worthwhile. + bool go_further() const { return true; } + + // Compute the explicit Hausdorff distance to the given primitive + void intersection(const Query& query, const Primitive& primitive) + { + // Have reached a single triangle + std::cout << "Reached Triangle " << primitive.id() << '\n'; + + /* TODO implement handling of a single triangle + / - Call Culling on B (First maybe don't cull, but only consider closest + / triangle), obtain local bounds for the triangle + / - Update global Hausdorff bounds according to the obtained local bounds + / - return the current best known global bounds + */ + } + + bool do_intersect(const Query& query, const Node& node) const + { + // Have reached a node, determine whether or not to enter it + + /* TODO implement processing of an AABB node + / - Determine distance of the node's bounding box (node.bbox()) to the + / closest point in B (First maybe any point) + / - If the distance is larger than the global lower bound, enter the + / node, i.e. return true. + */ + return true; + //return m_traits.do_intersect_object()(query, node.bbox()); + } + + private: + const AABBTraits& m_traits; + }; +} + +#endif //CGAL_PMP_INTERNAL_AABB_TRAVERSAL_TRAITS_WITH_HAUSDORFF_DISTANCE From 04c84efb3a5f383c87d7c8f0921a351954de9f98 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sun, 9 Jun 2019 18:08:07 +0200 Subject: [PATCH 07/48] Add second traversal trait for Hausdorff-based AABB traversal on the second triangle mesh. --- .../CGAL/Polygon_mesh_processing/distance.h | 10 +- ...traversal_traits_with_Hausdorff_distance.h | 101 +++++++++++++++++- 2 files changed, 102 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 85d4ff65c29..57a3df15b13 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -882,6 +882,7 @@ double bounded_error_Hausdorff_impl( typedef typename AABB_tree< AABB_traits >::AABB_traits Tree_traits; typedef typename Kernel::Point_3 Point_3; + typedef typename Kernel::Triangle_3 Triangle_3; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; @@ -914,16 +915,17 @@ double bounded_error_Hausdorff_impl( tm1_tree.build(); tm1_tree.accelerate_distance_queries(); - // Build traversal traits for tm1_tree - Hausdorff_primitive_traits traversal_traits( tm1_tree.traits() ); - tm1_tree.traversal( Point_3(0,0,0), traversal_traits ); - // Build an AABB tree on tm2 TM2_tree tm2_tree( faces(tm2).begin(), faces(tm2).end(), tm2, vpm2 ); tm2_tree.build(); tm2_tree.accelerate_distance_queries(); std::pair hint = tm2_tree.any_reference_point_and_id(); + // Build traversal traits for tm1_tree and tm2_tree + Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits() ); + Hausdorff_primitive_traits_tm2 traversal_traits_tm2( tm2_tree.traits() ); + tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); + // For each vertex in tm1, store the distance to the closest triangle of tm2 Vertex_closest_triangle_map vctm = get(Vertex_property_tag(), tm1); // For each triangle in tm1, sotre its respective local lower and upper bound 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 58bf1d85805..5ff34e3208d 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 @@ -27,18 +27,24 @@ namespace CGAL { + typedef std::pair Hausdorff_bounds; + /** - * @class Hausdorff_primitive_traits + * @class Hausdorff_primitive_traits_tm1 */ template - class Hausdorff_primitive_traits + class Hausdorff_primitive_traits_tm1 { typedef typename AABBTraits::Primitive Primitive; typedef ::CGAL::AABB_node Node; public: - Hausdorff_primitive_traits(const AABBTraits& traits) - : m_traits(traits) {} + Hausdorff_primitive_traits_tm1(const AABBTraits& traits) + : m_traits(traits) { + // Initialize the global bounds with 0., they will only grow. + h_lower = 0.; + h_upper = 0.; + } // Explore the whole tree, i.e. always enter children if the methods // do_intersect() below determines that it is worthwhile. @@ -48,7 +54,7 @@ namespace CGAL { void intersection(const Query& query, const Primitive& primitive) { // Have reached a single triangle - std::cout << "Reached Triangle " << primitive.id() << '\n'; + std::cout << "Reached Triangle in TM1: " << primitive.id() << '\n'; /* TODO implement handling of a single triangle / - Call Culling on B (First maybe don't cull, but only consider closest @@ -58,6 +64,8 @@ namespace CGAL { */ } + // Determine whether child nodes will still contribute to a larger + // Hausdorff distance and thus have to be entered bool do_intersect(const Query& query, const Node& node) const { // Have reached a node, determine whether or not to enter it @@ -74,6 +82,89 @@ namespace CGAL { private: const AABBTraits& m_traits; + // Global Hausdorff bounds to be taken track of during the traversal + double h_lower; + double h_upper; + }; + + + /** + * @class Hausdorff_primitive_traits_tm2 + */ + template + class Hausdorff_primitive_traits_tm2 + { + typedef typename AABBTraits::Primitive Primitive; + typedef ::CGAL::AABB_node Node; + + public: + Hausdorff_primitive_traits_tm2(const AABBTraits& traits) + : m_traits(traits) { + // Initialize the global bounds with 0., they will only grow. + h_local_upper = std::numeric_limits::infinity(); + h_local_lower_0 = std::numeric_limits::infinity(); + h_local_lower_1 = std::numeric_limits::infinity(); + h_local_lower_2 = std::numeric_limits::infinity(); + } + + // Explore the whole tree, i.e. always enter children if the methods + // do_intersect() below determines that it is worthwhile. + bool go_further() const { return true; } + + // Compute the explicit Hausdorff distance to the given primitive + void intersection(const Query& query, const Primitive& primitive) + { + // Have reached a single triangle + std::cout << "Reached Triangle in TM2:" << primitive.id() << '\n'; + + double distance = std::numeric_limits::infinity(); + /* + / TODO Determine the distance accroding to + / min_{b \in primitive} ( max_{vertex in query} ( d(vertex, b))) + */ + + // Update local upper bound + if ( distance < h_local_upper ) h_local_upper = distance; + + double vertex_distance = std::numeric_limits::infinity(); + /* + / TODO For all vertices v of the query triangle, determine the distance by + / min_{b \in primitive} ( d(v, b) ) + / Update h_local_lower_i accordingly + */ + + h_local_lower = std::max( std::max (h_local_lower_0, h_local_lower_1), h_local_lower_2 ); + } + + // Determine whether child nodes will still contribute to a smaller + // Hausdorff distance and thus have to be entered + bool do_intersect(const Query& query, const Node& node) const + { + // Have reached a node, determine whether or not to enter it + double distance = std::numeric_limits::infinity(); + + /* TODO Determine distance of the node's bounding box (node.bbox()) to + / the triangle given by the query object. + */ + + if (distance <= h_local_upper) return true; + else return false; + } + + Hausdorff_bounds get_local_bounds() + { + return Hausdorff_bounds( h_local_lower, h_local_upper ); + } + + private: + const AABBTraits& m_traits; + // Local Hausdorff bounds for the query triangle + double h_local_upper; + double h_local_lower; + // Local Hausdorff bounds for the query triangle's vertices + double h_local_lower_0; + double h_local_lower_1; + double h_local_lower_2; }; } From 14006623d98b9a756936ac19149093dafdeed28b Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sun, 9 Jun 2019 18:26:45 +0200 Subject: [PATCH 08/48] Fixed initialization of local Hausdorff bounds. --- .../internal/AABB_traversal_traits_with_Hausdorff_distance.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 5ff34e3208d..6461ea2b8e8 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 @@ -126,7 +126,7 @@ namespace CGAL { // Update local upper bound if ( distance < h_local_upper ) h_local_upper = distance; - double vertex_distance = std::numeric_limits::infinity(); + double vertex_distance = 0.; /* / TODO For all vertices v of the query triangle, determine the distance by / min_{b \in primitive} ( d(v, b) ) @@ -141,7 +141,7 @@ namespace CGAL { bool do_intersect(const Query& query, const Node& node) const { // Have reached a node, determine whether or not to enter it - double distance = std::numeric_limits::infinity(); + double distance = 0.; /* TODO Determine distance of the node's bounding box (node.bbox()) to / the triangle given by the query object. From d93aa73e7ef2069166bbb79b8c61d4883ee7ec57 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Tue, 18 Jun 2019 06:12:33 +0200 Subject: [PATCH 09/48] Pass tm2_tree instance to the traversal of tm1 --- .../CGAL/Polygon_mesh_processing/distance.h | 3 +- ...traversal_traits_with_Hausdorff_distance.h | 125 +++++++++--------- 2 files changed, 67 insertions(+), 61 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 57a3df15b13..31ff50ccff7 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -922,8 +922,7 @@ double bounded_error_Hausdorff_impl( std::pair hint = tm2_tree.any_reference_point_and_id(); // Build traversal traits for tm1_tree and tm2_tree - Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits() ); - Hausdorff_primitive_traits_tm2 traversal_traits_tm2( tm2_tree.traits() ); + Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree ); tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); // For each vertex in tm1, store the distance to the closest triangle of tm2 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 6461ea2b8e8..7684159ba4c 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 @@ -29,65 +29,6 @@ namespace CGAL { typedef std::pair Hausdorff_bounds; - /** - * @class Hausdorff_primitive_traits_tm1 - */ - template - class Hausdorff_primitive_traits_tm1 - { - typedef typename AABBTraits::Primitive Primitive; - typedef ::CGAL::AABB_node Node; - - public: - Hausdorff_primitive_traits_tm1(const AABBTraits& traits) - : m_traits(traits) { - // Initialize the global bounds with 0., they will only grow. - h_lower = 0.; - h_upper = 0.; - } - - // Explore the whole tree, i.e. always enter children if the methods - // do_intersect() below determines that it is worthwhile. - bool go_further() const { return true; } - - // Compute the explicit Hausdorff distance to the given primitive - void intersection(const Query& query, const Primitive& primitive) - { - // Have reached a single triangle - std::cout << "Reached Triangle in TM1: " << primitive.id() << '\n'; - - /* TODO implement handling of a single triangle - / - Call Culling on B (First maybe don't cull, but only consider closest - / triangle), obtain local bounds for the triangle - / - Update global Hausdorff bounds according to the obtained local bounds - / - return the current best known global bounds - */ - } - - // Determine whether child nodes will still contribute to a larger - // Hausdorff distance and thus have to be entered - bool do_intersect(const Query& query, const Node& node) const - { - // Have reached a node, determine whether or not to enter it - - /* TODO implement processing of an AABB node - / - Determine distance of the node's bounding box (node.bbox()) to the - / closest point in B (First maybe any point) - / - If the distance is larger than the global lower bound, enter the - / node, i.e. return true. - */ - return true; - //return m_traits.do_intersect_object()(query, node.bbox()); - } - - private: - const AABBTraits& m_traits; - // Global Hausdorff bounds to be taken track of during the traversal - double h_lower; - double h_upper; - }; - - /** * @class Hausdorff_primitive_traits_tm2 */ @@ -166,6 +107,72 @@ namespace CGAL { double h_local_lower_1; double h_local_lower_2; }; + + /** + * @class Hausdorff_primitive_traits_tm1 + */ + template + class Hausdorff_primitive_traits_tm1 + { + typedef AABB_face_graph_triangle_primitive TM2_primitive; + typedef typename AABB_tree< AABB_traits >::AABB_traits Tree_traits; + typedef typename AABBTraits::Primitive Primitive; + typedef ::CGAL::AABB_node Node; + typedef typename Kernel::Triangle_3 Triangle_3; + typedef AABB_tree< AABB_traits > TM2_tree; + + public: + Hausdorff_primitive_traits_tm1(const AABBTraits& traits, const TM2_tree& tree) + : m_traits(traits), tm2_tree(tree) { + // Initialize the global bounds with 0., they will only grow. + h_lower = 0.; + h_upper = 0.; + } + + // Explore the whole tree, i.e. always enter children if the methods + // do_intersect() below determines that it is worthwhile. + bool go_further() const { return true; } + + // Compute the explicit Hausdorff distance to the given primitive + void intersection(const Query& query, const Primitive& primitive) + { + // Have reached a single triangle + std::cout << "Reached Triangle in TM1: " << primitive.id() << '\n'; + + + Hausdorff_primitive_traits_tm2 traversal_traits_tm2( tm2_tree.traits() ); + /* TODO implement handling of a single triangle + / - Call Culling on B (First maybe don't cull, but only consider closest + / triangle), obtain local bounds for the triangle + / - Update global Hausdorff bounds according to the obtained local bounds + / - return the current best known global bounds + */ + } + + // Determine whether child nodes will still contribute to a larger + // Hausdorff distance and thus have to be entered + bool do_intersect(const Query& query, const Node& node) const + { + // Have reached a node, determine whether or not to enter it + + /* TODO implement processing of an AABB node + / - Determine distance of the node's bounding box (node.bbox()) to the + / closest point in B (First maybe any point) + / - If the distance is larger than the global lower bound, enter the + / node, i.e. return true. + */ + return true; + //return m_traits.do_intersect_object()(query, node.bbox()); + } + + private: + const AABBTraits& m_traits; + // AABB tree for the second triangle meshes + const TM2_tree& tm2_tree; + // Global Hausdorff bounds to be taken track of during the traversal + double h_lower; + double h_upper; + }; } #endif //CGAL_PMP_INTERNAL_AABB_TRAVERSAL_TRAITS_WITH_HAUSDORFF_DISTANCE From a4112c4c8b32a17a1b05e3fe07ca60243ebb3efe Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Tue, 18 Jun 2019 07:01:06 +0200 Subject: [PATCH 10/48] Implement calling of TM2_tree traversal. --- ..._traversal_traits_with_Hausdorff_distance.h | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) 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 7684159ba4c..c77a3a227cf 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 @@ -139,14 +139,18 @@ namespace CGAL { // Have reached a single triangle std::cout << "Reached Triangle in TM1: " << primitive.id() << '\n'; + // Call Culling on B with the single triangle found. + Hausdorff_primitive_traits_tm2 traversal_traits_tm2( tm2_tree.traits() ); + tm2_tree.traversal(primitive, traversal_traits_tm2); - Hausdorff_primitive_traits_tm2 traversal_traits_tm2( tm2_tree.traits() ); - /* TODO implement handling of a single triangle - / - Call Culling on B (First maybe don't cull, but only consider closest - / triangle), obtain local bounds for the triangle - / - Update global Hausdorff bounds according to the obtained local bounds - / - return the current best known global bounds - */ + // Update global Hausdorff bounds according to the obtained local bounds + Hausdorff_bounds local_bounds = traversal_traits_tm2.get_local_bounds(); + if (local_bounds.first > h_lower) { + h_lower = local_bounds.first; + } + if (local_bounds.second > h_upper) { + h_upper = local_bounds.second; + } } // Determine whether child nodes will still contribute to a larger From c76db1e3ab663de7a005c8378b4d4970b174e1b0 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Tue, 18 Jun 2019 07:27:33 +0200 Subject: [PATCH 11/48] Finish traversal traits of tm1, i.e. finish implementation of CullingOnA. --- ...traversal_traits_with_Hausdorff_distance.h | 34 ++++++++++++++----- 1 file changed, 25 insertions(+), 9 deletions(-) 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 c77a3a227cf..6be21480eb3 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 @@ -117,7 +117,9 @@ namespace CGAL { typedef AABB_face_graph_triangle_primitive TM2_primitive; typedef typename AABB_tree< AABB_traits >::AABB_traits Tree_traits; typedef typename AABBTraits::Primitive Primitive; + typedef typename AABBTraits::Bounding_box Bounding_box; typedef ::CGAL::AABB_node Node; + typedef typename Kernel::Point_3 Point_3; typedef typename Kernel::Triangle_3 Triangle_3; typedef AABB_tree< AABB_traits > TM2_tree; @@ -157,16 +159,30 @@ namespace CGAL { // Hausdorff distance and thus have to be entered bool do_intersect(const Query& query, const Node& node) const { - // Have reached a node, determine whether or not to enter it + /* Have reached a node, determine whether or not to enter it */ - /* TODO implement processing of an AABB node - / - Determine distance of the node's bounding box (node.bbox()) to the - / closest point in B (First maybe any point) - / - If the distance is larger than the global lower bound, enter the - / node, i.e. return true. - */ - return true; - //return m_traits.do_intersect_object()(query, node.bbox()); + // Get the bounding box of the nodes + Bounding_box bbox = node.bbox(); + // Compute its center + Point_3 center = Point_3( + (bbox.min(0) + bbox.max(0)) / 2, + (bbox.min(1) + bbox.max(1)) / 2, + (bbox.min(2) + bbox.max(2)) / 2); + // Find the point from TM2 closest to the center + Point_3 closest = tm2_tree.closest_point(center); + // Compute the distance of the center to the closest point in tm2 + double dist = approximate_sqrt(squared_distance(center, closest)); + // Compute the radius of the circumsphere of the bounding boxes + double radius = approximate_sqrt(squared_distance( + Point_3(bbox.min(0),bbox.min(1),bbox.min(2)), + Point_3(bbox.max(0),bbox.max(1),bbox.max(2))) + )/2.; + // If the distance is larger than the global lower bound, enter the node, i.e. return true. + if (dist + radius > h_lower) { + return true; + } else { + return false; + } } private: From 453ec1571ceb38dbbc283d3cfcea9de22a1db652 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Wed, 19 Jun 2019 18:47:41 +0200 Subject: [PATCH 12/48] Implement vertex (TM1) to triangle (TM2) distance in traversal processing. --- .../CGAL/Polygon_mesh_processing/distance.h | 2 +- ...traversal_traits_with_Hausdorff_distance.h | 62 ++++++++++++++++--- 2 files changed, 54 insertions(+), 10 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 31ff50ccff7..4c269ce9823 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -922,7 +922,7 @@ double bounded_error_Hausdorff_impl( std::pair hint = tm2_tree.any_reference_point_and_id(); // Build traversal traits for tm1_tree and tm2_tree - Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree ); + Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2 ); tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); // For each vertex in tm1, store the distance to the closest triangle of tm2 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 6be21480eb3..86cf640af02 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 @@ -32,15 +32,18 @@ namespace CGAL { /** * @class Hausdorff_primitive_traits_tm2 */ - template + template class Hausdorff_primitive_traits_tm2 { typedef typename AABBTraits::Primitive Primitive; typedef ::CGAL::AABB_node Node; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typename Kernel::Construct_projected_point_3 project_point; public: - Hausdorff_primitive_traits_tm2(const AABBTraits& traits) - : m_traits(traits) { + Hausdorff_primitive_traits_tm2(const AABBTraits& traits, const TriangleMesh& tm1, const TriangleMesh& tm2, const VPM1& vpm1, const VPM2& vpm2) + : m_traits(traits), m_tm1(tm1), m_tm2(tm2), m_vpm1(vpm1), m_vpm2(vpm2) { // Initialize the global bounds with 0., they will only grow. h_local_upper = std::numeric_limits::infinity(); h_local_lower_0 = std::numeric_limits::infinity(); @@ -58,12 +61,39 @@ namespace CGAL { // Have reached a single triangle std::cout << "Reached Triangle in TM2:" << primitive.id() << '\n'; - double distance = std::numeric_limits::infinity(); /* - / TODO Determine the distance accroding to + / Determine the distance accroding to / min_{b \in primitive} ( max_{vertex in query} ( d(vertex, b))) + / + / Here, we only have one triangle in B, i.e. tm2. Thus, it suffices to + / compute the distance of the vertices of the query triangles to the + / primitive triangle and use the maximum of the obtained distances. */ + // The query object is a triangle from TM1, get its vertices + halfedge_descriptor hd = halfedge(query.id(), m_tm1); + vertex_descriptor v0 = source(hd,m_tm1); + vertex_descriptor v1 = target(hd,m_tm1); + vertex_descriptor v2 = target(next(hd, m_tm1), m_tm1); + + // Compute distances of the vertices to the primitive triangle in TM2 + Triangle_from_face_descriptor_map face_to_triangle_map(&m_tm2, m_vpm2); + double v0_dist = squared_distance( + project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v0)), + get(m_vpm1, v0) + ); + double v1_dist = squared_distance( + project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v1)), + get(m_vpm1, v1) + ); + double v2_dist = squared_distance( + project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v2)), + get(m_vpm1, v2) + ); + + // Get the distance as maximizers over all vertices + double distance = approximate_sqrt(std::max(std::max(v0_dist,v1_dist),v2_dist)); + // Update local upper bound if ( distance < h_local_upper ) h_local_upper = distance; @@ -99,6 +129,12 @@ namespace CGAL { private: const AABBTraits& m_traits; + // The two triangle meshes + const TriangleMesh& m_tm1; + const TriangleMesh& m_tm2; + // Their respective vertex-point Maps + const VPM1& m_vpm1; + const VPM2& m_vpm2; // Local Hausdorff bounds for the query triangle double h_local_upper; double h_local_lower; @@ -111,7 +147,7 @@ namespace CGAL { /** * @class Hausdorff_primitive_traits_tm1 */ - template + template class Hausdorff_primitive_traits_tm1 { typedef AABB_face_graph_triangle_primitive TM2_primitive; @@ -122,10 +158,12 @@ namespace CGAL { typedef typename Kernel::Point_3 Point_3; typedef typename Kernel::Triangle_3 Triangle_3; typedef AABB_tree< AABB_traits > TM2_tree; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; public: - Hausdorff_primitive_traits_tm1(const AABBTraits& traits, const TM2_tree& tree) - : m_traits(traits), tm2_tree(tree) { + Hausdorff_primitive_traits_tm1(const AABBTraits& traits, const TM2_tree& tree, const TriangleMesh& tm1, const TriangleMesh& tm2 , const VPM1& vpm1, const VPM2& vpm2 ) + : m_traits(traits), tm2_tree(tree), m_tm1(tm1), m_tm2(tm2), m_vpm1(vpm1), m_vpm2(vpm2) { // Initialize the global bounds with 0., they will only grow. h_lower = 0.; h_upper = 0.; @@ -142,7 +180,7 @@ namespace CGAL { std::cout << "Reached Triangle in TM1: " << primitive.id() << '\n'; // Call Culling on B with the single triangle found. - Hausdorff_primitive_traits_tm2 traversal_traits_tm2( tm2_tree.traits() ); + Hausdorff_primitive_traits_tm2 traversal_traits_tm2( tm2_tree.traits(), m_tm1, m_tm2, m_vpm1, m_vpm2 ); tm2_tree.traversal(primitive, traversal_traits_tm2); // Update global Hausdorff bounds according to the obtained local bounds @@ -187,6 +225,12 @@ namespace CGAL { private: const AABBTraits& m_traits; + // The two triangle meshes + const TriangleMesh& m_tm1; + const TriangleMesh& m_tm2; + // Their vertex-point-maps + const VPM1 m_vpm1; + const VPM2 m_vpm2; // AABB tree for the second triangle meshes const TM2_tree& tm2_tree; // Global Hausdorff bounds to be taken track of during the traversal From 6221ebe148e95b4625650bf9495a1a2bb7777d1c Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 21 Jun 2019 07:22:04 +0200 Subject: [PATCH 13/48] Finish implementation of Culling on B (i.e. TM2) traversal. --- ...traversal_traits_with_Hausdorff_distance.h | 68 +++++++++++++------ 1 file changed, 49 insertions(+), 19 deletions(-) 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 86cf640af02..a140d9aa034 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 @@ -37,9 +37,11 @@ namespace CGAL { { typedef typename AABBTraits::Primitive Primitive; typedef ::CGAL::AABB_node Node; + typedef typename AABBTraits::Bounding_box Bounding_box; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typename Kernel::Construct_projected_point_3 project_point; + typedef typename Kernel::Point_3 Point_3; public: Hausdorff_primitive_traits_tm2(const AABBTraits& traits, const TriangleMesh& tm1, const TriangleMesh& tm2, const VPM1& vpm1, const VPM2& vpm2) @@ -72,8 +74,8 @@ namespace CGAL { // The query object is a triangle from TM1, get its vertices halfedge_descriptor hd = halfedge(query.id(), m_tm1); - vertex_descriptor v0 = source(hd,m_tm1); - vertex_descriptor v1 = target(hd,m_tm1); + vertex_descriptor v0 = source(hd, m_tm1); + vertex_descriptor v1 = target(hd, m_tm1); vertex_descriptor v2 = target(next(hd, m_tm1), m_tm1); // Compute distances of the vertices to the primitive triangle in TM2 @@ -94,32 +96,60 @@ namespace CGAL { // Get the distance as maximizers over all vertices double distance = approximate_sqrt(std::max(std::max(v0_dist,v1_dist),v2_dist)); - // Update local upper bound + // Since we are at the level of a single triangle in TM2, distance is + // actually the correct Hausdorff distance from the query triangle in + // TM1 to the primitive triangle in TM2 if ( distance < h_local_upper ) h_local_upper = distance; - - double vertex_distance = 0.; - /* - / TODO For all vertices v of the query triangle, determine the distance by - / min_{b \in primitive} ( d(v, b) ) - / Update h_local_lower_i accordingly - */ - - h_local_lower = std::max( std::max (h_local_lower_0, h_local_lower_1), h_local_lower_2 ); + if ( distance < h_local_lower ) h_local_lower = distance; } // Determine whether child nodes will still contribute to a smaller // Hausdorff distance and thus have to be entered bool do_intersect(const Query& query, const Node& node) const { - // Have reached a node, determine whether or not to enter it - double distance = 0.; + /* Have reached a node, determine whether or not to enter it */ - /* TODO Determine distance of the node's bounding box (node.bbox()) to - / the triangle given by the query object. - */ + // Get the bounding box of the nodes + Bounding_box bbox = node.bbox(); + // Compute its center + Point_3 bb_center = Point_3( + (bbox.min(0) + bbox.max(0)) / 2, + (bbox.min(1) + bbox.max(1)) / 2, + (bbox.min(2) + bbox.max(2)) / 2); - if (distance <= h_local_upper) return true; - else return false; + // Get the center of the query triangle + // The query object is a triangle from TM1, get its vertices + halfedge_descriptor hd = halfedge(query.id(), m_tm1); + Point_3 v0 = get(m_vpm1, source(hd, m_tm1)); + Point_3 v1 = get(m_vpm1, target(hd, m_tm1)); + Point_3 v2 = get(m_vpm1, target(next(hd, m_tm1), m_tm1)); + // Compute the barycenter of the triangle + Point_3 tri_center = Point_3( 0.3*(v0.x()+v1.x()+v2.x()), 0.3*(v0.y()+v1.y()+v2.y()), 0.3*(v0.z()+v1.z()+v2.z()) ); + + // Compute the distance of the center to the closest point in tm2 + double dist = approximate_sqrt(squared_distance(bb_center, tri_center)); + + // Compute the radius of the circumsphere of the bounding boxes + double bb_radius = approximate_sqrt(squared_distance( + Point_3(bbox.min(0),bbox.min(1),bbox.min(2)), + Point_3(bbox.max(0),bbox.max(1),bbox.max(2))) + )/2.; + + // Compute the radius of the circumcircle of the triangle + double tri_radius = approximate_sqrt( std::max( + std::max( + squared_distance(tri_center, v0), + squared_distance(tri_center, v1) + ), + squared_distance(tri_center, v2) + )); + + // If triangles in the node can still lower the upper bound, enter it + if (tri_radius + bb_radius >= dist) { + return true; + } else { + return false; + } } Hausdorff_bounds get_local_bounds() From 48e159ba3f48979ab6540b79b8f8bc2487ac16da Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 21 Jun 2019 22:55:56 +0200 Subject: [PATCH 14/48] Store candidate triangles in a set to be processed later. --- .../CGAL/Polygon_mesh_processing/distance.h | 37 +++++++----- ...traversal_traits_with_Hausdorff_distance.h | 56 +++++++++++++++---- 2 files changed, 67 insertions(+), 26 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 4c269ce9823..45ab6e038dc 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -896,20 +896,14 @@ double bounded_error_Hausdorff_impl( typedef typename boost::property_map::const_type Vertex_closest_triangle_map; typedef typename boost::property_map::const_type Triangle_hausdorff_bounds; + typedef std::pair Hausdorff_bounds; + typedef std::pair Candidate_triangle; + typedef typename std::vector Candidate_set; + typename Kernel::Compute_squared_distance_3 squared_distance; typename Kernel::Construct_projected_point_3 project_point; typename Kernel::FT dist; - // Store all vertices of tm1 in a vector - std::vector tm1_vertices; - tm1_vertices.reserve(num_vertices(tm1)); - tm1_vertices.insert(tm1_vertices.end(),vertices(tm1).begin(),vertices(tm1).end()); - - // Sort vertices along a Hilbert curve - spatial_sort( tm1_vertices.begin(), - tm1_vertices.end(), - Search_traits_3(vpm1) ); - // Build an AABB tree on tm1 TM1_tree tm1_tree( faces(tm1).begin(), faces(tm1).end(), tm1, vpm1 ); tm1_tree.build(); @@ -925,11 +919,8 @@ double bounded_error_Hausdorff_impl( Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2 ); tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); - // For each vertex in tm1, store the distance to the closest triangle of tm2 - Vertex_closest_triangle_map vctm = get(Vertex_property_tag(), tm1); - // For each triangle in tm1, sotre its respective local lower and upper bound - // on the Hausdorff measure - Triangle_hausdorff_bounds thb = get(Face_property_tag(), tm1); + Candidate_set candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); + Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); #if !defined(CGAL_LINKED_WITH_TBB) CGAL_static_assertion_msg (!(boost::is_convertible::value), @@ -948,6 +939,22 @@ double bounded_error_Hausdorff_impl( // else #endif { + // Store all vertices of tm1 in a vector + std::vector tm1_vertices; + tm1_vertices.reserve(num_vertices(tm1)); + tm1_vertices.insert(tm1_vertices.end(),vertices(tm1).begin(),vertices(tm1).end()); + + // Sort vertices along a Hilbert curve + spatial_sort( tm1_vertices.begin(), + tm1_vertices.end(), + Search_traits_3(vpm1) ); + + // For each vertex in tm1, store the distance to the closest triangle of tm2 + Vertex_closest_triangle_map vctm = get(Vertex_property_tag(), tm1); + // For each triangle in tm1, sotre its respective local lower and upper bound + // on the Hausdorff measure + Triangle_hausdorff_bounds thb = get(Face_property_tag(), tm1); + /* / For each vertex in the first mesh, find the closest triangle in the / second mesh, store it and also store the distance to this triangle 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 a140d9aa034..662399c4de1 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 @@ -44,13 +44,17 @@ namespace CGAL { typedef typename Kernel::Point_3 Point_3; public: - Hausdorff_primitive_traits_tm2(const AABBTraits& traits, const TriangleMesh& tm1, const TriangleMesh& tm2, const VPM1& vpm1, const VPM2& vpm2) + Hausdorff_primitive_traits_tm2( + const AABBTraits& traits, + const TriangleMesh& tm1, const TriangleMesh& tm2, + const VPM1& vpm1, const VPM2& vpm2, + const double h_lower_init, const double h_upper_init + ) : m_traits(traits), m_tm1(tm1), m_tm2(tm2), m_vpm1(vpm1), m_vpm2(vpm2) { - // Initialize the global bounds with 0., they will only grow. - h_local_upper = std::numeric_limits::infinity(); - h_local_lower_0 = std::numeric_limits::infinity(); - h_local_lower_1 = std::numeric_limits::infinity(); - h_local_lower_2 = std::numeric_limits::infinity(); + // Initialize the global bounds with infinity, triangles from TM2_tree + // try to minimize them. + h_local_lower = h_lower_init; + h_local_upper = h_upper_init; } // Explore the whole tree, i.e. always enter children if the methods @@ -152,6 +156,7 @@ namespace CGAL { } } + // Return the local Hausdorff bounds computed for the passed query triangle Hausdorff_bounds get_local_bounds() { return Hausdorff_bounds( h_local_lower, h_local_upper ); @@ -168,12 +173,9 @@ namespace CGAL { // Local Hausdorff bounds for the query triangle double h_local_upper; double h_local_lower; - // Local Hausdorff bounds for the query triangle's vertices - double h_local_lower_0; - double h_local_lower_1; - double h_local_lower_2; }; + /** * @class Hausdorff_primitive_traits_tm1 */ @@ -187,6 +189,8 @@ namespace CGAL { typedef ::CGAL::AABB_node Node; typedef typename Kernel::Point_3 Point_3; typedef typename Kernel::Triangle_3 Triangle_3; + typedef std::pair Candidate_triangle; + typedef typename std::vector Candidate_set; typedef AABB_tree< AABB_traits > TM2_tree; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; @@ -209,8 +213,18 @@ namespace CGAL { // Have reached a single triangle std::cout << "Reached Triangle in TM1: " << primitive.id() << '\n'; + // Map to transform faces of TM1 to actual triangles + Triangle_from_face_descriptor_map m_face_to_triangle_map( &m_tm1, m_vpm1 ); + Triangle_3 candidate_triangle = get(m_face_to_triangle_map, primitive.id()); + // Call Culling on B with the single triangle found. - Hausdorff_primitive_traits_tm2 traversal_traits_tm2( tm2_tree.traits(), m_tm1, m_tm2, m_vpm1, m_vpm2 ); + Hausdorff_primitive_traits_tm2< + Tree_traits, Primitive, Kernel, TriangleMesh, VPM1, VPM2 + > traversal_traits_tm2( + tm2_tree.traits(), m_tm1, m_tm2, m_vpm1, m_vpm2, + std::numeric_limits::infinity(), + std::numeric_limits::infinity() + ); tm2_tree.traversal(primitive, traversal_traits_tm2); // Update global Hausdorff bounds according to the obtained local bounds @@ -221,6 +235,12 @@ namespace CGAL { if (local_bounds.second > h_upper) { h_upper = local_bounds.second; } + + m_candidiate_triangles.push_back(Candidate_triangle(candidate_triangle, local_bounds)); + /* TODO Store the triangle given here is primitive as candidate triangle + together with the local bounds it optained to send it to + subdivision later. + */ } // Determine whether child nodes will still contribute to a larger @@ -253,6 +273,18 @@ namespace CGAL { } } + // Return those triangles from TM1 which are candidates for including a + // point realizing the Hausdorff distance + Candidate_set get_candidate_triangles() { + return m_candidiate_triangles; + } + + // Return the local Hausdorff bounds computed for the passed query triangle + Hausdorff_bounds get_global_bounds() + { + return Hausdorff_bounds( h_lower, h_upper ); + } + private: const AABBTraits& m_traits; // The two triangle meshes @@ -266,6 +298,8 @@ namespace CGAL { // Global Hausdorff bounds to be taken track of during the traversal double h_lower; double h_upper; + // List of candidate triangles with their Hausdorff bounds attached + Candidate_set m_candidiate_triangles; }; } From 7ec6c4ca328faefabf8faf766eec3aaccacd5735 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sat, 22 Jun 2019 09:02:27 +0200 Subject: [PATCH 15/48] Correct implementation of Culling on second mesh. --- ...traversal_traits_with_Hausdorff_distance.h | 42 +++++++++++++------ 1 file changed, 29 insertions(+), 13 deletions(-) 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 662399c4de1..16e9e49f39d 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 @@ -48,13 +48,16 @@ namespace CGAL { const AABBTraits& traits, const TriangleMesh& tm1, const TriangleMesh& tm2, const VPM1& vpm1, const VPM2& vpm2, - const double h_lower_init, const double h_upper_init + const double h_lower_init, const double h_upper_init, + const double h_v0_lower_init, const double h_v1_lower_init, const double h_v2_lower_init ) : m_traits(traits), m_tm1(tm1), m_tm2(tm2), m_vpm1(vpm1), m_vpm2(vpm2) { - // Initialize the global bounds with infinity, triangles from TM2_tree - // try to minimize them. + // Initialize the global and local bounds with the given values h_local_lower = h_lower_init; h_local_upper = h_upper_init; + h_v0_lower = h_v0_lower_init; + h_v1_lower = h_v1_lower_init; + h_v2_lower = h_v2_lower_init; } // Explore the whole tree, i.e. always enter children if the methods @@ -84,27 +87,33 @@ namespace CGAL { // Compute distances of the vertices to the primitive triangle in TM2 Triangle_from_face_descriptor_map face_to_triangle_map(&m_tm2, m_vpm2); - double v0_dist = squared_distance( + double v0_dist = approximate_sqrt(squared_distance( project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v0)), get(m_vpm1, v0) - ); - double v1_dist = squared_distance( + )); + if (v0_dist < h_v0_lower) h_v0_lower = v0_dist; + + double v1_dist = approximate_sqrt(squared_distance( project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v1)), get(m_vpm1, v1) - ); - double v2_dist = squared_distance( + )); + if (v1_dist < h_v1_lower) h_v1_lower = v1_dist; + + double v2_dist = approximate_sqrt(squared_distance( project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v2)), get(m_vpm1, v2) - ); + )); + if (v2_dist < h_v2_lower) h_v2_lower = v2_dist; // Get the distance as maximizers over all vertices - double distance = approximate_sqrt(std::max(std::max(v0_dist,v1_dist),v2_dist)); + double distance_upper = std::max(std::max(v0_dist,v1_dist),v2_dist); + double distance_lower = std::max(std::max(h_v0_lower,h_v1_lower),h_v2_lower); - // Since we are at the level of a single triangle in TM2, distance is + // Since we are at the level of a single triangle in TM2, distance_upper is // actually the correct Hausdorff distance from the query triangle in // TM1 to the primitive triangle in TM2 - if ( distance < h_local_upper ) h_local_upper = distance; - if ( distance < h_local_lower ) h_local_lower = distance; + if ( distance_upper < h_local_upper ) h_local_upper = distance_upper; + if ( distance_lower < h_local_lower ) h_local_lower = distance_lower; } // Determine whether child nodes will still contribute to a smaller @@ -173,6 +182,9 @@ namespace CGAL { // Local Hausdorff bounds for the query triangle double h_local_upper; double h_local_lower; + double h_v0_lower; + double h_v1_lower; + double h_v2_lower; }; @@ -223,6 +235,9 @@ namespace CGAL { > traversal_traits_tm2( tm2_tree.traits(), m_tm1, m_tm2, m_vpm1, m_vpm2, std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), std::numeric_limits::infinity() ); tm2_tree.traversal(primitive, traversal_traits_tm2); @@ -257,6 +272,7 @@ namespace CGAL { (bbox.min(1) + bbox.max(1)) / 2, (bbox.min(2) + bbox.max(2)) / 2); // Find the point from TM2 closest to the center + // TODO Insert a hint here to accelerate the query Point_3 closest = tm2_tree.closest_point(center); // Compute the distance of the center to the closest point in tm2 double dist = approximate_sqrt(squared_distance(center, closest)); From 2a523825a7502a5e6866e76ae4e2758a52989c17 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 24 Jun 2019 03:42:48 +0200 Subject: [PATCH 16/48] Implement a non-trivial toy example to test subdivision with. Intended solution is: Hausdorff distance will be attained at Point (0,0,1) and should be sqrt(3). --- ...usdorff_bounded_error_distance_example.cpp | 36 +++++++++++++++---- ...traversal_traits_with_Hausdorff_distance.h | 17 ++++----- 2 files changed, 39 insertions(+), 14 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index dd32a2b5e4e..eb68fec8e2b 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -10,21 +10,45 @@ #endif typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef K::Point_3 Point; -typedef CGAL::Surface_mesh Mesh; +typedef K::Point_3 Point_3; +typedef K::Triangle_3 Triangle_3; +typedef CGAL::Surface_mesh Mesh; +typedef Mesh::Vertex_index Vertex_index; namespace PMP = CGAL::Polygon_mesh_processing; int main() { Mesh tm1, tm2; - CGAL::make_tetrahedron(Point(.0,.0,.0), - Point(2,.0,.0), - Point(1,1,1), - Point(1,.0,2), +/* + CGAL::make_tetrahedron(Point_3(.0,.0,.0), + Point_3(2,.0,.0), + Point_3(1,1,1), + Point_3(1,.0,2), tm1); tm2=tm1; CGAL::Polygon_mesh_processing::isotropic_remeshing(tm2.faces(),.05, tm2); +*/ + tm1 = Mesh(); + tm1.add_vertex( Point_3(-1.,1.,1.) ); + tm1.add_vertex( Point_3(0.,-1.,1.) ); + tm1.add_vertex( Point_3(1.,1.,1.) ); + tm1.add_face( tm1.vertices() ); + + std::cout << "TM1 is valid: " << (tm1.is_valid() ? "true" : "false") << std::endl; + + tm2 = Mesh(); + Vertex_index w0 = tm2.add_vertex( Point_3(-1.,1.,0.) ); + Vertex_index w1 = tm2.add_vertex( Point_3(0.,-1.,0.) ); + Vertex_index w2 = tm2.add_vertex( Point_3(1.,1.,0.) ); + Vertex_index w3 = tm2.add_vertex( Point_3(0.,1.,-1.) ); + Vertex_index w4 = tm2.add_vertex( Point_3(-0.5,0.,-1.) ); + Vertex_index w5 = tm2.add_vertex( Point_3(0.5,0.,-1.) ); + tm2.add_face( w0, w3, w4 ); + tm2.add_face( w1, w4, w5 ); + tm2.add_face( w2, w5, w3 ); + + std::cout << "TM2 is valid: " << (tm2.is_valid() ? "true" : "false") << std::endl; std::cout << "Approximated Hausdorff distance: " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance 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 16e9e49f39d..2ffd206d9d1 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 @@ -67,8 +67,7 @@ namespace CGAL { // Compute the explicit Hausdorff distance to the given primitive void intersection(const Query& query, const Primitive& primitive) { - // Have reached a single triangle - std::cout << "Reached Triangle in TM2:" << primitive.id() << '\n'; + /* Have reached a single triangle, process it */ /* / Determine the distance accroding to @@ -222,8 +221,7 @@ namespace CGAL { // Compute the explicit Hausdorff distance to the given primitive void intersection(const Query& query, const Primitive& primitive) { - // Have reached a single triangle - std::cout << "Reached Triangle in TM1: " << primitive.id() << '\n'; + /* Have reached a single triangle, process it */ // Map to transform faces of TM1 to actual triangles Triangle_from_face_descriptor_map m_face_to_triangle_map( &m_tm1, m_vpm1 ); @@ -251,11 +249,14 @@ namespace CGAL { h_upper = local_bounds.second; } + std::cout << "Processed triangle: " << primitive.id() + << ", found bounds (low., up.): (" << local_bounds.first + << ", " << local_bounds.second << ")" << std::endl; + + // Store the triangle given as primitive here as candidate triangle + // together with the local bounds it obtained to sind it to subdivision + // later m_candidiate_triangles.push_back(Candidate_triangle(candidate_triangle, local_bounds)); - /* TODO Store the triangle given here is primitive as candidate triangle - together with the local bounds it optained to send it to - subdivision later. - */ } // Determine whether child nodes will still contribute to a larger From 68e9776af60865e723f746968b6d94be5472c180 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 24 Jun 2019 04:43:29 +0200 Subject: [PATCH 17/48] Implement subdivision of a triangle from the candidate set. --- .../CGAL/Polygon_mesh_processing/distance.h | 54 +++++++++++++++---- 1 file changed, 43 insertions(+), 11 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 45ab6e038dc..2a7a8f1e6dd 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -867,7 +867,7 @@ template traversal_traits_tm1( tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2 ); tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); + // TODO Implement the candidate_triangles set as Stack instead of Vector Candidate_set candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); + while ( (global_bounds.second - global_bounds.first > error_bound) && candidate_triangles.size() > 0 ) { + // Get the first triangle and its Hausdorff bounds from the candidate set + Candidate_triangle triangle_and_bound = candidate_triangles.front(); + // Remove it from the candidate set as it will be processed now + candidate_triangles.erase (candidate_triangles.begin(),candidate_triangles.begin()+1); + // Only process the triangle if it can contribute to the Hausdorff distance, + // i.e. if its Upper Bound is higher than the currently known best lower bound + if (triangle_and_bound.second.second > global_bounds.first) { + // Subdivide the triangle into four smaller triangles + Triangle_3 triangle_for_subdivision = triangle_and_bound.first; + Point_3 v0 = triangle_for_subdivision.vertex(0); + Point_3 v1 = triangle_for_subdivision.vertex(0); + Point_3 v2 = triangle_for_subdivision.vertex(0); + Point_3 v01 = Point_3( (v0.x()+ v1.x())/2., (v0.y()+v1.y())/2., (v0.z()+v1.z())/2. ); + Point_3 v02 = Point_3( (v0.x()+ v2.x())/2., (v0.y()+v2.y())/2., (v0.z()+v2.z())/2. ); + Point_3 v12 = Point_3( (v1.x()+ v2.x())/2., (v1.y()+v2.y())/2., (v1.z()+v2.z())/2. ); + Triangle_3 t0 = Triangle_3( v0, v01, v1); + Triangle_3 t1 = Triangle_3( v0, v02, v2); + Triangle_3 t2 = Triangle_3( v1, v12, v2); + Triangle_3 t3 = Triangle_3( v01, v02, v12); + + // TODO send each of the four triangles to Culling on B with the bounds of the parent triangle + // TODO add those triangles to the candidate set which still contribute to the Hausdorff distance + } + } + + // Print result found + std::cout << "Processing candidates finished, found distance (lower, upper): (" + << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; + + // Return linear interpolation between found upper and lower bound + return (global_bounds.first + global_bounds.second) / 2.; + #if !defined(CGAL_LINKED_WITH_TBB) CGAL_static_assertion_msg (!(boost::is_convertible::value), "Parallel_tag is enabled but TBB is unavailable."); @@ -938,6 +972,7 @@ double bounded_error_Hausdorff_impl( // } // else #endif +/* { // Store all vertices of tm1 in a vector std::vector tm1_vertices; @@ -955,11 +990,9 @@ double bounded_error_Hausdorff_impl( // on the Hausdorff measure Triangle_hausdorff_bounds thb = get(Face_property_tag(), tm1); - /* - / For each vertex in the first mesh, find the closest triangle in the - / second mesh, store it and also store the distance to this triangle - / in a dynamic vertex property - */ + // For each vertex in the first mesh, find the closest triangle in the + // second mesh, store it and also store the distance to this triangle + // in a dynamic vertex property for(vertex_descriptor vd : tm1_vertices) { // Get the point represented by the vertex @@ -982,11 +1015,9 @@ double bounded_error_Hausdorff_impl( // following std::vector candidate_triangles; - /* - / For each triangle in the first mesh, initialize its local upper and - / lower bound and store these in a dynamic face property for furture - / reference - */ + // For each triangle in the first mesh, initialize its local upper and + // lower bound and store these in a dynamic face property for furture + // reference for(face_descriptor fd : faces(tm1)) { // Initialize the local bounds for the current face fd @@ -1056,6 +1087,7 @@ double bounded_error_Hausdorff_impl( return (CGAL::approximate_sqrt(h_lower)+CGAL::approximate_sqrt(h_upper))/2.; } +*/ } } //end of namespace internal From 17b00036a90ceda4f66879b19d214e8d901cf763 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 24 Jun 2019 04:48:43 +0200 Subject: [PATCH 18/48] Fix typo in vertex lookup. --- .../include/CGAL/Polygon_mesh_processing/distance.h | 4 ++-- 1 file changed, 2 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 2a7a8f1e6dd..cc3759de006 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -934,8 +934,8 @@ double bounded_error_Hausdorff_impl( // Subdivide the triangle into four smaller triangles Triangle_3 triangle_for_subdivision = triangle_and_bound.first; Point_3 v0 = triangle_for_subdivision.vertex(0); - Point_3 v1 = triangle_for_subdivision.vertex(0); - Point_3 v2 = triangle_for_subdivision.vertex(0); + Point_3 v1 = triangle_for_subdivision.vertex(1); + Point_3 v2 = triangle_for_subdivision.vertex(2); Point_3 v01 = Point_3( (v0.x()+ v1.x())/2., (v0.y()+v1.y())/2., (v0.z()+v1.z())/2. ); Point_3 v02 = Point_3( (v0.x()+ v2.x())/2., (v0.y()+v2.y())/2., (v0.z()+v2.z())/2. ); Point_3 v12 = Point_3( (v1.x()+ v2.x())/2., (v1.y()+v2.y())/2., (v1.z()+v2.z())/2. ); From 3c73230272cbf9200bbca9f3ac019b93d7979946 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 24 Jun 2019 23:41:12 +0200 Subject: [PATCH 19/48] Finish first working prototype giving correct answer on the toy example. --- .../CGAL/Polygon_mesh_processing/distance.h | 77 +++++++++++++++++-- ...traversal_traits_with_Hausdorff_distance.h | 54 +++++-------- 2 files changed, 90 insertions(+), 41 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 cc3759de006..88cfd3d2f70 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -880,6 +880,7 @@ double bounded_error_Hausdorff_impl( typedef AABB_tree< AABB_traits > TM1_tree; typedef AABB_tree< AABB_traits > TM2_tree; typedef typename AABB_tree< AABB_traits >::AABB_traits Tree_traits; + typedef typename Tree_traits::Point_and_primitive_id Point_and_primitive_id; typedef typename Kernel::Point_3 Point_3; typedef typename Kernel::Triangle_3 Triangle_3; @@ -924,28 +925,88 @@ double bounded_error_Hausdorff_impl( Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); while ( (global_bounds.second - global_bounds.first > error_bound) && candidate_triangles.size() > 0 ) { + + // TODO Why does it only stop because of triangles, not because of the global bound condition? + std::cout << "There are currently " << candidate_triangles.size() + << " candidates. Global bounds are: " + << global_bounds.first << ", " << global_bounds.second << std::endl; + // Get the first triangle and its Hausdorff bounds from the candidate set Candidate_triangle triangle_and_bound = candidate_triangles.front(); // Remove it from the candidate set as it will be processed now candidate_triangles.erase (candidate_triangles.begin(),candidate_triangles.begin()+1); + // Only process the triangle if it can contribute to the Hausdorff distance, // i.e. if its Upper Bound is higher than the currently known best lower bound - if (triangle_and_bound.second.second > global_bounds.first) { - // Subdivide the triangle into four smaller triangles + // and the difference between the bounds to be obtained is larger than the + // user given error. + Hausdorff_bounds triangle_bounds = triangle_and_bound.second; + if ( (triangle_bounds.second > global_bounds.first) && (triangle_bounds.second - triangle_bounds.first > error_bound) ) { + // Get the triangle that is to be subdivided and read its vertices Triangle_3 triangle_for_subdivision = triangle_and_bound.first; Point_3 v0 = triangle_for_subdivision.vertex(0); Point_3 v1 = triangle_for_subdivision.vertex(1); Point_3 v2 = triangle_for_subdivision.vertex(2); + + // Check second stopping condition: All three vertices of the triangle + // are projected onto the same triangle in TM2 + Point_and_primitive_id closest_triangle_v0 = tm2_tree.closest_point_and_primitive(v0); + Point_and_primitive_id closest_triangle_v1 = tm2_tree.closest_point_and_primitive(v1); + Point_and_primitive_id closest_triangle_v2 = tm2_tree.closest_point_and_primitive(v2); + 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. + global_bounds.first = triangle_bounds.second; + } + + // Subdivide the triangle into four smaller triangles Point_3 v01 = Point_3( (v0.x()+ v1.x())/2., (v0.y()+v1.y())/2., (v0.z()+v1.z())/2. ); Point_3 v02 = Point_3( (v0.x()+ v2.x())/2., (v0.y()+v2.y())/2., (v0.z()+v2.z())/2. ); Point_3 v12 = Point_3( (v1.x()+ v2.x())/2., (v1.y()+v2.y())/2., (v1.z()+v2.z())/2. ); - Triangle_3 t0 = Triangle_3( v0, v01, v1); - Triangle_3 t1 = Triangle_3( v0, v02, v2); - Triangle_3 t2 = Triangle_3( v1, v12, v2); - Triangle_3 t3 = Triangle_3( v01, v02, v12); + std::array sub_triangles = { + Triangle_3( v0, v01, v02), Triangle_3( v1, v01, v12), + Triangle_3( v2, v02, v12), Triangle_3( v01, v02, v12) + }; - // TODO send each of the four triangles to Culling on B with the bounds of the parent triangle - // TODO add those triangles to the candidate set which still contribute to the Hausdorff distance + // Send each of the four triangles to Culling on B with the bounds of the parent triangle + for (int i=0; i<4; i++) { + // Call Culling on B with the single triangle found. + Hausdorff_primitive_traits_tm2 traversal_traits_tm2( + tm2_tree.traits(), tm2, vpm2, + triangle_bounds.first, + triangle_bounds.second, + std::numeric_limits::infinity(), + std::numeric_limits::infinity(), + std::numeric_limits::infinity() + ); + tm2_tree.traversal(sub_triangles[i], traversal_traits_tm2); + + // Get the highest current bound from all candidate triangles + double current_max = 0.; + for(auto&& ct: candidate_triangles) { + if (ct.second.second > current_max) { + current_max = ct.second.second; + } + } + + // Update global Hausdorff bounds according to the obtained local bounds + Hausdorff_bounds local_bounds = traversal_traits_tm2.get_local_bounds(); + if (local_bounds.first > global_bounds.first) { + global_bounds.first = local_bounds.first; + } + global_bounds.second = std::max(current_max, local_bounds.second); + + // Add the subtriangle to the candidate list + candidate_triangles.push_back(Candidate_triangle(sub_triangles[i], local_bounds)); + + // std::cout << "Split triangle (" << v0 << ", " << v1 << ", " << v2 + // << ") with bounds: (" << triangle_bounds.first << ", " + // << triangle_bounds.second << "), sub-triangle " << i + // << " (" << sub_triangles[i].vertex(0) << ", " << sub_triangles[i].vertex(1) << ", " << sub_triangles[i].vertex(2) + // << ") has bounds: (" + // << local_bounds.first << ", " << local_bounds.second << "), gobal bounds are: (" + // << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; + } } } 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 2ffd206d9d1..950d8ab1f38 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 @@ -32,26 +32,24 @@ namespace CGAL { /** * @class Hausdorff_primitive_traits_tm2 */ - template + template class Hausdorff_primitive_traits_tm2 { typedef typename AABBTraits::Primitive Primitive; typedef ::CGAL::AABB_node Node; typedef typename AABBTraits::Bounding_box Bounding_box; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typename Kernel::Construct_projected_point_3 project_point; typedef typename Kernel::Point_3 Point_3; public: Hausdorff_primitive_traits_tm2( const AABBTraits& traits, - const TriangleMesh& tm1, const TriangleMesh& tm2, - const VPM1& vpm1, const VPM2& vpm2, + const TriangleMesh& tm2, + const VPM2& vpm2, const double h_lower_init, const double h_upper_init, const double h_v0_lower_init, const double h_v1_lower_init, const double h_v2_lower_init ) - : m_traits(traits), m_tm1(tm1), m_tm2(tm2), m_vpm1(vpm1), m_vpm2(vpm2) { + : m_traits(traits), m_tm2(tm2), m_vpm2(vpm2) { // Initialize the global and local bounds with the given values h_local_lower = h_lower_init; h_local_upper = h_upper_init; @@ -79,29 +77,22 @@ namespace CGAL { */ // The query object is a triangle from TM1, get its vertices - halfedge_descriptor hd = halfedge(query.id(), m_tm1); - vertex_descriptor v0 = source(hd, m_tm1); - vertex_descriptor v1 = target(hd, m_tm1); - vertex_descriptor v2 = target(next(hd, m_tm1), m_tm1); + Point_3 v0 = query.vertex(0); + Point_3 v1 = query.vertex(1); + Point_3 v2 = query.vertex(2); // Compute distances of the vertices to the primitive triangle in TM2 Triangle_from_face_descriptor_map face_to_triangle_map(&m_tm2, m_vpm2); double v0_dist = approximate_sqrt(squared_distance( - project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v0)), - get(m_vpm1, v0) - )); + project_point( get(face_to_triangle_map, primitive.id()), v0), v0 ) ); if (v0_dist < h_v0_lower) h_v0_lower = v0_dist; double v1_dist = approximate_sqrt(squared_distance( - project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v1)), - get(m_vpm1, v1) - )); + project_point( get(face_to_triangle_map, primitive.id()), v1), v1 ) ); if (v1_dist < h_v1_lower) h_v1_lower = v1_dist; double v2_dist = approximate_sqrt(squared_distance( - project_point(get(face_to_triangle_map, primitive.id()), get(m_vpm1, v2)), - get(m_vpm1, v2) - )); + project_point( get(face_to_triangle_map, primitive.id()), v2), v2 ) ); if (v2_dist < h_v2_lower) h_v2_lower = v2_dist; // Get the distance as maximizers over all vertices @@ -131,10 +122,9 @@ namespace CGAL { // Get the center of the query triangle // The query object is a triangle from TM1, get its vertices - halfedge_descriptor hd = halfedge(query.id(), m_tm1); - Point_3 v0 = get(m_vpm1, source(hd, m_tm1)); - Point_3 v1 = get(m_vpm1, target(hd, m_tm1)); - Point_3 v2 = get(m_vpm1, target(next(hd, m_tm1), m_tm1)); + Point_3 v0 = query.vertex(0); + Point_3 v1 = query.vertex(1); + Point_3 v2 = query.vertex(2); // Compute the barycenter of the triangle Point_3 tri_center = Point_3( 0.3*(v0.x()+v1.x()+v2.x()), 0.3*(v0.y()+v1.y()+v2.y()), 0.3*(v0.z()+v1.z()+v2.z()) ); @@ -172,11 +162,9 @@ namespace CGAL { private: const AABBTraits& m_traits; - // The two triangle meshes - const TriangleMesh& m_tm1; + // the mesh of this tree const TriangleMesh& m_tm2; - // Their respective vertex-point Maps - const VPM1& m_vpm1; + // its vertex point map const VPM2& m_vpm2; // Local Hausdorff bounds for the query triangle double h_local_upper; @@ -208,7 +196,7 @@ namespace CGAL { public: Hausdorff_primitive_traits_tm1(const AABBTraits& traits, const TM2_tree& tree, const TriangleMesh& tm1, const TriangleMesh& tm2 , const VPM1& vpm1, const VPM2& vpm2 ) - : m_traits(traits), tm2_tree(tree), m_tm1(tm1), m_tm2(tm2), m_vpm1(vpm1), m_vpm2(vpm2) { + : m_traits(traits), m_tm2_tree(tree), m_tm1(tm1), m_tm2(tm2), m_vpm1(vpm1), m_vpm2(vpm2) { // Initialize the global bounds with 0., they will only grow. h_lower = 0.; h_upper = 0.; @@ -229,16 +217,16 @@ namespace CGAL { // Call Culling on B with the single triangle found. Hausdorff_primitive_traits_tm2< - Tree_traits, Primitive, Kernel, TriangleMesh, VPM1, VPM2 + Tree_traits, Triangle_3, Kernel, TriangleMesh, VPM2 > traversal_traits_tm2( - tm2_tree.traits(), m_tm1, m_tm2, m_vpm1, m_vpm2, + m_tm2_tree.traits(), m_tm2, m_vpm2, std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity() ); - tm2_tree.traversal(primitive, traversal_traits_tm2); + m_tm2_tree.traversal(candidate_triangle, traversal_traits_tm2); // Update global Hausdorff bounds according to the obtained local bounds Hausdorff_bounds local_bounds = traversal_traits_tm2.get_local_bounds(); @@ -274,7 +262,7 @@ namespace CGAL { (bbox.min(2) + bbox.max(2)) / 2); // Find the point from TM2 closest to the center // TODO Insert a hint here to accelerate the query - Point_3 closest = tm2_tree.closest_point(center); + Point_3 closest = m_tm2_tree.closest_point(center); // Compute the distance of the center to the closest point in tm2 double dist = approximate_sqrt(squared_distance(center, closest)); // Compute the radius of the circumsphere of the bounding boxes @@ -311,7 +299,7 @@ namespace CGAL { const VPM1 m_vpm1; const VPM2 m_vpm2; // AABB tree for the second triangle meshes - const TM2_tree& tm2_tree; + const TM2_tree& m_tm2_tree; // Global Hausdorff bounds to be taken track of during the traversal double h_lower; double h_upper; From d8b6e7dfb1ccea8746d53f6c25be2fb8efffa2ca Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 12 Jul 2019 18:04:08 +0200 Subject: [PATCH 20/48] Results of Code Review with Sebastien. --- ...usdorff_bounded_error_distance_example.cpp | 4 ++++ .../CGAL/Polygon_mesh_processing/distance.h | 20 +++++++++---------- ...traversal_traits_with_Hausdorff_distance.h | 9 +++++++-- 3 files changed, 21 insertions(+), 12 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index eb68fec8e2b..0c84ea03ca2 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -20,6 +20,10 @@ namespace PMP = CGAL::Polygon_mesh_processing; int main() { Mesh tm1, tm2; + // TODO Read real meshes and try method with it. + // std::ofstream(argv[1]) >> tm1; + // http://www.pmp-book.org/download/meshes/iphi_fullres.zip + // https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__meshing__grp.html#ga028a80dc84395650f67714fa7618ec53 /* CGAL::make_tetrahedron(Point_3(.0,.0,.0), Point_3(2,.0,.0), 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 88cfd3d2f70..1e74027ecf3 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -916,25 +916,23 @@ double bounded_error_Hausdorff_impl( tm2_tree.accelerate_distance_queries(); std::pair hint = tm2_tree.any_reference_point_and_id(); - // Build traversal traits for tm1_tree and tm2_tree + // Build traversal traits for tm1_tree Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2 ); - tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); + // Find candidate triangles in TM1 which might realise the Hausdorff bound + tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); // dummy point given as query as not needed // TODO Implement the candidate_triangles set as Stack instead of Vector + // check: https://www.boost.org/doc/libs/1_55_0/doc/html/heap.html + // Can already build a sorted structure while collecting the candidates Candidate_set candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); while ( (global_bounds.second - global_bounds.first > error_bound) && candidate_triangles.size() > 0 ) { - // TODO Why does it only stop because of triangles, not because of the global bound condition? - std::cout << "There are currently " << candidate_triangles.size() - << " candidates. Global bounds are: " - << global_bounds.first << ", " << global_bounds.second << std::endl; - // Get the first triangle and its Hausdorff bounds from the candidate set - Candidate_triangle triangle_and_bound = candidate_triangles.front(); + Candidate_triangle triangle_and_bound = candidate_triangles.back(); // Remove it from the candidate set as it will be processed now - candidate_triangles.erase (candidate_triangles.begin(),candidate_triangles.begin()+1); + candidate_triangles.pop_back(); // Only process the triangle if it can contribute to the Hausdorff distance, // i.e. if its Upper Bound is higher than the currently known best lower bound @@ -960,6 +958,7 @@ double bounded_error_Hausdorff_impl( } // Subdivide the triangle into four smaller triangles + // TODO Use CGAL::midpoint here. Point_3 v01 = Point_3( (v0.x()+ v1.x())/2., (v0.y()+v1.y())/2., (v0.z()+v1.z())/2. ); Point_3 v02 = Point_3( (v0.x()+ v2.x())/2., (v0.y()+v2.y())/2., (v0.z()+v2.z())/2. ); Point_3 v12 = Point_3( (v1.x()+ v2.x())/2., (v1.y()+v2.y())/2., (v1.z()+v2.z())/2. ); @@ -996,6 +995,7 @@ double bounded_error_Hausdorff_impl( } global_bounds.second = std::max(current_max, local_bounds.second); + // TODO Additionally store the face descriptor of the parent from TM1 in the Candidate_triangle. // Add the subtriangle to the candidate list candidate_triangles.push_back(Candidate_triangle(sub_triangles[i], local_bounds)); @@ -1212,7 +1212,7 @@ double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, Vpm2 vpm2 = choose_param(get_param(np2, internal_np::vertex_point), get_const_property_map(vertex_point, tm2)); - return internal::bounded_error_Hausdorff_impl(tm1, tm2, error_bound*error_bound, vpm1, vpm2); + return internal::bounded_error_Hausdorff_impl(tm1, tm2, error_bound, vpm1, vpm2); } template< class Concurrency_tag, 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 950d8ab1f38..818430dbeb4 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 @@ -125,7 +125,8 @@ namespace CGAL { Point_3 v0 = query.vertex(0); Point_3 v1 = query.vertex(1); Point_3 v2 = query.vertex(2); - // Compute the barycenter of the triangle + // Compute the centroid of the triangle + // TODO Use the centroid() method from the kernel here. Point_3 tri_center = Point_3( 0.3*(v0.x()+v1.x()+v2.x()), 0.3*(v0.y()+v1.y()+v2.y()), 0.3*(v0.z()+v1.z()+v2.z()) ); // Compute the distance of the center to the closest point in tm2 @@ -249,10 +250,14 @@ namespace CGAL { // Determine whether child nodes will still contribute to a larger // Hausdorff distance and thus have to be entered - bool do_intersect(const Query& query, const Node& node) const + // TODO Document Query object, explain why I don't need it here. + bool do_intersect(const Query& /*query*/, const Node& node) const { /* Have reached a node, determine whether or not to enter it */ + // TODO What's the closest distance of TM2 to the box given by node? + // Can we have a sharper bound on this than the one implemented below? + // Get the bounding box of the nodes Bounding_box bbox = node.bbox(); // Compute its center From 519d37d9d7bef230016c46a43a16ad9bf746c8f1 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 12 Jul 2019 18:17:24 +0200 Subject: [PATCH 21/48] Use Kernel method instead of computations by hand. --- .../include/CGAL/Polygon_mesh_processing/distance.h | 7 +++---- .../AABB_traversal_traits_with_Hausdorff_distance.h | 7 +------ 2 files changed, 4 insertions(+), 10 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 1e74027ecf3..2ef4798ff4c 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -958,10 +958,9 @@ double bounded_error_Hausdorff_impl( } // Subdivide the triangle into four smaller triangles - // TODO Use CGAL::midpoint here. - Point_3 v01 = Point_3( (v0.x()+ v1.x())/2., (v0.y()+v1.y())/2., (v0.z()+v1.z())/2. ); - Point_3 v02 = Point_3( (v0.x()+ v2.x())/2., (v0.y()+v2.y())/2., (v0.z()+v2.z())/2. ); - Point_3 v12 = Point_3( (v1.x()+ v2.x())/2., (v1.y()+v2.y())/2., (v1.z()+v2.z())/2. ); + Point_3 v01 = midpoint( v0, v1 ); + Point_3 v02 = midpoint( v0, v2 ); + Point_3 v12 = midpoint( v1, v2 ); std::array sub_triangles = { Triangle_3( v0, v01, v02), Triangle_3( v1, v01, v12), Triangle_3( v2, v02, v12), Triangle_3( v01, v02, v12) 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 818430dbeb4..6393c74a912 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 @@ -126,8 +126,7 @@ namespace CGAL { Point_3 v1 = query.vertex(1); Point_3 v2 = query.vertex(2); // Compute the centroid of the triangle - // TODO Use the centroid() method from the kernel here. - Point_3 tri_center = Point_3( 0.3*(v0.x()+v1.x()+v2.x()), 0.3*(v0.y()+v1.y()+v2.y()), 0.3*(v0.z()+v1.z()+v2.z()) ); + Point_3 tri_center = centroid( query ); // Compute the distance of the center to the closest point in tm2 double dist = approximate_sqrt(squared_distance(bb_center, tri_center)); @@ -238,10 +237,6 @@ namespace CGAL { h_upper = local_bounds.second; } - std::cout << "Processed triangle: " << primitive.id() - << ", found bounds (low., up.): (" << local_bounds.first - << ", " << local_bounds.second << ")" << std::endl; - // Store the triangle given as primitive here as candidate triangle // together with the local bounds it obtained to sind it to subdivision // later From 0332d9f00f02a275cf7378904287126e28e4d1b2 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 12 Jul 2019 18:55:55 +0200 Subject: [PATCH 22/48] Enable reading of real meshes and perturbation of them for distance computation. --- ...usdorff_bounded_error_distance_example.cpp | 25 +++++++++++++++---- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index 0c84ea03ca2..6f162f97058 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #if defined(CGAL_LINKED_WITH_TBB) #define TAG CGAL::Parallel_tag @@ -17,14 +18,11 @@ typedef Mesh::Vertex_index Vertex_index; namespace PMP = CGAL::Polygon_mesh_processing; -int main() +int main(int argc, char** argv) { Mesh tm1, tm2; - // TODO Read real meshes and try method with it. - // std::ofstream(argv[1]) >> tm1; - // http://www.pmp-book.org/download/meshes/iphi_fullres.zip - // https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__meshing__grp.html#ga028a80dc84395650f67714fa7618ec53 /* + // Easy Example CGAL::make_tetrahedron(Point_3(.0,.0,.0), Point_3(2,.0,.0), Point_3(1,1,1), @@ -33,6 +31,10 @@ int main() tm2=tm1; CGAL::Polygon_mesh_processing::isotropic_remeshing(tm2.faces(),.05, tm2); */ + +/* + // Example with point realizing the Hausdorff distance strictly lying in the + // interior of a triangle tm1 = Mesh(); tm1.add_vertex( Point_3(-1.,1.,1.) ); tm1.add_vertex( Point_3(0.,-1.,1.) ); @@ -53,6 +55,19 @@ int main() tm2.add_face( w2, w5, w3 ); std::cout << "TM2 is valid: " << (tm2.is_valid() ? "true" : "false") << std::endl; +*/ + + // Read a real mesh given by the user + std::ifstream input( argv[1] ); + input >> tm1; + std::cout << "Read a mesh with " << tm1.number_of_faces() << " triangles." << std::endl; + // Copy the mesh and perturb it slightly + tm2 = tm1; + bool do_project = false; + CGAL::Polygon_mesh_processing::random_perturbation( tm2.vertices(), tm2, 0.001, do_project ); + std::cout << "Perturbed the input mesh, now computing the Hausdorff distance." << std::endl; + +// https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__meshing__grp.html#ga028a80dc84395650f67714fa7618ec53 std::cout << "Approximated Hausdorff distance: " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance From e78cbff8a14cfe583a796e82d8758345f08b25ac Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 15 Jul 2019 14:12:01 +0200 Subject: [PATCH 23/48] Implement naive bounded Hausdoff computation by simple subdivision. --- ...usdorff_bounded_error_distance_example.cpp | 11 +- .../CGAL/Polygon_mesh_processing/distance.h | 158 ++++++++++++++++++ 2 files changed, 168 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index 6f162f97058..8280052edec 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -61,16 +61,25 @@ int main(int argc, char** argv) std::ifstream input( argv[1] ); input >> tm1; std::cout << "Read a mesh with " << tm1.number_of_faces() << " triangles." << std::endl; + std::ifstream input2( argv[2] ); + input2 >> tm2; + std::cout << "Read a mesh with " << tm2.number_of_faces() << " triangles." << std::endl; // Copy the mesh and perturb it slightly +/* tm2 = tm1; bool do_project = false; CGAL::Polygon_mesh_processing::random_perturbation( tm2.vertices(), tm2, 0.001, do_project ); std::cout << "Perturbed the input mesh, now computing the Hausdorff distance." << std::endl; +*/ // https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__meshing__grp.html#ga028a80dc84395650f67714fa7618ec53 std::cout << "Approximated Hausdorff distance: " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance - (tm1, tm2, 0.001) + (tm1, tm2, 0.01) + << std::endl; + std::cout << "Approximated Hausdorff distance (naive): " + << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance_naive + (tm1, tm2, 0.01) << std::endl; } 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 2ef4798ff4c..df1279220db 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1150,6 +1150,109 @@ double bounded_error_Hausdorff_impl( */ } +template +double recursive_hausdorff_subdivision( + const Point_3& v0, + const Point_3& v1, + const Point_3& v2, + const TM2_tree& tm2_tree, + const typename Kernel::FT& squared_error_bound) +{ +// std::cout << "Processing points " << v0 << ", " << v1 << ", " << v2 << std::endl; + + // If any edge length of the triangle is larger than (error_bound * 1/sqrt(3)), subdivide the triangle and proceed recursively + double max_edge_length = std::max( std::max( squared_distance( v0, v1 ), squared_distance( v0, v2 )), squared_distance( v1, v2 )); +// std::cout << "Maximum edge lenght: " << max_edge_length << "; Squared error bound: " << squared_error_bound << std::endl; + if ( max_edge_length < squared_error_bound ) { +// std::cout << "Maximum distance: " << std::max(std::max(squared_distance( v0, tm2_tree.closest_point(v0) ),squared_distance( v1, tm2_tree.closest_point(v1) ) ), squared_distance( v2, tm2_tree.closest_point(v2) )); + return std::max( + std::max( + squared_distance( v0, tm2_tree.closest_point(v0) ), + squared_distance( v1, tm2_tree.closest_point(v1) ) ), + squared_distance( v2, tm2_tree.closest_point(v2) ) + ); + } + + // Else return maximum of the distances of the three points to TM2 (via TM2_tree). + Point_3 v01 = midpoint( v0, v1 ); + Point_3 v02 = midpoint( v0, v2 ); + Point_3 v12 = midpoint( v1, v2 ); + + return std::max ( + std::max( + recursive_hausdorff_subdivision( v0,v01,v02,tm2_tree,squared_error_bound ), + recursive_hausdorff_subdivision( v1,v01,v12,tm2_tree,squared_error_bound ) + ), + std::max( + recursive_hausdorff_subdivision( v2,v02,v12,tm2_tree,squared_error_bound ), + recursive_hausdorff_subdivision( v01,v02,v12,tm2_tree,squared_error_bound ) + ) + ); +} + +template +double bounded_error_Hausdorff_naive_impl( + const TriangleMesh& tm1, + const TriangleMesh& tm2, + const typename Kernel::FT& error_bound, + VPM1 vpm1, + VPM2 vpm2) +{ + CGAL_assertion_code( bool is_triangle = is_triangle_mesh(tm1) && is_triangle_mesh(tm2) ); + CGAL_assertion_msg (is_triangle, + "One of the meshes is not triangulated. Distance computing impossible."); + + typedef AABB_face_graph_triangle_primitive TM2_primitive; + typedef AABB_tree< AABB_traits > TM2_tree; + + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename Kernel::Point_3 Point_3; + typedef typename Kernel::Triangle_3 Triangle_3; + + double squared_lower_bound = 0.; + double squared_error_bound = error_bound * error_bound; + + // Build an AABB tree on tm2 + TM2_tree tm2_tree( faces(tm2).begin(), faces(tm2).end(), tm2, vpm2 ); + tm2_tree.build(); + tm2_tree.accelerate_distance_queries(); + + Triangle_from_face_descriptor_map face_to_triangle_map( &tm1, vpm1 ); + + // Iterate over the triangles of TM1. + for(face_descriptor fd : faces(tm1)) + { + // Get the vertices of the face and pass them on to a recursive method. + Triangle_3 triangle = get(face_to_triangle_map, fd); + Point_3 v0 = triangle.vertex(0); + Point_3 v1 = triangle.vertex(1); + Point_3 v2 = triangle.vertex(2); + + double triangle_bound = recursive_hausdorff_subdivision( v0, v1, v2, tm2_tree, squared_error_bound ); + + if( triangle_bound > squared_lower_bound ) { + squared_lower_bound = triangle_bound; + } + } + + // Return linear interpolation between found upper and lower bound + return (approximate_sqrt( squared_lower_bound )); + +#if !defined(CGAL_LINKED_WITH_TBB) + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + // TODO implement parallelized version of the below here. +#endif +} + } //end of namespace internal /** @@ -1187,6 +1290,11 @@ double bounded_error_Hausdorff_impl( * The function `CGAL::parameters::all_default()` can be used to indicate to use the default values for * `np1` and specify custom values for `np2` */ + +/* + * Implementation of Bounded Hausdorff distance computation using AABBTree + * culling. + */ template< class Concurrency_tag, class TriangleMesh, class NamedParameters1, @@ -1234,6 +1342,56 @@ double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, return bounded_error_Hausdorff_distance(tm1, tm2, error_bound, parameters::all_default() ); } +/* + * Implementation of naive Bounded Hausdorff distance computation. + */ +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1, + class NamedParameters2> +double bounded_error_Hausdorff_distance_naive( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double error_bound, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + typedef typename GetGeomTraits::type Geom_traits; + + typedef typename GetVertexPointMap::const_type Vpm1; + typedef typename GetVertexPointMap::const_type Vpm2; + + using boost::choose_param; + using boost::get_param; + + Vpm1 vpm1 = choose_param(get_param(np1, internal_np::vertex_point), + get_const_property_map(vertex_point, tm1)); + Vpm2 vpm2 = choose_param(get_param(np2, internal_np::vertex_point), + get_const_property_map(vertex_point, tm2)); + + return internal::bounded_error_Hausdorff_naive_impl(tm1, tm2, error_bound, vpm1, vpm2); +} + +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1> +double bounded_error_Hausdorff_distance_naive( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double error_bound, + const NamedParameters1& np1) +{ + return bounded_error_Hausdorff_distance_naive(tm1, tm2, error_bound, np1, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh> +double bounded_error_Hausdorff_distance_naive( const TriangleMesh& tm1, + const TriangleMesh& tm2, + double error_bound) +{ + return bounded_error_Hausdorff_distance_naive(tm1, tm2, error_bound, parameters::all_default() ); +} + } } // end of namespace CGAL::Polygon_mesh_processing From 000604faa9a7641421f91030d6428cc1083fc830 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 15 Jul 2019 14:20:15 +0200 Subject: [PATCH 24/48] Comment naive comparison method. --- ...usdorff_bounded_error_distance_example.cpp | 7 +++-- .../CGAL/Polygon_mesh_processing/distance.h | 27 ++++++++++++------- 2 files changed, 23 insertions(+), 11 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index 8280052edec..0242dcfe4fd 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -74,12 +74,15 @@ int main(int argc, char** argv) // https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__meshing__grp.html#ga028a80dc84395650f67714fa7618ec53 + double error_bound = 0.01; + std::cout << "Approximated Hausdorff distance: " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance - (tm1, tm2, 0.01) + (tm1, tm2, error_bound) << std::endl; std::cout << "Approximated Hausdorff distance (naive): " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance_naive - (tm1, tm2, 0.01) + (tm1, tm2, error_bound) + << ", the actual distance is at most " << error_bound << " larger." << std::endl; } 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 df1279220db..25dba12f4b1 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1013,7 +1013,7 @@ double bounded_error_Hausdorff_impl( std::cout << "Processing candidates finished, found distance (lower, upper): (" << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; - // Return linear interpolation between found upper and lower bound + // Return linear interpolation between found lower and upper bound return (global_bounds.first + global_bounds.second) / 2.; #if !defined(CGAL_LINKED_WITH_TBB) @@ -1160,13 +1160,16 @@ double recursive_hausdorff_subdivision( const TM2_tree& tm2_tree, const typename Kernel::FT& squared_error_bound) { -// std::cout << "Processing points " << v0 << ", " << v1 << ", " << v2 << std::endl; - - // If any edge length of the triangle is larger than (error_bound * 1/sqrt(3)), subdivide the triangle and proceed recursively - double max_edge_length = std::max( std::max( squared_distance( v0, v1 ), squared_distance( v0, v2 )), squared_distance( v1, v2 )); -// std::cout << "Maximum edge lenght: " << max_edge_length << "; Squared error bound: " << squared_error_bound << std::endl; - if ( max_edge_length < squared_error_bound ) { -// std::cout << "Maximum distance: " << std::max(std::max(squared_distance( v0, tm2_tree.closest_point(v0) ),squared_distance( v1, tm2_tree.closest_point(v1) ) ), squared_distance( v2, tm2_tree.closest_point(v2) )); + // If all edge lengths of the triangle are below the error_bound, + // return maximum of the distances of the three points to TM2 (via TM2_tree). + double max_squared_edge_length = + std::max( + std::max( + squared_distance( v0, v1 ), + squared_distance( v0, v2 )), + squared_distance( v1, v2 ) + ); + if ( max_squared_edge_length < squared_error_bound ) { return std::max( std::max( squared_distance( v0, tm2_tree.closest_point(v0) ), @@ -1175,7 +1178,7 @@ double recursive_hausdorff_subdivision( ); } - // Else return maximum of the distances of the three points to TM2 (via TM2_tree). + // Else subdivide the triangle and proceed recursively Point_3 v01 = midpoint( v0, v1 ); Point_3 v02 = midpoint( v0, v2 ); Point_3 v12 = midpoint( v1, v2 ); @@ -1216,7 +1219,9 @@ double bounded_error_Hausdorff_naive_impl( typedef typename Kernel::Point_3 Point_3; typedef typename Kernel::Triangle_3 Triangle_3; + // Initially, no lower bound is known double squared_lower_bound = 0.; + // Work with squares in the following, only draw sqrt at the very end double squared_error_bound = error_bound * error_bound; // Build an AABB tree on tm2 @@ -1224,6 +1229,7 @@ double bounded_error_Hausdorff_naive_impl( tm2_tree.build(); tm2_tree.accelerate_distance_queries(); + // Build a map to obtain actual triangles from the face descriptors of tm1. Triangle_from_face_descriptor_map face_to_triangle_map( &tm1, vpm1 ); // Iterate over the triangles of TM1. @@ -1235,8 +1241,11 @@ double bounded_error_Hausdorff_naive_impl( Point_3 v1 = triangle.vertex(1); Point_3 v2 = triangle.vertex(2); + // Recursively process the current triangle to obtain a lower bound on + // its Hausdorff distance. double triangle_bound = recursive_hausdorff_subdivision( v0, v1, v2, tm2_tree, squared_error_bound ); + // Store the largest lower bound. if( triangle_bound > squared_lower_bound ) { squared_lower_bound = triangle_bound; } From 313402589a3fc3cdecdb32615119063e55826f58 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 15 Jul 2019 14:54:52 +0200 Subject: [PATCH 25/48] Bugfix computation of upper bound, should not go below lower bound. --- .../include/CGAL/Polygon_mesh_processing/distance.h | 7 ++++++- 1 file changed, 6 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 25dba12f4b1..d2957a78576 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -954,7 +954,9 @@ double bounded_error_Hausdorff_impl( 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. + // TODO Update the reference to the realizing triangle here as this is the best current guess. global_bounds.first = triangle_bounds.second; + continue; } // Subdivide the triangle into four smaller triangles @@ -992,7 +994,10 @@ double bounded_error_Hausdorff_impl( if (local_bounds.first > global_bounds.first) { global_bounds.first = local_bounds.first; } - global_bounds.second = std::max(current_max, local_bounds.second); + global_bounds.second = std::max( + std::max(current_max, local_bounds.second), + global_bounds.first + ); // TODO Additionally store the face descriptor of the parent from TM1 in the Candidate_triangle. // Add the subtriangle to the candidate list From 30b6c70b3b84334fe9ec94458ef5fb88f4e667ff Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Thu, 18 Jul 2019 09:56:13 +0200 Subject: [PATCH 26/48] Add naive and optimized Hausdorff implementation to test file. --- .../CGAL/Polygon_mesh_processing/distance.h | 2 ++ .../test_pmp_distance.cpp | 23 +++++++++++++++++-- 2 files changed, 23 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 d2957a78576..f81ab7c7874 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1015,8 +1015,10 @@ double bounded_error_Hausdorff_impl( } // Print result found +/* std::cout << "Processing candidates finished, found distance (lower, upper): (" << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; +*/ // Return linear interpolation between found lower and upper bound return (global_bounds.first + global_bounds.second) / 2.; diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp index 08936586a6f..e1bc4c828dd 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp @@ -283,7 +283,7 @@ int main(int, char** argv) CGAL::Real_timer time; #if defined(CGAL_LINKED_WITH_TBB) time.start(); - std::cout << "Distance between meshes (parallel) " + std::cout << "Hausdorff distance approximation between meshes (parallel) " << PMP::approximate_Hausdorff_distance( m1,m2,PMP::parameters::number_of_points_per_area_unit(4000)) << "\n"; @@ -293,13 +293,32 @@ int main(int, char** argv) time.reset(); time.start(); - std::cout << "Distance between meshes (sequential) " + std::cout << "Hausdorff distance approximation between meshes (sequential) " << PMP::approximate_Hausdorff_distance( m1,m2,PMP::parameters::number_of_points_per_area_unit(4000)) << "\n"; time.stop(); std::cout << "done in " << time.time() << "s.\n"; + time.reset(); + time.start(); + std::cout << "Lower bound on Hausdorff distance between meshes (sequential), " + << "the actual distance is at most " << 0.01 << " larger " + << PMP::bounded_error_Hausdorff_distance_naive( + m1,m2,0.01) + << "\n"; + time.stop(); + std::cout << "done in " << time.time() << "s.\n"; + + time.reset(); + time.start(); + std::cout << "Bounded Hausdorff distance between meshes (sequential), optimized implementation " + << PMP::bounded_error_Hausdorff_distance( + m1,m2,0.01) + << "\n"; + time.stop(); + std::cout << "done in " << time.time() << "s.\n"; + general_tests(m1,m2); test_concept(); From 98e92a973bdb6771d39ad80cd0439323a6ba911c Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Thu, 18 Jul 2019 11:42:16 +0200 Subject: [PATCH 27/48] Introduce an additional stopping criterion to stop in degenerate cases where a triangle is never projected onto a single triangle in the other mesh. --- .../CGAL/Polygon_mesh_processing/distance.h | 29 +++++++++++++++++++ ...traversal_traits_with_Hausdorff_distance.h | 4 +++ .../test_pmp_distance.cpp | 2 +- 3 files changed, 34 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 f81ab7c7874..2b290837bed 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -927,8 +927,22 @@ double bounded_error_Hausdorff_impl( Candidate_set candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); +/* + std::cout << "Found " << candidate_triangles.size() << " candidates." << std::endl; + for (int i=0; i error_bound) && candidate_triangles.size() > 0 ) { + // std::cout << "Current number candidates: " << candidate_triangles.size() << std::endl; + // std::cout << "Current global bounds: (" << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; + // Get the first triangle and its Hausdorff bounds from the candidate set Candidate_triangle triangle_and_bound = candidate_triangles.back(); // Remove it from the candidate set as it will be processed now @@ -939,6 +953,9 @@ double bounded_error_Hausdorff_impl( // and the difference between the bounds to be obtained is larger than the // user given error. Hausdorff_bounds triangle_bounds = triangle_and_bound.second; + + // std::cout << "Current triangle bounds: (" << triangle_bounds.first << ", " << triangle_bounds.second << ")" << std::endl; + if ( (triangle_bounds.second > global_bounds.first) && (triangle_bounds.second - triangle_bounds.first > error_bound) ) { // Get the triangle that is to be subdivided and read its vertices Triangle_3 triangle_for_subdivision = triangle_and_bound.first; @@ -959,6 +976,18 @@ double bounded_error_Hausdorff_impl( continue; } + // Check third stopping condition: All edge lengths of the triangle are + // smaller than the given error bound, cannot get results beyond this + // bound. + if ( squared_distance( v0, v1 ) < squared_error_bound + && squared_distance( v0, v2 ) < squared_error_bound + && squared_distance( v1, v2 ) < squared_error_bound ) { + // The upper bound of this triangle is within error tolerance of + // the actual upper bound, use it. + global_bounds.first = triangle_bounds.second; + continue; + } + // Subdivide the triangle into four smaller triangles Point_3 v01 = midpoint( v0, v1 ); Point_3 v02 = midpoint( v0, v2 ); 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 6393c74a912..ee93fff07eb 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 @@ -230,6 +230,10 @@ namespace CGAL { // Update global Hausdorff bounds according to the obtained local bounds Hausdorff_bounds local_bounds = traversal_traits_tm2.get_local_bounds(); +/* + std::cout << "Processed triangle " << primitive.id() << " with bounds (" + << local_bounds.first << ", " << local_bounds.second << ")" << std::endl; +*/ if (local_bounds.first > h_lower) { h_lower = local_bounds.first; } diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp index e1bc4c828dd..33b17e54de2 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp @@ -303,7 +303,7 @@ int main(int, char** argv) time.reset(); time.start(); std::cout << "Lower bound on Hausdorff distance between meshes (sequential), " - << "the actual distance is at most " << 0.01 << " larger " + << "the actual distance is at most " << 0.01 << " larger than " << PMP::bounded_error_Hausdorff_distance_naive( m1,m2,0.01) << "\n"; From 130dcbc1580a33e31fd8559b2c97e5cbe501075a Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 19 Jul 2019 11:07:06 +0200 Subject: [PATCH 28/48] Bugfix entering nodes of TM2 when traversing it. --- ...BB_traversal_traits_with_Hausdorff_distance.h | 16 ++++++++++------ .../test_pmp_distance.cpp | 4 ++-- 2 files changed, 12 insertions(+), 8 deletions(-) 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 ee93fff07eb..6cf6c307146 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 @@ -146,8 +146,15 @@ namespace CGAL { squared_distance(tri_center, v2) )); - // If triangles in the node can still lower the upper bound, enter it - if (tri_radius + bb_radius >= dist) { + // Find a lower bound on the distance of the triangle to the bounding box + // by taking the distance from the triangle center to the bbox center + // and subtracting the radii of both elements. Thereby, each triangle in + // the bounding box has at least distance + // ( dist - triangle_radius - bbox_radius ) + // to the query triangle. Only if this lower bound is lower than the + // current best known lower bound, we enter the node, trying to + // lower the lower bound further. + if ( (dist - tri_radius - bb_radius) <= h_local_lower) { return true; } else { return false; @@ -230,10 +237,7 @@ namespace CGAL { // Update global Hausdorff bounds according to the obtained local bounds Hausdorff_bounds local_bounds = traversal_traits_tm2.get_local_bounds(); -/* - std::cout << "Processed triangle " << primitive.id() << " with bounds (" - << local_bounds.first << ", " << local_bounds.second << ")" << std::endl; -*/ + if (local_bounds.first > h_lower) { h_lower = local_bounds.first; } diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp index 33b17e54de2..6a6ac16db0a 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp @@ -290,7 +290,7 @@ int main(int, char** argv) time.stop(); std::cout << "done in " << time.time() << "s.\n"; #endif - +/* time.reset(); time.start(); std::cout << "Hausdorff distance approximation between meshes (sequential) " @@ -309,7 +309,7 @@ int main(int, char** argv) << "\n"; time.stop(); std::cout << "done in " << time.time() << "s.\n"; - +*/ time.reset(); time.start(); std::cout << "Bounded Hausdorff distance between meshes (sequential), optimized implementation " From c91c780cae4d1e26e178578a5c8cafeb884ab91c Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 19 Jul 2019 11:27:59 +0200 Subject: [PATCH 29/48] Add timing to the example file and restore all tests in the test file. --- ...hausdorff_bounded_error_distance_example.cpp | 17 ++++++++++++----- .../test_pmp_distance.cpp | 4 ++-- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index 0242dcfe4fd..01ec8a43329 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #if defined(CGAL_LINKED_WITH_TBB) #define TAG CGAL::Parallel_tag @@ -61,28 +62,34 @@ int main(int argc, char** argv) std::ifstream input( argv[1] ); input >> tm1; std::cout << "Read a mesh with " << tm1.number_of_faces() << " triangles." << std::endl; +/* std::ifstream input2( argv[2] ); input2 >> tm2; std::cout << "Read a mesh with " << tm2.number_of_faces() << " triangles." << std::endl; +*/ // Copy the mesh and perturb it slightly -/* tm2 = tm1; bool do_project = false; - CGAL::Polygon_mesh_processing::random_perturbation( tm2.vertices(), tm2, 0.001, do_project ); + CGAL::Polygon_mesh_processing::random_perturbation( tm2.vertices(), tm2, 0.1, do_project ); std::cout << "Perturbed the input mesh, now computing the Hausdorff distance." << std::endl; -*/ - -// https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__meshing__grp.html#ga028a80dc84395650f67714fa7618ec53 + // Set an error bound double error_bound = 0.01; + CGAL::Real_timer time; + + time.start(); std::cout << "Approximated Hausdorff distance: " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance (tm1, tm2, error_bound) << std::endl; + time.stop(); + std::cout << "Processing took " << time.time() << "s." << std::endl; +/* std::cout << "Approximated Hausdorff distance (naive): " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance_naive (tm1, tm2, error_bound) << ", the actual distance is at most " << error_bound << " larger." << std::endl; +*/ } diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp index 6a6ac16db0a..33b17e54de2 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp @@ -290,7 +290,7 @@ int main(int, char** argv) time.stop(); std::cout << "done in " << time.time() << "s.\n"; #endif -/* + time.reset(); time.start(); std::cout << "Hausdorff distance approximation between meshes (sequential) " @@ -309,7 +309,7 @@ int main(int, char** argv) << "\n"; time.stop(); std::cout << "done in " << time.time() << "s.\n"; -*/ + time.reset(); time.start(); std::cout << "Bounded Hausdorff distance between meshes (sequential), optimized implementation " From 156cac510777f3bb9889d5f3d3c8a5e16eb16000 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sat, 20 Jul 2019 19:05:55 +0200 Subject: [PATCH 30/48] Implement benchmarking on both time and number of culled triangles. --- ...usdorff_bounded_error_distance_example.cpp | 63 ++++++++++++++++--- .../CGAL/Polygon_mesh_processing/distance.h | 2 + ...traversal_traits_with_Hausdorff_distance.h | 13 ++++ 3 files changed, 71 insertions(+), 7 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index 01ec8a43329..f4146d98877 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -21,7 +21,15 @@ namespace PMP = CGAL::Polygon_mesh_processing; int main(int argc, char** argv) { + // Objects to hold the meshes Mesh tm1, tm2; + + // Timer to measure the runtime of h-distance Distance_computation + CGAL::Real_timer time; + + // Set an error bound + double error_bound = 0.01; + /* // Easy Example CGAL::make_tetrahedron(Point_3(.0,.0,.0), @@ -32,7 +40,7 @@ int main(int argc, char** argv) tm2=tm1; CGAL::Polygon_mesh_processing::isotropic_remeshing(tm2.faces(),.05, tm2); */ - +// ---------------------------------------------------------------------------- /* // Example with point realizing the Hausdorff distance strictly lying in the // interior of a triangle @@ -57,27 +65,67 @@ int main(int argc, char** argv) std::cout << "TM2 is valid: " << (tm2.is_valid() ? "true" : "false") << std::endl; */ - +// ---------------------------------------------------------------------------- +/* // Read a real mesh given by the user std::ifstream input( argv[1] ); input >> tm1; std::cout << "Read a mesh with " << tm1.number_of_faces() << " triangles." << std::endl; -/* + std::ifstream input2( argv[2] ); input2 >> tm2; std::cout << "Read a mesh with " << tm2.number_of_faces() << " triangles." << std::endl; -*/ + // Copy the mesh and perturb it slightly tm2 = tm1; bool do_project = false; CGAL::Polygon_mesh_processing::random_perturbation( tm2.vertices(), tm2, 0.1, do_project ); std::cout << "Perturbed the input mesh, now computing the Hausdorff distance." << std::endl; +*/ +// ---------------------------------------------------------------------------- - // Set an error bound - double error_bound = 0.01; + // Pairwise computation on a set of benchmark models + std::vector models = {"80","102","128","162","204","256","320","402","504","632","792","992","1242"}; + int num_models = models.size(); + std::string prefix = "/home/martin/Downloads/bunnies/bunny_"; + std::string postfix = ".off"; + std::vector error_bounds = { 0.1, 0.01, 0.001, 0.0001 }; - CGAL::Real_timer time; + std::ifstream input(prefix + models[0] + postfix); + // input >> tm1; + // input >> tm2; + // std::cout << "Initialized with a mesh at " << (prefix + models[0] + postfix) << " with " << tm1.number_of_faces() << " triangles." << std::endl; + for(int i=0; i> tm1; + // std::cout << "Read a mesh at " << (prefix + models[i] + postfix) << " with " << tm1.number_of_faces() << " faces." << std::endl; + + for(int j=0; j> tm2; + // std::cout << "Read a mesh at " << (prefix + models[j] + postfix) << " with " << tm2.number_of_faces() << " faces." << std::endl; + // + // std::cout << "Read two meshes with " << tm1.number_of_faces() << ", " << tm2.number_of_faces() << " triangles respectively." << std::endl; + + for (int k=0; k(tm1, tm2, error_bounds[k]); + time.stop(); + std::cout << models[i] << " " << models[j] << " " << time.time() << " " << error_bounds[k] << " " << h_dist << std::endl; + } + } + } + + +/* time.start(); std::cout << "Approximated Hausdorff distance: " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance @@ -85,6 +133,7 @@ int main(int argc, char** argv) << std::endl; time.stop(); std::cout << "Processing took " << time.time() << "s." << std::endl; +*/ /* std::cout << "Approximated Hausdorff distance (naive): " << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance_naive 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 2b290837bed..2123aada381 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -927,6 +927,8 @@ double bounded_error_Hausdorff_impl( Candidate_set candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); + // std::cout << "Culled " << traversal_traits_tm1.get_num_culled_triangles() << " out of " << tm1.num_faces() << std::endl; + /* std::cout << "Found " << candidate_triangles.size() << " candidates." << std::endl; for (int i=0; i m_face_to_triangle_map( &m_tm1, m_vpm1 ); Triangle_3 candidate_triangle = get(m_face_to_triangle_map, primitive.id()); // Call Culling on B with the single triangle found. + // TODO If we project the vertices of the triangle onto B and consider + // the edge length of the triangle, can this give us better bounds + // to initialize with here? Hausdorff_primitive_traits_tm2< Tree_traits, Triangle_3, Kernel, TriangleMesh, VPM2 > traversal_traits_tm2( @@ -298,6 +304,11 @@ namespace CGAL { return Hausdorff_bounds( h_lower, h_upper ); } + int get_num_culled_triangles() + { + return (m_tm1.num_faces() - m_investigated_on_tm1); + } + private: const AABBTraits& m_traits; // The two triangle meshes @@ -313,6 +324,8 @@ namespace CGAL { double h_upper; // List of candidate triangles with their Hausdorff bounds attached Candidate_set m_candidiate_triangles; + // Number of triangles investigated in the procedure + int m_investigated_on_tm1; }; } From 680edc46b1d7de2811e1e81f12fc8d99ff7c041e Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sun, 21 Jul 2019 11:07:21 +0200 Subject: [PATCH 31/48] Modify triangle bounding box distance to have sharper bounds. --- ...traversal_traits_with_Hausdorff_distance.h | 95 ++++++++++++++++++- 1 file changed, 91 insertions(+), 4 deletions(-) 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 2606c1bd5b3..9cd60db374a 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 @@ -40,6 +40,7 @@ namespace CGAL { typedef typename AABBTraits::Bounding_box Bounding_box; typename Kernel::Construct_projected_point_3 project_point; typedef typename Kernel::Point_3 Point_3; + typedef typename Kernel::Vector_3 Vector_3; public: Hausdorff_primitive_traits_tm2( @@ -112,6 +113,7 @@ namespace CGAL { { /* Have reached a node, determine whether or not to enter it */ + /* // Get the bounding box of the nodes Bounding_box bbox = node.bbox(); // Compute its center @@ -159,6 +161,60 @@ namespace CGAL { } else { return false; } + */ + + // Get the bounding box of the nodes + Bounding_box bbox = node.bbox(); + // Get the vertices of the query triangle + Point_3 v0 = query.vertex(0); + Point_3 v1 = query.vertex(1); + Point_3 v2 = query.vertex(2); + // Find the axis aligned bbox of the triangle + Point_3 tri_min = Point_3 ( + std::min(std::min( v0.x(), v1.x()), v2.x() ), + std::min(std::min( v0.y(), v1.y()), v2.y() ), + std::min(std::min( v0.z(), v1.z()), v2.z() ) + ); + Point_3 tri_max = Point_3 ( + std::max(std::max( v0.x(), v1.x()), v2.x() ), + std::max(std::max( v0.y(), v1.y()), v2.y() ), + std::max(std::max( v0.z(), v1.z()), v2.z() ) + ); + + // Compute distance of the bounding boxes + // Distance along the x-axis + double dist_x = 0.; + if ( tri_max.x() < bbox.min(0) ) { + dist_x = bbox.min(0) - tri_max.x(); + } else if ( bbox.max(0) < tri_min.x() ) { + dist_x = tri_min.x() - bbox.max(0); + } + // Distance along the y-axis + double dist_y = 0.; + if ( tri_max.y() < bbox.min(1) ) { + dist_y = bbox.min(1) - tri_max.y(); + } else if ( bbox.max(1) < tri_min.y() ) { + dist_y = tri_min.y() - bbox.max(1); + } + // Distance along the y-axis + double dist_z = 0.; + if ( tri_max.z() < bbox.min(2) ) { + dist_z = bbox.min(2) - tri_max.z(); + } else if ( bbox.max(2) < tri_min.z() ) { + dist_z = tri_min.z() - bbox.max(2); + } + + // Lower bound on the distance between the two bounding boxes is given + // as the length of the diagonal of the bounding box between them + double dist = approximate_sqrt( Vector_3(dist_x,dist_y,dist_z).squared_length() ); + + // Check whether investigating the bbox can still lower the Hausdorff + // distance. If so, enter the box. + if ( dist <= h_local_lower) { + return true; + } else { + return false; + } } // Return the local Hausdorff bounds computed for the passed query triangle @@ -194,6 +250,7 @@ namespace CGAL { typedef typename AABBTraits::Bounding_box Bounding_box; typedef ::CGAL::AABB_node Node; typedef typename Kernel::Point_3 Point_3; + typedef typename Kernel::Vector_3 Vector_3; typedef typename Kernel::Triangle_3 Triangle_3; typedef std::pair Candidate_triangle; typedef typename std::vector Candidate_set; @@ -263,10 +320,7 @@ namespace CGAL { bool do_intersect(const Query& /*query*/, const Node& node) const { /* Have reached a node, determine whether or not to enter it */ - - // TODO What's the closest distance of TM2 to the box given by node? - // Can we have a sharper bound on this than the one implemented below? - +/* // Get the bounding box of the nodes Bounding_box bbox = node.bbox(); // Compute its center @@ -290,6 +344,39 @@ namespace CGAL { } else { return false; } +*/ + // Get the bounding box of the nodes + Bounding_box bbox = node.bbox(); + // Compute its center + Point_3 center = Point_3( + (bbox.min(0) + bbox.max(0)) / 2, + (bbox.min(1) + bbox.max(1)) / 2, + (bbox.min(2) + bbox.max(2)) / 2); + // Find the point from TM2 closest to the center + // TODO Insert a hint here to accelerate the query + Point_3 closest = m_tm2_tree.closest_point(center); + // Compute the difference vector between the bbox center and the closest + // point in tm2 + Vector_3 difference = Vector_3( closest, center ); + // Shift the vector to be the difference between the farthest corner + // of the bounding box away from the closest point on TM2 + double diff_x = (bbox.max(0) - bbox.min(0)) / 2; + if (difference.x() < 0) diff_x = diff_x * -1.; + double diff_y = (bbox.max(1) - bbox.min(1)) / 2; + if (difference.y() < 0) diff_y = diff_y * -1.; + double diff_z = (bbox.max(2) - bbox.min(2)) / 2; + if (difference.z() < 0) diff_z = diff_z * -1.; + difference = difference + Vector_3( diff_x, diff_y, diff_z ); + // Compute distance from the farthest corner of the bbox to the closest + // point in TM2 + double dist = approximate_sqrt( difference.squared_length() ); + + // If the distance is larger than the global lower bound, enter the node, i.e. return true. + if (dist > h_lower) { + return true; + } else { + return false; + } } // Return those triangles from TM1 which are candidates for including a From 405cf122493c4e9c98a9164d8424a1ba287df0b2 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sun, 21 Jul 2019 14:18:35 +0200 Subject: [PATCH 32/48] Implement Priority Queue for the set of candidate triangles. --- .../CGAL/Polygon_mesh_processing/distance.h | 44 +++++++++---------- ...traversal_traits_with_Hausdorff_distance.h | 41 ++++++++++++++--- .../test_pmp_distance.cpp | 8 ++-- 3 files changed, 63 insertions(+), 30 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 2123aada381..3293fbbc116 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -898,13 +898,19 @@ double bounded_error_Hausdorff_impl( typedef typename boost::property_map::const_type Triangle_hausdorff_bounds; typedef std::pair Hausdorff_bounds; - typedef std::pair Candidate_triangle; - typedef typename std::vector Candidate_set; typename Kernel::Compute_squared_distance_3 squared_distance; typename Kernel::Construct_projected_point_3 project_point; typename Kernel::FT dist; + typedef + #if BOOST_VERSION >= 105000 + boost::heap::priority_queue< Candidate_triangle, boost::heap::compare< std::greater > > > + #else + std::priority_queue< Candidate_triangle > + #endif + Heap_type; + // Build an AABB tree on tm1 TM1_tree tm1_tree( faces(tm1).begin(), faces(tm1).end(), tm1, vpm1 ); tm1_tree.build(); @@ -924,7 +930,7 @@ double bounded_error_Hausdorff_impl( // TODO Implement the candidate_triangles set as Stack instead of Vector // check: https://www.boost.org/doc/libs/1_55_0/doc/html/heap.html // Can already build a sorted structure while collecting the candidates - Candidate_set candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); + Heap_type candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); // std::cout << "Culled " << traversal_traits_tm1.get_num_culled_triangles() << " out of " << tm1.num_faces() << std::endl; @@ -940,27 +946,27 @@ double bounded_error_Hausdorff_impl( double squared_error_bound = error_bound * error_bound; - while ( (global_bounds.second - global_bounds.first > error_bound) && candidate_triangles.size() > 0 ) { + while ( (global_bounds.second - global_bounds.first > error_bound) && !candidate_triangles.empty() ) { // std::cout << "Current number candidates: " << candidate_triangles.size() << std::endl; // std::cout << "Current global bounds: (" << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; // Get the first triangle and its Hausdorff bounds from the candidate set - Candidate_triangle triangle_and_bound = candidate_triangles.back(); + Candidate_triangle triangle_and_bound = candidate_triangles.top(); // Remove it from the candidate set as it will be processed now - candidate_triangles.pop_back(); + candidate_triangles.pop(); // Only process the triangle if it can contribute to the Hausdorff distance, // i.e. if its Upper Bound is higher than the currently known best lower bound // and the difference between the bounds to be obtained is larger than the // user given error. - Hausdorff_bounds triangle_bounds = triangle_and_bound.second; + Hausdorff_bounds triangle_bounds = triangle_and_bound.m_bounds; // std::cout << "Current triangle bounds: (" << triangle_bounds.first << ", " << triangle_bounds.second << ")" << std::endl; if ( (triangle_bounds.second > global_bounds.first) && (triangle_bounds.second - triangle_bounds.first > error_bound) ) { // Get the triangle that is to be subdivided and read its vertices - Triangle_3 triangle_for_subdivision = triangle_and_bound.first; + Triangle_3 triangle_for_subdivision = triangle_and_bound.m_triangle; Point_3 v0 = triangle_for_subdivision.vertex(0); Point_3 v1 = triangle_for_subdivision.vertex(1); Point_3 v2 = triangle_for_subdivision.vertex(2); @@ -1012,27 +1018,17 @@ double bounded_error_Hausdorff_impl( ); tm2_tree.traversal(sub_triangles[i], traversal_traits_tm2); - // Get the highest current bound from all candidate triangles - double current_max = 0.; - for(auto&& ct: candidate_triangles) { - if (ct.second.second > current_max) { - current_max = ct.second.second; - } - } - - // Update global Hausdorff bounds according to the obtained local bounds + // Update global lower Hausdorff bound according to the obtained local bounds Hausdorff_bounds local_bounds = traversal_traits_tm2.get_local_bounds(); if (local_bounds.first > global_bounds.first) { global_bounds.first = local_bounds.first; } - global_bounds.second = std::max( - std::max(current_max, local_bounds.second), - global_bounds.first - ); // TODO Additionally store the face descriptor of the parent from TM1 in the Candidate_triangle. // Add the subtriangle to the candidate list - candidate_triangles.push_back(Candidate_triangle(sub_triangles[i], local_bounds)); + candidate_triangles.push( + Candidate_triangle(sub_triangles[i], local_bounds) + ); // std::cout << "Split triangle (" << v0 << ", " << v1 << ", " << v2 // << ") with bounds: (" << triangle_bounds.first << ", " @@ -1042,6 +1038,10 @@ double bounded_error_Hausdorff_impl( // << local_bounds.first << ", " << local_bounds.second << "), gobal bounds are: (" // << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; } + + // Update global upper Hausdorff bound after subdivision + double current_max = candidate_triangles.top().m_bounds.second; + global_bounds.second = std::max(current_max, global_bounds.first); } } 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 9cd60db374a..ccafde7ecc8 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 @@ -29,6 +29,28 @@ namespace CGAL { typedef std::pair Hausdorff_bounds; + /** + * @struct Candidate_triangle + */ + template + struct Candidate_triangle { + typedef typename Kernel::Triangle_3 Triangle_3; + + Candidate_triangle(const Triangle_3& triangle, const Hausdorff_bounds& bounds) + : m_triangle(triangle), m_bounds(bounds) {} + + Triangle_3 m_triangle; + Hausdorff_bounds m_bounds; + + #if BOOST_VERSION >= 105000 + bool operator<(const Candidate_triangle& other) const { return m_bounds.second < other.m_bounds.second; } + bool operator>(const Candidate_triangle& other) const { return m_bounds.second > other.m_bounds.second; } + #else + bool operator>(const Candidate_triangle& other) const { return m_bounds.second < other.m_bounds.second; } + bool operator<(const Candidate_triangle& other) const { return m_bounds.second > other.m_bounds.second; } + #endif + }; + /** * @class Hausdorff_primitive_traits_tm2 */ @@ -252,11 +274,17 @@ namespace CGAL { typedef typename Kernel::Point_3 Point_3; typedef typename Kernel::Vector_3 Vector_3; typedef typename Kernel::Triangle_3 Triangle_3; - typedef std::pair Candidate_triangle; - typedef typename std::vector Candidate_set; + typedef typename std::vector> Candidate_set; typedef AABB_tree< AABB_traits > TM2_tree; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef + #if BOOST_VERSION >= 105000 + boost::heap::priority_queue< Candidate_triangle, boost::heap::compare< std::greater > > > + #else + std::priority_queue< Candidate_triangle > + #endif + Heap_type; public: Hausdorff_primitive_traits_tm1(const AABBTraits& traits, const TM2_tree& tree, const TriangleMesh& tm1, const TriangleMesh& tm2 , const VPM1& vpm1, const VPM2& vpm2 ) @@ -311,7 +339,8 @@ namespace CGAL { // Store the triangle given as primitive here as candidate triangle // together with the local bounds it obtained to sind it to subdivision // later - m_candidiate_triangles.push_back(Candidate_triangle(candidate_triangle, local_bounds)); + m_candidiate_triangles.push( Candidate_triangle(candidate_triangle, local_bounds) ); + pq.push( Candidate_triangle(candidate_triangle, local_bounds) ); } // Determine whether child nodes will still contribute to a larger @@ -381,7 +410,7 @@ namespace CGAL { // Return those triangles from TM1 which are candidates for including a // point realizing the Hausdorff distance - Candidate_set get_candidate_triangles() { + Heap_type get_candidate_triangles() { return m_candidiate_triangles; } @@ -410,7 +439,9 @@ namespace CGAL { double h_lower; double h_upper; // List of candidate triangles with their Hausdorff bounds attached - Candidate_set m_candidiate_triangles; + Heap_type m_candidiate_triangles; + // Heap of candidate triangles with their Hausdorff bounds attached + Heap_type pq; // Number of triangles investigated in the procedure int m_investigated_on_tm1; }; diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp index 33b17e54de2..6a499023c7b 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp @@ -300,12 +300,14 @@ int main(int, char** argv) time.stop(); std::cout << "done in " << time.time() << "s.\n"; + double error_bound = 0.001; + time.reset(); time.start(); std::cout << "Lower bound on Hausdorff distance between meshes (sequential), " - << "the actual distance is at most " << 0.01 << " larger than " + << "the actual distance is at most " << error_bound << " larger than " << PMP::bounded_error_Hausdorff_distance_naive( - m1,m2,0.01) + m1,m2,error_bound) << "\n"; time.stop(); std::cout << "done in " << time.time() << "s.\n"; @@ -314,7 +316,7 @@ int main(int, char** argv) time.start(); std::cout << "Bounded Hausdorff distance between meshes (sequential), optimized implementation " << PMP::bounded_error_Hausdorff_distance( - m1,m2,0.01) + m1,m2,error_bound) << "\n"; time.stop(); std::cout << "done in " << time.time() << "s.\n"; From e7e724e4f9eef0e0d2adac11b575a7a7a41b96bd Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 23 Aug 2019 07:44:08 +0200 Subject: [PATCH 33/48] Include Benchmark examples in the Hausdorff examples file. --- ...usdorff_bounded_error_distance_example.cpp | 157 +++++++++--------- 1 file changed, 82 insertions(+), 75 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index f4146d98877..c0a18953ac2 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -1,8 +1,13 @@ +#include +#include +#include #include #include +#include #include #include #include +#include #include #if defined(CGAL_LINKED_WITH_TBB) @@ -14,6 +19,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef K::Point_3 Point_3; typedef K::Triangle_3 Triangle_3; +typedef K::Vector_3 Vector_3; typedef CGAL::Surface_mesh Mesh; typedef Mesh::Vertex_index Vertex_index; @@ -28,21 +34,35 @@ int main(int argc, char** argv) CGAL::Real_timer time; // Set an error bound - double error_bound = 0.01; + double error_bound = 0.0001; -/* - // Easy Example +// ---------------------------------------------------------------------------- + + // Easy Example of a tetrahedron and a remeshed version of itself + + // Create the Tetrahedron CGAL::make_tetrahedron(Point_3(.0,.0,.0), Point_3(2,.0,.0), Point_3(1,1,1), Point_3(1,.0,2), tm1); + // Copy it and remesh it tm2=tm1; - CGAL::Polygon_mesh_processing::isotropic_remeshing(tm2.faces(),.05, tm2); -*/ + PMP::isotropic_remeshing(tm2.faces(),.05, tm2); + // Compute the Hausdorff distance + time.reset(); + time.start(); + std::cout << "Approximated Hausdorff distance: " + << PMP::bounded_error_Hausdorff_distance + (tm1, tm2, error_bound) + << std::endl; + time.stop(); + std::cout << "Processing took " << time.time() << "s." << std::endl; + // ---------------------------------------------------------------------------- -/* + // Example with point realizing the Hausdorff distance strictly lying in the + // interior of a triangle tm1 = Mesh(); tm1.add_vertex( Point_3(-1.,1.,1.) ); @@ -50,8 +70,6 @@ int main(int argc, char** argv) tm1.add_vertex( Point_3(1.,1.,1.) ); tm1.add_face( tm1.vertices() ); - std::cout << "TM1 is valid: " << (tm1.is_valid() ? "true" : "false") << std::endl; - tm2 = Mesh(); Vertex_index w0 = tm2.add_vertex( Point_3(-1.,1.,0.) ); Vertex_index w1 = tm2.add_vertex( Point_3(0.,-1.,0.) ); @@ -63,82 +81,71 @@ int main(int argc, char** argv) tm2.add_face( w1, w4, w5 ); tm2.add_face( w2, w5, w3 ); - std::cout << "TM2 is valid: " << (tm2.is_valid() ? "true" : "false") << std::endl; -*/ + // Compute the Hausdorff distance + time.reset(); + time.start(); + std::cout << "Approximated Hausdorff distance: " + << PMP::bounded_error_Hausdorff_distance + (tm1, tm2, error_bound) + << std::endl; + time.stop(); + std::cout << "Processing took " << time.time() << "s." << std::endl; + // ---------------------------------------------------------------------------- -/* - // Read a real mesh given by the user + + // Read a real meshes given by the user, perturb it slightly and compute the + // Hausdorff distance between the original mesh and its pertubation + std::ifstream input( argv[1] ); input >> tm1; std::cout << "Read a mesh with " << tm1.number_of_faces() << " triangles." << std::endl; + // Copy the mesh and perturb it slightly + tm2 = tm1; + bool do_project = false; + PMP::random_perturbation( tm2.vertices(), tm2, 0.1, do_project ); + std::cout << "Perturbed the input mesh, now computing the Hausdorff distance." << std::endl; + + // Compute the Hausdorff distance + time.reset(); + time.start(); + std::cout << "Approximated Hausdorff distance: " + << PMP::bounded_error_Hausdorff_distance + (tm1, tm2, error_bound) + << std::endl; + time.stop(); + std::cout << "Processing took " << time.time() << "s." << std::endl; + +// ---------------------------------------------------------------------------- + + // Read two meshes given by the user, initially place them at their originally + // given position. Move the second mesh in 300 steps away from the first one. + // Print how the Hausdorff distance changes. + + std::ifstream input1( argv[1] ); + input1 >> tm1; + std::cout << "Read a mesh with " << tm1.number_of_faces() << " triangles." << std::endl; + std::ifstream input2( argv[2] ); input2 >> tm2; std::cout << "Read a mesh with " << tm2.number_of_faces() << " triangles." << std::endl; - // Copy the mesh and perturb it slightly - tm2 = tm1; - bool do_project = false; - CGAL::Polygon_mesh_processing::random_perturbation( tm2.vertices(), tm2, 0.1, do_project ); - std::cout << "Perturbed the input mesh, now computing the Hausdorff distance." << std::endl; -*/ -// ---------------------------------------------------------------------------- + CGAL::Bbox_3 bb = PMP::bbox(tm2); + double dist = CGAL::approximate_sqrt( Vector_3(bb.xmax() - bb.xmin(), bb.ymax() - bb.ymin(), bb.zmax() - bb.zmin()).squared_length() ); - // Pairwise computation on a set of benchmark models - std::vector models = {"80","102","128","162","204","256","320","402","504","632","792","992","1242"}; - int num_models = models.size(); - std::string prefix = "/home/martin/Downloads/bunnies/bunny_"; - std::string postfix = ".off"; - std::vector error_bounds = { 0.1, 0.01, 0.001, 0.0001 }; - - std::ifstream input(prefix + models[0] + postfix); - // input >> tm1; - // input >> tm2; - // std::cout << "Initialized with a mesh at " << (prefix + models[0] + postfix) << " with " << tm1.number_of_faces() << " triangles." << std::endl; - - for(int i=0; i> tm1; - // std::cout << "Read a mesh at " << (prefix + models[i] + postfix) << " with " << tm1.number_of_faces() << " faces." << std::endl; - - for(int j=0; j> tm2; - // std::cout << "Read a mesh at " << (prefix + models[j] + postfix) << " with " << tm2.number_of_faces() << " faces." << std::endl; - // - // std::cout << "Read two meshes with " << tm1.number_of_faces() << ", " << tm2.number_of_faces() << " triangles respectively." << std::endl; - - for (int k=0; k(tm1, tm2, error_bounds[k]); - time.stop(); - std::cout << models[i] << " " << models[j] << " " << time.time() << " " << error_bounds[k] << " " << h_dist << std::endl; - } - } + for (int i=0; i<300; i++) { + PMP::transform( + CGAL::Aff_transformation_3 ( CGAL::Translation(), Vector_3( 0.01*dist, 0.01*dist, 0.01*dist ) ), + tm2 + ); + time.reset(); + time.start(); + std::cout << "Position: " << i << std::endl; + std::cout << "Approximated Hausdorff distance: " + << PMP::bounded_error_Hausdorff_distance + (tm1, tm2, error_bound) + << std::endl; + time.stop(); + std::cout << "Processing took " << time.time() << "s." << std::endl; } - - -/* - time.start(); - std::cout << "Approximated Hausdorff distance: " - << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance - (tm1, tm2, error_bound) - << std::endl; - time.stop(); - std::cout << "Processing took " << time.time() << "s." << std::endl; -*/ -/* - std::cout << "Approximated Hausdorff distance (naive): " - << CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance_naive - (tm1, tm2, error_bound) - << ", the actual distance is at most " << error_bound << " larger." - << std::endl; -*/ } From 1e37c32d6e1dc29398403f6412f1466699928cf7 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 23 Aug 2019 07:46:52 +0200 Subject: [PATCH 34/48] Write first lines of documentation to test linking and inclusion of example. --- .../Polygon_mesh_processing.txt | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index 004aad17948..47f486bbc62 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt @@ -729,10 +729,12 @@ which enables to treat one or several connected components as a face graph. \cgalExample{Polygon_mesh_processing/face_filtered_graph_example.cpp} -\section PMPDistance Approximate Hausdorff Distance +\section PMPDistance (Approximate) Hausdorff Distance This package provides methods to compute (approximate) distances between meshes and point sets. +\subsection Approximate Hausdorff Distance + The function \link approximate_Hausdorff_distance() `approximate_Hausdorff_distance()`\endlink computes an approximation of the Hausdorff distance from a mesh `tm1` to a mesh `tm2`. Given a a sampling of `tm1`, it computes the distance to `tm2` of the farthest sample point to `tm2` \cgalCite{cignoni1998metro}. @@ -757,12 +759,12 @@ computes an approximation of the Hausdorff distance from a mesh to a point set. For each triangle, a lower and upper bound of the Hausdorff distance to the point set are computed. Triangles are refined until the difference between the bounds is lower than a user-defined precision threshold. -\subsection AHDExample Approximate Hausdorff Distance Example +\subsubsection AHDExample Approximate Hausdorff Distance Example In the following example, a mesh is isotropically remeshed and the approximate distance between the input and the output is computed. \cgalExample{Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp} -\subsection PoissonDistanceExample Max Distance Between Point Set and Surface Example +\subsubsection PoissonDistanceExample Max Distance Between Point Set and Surface Example In \ref Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp, a triangulated surface mesh is constructed from a point set using the \link PkgPoissonSurfaceReconstruction3 Poisson reconstruction algorithm \endlink, @@ -771,6 +773,14 @@ with the following code: \snippet Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp PMP_distance_snippet +\subsection Bounded Hausdorff Distance + +The function \link bounded_error_Hausdorff_distance() `bounded_error_Hausdorff_distance()`\endlink +computes an estimate of the Hausdorff distance of two triangle meshes which is +bounded by a user-given error bound. + +\cgalExample{POlygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp} + \section PMPDetectFeatures Feature Detection This package provides methods to detect some features of a polygon mesh. From ba0a67190f02f55209c59a7120aa847ecaafab44 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 23 Aug 2019 07:47:46 +0200 Subject: [PATCH 35/48] Clean up Code. --- .../CGAL/Polygon_mesh_processing/distance.h | 238 +----------------- ...traversal_traits_with_Hausdorff_distance.h | 93 +------ 2 files changed, 15 insertions(+), 316 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 3293fbbc116..1bd5f39fe85 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -816,49 +816,7 @@ double approximate_symmetric_Hausdorff_distance(const TriangleMesh& tm1, //////////////////////////////////////////////////////////////////////// namespace internal { -/* -#if defined(CGAL_LINKED_WITH_TBB) -template -struct Distance_computation{ - const AABB_tree& tree; - const std::vector& sample_points; - Point_3 initial_hint; - tbb::atomic* distance; - Distance_computation( - const AABB_tree& tree, - const Point_3& p, - const std::vector& sample_points, - tbb::atomic* d) - : tree(tree) - , sample_points(sample_points) - , initial_hint(p) - , distance(d) - {} - - void - operator()(const tbb::blocked_range& range) const - { - Point_3 hint = initial_hint; - double hdist = 0; - for( std::size_t i = range.begin(); i != range.end(); ++i) - { - hint = tree.closest_point(sample_points[i], hint); - typename Kernel_traits::Kernel::Compute_squared_distance_3 squared_distance; - double d = to_double(CGAL::approximate_sqrt( squared_distance(hint,sample_points[i]) )); - if (d>hdist) hdist=d; - } - - // update max value stored in distance - double current_value = *distance; - while( current_value < hdist ) - { - current_value = distance->compare_and_swap(hdist, current_value); - } - } -}; -#endif -*/ template hint = tm2_tree.any_reference_point_and_id(); // Build traversal traits for tm1_tree - Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2 ); + Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2, hint.first ); // Find candidate triangles in TM1 which might realise the Hausdorff bound tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); // dummy point given as query as not needed - // TODO Implement the candidate_triangles set as Stack instead of Vector - // check: https://www.boost.org/doc/libs/1_55_0/doc/html/heap.html - // Can already build a sorted structure while collecting the candidates + // TODO Is there a better/faster data structure than the Heap used here? + // Can already build a sorted structure while collecting the candidates Heap_type candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); - // std::cout << "Culled " << traversal_traits_tm1.get_num_culled_triangles() << " out of " << tm1.num_faces() << std::endl; - -/* - std::cout << "Found " << candidate_triangles.size() << " candidates." << std::endl; - for (int i=0; i(sub_triangles[i], local_bounds) ); - - // std::cout << "Split triangle (" << v0 << ", " << v1 << ", " << v2 - // << ") with bounds: (" << triangle_bounds.first << ", " - // << triangle_bounds.second << "), sub-triangle " << i - // << " (" << sub_triangles[i].vertex(0) << ", " << sub_triangles[i].vertex(1) << ", " << sub_triangles[i].vertex(2) - // << ") has bounds: (" - // << local_bounds.first << ", " << local_bounds.second << "), gobal bounds are: (" - // << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; } // Update global upper Hausdorff bound after subdivision @@ -1045,12 +985,6 @@ double bounded_error_Hausdorff_impl( } } - // Print result found -/* - std::cout << "Processing candidates finished, found distance (lower, upper): (" - << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; -*/ - // Return linear interpolation between found lower and upper bound return (global_bounds.first + global_bounds.second) / 2.; @@ -1058,134 +992,8 @@ double bounded_error_Hausdorff_impl( CGAL_static_assertion_msg (!(boost::is_convertible::value), "Parallel_tag is enabled but TBB is unavailable."); #else - // TODO implement parallelized version of the below here. - // if (boost::is_convertible::value) - // { - // tbb::atomic distance; - // distance=0; - // Distance_computation f(tm2_tree - // , hint, sample_points, &distance); - // tbb::parallel_for(tbb::blocked_range(0, sample_points.size()), f); - // return distance; - // } - // else + // TODO implement parallelized version of the above here. #endif -/* - { - // Store all vertices of tm1 in a vector - std::vector tm1_vertices; - tm1_vertices.reserve(num_vertices(tm1)); - tm1_vertices.insert(tm1_vertices.end(),vertices(tm1).begin(),vertices(tm1).end()); - - // Sort vertices along a Hilbert curve - spatial_sort( tm1_vertices.begin(), - tm1_vertices.end(), - Search_traits_3(vpm1) ); - - // For each vertex in tm1, store the distance to the closest triangle of tm2 - Vertex_closest_triangle_map vctm = get(Vertex_property_tag(), tm1); - // For each triangle in tm1, sotre its respective local lower and upper bound - // on the Hausdorff measure - Triangle_hausdorff_bounds thb = get(Face_property_tag(), tm1); - - // For each vertex in the first mesh, find the closest triangle in the - // second mesh, store it and also store the distance to this triangle - // in a dynamic vertex property - for(vertex_descriptor vd : tm1_vertices) - { - // Get the point represented by the vertex - typename boost::property_traits::reference pt = get(vpm1, vd); - // Use the AABB tree to find the closest point and face in tm2 - hint = tm2_tree.closest_point_and_primitive(pt, hint); - // Compute the distance of the point to the closest point in tm2 - dist = squared_distance(hint.first, pt); - double d = to_double(dist); - // Store the distance and the closest triangle in the corresponding map - put(vctm, vd, std::make_pair(d, hint.second)); - } - - // Maps the faces of tm2 to actual triangles - Triangle_from_face_descriptor_map face_to_triangle_map(&tm2, vpm2); - // Initialize global bounds on the Hausdorff measure - double h_lower = 0.; - double h_upper = 0.; - // Initialize an array of candidate triangles in A to be procesed in the - // following - std::vector candidate_triangles; - - // For each triangle in the first mesh, initialize its local upper and - // lower bound and store these in a dynamic face property for furture - // reference - for(face_descriptor fd : faces(tm1)) - { - // Initialize the local bounds for the current face fd - double h_triangle_lower = 0.; - double h_triangle_upper = std::numeric_limits::infinity(); - - // Create a halfedge descriptor for the current face and store the vertices - halfedge_descriptor hd = halfedge(fd, tm1); - std::array face_vertices = {source(hd,tm1), target(hd,tm1), target(next(hd, tm1),tm1)}; - - // Get the distance and closest triangle in tm2 for each vertex of fd - std::array,3> vertex_properties = { - get(vctm, face_vertices[0]), - get(vctm, face_vertices[1]), - get(vctm, face_vertices[2])}; - - // Convert the closest faces of tm2 to triangles - std::array triangles_in_B = { - get(face_to_triangle_map, vertex_properties[0].second), - get(face_to_triangle_map, vertex_properties[1].second), - get(face_to_triangle_map, vertex_properties[2].second)}; - - for(int i=0; i<3; ++i) - { - // Iterate over the vertices by i, the distance to the closest point in - // tm2 computed above is a lower bound for the local triangle - h_triangle_lower = (std::max)(h_triangle_lower, vertex_properties[i].first); - - // Iterate over the triangles by i, if the triangles are the same, we do - // not need to compute the distance, only compute it if the triangles - // differ - double face_distance_1 = vertex_properties[i].second==vertex_properties[(i+1)%3].second - ? vertex_properties[(i+1)%3].first - : squared_distance(project_point(triangles_in_B[i], get(vpm1, face_vertices[(i+1)%3])), get(vpm1, face_vertices[(i+1)%3])); - double face_distance_2 = vertex_properties[i].second==vertex_properties[(i+2)%3].second - ? vertex_properties[(i+2)%3].first - : squared_distance(project_point(triangles_in_B[i], get(vpm1, face_vertices[(i+2)%3])), get(vpm1, face_vertices[(i+2)%3])); - - // Update the local lower bound of the triangle - h_triangle_upper = (std::min)( - (std::max)( - (std::max)(face_distance_1, face_distance_2), - vertex_properties[i].first), - h_triangle_upper); - } - - // Store the computed lower and upper bound in a dynamic face property - put(thb, fd, Hausdorff_bounds(h_triangle_lower, h_triangle_upper)); - h_lower = (std::max)(h_lower, h_triangle_lower); - h_upper = (std::max)(h_upper, h_triangle_upper); - - // Only process the triangle further if it can still contribute to a - // Hausdorff distance - if (h_triangle_upper > h_lower) { - - // TODO culling on B - - candidate_triangles.push_back(fd); - } - } - - - - // TODO Iterate over candidate_triangles and kill those which cannot contribute anymore - - // TODO Send the remaining triangles to the Subdivision - - return (CGAL::approximate_sqrt(h_lower)+CGAL::approximate_sqrt(h_upper))/2.; - } -*/ } template Intel TBB library. - * To control the number of threads used, the user may use the `tbb::task_scheduler_init` class. - * See the TBB documentation - * for more details. - * - * @tparam Concurrency_tag enables sequential versus parallel algorithm. - * Possible values are `Sequential_tag` - * and `Parallel_tag`. - * @tparam TriangleMesh a model of the concept `FaceListGraph` - * @tparam NamedParameters1 a sequence of \ref pmp_namedparameters "Named Parameters" for `tm1` - * @tparam NamedParameters2 a sequence of \ref pmp_namedparameters "Named Parameters" for `tm2` - * - * @param tm1 the triangle mesh that will be sampled - * @param tm2 the triangle mesh to compute the distance to - * @param np1 optional sequence of \ref pmp_namedparameters "Named Parameters" for `tm1` passed to `sample_triangle_mesh()`. - * - * @param np2 optional sequence of \ref pmp_namedparameters "Named Parameters" for `tm2` among the ones listed below - * - * \cgalNamedParamsBegin - * \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm2` - * If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` must be available in `TriangleMesh` - * and in all places where `vertex_point_map` is used. - * \cgalParamEnd - * \cgalNamedParamsEnd - * The function `CGAL::parameters::all_default()` can be used to indicate to use the default values for - * `np1` and specify custom values for `np2` - */ - /* * Implementation of Bounded Hausdorff distance computation using AABBTree * culling. 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 ccafde7ecc8..bfc57ee5895 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 @@ -133,58 +133,6 @@ namespace CGAL { // Hausdorff distance and thus have to be entered bool do_intersect(const Query& query, const Node& node) const { - /* Have reached a node, determine whether or not to enter it */ - - /* - // Get the bounding box of the nodes - Bounding_box bbox = node.bbox(); - // Compute its center - Point_3 bb_center = Point_3( - (bbox.min(0) + bbox.max(0)) / 2, - (bbox.min(1) + bbox.max(1)) / 2, - (bbox.min(2) + bbox.max(2)) / 2); - - // Get the center of the query triangle - // The query object is a triangle from TM1, get its vertices - Point_3 v0 = query.vertex(0); - Point_3 v1 = query.vertex(1); - Point_3 v2 = query.vertex(2); - // Compute the centroid of the triangle - Point_3 tri_center = centroid( query ); - - // Compute the distance of the center to the closest point in tm2 - double dist = approximate_sqrt(squared_distance(bb_center, tri_center)); - - // Compute the radius of the circumsphere of the bounding boxes - double bb_radius = approximate_sqrt(squared_distance( - Point_3(bbox.min(0),bbox.min(1),bbox.min(2)), - Point_3(bbox.max(0),bbox.max(1),bbox.max(2))) - )/2.; - - // Compute the radius of the circumcircle of the triangle - double tri_radius = approximate_sqrt( std::max( - std::max( - squared_distance(tri_center, v0), - squared_distance(tri_center, v1) - ), - squared_distance(tri_center, v2) - )); - - // Find a lower bound on the distance of the triangle to the bounding box - // by taking the distance from the triangle center to the bbox center - // and subtracting the radii of both elements. Thereby, each triangle in - // the bounding box has at least distance - // ( dist - triangle_radius - bbox_radius ) - // to the query triangle. Only if this lower bound is lower than the - // current best known lower bound, we enter the node, trying to - // lower the lower bound further. - if ( (dist - tri_radius - bb_radius) <= h_local_lower) { - return true; - } else { - return false; - } - */ - // Get the bounding box of the nodes Bounding_box bbox = node.bbox(); // Get the vertices of the query triangle @@ -287,8 +235,12 @@ namespace CGAL { Heap_type; public: - Hausdorff_primitive_traits_tm1(const AABBTraits& traits, const TM2_tree& tree, const TriangleMesh& tm1, const TriangleMesh& tm2 , const VPM1& vpm1, const VPM2& vpm2 ) - : m_traits(traits), m_tm2_tree(tree), m_tm1(tm1), m_tm2(tm2), m_vpm1(vpm1), m_vpm2(vpm2) { + Hausdorff_primitive_traits_tm1( + const AABBTraits& traits, const TM2_tree& tree, const TriangleMesh& tm1, + const TriangleMesh& tm2 , const VPM1& vpm1, const VPM2& vpm2, + const Point_3 hint ) + : m_traits(traits), m_tm2_tree(tree), m_tm1(tm1), m_tm2(tm2), + m_vpm1(vpm1), m_vpm2(vpm2) { // Initialize the global bounds with 0., they will only grow. h_lower = 0.; h_upper = 0.; @@ -310,10 +262,11 @@ namespace CGAL { Triangle_from_face_descriptor_map m_face_to_triangle_map( &m_tm1, m_vpm1 ); Triangle_3 candidate_triangle = get(m_face_to_triangle_map, primitive.id()); + // TODO Can we initialize the bounds here, s.t. we don't start with infty? + // Can we initialize the bounds depending on the closest points in tm2 + // seen from the three vertices? + // Call Culling on B with the single triangle found. - // TODO If we project the vertices of the triangle onto B and consider - // the edge length of the triangle, can this give us better bounds - // to initialize with here? Hausdorff_primitive_traits_tm2< Tree_traits, Triangle_3, Kernel, TriangleMesh, VPM2 > traversal_traits_tm2( @@ -349,7 +302,6 @@ namespace CGAL { bool do_intersect(const Query& /*query*/, const Node& node) const { /* Have reached a node, determine whether or not to enter it */ -/* // Get the bounding box of the nodes Bounding_box bbox = node.bbox(); // Compute its center @@ -358,31 +310,6 @@ namespace CGAL { (bbox.min(1) + bbox.max(1)) / 2, (bbox.min(2) + bbox.max(2)) / 2); // Find the point from TM2 closest to the center - // TODO Insert a hint here to accelerate the query - Point_3 closest = m_tm2_tree.closest_point(center); - // Compute the distance of the center to the closest point in tm2 - double dist = approximate_sqrt(squared_distance(center, closest)); - // Compute the radius of the circumsphere of the bounding boxes - double radius = approximate_sqrt(squared_distance( - Point_3(bbox.min(0),bbox.min(1),bbox.min(2)), - Point_3(bbox.max(0),bbox.max(1),bbox.max(2))) - )/2.; - // If the distance is larger than the global lower bound, enter the node, i.e. return true. - if (dist + radius > h_lower) { - return true; - } else { - return false; - } -*/ - // Get the bounding box of the nodes - Bounding_box bbox = node.bbox(); - // Compute its center - Point_3 center = Point_3( - (bbox.min(0) + bbox.max(0)) / 2, - (bbox.min(1) + bbox.max(1)) / 2, - (bbox.min(2) + bbox.max(2)) / 2); - // Find the point from TM2 closest to the center - // TODO Insert a hint here to accelerate the query Point_3 closest = m_tm2_tree.closest_point(center); // Compute the difference vector between the bbox center and the closest // point in tm2 From 68f7fac25435fb6957f9e865375056d48d71511a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 23 Aug 2019 09:49:08 +0200 Subject: [PATCH 36/48] update doc to make links working --- .../doc/Polygon_mesh_processing/PackageDescription.txt | 1 + .../doc/Polygon_mesh_processing/Polygon_mesh_processing.txt | 4 ++-- .../doc/Polygon_mesh_processing/examples.txt | 1 + .../include/CGAL/Polygon_mesh_processing/distance.h | 3 ++- 4 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index 032167d48e0..9051cce6097 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -176,6 +176,7 @@ and provides a list of the parameters that are used in this package. - \link measure_grp `CGAL::Polygon_mesh_processing::centroid()` \endlink \cgalCRPSection{Distance Functions} +- `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` - `CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance()` - `CGAL::Polygon_mesh_processing::approximate_symmetric_Hausdorff_distance()` - `CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set()` diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index 47f486bbc62..214b290b46d 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt @@ -775,11 +775,11 @@ with the following code: \subsection Bounded Hausdorff Distance -The function \link bounded_error_Hausdorff_distance() `bounded_error_Hausdorff_distance()`\endlink +The function `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` computes an estimate of the Hausdorff distance of two triangle meshes which is bounded by a user-given error bound. -\cgalExample{POlygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp} +\cgalExample{Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp} \section PMPDetectFeatures Feature Detection diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt index f5aa23b864b..a8bde70c86b 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt @@ -23,4 +23,5 @@ \example Polygon_mesh_processing/detect_features_example.cpp \example Polygon_mesh_processing/manifoldness_repair_example.cpp \example Polygon_mesh_processing/repair_polygon_soup_example.cpp +\example Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp */ 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 1bd5f39fe85..409f7cadd5f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1110,7 +1110,8 @@ double bounded_error_Hausdorff_naive_impl( } //end of namespace internal -/* +/** + * \ingroup PMP_distance_grp * Implementation of Bounded Hausdorff distance computation using AABBTree * culling. */ From 316cae0dd95f5f819727b05873c6b370e6a23443 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Fri, 23 Aug 2019 17:30:07 +0200 Subject: [PATCH 37/48] Finalize documentation of bounded Hausdorff Distance computation and make print outs optional via a debug parameter. --- Documentation/doc/biblio/geom.bib | 11 ++++ .../Polygon_mesh_processing.txt | 31 +++++++++- .../CGAL/Polygon_mesh_processing/distance.h | 59 ++++++++++++++++--- 3 files changed, 90 insertions(+), 11 deletions(-) diff --git a/Documentation/doc/biblio/geom.bib b/Documentation/doc/biblio/geom.bib index d4bf473db7f..56f88ba33b6 100644 --- a/Documentation/doc/biblio/geom.bib +++ b/Documentation/doc/biblio/geom.bib @@ -152016,3 +152016,14 @@ pages = {179--189} HAL_ID = {inria-00412437}, HAL_VERSION = {v1}, } + +@inproceedings{tang2009interactive, + title={Interactive Hausdorff distance computation for general polygonal models}, + author={Tang, Min and Lee, Minkyoung and Kim, Young J}, + booktitle={ACM Transactions on Graphics (TOG)}, + volume={28}, + number={3}, + pages={74}, + year={2009}, + organization={ACM} +} diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index 214b290b46d..bf0ab568159 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt @@ -733,7 +733,7 @@ which enables to treat one or several connected components as a face graph. This package provides methods to compute (approximate) distances between meshes and point sets. -\subsection Approximate Hausdorff Distance +\subsection ApproxHD Approximate Hausdorff Distance The function \link approximate_Hausdorff_distance() `approximate_Hausdorff_distance()`\endlink computes an approximation of the Hausdorff distance from a mesh `tm1` to a mesh `tm2`. Given a @@ -773,11 +773,36 @@ with the following code: \snippet Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp PMP_distance_snippet -\subsection Bounded Hausdorff Distance +\subsection BoundedHD Bounded Hausdorff Distance The function `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()` computes an estimate of the Hausdorff distance of two triangle meshes which is -bounded by a user-given error bound. +bounded by a user-given error bound. Given two meshes tm1 and tm2, it follows +the procedure given by \cgalCite{tang2009interactive}. Namely, an AABB tree is +built on tm1 and tm2 respectively. The tree on tm1 is used to iterate over all +triangles in tm1. For each triangle t in tm1, by traversing the tree on tm2, +by keeping track of current global bounds on the distance, it +is estimated whether this triangle can still contribute to the actual +Hausdorff distance. From this process, a set of candidate triangles is selected. + +The candidate triangles are subsequently subdivided and for each smaller +triangle, the tree on tm2 is traversed again. This is repeated until the +triangle is smaller than the user-given error bound, all vertices of the +triangle are projected onto the same triangle in tm2, or the triangle's upper +bound is lower than the global lower bound. After creation, the subdivided +triangles are added to the list of candidate triangles. Thereby, all candidate +triangles are processed until a triangle is found in which the Hausdorff +distance is realized or in which it is guaranteed to be realized within the +user-given error bound. + +\subsubsection BHDExample Bounded Hausdorff Distance Example + +In the following examples: (a) the distance of a tetrahedron to a remeshed +version of itself is computed, (b) the distance of two geometries is computed +which is realized strictly in the interior of a triangle of the first geometry, +(c) a perturbation of a user-given mesh is compared to the original user-given +mesh, (d) two user-given meshes are compared, where the second mesh is gradually +moved away from the first one. \cgalExample{Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp} 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 409f7cadd5f..111e60cfd16 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -890,15 +890,14 @@ double bounded_error_Hausdorff_impl( Heap_type candidate_triangles = traversal_traits_tm1.get_candidate_triangles(); Hausdorff_bounds global_bounds = traversal_traits_tm1.get_global_bounds(); - std::cout << "Culled " << traversal_traits_tm1.get_num_culled_triangles() << " out of " << tm1.num_faces() << std::endl; + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "Culled " << traversal_traits_tm1.get_num_culled_triangles() << " out of " << tm1.num_faces() << std::endl; + #endif double squared_error_bound = error_bound * error_bound; while ( (global_bounds.second - global_bounds.first > error_bound) && !candidate_triangles.empty() ) { - // std::cout << "Current number candidates: " << candidate_triangles.size() << std::endl; - // std::cout << "Current global bounds: (" << global_bounds.first << ", " << global_bounds.second << ")" << std::endl; - // Get the first triangle and its Hausdorff bounds from the candidate set Candidate_triangle triangle_and_bound = candidate_triangles.top(); // Remove it from the candidate set as it will be processed now @@ -910,8 +909,6 @@ double bounded_error_Hausdorff_impl( // user given error. Hausdorff_bounds triangle_bounds = triangle_and_bound.m_bounds; - // std::cout << "Current triangle bounds: (" << triangle_bounds.first << ", " << triangle_bounds.second << ")" << std::endl; - if ( (triangle_bounds.second > global_bounds.first) && (triangle_bounds.second - triangle_bounds.first > error_bound) ) { // Get the triangle that is to be subdivided and read its vertices Triangle_3 triangle_for_subdivision = triangle_and_bound.m_triangle; @@ -1112,8 +1109,21 @@ double bounded_error_Hausdorff_naive_impl( /** * \ingroup PMP_distance_grp - * Implementation of Bounded Hausdorff distance computation using AABBTree - * culling. + * returns an estimate on the Hausdorff distance between `tm1` and `tm2` that + * is at most `error_bound` away from the actual Hausdorff distance between + * the two given meshes. + * @tparam Concurrency_tag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag` + * and `Parallel_tag`. Currently, parall computation is + * not implemented, though. + * @tparam TriangleMesh a model of the concept `FaceListGraph` + * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" + * @param tm1 a triangle mesh + * @param tm2 a second triangle mesh + * @param error_bound Maximum bound by which the Hausdorff distance estimate is + * allowed to deviate from the actual Hausdorff distance. + * @param np1 an optional sequence of \ref pmp_namedparameters "Named Parameters" + * @param np2 an optional sequence of \ref pmp_namedparameters "Named Parameters" */ template< class Concurrency_tag, class TriangleMesh, @@ -1142,6 +1152,23 @@ double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, return internal::bounded_error_Hausdorff_impl(tm1, tm2, error_bound, vpm1, vpm2); } +/** + * \ingroup PMP_distance_grp + * returns an estimate on the Hausdorff distance between `tm1` and `tm2` that + * is at most `error_bound` away from the actual Hausdorff distance between + * the two given meshes. + * @tparam Concurrency_tag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag` + * and `Parallel_tag`. Currently, parall computation is + * not implemented, though. + * @tparam TriangleMesh a model of the concept `FaceListGraph` + * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" + * @param tm1 a triangle mesh + * @param tm2 a second triangle mesh + * @param error_bound Maximum bound by which the Hausdorff distance estimate is + * allowed to deviate from the actual Hausdorff distance. + * @param np1 an optional sequence of \ref pmp_namedparameters "Named Parameters" + */ template< class Concurrency_tag, class TriangleMesh, class NamedParameters1> @@ -1153,6 +1180,22 @@ double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, return bounded_error_Hausdorff_distance(tm1, tm2, error_bound, np1, parameters::all_default()); } +/** + * \ingroup PMP_distance_grp + * returns an estimate on the Hausdorff distance between `tm1` and `tm2` that + * is at most `error_bound` away from the actual Hausdorff distance between + * the two given meshes. + * @tparam Concurrency_tag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag` + * and `Parallel_tag`. Currently, parall computation is + * not implemented, though. + * @tparam TriangleMesh a model of the concept `FaceListGraph` + * @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" + * @param tm1 a triangle mesh + * @param tm2 a second triangle mesh + * @param error_bound Maximum bound by which the Hausdorff distance estimate is + * allowed to deviate from the actual Hausdorff distance. + */ template< class Concurrency_tag, class TriangleMesh> double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, From c5c6ef03baef04fe0cb010c4cd8c0ff2897d5479 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sat, 24 Aug 2019 11:09:25 +0200 Subject: [PATCH 38/48] Improve text in documentation. --- .../Polygon_mesh_processing/Polygon_mesh_processing.txt | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index bf0ab568159..5bf0d10fae8 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt @@ -780,10 +780,11 @@ computes an estimate of the Hausdorff distance of two triangle meshes which is bounded by a user-given error bound. Given two meshes tm1 and tm2, it follows the procedure given by \cgalCite{tang2009interactive}. Namely, an AABB tree is built on tm1 and tm2 respectively. The tree on tm1 is used to iterate over all -triangles in tm1. For each triangle t in tm1, by traversing the tree on tm2, -by keeping track of current global bounds on the distance, it -is estimated whether this triangle can still contribute to the actual -Hausdorff distance. From this process, a set of candidate triangles is selected. +triangles in tm1. Throughout the traversal, the procedure keeps track of a global +lower and upper bound on the Hausdorff distance respectively. For each triangle +t in tm1, by traversing the tree on tm2, it is estimated via the global bounds +whether t can still contribute to the actual Hausdorff distance. From this +process, a set of candidate triangles is selected. The candidate triangles are subsequently subdivided and for each smaller triangle, the tree on tm2 is traversed again. This is repeated until the From 5f2069a6b3618f2526fcfa083cc7315abe8699d8 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Sat, 24 Aug 2019 11:09:36 +0200 Subject: [PATCH 39/48] Document named parameters. --- .../CGAL/Polygon_mesh_processing/distance.h | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 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 111e60cfd16..fa3a8ae912e 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -1122,8 +1122,15 @@ double bounded_error_Hausdorff_naive_impl( * @param tm2 a second triangle mesh * @param error_bound Maximum bound by which the Hausdorff distance estimate is * allowed to deviate from the actual Hausdorff distance. - * @param np1 an optional sequence of \ref pmp_namedparameters "Named Parameters" - * @param np2 an optional sequence of \ref pmp_namedparameters "Named Parameters" + * @param np1 an optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below + * @param np2 an optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below + * \cgalNamedParamsBegin + * \cgalParamBegin{vertex_point_map} the property map with the points + * associated to the vertices of `tm`. If this parameter is omitted, + * an internal property map for `CGAL::vertex_point_t` + * must be available for `TriangleMesh`. + * \cgalParamEnd + * \cgalNamedParamsEnd */ template< class Concurrency_tag, class TriangleMesh, @@ -1167,7 +1174,14 @@ double bounded_error_Hausdorff_distance( const TriangleMesh& tm1, * @param tm2 a second triangle mesh * @param error_bound Maximum bound by which the Hausdorff distance estimate is * allowed to deviate from the actual Hausdorff distance. - * @param np1 an optional sequence of \ref pmp_namedparameters "Named Parameters" + * @param np1 an optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below + * \cgalNamedParamsBegin + * \cgalParamBegin{vertex_point_map} the property map with the points + * associated to the vertices of `tm`. If this parameter is omitted, + * an internal property map for `CGAL::vertex_point_t` + * must be available for `TriangleMesh`. + * \cgalParamEnd + * \cgalNamedParamsEnd */ template< class Concurrency_tag, class TriangleMesh, From 677908e405b024e01e4495b75ef581a6dca5e024 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 30 Aug 2019 13:57:56 +0200 Subject: [PATCH 40/48] Fix NP usage --- .../hausdorff_bounded_error_distance_example.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp index c0a18953ac2..491573b68f2 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -103,7 +103,7 @@ int main(int argc, char** argv) // Copy the mesh and perturb it slightly tm2 = tm1; bool do_project = false; - PMP::random_perturbation( tm2.vertices(), tm2, 0.1, do_project ); + PMP::random_perturbation( tm2.vertices(), tm2, 0.1, CGAL::parameters::do_project(do_project)); std::cout << "Perturbed the input mesh, now computing the Hausdorff distance." << std::endl; // Compute the Hausdorff distance From 6027b2d5f78f3e1847b9daeffd9c0a33a8eaab55 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 26 Sep 2019 15:15:50 +0200 Subject: [PATCH 41/48] Add Traits::intersection_with_priority() --- AABB_tree/include/CGAL/AABB_tree.h | 16 +++++ .../CGAL/internal/AABB_tree/AABB_node.h | 68 +++++++++++++++++++ .../CGAL/Polygon_mesh_processing/distance.h | 4 +- ...traversal_traits_with_Hausdorff_distance.h | 15 ++++ 4 files changed, 101 insertions(+), 2 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_tree.h b/AABB_tree/include/CGAL/AABB_tree.h index b57f0260a05..f7e7e131434 100644 --- a/AABB_tree/include/CGAL/AABB_tree.h +++ b/AABB_tree/include/CGAL/AABB_tree.h @@ -520,6 +520,22 @@ public: } } + + template + void traversal_with_priority(const Query& query, Traversal_traits& traits) const + { + switch(size()) + { + case 0: + break; + case 1: + traits.intersection(query, singleton_data()); + break; + default: // if(size() >= 2) + root_node()->template traversal_with_priority(query, traits, m_primitives.size()); + } + } + private: typedef AABB_node Node; diff --git a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h index 5d0cf5ce4d2..d4868b80157 100644 --- a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h +++ b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h @@ -86,6 +86,11 @@ public: Traversal_traits& traits, const std::size_t nb_primitives) const; + template + void traversal_with_priority(const Query& query, + Traversal_traits& traits, + const std::size_t nb_primitives) const; + private: typedef AABBTraits AABB_traits; typedef AABB_node Node; @@ -201,6 +206,69 @@ AABB_node::traversal(const Query& query, } } + + +template +template +void +AABB_node::traversal_with_priority(const Query& query, + Traversal_traits& traits, + const std::size_t nb_primitives) const +{ + // Recursive traversal + switch(nb_primitives) + { + case 2: + traits.intersection(query, left_data()); + if( traits.go_further() ) + { + traits.intersection(query, right_data()); + } + break; + case 3: + traits.intersection(query, left_data()); + if( traits.go_further() && traits.do_intersect(query, right_child()) ) + { + right_child().traversal_with_priority(query, traits, 2); + } + break; + default: + bool ileft, iright; + typename Traversal_traits::Priority pleft, pright; + std::tie(ileft,pleft) = traits.do_intersect_with_priority(query, left_child()); + std::tie(iright,pright) = traits.do_intersect_with_priority(query, right_child()); + + if(pleft >= pright) + { + if( ileft ) + { + left_child().traversal_with_priority(query, traits, nb_primitives/2); + if( traits.go_further() && traits.do_intersect(query, right_child()) ) + { + right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + } + } + else if( iright ) + { + right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + } + }else{ + if( iright ) + { + right_child().traversal_with_priority(query, traits, nb_primitives/2); + if( traits.go_further() && traits.do_intersect(query, left_child()) ) + { + left_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + } + } + else if( ileft ) + { + left_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + } + } + } +} + } // end namespace CGAL #endif // CGAL_AABB_NODE_H 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 fa3a8ae912e..f9307b1e8c4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -883,7 +883,7 @@ double bounded_error_Hausdorff_impl( // Build traversal traits for tm1_tree Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2, hint.first ); // Find candidate triangles in TM1 which might realise the Hausdorff bound - tm1_tree.traversal( Point_3(0,0,0), traversal_traits_tm1 ); // dummy point given as query as not needed + tm1_tree.traversal_with_priority( Point_3(0,0,0), traversal_traits_tm1 ); // dummy point given as query as not needed // TODO Is there a better/faster data structure than the Heap used here? // Can already build a sorted structure while collecting the candidates @@ -961,7 +961,7 @@ double bounded_error_Hausdorff_impl( std::numeric_limits::infinity(), std::numeric_limits::infinity() ); - tm2_tree.traversal(sub_triangles[i], traversal_traits_tm2); + tm2_tree.traversal_with_priority(sub_triangles[i], traversal_traits_tm2); // Update global lower Hausdorff bound according to the obtained local bounds Hausdorff_bounds local_bounds = traversal_traits_tm2.get_local_bounds(); 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 bfc57ee5895..b6a7e0804ef 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 @@ -65,6 +65,8 @@ namespace CGAL { typedef typename Kernel::Vector_3 Vector_3; public: + typedef int Priority; + Hausdorff_primitive_traits_tm2( const AABBTraits& traits, const TriangleMesh& tm2, @@ -187,6 +189,12 @@ namespace CGAL { } } + std::pair do_intersect_with_priority(const Query& query, const Node& node) const + { + bool b = do_intersect(query, node); + return std::make_pair(b, Priority(0)); + } + // Return the local Hausdorff bounds computed for the passed query triangle Hausdorff_bounds get_local_bounds() { @@ -235,6 +243,7 @@ namespace CGAL { Heap_type; public: + typedef int Priority; Hausdorff_primitive_traits_tm1( const AABBTraits& traits, const TM2_tree& tree, const TriangleMesh& tm1, const TriangleMesh& tm2 , const VPM1& vpm1, const VPM2& vpm2, @@ -335,6 +344,12 @@ namespace CGAL { } } + std::pair do_intersect_with_priority(const Query& query, const Node& node) const + { + bool b = do_intersect(query, node); + return std::make_pair(b, Priority(0)); + } + // Return those triangles from TM1 which are candidates for including a // point realizing the Hausdorff distance Heap_type get_candidate_triangles() { From a4aa6a7b770b58ad9be00a651293d59f53e1b967 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 26 Sep 2019 15:54:43 +0200 Subject: [PATCH 42/48] rework expression --- .../CGAL/internal/AABB_tree/AABB_node.h | 35 +++++++------------ 1 file changed, 13 insertions(+), 22 deletions(-) diff --git a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h index d4868b80157..94eacd277d7 100644 --- a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h +++ b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h @@ -206,8 +206,6 @@ AABB_node::traversal(const Query& query, } } - - template template void @@ -238,34 +236,27 @@ AABB_node::traversal_with_priority(const Query& query, std::tie(ileft,pleft) = traits.do_intersect_with_priority(query, left_child()); std::tie(iright,pright) = traits.do_intersect_with_priority(query, right_child()); - if(pleft >= pright) + if (ileft) { - if( ileft ) + if (iright) + { + if(pleft >= pright) { left_child().traversal_with_priority(query, traits, nb_primitives/2); - if( traits.go_further() && traits.do_intersect(query, right_child()) ) - { - right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); - } + if( traits.go_further() /* && traits.do_intersect(query, right_child()) */ ) // TODO shall we call again do_intersect? + right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); } - else if( iright ) - { - right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); - } - }else{ - if( iright ) + else { right_child().traversal_with_priority(query, traits, nb_primitives/2); - if( traits.go_further() && traits.do_intersect(query, left_child()) ) - { - left_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); - } - } - else if( ileft ) - { - left_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + if( traits.go_further() /* && traits.do_intersect(query, left_child()) */ ) // TODO shall we call again do_intersect? + left_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); } + } } + else + if (iright) + right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); } } From 3e5e721542ac6ad7852251818e68ff248670f207 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 26 Sep 2019 16:06:24 +0200 Subject: [PATCH 43/48] plug do_intersect_with_priority --- ...traversal_traits_with_Hausdorff_distance.h | 34 +++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) 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 b6a7e0804ef..0683688cf67 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 @@ -65,8 +65,8 @@ namespace CGAL { typedef typename Kernel::Vector_3 Vector_3; public: - typedef int Priority; - + typedef double Priority; + Hausdorff_primitive_traits_tm2( const AABBTraits& traits, const TriangleMesh& tm2, @@ -133,7 +133,7 @@ namespace CGAL { // Determine whether child nodes will still contribute to a smaller // Hausdorff distance and thus have to be entered - bool do_intersect(const Query& query, const Node& node) const + std::pair do_intersect_with_priority(const Query& query, const Node& node) const { // Get the bounding box of the nodes Bounding_box bbox = node.bbox(); @@ -183,18 +183,18 @@ namespace CGAL { // Check whether investigating the bbox can still lower the Hausdorff // distance. If so, enter the box. if ( dist <= h_local_lower) { - return true; + return std::make_pair(true, -dist); } else { - return false; + return std::make_pair(false, 0); } } - std::pair do_intersect_with_priority(const Query& query, const Node& node) const + bool do_intersect(const Query& /*query*/, const Node& /* node */) const { - bool b = do_intersect(query, node); - return std::make_pair(b, Priority(0)); + CGAL_assertion(false); + return false; } - + // Return the local Hausdorff bounds computed for the passed query triangle Hausdorff_bounds get_local_bounds() { @@ -243,7 +243,7 @@ namespace CGAL { Heap_type; public: - typedef int Priority; + typedef double Priority; Hausdorff_primitive_traits_tm1( const AABBTraits& traits, const TM2_tree& tree, const TriangleMesh& tm1, const TriangleMesh& tm2 , const VPM1& vpm1, const VPM2& vpm2, @@ -308,7 +308,7 @@ namespace CGAL { // Determine whether child nodes will still contribute to a larger // Hausdorff distance and thus have to be entered // TODO Document Query object, explain why I don't need it here. - bool do_intersect(const Query& /*query*/, const Node& node) const + std::pair do_intersect_with_priority(const Query& /*query*/, const Node& node) const { /* Have reached a node, determine whether or not to enter it */ // Get the bounding box of the nodes @@ -338,18 +338,18 @@ namespace CGAL { // If the distance is larger than the global lower bound, enter the node, i.e. return true. if (dist > h_lower) { - return true; + return std::make_pair(true, dist); } else { - return false; + return std::make_pair(false, 0); } } - std::pair do_intersect_with_priority(const Query& query, const Node& node) const + bool do_intersect(const Query& /*query*/, const Node& /* node */) const { - bool b = do_intersect(query, node); - return std::make_pair(b, Priority(0)); + CGAL_assertion(false); + return false; } - + // Return those triangles from TM1 which are candidates for including a // point realizing the Hausdorff distance Heap_type get_candidate_triangles() { From 89094af48d1c5b02c6c70667c144cc55ffca5d48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Thu, 26 Sep 2019 17:59:54 +0200 Subject: [PATCH 44/48] use do_intersect correctly as it is still needed --- .../AABB_traversal_traits_with_Hausdorff_distance.h | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) 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 0683688cf67..2cdc3ae9793 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 @@ -189,10 +189,9 @@ namespace CGAL { } } - bool do_intersect(const Query& /*query*/, const Node& /* node */) const + bool do_intersect(const Query& query, const Node& node ) const { - CGAL_assertion(false); - return false; + return this->do_intersect_with_priority(query, node).first; } // Return the local Hausdorff bounds computed for the passed query triangle @@ -286,7 +285,7 @@ namespace CGAL { std::numeric_limits::infinity(), std::numeric_limits::infinity() ); - m_tm2_tree.traversal(candidate_triangle, traversal_traits_tm2); + m_tm2_tree.traversal_with_priority(candidate_triangle, traversal_traits_tm2); // Update global Hausdorff bounds according to the obtained local bounds Hausdorff_bounds local_bounds = traversal_traits_tm2.get_local_bounds(); @@ -344,10 +343,9 @@ namespace CGAL { } } - bool do_intersect(const Query& /*query*/, const Node& /* node */) const + bool do_intersect(const Query& query, const Node& node ) const { - CGAL_assertion(false); - return false; + return this->do_intersect_with_priority(query, node).first; } // Return those triangles from TM1 which are candidates for including a From 6e73f92fd21c01a71944d2d503a58c72d4d6c197 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 27 Sep 2019 07:15:17 +0200 Subject: [PATCH 45/48] fix nb of primitives --- AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h index 94eacd277d7..2bbac2230c6 100644 --- a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h +++ b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h @@ -248,9 +248,9 @@ AABB_node::traversal_with_priority(const Query& query, } else { - right_child().traversal_with_priority(query, traits, nb_primitives/2); + right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); if( traits.go_further() /* && traits.do_intersect(query, left_child()) */ ) // TODO shall we call again do_intersect? - left_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); + left_child().traversal_with_priority(query, traits, nb_primitives/2); } } } From 16e0a475f40b81200907de3cdcf253f36f84a093 Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 14 Oct 2019 00:28:08 +0200 Subject: [PATCH 46/48] Fix traversal order, add comment on benchmark results. --- AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h index 2bbac2230c6..87e74fb18f0 100644 --- a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h +++ b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h @@ -240,22 +240,29 @@ AABB_node::traversal_with_priority(const Query& query, { if (iright) { + // Both children have to be inspected if(pleft >= pright) { + // Inspect left first, has higher priority left_child().traversal_with_priority(query, traits, nb_primitives/2); - if( traits.go_further() /* && traits.do_intersect(query, right_child()) */ ) // TODO shall we call again do_intersect? + if( traits.go_further() ) //&& traits.do_intersect(query, right_child()) ) // TODO shall we call again do_intersect? (Benchmarks show it slows down the Hausdorff Distance) right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); } else { + // Inspect right first, has higher priority right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); - if( traits.go_further() /* && traits.do_intersect(query, left_child()) */ ) // TODO shall we call again do_intersect? + if( traits.go_further() ) //&& traits.do_intersect(query, left_child()) ) // TODO shall we call again do_intersect? (Benchmarks show it slows down the Hausdorff Distance) left_child().traversal_with_priority(query, traits, nb_primitives/2); } } + else + // Only the left child has to be inspected + left_child().traversal_with_priority(query, traits, nb_primitives/2); } else if (iright) + // Only the right child has to be inspected right_child().traversal_with_priority(query, traits, nb_primitives-nb_primitives/2); } } From fed345c07a479f3887bcb2ef1b13bff69b62701e Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 14 Oct 2019 00:32:55 +0200 Subject: [PATCH 47/48] Add notes from skype meeting. --- .../include/CGAL/Polygon_mesh_processing/distance.h | 5 +++++ .../internal/AABB_traversal_traits_with_Hausdorff_distance.h | 2 ++ 2 files changed, 7 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 f9307b1e8c4..1e2d7ce0b0f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -34,6 +34,7 @@ #include #include #include +#include #include #include @@ -882,7 +883,11 @@ double bounded_error_Hausdorff_impl( // Build traversal traits for tm1_tree Hausdorff_primitive_traits_tm1 traversal_traits_tm1( tm1_tree.traits(), tm2_tree, tm1, tm2, vpm1, vpm2, hint.first ); + // Find candidate triangles in TM1 which might realise the Hausdorff bound +// TODO Initialize the distances on all the vertices first and store those. +// TODO Do not traverse TM1, but only TM2, i.e. reduce to Culling on TM2 (Can do this for all triangles in TM1 in parallel) + tm1_tree.traversal_with_priority( Point_3(0,0,0), traversal_traits_tm1 ); // dummy point given as query as not needed // TODO Is there a better/faster data structure than the Heap used here? 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 2cdc3ae9793..9e03cb17942 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 @@ -91,6 +91,7 @@ namespace CGAL { void intersection(const Query& query, const Primitive& primitive) { /* Have reached a single triangle, process it */ + // TODO Already perform these computations once we have <=k /* / Determine the distance accroding to @@ -285,6 +286,7 @@ namespace CGAL { std::numeric_limits::infinity(), std::numeric_limits::infinity() ); + // TODO Pass on the current global bounds to the TM2 tree traversal. There, only enter subtrees that can still be better than the current global bound. m_tm2_tree.traversal_with_priority(candidate_triangle, traversal_traits_tm2); // Update global Hausdorff bounds according to the obtained local bounds From 2d4c2543627c276e98e278af5a636b3e313c1d7d Mon Sep 17 00:00:00 2001 From: Martin Skrodzki Date: Mon, 14 Oct 2019 00:51:40 +0200 Subject: [PATCH 48/48] Keep current global bounds in mind when traversing tm2. --- .../AABB_traversal_traits_with_Hausdorff_distance.h | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) 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 9e03cb17942..a1d7e6ab1a2 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 @@ -71,11 +71,13 @@ namespace CGAL { const AABBTraits& traits, const TriangleMesh& tm2, const VPM2& vpm2, + const double h_upper_current_global, const double h_lower_init, const double h_upper_init, const double h_v0_lower_init, const double h_v1_lower_init, const double h_v2_lower_init ) : m_traits(traits), m_tm2(tm2), m_vpm2(vpm2) { // Initialize the global and local bounds with the given values + h_upper_global = h_upper_current_global; h_local_lower = h_lower_init; h_local_upper = h_upper_init; h_v0_lower = h_v0_lower_init; @@ -182,8 +184,8 @@ namespace CGAL { double dist = approximate_sqrt( Vector_3(dist_x,dist_y,dist_z).squared_length() ); // Check whether investigating the bbox can still lower the Hausdorff - // distance. If so, enter the box. - if ( dist <= h_local_lower) { + // distance and improve the current global bound. If so, enter the box. + if ( dist <= std::min(h_local_lower, h_upper_global) ) { return std::make_pair(true, -dist); } else { return std::make_pair(false, 0); @@ -207,6 +209,8 @@ namespace CGAL { const TriangleMesh& m_tm2; // its vertex point map const VPM2& m_vpm2; + // Current global upper bound on the Hausdorff distance + double h_upper_global; // Local Hausdorff bounds for the query triangle double h_local_upper; double h_local_lower; @@ -280,13 +284,13 @@ namespace CGAL { Tree_traits, Triangle_3, Kernel, TriangleMesh, VPM2 > traversal_traits_tm2( m_tm2_tree.traits(), m_tm2, m_vpm2, + (h_upper == 0.) ? std::numeric_limits::infinity() : h_upper, // Only pass current global bounds if they have been established yet std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity(), std::numeric_limits::infinity() ); - // TODO Pass on the current global bounds to the TM2 tree traversal. There, only enter subtrees that can still be better than the current global bound. m_tm2_tree.traversal_with_priority(candidate_triangle, traversal_traits_tm2); // Update global Hausdorff bounds according to the obtained local bounds