sebastien review

This commit is contained in:
Dmitry Anisimov 2021-06-18 14:44:09 +02:00
parent 6ae21a4379
commit 859bae9036
5 changed files with 448 additions and 280 deletions

View File

@ -202,6 +202,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
\cgalCRPSection{Distance Functions}
- `CGAL::Polygon_mesh_processing::bounded_error_Hausdorff_distance()`
- `CGAL::Polygon_mesh_processing::bounded_error_symmetric_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()`

View File

@ -959,7 +959,7 @@ with the following code:
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
the procedure given by \cgalCite{tang2009interactive}. Namely, bounded volume hierarchy (BVH) 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
@ -977,6 +977,10 @@ 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.
In the current implementation, the BVH used is an AABB-tree and not the swept sphere
volumes as used in the original implementation. This should explain the runtime difference
observed with the original implementation.
\subsubsection BHDExample Bounded Hausdorff Distance Example
In the following examples: (a) the distance of a tetrahedron to a remeshed

View File

@ -1,15 +1,11 @@
#include <CGAL/Bbox_3.h>
#include <CGAL/Real_timer.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Aff_transformation_3.h>
#include <CGAL/IO/PLY.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Polygon_mesh_processing/random_perturbation.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/Polygon_mesh_processing/transform.h>
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
@ -21,43 +17,19 @@ using TAG = CGAL::Sequential_tag;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
using Polyhedron = CGAL::Polyhedron_3<Kernel>;
using Affine_transformation_3 = CGAL::Aff_transformation_3<Kernel>;
using Timer = CGAL::Real_timer;
namespace PMP = CGAL::Polygon_mesh_processing;
template<class Mesh>
void get_mesh(const std::string filepath, Mesh& mesh) {
int main(int argc, char** argv) {
mesh.clear();
std::ifstream input(filepath);
input >> mesh;
std::cout << "* getting mesh with " << num_faces(mesh) << " faces" << std::endl;
}
const double error_bound = 1e-4;
const std::string filepath = (argc > 1 ? argv[1] : "data/blobby.off");
void get_meshes(
const std::string filepath1, const std::string filepath2,
Surface_mesh& mesh1, Surface_mesh& mesh2) {
// We create a tetrahedron, remesh it, and compute the distance.
// The expected distance is error_bound.
std::cout << std::endl << "* remeshing tetrahedron example:" << std::endl;
get_mesh(filepath1, mesh1);
get_mesh(filepath2, mesh2);
}
void save_mesh(const Surface_mesh& mesh, const std::string filepath) {
if (!CGAL::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) {
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);
@ -72,154 +44,26 @@ void remeshing_tetrahedon_example(
mesh2.faces(), target_edge_length, mesh2,
PMP::parameters::edge_is_constrained_map(is_constrained_map));
if (save) save_mesh(mesh1, "mesh1");
if (save) save_mesh(mesh2, "mesh2");
timer.reset();
timer.start();
std::cout << "* bounded Hausdorff distance: " <<
std::cout << "* one-sided bounded-error Hausdorff distance: " <<
PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound) << std::endl;
timer.stop();
std::cout << "* processing time: " << timer.time() << " s." << std::endl;
}
// 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) {
// We load a mesh, save it in two different containers, and
// translate the second mesh by 1 unit. The expected distance is 1.
std::cout << std::endl << "* moving mesh example:" << std::endl;
Timer timer;
Surface_mesh mesh1, mesh2;
std::cout << std::endl << "* (E2) interior triangle example:" << std::endl;
Surface_mesh surface_mesh;
std::ifstream sm_input(filepath);
sm_input >> surface_mesh;
mesh1.add_vertex(Point_3(-1, 1, 1));
mesh1.add_vertex(Point_3( 0, -1, 1));
mesh1.add_vertex(Point_3( 1, 1, 1));
mesh1.add_face(mesh1.vertices());
Polyhedron polyhedron;
std::ifstream poly_input(filepath);
poly_input >> polyhedron;
auto v0 = mesh2.add_vertex(Point_3(-1.0, 1, 0));
auto v1 = mesh2.add_vertex(Point_3( 0.0, -1, 0));
auto v2 = mesh2.add_vertex(Point_3( 1.0, 1, 0));
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);
PMP::transform(Affine_transformation_3(CGAL::Translation(),
Vector_3(FT(0), FT(0), FT(1))), polyhedron);
if (save) save_mesh(mesh1, "mesh1");
if (save) save_mesh(mesh2, "mesh2");
timer.reset();
timer.start();
std::cout << "* bounded Hausdorff distance: " <<
PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound) << std::endl;
timer.stop();
std::cout << "* processing time: " << timer.time() << " s." << std::endl;
}
// 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_mesh_example(
const std::string filepath, const double error_bound, const bool save = false) {
Timer timer;
std::cout << std::endl << "* (E3) perturbing mesh example (Surface Mesh):" << 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;
if (save) save_mesh(mesh2, "mesh2");
timer.reset();
timer.start();
std::cout << "* bounded Hausdorff distance: " <<
PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound) << std::endl;
timer.stop();
std::cout << "* processing time: " << timer.time() << " s." << std::endl;
}
// Read two meshes and store them in two different face graph containers,
// perturb the second mesh, and compute the Hausdorff distance.
void perturbing_mesh_example_with_polyhedron(
const std::string filepath, const double error_bound) {
Timer timer;
std::cout << std::endl << "* (E3) perturbing mesh example (Polyhedron):" << std::endl;
Surface_mesh mesh1;
Polyhedron mesh2;
get_mesh(filepath, mesh1);
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;
timer.reset();
timer.start();
std::cout << "* bounded Hausdorff distance 1->2: " <<
PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound) << std::endl;
std::cout << "* bounded Hausdorff distance 2->1: " <<
PMP::bounded_error_Hausdorff_distance<TAG>(mesh2, mesh1, error_bound) << std::endl;
timer.stop();
std::cout << "* processing time: " << timer.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.
void moving_mesh_example(
const std::string filepath1, const std::string filepath2,
const std::size_t n, const double error_bound, const bool save = false) {
Timer timer;
std::cout << std::endl << "* (E4) moving 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 FT t = FT(1) / FT(100);
if (save) save_mesh(mesh2, "mesh-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);
timer.reset();
timer.start();
std::cout << "* position: " << i << std::endl;
std::cout << "* bounded Hausdorff distance: " <<
PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound) << std::endl;
timer.stop();
std::cout << "* processing time: " << timer.time() << " s." << std::endl;
if (save) save_mesh(mesh2, "mesh-" + std::to_string(i + 1));
}
}
int main(int argc, char** argv) {
const double error_bound = 1e-4;
std::cout << std::endl << "* error bound: " << error_bound << std::endl;
std::string filepath = (argc > 1 ? argv[1] : "data/blobby.off");
remeshing_tetrahedon_example(error_bound);
interior_triangle_example(error_bound);
perturbing_mesh_example(filepath, error_bound);
perturbing_mesh_example_with_polyhedron(filepath, error_bound);
moving_mesh_example(filepath, filepath, 5, error_bound);
std::cout << std::endl;
std::cout << "* symmetric bounded-error Hausdorff distance: " <<
PMP::bounded_error_symmetric_Hausdorff_distance<TAG>(surface_mesh, polyhedron, error_bound)
<< std::endl << std::endl;
return EXIT_SUCCESS;
}

