Mesh_3: Add weighted images

BGL/include/CGAL/boost/parameter.h

@@ -110,6 +110,7 @@ BOOST_PARAMETER_NAME( (pointer_to_stop_atomic_boolean, tag ) pointer_to_stop_ato
 BOOST_PARAMETER_NAME( (function, tag ) function_)
 BOOST_PARAMETER_NAME( (bounding_object, tag ) bounding_object_)
 BOOST_PARAMETER_NAME( (relative_error_bound, tag ) relative_error_bound_)
+BOOST_PARAMETER_NAME( (weights, tag) weights_)
 BOOST_PARAMETER_NAME( (p_rng, tag ) p_rng_)
 BOOST_PARAMETER_NAME( (null_subdomain_index, tag ) null_subdomain_index_)
 BOOST_PARAMETER_NAME( (construct_surface_patch_index, tag ) construct_surface_patch_index_)

CGAL_ImageIO/include/CGAL/Image_3.h

@@ -102,7 +102,7 @@ class CGAL_IMAGEIO_EXPORT Image_3
 //     std::cerr << "Image_3::copy_constructor\n";
   }
 
-  Image_3(_image* im, Own own_the_data = OWN_THE_DATA)
+  explicit Image_3(_image* im, Own own_the_data = OWN_THE_DATA)
   {
     private_read(im, own_the_data);
   }
@@ -146,9 +146,9 @@ class CGAL_IMAGEIO_EXPORT Image_3
   double vy() const { return image_ptr->vy; }
   double vz() const { return image_ptr->vz; }
 
-  double tx() const { return image_ptr->tx; }
-  double ty() const { return image_ptr->ty; }
-  double tz() const { return image_ptr->tz; }
+  float tx() const { return image_ptr->tx; }
+  float ty() const { return image_ptr->ty; }
+  float tz() const { return image_ptr->tz; }
 
   float value(const std::size_t i,
               const std::size_t j,
@@ -157,6 +157,11 @@ class CGAL_IMAGEIO_EXPORT Image_3
     return ::evaluate(image(),i,j,k);
   }
 
+  bool is_valid() const
+  {
+    return image_ptr.get() != nullptr;
+  }
+
 public:
 
   bool read(const char* file)

Documentation/doc/biblio/cgal_manual.bib

@@ -3256,6 +3256,12 @@ @article{cgal:Wwshap-eepec-20
      publisher = {ACM}
 }
 
