diff --git a/AABB_tree/include/CGAL/AABB_tree.h b/AABB_tree/include/CGAL/AABB_tree.h index 5f0521e219f..a0d35a0a99c 100644 --- a/AABB_tree/include/CGAL/AABB_tree.h +++ b/AABB_tree/include/CGAL/AABB_tree.h @@ -527,6 +527,21 @@ 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 d8edaa91d35..209da8536e6 100644 --- a/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h +++ b/AABB_tree/include/CGAL/internal/AABB_tree/AABB_node.h @@ -68,6 +68,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; @@ -152,6 +157,67 @@ 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 (ileft) + { + 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? (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? (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); + } +} + } // end namespace CGAL #endif // CGAL_AABB_NODE_H diff --git a/Documentation/doc/biblio/geom.bib b/Documentation/doc/biblio/geom.bib index ca708451d71..4b7a0483d27 100644 --- a/Documentation/doc/biblio/geom.bib +++ b/Documentation/doc/biblio/geom.bib @@ -152050,4 +152050,14 @@ pages = {179--189} Pages = {215--224}, Year = {2012}, Url = {http://monge.univ-mlv.fr/~colinde/pub/09edgewidth.pdf} + +@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/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index 860b966cdcd..a1d41428909 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -200,6 +200,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - \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 0a5ed04fdd6..0505ebe38a3 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 @@ -910,10 +910,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 Hausdorff Distance This package provides methods to compute (approximate) distances between meshes and point sets. +\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 a sampling of `tm1`, it computes the distance to `tm2` of the farthest sample point to `tm2` \cgalCite{cignoni1998metro}. @@ -938,12 +940,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, @@ -952,6 +954,40 @@ with the following code: \snippet Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp PMP_distance_snippet +\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. 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. 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 +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} + \section PMPDetectFeatures Feature Detection This package provides methods to detect some features of a polygon mesh. diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt index 6c19cbdbb15..1a5d1662b69 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt @@ -33,4 +33,5 @@ \example Polygon_mesh_processing/locate_example.cpp \example Polygon_mesh_processing/orientation_pipeline_example.cpp \example Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp +\example 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 e34c91b8892..a37ff600f70 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -37,6 +37,7 @@ include(CGAL_Eigen3_support) # ########################################################## create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp") +create_single_source_cgal_program( "hausdorff_bounded_error_distance_example.cpp") if(TARGET CGAL::Eigen3_support) create_single_source_cgal_program("hole_filling_example.cpp") @@ -123,6 +124,9 @@ if(TARGET CGAL::TBB_support) target_link_libraries(self_intersections_example PUBLIC CGAL::TBB_support) target_link_libraries(hausdorff_distance_remeshing_example PUBLIC CGAL::TBB_support) + target_link_libraries(hausdorff_bounded_error_distance_example + PUBLIC CGAL::TBB_support) + create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp") target_link_libraries(corefinement_parallel_union_meshes PUBLIC CGAL::TBB_support) 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..491573b68f2 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_bounded_error_distance_example.cpp @@ -0,0 +1,151 @@ +#include +#include +#include +#include +#include +#include +#include +#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_3; +typedef K::Triangle_3 Triangle_3; +typedef K::Vector_3 Vector_3; +typedef CGAL::Surface_mesh Mesh; +typedef Mesh::Vertex_index Vertex_index; + +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.0001; + +// ---------------------------------------------------------------------------- + + // 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; + 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.) ); + tm1.add_vertex( Point_3(0.,-1.,1.) ); + tm1.add_vertex( Point_3(1.,1.,1.) ); + tm1.add_face( tm1.vertices() ); + + 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 ); + + // 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 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, CGAL::parameters::do_project(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; + + 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() ); + + 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; + } +} 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 b945dccf702..1b42bf08a05 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -8,7 +8,7 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // -// 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 @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -26,6 +27,7 @@ #include #include #include +#include #include #ifdef CGAL_LINKED_WITH_TBB @@ -39,6 +41,7 @@ #include #include #include +#include namespace CGAL { namespace Polygon_mesh_processing { @@ -1310,8 +1313,467 @@ double approximate_symmetric_Hausdorff_distance(const TriangleMesh& tm1, tm1, tm2, parameters::all_default(), parameters::all_default()); } +//////////////////////////////////////////////////////////////////////// + +namespace internal { + +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 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 Tree_traits::Point_and_primitive_id Point_and_primitive_id; + + 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; + 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; + + typedef std::pair Hausdorff_bounds; + + 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(); + tm1_tree.accelerate_distance_queries(); + + // 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 + 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? + // 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(); + + #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() ) { + + // 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 + 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.m_bounds; + + 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; + 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. + // TODO Update the reference to the realizing triangle here as this is the best current guess. + global_bounds.first = triangle_bounds.second; + 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 ); + 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) + }; + + // 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_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(); + if (local_bounds.first > global_bounds.first) { + global_bounds.first = local_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( + Candidate_triangle(sub_triangles[i], local_bounds) + ); + } + + // 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); + } + } + + // Return linear interpolation between found lower and upper 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."); +#else + // TODO implement parallelized version of the above here. +#endif } -} // end of namespace CGAL::Polygon_mesh_processing + +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) +{ + // 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) ), + squared_distance( v1, tm2_tree.closest_point(v1) ) ), + squared_distance( v2, tm2_tree.closest_point(v2) ) + ); + } + + // 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 ); + + 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; + + // 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 + TM2_tree tm2_tree( faces(tm2).begin(), faces(tm2).end(), tm2, vpm2 ); + 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. + 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); + + // 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; + } + } + + // 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 + +/** + * \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" 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, + 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); +} + +/** + * \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" 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, + 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()); +} + +/** + * \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, + const TriangleMesh& tm2, + double error_bound) +{ + 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 #endif //CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H 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..a1d7e6ab1a2 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/AABB_traversal_traits_with_Hausdorff_distance.h @@ -0,0 +1,396 @@ +// 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 { + + 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 + */ + template + class Hausdorff_primitive_traits_tm2 + { + typedef typename AABBTraits::Primitive Primitive; + typedef ::CGAL::AABB_node Node; + 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: + typedef double Priority; + + Hausdorff_primitive_traits_tm2( + 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; + 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 + // 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, process it */ + // TODO Already perform these computations once we have <=k + + /* + / 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 + 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()), 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()), 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()), v2), v2 ) ); + if (v2_dist < h_v2_lower) h_v2_lower = v2_dist; + + // Get the distance as maximizers over all vertices + 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_upper is + // actually the correct Hausdorff distance from the query triangle in + // TM1 to the primitive triangle in TM2 + 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 + // Hausdorff distance and thus have to be entered + 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(); + // 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 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); + } + } + + bool do_intersect(const Query& query, const Node& node ) const + { + return this->do_intersect_with_priority(query, node).first; + } + + // 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 ); + } + + private: + const AABBTraits& m_traits; + // the mesh of this tree + 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; + double h_v0_lower; + double h_v1_lower; + double h_v2_lower; + }; + + + /** + * @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 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 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: + 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, + 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.; + // Initialize number of culled triangles + m_investigated_on_tm1 = 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, process it */ + m_investigated_on_tm1++; + + // 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()); + + // 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. + Hausdorff_primitive_traits_tm2< + 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() + ); + 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(); + + if (local_bounds.first > h_lower) { + h_lower = local_bounds.first; + } + if (local_bounds.second > h_upper) { + h_upper = local_bounds.second; + } + + // 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( Candidate_triangle(candidate_triangle, local_bounds) ); + pq.push( Candidate_triangle(candidate_triangle, local_bounds) ); + } + + // 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. + 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 + 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 = 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 std::make_pair(true, dist); + } else { + return std::make_pair(false, 0); + } + } + + bool do_intersect(const Query& query, const Node& node ) const + { + return this->do_intersect_with_priority(query, node).first; + } + + // Return those triangles from TM1 which are candidates for including a + // point realizing the Hausdorff distance + Heap_type 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 ); + } + + int get_num_culled_triangles() + { + return (m_tm1.num_faces() - m_investigated_on_tm1); + } + + 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& m_tm2_tree; + // 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 + 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; + }; +} + +#endif //CGAL_PMP_INTERNAL_AABB_TRAVERSAL_TRAITS_WITH_HAUSDORFF_DISTANCE 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 e88bfe425e9..83f7dc1bd86 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 @@ -294,7 +294,7 @@ int main(int argc, 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"; @@ -304,13 +304,34 @@ int main(int argc, 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"; + 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 " << error_bound << " larger than " + << PMP::bounded_error_Hausdorff_distance_naive( + m1,m2,error_bound) + << "\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,error_bound) + << "\n"; + time.stop(); + std::cout << "done in " << time.time() << "s.\n"; + general_tests(m1,m2); test_concept();