View File

@ -1031,12 +1031,12 @@ double approximate_Hausdorff_distance(
CGAL_assertion_code( bool is_triangle = is_triangle_mesh(tm) );
CGAL_assertion_msg (is_triangle,
"Mesh is not triangulated. Distance computing impossible.");
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "Nb sample points " << sample_points.size() << "\n";
#endif
typedef typename Kernel::Point_3 Point_3;
std::vector<Point_3> sample_points
(boost::begin(original_sample_points), boost::end(original_sample_points) );
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "Nb sample points " << sample_points.size() << "\n";
#endif
spatial_sort(sample_points.begin(), sample_points.end());
@ -1327,6 +1327,9 @@ double approximate_symmetric_Hausdorff_distance(const TriangleMesh& tm1,
// Use this def in order to get back the parallel version of the one-sided Hausdorff code!
// #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
namespace internal {
template< class Kernel,
@ -1356,13 +1359,14 @@ std::pair<typename Kernel::FT, bool> preprocess_bounded_error_Hausdorff_impl(
{
using FT = typename Kernel::FT;
using Point_3 = typename Kernel::Point_3;
using Timer = CGAL::Real_timer;
#ifdef CGAL_HAUSDORFF_DEBUG
using Timer = CGAL::Real_timer;
Timer timer;
timer.start();
// std::cout << "* preprocessing begin ...." << std::endl;
std::cout.precision(20);
std::cout << "* preprocessing begin ...." << std::endl;
std::cout.precision(17);
#endif
// Compute the max value that is used as infinity value for the given meshes.
// In our case, it is twice the length of the diagonal of the bbox of two input meshes.
@ -1392,9 +1396,11 @@ std::pair<typename Kernel::FT, bool> preprocess_bounded_error_Hausdorff_impl(
match_faces(tm1, tm2, std::back_inserter(common),
std::back_inserter(tm1_only), std::back_inserter(tm2_only), np1, np2);
// std::cout << "- common: " << common.size() << std::endl;
// std::cout << "- tm1 only: " << tm1_only.size() << std::endl;
// std::cout << "- tm2 only: " << tm2_only.size() << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "- common: " << common.size() << std::endl;
std::cout << "- tm1 only: " << tm1_only.size() << std::endl;
std::cout << "- tm2 only: " << tm2_only.size() << std::endl;
#endif
if (is_one_sided_distance) { // one-sided distance
@ -1436,9 +1442,11 @@ std::pair<typename Kernel::FT, bool> preprocess_bounded_error_Hausdorff_impl(
tm2_tree.insert(faces2.begin(), faces2.end(), tm2, vpm2);
}
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// std::cout << "* .... end preprocessing" << std::endl;
// std::cout << "* preprocessing time (sec.): " << timer.time() << std::endl;
std::cout << "* .... end preprocessing" << std::endl;
std::cout << "* preprocessing time (sec.): " << timer.time() << std::endl;
#endif
return std::make_pair(infinity_value, rebuild);
}
@ -1480,18 +1488,19 @@ double bounded_error_Hausdorff_impl(
using Face_handle_2 = typename boost::graph_traits<TriangleMesh2>::face_descriptor;
using Candidate = Candidate_triangle<Kernel, Face_handle_1, Face_handle_2>;
using Timer = CGAL::Real_timer;
std::cout.precision(20);
CGAL_precondition(error_bound >= FT(0));
CGAL_precondition(tm1_tree.size() > 0);
CGAL_precondition(tm2_tree.size() > 0);
// First, we apply culling.
// std::cout << "- applying culling" << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
using Timer = CGAL::Real_timer;
Timer timer;
timer.start();
std::cout << "- applying culling" << std::endl;
std::cout.precision(17);
#endif
// Build traversal traits for tm1_tree.
TM1_hd_traits traversal_traits_tm1(
@ -1506,12 +1515,16 @@ double bounded_error_Hausdorff_impl(
auto& candidate_triangles = traversal_traits_tm1.get_candidate_triangles();
auto global_bounds = traversal_traits_tm1.get_global_bounds();
// std::cout << "- number of candidate triangles: " << candidate_triangles.size() << std::endl;
// const FT culling_rate = FT(100) - (FT(candidate_triangles.size()) / FT(tm1_tree.size()) * FT(100));
// std::cout << "- culling rate: " << culling_rate << "%" << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "- number of candidate triangles: " << candidate_triangles.size() << std::endl;
const FT culling_rate = FT(100) - (FT(candidate_triangles.size()) / FT(tm1_tree.size()) * FT(100));
std::cout << "- culling rate: " << culling_rate << "%" << std::endl;
#endif
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// std::cout << "* culling (sec.): " << timer.time() << std::endl;
std::cout << "* culling (sec.): " << timer.time() << std::endl;
#endif
CGAL_assertion(global_bounds.lower >= FT(0));
CGAL_assertion(global_bounds.upper >= FT(0));
@ -1526,10 +1539,11 @@ double bounded_error_Hausdorff_impl(
CGAL_assertion(!traversal_traits_tm1.early_quit());
// Second, we apply subdivision.
// std::cout << "- applying subdivision" << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
timer.reset();
timer.start();
std::cout << "- applying subdivision" << std::endl;
#endif
// See Section 5.1 in the paper.
const FT squared_error_bound = error_bound * error_bound;
@ -1657,8 +1671,10 @@ double bounded_error_Hausdorff_impl(
}
}
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// std::cout << "* subdivision (sec.): " << timer.time() << std::endl;
std::cout << "* subdivision (sec.): " << timer.time() << std::endl;
#endif
// Compute linear interpolation between the found lower and upper bounds.
CGAL_assertion(global_bounds.lower >= FT(0));
@ -1667,12 +1683,6 @@ double bounded_error_Hausdorff_impl(
const double hdist = CGAL::to_double((global_bounds.lower + global_bounds.upper) / FT(2));
// Get realizing triangles.
// std::cout << "- found lface 1:" << static_cast<int>(global_bounds.lpair.first) << std::endl;
// std::cout << "- found lface 2:" << static_cast<int>(global_bounds.lpair.second) << std::endl;
// std::cout << "- found uface 1:" << static_cast<int>(global_bounds.upair.first) << std::endl;
// std::cout << "- found uface 2:" << static_cast<int>(global_bounds.upair.second) << std::endl;
CGAL_assertion(global_bounds.lpair.first != boost::graph_traits<TriangleMesh1>::null_face());
CGAL_assertion(global_bounds.lpair.second != boost::graph_traits<TriangleMesh2>::null_face());
CGAL_assertion(global_bounds.upair.first != boost::graph_traits<TriangleMesh1>::null_face());
@ -1708,7 +1718,9 @@ struct Triangle_mesh_wrapper {
template<class TM1Wrapper, class TM2Wrapper>
struct Bounded_error_preprocessing {
#ifdef CGAL_HAUSDORFF_DEBUG
using Timer = CGAL::Real_timer;
#endif
std::vector<boost::any>& tm_wrappers;
// Constructor.
@ -1731,9 +1743,14 @@ struct Bounded_error_preprocessing {
// TODO: make AABB tree build parallel!
void operator()(const tbb::blocked_range<std::size_t>& range) {
#ifdef CGAL_HAUSDORFF_DEBUG
Timer timer;
timer.reset();
timer.start();
std::cout.precision(17);
#endif
for (std::size_t i = range.begin(); i != range.end(); ++i) {
CGAL_assertion(i < tm_wrappers.size());
auto& tm_wrapper = tm_wrappers[i];
@ -1747,8 +1764,11 @@ struct Bounded_error_preprocessing {
CGAL_assertion_msg(false, "Error: wrong boost any type!");
}
}
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// std::cout << "* time operator() preprocessing (sec.): " << timer.time() << std::endl;
std::cout << "* time operator() preprocessing (sec.): " << timer.time() << std::endl;
#endif
}
void join(Bounded_error_preprocessing&) { }
@ -1764,7 +1784,9 @@ template< class TriangleMesh1,
struct Bounded_error_distance_computation {
using FT = typename Kernel::FT;
#ifdef CGAL_HAUSDORFF_DEBUG
using Timer = CGAL::Real_timer;
#endif
const std::vector<TriangleMesh1>& tm1_parts; const TriangleMesh2& tm2;
const FT error_bound; const VPM1& vpm1; const VPM2& vpm2;
@ -1796,9 +1818,14 @@ struct Bounded_error_distance_computation {
}
void operator()(const tbb::blocked_range<std::size_t>& range) {
#ifdef CGAL_HAUSDORFF_DEBUG
Timer timer;
timer.reset();
timer.start();
std::cout.precision(17);
#endif
double hdist = -1.0;
auto stub = CGAL::Emptyset_iterator();
@ -1816,8 +1843,11 @@ struct Bounded_error_distance_computation {
if (dist > hdist) hdist = dist;
}
if (hdist > distance) distance = hdist;
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// std::cout << "* time operator() computation (sec.): " << timer.time() << std::endl;
std::cout << "* time operator() computation (sec.): " << timer.time() << std::endl;
#endif
}
void join(Bounded_error_distance_computation& rhs) {
@ -1871,8 +1901,11 @@ double bounded_error_one_sided_Hausdorff_impl(
using Face_handle_1 = typename boost::graph_traits<TM1>::face_descriptor;
using Face_handle_2 = typename boost::graph_traits<TM2>::face_descriptor;
using Timer = CGAL::Real_timer;
// This is parallel version: we split the tm1 into parts, build trees for all parts, and
// run in parallel all BHD computations. The final distance is obtained by taking the max
// between BHDs computed for these parts with respect to tm2.
// This is off by default because the parallel version does not show much of runtime improvement.
// The slowest part is building AABB trees and this is what should be accelerated in the future.
#if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD)
using Point_3 = typename Kernel::Point_3;
using TMF = CGAL::Face_filtered_graph<TM1>;
@ -1887,8 +1920,11 @@ double bounded_error_one_sided_Hausdorff_impl(
std::vector<boost::any> tm_wrappers;
#endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED)
#ifdef CGAL_HAUSDORFF_DEBUG
using Timer = CGAL::Real_timer;
Timer timer;
std::cout.precision(20);
std::cout.precision(17);
#endif
TM1_tree tm1_tree;
TM2_tree tm2_tree;
@ -1898,15 +1934,19 @@ double bounded_error_one_sided_Hausdorff_impl(
// TODO: add to NP!
const int nb_cores = 4;
const std::size_t min_nb_faces_to_split = 100; // TODO: increase this number?
// std::cout << "* num cores: " << nb_cores << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "* num cores: " << nb_cores << std::endl;
#endif
if (
boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value &&
nb_cores > 1 && faces(tm1).size() >= min_nb_faces_to_split) {
// (0) -- Compute infinity value.
#ifdef CGAL_HAUSDORFF_DEBUG
timer.reset();
timer.start();
#endif
const auto bbox1 = bbox(tm1);
const auto bbox2 = bbox(tm2);
const auto bb = bbox1 + bbox2;
@ -1915,39 +1955,54 @@ double bounded_error_one_sided_Hausdorff_impl(
Point_3(bb.xmax(), bb.ymax(), bb.zmax()));
infinity_value = CGAL::approximate_sqrt(sq_dist) * FT(2);
CGAL_assertion(infinity_value >= FT(0));
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// const double time0 = timer.time();
// std::cout << "- computing infinity (sec.): " << time0 << std::endl;
const double time0 = timer.time();
std::cout << "- computing infinity (sec.): " << time0 << std::endl;
#endif
// (1) -- Create partition of tm1.
#ifdef CGAL_HAUSDORFF_DEBUG
timer.reset();
timer.start();
#endif
using Face_property_tag = CGAL::dynamic_face_property_t<int>;
auto face_pid_map = get(Face_property_tag(), tm1);
CGAL::METIS::partition_graph(
tm1, nb_cores, CGAL::parameters::
face_partition_id_map(face_pid_map));
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// const double time1 = timer.time();
// std::cout << "- computing partition time (sec.): " << time1 << std::endl;
const double time1 = timer.time();
std::cout << "- computing partition time (sec.): " << time1 << std::endl;
#endif
// (2) -- Create a filtered face graph for each part.
#ifdef CGAL_HAUSDORFF_DEBUG
timer.reset();
timer.start();
#endif
tm1_parts.reserve(nb_cores);
for (int i = 0; i < nb_cores; ++i) {
tm1_parts.emplace_back(tm1, i, face_pid_map);
// CGAL_assertion(tm1_parts.back().is_selection_valid()); // TODO: why is it triggered sometimes?
// std::cout << "- part " << i << " size: " << tm1_parts.back().number_of_faces() << std::endl;
// TODO: why is it triggered sometimes?
// CGAL_assertion(tm1_parts.back().is_selection_valid());
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "- part " << i << " size: " << tm1_parts.back().number_of_faces() << std::endl;
#endif
}
CGAL_assertion(tm1_parts.size() == nb_cores);
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// const double time2 = timer.time();
// std::cout << "- creating graphs time (sec.): " << time2 << std::endl;
const double time2 = timer.time();
std::cout << "- creating graphs time (sec.): " << time2 << std::endl;
#endif
// (3) -- Preprocess all input data.
#ifdef CGAL_HAUSDORFF_DEBUG
timer.reset();
timer.start();
#endif
tm1_trees.resize(tm1_parts.size());
tm_wrappers.reserve(tm1_parts.size() + 1);
for (std::size_t i = 0; i < tm1_parts.size(); ++i) {
@ -1957,20 +2012,26 @@ double bounded_error_one_sided_Hausdorff_impl(
CGAL_assertion(tm_wrappers.size() == tm1_parts.size() + 1);
Bounded_error_preprocessing<TM1_wrapper, TM2_wrapper> bep(tm_wrappers);
tbb::parallel_reduce(tbb::blocked_range<std::size_t>(0, tm_wrappers.size()), bep);
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// const double time3 = timer.time();
// std::cout << "- creating trees time (sec.) " << time3 << std::endl;
const double time3 = timer.time();
std::cout << "- creating trees time (sec.) " << time3 << std::endl;
#endif
// Final timing.
// std::cout << "* preprocessing parallel time (sec.) " <<
// time0 + time1 + time2 + time3 << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "* preprocessing parallel time (sec.) " <<
time0 + time1 + time2 + time3 << std::endl;
#endif
} else // sequential version
#endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED)
{
// std::cout << "* preprocessing sequential version " << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
timer.reset();
timer.start();
std::cout << "* preprocessing sequential version " << std::endl;
#endif
bool rebuild = false;
std::vector<Face_handle_1> tm1_only;
std::vector<Face_handle_2> tm2_only;
@ -1984,13 +2045,19 @@ double bounded_error_one_sided_Hausdorff_impl(
tm1_tree.do_not_accelerate_distance_queries();
tm2_tree.accelerate_distance_queries();
}
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// std::cout << "* preprocessing sequential time (sec.) " << timer.time() << std::endl;
std::cout << "* preprocessing sequential time (sec.) " << timer.time() << std::endl;
#endif
}
// std::cout << "* infinity_value: " << infinity_value << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "* infinity_value: " << infinity_value << std::endl;
#endif
if (infinity_value < FT(0)) {
// std::cout << "* culling rate: 100%" << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "* culling rate: 100%" << std::endl;
#endif
const auto face1 = *(faces(tm1).begin());
const auto face2 = *(faces(tm2).begin());
*out++ = std::make_pair(face1, face2);
@ -2002,15 +2069,19 @@ double bounded_error_one_sided_Hausdorff_impl(
const FT initial_bound = error_bound;
std::atomic<double> hdist;
#ifdef CGAL_HAUSDORFF_DEBUG
timer.reset();
timer.start();
#endif
#if defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED) && defined(USE_PARALLEL_BEHD)
if (
boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value &&
nb_cores > 1 && faces(tm1).size() >= min_nb_faces_to_split) {
// std::cout << "* executing parallel version " << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "* executing parallel version " << std::endl;
#endif
Bounded_error_distance_computation<TMF, TM2, VPM1, VPM2, TMF_tree, TM2_tree, Kernel> bedc(
tm1_parts, tm2, error_bound, vpm1, vpm2,
infinity_value, initial_bound, tm1_trees, tm2_tree);
@ -2020,21 +2091,24 @@ double bounded_error_one_sided_Hausdorff_impl(
} else // sequential version
#endif // defined(CGAL_LINKED_WITH_TBB) && defined(CGAL_METIS_ENABLED)
{
// std::cout << "* executing sequential version " << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout << "* executing sequential version " << std::endl;
#endif
hdist = bounded_error_Hausdorff_impl<Kernel>(
tm1, tm2, error_bound, vpm1, vpm2,
infinity_value, initial_bound, distance_bound,
tm1_tree, tm2_tree, out);
}
#ifdef CGAL_HAUSDORFF_DEBUG
timer.stop();
// std::cout << "* computation time (sec.) " << timer.time() << std::endl;
std::cout << "* computation time (sec.) " << timer.time() << std::endl;
#endif
CGAL_assertion(hdist >= 0.0);
return hdist;
}
// TODO: should we add a parallel version here?
template< class Concurrency_tag,
class Kernel,
class TriangleMesh1,
@ -2064,14 +2138,6 @@ double bounded_error_symmetric_Hausdorff_impl(
"Parallel_tag is enabled but at least TBB or METIS is unavailable.");
#endif
// Naive version.
// TODO: we can use it for parallel version to simplify our life.
// const double hdist1 = bounded_error_one_sided_Hausdorff_impl<Concurrency_tag, Kernel>(
// tm1, tm2, error_bound, distance_bound, compare_meshes, vpm1, vpm2, np1, np2, out1);
// const double hdist2 = bounded_error_one_sided_Hausdorff_impl<Concurrency_tag, Kernel>(
// tm2, tm1, error_bound, distance_bound, compare_meshes, vpm2, vpm1, np1, np2, out2);
// return (CGAL::max)(hdist1, hdist2);
// Optimized version.
// -- We compare meshes only if it is required.
// -- We first build trees and rebuild them only if it is required.
@ -2090,7 +2156,6 @@ double bounded_error_symmetric_Hausdorff_impl(
using Face_handle_1 = typename boost::graph_traits<TriangleMesh1>::face_descriptor;
using Face_handle_2 = typename boost::graph_traits<TriangleMesh2>::face_descriptor;
std::cout.precision(20);
std::vector<Face_handle_1> tm1_only;
std::vector<Face_handle_2> tm2_only;
@ -2104,7 +2169,10 @@ double bounded_error_symmetric_Hausdorff_impl(
tm1_tree, tm2_tree, tm1_only, tm2_only);
if (infinity_value < FT(0)) {
// std::cout << "* culling rate: 100%" << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout.precision(17);
std::cout << "* culling rate: 100%" << std::endl;
#endif
const auto face1 = *(faces(tm1).begin());
const auto face2 = *(faces(tm2).begin());
*out1++ = std::make_pair(face1, face2);
@ -2470,7 +2538,7 @@ template< class Concurrency_tag,
class TriangleMesh2,
class NamedParameters1,
class NamedParameters2 >
bool is_further_than(
bool is_Hausdorff_distance_larger(
const TriangleMesh1& tm1,
const TriangleMesh2& tm2,
const double distance_bound,
@ -2516,8 +2584,11 @@ bool is_further_than(
}
CGAL_assertion(hdist >= 0.0);
// std::cout << "- fin distance: " << hdist << std::endl;
// std::cout << "- max distance: " << distance_bound << std::endl;
#ifdef CGAL_HAUSDORFF_DEBUG
std::cout.precision(17);
std::cout << "- fin distance: " << hdist << std::endl;
std::cout << "- max distance: " << distance_bound << std::endl;
#endif
return hdist > distance_bound;
}
@ -2525,27 +2596,27 @@ template< class Concurrency_tag,
class TriangleMesh1,
class TriangleMesh2,
class NamedParameters1 >
double is_further_than(
double is_Hausdorff_distance_larger(
const TriangleMesh1& tm1,
const TriangleMesh2& tm2,
const double max_distance,
const double error_bound,
const NamedParameters1& np1)
{
return is_further_than<Concurrency_tag>(
return is_Hausdorff_distance_larger<Concurrency_tag>(
tm1, tm2, max_distance, error_bound, np1, parameters::all_default());
}
template< class Concurrency_tag,
class TriangleMesh1,
class TriangleMesh2 >
double is_further_than(
double is_Hausdorff_distance_larger(
const TriangleMesh1& tm1,
const TriangleMesh2& tm2,
const double max_distance = 1.0,
const double error_bound = 0.0001)
{
return is_further_than<Concurrency_tag>(
return is_Hausdorff_distance_larger<Concurrency_tag>(
tm1, tm2, max_distance, error_bound, parameters::all_default());
}

View File

@ -1,6 +1,9 @@
// Use it to include parallel computations in the bounded error Hausdorff distance.
// #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
#include <CGAL/Bbox_3.h>
#include <CGAL/Real_timer.h>
#include <CGAL/Simple_cartesian.h>
@ -10,6 +13,7 @@
#include <CGAL/IO/PLY.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
@ -28,6 +32,7 @@ using Triangle_3 = typename Kernel::Triangle_3;
using TAG = CGAL::Sequential_tag;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
using Polyhedron = CGAL::Polyhedron_3<Kernel>;
using Affine_transformation_3 = CGAL::Aff_transformation_3<Kernel>;
using Timer = CGAL::Real_timer;
@ -76,23 +81,26 @@ struct Bounded_error_hd_wrapper {
}
};
void get_mesh(const std::string filepath, Surface_mesh& mesh) {
template<typename PolygonMesh>
void get_mesh(const std::string filepath, PolygonMesh& mesh) {
mesh.clear();
std::ifstream input(filepath);
input >> mesh;
std::cout << "* getting mesh with " << mesh.number_of_faces() << " faces" << std::endl;
std::cout << "* getting mesh with " << faces(mesh).size() << " faces" << std::endl;
}
template<typename PolygonMesh1, typename PolygonMesh2>
void get_meshes(
const std::string filepath1, const std::string filepath2,
Surface_mesh& mesh1, Surface_mesh& mesh2) {
PolygonMesh1& mesh1, PolygonMesh2& mesh2) {
get_mesh(filepath1, mesh1);
get_mesh(filepath2, mesh2);
}
void save_mesh(const Surface_mesh& mesh, const std::string filepath) {
template<typename PolygonMesh>
void save_mesh(const PolygonMesh& mesh, const std::string filepath) {
if (!CGAL::write_PLY(filepath + ".ply", mesh)) {
std::cerr << "ERROR: cannot save this file: " << filepath << std::endl;
@ -100,6 +108,176 @@ void save_mesh(const Surface_mesh& mesh, const std::string filepath) {
}
}
// An easy example of a tetrahedron and its remeshed version.
void remeshing_tetrahedon_example(
const double error_bound, const bool save = false) {
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);
mesh2 = mesh1;
using edge_descriptor = typename boost::graph_traits<Surface_mesh>::edge_descriptor;
Surface_mesh::Property_map<edge_descriptor, bool> is_constrained_map =
mesh2.add_property_map<edge_descriptor, bool>("e:is_constrained", true).first;
const double target_edge_length = 0.05;
PMP::isotropic_remeshing(
mesh2.faces(), target_edge_length, mesh2,
PMP::parameters::edge_is_constrained_map(is_constrained_map));
if (save) save_mesh(mesh1, "mesh1");
if (save) save_mesh(mesh2, "mesh2");
timer.reset();
timer.start();
const double hdist = PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound);
std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl;
timer.stop();
std::cout << "* processing time: " << timer.time() << " s." << std::endl;
assert(hdist == error_bound);
}
// 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) {
Timer timer;
Surface_mesh mesh1, mesh2;
std::cout << std::endl << "* (E2) interior Triangle example:" << std::endl;
mesh1.add_vertex(Point_3(-1, 1, 1));
mesh1.add_vertex(Point_3( 0, -1, 1));
mesh1.add_vertex(Point_3( 1, 1, 1));
mesh1.add_face(mesh1.vertices());
auto v0 = mesh2.add_vertex(Point_3(-1.0, 1, 0));
auto v1 = mesh2.add_vertex(Point_3( 0.0, -1, 0));
auto v2 = mesh2.add_vertex(Point_3( 1.0, 1, 0));
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);
if (save) save_mesh(mesh1, "mesh1");
if (save) save_mesh(mesh2, "mesh2");
timer.reset();
timer.start();
const double hdist = PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound);
std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl;
timer.stop();
std::cout << "* processing time: " << timer.time() << " s." << std::endl;
assert(hdist >= 1.0);
}
// 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;
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;
if (save) save_mesh(mesh2, "mesh2");
timer.reset();
timer.start();
const double hdist = PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound);
std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl;
timer.stop();
std::cout << "* processing time: " << timer.time() << " s." << std::endl;
assert(hdist > 0.0);
}
// 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;
Surface_mesh mesh1;
Polyhedron mesh2;
get_mesh(filepath, mesh1);
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;
timer.reset();
timer.start();
const double hdist1 = PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound);
const double hdist2 = PMP::bounded_error_Hausdorff_distance<TAG>(mesh2, mesh1, error_bound);
std::cout << "* bounded-error Hausdorff distance 1->2: " << hdist1 << std::endl;
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(hdist2 > hdist1);
const double hdist3 = PMP::bounded_error_symmetric_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound);
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) {
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 FT t = FT(1) / FT(100);
if (save) save_mesh(mesh2, "mesh-0");
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);
timer.reset();
timer.start();
const double hdist = PMP::bounded_error_Hausdorff_distance<TAG>(mesh1, mesh2, error_bound);
std::cout << "* position: " << i << std::endl;
std::cout << "* bounded-error Hausdorff distance: " << hdist << std::endl;
timer.stop();
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;
}
}
template<class FunctionWrapper>
void test_0(const FunctionWrapper& functor, const bool save = false) {
@ -121,14 +299,16 @@ void test_0(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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);
}
template<class FunctionWrapper>
@ -155,14 +335,16 @@ void test_1(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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);
}
template<class FunctionWrapper>
@ -189,14 +371,16 @@ void test_2(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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);
}
template<class FunctionWrapper>
@ -224,14 +408,16 @@ void test_3(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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;
assert(CGAL::abs(dista - CGAL::sqrt(2.0)) < 1e-5); // error bound is 1e-4
assert(distb == 2.0);
assert(distc == naive);
}
template<class FunctionWrapper>
@ -258,14 +444,17 @@ void test_4(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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;
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);
}
template<class FunctionWrapper>
@ -293,14 +482,16 @@ void test_5(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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;
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);
}
template<class FunctionWrapper>
@ -333,14 +524,16 @@ void test_6(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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;
assert(dista == 0.0);
assert(CGAL::abs(distb - 0.707106) < 1e-5); // error bound is 1e-4
assert(distc == naive);
}
template<class FunctionWrapper>
@ -373,14 +566,16 @@ void test_7(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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;
assert(dista == 0.5);
assert(CGAL::abs(distb - 0.866025) < 1e-5); // error bound is 1e-4
assert(distc == naive);
}
template<class FunctionWrapper>
@ -413,14 +608,16 @@ void test_8(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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);
}
template<class FunctionWrapper>
@ -460,14 +657,16 @@ void test_9(const FunctionWrapper& functor, const bool save = false) {
const double dista = functor(mesh1, mesh2);
const double distb = functor(mesh2, mesh1);
const double naive = (CGAL::max)(dista, distb);
const double distc = functor.symmetric(mesh1, mesh2);
assert(distc == naive);
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);
}
template<class FunctionWrapper>
@ -527,6 +726,9 @@ void test_one_versus_another(
const double distb1 = functor2(mesh1, mesh2);
std::cout << "* Hausdorff distance1: " << dista1 << std::endl;
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
}
template<
@ -563,6 +765,9 @@ void test_real_meshes(
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<class FunctionWrapper>
@ -598,6 +803,10 @@ void test_timings(const std::string filepath, const FunctionWrapper& functor) {
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);
assert(timeab < timea + timeb);
PMP::transform(Affine_transformation_3(CGAL::Translation(),
Vector_3(FT(0), FT(0), FT(1))), mesh2);
@ -624,9 +833,16 @@ void test_timings(const std::string filepath, const FunctionWrapper& functor) {
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);
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<class FunctionWrapper>
@ -672,10 +888,12 @@ void test_bunny(
times.push_back(timer.time());
std::cout << "* distance / Hausdorff / time (sec.) : " <<
distance << " / " << distc << " / " << times.back() << std::endl;
assert(distc > 0.0);
} else {
// t is the step where n is the number of steps.
double curr_dist = 1.0;
const FT t = FT(1) / static_cast<FT>(n);
for (int k = n; k >= 0; --k) {
auto mesh = mesh2;
@ -693,6 +911,8 @@ void test_bunny(
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;
}
}
@ -710,6 +930,10 @@ void test_bunny(
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);
}
Triangle_3 get_triangle(const int index, const Surface_mesh& mesh) {
@ -721,6 +945,7 @@ Triangle_3 get_triangle(const int index, const Surface_mesh& mesh) {
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;
@ -744,6 +969,9 @@ void compute_realizing_triangles(
std::cout << "* Hausdorff: " << hdist << std::endl;
std::cout << "* f1 / f2: " << f1 << " / " << f2 << std::endl;
assert(f1 == 0);
assert(f2 == 0 || f2 == 161);
if (f1 != -1 && save) {
const auto triangle = get_triangle(f1, mesh1);
Surface_mesh sm1;
@ -864,6 +1092,10 @@ 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)
@ -880,7 +1112,7 @@ void test_early_quit(const std::string filepath) {
std::cout << " ---- distance 0.0 = 0.0 ---- " << std::endl;
timer.reset();
timer.start();
assert(!PMP::is_further_than<TAG>(mesh1, mesh2, 0.0));
assert(!PMP::is_Hausdorff_distance_larger<TAG>(mesh1, mesh2, 0.0));
timer.stop();
const double timea = timer.time();
@ -890,25 +1122,25 @@ void test_early_quit(const std::string filepath) {
std::cout << " ---- distance 0.5 < 1.0 ---- " << std::endl;
timer.reset();
timer.start();
assert(!PMP::is_further_than<TAG>(mesh1, mesh2, 1.0));
assert(!PMP::is_Hausdorff_distance_larger<TAG>(mesh1, mesh2, 1.0));
timer.stop();
const double timeb = timer.time();
std::cout << " ---- distance 0.5 < 0.6 ---- " << std::endl;
timer.reset();
timer.start();
assert(!PMP::is_further_than<TAG>(mesh1, mesh2, 0.6));
assert(!PMP::is_Hausdorff_distance_larger<TAG>(mesh1, mesh2, 0.6));
timer.stop();
const double timec = timer.time();
std::cout << " ---- distance 0.5 > 0.4 ---- " << std::endl;
timer.reset();
timer.start();
assert(PMP::is_further_than<TAG>(mesh1, mesh2, 0.4));
assert(PMP::is_Hausdorff_distance_larger<TAG>(mesh1, mesh2, 0.4));
timer.stop();
const double timed = timer.time();
std::cout << " ---- distance 0.5 > 0.0 ---- " << std::endl;
timer.reset();
timer.start();
assert(PMP::is_further_than<TAG>(mesh1, mesh2, 0.0));
assert(PMP::is_Hausdorff_distance_larger<TAG>(mesh1, mesh2, 0.0));
timer.stop();
const double timee = timer.time();
@ -917,6 +1149,21 @@ void test_early_quit(const std::string filepath) {
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);
}
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);
perturbing_polyhedron_mesh_example(filepath, error_bound);
moving_surface_mesh_example(filepath, filepath, 5, error_bound);
}
int main(int argc, char** argv) {
@ -929,6 +1176,7 @@ int main(int argc, char** argv) {
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] : "data/blobby.off");
run_examples(error_bound, filepath);
// ------------------------------------------------------------------------ //
// Tests.