Merge pull request #6978 from janetournois/Mesh_3-detect_cc_in_labeled_images-GF

Mesh 3 for labeled images - avoid vertex clusters on surfaces

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<class C3T3, class MeshDomain, class MeshCriteria>
@@ -114,7 +127,13 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3,
                                              image.vy()),
                                   image.vz());
 
-  typedef std::vector<std::pair<Bare_point, std::size_t> > Seeds;
+  struct Seed {
+    std::size_t i, j, k;
+    std::size_t radius;
+  };
+  using Seeds = std::vector<Seed>;
+  using Subdomain = typename Mesh_domain::Subdomain;
+
   Seeds seeds;
   Get_point<Bare_point> get_point(&image);
   std::cout << "Searching for connected components..." << std::endl;
@@ -122,20 +141,34 @@ 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 Bare_point seed_point = get_point(seed.i, seed.j, seed.k);
+    Cell_handle seed_cell = tr.locate(cwp(seed_point));
+
+    const Subdomain seed_label
+      = domain.is_in_domain_object()(seed_point);
+    const Subdomain seed_cell_label
+      = (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()));
+
+    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
+
+    const double radius = double(seed.radius + 1)* max_v;
     CGAL::Random_points_on_sphere_3<Bare_point> points_on_sphere_3(radius);
     typename Mesh_domain::Construct_intersection construct_intersection =
       domain.construct_intersection_object();
 
     std::vector<Vector_3> directions;
-    if(it->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));
@@ -153,9 +186,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);

Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h

@@ -34,14 +34,12 @@
 template <typename PointsOutputIterator,
           typename DomainsOutputIterator,
           typename TransformOperator,
-          typename Construct_point,
           typename Image_word_type>
 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,8 +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),
-                                       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")