Merge pull request #5463 from janetournois/Mesh_3-example_initialization_cc_in_gray_images-jtournois

Mesh 3 - connected components initialization for 3D gray images

Mesh_3/doc/Mesh_3/Mesh_3.txt

@@ -573,13 +573,6 @@ constructor of the `Mesh_criteria` instance.
 Cut view of a 3D mesh produced from an implicit domain
 \cgalFigureEnd
 
-\subsubsection Mesh_33DDomainsGrayImageIsosurfaces 3D Domains Bounded by Isosurfaces in 3D Gray-Level Images
-
-The following example produces a 3D mesh for a domain whose boundary surface
-is the isosurface associated to an isovalue inside the input gray-level
-3D image. In the distribution you can also find the example \ref Mesh_3/mesh_3D_gray_vtk_image.cpp which can deal with DICOM files as input.
-
-\cgalExample{Mesh_3/mesh_3D_gray_image.cpp}
 
 \subsection Mesh_3MeshingMultipleDomains Meshing Multiple Domains
 
@@ -682,8 +675,18 @@ This allows to remesh a surface, and is equivalent to the function `make_surface
 View of a remeshed surface. (Left) input mesh (Right) output mesh. Code from subsection \ref Mesh_3RemeshingPolyhedralSurface generates the file.
 \cgalFigureEnd
 
-\subsection Mesh_3DomainsFromSegmented3DImages Domains From Segmented 3D Images
+\subsection Mesh_3DomainsFrom3DImages Domains From 3D Images
 
+\subsubsection Mesh_33DDomainsGrayImageIsosurfaces 3D Domains Bounded by Isosurfaces in 3D Gray-Level Images
+
+The following example produces a 3D mesh for a domain whose boundary surface
+is the isosurface associated to an isovalue inside the input gray-level
+3D image. In the distribution you can also find the example \ref Mesh_3/mesh_3D_gray_vtk_image.cpp which can deal with DICOM files as input.
+
+\cgalExample{Mesh_3/mesh_3D_gray_image.cpp}
+
+
+\subsubsection Mesh_3DomainsFromSegmented3DImages Domains From Segmented 3D Images
 \anchor Mesh_3_subsection_examples_3d_image
 The following code produces a 3D mesh from
 a 3D image. The image is a segmented medical image in which each
@@ -703,7 +706,7 @@ The resulting mesh is shown in \cgalFigureRef{figureliver_3d_image_mesh}.
 Cut view of a 3D mesh produced from a segmented liver image. Code from subsection \ref Mesh_3_subsection_examples_3d_image generates this file.
 \cgalFigureEnd
 
-\subsubsection Mesh_3DomainsFromSegmented3DImagesWithCustomInitialization Domains From Segmented 3D Images, with a Custom Initialization
+\subsubsection Mesh_3DomainsFrom3DImagesWithCustomInitialization Domains From 3D Images, with a Custom Initialization
 
 The example \ref Mesh_3/mesh_3D_image_with_custom_initialization.cpp is a modification
 of \ref Mesh_3/mesh_3D_image.cpp. The goal of that example is to show how
@@ -788,6 +791,25 @@ create a 3D image using the undocumented API of CGAL_ImageIO.
 The code of the function `%random_labeled_image()` is in the header file \ref
 Mesh_3/random_labeled_image.h\.
 
+
+The example \ref Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp is another
+custom initialization example, for meshing of 3D gray-level images. Similarly to
+the segmented image example above, the code consists in:
+  -# the creation of an empty `%c3t3` object,
+  -# a call to a non-documented function
+    `initialize_triangulation_from_gray_image()` that inserts points in
+    the triangulation,
+  -# then the call to `refine_mesh_3()`.
+
+\snippet Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp Meshing
+
+The code of the function `initialize_triangulation_from_gray_image()` is
+in the non-documented header \ref
+CGAL/Mesh_3/initialize_triangulation_from_gray_image.h\. As it is
+undocumented and may be removed or modified at any time, if you wish to
+use it then you should copy-paste it to your user code.
+
+
 \subsection Mesh_3UsingVariableSizingField Using Variable Sizing Field
 
 \subsubsection Mesh_3SizingFieldasanAnalyticalFunction Sizing Field as an Analytical Function

Mesh_3/doc/Mesh_3/examples.txt

@@ -1,9 +1,11 @@
 /*!
 \example Mesh_3/implicit_functions.cpp
 \example Mesh_3/mesh_3D_image.cpp
+\example Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp
 \example Mesh_3/mesh_3D_image_with_features.cpp
 \example Mesh_3/mesh_3D_image_with_custom_initialization.cpp
 \example Mesh_3/random_labeled_image.h
+\example CGAL/Mesh_3/initialize_triangulation_from_gray_image.h
 \example CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h
 \example Mesh_3/mesh_3D_image_variable_size.cpp
 \example Mesh_3/mesh_hybrid_mesh_domain.cpp

Mesh_3/examples/Mesh_3/CMakeLists.txt

@@ -149,6 +149,11 @@ if(TARGET CGAL::CGAL_ImageIO)
     target_link_libraries(mesh_3D_image_with_custom_initialization
                           PUBLIC CGAL::Eigen3_support)
 
+    create_single_source_cgal_program(
+      "mesh_3D_gray_image_with_custom_initialization.cpp")
+    target_link_libraries(mesh_3D_gray_image_with_custom_initialization
+                          PUBLIC CGAL::Eigen3_support)
+
     create_single_source_cgal_program("mesh_3D_image_variable_size.cpp")
     target_link_libraries(mesh_3D_image_variable_size
                           PUBLIC CGAL::Eigen3_support)
@@ -172,6 +177,7 @@ if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support)
     mesh_3D_image
     mesh_3D_image_variable_size
     mesh_3D_image_with_custom_initialization
+    mesh_3D_gray_image_with_custom_initialization
     mesh_3D_image_with_features
     mesh_implicit_domains
     mesh_implicit_sphere

Mesh_3/examples/Mesh_3/mesh_3D_gray_image_with_custom_initialization.cpp

@@ -0,0 +1,71 @@
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Mesh_3/initialize_triangulation_from_gray_image.h>
+
+#include <CGAL/Labeled_mesh_domain_3.h>
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/Image_3.h>
+#include <functional>
+
+typedef float Image_word_type;
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Labeled_mesh_domain_3<K> Mesh_domain;
+
+// Parallel tag
+#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Parallel_tag Concurrency_tag;
+#else
+typedef CGAL::Sequential_tag Concurrency_tag;
+#endif
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+int main(int argc, char* argv[])
+{
+  const char* fname = (argc > 1) ? argv[1] : "data/skull_2.9.inr";
+  /// [Load image]
+  CGAL::Image_3 image;
+  if (!image.read(fname)) {
+    std::cerr << "Error: Cannot read file " << fname << std::endl;
+    return EXIT_FAILURE;
+  }
+  /// [Domain creation]
+  Mesh_domain domain =
+    Mesh_domain::create_gray_image_mesh_domain(image, 2.9f, 0.f);
+  /// [Domain creation]
+
+  /// [Mesh criteria]
+  Mesh_criteria criteria(facet_angle = 30, facet_size = 6, facet_distance = 2,
+    cell_radius_edge_ratio = 3, cell_size = 8);
+
+  /// [Meshing]
+  C3t3 c3t3;
+  initialize_triangulation_from_gray_image(c3t3,
+                                           domain,
+                                           image,
+                                           criteria,
+                                           2.9f,//isolevel
+                                           Image_word_type(0));
+  CGAL::refine_mesh_3(c3t3, domain, criteria);
+  /// [Meshing]
+
+  /// Output
+  CGAL::dump_c3t3(c3t3, "out");
+
+  return 0;
+}

Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_gray_image.h

@@ -0,0 +1,50 @@
+// Copyright (c) 2015,2016 GeometryFactory
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Laurent Rineau, Jane Tournois
+
+#ifndef CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_GRAY_IMAGE_H
+#define CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_GRAY_IMAGE_H
+
+#include <CGAL/license/Mesh_3.h>
+
+#include <CGAL/Labeled_mesh_domain_3.h>
+#include <CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h>
+
+#include <CGAL/tags.h>
+
+template<class C3T3, class MeshDomain, class MeshCriteria,
+         typename FT,
+         typename Image_word_type,
+         typename Functor = CGAL::Null_functor>
+void initialize_triangulation_from_gray_image(C3T3& c3t3,
+         const MeshDomain& domain,
+         const CGAL::Image_3& image,
+         const MeshCriteria& criteria,
+         const FT& iso_value,
+         Image_word_type,
+         const Functor image_values_to_subdomain_indices = CGAL::Null_functor(),
+         bool protect_features = false)
+{
+  typedef typename CGAL::Default::Get<Functor, CGAL::Null_functor>::type Functor_;
+
+  using CGAL::Mesh_3::internal::Create_gray_image_values_to_subdomain_indices;
+  typedef Create_gray_image_values_to_subdomain_indices<Functor_> C_i_v_t_s_i;
+  typedef typename C_i_v_t_s_i::type Image_values_to_subdomain_indices;
+  Image_values_to_subdomain_indices transform_fct =
+    C_i_v_t_s_i()(image_values_to_subdomain_indices, iso_value);
+
+  initialize_triangulation_from_labeled_image(c3t3, domain, image, criteria,
+                                                                  Image_word_type(),
+                                                                  protect_features,
+                                                                  transform_fct);
+}
+
+#endif // CGAL_MESH_3_INITIALIZE_TRIANGULATION_FROM_GRAY_IMAGE_H

Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h

@@ -17,6 +17,7 @@
 
 #include <CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h>
 #include <CGAL/Mesh_3/squared_distance_Point_3_Triangle_3.h>
+#include <CGAL/Labeled_mesh_domain_3.h>
 #include <CGAL/make_mesh_3.h>
 
 #include <CGAL/enum.h>
@@ -31,17 +32,23 @@ template <typename Point>
 struct Get_point
 {
   const double vx, vy, vz;
+  const double tx, ty, tz;
   Get_point(const CGAL::Image_3* image)
     : vx(image->vx())
     , vy(image->vy())
     , vz(image->vz())
+    , tx(image->tx())
+    , ty(image->ty())
+    , tz(image->tz())
   {}
 
   Point operator()(const std::size_t i,
                    const std::size_t j,
                    const std::size_t k) const
   {
-    return Point(double(i) * vx, double(j) * vy, double(k) * vz);
+    return Point(double(i) * vx + tx,
+                 double(j) * vy + ty,
+                 double(k) * vz + tz);
   }
 };
 template<class C3T3, class MeshDomain, class MeshCriteria>
@@ -64,16 +71,16 @@ void init_tr_from_labeled_image_call_init_features(C3T3& c3t3,
             << " initial points on 1D-features" << std::endl;
 }
 
-
 template<class C3T3, class MeshDomain, class MeshCriteria,
-         typename Image_word_type>
+         typename Image_word_type,
+         typename TransformOperator = CGAL::Identity<Image_word_type> >
 void initialize_triangulation_from_labeled_image(C3T3& c3t3,
       const MeshDomain&   domain,
       const CGAL::Image_3& image,
       const MeshCriteria& criteria,
       Image_word_type,
-                                                 bool protect_features = false
+      bool protect_features = false,
-                                                 )
+      TransformOperator transform = CGAL::Identity<Image_word_type>())
 {
   typedef typename C3T3::Triangulation       Tr;
   typedef typename Tr::Geom_traits           Gt;
@@ -111,11 +118,10 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3,
   Seeds seeds;
   Get_point<Bare_point> get_point(&image);
   std::cout << "Searching for connected components..." << std::endl;
-  CGAL::Identity<Image_word_type> no_transformation;
   search_for_connected_components_in_labeled_image(image,
                                                    std::back_inserter(seeds),
                                                    CGAL::Emptyset_iterator(),
-                                                   no_transformation,
+                                                   transform,
                                                    get_point,
                                                    Image_word_type());
   std::cout << "  " << seeds.size() << " components were found." << std::endl;