From db6b0519737a6e685769f837f09abbd00469ae72 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 20 Oct 2022 12:07:27 +0200 Subject: [PATCH 1/8] on image boundaries, construct point away from the surface this avoids creating point clusters on surfaces --- ...tialize_triangulation_from_labeled_image.h | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index 74bb943f76e..c7476e96c67 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -33,6 +33,7 @@ struct Get_point { const double vx, vy, vz; const double tx, ty, tz; + const std::size_t xdim, ydim, zdim; Get_point(const CGAL::Image_3* image) : vx(image->vx()) , vy(image->vy()) @@ -40,15 +41,27 @@ struct Get_point , tx(image->tx()) , ty(image->ty()) , tz(image->tz()) + , xdim(image->xdim()) + , ydim(image->ydim()) + , zdim(image->zdim()) {} Point operator()(const std::size_t i, const std::size_t j, const std::size_t k) const { - return Point(double(i) * vx + tx, - double(j) * vy + ty, - double(k) * vz + tz); + double x = double(i) * vx + tx; + double y = double(j) * vy + ty; + double z = double(k) * vz + tz; + + if (i == 0) x += 1. / 6. * vx; + else if (i == xdim - 1) x -= 1. / 6. * vx; + if (j == 0) y += 1. / 6. * vy; + else if (j == ydim - 1) y -= 1. / 6. * vy; + if (k == 0) z += 1. / 6. * vz; + else if (k == zdim - 1) z -= 1. / 6. * vz; + + return Point(x, y, z); } }; template From 28a69460706803ababf0059057122e96f759de40 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 21 Oct 2022 15:05:30 +0200 Subject: [PATCH 2/8] do not run random shooting in connected components that already are represented during feature protection, most of the connected components, in particular on the bbox boundary, are already represented and ready to start/hint Delaunay refinement if the chosen seed lies in a cell that already belongs (wrt Is_in_domain(cc)) to the right connected component, random shooting is canceled and the loop continues to next seed --- ...tialize_triangulation_from_labeled_image.h | 38 +++++++++++++++---- ...or_connected_components_in_labeled_image.h | 4 +- 2 files changed, 31 insertions(+), 11 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index c7476e96c67..ba3cad5f699 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -127,7 +128,11 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, image.vy()), image.vz()); - typedef std::vector > Seeds; + using Point_indices = std::tuple; + using Seed = std::pair; + using Seeds = std::vector; + using Subdomain_index = typename Mesh_domain::Subdomain_index; + Seeds seeds; Get_point get_point(&image); std::cout << "Searching for connected components..." << std::endl; @@ -135,20 +140,36 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, std::back_inserter(seeds), CGAL::Emptyset_iterator(), transform, - get_point, Image_word_type()); std::cout << " " << seeds.size() << " components were found." << std::endl; std::cout << "Construct initial points..." << std::endl; - for(typename Seeds::const_iterator it = seeds.begin(), end = seeds.end(); - it != end; ++it) + for(const Seed seed : seeds) { - const double radius = double(it->second + 1)* max_v; + const std::size_t i = get<0>(seed.first); + const std::size_t j = get<1>(seed.first); + const std::size_t k = get<2>(seed.first); + + const Bare_point seed_point = get_point(i, j, k); + Cell_handle seed_cell = tr.locate(cwp(seed_point)); + + const boost::optional seed_label + = domain.is_in_domain_object()(seed_point); + const boost::optional seed_cell_label + = domain.is_in_domain_object()( + seed_cell->weighted_circumcenter(tr.geom_traits())); + + if ( seed_label != boost::none + && seed_cell_label != boost::none + && seed_label.get() == seed_cell_label.get()) + continue; //this means the connected component has already been initialized + + const double radius = double(seed.second + 1)* max_v; CGAL::Random_points_on_sphere_3 points_on_sphere_3(radius); typename Mesh_domain::Construct_intersection construct_intersection = domain.construct_intersection_object(); std::vector directions; - if(it->second < 2) { + if(seed.second < 2) { // shoot in six directions directions.push_back(Vector_3(-radius, 0, 0)); directions.push_back(Vector_3(+radius, 0, 0)); @@ -166,9 +187,10 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, for(const Vector_3& v : directions) { - const Bare_point test = it->first + v; + const Bare_point test = seed_point + v; + const typename Mesh_domain::Intersection intersect = - construct_intersection(Segment_3(it->first, test)); + construct_intersection(Segment_3(seed_point, test)); if (std::get<2>(intersect) != 0) { const Bare_point& bpi = std::get<0>(intersect); diff --git a/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h index 56c5d7bc6d5..23ca8cb5cc3 100644 --- a/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h @@ -34,14 +34,12 @@ template void search_for_connected_components_in_labeled_image(const CGAL::Image_3& image, PointsOutputIterator it, DomainsOutputIterator dom_it, TransformOperator transform, - Construct_point point, Image_word_type) { const std::size_t nx = image.xdim(); @@ -210,7 +208,7 @@ search_for_connected_components_in_labeled_image(const CGAL::Image_3& image, { // if(nb_voxels >= 100) { - *it++ = std::make_pair(point(i, j, k), + *it++ = std::make_pair(std::make_tuple(i, j, k), depth+1); #if CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE > 1 std::cerr << boost::format("Found seed %5%, which is voxel " From 8b76b55b82d66b3b65b1048f2af738420a294224 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 21 Oct 2022 15:48:45 +0200 Subject: [PATCH 3/8] remove useless include --- .../CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index ba3cad5f699..e4be4b323f9 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -24,7 +24,6 @@ #include #include #include -#include #include #include From 7faa73bb27df4e0cee95183e942964e8c15af609 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Mon, 24 Oct 2022 12:47:04 +0200 Subject: [PATCH 4/8] locate of seed that is outside affine hull returns Cell_handle() accessing the circumcenter was then failing --- .../Mesh_3/initialize_triangulation_from_labeled_image.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index e4be4b323f9..491ade42e05 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -154,8 +154,10 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, const boost::optional seed_label = domain.is_in_domain_object()(seed_point); const boost::optional seed_cell_label - = domain.is_in_domain_object()( - seed_cell->weighted_circumcenter(tr.geom_traits())); + = (seed_cell == Cell_handle()) //seed_point is OUTSIDE_AFFINE_HULL + ? boost::none + : domain.is_in_domain_object()( + seed_cell->weighted_circumcenter(tr.geom_traits())); if ( seed_label != boost::none && seed_cell_label != boost::none From a8877c5a441f309ef842f58b8625bf97fea55f6e Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 25 Oct 2022 17:24:57 +0200 Subject: [PATCH 5/8] add struct Seed for better readability and avoid using boost::none to be ready for std::optional (C++17) --- ...tialize_triangulation_from_labeled_image.h | 29 +++++++++---------- ...or_connected_components_in_labeled_image.h | 3 +- 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index 491ade42e05..5c7314048f4 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -127,10 +127,13 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, image.vy()), image.vz()); - using Point_indices = std::tuple; - using Seed = std::pair; + struct Seed { + std::size_t i, j, k; + std::size_t radius; + }; using Seeds = std::vector; using Subdomain_index = typename Mesh_domain::Subdomain_index; + using Subdomain = typename Mesh_domain::Subdomain; Seeds seeds; Get_point get_point(&image); @@ -144,33 +147,29 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, std::cout << "Construct initial points..." << std::endl; for(const Seed seed : seeds) { - const std::size_t i = get<0>(seed.first); - const std::size_t j = get<1>(seed.first); - const std::size_t k = get<2>(seed.first); - - const Bare_point seed_point = get_point(i, j, k); + const Bare_point seed_point = get_point(seed.i, seed.j, seed.k); Cell_handle seed_cell = tr.locate(cwp(seed_point)); - const boost::optional seed_label + const Subdomain seed_label = domain.is_in_domain_object()(seed_point); - const boost::optional seed_cell_label + const Subdomain seed_cell_label = (seed_cell == Cell_handle()) //seed_point is OUTSIDE_AFFINE_HULL - ? boost::none + ? Subdomain() : domain.is_in_domain_object()( seed_cell->weighted_circumcenter(tr.geom_traits())); - if ( seed_label != boost::none - && seed_cell_label != boost::none - && seed_label.get() == seed_cell_label.get()) + if ( seed_label.has_value() + && seed_cell_label.has_value() + && *seed_label == *seed_cell_label) continue; //this means the connected component has already been initialized - const double radius = double(seed.second + 1)* max_v; + const double radius = double(seed.radius + 1)* max_v; CGAL::Random_points_on_sphere_3 points_on_sphere_3(radius); typename Mesh_domain::Construct_intersection construct_intersection = domain.construct_intersection_object(); std::vector directions; - if(seed.second < 2) { + if(seed.radius < 2) { // shoot in six directions directions.push_back(Vector_3(-radius, 0, 0)); directions.push_back(Vector_3(+radius, 0, 0)); diff --git a/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h index 23ca8cb5cc3..26849514840 100644 --- a/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h @@ -208,8 +208,7 @@ search_for_connected_components_in_labeled_image(const CGAL::Image_3& image, { // if(nb_voxels >= 100) { - *it++ = std::make_pair(std::make_tuple(i, j, k), - depth+1); + *it++ = { i, j, k, std::size_t(depth + 1) }; #if CGAL_MESH_3_SEARCH_FOR_CONNECTED_COMPONENTS_IN_LABELED_IMAGE_VERBOSE > 1 std::cerr << boost::format("Found seed %5%, which is voxel " "(%1%, %2%, %3%), value=%4%\n") From 168825a7f116f2ea2db43b98e34e31303e25307f Mon Sep 17 00:00:00 2001 From: Sebastien Loriot Date: Wed, 2 Nov 2022 07:49:59 +0100 Subject: [PATCH 6/8] fix warning --- .../CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index 5c7314048f4..6ea46d43663 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -132,7 +132,6 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, std::size_t radius; }; using Seeds = std::vector; - using Subdomain_index = typename Mesh_domain::Subdomain_index; using Subdomain = typename Mesh_domain::Subdomain; Seeds seeds; From 7f70d48ab9d1797a6565a5be67bbafd04320bede Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 3 Nov 2022 12:52:05 +0100 Subject: [PATCH 7/8] seed_cell cannot be infinite --- .../CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index 6ea46d43663..1d6d26718d4 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -152,8 +152,8 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, const Subdomain seed_label = domain.is_in_domain_object()(seed_point); const Subdomain seed_cell_label - = (seed_cell == Cell_handle()) //seed_point is OUTSIDE_AFFINE_HULL - ? Subdomain() + = (seed_cell == Cell_handle() || tr.is_infinite(seed_cell)) + ? Subdomain() //seed_point is OUTSIDE_AFFINE_HULL : domain.is_in_domain_object()( seed_cell->weighted_circumcenter(tr.geom_traits())); From 3816d8c2af1836a7d2c5d378f6fe06ea82158e7b Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 3 Nov 2022 12:53:38 +0100 Subject: [PATCH 8/8] boost::optional has_value() is available only since boost 1.68 and CGAL supports boost >= 1.66 --- .../CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index 1d6d26718d4..5a7b52b6685 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -157,8 +157,8 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, : domain.is_in_domain_object()( seed_cell->weighted_circumcenter(tr.geom_traits())); - if ( seed_label.has_value() - && seed_cell_label.has_value() + if ( seed_label != boost::none + && seed_cell_label != boost::none && *seed_label == *seed_cell_label) continue; //this means the connected component has already been initialized