diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_hausdorff_bounded_error_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_hausdorff_bounded_error_distance.cpp index 2068553b08d..7e4799bbeb1 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_hausdorff_bounded_error_distance.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_hausdorff_bounded_error_distance.cpp @@ -1,24 +1,24 @@ // Use it to include parallel computations in the bounded error Hausdorff distance. -// #define USE_PARALLEL_BEHD +#define USE_PARALLEL_BEHD // Use this def in order to get all DEBUG info related to the bounded-error Hausdorff code! -// #define CGAL_HAUSDORFF_DEBUG +#define CGAL_HAUSDORFF_DEBUG -#include -#include #include #include #include -#include -#include - #include #include + +#include +#include +#include #include #include #include #include #include +#include using SCK = CGAL::Simple_cartesian; using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel; @@ -32,106 +32,164 @@ using Triangle_3 = typename Kernel::Triangle_3; using TAG = CGAL::Sequential_tag; using Surface_mesh = CGAL::Surface_mesh; -using Polyhedron = CGAL::Polyhedron_3; +using Polyhedron = CGAL::Surface_mesh; using Affine_transformation_3 = CGAL::Aff_transformation_3; using Timer = CGAL::Real_timer; using Face_handle = typename boost::graph_traits::face_descriptor; namespace PMP = CGAL::Polygon_mesh_processing; -struct Approximate_hd_wrapper { +struct Approximate_hd_wrapper +{ const double m_num_samples; + const double m_error_bound; + std::string name() const { return "approximate"; } - Approximate_hd_wrapper(const double num_samples) : m_num_samples(num_samples) { } - double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + + explicit Approximate_hd_wrapper(const double num_samples, + const double error_bound) + : m_num_samples(num_samples), m_error_bound(error_bound) + { } + + double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const + { return PMP::approximate_Hausdorff_distance(tm1, tm2, - CGAL::parameters::number_of_points_per_area_unit(m_num_samples), - CGAL::parameters::number_of_points_per_area_unit(m_num_samples)); + CGAL::parameters::number_of_points_on_faces(m_num_samples).random_seed(0), + CGAL::parameters::number_of_points_on_faces(m_num_samples).random_seed(0)); } - double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + + double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const + { return PMP::approximate_symmetric_Hausdorff_distance(tm1, tm2, - CGAL::parameters::number_of_points_per_area_unit(m_num_samples), - CGAL::parameters::number_of_points_per_area_unit(m_num_samples)); + CGAL::parameters::number_of_points_on_faces(m_num_samples).random_seed(0), + CGAL::parameters::number_of_points_on_faces(m_num_samples).random_seed(0)); + } + + bool check_value(const double val, const double expected) const + { + const bool res = ((CGAL::abs(val - expected)) <= m_error_bound); + if(!res) + std::cerr << "Expected " << expected << " but got " << val << std::endl; + + return res; } }; -struct Naive_bounded_error_hd_wrapper { +struct Naive_bounded_error_hd_wrapper +{ const double m_error_bound; + std::string name() const { return "naive bounded error"; } - Naive_bounded_error_hd_wrapper(const double error_bound) : m_error_bound(error_bound) { } - double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + + explicit Naive_bounded_error_hd_wrapper(const double error_bound) : m_error_bound(error_bound) { } + + double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const + { return PMP::bounded_error_Hausdorff_distance_naive(tm1, tm2, m_error_bound); } - double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + + double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const + { const double dista = operator()(tm1, tm2); const double distb = operator()(tm2, tm1); return (CGAL::max)(dista, distb); } + + bool check_value(const double val, const double expected) const + { + const bool res = ((CGAL::abs(val - expected)) <= m_error_bound); + if(!res) + std::cerr << "Expected " << expected << " but got " << val << std::endl; + + return res; + } }; -struct Bounded_error_hd_wrapper { +struct Bounded_error_hd_wrapper +{ const double m_error_bound; + std::string name() const { return "bounded error"; } - Bounded_error_hd_wrapper(const double error_bound) : m_error_bound(error_bound) { } - double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + + explicit Bounded_error_hd_wrapper(const double error_bound) : m_error_bound(error_bound) { } + + double operator()(const Surface_mesh& tm1, const Surface_mesh& tm2) const + { return PMP::bounded_error_Hausdorff_distance(tm1, tm2, m_error_bound); } - double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const { + + double symmetric(const Surface_mesh& tm1, const Surface_mesh& tm2) const + { return PMP::bounded_error_symmetric_Hausdorff_distance(tm1, tm2, m_error_bound); } + + bool check_value(const double val, const double expected) const + { + const bool res = ((CGAL::abs(val - expected)) <= m_error_bound); + if(!res) + std::cerr << "Expected " << expected << " but got " << val << std::endl; + + return res; + } }; template -void get_mesh(const std::string filepath, PolygonMesh& mesh) { - +void get_mesh(const std::string& filepath, PolygonMesh& mesh) +{ mesh.clear(); - std::ifstream input(filepath); - input >> mesh; + CGAL::IO::read_polygon_mesh(filepath, mesh); std::cout << "* getting mesh with " << faces(mesh).size() << " faces" << std::endl; } template -void get_meshes( - const std::string filepath1, const std::string filepath2, - PolygonMesh1& mesh1, PolygonMesh2& mesh2) { - +void get_meshes(const std::string& filepath1, + const std::string& filepath2, + PolygonMesh1& mesh1, + PolygonMesh2& mesh2) +{ get_mesh(filepath1, mesh1); get_mesh(filepath2, mesh2); } template -void save_mesh(const PolygonMesh& mesh, const std::string filepath) { - - if (!CGAL::IO::write_PLY(filepath + ".ply", mesh)) { +void save_mesh(const PolygonMesh& mesh, + const std::string& filepath) +{ + if(!CGAL::IO::write_PLY(filepath + ".ply", mesh)) + { std::cerr << "ERROR: cannot save this file: " << filepath << std::endl; exit(EXIT_FAILURE); } } // An easy example of a tetrahedron and its remeshed version. -void remeshing_tetrahedon_example( - const double error_bound, const bool save = false) { - +void remeshing_tetrahedon_example(const double error_bound, + const bool save = true) +{ Timer timer; Surface_mesh mesh1, mesh2; std::cout << std::endl << "* (E1) remeshing Tetrahedron example:" << std::endl; - CGAL::make_tetrahedron( - Point_3(0, 0, 0), Point_3(2, 0, 0), - Point_3(1, 1, 1), Point_3(1, 0, 2), mesh1); + CGAL::make_tetrahedron(Point_3(0, 0, 0), Point_3(2, 0, 0), + Point_3(1, 1, 1), Point_3(1, 0, 2), mesh1); mesh2 = mesh1; using edge_descriptor = typename boost::graph_traits::edge_descriptor; Surface_mesh::Property_map is_constrained_map = mesh2.add_property_map("e:is_constrained", true).first; - const double target_edge_length = 0.05; - PMP::isotropic_remeshing( - mesh2.faces(), target_edge_length, mesh2, - CGAL::parameters::edge_is_constrained_map(is_constrained_map)); + const double target_edge_length = 0.1; + PMP::isotropic_remeshing(mesh2.faces(), target_edge_length, mesh2, + CGAL::parameters::edge_is_constrained_map(is_constrained_map)); - if (save) save_mesh(mesh1, "mesh1"); - if (save) save_mesh(mesh2, "mesh2"); + std::cout << "Mesh1: " << num_vertices(mesh1) << " nv " << num_faces(mesh1) << " nf" << std::endl; + std::cout << "Mesh2: " << vertices(mesh2).size() << " nv " << faces(mesh2).size() << " nf" << std::endl; + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } timer.reset(); timer.start(); @@ -139,14 +197,13 @@ void remeshing_tetrahedon_example( std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; timer.stop(); std::cout << "* processing time: " << timer.time() << " s." << std::endl; - assert(hdist == error_bound); + assert(hdist <= error_bound); // real theoretical distance is 0 } -// Example with a point realizing the Hausdorff distance strictly -// lying in the interior of a triangle. -void interior_triangle_example( - const double error_bound, const bool save = false) { - +// Example with a point realizing the Hausdorff distance strictly lying in the interior of a triangle. +void interior_triangle_example(const double error_bound, + const bool save = true) +{ Timer timer; Surface_mesh mesh1, mesh2; std::cout << std::endl << "* (E2) interior Triangle example:" << std::endl; @@ -162,56 +219,66 @@ void interior_triangle_example( auto v3 = mesh2.add_vertex(Point_3( 0.0, 1, -1)); auto v4 = mesh2.add_vertex(Point_3(-0.5, 0, -1)); auto v5 = mesh2.add_vertex(Point_3( 0.5, 0, -1)); - mesh2.add_face(v0, v3, v4); - mesh2.add_face(v1, v4, v5); - mesh2.add_face(v2, v5, v3); + mesh2.add_face(v3, v4, v5); + mesh2.add_face(v0, v4, v3); + mesh2.add_face(v1, v5, v4); + mesh2.add_face(v2, v3, v5); - if (save) save_mesh(mesh1, "mesh1"); - if (save) save_mesh(mesh2, "mesh2"); + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } timer.reset(); timer.start(); const double hdist = PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound); - std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; timer.stop(); + std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; std::cout << "* processing time: " << timer.time() << " s." << std::endl; - assert(hdist >= 1.0); + + assert(1.59168 < hdist && hdist < 1.59169); } // Read a real mesh given by the user, perturb it slightly, and compute the // Hausdorff distance between the original mesh and its pertubation. -void perturbing_surface_mesh_example( - const std::string filepath, const double error_bound, const bool save = false) { - - Timer timer; +void perturbing_surface_mesh_example(const std::string& filepath, + const double error_bound, + const bool save = true) +{ std::cout << std::endl << "* (E3) perturbing Surface Mesh example:" << std::endl; Surface_mesh mesh1, mesh2; get_meshes(filepath, filepath, mesh1, mesh2); const double max_size = 0.1; - PMP::random_perturbation( - mesh2.vertices(), mesh2, max_size, CGAL::parameters::do_project(false)); - std::cout << "* perturbing the second mesh" << std::endl; + std::cout << "* perturbing the second mesh; max_size: " << max_size << std::endl; + PMP::random_perturbation(mesh2.vertices(), mesh2, max_size, CGAL::parameters::do_project(false) + .random_seed(0)); - if (save) save_mesh(mesh2, "mesh2"); + if(save) + { + save_mesh(mesh2, "mesh1"); + save_mesh(mesh1, "mesh2"); + } - timer.reset(); + Timer timer; timer.start(); const double hdist = PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound); - std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; timer.stop(); + std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; std::cout << "* processing time: " << timer.time() << " s." << std::endl; - assert(hdist > 0.0); + + assert(hdist > 0. && hdist <= max_size); } // Read two meshes and store them in two different face graph containers, // perturb the second mesh, and compute the Hausdorff distance. -void perturbing_polyhedron_mesh_example( - const std::string filepath, const double error_bound) { - - Timer timer; - std::cout << std::endl << "* (E3) perturbing Polyhedron mesh example:" << std::endl; +void perturbing_polyhedron_mesh_example(const std::string& filepath, + const double error_bound, + const bool save = true) +{ + std::cout << std::endl << "* (E3 bis) perturbing Polyhedron mesh example:" << std::endl; Surface_mesh mesh1; Polyhedron mesh2; @@ -219,72 +286,88 @@ void perturbing_polyhedron_mesh_example( get_mesh(filepath, mesh2); const double max_size = 0.1; - PMP::random_perturbation( - vertices(mesh2), mesh2, max_size, CGAL::parameters::do_project(false)); - std::cout << "* perturbing the second mesh" << std::endl; + std::cout << "* perturbing the second mesh; max_size: " << max_size << std::endl; + PMP::random_perturbation(vertices(mesh2), mesh2, max_size, CGAL::parameters::do_project(false) + .random_seed(0)); - timer.reset(); + if(save) + { + save_mesh(mesh2, "mesh1"); + save_mesh(mesh1, "mesh2"); + } + + Timer timer; timer.start(); const double hdist1 = PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound); - const double hdist2 = PMP::bounded_error_Hausdorff_distance(mesh2, mesh1, error_bound); std::cout << "* bounded-error Hausdorff distance 1->2: " << hdist1 << std::endl; + const double hdist2 = PMP::bounded_error_Hausdorff_distance(mesh2, mesh1, error_bound); std::cout << "* bounded-error Hausdorff distance 2->1: " << hdist2 << std::endl; timer.stop(); std::cout << "* processing time: " << timer.time() << " s." << std::endl; - assert(hdist1 > 0.0); - assert(hdist2 > 0.0); + assert(hdist1 > 0. && hdist1 <= max_size); + assert(hdist2 > 0.); assert(hdist2 > hdist1); const double hdist3 = PMP::bounded_error_symmetric_Hausdorff_distance(mesh1, mesh2, error_bound); + std::cout << "* symmetric bounded-error Hausdorff distance: " << hdist3 << std::endl; assert(hdist3 == (CGAL::max)(hdist1, hdist2)); } // 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. -void moving_surface_mesh_example( - const std::string filepath1, const std::string filepath2, - const std::size_t n, const double error_bound, const bool save = false) { - +void moving_surface_mesh_example(const std::string& filepath1, + const std::string& filepath2, + const std::size_t n, + const double error_bound, + const bool save = true) +{ Timer timer; std::cout << std::endl << "* (E4) moving Surface Mesh example:" << std::endl; Surface_mesh mesh1, mesh2; get_meshes(filepath1, filepath2, mesh1, mesh2); - const auto bbox = PMP::bbox(mesh2); - const FT distance = CGAL::approximate_sqrt(CGAL::squared_distance( - Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin()), - Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()))); + const CGAL::Bbox_3 bbox = PMP::bbox(mesh2); + const FT distance = CGAL::approximate_sqrt(CGAL::squared_distance(Point_3(bbox.xmin(), bbox.ymin(), bbox.zmin()), + Point_3(bbox.xmax(), bbox.ymax(), bbox.zmax()))); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const FT t = FT(1) / FT(100); - if (save) save_mesh(mesh2, "mesh-0"); + const FT a = t * distance; + const Vector_3 v(a, a, a); + const FT vl = std::sqrt(3) * t * distance; - double curr_dist = 0.0; - for (std::size_t i = 0; i < n; ++i) { - PMP::transform(Affine_transformation_3(CGAL::Translation(), - Vector_3(t * distance, t * distance, t * distance)), mesh2); + for(std::size_t i=0; i(mesh1, mesh2, error_bound); - std::cout << "* position: " << i << std::endl; - std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; timer.stop(); + std::cout << "* position: " << i << " translation distance: " << (i+1) * vl << std::endl; + std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl; std::cout << "* processing time: " << timer.time() << " s." << std::endl; - if (save) save_mesh(mesh2, "mesh-" + std::to_string(i + 1)); - assert(hdist > curr_dist); - curr_dist = hdist; + const FT l = (i+1) * vl - error_bound, u = (i+1) * vl + error_bound; + assert(hdist >= l && hdist <= u); } } template -void test_0(const FunctionWrapper& functor, const bool save = false) { - +void test_0(const FunctionWrapper& functor, + const bool save = true) +{ // The same triangle. // Expected distance is 0. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 0 ---- " << std::endl; @@ -292,32 +375,35 @@ void test_0(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); - mesh2 = mesh1; - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); std::cout << "* Hausdorff distance (expected 0): " << dista << std::endl; std::cout << "* HInverted distance (expected 0): " << distb << std::endl; std::cout << "* Symmetric distance (expected 0): " << distc << std::endl; - assert(dista == 0.0); - assert(distb == 0.0); - assert(distc == naive); + assert(functor.check_value(dista, 0)); + assert(functor.check_value(distb, 0)); + assert(functor.check_value(disto, distc)); } template -void test_1(const FunctionWrapper& functor, const bool save = false) { - - // Two triangles are parallel and 1 unit distance away from each other. +void test_1(const FunctionWrapper& functor, + const bool save = true) +{ + // Two parallel triangles, 1 unit distance away from each other. // Expected distance is 1. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 1 ---- " << std::endl; @@ -325,35 +411,39 @@ void test_1(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); mesh2.add_vertex(Point_3(0, 0, 1)); mesh2.add_vertex(Point_3(2, 0, 1)); mesh2.add_vertex(Point_3(1, 1, 1)); mesh2.add_face(mesh2.vertices()); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); std::cout << "* Hausdorff distance (expected 1): " << dista << std::endl; std::cout << "* HInverted distance (expected 1): " << distb << std::endl; std::cout << "* Symmetric distance (expected 1): " << distc << std::endl; - assert(dista == 1.0); - assert(distb == 1.0); - assert(distc == naive); + assert(functor.check_value(dista, 1.)); + assert(functor.check_value(distb, 1.)); + assert(functor.check_value(disto, distc)); } template -void test_2(const FunctionWrapper& functor, const bool save = false) { - +void test_2(const FunctionWrapper& functor, + const bool save = true) +{ // One triangle is orthogonal to the other one and shares a common edge. // Expected distance is 1. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 2 ---- " << std::endl; @@ -361,36 +451,40 @@ void test_2(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); mesh2.add_vertex(Point_3(0, 0, 0)); mesh2.add_vertex(Point_3(2, 0, 0)); mesh2.add_vertex(Point_3(1, 0, 1)); mesh2.add_face(mesh2.vertices()); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); std::cout << "* Hausdorff distance (expected 1): " << dista << std::endl; std::cout << "* HInverted distance (expected 1): " << distb << std::endl; std::cout << "* Symmetric distance (expected 1): " << distc << std::endl; - assert(dista == 1.0); - assert(distb == 1.0); - assert(distc == naive); + assert(functor.check_value(dista, 1.)); + assert(functor.check_value(distb, 1.)); + assert(functor.check_value(disto, distc)); } template -void test_3(const FunctionWrapper& functor, const bool save = false) { - +void test_3(const FunctionWrapper& functor, + const bool save = true) +{ // One triangle is orthogonal to the other one and shares a common edge // that is moved 1 unit distance away. - // Expected distances are sqrt(2) and 2. + // Expected distances are 2 and sqrt(2) ~= 1.4142135623730951. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 3 ---- " << std::endl; @@ -398,35 +492,39 @@ void test_3(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); mesh2.add_vertex(Point_3(0, 0, 1)); mesh2.add_vertex(Point_3(2, 0, 1)); mesh2.add_vertex(Point_3(1, 0, 2)); mesh2.add_face(mesh2.vertices()); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); - std::cout << "* Hausdorff distance (expected sqrt(2)): " << dista << std::endl; - std::cout << "* HInverted distance (expected 2 ): " << distb << std::endl; - std::cout << "* Symmetric distance (expected 2 ): " << distc << std::endl; + std::cout << "* Hausdorff distance (expected 1.4142135623730951): " << dista << std::endl; + std::cout << "* HInverted distance (expected 2): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 2): " << distc << std::endl; - assert(CGAL::abs(dista - CGAL::sqrt(2.0)) < 1e-5); // error bound is 1e-4 - assert(distb == 2.0); - assert(distc == naive); + assert(functor.check_value(dista, std::sqrt(2.))); + assert(functor.check_value(distb, 2.)); + assert(functor.check_value(disto, distc)); } template -void test_4(const FunctionWrapper& functor, const bool save = false) { - +void test_4(const FunctionWrapper& functor, + const bool save = true) +{ // One triangle is orthogonal to the other one and shares a common vertex. - // Expected distance is 1.2247448713915889407. + // Expected distance is std::sqrt(2) * std::sin(60°) ~= 1.2247448713915892 - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 4 ---- " << std::endl; @@ -434,37 +532,39 @@ void test_4(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); mesh2.add_vertex(Point_3(0, 1, 1)); mesh2.add_vertex(Point_3(2, 1, 1)); mesh2.add_vertex(Point_3(1, 1, 0)); mesh2.add_face(mesh2.vertices()); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); - std::cout << "* Hausdorff distance (expected 1.22): " << dista << std::endl; - std::cout << "* HInverted distance (expected 1.22): " << distb << std::endl; - std::cout << "* Symmetric distance (expected 1.22): " << distc << std::endl; + std::cout << "* Hausdorff distance (expected 1.2247448713915892): " << dista << std::endl; + std::cout << "* HInverted distance (expected 1.2247448713915892): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 1.2247448713915892): " << distc << std::endl; - assert(CGAL::abs(dista - 1.224744) < 1e-5); // error bound is 1e-4 - assert(CGAL::abs(distb - 1.224744) < 1e-5); // error bound is 1e-4 - assert(dista == distb); - assert(distc == naive); + assert(functor.check_value(dista, 1.2247448713915892)); + assert(functor.check_value(distb, 1.2247448713915892)); + assert(functor.check_value(disto, distc)); } template -void test_5(const FunctionWrapper& functor, const bool save = false) { +void test_5(const FunctionWrapper& functor, const bool save = true) { // One triangle is orthogonal to the other one and shares a common vertex // that is moved 1 unit distance away. - // Expected distances are 1.7320508075688771932 and 2.1213203435596423851. + // Expected distances are sqrt(3) ~= 1.7320508075688771 and sqrt(4.5) ~= 2.1213203435596423. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 5 ---- " << std::endl; @@ -472,36 +572,40 @@ void test_5(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); mesh2.add_vertex(Point_3(0, 1, 2)); mesh2.add_vertex(Point_3(2, 1, 2)); mesh2.add_vertex(Point_3(1, 1, 1)); mesh2.add_face(mesh2.vertices()); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); - std::cout << "* Hausdorff distance (expected 1.73): " << dista << std::endl; - std::cout << "* HInverted distance (expected 2.12): " << distb << std::endl; - std::cout << "* Symmetric distance (expected 2.12): " << distc << std::endl; + std::cout << "* Hausdorff distance (expected 1.7320508075688771): " << dista << std::endl; + std::cout << "* HInverted distance (expected 2.1213203435596423): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 2.1213203435596423): " << distc << std::endl; - assert(CGAL::abs(dista - 1.732050) < 1e-5); // error bound is 1e-4 - assert(CGAL::abs(distb - 2.121320) < 1e-5); // error bound is 1e-4 - assert(distc == naive); + assert(functor.check_value(dista, 1.7320508075688771)); + assert(functor.check_value(distb, 2.1213203435596423)); + assert(functor.check_value(disto, distc)); } template -void test_6(const FunctionWrapper& functor, const bool save = false) { - +void test_6(const FunctionWrapper& functor, + const bool save = true) +{ // The first and second mesh have different number of triangles. - // They are parallel and lie at the same plane. The middle triangle is overlapping. - // Expected distances are 0 and 0.70710678118654757274. + // They are parallel and lie in the same plane. The middle triangle is common to both meshes. + // Expected distances are 0 and 1/sqrt(2) ~=0.70710678118654746. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 6 ---- " << std::endl; @@ -509,7 +613,6 @@ void test_6(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); const auto v0 = mesh2.add_vertex(Point_3(0, 0, 0)); const auto v1 = mesh2.add_vertex(Point_3(2, 0, 0)); @@ -520,30 +623,35 @@ void test_6(const FunctionWrapper& functor, const bool save = false) { mesh2.add_face(v0, v1, v2); mesh2.add_face(v2, v1, v3); mesh2.add_face(v0, v2, v4); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); - std::cout << "* Hausdorff distance (expected 0.0): " << dista << std::endl; - std::cout << "* HInverted distance (expected 0.7): " << distb << std::endl; - std::cout << "* Symmetric distance (expected 0.7): " << distc << std::endl; + std::cout << "* Hausdorff distance (expected 0): " << dista << std::endl; + std::cout << "* HInverted distance (expected 0.70710678118654746): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 0.70710678118654746): " << distc << std::endl; - assert(dista == 0.0); - assert(CGAL::abs(distb - 0.707106) < 1e-5); // error bound is 1e-4 - assert(distc == naive); + assert(functor.check_value(dista, 0.)); + assert(functor.check_value(distb, 0.70710678118654746)); + assert(functor.check_value(disto, distc)); } template -void test_7(const FunctionWrapper& functor, const bool save = false) { - +void test_7(const FunctionWrapper& functor, + const bool save = true) +{ // One triangle is moved to 0.5 unit distance away from 3 other triangles. // The first and second meshes are parallel. - // Expected distances are 0.5 and 0.86602540378443859659. + // Expected distances are 0.5 and sqrt(0.75) ~= 0.8660254037844386. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 7 ---- " << std::endl; @@ -551,41 +659,44 @@ void test_7(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); const auto v0 = mesh2.add_vertex(Point_3(0, 0, 0.5)); const auto v1 = mesh2.add_vertex(Point_3(2, 0, 0.5)); const auto v2 = mesh2.add_vertex(Point_3(1, 1, 0.5)); const auto v3 = mesh2.add_vertex(Point_3(2, 1, 0.5)); const auto v4 = mesh2.add_vertex(Point_3(0, 1, 0.5)); - mesh2.add_face(v0, v1, v2); mesh2.add_face(v2, v1, v3); mesh2.add_face(v0, v2, v4); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); - std::cout << "* Hausdorff distance (expected 0.50): " << dista << std::endl; - std::cout << "* HInverted distance (expected 0.86): " << distb << std::endl; - std::cout << "* Symmetric distance (expected 0.86): " << distc << std::endl; + std::cout << "* Hausdorff distance (expected 0.5): " << dista << std::endl; + std::cout << "* HInverted distance (expected 0.8660254037844386): " << distb << std::endl; + std::cout << "* Symmetric distance (expected 0.8660254037844386): " << distc << std::endl; - assert(dista == 0.5); - assert(CGAL::abs(distb - 0.866025) < 1e-5); // error bound is 1e-4 - assert(distc == naive); + assert(functor.check_value(dista, 0.5)); + assert(functor.check_value(distb, 0.8660254037844386)); + assert(functor.check_value(disto, distc)); } template -void test_8(const FunctionWrapper& functor, const bool save = false) { - +void test_8(const FunctionWrapper& functor, + const bool save = true) +{ // One mesh has one triangle at zero level, another mesh has two triangles // where the first one is at level 1 and the second one is at level 2. // Expected distances are 1 and 2. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 8 ---- " << std::endl; @@ -593,7 +704,6 @@ void test_8(const FunctionWrapper& functor, const bool save = false) { mesh1.add_vertex(Point_3(2, 0, 0)); mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - if (save) save_mesh(mesh1, "mesh1"); const auto v0 = mesh2.add_vertex(Point_3(0, 0, 1)); const auto v1 = mesh2.add_vertex(Point_3(2, 0, 1)); @@ -601,33 +711,37 @@ void test_8(const FunctionWrapper& functor, const bool save = false) { const auto v3 = mesh2.add_vertex(Point_3(0, 0, 2)); const auto v4 = mesh2.add_vertex(Point_3(2, 0, 2)); const auto v5 = mesh2.add_vertex(Point_3(1, 1, 2)); - mesh2.add_face(v0, v1, v2); mesh2.add_face(v3, v4, v5); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); std::cout << "* Hausdorff distance (expected 1): " << dista << std::endl; std::cout << "* HInverted distance (expected 2): " << distb << std::endl; std::cout << "* Symmetric distance (expected 2): " << distc << std::endl; - assert(dista == 1.0); - assert(distb == 2.0); - assert(distc == naive); + assert(functor.check_value(dista, 1.)); + assert(functor.check_value(distb, 2.)); + assert(functor.check_value(disto, distc)); } template -void test_9(const FunctionWrapper& functor, const bool save = false) { - +void test_9(const FunctionWrapper& functor, + const bool save = true) +{ // Two meshes partially overlap, have 2 triangles in common and each one has // two its own trianles. All triangles form a Z shape where the height is 1. // The expected result is 1. - std::cout.precision(20); Surface_mesh mesh1, mesh2; std::cout << " ---- testing 9 ---- " << std::endl; @@ -641,7 +755,6 @@ void test_9(const FunctionWrapper& functor, const bool save = false) { mesh1.add_face(v2, v1, v3); mesh1.add_face(v1, v4, v3); mesh1.add_face(v3, v4, v5); - if (save) save_mesh(mesh1, "mesh1"); v0 = mesh2.add_vertex(Point_3(2, 0, 1)); v1 = mesh2.add_vertex(Point_3(1, 0, 0)); @@ -653,25 +766,30 @@ void test_9(const FunctionWrapper& functor, const bool save = false) { mesh2.add_face(v3, v4, v5); mesh2.add_face(v4, v0, v5); mesh2.add_face(v5, v0, v2); - if (save) save_mesh(mesh2, "mesh2"); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } const double dista = functor(mesh1, mesh2); const double distb = functor(mesh2, mesh1); - const double naive = (CGAL::max)(dista, distb); + const double disto = (CGAL::max)(dista, distb); const double distc = functor.symmetric(mesh1, mesh2); std::cout << "* Hausdorff distance (expected 1): " << dista << std::endl; std::cout << "* HInverted distance (expected 1): " << distb << std::endl; std::cout << "* Symmetric distance (expected 1): " << distc << std::endl; - assert(dista == 1.0); - assert(distb == 1.0); - assert(distc == naive); + assert(functor.check_value(dista, 1.)); + assert(functor.check_value(distb, 1.)); + assert(functor.check_value(disto, distc)); } template -void test_synthetic_data(const FunctionWrapper& functor) { - +void test_synthetic_data(const FunctionWrapper& functor) +{ std::cout << std::endl << "-- test synthetic data:" << std::endl << std::endl; std::cout << "* name -> " << functor.name() << std::endl; @@ -687,14 +805,12 @@ void test_synthetic_data(const FunctionWrapper& functor) { test_9(functor); // overlapping } -template< -typename FunctionWrapper1, -typename FunctionWrapper2> -void test_one_versus_another( - const FunctionWrapper1& functor1, - const FunctionWrapper2& functor2) { - - std::cout.precision(20); +template +void test_one_versus_another(const FunctionWrapper1& functor1, + const FunctionWrapper2& functor2, + const double error_bound) +{ const std::string filepath1 = "data/tetrahedron.off"; const std::string filepath2 = "data/tetrahedron-remeshed.off"; @@ -710,287 +826,126 @@ void test_one_versus_another( // The expected distance is 0. std::cout << " ---- testing 0 ---- " << std::endl; const double dista0 = functor1(mesh1, mesh2); - const double distb0 = functor2(mesh1, mesh2); std::cout << "* Hausdorff distance1: " << dista0 << std::endl; + const double distb0 = functor2(mesh1, mesh2); std::cout << "* Hausdorff distance2: " << distb0 << std::endl; + assert(dista0 <= error_bound); + assert(distb0 <= error_bound); - std::cout << "* traslating by 1 unit ..." << std::endl; - PMP::transform(Affine_transformation_3(CGAL::Translation(), - Vector_3(FT(0), FT(0), FT(1))), mesh2); + std::cout << "* translating by 1 unit ..." << std::endl; + PMP::transform(Affine_transformation_3(CGAL::Translation(), Vector_3(FT(0), FT(0), FT(1))), mesh2); // TEST 1. // Translate by 1 unit distance and compare. // The expected distance is 1. std::cout << " ---- testing 1 ---- " << std::endl; const double dista1 = functor1(mesh1, mesh2); - const double distb1 = functor2(mesh1, mesh2); std::cout << "* Hausdorff distance1: " << dista1 << std::endl; + const double distb1 = functor2(mesh1, mesh2); std::cout << "* Hausdorff distance2: " << distb1 << std::endl; - assert(CGAL::abs(dista0 - distb0) < 1e-3); // error bound is 1e-4 - assert(CGAL::abs(dista1 - distb1) < 1e-3); // error bound is 1e-4 + assert(std::abs(dista0 - distb0) <= 2. * error_bound); + assert(std::abs(dista1 - distb1) <= 2. * error_bound); } -template< -typename FunctionWrapper1, -typename FunctionWrapper2> -void test_real_meshes( - const std::string filepath1, - const std::string filepath2, - const FunctionWrapper1& functor1, - const FunctionWrapper2& functor2) { - - std::cout.precision(20); +template +void test_real_meshes(const std::string& filepath1, + const std::string& filepath2, + const FunctionWrapper1& functor1, + const FunctionWrapper2& functor2, + const double error_bound) +{ std::cout << std::endl << "-- test real meshes:" << std::endl << std::endl; + std::cout << "* name 1 -> " << functor1.name() << std::endl; + std::cout << "* name 2 -> " << functor2.name() << std::endl; std::cout << "* input path 1: " << filepath1 << std::endl; std::cout << "* input path 2: " << filepath2 << std::endl; Surface_mesh mesh1, mesh2; get_meshes(filepath1, filepath2, mesh1, mesh2); - std::cout << std::endl; - std::cout << "* name 1 -> " << functor1.name() << std::endl; - std::cout << "* name 2 -> " << functor2.name() << std::endl; - // Load and compare. - std::cout << std::endl; - std::cout << " ---- testing ---- " << std::endl; - const double dista0 = functor1(mesh1, mesh2); - const double dista1 = functor1(mesh2, mesh1); - const double distb0 = functor2(mesh1, mesh2); - const double distb1 = functor2(mesh2, mesh1); - std::cout << std::endl; - std::cout << "* Hausdorff distance1 f: " << dista0 << std::endl; - std::cout << "* Hausdorff distance1 b: " << dista1 << std::endl; - std::cout << std::endl; - std::cout << "* Hausdorff distance2 f: " << distb0 << std::endl; - std::cout << "* Hausdorff distance2 b: " << distb1 << std::endl; - - assert(CGAL::abs(dista0 - distb0) < 1e-3); // error bound is 1e-4 - assert(CGAL::abs(dista1 - distb1) < 1e-3); // error bound is 1e-4 -} - -template -void test_timings(const std::string filepath, const FunctionWrapper& functor) { - - std::cout.precision(20); - std::cout << std::endl << "-- test timings: " << functor.name() << std::endl << std::endl; - - Timer timer; - Surface_mesh mesh1, mesh2; - get_mesh(filepath, mesh1); - PMP::isotropic_remeshing(faces(mesh1), 0.005, mesh1); - mesh1.collect_garbage(); - mesh2 = mesh1; - - timer.reset(); - timer.start(); - const double dista1 = functor(mesh1, mesh2); - timer.stop(); - double timea = timer.time(); - - timer.reset(); - timer.start(); - const double distb1 = functor(mesh2, mesh1); - timer.stop(); - double timeb = timer.time(); - - timer.reset(); - timer.start(); - const double distc1 = functor.symmetric(mesh1, mesh2); - timer.stop(); - double timeab = timer.time(); - - std::cout << "* time a1 (sec.): " << timea << std::endl; - std::cout << "* time b1 (sec.): " << timeb << std::endl; - std::cout << "* time ab1 naive (sec.): " << timea + timeb << std::endl; - std::cout << "* time ab1 optimized (sec.): " << timeab << std::endl; - - assert(timea > 0.0); - assert(timeb > 0.0); - - // If the machine is loaded, this assert may go false. - // assert(timeab < timea + timeb); - - PMP::transform(Affine_transformation_3(CGAL::Translation(), - Vector_3(FT(0), FT(0), FT(1))), mesh2); - - timer.reset(); - timer.start(); - const double dista2 = functor(mesh1, mesh2); - timer.stop(); - timea = timer.time(); - - timer.reset(); - timer.start(); - const double distb2 = functor(mesh2, mesh1); - timer.stop(); - timeb = timer.time(); - - timer.reset(); - timer.start(); - const double distc2 = functor.symmetric(mesh1, mesh2); - timer.stop(); - timeab = timer.time(); - - std::cout << "* time a2 (sec.): " << timea << std::endl; - std::cout << "* time b2 (sec.): " << timeb << std::endl; - std::cout << "* time ab2 naive (sec.): " << timea + timeb << std::endl; - std::cout << "* time ab2 optimized (sec.): " << timeab << std::endl; - - assert(timea > 0.0); - assert(timeb > 0.0); - - // If the machine is loaded, this assert may go false. - // assert(timeab < timea + timeb); - - std::cout << "* dista = " << dista1 << " , " << dista2 << std::endl; - std::cout << "* distb = " << distb1 << " , " << distb2 << std::endl; - std::cout << "* distab = " << distc1 << " , " << distc2 << std::endl; - - assert(dista1 == distb1 && distb1 == distc1); - assert(dista2 == distb2 && distb2 == distc2); -} - -template -void test_bunny( - const FunctionWrapper& functor, const int n = 5, const bool save = false) { - - std::cout.precision(20); - const std::string filepath1 = "data/bunny_16300.off"; // approx 16.3K - const std::string filepath2 = "data/bunny_69400.off"; // approx 69.4K - - std::cout << std::endl << "-- test bunny:" << std::endl << std::endl; - std::cout << "* name -> " << functor.name() << std::endl; - - Surface_mesh mesh1, mesh2; - get_meshes(filepath1, filepath2, mesh1, mesh2); - if (save) save_mesh(mesh1, "mesh1"); - if (save) save_mesh(mesh2, "mesh2"); - - // Get 3 times longest dimension. - const auto bbox = PMP::bbox(mesh2); - const FT dist1 = static_cast(CGAL::abs(bbox.xmax() - bbox.xmin())); - const FT dist2 = static_cast(CGAL::abs(bbox.ymax() - bbox.ymin())); - const FT dist3 = static_cast(CGAL::abs(bbox.zmax() - bbox.zmin())); - const FT dim = FT(3) * (CGAL::max)((CGAL::max)(dist1, dist2), dist3); - - // Get timings. - Timer timer; - std::vector times; - times.reserve(n); - - if (n == 0) { - - const FT distance = dim; - PMP::transform(Affine_transformation_3(CGAL::Translation(), - Vector_3(distance, distance, distance)), mesh2); - if (save) save_mesh(mesh2, "mesh2"); - timer.reset(); + auto run = [&](auto f, const Surface_mesh& m1, const Surface_mesh& m2) -> std::pair + { + Timer timer; timer.start(); - // const double dista = functor(mesh1, mesh2); - const double distb = functor(mesh2, mesh1); - const double distc = distb; // (CGAL::max)(dista, distb); + const double dist = f(m1, m2); timer.stop(); - times.push_back(timer.time()); - std::cout << "* distance / Hausdorff / time (sec.) : " << - distance << " / " << distc << " / " << times.back() << std::endl; - assert(distc > 0.0); + return std::make_pair(dist, timer.time()); + }; - } else { + double timea0, timeb0, timea1, timeb1; + double dista0, distb0, dista1, distb1; + std::tie(dista0, timea0) = run(functor1, mesh1, mesh2); + std::tie(distb0, timeb0) = run(functor2, mesh1, mesh2); + std::tie(dista1, timea1) = run(functor1, mesh2, mesh1); + std::tie(distb1, timeb1) = run(functor2, mesh2, mesh1); - // t is the step where n is the number of steps. - double curr_dist = 1.0; - const FT t = FT(1) / static_cast(n); - for (int k = n; k >= 0; --k) { - auto mesh = mesh2; - const FT distance = k * t * dim; - PMP::transform(Affine_transformation_3(CGAL::Translation(), - Vector_3(distance, distance, distance)), mesh); - if (save) save_mesh(mesh, "mesh2-" + std::to_string(k)); + std::cout << std::endl; + std::cout << "* Hausdorff distance1-2 " << functor1.name() << " dist " << dista0 << " runtime: " << timea0 << std::endl; + std::cout << "* Hausdorff distance1-2 " << functor2.name() << " dist " << distb0 << " runtime: " << timeb0 << std::endl; + std::cout << std::endl; + std::cout << "* Hausdorff distance2-1 " << functor1.name() << " dist " << dista1 << " runtime: " << timea1 << std::endl; + std::cout << "* Hausdorff distance2-1 " << functor2.name() << " dist " << distb1 << " runtime: " << timeb1 << std::endl; - timer.reset(); - timer.start(); - const double dista = functor(mesh1, mesh); - const double distb = functor(mesh, mesh1); - const double distc = (CGAL::max)(dista, distb); - timer.stop(); - times.push_back(timer.time()); - std::cout << "* distance / Hausdorff / time (sec.) " << k << " : " << - distance << " / " << distc << " / " << times.back() << std::endl; - assert(distc < curr_dist); - curr_dist = distc; - } - } - - double min_time = +std::numeric_limits::infinity(); - double max_time = -std::numeric_limits::infinity(); - double avg_time = 0.0; - for (const double time : times) { - min_time = (CGAL::min)(min_time, time); - max_time = (CGAL::max)(max_time, time); - avg_time += time; - } - avg_time /= static_cast(times.size()); - - std::cout << std::endl << "* timings (msec.): " << std::endl; - std::cout << "* avg time: " << avg_time * 1000.0 << std::endl; - std::cout << "* min time: " << min_time * 1000.0 << std::endl; - std::cout << "* max time: " << max_time * 1000.0 << std::endl; - - assert(min_time <= max_time); - assert(min_time <= avg_time); - assert(avg_time <= max_time); + assert(std::abs(dista0 - distb0) <= 2. * error_bound); + assert(std::abs(dista1 - distb1) <= 2. * error_bound); } -Triangle_3 get_triangle(const int index, const Surface_mesh& mesh) { - +Triangle_3 get_triangle(const int index, + const Surface_mesh& mesh) +{ const auto mfaces = faces(mesh); assert(index >= 0); assert(index < static_cast(mfaces.size())); - const auto face = *(mfaces.begin() + index); + const auto face = *(std::next(mfaces.begin(), index)); const auto he = halfedge(face, mesh); const auto vertices = vertices_around_face(he, mesh); assert(vertices.size() >= 3); + auto vit = vertices.begin(); - const auto& p0 = mesh.point(*vit); ++vit; - const auto& p1 = mesh.point(*vit); ++vit; - const auto& p2 = mesh.point(*vit); + const Point_3& p0 = mesh.point(*vit); ++vit; + const Point_3& p1 = mesh.point(*vit); ++vit; + const Point_3& p2 = mesh.point(*vit); + return Triangle_3(p0, p1, p2); } -void compute_realizing_triangles( - const Surface_mesh& mesh1, const Surface_mesh& mesh2, - const double error_bound, const std::string prefix, const bool save = false) { - +void compute_realizing_triangles(const Surface_mesh& mesh1, + const Surface_mesh& mesh2, + const double error_bound, + const int expected_f1, + const int expected_f2, + const std::string& prefix, + const bool save = true) +{ std::cout << "* getting realizing triangles: " << std::endl; std::vector< std::pair > fpairs; - const double hdist = - PMP::bounded_error_Hausdorff_distance(mesh1, mesh2, error_bound, - CGAL::parameters::output_iterator(std::back_inserter(fpairs))); + const double hdist = PMP::bounded_error_Hausdorff_distance( + mesh1, mesh2, error_bound, + CGAL::parameters::output_iterator(std::back_inserter(fpairs))); assert(fpairs.size() == 2); // lower bound face pair + upper bound face pair + const int f1 = static_cast(fpairs.back().first); // f1 / f2: upper bound const int f2 = static_cast(fpairs.back().second); std::cout << "* Hausdorff: " << hdist << std::endl; std::cout << "* f1 / f2: " << f1 << " / " << f2 << std::endl; + std::cout << "* expected f1 / f2: " << expected_f1 << " / " << expected_f2 << std::endl; - assert(f1 == 0); - assert(f2 == 0 || f2 == 161); - - if (f1 != -1 && save) { - const auto triangle = get_triangle(f1, mesh1); + if(save) + { + auto triangle = get_triangle(f1, mesh1); Surface_mesh sm1; sm1.add_vertex(triangle[0]); sm1.add_vertex(triangle[1]); sm1.add_vertex(triangle[2]); sm1.add_face(sm1.vertices()); save_mesh(sm1, prefix + "triangle1"); - } - if (f2 != -1 && save) { - const auto triangle = get_triangle(f2, mesh2); + triangle = get_triangle(f2, mesh2); Surface_mesh sm2; sm2.add_vertex(triangle[0]); sm2.add_vertex(triangle[1]); @@ -998,12 +953,14 @@ void compute_realizing_triangles( sm2.add_face(sm2.vertices()); save_mesh(sm2, prefix + "triangle2"); } + + assert(f1 == expected_f1); + assert(f2 == expected_f2); } -void test_realizing_triangles( - const double error_bound, const bool save = false) { - - std::cout.precision(20); +void test_realizing_triangles(const double error_bound, + const bool save = true) +{ std::cout << std::endl << "-- test realizing triangles:" << std::endl << std::endl; // Basic test. @@ -1015,66 +972,63 @@ void test_realizing_triangles( mesh1.add_vertex(Point_3(1, 1, 0)); mesh1.add_face(mesh1.vertices()); - mesh2.add_vertex(Point_3(0, 0, error_bound / 2.0)); - mesh2.add_vertex(Point_3(2, 0, error_bound / 2.0)); - mesh2.add_vertex(Point_3(1, 1, error_bound / 2.0)); + mesh2.add_vertex(Point_3(0, 0, error_bound / 2.)); + mesh2.add_vertex(Point_3(2, 0, error_bound / 2.)); + mesh2.add_vertex(Point_3(1, 1, error_bound / 2.)); mesh2.add_face(mesh2.vertices()); - if (save) save_mesh(mesh1, "1mesh1"); - if (save) save_mesh(mesh2, "1mesh2"); + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } - compute_realizing_triangles(mesh1, mesh2, error_bound, "1", save); + compute_realizing_triangles(mesh1, mesh2, error_bound, 0, 0, "1", save); - mesh2.clear(); - mesh2.add_vertex(Point_3(0, 0, 1)); - mesh2.add_vertex(Point_3(2, 0, 1)); - mesh2.add_vertex(Point_3(1, 1, 1)); - mesh2.add_face(mesh2.vertices()); - - if (save) save_mesh(mesh2, "1mesh2"); - compute_realizing_triangles(mesh1, mesh2, error_bound, "1", save); - - // Complex test. std::cout << std::endl << " ---- complex test ---- " << std::endl; - mesh1.clear(); - mesh2.clear(); - const std::string filepath1 = "data/tetrahedron.off"; - const std::string filepath2 = "data/tetrahedron-remeshed.off"; - - std::array vhs1 = { - mesh1.add_vertex(Point_3(0, 1, 3)), - mesh1.add_vertex(Point_3(1, 1, 3)), - mesh1.add_vertex(Point_3(1, 0, 3)) - }; - mesh1.add_face(vhs1); - - std::array vhs2 = { - mesh2.add_vertex(Point_3(0, 1, 3.5)), - mesh2.add_vertex(Point_3(1, 1, 3.5)), - mesh2.add_vertex(Point_3(1, 0, 3.5)) - }; - mesh2.add_face(vhs2); Surface_mesh tmp1,tmp2; - get_meshes(filepath1, filepath2, tmp1, tmp2); + std::array vhs1 = + { + tmp1.add_vertex(Point_3(0, 1, 3)), + tmp1.add_vertex(Point_3(1, 1, 3)), + tmp1.add_vertex(Point_3(1, 0, 3)) + }; + tmp1.add_face(vhs1); - PMP::transform(Affine_transformation_3( - CGAL::Translation(), Vector_3(0, 0, 10 * error_bound)), tmp2); + std::array vhs2 = + { + tmp2.add_vertex(Point_3(0, 1, 3.5)), + tmp2.add_vertex(Point_3(1, 1, 3.5)), + tmp2.add_vertex(Point_3(1, 0, 3.5)) + }; + tmp2.add_face(vhs2); + + const std::string filepath1 = "data/tetrahedron.off"; + const std::string filepath2 = "data/tetrahedron-remeshed.off"; + mesh1.clear(); + mesh2.clear(); + get_meshes(filepath1, filepath2, mesh1, mesh2); + PMP::transform(Affine_transformation_3(CGAL::Translation(), Vector_3(0, 0, 10 * error_bound)), tmp2); mesh1.join(tmp1); mesh2.join(tmp2); - if (save) save_mesh(mesh1, "2mesh1"); - if (save) save_mesh(mesh2, "2mesh2"); + if(save) + { + save_mesh(mesh1, "mesh1-2"); + save_mesh(mesh2, "mesh2-2"); + } - compute_realizing_triangles(mesh1, mesh2, error_bound, "2", save); + const int expected_f1 = static_cast(faces(mesh1).size() - 1); + const int expected_f2 = static_cast(faces(mesh2).size() - 1); + compute_realizing_triangles(mesh1, mesh2, error_bound, expected_f1, expected_f2, "2", save); } #if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) -void test_parallel_version( - const std::string filepath, const double error_bound) { - - std::cout.precision(20); +void test_parallel_version(const std::string& filepath, + const double error_bound) +{ std::cout << std::endl << "-- test parallel version:" << std::endl << std::endl; Timer timer; @@ -1085,7 +1039,7 @@ void test_parallel_version( std::cout << " ---- one-sided distance ---- " << std::endl; PMP::transform(Affine_transformation_3(CGAL::Translation(), - Vector_3(FT(0), FT(0), FT(1))), mesh2); + Vector_3(FT(0), FT(0), FT(1))), mesh2); std::cout << " ---- SEQUENTIAL ---- " << std::endl; timer.reset(); @@ -1113,31 +1067,27 @@ void test_parallel_version( std::cout << "* dista seq = " << dista << std::endl; std::cout << "* distb par = " << distb << std::endl; - assert(timea > 0.0); - assert(timeb > 0.0); assert(dista == distb); } #endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) -void test_early_quit(const std::string filepath) { - - std::cout.precision(20); +void test_early_quit(const std::string& filepath, + const bool save = true) +{ std::cout << std::endl << "-- test early quit:" << std::endl << std::endl; Timer timer; Surface_mesh mesh1, mesh2; get_meshes(filepath, filepath, mesh1, mesh2); - std::cout << std::endl; - std::cout << " ---- distance 0.0 = 0.0 ---- " << std::endl; - timer.reset(); - timer.start(); - assert(!PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 0.0, 0.0001)); - timer.stop(); - const double timea = timer.time(); - PMP::transform(Affine_transformation_3(CGAL::Translation(), - Vector_3(0.0, 0.0, 0.5)), mesh2); + Vector_3(0.0, 0.0, 0.5)), mesh2); + + if(save) + { + save_mesh(mesh1, "mesh1"); + save_mesh(mesh2, "mesh2"); + } std::cout << " ---- distance 0.5 < 1.0 ---- " << std::endl; timer.reset(); @@ -1145,40 +1095,29 @@ void test_early_quit(const std::string filepath) { assert(!PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 1.0, 0.0001)); timer.stop(); const double timeb = timer.time(); - std::cout << " ---- distance 0.5 < 0.6 ---- " << std::endl; + + std::cout << std::endl << " ---- distance 0.5 < 0.6 ---- " << std::endl; timer.reset(); timer.start(); assert(!PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 0.6, 0.0001)); timer.stop(); const double timec = timer.time(); - std::cout << " ---- distance 0.5 > 0.4 ---- " << std::endl; + + std::cout << std::endl << " ---- distance 0.5 > 0.4 ---- " << std::endl; timer.reset(); timer.start(); assert(PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 0.4, 0.0001)); timer.stop(); const double timed = timer.time(); - std::cout << " ---- distance 0.5 > 0.0 ---- " << std::endl; - timer.reset(); - timer.start(); - assert(PMP::is_Hausdorff_distance_larger(mesh1, mesh2, 0.0, 0.0001)); - timer.stop(); - const double timee = timer.time(); - std::cout << "* timea 0.0 = 0.0 = " << timea << std::endl; - std::cout << "* timeb 0.5 < 1.0 = " << timeb << std::endl; - std::cout << "* timec 0.5 < 0.6 = " << timec << std::endl; - std::cout << "* timed 0.5 > 0.4 = " << timed << std::endl; - std::cout << "* timee 0.5 > 0.0 = " << timee << std::endl; - - assert(timea > 0.0); - assert(timeb > 0.0); - assert(timec > 0.0); - assert(timed > 0.0); - assert(timee > 0.0); + std::cout << "* timeb 0.5 < 1.0; time: " << timeb << std::endl; + std::cout << "* timec 0.5 < 0.6; time: " << timec << std::endl; + std::cout << "* timed 0.5 > 0.4; time: " << timed << std::endl; } -void run_examples(const double error_bound, const std::string filepath) { - +void run_examples(const double error_bound, + const std::string& filepath) +{ remeshing_tetrahedon_example(error_bound); interior_triangle_example(error_bound); perturbing_surface_mesh_example(filepath, error_bound); @@ -1186,73 +1125,58 @@ void run_examples(const double error_bound, const std::string filepath) { moving_surface_mesh_example(filepath, filepath, 5, error_bound); } -int main(int argc, char** argv) { - - // std::string name; - // std::cin >> name; +int main(int argc, char** argv) +{ + std::cout.precision(17); + std::cerr.precision(17); const double error_bound = 1e-4; - const double num_samples = 10.; + const double num_samples = 1000.; std::cout << std::endl << "* error bound: " << error_bound << std::endl; - // std::cout << std::endl << "* number of samples: " << num_samples << std::endl; - std::string filepath = (argc > 1 ? argv[1] : CGAL::data_file_path("meshes/blobby.off")); + std::cout << std::endl << "* number of samples: " << num_samples << std::endl; + + const std::string filepath = (argc > 1 ? argv[1] : CGAL::data_file_path("meshes/blobby.off")); run_examples(error_bound, filepath); // ------------------------------------------------------------------------ // // Tests. - // Approximate_hd_wrapper does not work with EPECK! - Approximate_hd_wrapper apprx_hd(num_samples); + Approximate_hd_wrapper apprx_hd(num_samples, error_bound); Naive_bounded_error_hd_wrapper naive_hd(error_bound); Bounded_error_hd_wrapper bound_hd(error_bound); // --- Testing basic properties. - - // test_synthetic_data(apprx_hd); - // test_synthetic_data(naive_hd); + test_synthetic_data(apprx_hd); +// test_synthetic_data(naive_hd); // naive takes too long test_synthetic_data(bound_hd); // --- Compare on common meshes. - - // test_one_versus_another(apprx_hd, naive_hd); - // test_one_versus_another(naive_hd, bound_hd); - test_one_versus_another(bound_hd, apprx_hd); +// test_one_versus_another(apprx_hd, naive_hd, error_bound); // naive takes too long +// test_one_versus_another(naive_hd, bound_hd, error_bound); // naive takes too long + test_one_versus_another(bound_hd, apprx_hd, error_bound); // --- Compare on real meshes. - const std::string filepath1 = (argc > 1 ? argv[1] : CGAL::data_file_path("meshes/blobby.off")); const std::string filepath2 = (argc > 2 ? argv[2] : "data/tetrahedron-remeshed.off"); - // test_real_meshes(filepath1, filepath2, apprx_hd, naive_hd); - // test_real_meshes(filepath1, filepath2, naive_hd, bound_hd); - test_real_meshes(filepath1, filepath2, bound_hd, apprx_hd); - - // --- Compare timings. - - filepath = (argc > 1 ? argv[1] : CGAL::data_file_path("meshes/blobby.off")); - // test_timings(filepath, apprx_hd); - // test_timings(filepath, naive_hd); - test_timings(filepath, bound_hd); - - // --- Compare with the paper. - - // test_bunny(apprx_hd); - // test_bunny(naive_hd); - // test_bunny(bound_hd, 3); +// test_real_meshes(filepath1, filepath2, apprx_hd, naive_hd, error_bound); +// test_real_meshes(filepath1, filepath2, naive_hd, bound_hd, error_bound); + test_real_meshes(filepath1, filepath2, bound_hd, apprx_hd, error_bound); +// test_real_meshes("data/elephant_concave_hole.off", CGAL::data_file_path("meshes/mech-holes-shark.off"), bound_hd, apprx_hd, error_bound); // commenting because approx hausdorff is annoyingly rough + test_real_meshes("data/small_spheres.off", "data/overlapping_triangles.off", bound_hd, apprx_hd, error_bound); // --- Test realizing triangles. test_realizing_triangles(error_bound); - #if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) +#if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD) // --- Test parallelization. test_parallel_version(filepath, error_bound); - #endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) +#endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) // --- Test early quit. test_early_quit(filepath); - // ------------------------------------------------------------------------ // + std::cout << "Done!" << std::endl; - std::cout << std::endl; return EXIT_SUCCESS; }