+@article{stalling1998weighted,
+  title={Weighted labels for 3D image segmentation},
+  author={Stalling, Detlev and Z{\"o}ckler, Malte and Sander, Oliver and Hege, Hans-Christian},
+  year={1998}
+}
+
 % ----------------------------------------------------------------------------
 % END OF BIBFILE
 % ----------------------------------------------------------------------------

Installation/cmake/modules/CGAL_ITK_support.cmake

@@ -0,0 +1,7 @@
+if(ITK_FOUND AND NOT TARGET CGAL::ITK_support)
+  add_library(CGAL::ITK_support INTERFACE IMPORTED)
+  set_target_properties(CGAL::ITK_support PROPERTIES
+    INTERFACE_COMPILE_DEFINITIONS "CGAL_USE_ITK"
+    INTERFACE_INCLUDE_DIRECTORIES "${ITK_INCLUDE_DIRS}"
+    INTERFACE_LINK_LIBRARIES "${ITK_LIBRARIES}")
+endif()

Installation/cmake/modules/list_of_whitelisted_headers.cmake

@@ -41,7 +41,7 @@ set(list_of_whitelisted_headers_txt [=[
   CGAL/IO/Triangulation_geomview_ostream_3.h
   CGAL/IO/Triangulation_geomview_ostream_2.h
   CGAL/IO/Polyhedron_geomview_ostream.h
-
+  CGAL/Mesh_3/generate_label_weights.h
 ]=])
 
 separate_arguments(list_of_whitelisted_headers UNIX_COMMAND ${list_of_whitelisted_headers_txt})

Mesh_3/doc/Mesh_3/CGAL/Labeled_mesh_domain_3.h

@@ -233,6 +233,12 @@ The parameters are optional unless otherwise specified.
 <li> <b>`parameters::image` (mandatory)</b> the input 3D image. Must
 be a `CGAL::Image_3` object.
 
+<li> <b>`parameters::weights`</b> an input 3D image that provides
+weights associated to each voxel (the word type is `unsigned char`,
+and the voxels values are integers between 0 and 255).
+The weights image can be generated with `CGAL::Mesh_3::generate_label_weights()`.
+Its dimensions must be the same as the dimensions of `parameters::image`.
+
 <li><b>`parameter::value_outside`</b> the value attached to voxels
  outside of the domain to be meshed. Its default value is `0`.
 
@@ -246,6 +252,11 @@ From the example (\ref Mesh_3/mesh_3D_image.cpp):
 
 \snippet Mesh_3/mesh_3D_image.cpp Domain creation
 
+From the example (\ref Mesh_3/mesh_3D_weighted_image.cpp),
+where the labeled image is used with a precomputed 3D image of weights :
+
+\snippet Mesh_3/mesh_3D_weighted_image.cpp Domain creation
+
  */
 template <typename ... A_i>
 static

Mesh_3/doc/Mesh_3/Doxyfile.in

@@ -5,7 +5,8 @@ ALIASES += "cgalDescribePolylineType=A polyline is defined as a sequence of poin
 
 INPUT += \
       ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Polyhedral_complex_mesh_domain_3.h \
-      ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_domain_with_polyline_features_3.h
+      ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_domain_with_polyline_features_3.h \
+      ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/generate_label_weights.h
 
 PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Mesh Generation"
 HTML_EXTRA_FILES           =  ${CGAL_PACKAGE_DOC_DIR}/fig/implicit_domain_3.jpg \

Mesh_3/doc/Mesh_3/Mesh_3.txt

@@ -706,6 +706,24 @@ 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_3DomainsFrom3DImagesWithWeights Domains From Segmented 3D Images, with Weights
+
+When a segmented image is given as input, the generated mesh surface sometimes sticks too closely
+to the voxels surface, causing an aliasing effect.
+A solution to generate a smooth and accurate output surface was described by Stalling et al in
+\cgalCite{stalling1998weighted}. It consists in generating a second input image, made
+of integer coefficients called *weights*, and use those weights to define smoother domain boundaries.
+The 3D image of weights can be generated using `CGAL::Mesh_3::generate_weights()` as shown in
+the following example.
+
+\cgalExample{Mesh_3/mesh_3D_weighted_image.cpp}
+
+\cgalFigureBegin{figure_weightedImage, weighted_images.jpg}
+Surface of the output mesh generated with a very small `facet_distance`
+without the weights (left, 25563 vertices) and with the weights (right, 19936 vertices).
+\cgalFigureEnd
+
+
 \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

Mesh_3/doc/Mesh_3/PackageDescription.txt

@@ -115,6 +115,7 @@ and their associated classes:
 - `CGAL::lloyd_optimize_mesh_3()`
 - `CGAL::odt_optimize_mesh_3()`
 - `CGAL::facets_in_complex_3_to_triangle_mesh()`
+- `CGAL::Mesh_3::generate_label_weights()`
 
 \cgalCRPSection{CGAL::parameters Functions}
 

Mesh_3/doc/Mesh_3/examples.txt

@@ -4,6 +4,7 @@
 \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/mesh_3D_weighted_image.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

Mesh_3/examples/Mesh_3/CMakeLists.txt

@@ -161,10 +161,21 @@ if(TARGET CGAL::CGAL_ImageIO)
     create_single_source_cgal_program("mesh_3D_image_variable_size.cpp")
     target_link_libraries(mesh_3D_image_variable_size
                           PUBLIC CGAL::Eigen3_support)
+
+    find_package(ITK QUIET COMPONENTS ITKCommon ITKThresholding ITKSmoothing ITKImageIntensity)
+    if(ITK_FOUND)
+      include(CGAL_ITK_support)
+      create_single_source_cgal_program("mesh_3D_weighted_image.cpp")
+      target_link_libraries(mesh_3D_weighted_image
+                            PUBLIC CGAL::Eigen3_support CGAL::ITK_support)
+    else(ITK_FOUND)
+        message(STATUS "NOTICE: The examples that need ITK will not be compiled.")
+    endif(ITK_FOUND)
+
   else()
     message(
       STATUS
-        "NOTICE: The examples mesh_3D_image.cpp, mesh_3D_image_variable_size.cpp, mesh_optimization_example.cpp and mesh_optimization_lloyd_example.cpp need CGAL_ImageIO to be configured with ZLIB support, and will not be compiled."
+        "NOTICE: The examples mesh_3D_image.cpp, mesh_3D_weighted_image.cpp, mesh_3D_image_variable_size.cpp, mesh_optimization_example.cpp and mesh_optimization_lloyd_example.cpp need CGAL_ImageIO to be configured with ZLIB support, and will not be compiled."
     )
   endif()
 
@@ -179,6 +190,7 @@ if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support)
   foreach(
     target
     mesh_3D_image
+    mesh_3D_weighted_image
     mesh_3D_image_variable_size
     mesh_3D_image_with_custom_initialization
     mesh_3D_gray_image_with_custom_initialization

Mesh_3/examples/Mesh_3/mesh_3D_weighted_image.cpp

@@ -0,0 +1,67 @@
+#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/Labeled_mesh_domain_3.h>
+#include <CGAL/Mesh_3/generate_label_weights.h>
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/Image_3.h>
+#include <CGAL/IO/File_binary_mesh_3.h>
+#include <CGAL/tags.h>
+
+// Domain
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
+
+// Triangulation
+using Tr   = CGAL::Mesh_triangulation_3<Mesh_domain,
+                                        CGAL::Default,
+                                        CGAL::Parallel_if_available_tag>::type;
+using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
+
+// Criteria
+using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+int main(int argc, char* argv[])
+{
+  /// [Loads image]
+  const char* fname = (argc > 1) ? argv[1] : "data/liver.inr.gz";
+  CGAL::Image_3 image;
+  if(!image.read(fname)){
+    std::cerr << "Error: Cannot read file " <<  fname << std::endl;
+    return EXIT_FAILURE;
+  }
+  /// [Loads image]
+
+  /// [Domain creation]
+  const float sigma = 10.f;
+  CGAL::Image_3 img_weights =
+    CGAL::Mesh_3::generate_label_weights(image, sigma);
+
+  Mesh_domain domain
+    = Mesh_domain::create_labeled_image_mesh_domain(image,
+                                                    weights = img_weights,
+                                                    relative_error_bound = 1e-6);
+  /// [Domain creation]
+
+  // Mesh criteria
+  Mesh_criteria criteria(facet_angle=30, facet_size=6, facet_distance=0.5,
+                         cell_radius_edge_ratio=3, cell_size=8);
+
+  /// [Meshing]
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+  /// [Meshing]
+
+  // Output
+  std::ofstream medit_file("out.mesh");
+  c3t3.output_to_medit(medit_file);
+  std::ofstream bin_file("out.binary.cgal", std::ios_base::binary);
+  CGAL::IO::save_binary_file(bin_file, c3t3);
+
+  return 0;
+}

Mesh_3/include/CGAL/Labeled_mesh_domain_3.h

@@ -39,6 +39,7 @@
 // support for `CGAL::Image_3`
 #include <CGAL/Image_3.h>
 #include <CGAL/Mesh_3/Image_to_labeled_function_wrapper.h>
+#include <CGAL/Mesh_3/Image_plus_weights_to_labeled_function_wrapper.h>
 
 // support for implicit functions
 #include <CGAL/Implicit_to_labeling_function_wrapper.h>
@@ -418,6 +419,7 @@ class Labeled_mesh_domain_3 :
                                   (optional
                                    (relative_error_bound_, (const FT&),
                                     FT(1e-3))
+                                   (weights_, (const CGAL::Image_3&), CGAL::Image_3())
                                    (value_outside_, *, 0)
                                    (p_rng_, (CGAL::Random*), (CGAL::Random*)(0))
                                    (image_values_to_subdomain_indices_, *,
@@ -429,18 +431,37 @@ class Labeled_mesh_domain_3 :
                                   )
   {
     namespace p = CGAL::parameters;
-    return Labeled_mesh_domain_3
-      (create_labeled_image_wrapper
+    if (weights_.is_valid())
+    {
+      return Labeled_mesh_domain_3
+      (create_weighted_labeled_image_wrapper
        (image_,
+        weights_,
         image_values_to_subdomain_indices_,
         value_outside_),
-       Mesh_3::internal::compute_bounding_box(image_),
-       p::relative_error_bound = relative_error_bound_,
-       p::p_rng = p_rng_,
-       p::null_subdomain_index =
-         create_null_subdomain_index(null_subdomain_index_),
-       p::construct_surface_patch_index =
-         create_construct_surface_patch_index(construct_surface_patch_index_));
+        Mesh_3::internal::compute_bounding_box(image_),
+        p::relative_error_bound = relative_error_bound_,
+        p::p_rng = p_rng_,
+        p::null_subdomain_index =
+        create_null_subdomain_index(null_subdomain_index_),
+        p::construct_surface_patch_index =
+        create_construct_surface_patch_index(construct_surface_patch_index_));
+    }
+    else
+    {
+      return Labeled_mesh_domain_3
+      (create_labeled_image_wrapper
+      (image_,
+        image_values_to_subdomain_indices_,
+        value_outside_),
+        Mesh_3::internal::compute_bounding_box(image_),
+        p::relative_error_bound = relative_error_bound_,
+        p::p_rng = p_rng_,
+        p::null_subdomain_index =
+        create_null_subdomain_index(null_subdomain_index_),
+        p::construct_surface_patch_index =
+        create_construct_surface_patch_index(construct_surface_patch_index_));
+    }
   }
 
   BOOST_PARAMETER_MEMBER_FUNCTION(
@@ -866,6 +887,33 @@ class Labeled_mesh_domain_3 :
                    transform_fct(value_outside));
   }
 
+  template <typename Image_word_type,
+            typename FT, typename Functor>
+  static
+  Function
+  create_weighted_labeled_image_wrapper_with_know_word_type
+  (const CGAL::Image_3& image,
+   const CGAL::Image_3& weights,
+   const Functor& image_values_to_subdomain_indices,
+   const FT& value_outside)
+  {
+    using Mesh_3::internal::Create_labeled_image_values_to_subdomain_indices;
+    typedef Create_labeled_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);
+
+    typedef Mesh_3::Image_plus_weights_to_labeled_function_wrapper<
+      Image_word_type,
+      int, //interpolation type
+      unsigned char, // Weights_type,
+      Subdomain_index> Wrapper;
+    return Wrapper(image,
+                   weights,
+                   transform_fct,
+                   transform_fct(value_outside));
+  }
+
   template <typename FT, typename Functor>
   static
   Function
@@ -885,6 +933,27 @@ class Labeled_mesh_domain_3 :
     return Function();
   }
 
+  template <typename FT, typename Functor>
+  static
+  Function
+  create_weighted_labeled_image_wrapper(const CGAL::Image_3& image,
+                                        const CGAL::Image_3& weights,
+                                        const Functor& image_values_to_subdomain_indices,
+                                        const FT& value_outside)
+  {
+    CGAL_IMAGE_IO_CASE(image.image(),
+      return create_weighted_labeled_image_wrapper_with_know_word_type<Word>
+                        (image,
+                         weights,
+                         image_values_to_subdomain_indices,
+                         value_outside);
+                        );
+    CGAL_error_msg("This place should never be reached, because it would mean "
+      "the image word type is a type that is not handled by "
+      "CGAL_ImageIO.");
+    return Function();
+  }
+
   static
   Construct_surface_patch_index
   create_construct_surface_patch_index(const Null_functor&) {