From 26d1a68878e1bf0f6b17968a42018679380eff91 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 17 Jul 2023 17:31:48 +0200 Subject: [PATCH 1/7] Add a tentative new initialization scheme for Mesh_3's polyhedral domains --- .../include/CGAL/Polyhedral_mesh_domain_3.h | 118 +++++++++++++----- 1 file changed, 88 insertions(+), 30 deletions(-) diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index 2d9bf625f57..75ab6ec89d7 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -687,54 +687,112 @@ Polyhedral_mesh_domain_3:: Construct_initial_points::operator()(OutputIterator pts, const int n) const { - typename IGT::Construct_ray_3 ray = IGT().construct_ray_3_object(); - typename IGT::Construct_vector_3 vector = IGT().construct_vector_3_object(); + typedef typename std::pair Point_with_index; - const Bounding_box bbox = r_domain_.tree_.bbox(); - const Point_3 center( FT( (bbox.xmin() + bbox.xmax()) / 2), - FT( (bbox.ymin() + bbox.ymax()) / 2), - FT( (bbox.zmin() + bbox.zmax()) / 2) ); - - CGAL::Random& rng = *(r_domain_.p_rng_ != 0 ? - r_domain_.p_rng_ : - new Random(0)); - Random_points_on_sphere_3 random_point(1., rng); - - int i = n; # ifdef CGAL_MESH_3_VERBOSE std::cerr << "construct initial points:" << std::endl; # endif - // Point construction by ray shooting from the center of the enclosing bbox - while ( i > 0 ) + + CGAL_precondition(!r_domain_.tree_.empty()); + + const Bounding_box bbox = r_domain_.tree_.bbox(); + + // The initialization proceeds as follows: + // - Generate a grid of n' points with n' > n, to get good candidates + // - Project them onto the surface + // - Grid snap to avoid very close points (ideally that would depend on the sizing field) + // - Keep the n closest points + + const int n_prime = 10 * n; + const int n_dir = std::ceil(std::cbrt(n_prime)); + + std::vector grid_points; + grid_points.reserve(n_dir * n_dir * n_dir); + + // generate points in the bounding bbox + for(int i=0; i projected_points; + + // @todo trivial to parallelize + for(const Point_3& gp : grid_points) { - const Ray_3 ray_shot = ray(center, vector(CGAL::ORIGIN,*random_point)); + typename AABB_tree_::Point_and_primitive_id cp_and_p = r_domain_.tree_.closest_point_and_primitive(gp); + auto index = r_domain_.index_from_surface_patch_index( + r_domain_.make_surface_index(cp_and_p.second /*primitive ID*/)); + projected_points.emplace(cp_and_p.first, index); + } -#ifdef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3 - Intersection intersection = r_domain_.construct_intersection_object()(ray_shot); - if(std::get<2>(intersection) != 0) { -#else - if(r_domain_.do_intersect_surface_object()(ray_shot)) { - Intersection intersection = r_domain_.construct_intersection_object()(ray_shot); -#endif - *pts++ = std::make_pair(std::get<0>(intersection), - std::get<1>(intersection)); + CGAL_warning(projected_points.size() > n); - --i; + // n farthest points + std::vector initial_points; + + // start with a random point + CGAL::Random& rng = *(r_domain_.p_rng_ != 0 ? r_domain_.p_rng_ : new Random(0)); + + auto it = projected_points.begin(); + std::advance(it, rng.get_int(0, projected_points.size())); + initial_points.push_back(*it); + + // for each next initial point, take the point farthest from all existing initial points + while(initial_points.size() < static_cast(n)) + { + const Point_with_index* next_ip = nullptr; + FT farthest_distance = FT(0); + + for(const Point_with_index& pp : projected_points) + { + FT min_distance = FT(std::numeric_limits::max()); + for(const Point_with_index& ip : initial_points) + min_distance = (std::min)(min_distance, CGAL::squared_distance(pp.first, ip.first)); + + if(min_distance > farthest_distance) + { + farthest_distance = min_distance; + next_ip = &pp; + } + } + + if(is_zero(farthest_distance)) + break; + + initial_points.emplace_back(*next_ip); #ifdef CGAL_MESH_3_VERBOSE std::cerr << boost::format("\r \r" "%1%/%2% initial point(s) found...") - % (n - i) + % initial_points.size() % n; # endif - } - ++random_point; } #ifdef CGAL_MESH_3_VERBOSE std::cerr << std::endl; + + for(const Point_with_index& ip : initial_points) + std::cerr << "\t" << ip.first << std::endl; + + if(initial_points.size() != n) + std::cerr << "Warning: failed to construct " << n << " initial points" << std::endl; #endif - if(r_domain_.p_rng_ == 0) delete &rng; + + for(const Point_with_index& ip : initial_points) + *pts++ = ip; + + if(r_domain_.p_rng_ == nullptr) + delete &rng; + return pts; } From 370070ea2c9545cc7c399b158d2ffc12cd8f06bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 18 Jul 2023 11:39:18 +0200 Subject: [PATCH 2/7] Add a bit of randomness in the grid point positions --- .../include/CGAL/Polyhedral_mesh_domain_3.h | 21 ++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index 75ab6ec89d7..3e5a4d7de1f 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -697,6 +697,8 @@ Construct_initial_points::operator()(OutputIterator pts, const Bounding_box bbox = r_domain_.tree_.bbox(); + CGAL::Random& rng = *(r_domain_.p_rng_ != 0 ? r_domain_.p_rng_ : new Random(0)); + // The initialization proceeds as follows: // - Generate a grid of n' points with n' > n, to get good candidates // - Project them onto the surface @@ -714,9 +716,20 @@ Construct_initial_points::operator()(OutputIterator pts, for(int j=0; j initial_points; // start with a random point - CGAL::Random& rng = *(r_domain_.p_rng_ != 0 ? r_domain_.p_rng_ : new Random(0)); - auto it = projected_points.begin(); std::advance(it, rng.get_int(0, projected_points.size())); initial_points.push_back(*it); From fb8fc437beccf906c65a391eb46a9ee45ad6c2dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Fri, 4 Aug 2023 14:45:53 +0200 Subject: [PATCH 3/7] Fix comment --- Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index b5bd7a4a76c..ca9fd32aa4b 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -704,7 +704,7 @@ Construct_initial_points::operator()(OutputIterator pts, // - Generate a grid of n' points with n' > n, to get good candidates // - Project them onto the surface // - Grid snap to avoid very close points (ideally that would depend on the sizing field) - // - Keep the n closest points + // - Keep the n farthest points const int n_prime = 10 * n; const int n_dir = std::ceil(std::cbrt(n_prime)); From de957d430fd9d48c365b263a81885c26c5189808 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 5 Sep 2023 14:29:21 +0200 Subject: [PATCH 4/7] Fix warning --- Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index ca9fd32aa4b..2371dd94d75 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -795,7 +795,7 @@ Construct_initial_points::operator()(OutputIterator pts, for(const Point_with_index& ip : initial_points) std::cerr << "\t" << ip.first << std::endl; - if(initial_points.size() != n) + if(initial_points.size() != static_cast(n)) std::cerr << "Warning: failed to construct " << n << " initial points" << std::endl; #endif From ed1db1140c1338b116c789bce2196d953d9faad8 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 17 Oct 2023 17:32:51 +0200 Subject: [PATCH 5/7] protect max with () for msvc --- Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index 2371dd94d75..9eb2fad6141 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -765,7 +765,7 @@ Construct_initial_points::operator()(OutputIterator pts, for(const Point_with_index& pp : projected_points) { - FT min_distance = FT(std::numeric_limits::max()); + FT min_distance = FT((std::numeric_limits::max)()); for(const Point_with_index& ip : initial_points) min_distance = (std::min)(min_distance, CGAL::squared_distance(pp.first, ip.first)); From 27eaea6c4a01a9f51260b3f5740db6af7abaf8dc Mon Sep 17 00:00:00 2001 From: Sebastien Loriot Date: Wed, 18 Oct 2023 11:54:20 +0200 Subject: [PATCH 6/7] fix warning --- Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index 9eb2fad6141..2e1e7b2e124 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -747,7 +747,7 @@ Construct_initial_points::operator()(OutputIterator pts, projected_points.emplace(cp_and_p.first, index); } - CGAL_warning(projected_points.size() > n); + CGAL_warning(projected_points.size() > static_cast(n)); // n farthest points std::vector initial_points; From 02214aedae2b883fe4b84b480c4dc5d7806c78bc Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 19 Oct 2023 08:43:57 +0100 Subject: [PATCH 7/7] Add static_cast to quiet VC++ --- Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index 2e1e7b2e124..fadc59a10c0 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -707,7 +707,7 @@ Construct_initial_points::operator()(OutputIterator pts, // - Keep the n farthest points const int n_prime = 10 * n; - const int n_dir = std::ceil(std::cbrt(n_prime)); + const int n_dir = static_cast(std::ceil(std::cbrt(n_prime))); std::vector grid_points; grid_points.reserve(n_dir * n_dir * n_dir); @@ -754,7 +754,7 @@ Construct_initial_points::operator()(OutputIterator pts, // start with a random point auto it = projected_points.begin(); - std::advance(it, rng.get_int(0, projected_points.size())); + std::advance(it, rng.get_int(0, static_cast(projected_points.size()))); initial_points.push_back(*it); // for each next initial point, take the point farthest from all existing initial points