Merge pull request #6129 from janetournois/Mesh_3-triple_line_extraction-GF

Mesh_3 - add detection of intersection lines from labeled images

# Conflicts:
#	Installation/CHANGES.md

CGAL_ImageIO/include/CGAL/ImageIO.h

@@ -20,7 +20,7 @@
 
 #include <stdlib.h>
 #include <stdio.h>
-#include <boost/cstdint.hpp> // for uint32_t, etc.
+#include <cstdint> // for uint32_t, etc.
 
 #ifdef CGAL_USE_ZLIB
 #include <zlib.h>
@@ -562,38 +562,50 @@ struct Word_type_generator<WK_FLOAT, sign, 8>
 template <>
 struct Word_type_generator<WK_FIXED, SGN_SIGNED, 1>
 {
-//   typedef boost::int8_t type;
+//   typedef std::int8_t type;
   typedef char type;
 };
 
 template <>
 struct Word_type_generator<WK_FIXED, SGN_UNSIGNED, 1>
 {
-  typedef boost::uint8_t type;
+  typedef std::uint8_t type;
 };
 
 template <>
 struct Word_type_generator<WK_FIXED, SGN_SIGNED, 2>
 {
-  typedef boost::int16_t type;
+  typedef std::int16_t type;
 };
 
 template <>
 struct Word_type_generator<WK_FIXED, SGN_UNSIGNED, 2>
 {
-  typedef boost::uint16_t type;
+  typedef std::uint16_t type;
 };
 
 template <>
 struct Word_type_generator<WK_FIXED, SGN_SIGNED, 4>
 {
-  typedef boost::int32_t type;
+  typedef std::int32_t type;
 };
 
 template <>
 struct Word_type_generator<WK_FIXED, SGN_UNSIGNED, 4>
 {
-  typedef boost::uint32_t type;
+  typedef std::uint32_t type;
+};
+
+template <>
+struct Word_type_generator<WK_FIXED, SGN_SIGNED, 8>
+{
+  typedef std::int64_t type;
+};
+
+template <>
+struct Word_type_generator<WK_FIXED, SGN_UNSIGNED, 8>
+{
+  typedef std::uint64_t type;
 };
 
 template <WORD_KIND wordKind, SIGN sign, std::size_t wdim>

Documentation/doc/biblio/cgal_manual.bib

@@ -3270,6 +3270,14 @@ pages = "207--221"
   year={1998}
 }
 
+
+@techreport{cgal:hssz-gmcabonbc-97,
+  title={A generalized marching cubes algorithm based on non-binary classifications},
+  author={H-C. Hege and M. Seebass and D. Stalling and M. Zöckler},
+  number={SC 97-05},
+  year={1997}
+}
+
 % ----------------------------------------------------------------------------
 % END OF BIBFILE
 % ----------------------------------------------------------------------------

Installation/CHANGES.md

@@ -38,6 +38,13 @@ CGAL tetrahedral Delaunay refinement algorithm.
 
 -   This new package wraps all the existing code that deals with a `MeshComplex_3InTriangulation_3` to describe 3D simplicial meshes, and makes the data structure independent from the tetrahedral mesh generation package.
 
+### [Tetrahedral Mesh Generation](https://doc.cgal.org/5.6/Manual/packages.html#PkgMesh3)
+
+- Added two new named parameters to the named constructor `CGAL::create_labeled_image_mesh_domain()`
+for automatic detection and protection
+of 1D-curves that lie at the intersection of three or more subdomains,
+extracted from labeled images.
+
 ### [Shape Detection](https://doc.cgal.org/5.6/Manual/packages.html#PkgShapeDetection) (breaking change, major changes)
 
 -   **Breaking change**: The region growing part of the package have been reworked to fix design issues introduced with the handling `FaceGraph` models.
@@ -59,6 +66,7 @@ CGAL tetrahedral Delaunay refinement algorithm.
     and cylinders in 3D.
 
 ### [2D Arrangements](https://doc.cgal.org/5.6/Manual/packages.html#PkgArrangementOnSurface2)
+
 -   Fixed some code that handles geodesic-curves on spheres that compare x- and y-coordinates on the boundary of the parameter space. It mainly effected the naive point-location.
 
 ### [2D Convex Hulls](https://doc.cgal.org/5.6/Manual/packages.html#PkgConvexHull2)

Mesh_3/doc/Mesh_3/Doxyfile.in

@@ -18,7 +18,9 @@ INPUT += \
       ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_facet_topology.h \
       ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_vertex_base_3.h \
       ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_cell_base_3.h \
-      ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Compact_mesh_cell_base_3.h
+      ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Compact_mesh_cell_base_3.h \
+      ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/Detect_features_in_image.h \
+      ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/Detect_features_on_image_bbox.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

@@ -735,7 +735,6 @@ 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
@@ -926,23 +925,27 @@ The first modification is the type of the mesh domain. Instead of being
 
 \snippet Mesh_3/mesh_3D_image_with_features.cpp Domain definition
 
-Then, in the function `%main()`, after the `%domain` object has been created,
+In the %main() function, the domain is created with an additional argument - a
-a dedicated function computes the 1D-features, and adds them to the domain.
+dedicated functor that computes the one-dimensional features, which are then
+added to the domain.
 
-\snippet Mesh_3/mesh_3D_image_with_features.cpp Call add_1D_features
+\snippet Mesh_3/mesh_3D_image_with_detection_of_features.cpp Domain creation
 
-The function template `%add_1D_features()` is defined in the example
+The `CGAL::Mesh_3::Detect_features_in_image` functor is defined in its own
-file. It uses non-documented code from \cgal, that should be copy-pasted in
+header file. It computes the one-dimensional features that correspond to the
-any user-code willing to use similar code. It uses the undocumented
+intersections of the bounding box of the image with the surfaces defined by the
-function template `%CGAL::polylines_to_protect` that computes the
+image, as well as polylines that lie at the intersection of three or more
-1D-features that correspond to the intersection of the bounding box of the
+subdomains (including the outside). It then constructs a graph of these polyline
-image with the surfaces defined by the image. At the same time, a few other
+features. The named constructor adds this feature graph to the domain for later
-polylines are added as 1D-features, to protect 1D curves in the interior of
+feature protection. The original feature detection algorithm was described in
-the image. Then, the method
+\cgalCite{cgal:hssz-gmcabonbc-97}, which provides a list of possible voxel
-`CGAL::Mesh_domain_with_polyline_features_3::add_features` is called twice
+configurations. The feature detection implemented in \cgal generalizes this
-to add the computed 1D-features to the mesh domain.
+description.
 
-\snippet Mesh_3/mesh_3D_image_with_features.cpp Add 1D features
+The example \ref Mesh_3/mesh_3D_image_with_features.cpp shows how
+user-specified input polylines can further be added as 1D features to the mesh domain.
+
+\snippet Mesh_3/mesh_3D_image_with_features.cpp Domain creation
 
 In the meshing criteria, if 1D features are added to the domain, the user
 can define the parameter `edge_size` of the criteria class
@@ -970,7 +973,7 @@ protection of the 1D-features.
 Left: the mesh without any call to `%add_features()`.
 
 Middle: the mesh with only the 1D-features computed by
-`%CGAL::polylines_to_protect()`.
+`CGAL::Mesh_3::Detect_features_on_image_bbox()`.
 
 Right: the mesh with added 1D-features in the interior of the bounding box
 of the image.

Mesh_3/doc/Mesh_3/PackageDescription.txt

@@ -16,6 +16,10 @@
 /// \ingroup PkgMesh3Ref
 /// The classes in this group are models of domain concepts and their associated classes.
 
+/// \defgroup PkgMesh3FeatureDetection Feature Detection
+/// \ingroup PkgMesh3Ref
+/// The functors in this group perform polyline features detection in input domains.
+
 /// \defgroup PkgMesh3Functions Mesh Generation Functions
 /// \ingroup PkgMesh3Ref
 /// The two main functions to generate a mesh are `make_mesh_3()` and `refine_mesh_3()`.
@@ -103,6 +107,11 @@ and their associated classes:
 - `CGAL::Labeled_image_mesh_domain_3<Image,BGT>` (deprecated)
 - `CGAL::Gray_image_mesh_domain_3<Image,BGT,Image_word_type>` (deprecated)
 
+The following functors are available for feature detection:
+
+- `CGAL::Mesh_3::Detect_features_in_image`
+- `CGAL::Mesh_3::Detect_features_on_image_bbox`
+
 \cgalCRPSection{Function Templates}
 
 - `CGAL::make_mesh_3()`

Mesh_3/doc/Mesh_3/examples.txt

@@ -4,6 +4,8 @@
 \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_image_with_detection_of_features.cpp
+\example Mesh_3/mesh_3D_image_with_input_features.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

Mesh_3/examples/Mesh_3/CMakeLists.txt

@@ -133,6 +133,12 @@ if(TARGET CGAL::CGAL_ImageIO)
   create_single_source_cgal_program("mesh_3D_image_with_features.cpp")
   target_link_libraries(mesh_3D_image_with_features PUBLIC CGAL::Eigen3_support)
 
+  create_single_source_cgal_program("mesh_3D_image_with_input_features.cpp")
+  target_link_libraries(mesh_3D_image_with_input_features PUBLIC CGAL::Eigen3_support)
+
+  create_single_source_cgal_program("mesh_3D_image_with_detection_of_features.cpp")
+  target_link_libraries(mesh_3D_image_with_detection_of_features PUBLIC CGAL::Eigen3_support)
+
   if(CGAL_ImageIO_USE_ZLIB)
     create_single_source_cgal_program("mesh_optimization_example.cpp")
     target_link_libraries(mesh_optimization_example PUBLIC CGAL::Eigen3_support)
@@ -185,6 +191,8 @@ if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support)
     mesh_3D_image_with_custom_initialization
     mesh_3D_gray_image_with_custom_initialization
     mesh_3D_image_with_features
+    mesh_3D_image_with_detection_of_features
+    mesh_3D_image_with_input_features
     mesh_implicit_domains
     mesh_implicit_sphere
     mesh_implicit_sphere_variable_size

Mesh_3/examples/Mesh_3/mesh_3D_image_variable_size.cpp

@@ -1,5 +1,3 @@
-#define CGAL_MESH_3_VERBOSE
-
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
 #include <CGAL/Mesh_triangulation_3.h>

Mesh_3/examples/Mesh_3/mesh_3D_image_with_detection_of_features.cpp

@@ -0,0 +1,86 @@
+#include <vector>
+#include <iostream>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/Image_3.h>
+
+/// [Domain definition]
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Mesh_domain_with_polyline_features_3.h>
+#include <CGAL/Labeled_mesh_domain_3.h>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Labeled_mesh_domain_3<K> Image_domain;
+typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Mesh_domain;
+/// [Domain definition]
+
+#include <CGAL/Mesh_3/Detect_features_in_image.h>
+
+#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
+namespace params = CGAL::parameters;
+
+int main(int argc, char* argv[])
+{
+  const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/420.inr");
+
+  /// [Loads image]
+  CGAL::Image_3 image;
+  if(!image.read(fname)){
+    std::cerr << "Error: Cannot read file " <<  fname << std::endl;
+    return EXIT_FAILURE;
+  }
+  /// [Loads image]
+
+  /// [Domain creation]
+  Mesh_domain domain
+    = Mesh_domain::create_labeled_image_mesh_domain(image,
+         params::features_detector = CGAL::Mesh_3::Detect_features_in_image());
+  /// [Domain creation]
+
+  CGAL::Bbox_3 bbox = domain.bbox();
+  double diag = CGAL::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                           CGAL::square(bbox.ymax() - bbox.ymin()) +
+                           CGAL::square(bbox.zmax() - bbox.zmin()));
+  double sizing_default = diag * 0.05;
+
+  /// [Mesh criteria]
+  /// Note that `edge_size` is needed with 1D-features
+  Mesh_criteria criteria(params::edge_size = sizing_default,
+    params::facet_angle = 30,
+    params::facet_size = sizing_default,
+    params::facet_distance = sizing_default / 10,
+    params::facet_topology = CGAL::FACET_VERTICES_ON_SAME_SURFACE_PATCH,
+    params::cell_radius_edge_ratio = 0,
+    params::cell_size = 0
+  );
+  /// [Mesh criteria]
+
+  /// [Meshing]
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
+                                      params::no_exude(),
+                                      params::no_perturb());
+  /// [Meshing]
+
+  // Output
+  CGAL::dump_c3t3(c3t3, "out");
+
+  return 0;
+}

Mesh_3/examples/Mesh_3/mesh_3D_image_with_features.cpp

@@ -18,6 +18,8 @@ typedef CGAL::Labeled_mesh_domain_3<K> Image_domain;
 typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Mesh_domain;
 /// [Domain definition]
 
+#include <CGAL/Mesh_3/Detect_features_on_image_bbox.h>
+
 #ifdef CGAL_CONCURRENT_MESH_3
 typedef CGAL::Parallel_tag Concurrency_tag;
 #else
@@ -34,40 +36,8 @@ typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
 
 namespace params = CGAL::parameters;
 
-/// [Add 1D features]
+// Read input features
 #include "read_polylines.h"
-#include <CGAL/Mesh_3/polylines_to_protect.h> // undocumented header
-
-// Protect the intersection of the object with the box of the image,
-// by declaring 1D-features. Note that `CGAL::polylines_to_protect` is
-// not documented.
-bool add_1D_features(const CGAL::Image_3& image,
-                     Mesh_domain& domain,
-                     const std::string lines_fname)
-{
-  typedef K::Point_3 Point_3;
-  typedef unsigned char Word_type;
-
-  std::vector<std::vector<Point_3> > features_inside;
-  if(!read_polylines(lines_fname, features_inside)) // see file "read_polylines.h"
-  {
-    std::cerr << "Error: Cannot read file " <<  lines_fname << std::endl;
-    return false;
-  }
-
-  std::vector<std::vector<Point_3> > polylines_on_bbox;
-  CGAL::polylines_to_protect<Point_3, Word_type>(image, polylines_on_bbox,
-                                                 features_inside.begin(),
-                                                 features_inside.end());
-
-  domain.add_features(polylines_on_bbox.begin(), polylines_on_bbox.end());
-
-  // It is very important that the polylines from the file `lines_fname`
-  // contain only polylines in the inside of the box of the image.
-  domain.add_features(features_inside.begin(), features_inside.end());
-  return true;
-}
-/// [Add 1D features]
 
 int main(int argc, char* argv[])
 {
@@ -79,21 +49,28 @@ int main(int argc, char* argv[])
     return EXIT_FAILURE;
   }
 
-  // Domain
+  /// Load 1D-features
-  Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image);
+  const std::string lines_fname = (argc > 2) ? argv[2] : CGAL::data_file_path("images/420.polylines.txt");
-
+  std::vector<std::vector<K::Point_3> > features_inside;
-  /// Declare 1D-features, see above [Call add_1D_features]
+  if (!read_polylines(lines_fname, features_inside)) // see file "read_polylines.h"
-  const std::string lines_fname = (argc>2)?argv[2]:CGAL::data_file_path("images/420.polylines.txt");
+  {
-
+    std::cerr << "Error: Cannot read file " << lines_fname << std::endl;
-  if(!add_1D_features(image, domain, lines_fname)) {
     return EXIT_FAILURE;
   }
-  /// [Call add_1D_features]
+
+  /// [Domain creation]
+  Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image,
+    params::features_detector = CGAL::Mesh_3::Detect_features_on_image_bbox(),
+    params::input_features  = std::cref(features_inside));//use std::cref to avoid a copy
+  /// [Domain creation]
 
   /// Note that `edge_size` is needed with 1D-features [Mesh criteria]
-  Mesh_criteria criteria(params::edge_size(6).
+  Mesh_criteria criteria(params::edge_size = 6.,
-                                 facet_angle(30).facet_size(6).facet_distance(4).
+    params::facet_angle = 30,
-                                 cell_radius_edge_ratio(3).cell_size(8));
+    params::facet_size = 6,
+    params::facet_distance = 4,
+    params::cell_radius_edge_ratio = 3,
+    params::cell_size = 8);
   /// [Mesh criteria]
 
   // Meshing
@@ -104,5 +81,5 @@ int main(int argc, char* argv[])
   CGAL::IO::write_MEDIT(medit_file, c3t3);
   medit_file.close();
 
-  return 0;
+  return EXIT_SUCCESS;
 }

Mesh_3/examples/Mesh_3/mesh_3D_image_with_input_features.cpp

@@ -0,0 +1,82 @@
+#include <vector>
+#include <iostream>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/Image_3.h>
+
+/// [Domain definition]
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Mesh_domain_with_polyline_features_3.h>
+#include <CGAL/Labeled_mesh_domain_3.h>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Labeled_mesh_domain_3<K> Image_domain;
+typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Mesh_domain;
+/// [Domain definition]
+
+#ifdef CGAL_CONCURRENT_MESH_3
+using Concurrency_tag = CGAL::Parallel_tag;
+#else
+using Concurrency_tag = CGAL::Sequential_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;
+
+namespace params = CGAL::parameters;
+
+#include "read_polylines.h"
+
+int main(int argc, char* argv[])
+{
+  const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("images/420.inr");
+  // Loads image
+  CGAL::Image_3 image;
+  if(!image.read(fname)){
+    std::cerr << "Error: Cannot read file " <<  fname << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  /// Declare 1D-features
+  const std::string lines_fname = (argc>2)?argv[2]:CGAL::data_file_path("images/420.polylines.txt");
+  using Point_3 = K::Point_3;
+  std::vector<std::vector<Point_3> > features_inside;
+  if (!read_polylines(lines_fname, features_inside)) // see file "read_polylines.h"
+  {
+    std::cerr << "Error: Cannot read file " << lines_fname << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  /// [Domain creation]
+  Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image,
+    params::input_features = std::cref(features_inside));//use std::cref to avoid a copy
+  /// [Domain creation]
+
+  /// Note that `edge_size` is needed with 1D-features [Mesh criteria]
+  Mesh_criteria criteria(params::edge_size = 6,
+    params::facet_angle = 30,
+    params::facet_size = 6,
+    params::facet_distance = 4,
+    params::cell_radius_edge_ratio = 3,
+    params::cell_size = 8);
+  /// [Mesh criteria]
+
+  // Meshing
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+
+  // Output
+  std::ofstream medit_file("out.mesh");
+  CGAL::IO::write_MEDIT(medit_file, c3t3);
+  medit_file.close();
+
+  return EXIT_SUCCESS;
+}

Mesh_3/include/CGAL/Implicit_mesh_domain_3.h

@@ -24,6 +24,7 @@
 
 #include <CGAL/Labeled_mesh_domain_3.h>
 #include <CGAL/Implicit_to_labeling_function_wrapper.h>
+#include <CGAL/Mesh_3/Null_subdomain_index.h>
 #include <CGAL/Random.h>
 
 namespace CGAL {

Mesh_3/include/CGAL/Labeled_image_mesh_domain_3.h

@@ -26,6 +26,7 @@
 #include <CGAL/Random.h>
 #include <CGAL/Labeled_mesh_domain_3.h>
 #include <CGAL/Mesh_3/Image_to_labeled_function_wrapper.h>
+#include <CGAL/Mesh_3/Null_subdomain_index.h>
 #include <CGAL/Bbox_3.h>
 #include <CGAL/Default.h>
 

Mesh_3/include/CGAL/Labeled_mesh_domain_3.h

@@ -31,6 +31,7 @@
 #include <CGAL/Origin.h>
 
 #include <functional>
+#include <type_traits>
 
 #include <CGAL/SMDS_3/internal/Handle_IO_for_pair_of_int.h>
 #include <CGAL/SMDS_3/internal/indices_management.h>
@@ -48,6 +49,11 @@
 #endif
 #include <boost/optional.hpp>
 
+#include <CGAL/Mesh_3/Null_subdomain_index.h>
+#include <CGAL/Mesh_domain_with_polyline_features_3.h>
+
+#include <CGAL/Mesh_3/polylines_to_protect.h>
+
 namespace CGAL {
 namespace Mesh_3 {
 namespace internal {
@@ -117,15 +123,67 @@ namespace internal {
     }
   };
 
+  // Detect_features_in_domain
+  template<typename Point, typename DetectFunctor>
+  struct Detect_features_in_domain {
+    std::vector<std::vector<Point>>
+    operator()(const CGAL::Image_3& image, DetectFunctor functor) const {
+#if defined(BOOST_MSVC) && (BOOST_MSVC < 1910) //before msvc2017
+      return functor.operator()<Point>(image);
+#else
+      return functor.template operator()<Point>(image);
+#endif
+    }
+  };
+  // specialization for `Null_functor`: create the default functor
+  template<typename Point>
+  struct Detect_features_in_domain<Point, Null_functor> {
+    std::vector<std::vector<Point>>
+    operator()(const CGAL::Image_3&, Null_functor) const {
+      return std::vector<std::vector<Point>>();
+    }
+  };
+
+  template<typename Point, typename DetectFunctor>
+  std::vector<std::vector<Point>>
+    detect_features(const CGAL::Image_3& image, DetectFunctor functor)
+  {
+    Detect_features_in_domain<Point, DetectFunctor> detector;
+    return detector(image, functor);
+  }
+
+  template<bool WithFeatures>
+  struct Add_features_in_domain {
+    template<typename MeshDomain, typename InputFeatureRange, typename DetectFunctor>
+    void operator()(const CGAL::Image_3&, MeshDomain&, const InputFeatureRange&, DetectFunctor)
+    {}
+  };
+
+  template<>
+  struct Add_features_in_domain<true>
+  {
+    template<typename MeshDomain, typename InputFeatureRange, typename DetectFunctor>
+    void operator()(const CGAL::Image_3& image,
+                    MeshDomain& domain,
+                    const InputFeatureRange& input_features,
+                    DetectFunctor functor)
+    {
+      using P = typename MeshDomain::Point_3;
+      auto detected_feature_range
+        = CGAL::Mesh_3::internal::detect_features<P>(image, functor);
+
+      CGAL::merge_and_snap_polylines(image, detected_feature_range, input_features);
+
+      if (!input_features.empty())
+        domain.add_features(input_features.begin(), input_features.end());
+      domain.add_features(detected_feature_range.begin(), detected_feature_range.end());
+    }
+  };
+
 } // end namespace CGAL::Mesh_3::internal
 } // end namespace CGAL::Mesh_3
 
 #ifndef DOXYGEN_RUNNING
-struct Null_subdomain_index {
-  template <typename T>
-  bool operator()(const T& x) const { return 0 == x; }
-};
-
 template <typename Subdomain_index>
 struct Construct_pair_from_subdomain_indices {
   typedef std::pair<Subdomain_index, Subdomain_index> result_type;
@@ -539,6 +597,11 @@ public:
    * domain to be discretized is the union of voxels that have non-zero
    * values.
    *
+   * \returns either a `Labeled_mesh_domain_3`,
+   *   or a `Mesh_domain_with_polyline_features_3<Labeled_mesh_domain_3>`
+   *   depending on whether one or more of the named parameters
+   *   `features_detector` and `input_features` are provided.
+   *
    * \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
    * \param image_ the input 3D image.
    * \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below:
@@ -565,6 +628,34 @@ public:
    *                            bound, relative to the diameter of the box of the image.}
    *     \cgalParamDefault{FT(1e-3)}
    *   \cgalParamNEnd
+   *
+   *   \cgalParamNBegin{features_detector}
+   *    \cgalParamDescription{ a functor that implements
+   *      `std::vector<std::vector<Point>> operator()(const Image_3& img) const`,
+   *      where `%Point` matches the mesh domain point type.
+   *      It returns a range of detected polyline features, which are added
+   *      to the domain for feature protection.
+   *      See \ref PkgMesh3FeatureDetection for a list of available functors.}
+   *    \cgalParamDefault{CGAL::Null_functor()}
+   *    \cgalParamExtra{The return type of the function depends on whether this parameter
+                        or `input_features` are provided or not.}
+   *    \cgalParamExtra{If `weights` is provided, this parameter is ignored}
+   *   \cgalParamNEnd
+   *
+   *   \cgalParamNBegin{input_features}
+   *    \cgalParamDescription{ a `Range` of polyline features, represented as `Range`s of `Point_3`.
+   *         Polyline features are added to the domain for further feature protection.
+   *         Input polyline features must be different from the detected features
+   *         and can intersect only at vertices, if they do. Otherwise,
+   *         the meshing process may not terminate.}
+   *    \cgalParamDefault{`std::vector<std::vector<Point_3>>()`}
+   *    \cgalParamExtra{The return type of the function depends on whether this parameter
+                        or `input_features` are provided or not.}
+   *    \cgalParamExtra{It is recommended to pass a const-reference for this parameter,
+   *                    possibly using `std::cref(polylines_range)` to avoid useless copies.}
+   *    \cgalParamExtra{If `weights` is provided, this parameter is ignored}
+   *   \cgalParamNEnd
+   *
    * \cgalNamedParamsEnd
    *
    * \cgalHeading{Example}
@@ -578,13 +669,24 @@ public:
    *
    * \snippet Mesh_3/mesh_3D_weighted_image.cpp Domain creation
    *
+   * From the example (\ref Mesh_3/mesh_3D_image_with_detection_of_features.cpp)
+   * where the features are detected in `image`:
+   *
+   * \snippet Mesh_3/mesh_3D_image_with_detection_of_features.cpp Domain creation
+   *
+   * From the example (\ref Mesh_3/mesh_3D_image_with_input_features.cpp)
+   * where the features are provided by the user:
+   *
+   * \snippet Mesh_3/mesh_3D_image_with_input_features.cpp Domain creation
    */
   template<typename CGAL_NP_TEMPLATE_PARAMETERS>
-  static Labeled_mesh_domain_3 create_labeled_image_mesh_domain(const CGAL::Image_3& image_, const CGAL_NP_CLASS& np = parameters::default_values())
+  static auto
+  create_labeled_image_mesh_domain(const CGAL::Image_3& image_, const CGAL_NP_CLASS& np = parameters::default_values())
   {
     using parameters::get_parameter;
     using parameters::get_parameter_reference;
     using parameters::choose_parameter;
+
     auto iso_value_ = choose_parameter(get_parameter(np, internal_np::iso_value_param), 0);
     auto value_outside_ = choose_parameter(get_parameter(np, internal_np::voxel_value), 0);
     FT relative_error_bound_ = choose_parameter(get_parameter(np, internal_np::error_bound), FT(1e-3));
@@ -592,41 +694,62 @@ public:
     CGAL::Random* p_rng_ = choose_parameter(get_parameter(np, internal_np::rng), nullptr);
     auto null_subdomain_index_ = choose_parameter(get_parameter(np, internal_np::null_subdomain_index_param), Null_functor());
     auto construct_surface_patch_index_ = choose_parameter(get_parameter(np, internal_np::surface_patch_index), Null_functor());
-    const CGAL::Image_3& weights_ = choose_parameter(get_parameter_reference(np, internal_np::weights_param), CGAL::Image_3());
+
+    using Image_ref_type = typename internal_np::Lookup_named_param_def<internal_np::weights_param_t,
+                                                                        CGAL_NP_CLASS,
+                                                                        CGAL::Image_3>::reference;
+    CGAL::Image_3 no_weights;
+    const Image_ref_type weights_ = choose_parameter(get_parameter_reference(np, internal_np::weights_param), no_weights);
+    auto features_detector_ = choose_parameter(get_parameter(np, internal_np::features_detector_param), Null_functor());
+
+    using Default_input_features = std::vector<std::vector<typename Labeled_mesh_domain_3::Point_3>>;
+    using Input_features_ref_type = typename internal_np::Lookup_named_param_def<internal_np::input_features_param_t,
+                                                                                 CGAL_NP_CLASS,
+                                                                                 Default_input_features>::reference;
+    Default_input_features empty_vec;
+    Input_features_ref_type input_features_
+          = choose_parameter(get_parameter_reference(np, internal_np::input_features_param), empty_vec);
+
     CGAL_USE(iso_value_);
     namespace p = CGAL::parameters;
 
+    auto image_wrapper = weights_.is_valid()
+      ? create_weighted_labeled_image_wrapper(image_,
+                                              weights_,
+                                              image_values_to_subdomain_indices_,
+                                              value_outside_)
+      : create_labeled_image_wrapper(image_,
+                                     image_values_to_subdomain_indices_,
+                                     value_outside_);
+
+    // warning : keep Return_type consistent with actual return type
+    const bool no_features
+      =  CGAL::parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::features_detector_param_t>::value
+      && CGAL::parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::input_features_param_t>::value;
+    using Return_type = std::conditional_t <
+      no_features,
+      Labeled_mesh_domain_3,
+      Mesh_domain_with_polyline_features_3<Labeled_mesh_domain_3>
+    >;
+
+    Return_type domain
+      (p::function = image_wrapper,
+       p::bounding_object = 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_));
+
     if (weights_.is_valid())
-    {
+      return domain;
-      return Labeled_mesh_domain_3
+
-              (p::function = create_weighted_labeled_image_wrapper
+    // features
-                       (image_,
+    Mesh_3::internal::Add_features_in_domain<!no_features>()
-                        weights_,
+      (image_, domain, input_features_, features_detector_);
-                        image_values_to_subdomain_indices_,
+
-                        value_outside_),
+    return domain;
-               p::bounding_object = 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
-              (p::function = create_labeled_image_wrapper
-                       (image_,
-                        image_values_to_subdomain_indices_,
-                        value_outside_),
-               p::bounding_object = 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_));
-    }
   }
 /// @}
 
@@ -686,7 +809,7 @@ public:
   }
 
   template<typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT>
-  static Labeled_mesh_domain_3 create_labeled_image_mesh_domain(const CGAL_NP_CLASS& np)
+  static auto create_labeled_image_mesh_domain(const CGAL_NP_CLASS& np)
   {
     static_assert(!parameters::is_default_parameter<CGAL_NP_CLASS, internal_np::image_3_param_t>::value, "Value for required parameter not found");
     using parameters::get_parameter_reference;
@@ -698,10 +821,10 @@ public:
   template<typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_1,
            typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_2,
            typename ... NP>
-  static Labeled_mesh_domain_3 create_labeled_image_mesh_domain(const CGAL::Image_3& image_,
+  static auto create_labeled_image_mesh_domain(const CGAL::Image_3& image_,
-                                                                const CGAL_NP_CLASS_1&  np1,
+                                               const CGAL_NP_CLASS_1&  np1,
-                                                                const CGAL_NP_CLASS_2&  np2,
+                                               const CGAL_NP_CLASS_2&  np2,
-                                                                const NP& ... nps)
+                                               const NP& ... nps)
   {
     return create_labeled_image_mesh_domain(image_, internal_np::combine_named_parameters(np1, np2, nps...));
   }
@@ -709,9 +832,9 @@ public:
   template<typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_1,
            typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_2,
            typename ... NP>
-  static Labeled_mesh_domain_3 create_labeled_image_mesh_domain(const CGAL_NP_CLASS_1&  np1,
+  static auto create_labeled_image_mesh_domain(const CGAL_NP_CLASS_1& np1,
-                                                                const CGAL_NP_CLASS_2&  np2,
+                                               const CGAL_NP_CLASS_2& np2,
-                                                                const NP& ... nps)
+                                               const NP& ... nps)
   {
     return create_labeled_image_mesh_domain(internal_np::combine_named_parameters(np1, np2, nps...));
   }
@@ -904,9 +1027,9 @@ public:
       // null(f(p)) means p is outside the domain
       Subdomain_index index = (r_domain_.function_)(p);
       if ( r_domain_.null(index) )
-        return Subdomain();
+        return Subdomain{};
       else
-        return Subdomain(index);
+        return Subdomain{ index };
     }
   private:
     const Labeled_mesh_domain_3& r_domain_;

Mesh_3/include/CGAL/Mesh_3/Detect_features_in_image.h

@@ -0,0 +1,290 @@
+// Copyright (c) 2022 GeometryFactory (France).
+// 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)     : Sebastien Loriot, Jane Tournois
+//
+//******************************************************************************
+//
+//******************************************************************************
+
+#ifndef CGAL_MESH_3_DETECT_FEATURES_IN_IMAGE_H
+#define CGAL_MESH_3_DETECT_FEATURES_IN_IMAGE_H
+
+#include <CGAL/license/Mesh_3.h>
+
+#include <CGAL/ImageIO.h>
+
+#include <CGAL/Delaunay_triangulation_3.h>
+
+#include <CGAL/Mesh_3/features_detection/features_detection.h>
+#include <CGAL/Mesh_3/features_detection/coordinates.h>
+#include <CGAL/Mesh_3/features_detection/combinations.h>
+#include <CGAL/Mesh_3/features_detection/cases_table.h>
+#include <CGAL/Mesh_3/features_detection/cube_isometries.h>
+#include <CGAL/Mesh_3/features_detection/features_detection_helpers.h>
+
+#include <CGAL/Mesh_3/polylines_to_protect.h>
+#include <CGAL/Kernel_traits.h>
+
+#include <vector>
+#include <array>
+
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+#include <boost/range/join.hpp>
+#endif
+
+namespace CGAL
+{
+namespace Mesh_3
+{
+namespace internal
+{
+
+// Protect the intersection of the object with the box of the image,
+// by declaring 1D-features. Note that `CGAL::polylines_to_protect` is
+// not documented.
+template<typename Word_type, typename P>
+std::vector<std::vector<P>>
+detect_features_in_image_with_know_word_type(const CGAL::Image_3& image)
+{
+  using Gt = typename CGAL::Kernel_traits<P>::Kernel;
+  using Point_3 = P;
+  using Vector_3 = typename Gt::Vector_3;
+  using Polyline_type = std::vector<Point_3>;
+  using Polylines = std::vector<Polyline_type>;
+
+  CGAL::Mesh_3::Triple_line_extractor<Point_3> lines;
+
+  Polylines features_inside;
+
+  const double vx = image.vx();
+  const double vy = image.vy();
+  const double vz = image.vz();
+  const double dist_bound = (std::min)(vx, (std::min)(vy, vz)) / 256;
+  const double sq_dist_bound = dist_bound * dist_bound;
+
+  const std::size_t xdim = image.xdim();
+  const std::size_t ydim = image.ydim();
+  const std::size_t zdim = image.zdim();
+
+  const float tx = image.tx();
+  const float ty = image.ty();
+  const float tz = image.tz();
+
+  using CGAL::IMAGEIO::static_evaluate;
+
+  using Del = CGAL::Delaunay_triangulation_3<Gt>;
+  using Cell_handle = typename Del::Cell_handle;
+  using Vertex_handle = typename Del::Vertex_handle;
+  Del triangulation;
+  Cell_handle start_cell;
+
+  using Word //use unsigned integral Word type to use it as an index
+    = typename CGAL::IMAGEIO::Word_type_generator<WK_FIXED, SGN_UNSIGNED, sizeof(Word_type)>::type;
+
+  using Color_transform = internal::Color_transformation_helper<Word>;
+  typename Color_transform::type color_transformation;
+  std::array<Word, 8>            inv_color_transformation;
+
+  using Permutation = internal::Permutation;
+  using Coord = internal::Coordinates;
+
+  for (std::size_t k = 0, end_k = zdim - 1; k < end_k; ++k)
+    for (std::size_t j = 0, end_j = ydim - 1; j < end_j; ++j)
+      for (std::size_t i = 0, end_i = xdim - 1; i < end_i; ++i)
+      {
+        Vector_3 translation{ i * vx + tx,
+                              j * vy + ty,
+                              k * vz + tz };
+
+        const std::array<Word, 8> cube = {
+          static_evaluate<Word>(image.image(), i  , j  , k),
+          static_evaluate<Word>(image.image(), i + 1, j  , k),
+          static_evaluate<Word>(image.image(), i  , j + 1, k),
+          static_evaluate<Word>(image.image(), i + 1, j + 1, k),
+          static_evaluate<Word>(image.image(), i  , j  , k + 1),
+          static_evaluate<Word>(image.image(), i + 1, j  , k + 1),
+          static_evaluate<Word>(image.image(), i  , j + 1, k + 1),
+          static_evaluate<Word>(image.image(), i + 1, j + 1, k + 1),
+        }; /// TODO: optimize the access to the image data
+        bool monocolor = (cube[0] == cube[1]);
+        for (int i = 2; i < 8; ++i) monocolor = monocolor && (cube[0] == cube[i]);
+        if (monocolor) continue;
+
+        Color_transform::reset(color_transformation);
+
+        std::uint8_t nb_color = 0;
+        for (int i = 0; i < 8; ++i) {
+          if (!Color_transform::is_valid(color_transformation, cube[i]))
+          {
+            color_transformation[cube[i]] = nb_color;
+            inv_color_transformation[nb_color] = cube[i];
+            ++nb_color;
+          }
+        }
+        std::array<std::uint8_t, 8> reference_cube = {
+            color_transformation[cube[0]],
+            color_transformation[cube[1]],
+            color_transformation[cube[2]],
+            color_transformation[cube[3]],
+            color_transformation[cube[4]],
+            color_transformation[cube[5]],
+            color_transformation[cube[6]],
+            color_transformation[cube[7]]
+        };
+        auto case_it = internal::find_case(internal::cases, reference_cube);
+        const bool case_found = (case_it != std::end(internal::cases));
+        if (case_found) reference_cube = internal::combinations[(*case_it)[8]];
+        else {
+          //std::cerr << "Warning: case not found: " << reference_cube << '\n';
+          CGAL_error();
+        };
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+        CGAL::Mesh_3::internal::debug_cerr("Cube", cube);
+        CGAL::Mesh_3::internal::debug_cerr("reference cube", reference_cube);
+        CGAL::Mesh_3::internal::debug_cerr("with transformation", internal::cube_isometries[(*case_it)[9]]);
+#endif // CGAL_DEBUG_TRIPLE_LINES
+
+        auto fct_it = lines.create_polylines_fcts.find(reference_cube);
+        if (fct_it != lines.create_polylines_fcts.end())
+        {
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+          CGAL::Mesh_3::internal::debug_cerr("Using the function of", Cube(fct_it->first));
+#endif // CGAL_DEBUG_TRIPLE_LINES
+
+          Polylines cube_features = (fct_it->second)(10);
+          if (case_found)
+          {
+            const Permutation& transformation = internal::cube_isometries[(*case_it)[9]];
+
+            Coord a1 = internal::coordinates[transformation[0]];
+            Coord u = internal::minus(internal::coordinates[transformation[1]], a1);
+            Coord v = internal::minus(internal::coordinates[transformation[2]], a1);
+            Coord w = internal::minus(internal::coordinates[transformation[4]], a1);
+
+            const Point_3  pa{ a1[0], a1[1], a1[2] };
+            const Vector_3 vu{ u[0], u[1], u[2] };
+            const Vector_3 vv{ v[0], v[1], v[2] };
+            const Vector_3 vw{ w[0], w[1], w[2] };
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+            std::cerr << "pa: " << pa << "\n";
+            std::cerr << "vu: " << vu << "\n";
+            std::cerr << "vv: " << vv << "\n";
+            std::cerr << "vw: " << vw << "\n";
+#endif // CGAL_DEBUG_TRIPLE_LINES
+            for (auto& polyline : cube_features) {
+              for (auto& point : polyline) {
+                point = pa
+                  + point.x() * vu
+                  + point.y() * vv
+                  + point.z() * vw;
+                point = { vx * point.x(),
+                          vy * point.y(),
+                          vz * point.z(), };
+                point = point + translation;
+              }
+              for (int i = 0; i < 2; ++i) {
+                Point_3& extremity = (i == 0) ? polyline.front() : polyline.back();
+                Vertex_handle vh = triangulation.nearest_vertex(extremity, start_cell);
+                if (Vertex_handle() != vh) {
+                  if (squared_distance(vh->point(), extremity) < sq_dist_bound) {
+                    extremity = vh->point();
+                  }
+                }
+                vh = triangulation.insert(extremity, start_cell);
+                start_cell = vh->cell();
+              }
+              features_inside.push_back(std::move(polyline));
+            } // end loop on polylines
+          } // end case where the transformation is not the identity
+        } // end if the reference_cube has polylines
+      }
+
+  // call the split_graph_into_polylines, to create long polylines from the
+  // short polylines that were generated per voxel.
+  Polylines polylines_inside;
+  CGAL::polylines_to_protect(polylines_inside,
+                             features_inside.begin(),
+                             features_inside.end());
+
+  Polylines polylines_on_bbox;
+  CGAL::polylines_to_protect<Point_3, Word_type>(image, polylines_on_bbox,
+                                                 polylines_inside.begin(),
+                                                 polylines_inside.end());
+
+  polylines_inside.insert(polylines_inside.end(),
+                          polylines_on_bbox.begin(),
+                          polylines_on_bbox.end());
+
+#ifdef CGAL_DEBUG_TRIPLE_LINES
+  std::ofstream output_polylines("out-generated.polylines.txt");
+  output_polylines.precision(17);
+  for (auto poly : boost::range::join(polylines_on_bbox, polylines_inside)) {
+    output_polylines << poly.size();
+    for (auto p : poly) output_polylines << " " << p;
+    output_polylines << std::endl;
+  }
+#endif
+
+  return polylines_inside;
+}
+
+
+}// namespace internal
+
+/*!
+* \ingroup PkgMesh3FeatureDetection
+*
+* Functor for feature detection in labeled images.
+*/
+struct Detect_features_in_image
+{
+public:
+  /*!
+  * detects and constructs the polylines that lie at the
+  * intersection of three or more subdomains.
+  *
+  * Each subdomain inside the bounding box
+  * of the input labeled image is defined as the set of voxels
+  * with the same value. The outside of the bounding box
+  * of the image is considered as a subdomain with voxel value
+  * `value_outside` (see \link CGAL::Labeled_mesh_domain_3::create_labeled_image_mesh_domain `create_labeled_image_mesh_domain()` \endlink
+  * parameters description). Hence, this function also computes
+  * intersections with the image bounding box.
+  *
+  * \tparam Point class model of `Kernel::Point_3`. It
+  * must match the triangulation point type.
+  *
+  * \param image the input image
+  *
+  * \returns a `std::vector<std::vector<Point>>`
+  * containing the constructed polylines for later feature protection.
+  */
+  template<typename Point>
+  std::vector<std::vector<Point>>
+  operator()(const CGAL::Image_3& image) const
+  {
+    CGAL_IMAGE_IO_CASE(image.image(),
+      return (internal::detect_features_in_image_with_know_word_type<Word, Point>(image));
+    );
+    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 std::vector<std::vector<Point>>();
+  }
+};
+
+
+}//end namespace Mesh_3
+}//end namespace CGAL
+
+
+#endif //CGAL_MESH_3_DETECT_FEATURES_IN_IMAGE_H

Mesh_3/include/CGAL/Mesh_3/Detect_features_on_image_bbox.h

@@ -0,0 +1,100 @@
+// Copyright (c) 2022 GeometryFactory (France).
+// 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_DETECT_FEATURES_ON_IMAGE_BBOX_H
+#define CGAL_MESH_3_DETECT_FEATURES_ON_IMAGE_BBOX_H
+
+#include <CGAL/license/Mesh_3.h>
+
+#include <CGAL/Mesh_3/polylines_to_protect.h>
+
+#include <vector>
+
+
+namespace CGAL
+{
+namespace Mesh_3
+{
+namespace internal
+{
+
+template<typename Point>
+std::vector<std::vector<Point>>
+detect_features_on_bbox(const CGAL::Image_3& image)
+{
+  using Point_3 = Point;
+  using Polyline_type = std::vector<Point_3>;
+  using Polylines = std::vector<Polyline_type>;
+
+  Polylines polylines_on_bbox;
+
+  CGAL_IMAGE_IO_CASE(image.image(),
+    {
+      (CGAL::polylines_to_protect<Point_3, Word>(image, polylines_on_bbox));
+      return polylines_on_bbox;
+    }
+  );
+  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 Polylines();
+}
+
+}// namespace internal
+
+/*!
+* \ingroup PkgMesh3FeatureDetection
+*
+* Functor for feature detection in labeled images.
+*/
+struct Detect_features_on_image_bbox
+{
+public:
+
+  /*!
+  * detects and constructs the polylines that lie at the
+  * intersection of two or more subdomains and the bounding box
+  * of the input labeled image.
+  *
+  * Each subdomain inside the bounding box
+  * of the input labeled image is defined as the set of voxels
+  * with the same value. The outside of the bounding box
+  * of the image is considered as a subdomain with voxel value
+  * `value_outside` (see \link CGAL::Labeled_mesh_domain_3::create_labeled_image_mesh_domain `create_labeled_image_mesh_domain()` \endlink
+  * parameters description). Hence, this function computes
+  * intersections of "internal" subdomains with the image bounding box.
+  *
+  * \tparam Point class model of `Kernel::Point_3`. The point type
+  * must match the triangulation point type.
+  *
+  * \param image the input image
+  *
+  * \returns a `std::vector<std::vector<Point>>`
+  * containing the constructed polylines for later feature protection.
+  */
+  template<typename Point>
+  std::vector<std::vector<Point>>
+  operator()(const CGAL::Image_3& image) const
+  {
+    return internal::detect_features_on_bbox<Point>(image);
+  }
+};
+
+}//end namespace Mesh_3
+}//end namespace CGAL
+
+
+#endif //CGAL_MESH_3_DETECT_FEATURES_ON_IMAGE_BBOX_H

Mesh_3/include/CGAL/Mesh_3/Null_subdomain_index.h

@@ -0,0 +1,29 @@
+// Copyright (c) 2015 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
+
+#ifndef CGAL_MESH_3_NULL_SUBDOMAIN_INDEX
+#define CGAL_MESH_3_NULL_SUBDOMAIN_INDEX
+
+#include <CGAL/license/Mesh_3.h>
+
+#ifndef DOXYGEN_RUNNING
+
+namespace CGAL {
+  struct Null_subdomain_index {
+    template <typename T>
+    bool operator()(const T& x) const { return 0 == x; }
+  };
+}
+
+#endif
+
+#endif //CGAL_MESH_3_NULL_SUBDOMAIN_INDEX

Mesh_3/include/CGAL/Mesh_3/features_detection/combinations.h

@@ -0,0 +1,235 @@
+// Copyright (c) 2022 GeometryFactory (France).
+// 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)     : Sebastien Loriot
+//
+//******************************************************************************
+//
+//******************************************************************************
+
+
+#ifndef CGAL_MESH_3_FEATURES_DETECTION_COMBINATIONS_H
+#define CGAL_MESH_3_FEATURES_DETECTION_COMBINATIONS_H
+
+#include <CGAL/license/Mesh_3.h>
+
+#include <array>
+
+namespace CGAL
+{
+namespace Mesh_3
+{
+namespace internal
+{
+  const std::array<std::uint8_t, 8> combinations[]
+                            = {
+                                // 1 color
+                                {0,0,0,0,0,0,0,0}, //0
+                                // 2 colors
+                                {0,0,0,0,0,0,0,1}, //1
+                                {0,0,0,0,0,0,1,1}, //2
+                                // 3 colors
+                                {0,0,0,0,0,0,1,2}, //3
+                                {0,0,0,0,0,1,1,0},
+                                {0,0,0,0,0,1,1,1},
+                                {0,0,0,0,0,1,1,2},
+                                {0,0,0,0,0,1,2,0},
+                                {0,0,0,0,0,1,2,1},
+                                {0,0,0,0,1,1,1,1},
+                                {0,0,0,0,1,1,1,2},
+                                {0,0,0,0,1,1,2,2},
+                                {0,0,0,0,1,2,2,1},
+                                {0,0,0,1,0,1,1,0},
+                                {0,0,0,1,0,1,1,1},
+                                {0,0,0,1,0,1,1,2},
+                                {0,0,0,1,0,1,2,0},
+                                {0,0,0,1,0,1,2,1},
+                                {0,0,0,1,0,1,2,2},
+                                {0,0,0,1,1,0,0,0},
+                                {0,0,0,1,1,0,0,1},
+                                {0,0,0,1,1,0,0,2},
+                                {0,0,0,1,1,0,1,1},
+                                {0,0,0,1,1,0,1,2},
+                                {0,0,0,1,1,0,2,2},
+                                {0,0,0,1,1,1,1,0},
+                                {0,0,0,1,1,1,1,2},
+                                {0,0,0,1,1,1,2,0},
+                                {0,0,0,1,1,1,2,1},
+                                {0,0,0,1,1,1,2,2},
+                                {0,0,0,1,1,2,2,0},
+                                {0,0,0,1,1,2,2,1},
+                                {0,0,0,1,1,2,2,2},
+                                {0,0,0,1,2,0,0,0},
+                                {0,0,0,1,2,0,0,1},
+                                {0,0,0,1,2,0,0,2},
+                                {0,0,0,1,2,0,1,2},
+                                {0,0,0,1,2,0,2,1},
+                                {0,0,0,1,2,1,1,0},
+                                {0,0,0,1,2,1,1,2},
+                                {0,0,0,1,2,1,2,0},
+                                {0,0,0,1,2,1,2,1},
+                                {0,0,0,1,2,1,2,2},
+                                {0,0,0,1,2,2,2,1},
+                                {0,0,1,1,1,1,0,0},
+                                {0,0,1,1,1,1,0,2},
+                                {0,0,1,1,1,1,2,2},
+                                {0,0,1,1,1,2,0,2},
+                                {0,0,1,1,1,2,2,0},
+                                {0,0,1,2,1,2,0,0},
+                                {0,0,1,2,1,2,0,1},
+                                {0,0,1,2,1,2,2,1},
+                                {0,0,1,2,2,1,0,0},
+                                {0,0,1,2,2,1,0,1},
+                                {0,1,1,0,1,0,0,1},
+                                {0,1,1,0,1,0,0,2},
+                                {0,1,1,0,1,2,2,1},
+                                {0,1,1,2,1,2,2,0}, //57
+                                // 4 colors
+                                {0,0,0,0,0,1,2,3}, //58
+                                {0,0,0,0,1,1,2,3},
+                                {0,0,0,0,1,2,2,3},
+                                {0,0,0,1,0,1,2,3},
+                                {0,0,0,1,0,2,3,0},
+                                {0,0,0,1,0,2,3,1},
+                                {0,0,0,1,1,0,2,3},
+                                {0,0,0,1,1,1,2,3},
+                                {0,0,0,1,1,2,2,3},
+                                {0,0,0,1,1,2,3,0},
+                                {0,0,0,1,1,2,3,1},
+                                {0,0,0,1,1,2,3,2},
+                                {0,0,0,1,2,0,0,3},
+                                {0,0,0,1,2,0,1,3},
+                                {0,0,0,1,2,0,2,3},
+                                {0,0,0,1,2,0,3,3},
+                                {0,0,0,1,2,1,1,3},
+                                {0,0,0,1,2,1,2,3},
+                                {0,0,0,1,2,1,3,0},
+                                {0,0,0,1,2,1,3,1},
+                                {0,0,0,1,2,1,3,2},
+                                {0,0,0,1,2,1,3,3},
+                                {0,0,0,1,2,2,2,3},
+                                {0,0,0,1,2,2,3,0},
+                                {0,0,0,1,2,2,3,1},
+                                {0,0,0,1,2,2,3,2},
+                                {0,0,0,1,2,2,3,3},
+                                {0,0,0,1,2,3,3,0},
+                                {0,0,0,1,2,3,3,1},
+                                {0,0,0,1,2,3,3,2},
+                                {0,0,0,1,2,3,3,3},
+                                {0,0,1,1,1,1,2,3},
+                                {0,0,1,1,1,2,0,3},
+                                {0,0,1,1,1,2,2,3},
+                                {0,0,1,1,1,2,3,0},
+                                {0,0,1,1,1,2,3,2},
+                                {0,0,1,1,1,2,3,3},
+                                {0,0,1,1,2,2,3,3},
+                                {0,0,1,1,2,3,2,3},
+                                {0,0,1,1,2,3,3,2},
+                                {0,0,1,2,1,2,0,3},
+                                {0,0,1,2,1,2,2,3},
+                                {0,0,1,2,1,2,3,3},
+                                {0,0,1,2,1,3,0,0},
+                                {0,0,1,2,1,3,0,1},
+                                {0,0,1,2,1,3,0,2},
+                                {0,0,1,2,1,3,2,0},
+                                {0,0,1,2,1,3,2,1},
+                                {0,0,1,2,1,3,2,3},
+                                {0,0,1,2,2,1,0,3},
+                                {0,0,1,2,2,1,1,3},
+                                {0,0,1,2,2,1,3,3},
+                                {0,0,1,2,2,3,0,0},
+                                {0,0,1,2,2,3,0,1},
+                                {0,0,1,2,2,3,0,2},
+                                {0,0,1,2,2,3,1,3},
+                                {0,0,1,2,2,3,3,1},
+                                {0,1,1,0,1,0,2,3},
+                                {0,1,1,0,1,2,2,3},
+                                {0,1,1,0,1,2,3,1},
+                                {0,1,1,0,2,3,3,2},
+                                {0,1,1,2,1,2,2,3},
+                                {0,1,1,2,1,2,3,0},
+                                {0,1,1,2,2,3,3,0},
+                                {0,1,1,2,3,0,2,3},
+                                {0,1,2,3,3,2,1,0}, //123
+                                // 5 colors
+                                {0,0,0,0,1,2,3,4}, //124
+                                {0,0,0,1,0,2,3,4},
+                                {0,0,0,1,1,2,3,4},
+                                {0,0,0,1,2,0,3,4},
+                                {0,0,0,1,2,1,3,4},
+                                {0,0,0,1,2,2,3,4},
+                                {0,0,0,1,2,3,3,4},
+                                {0,0,0,1,2,3,4,0},
+                                {0,0,0,1,2,3,4,1},
+                                {0,0,0,1,2,3,4,2},
+                                {0,0,0,1,2,3,4,3},
+                                {0,0,1,1,1,2,3,4},
+                                {0,0,1,1,2,2,3,4},
+                                {0,0,1,1,2,3,2,4},
+                                {0,0,1,1,2,3,3,4},
+                                {0,0,1,2,1,2,3,4},
+                                {0,0,1,2,1,3,0,4},
+                                {0,0,1,2,1,3,2,4},
+                                {0,0,1,2,1,3,4,0},
+                                {0,0,1,2,1,3,4,1},
+                                {0,0,1,2,1,3,4,2},
+                                {0,0,1,2,1,3,4,4},
+                                {0,0,1,2,2,1,3,4},
+                                {0,0,1,2,2,3,0,4},
+                                {0,0,1,2,2,3,1,4},
+                                {0,0,1,2,2,3,2,4},
+                                {0,0,1,2,2,3,3,4},
+                                {0,0,1,2,2,3,4,4},
+                                {0,0,1,2,3,4,0,0},
+                                {0,0,1,2,3,4,0,1},
+                                {0,0,1,2,3,4,1,4},
+                                {0,0,1,2,3,4,2,1},
+                                {0,0,1,2,3,4,2,3},
+                                {0,1,1,0,1,2,3,4},
+                                {0,1,1,0,2,3,3,4},
+                                {0,1,1,2,1,2,3,4},
+                                {0,1,1,2,1,3,4,0},
+                                {0,1,1,2,1,3,4,1},
+                                {0,1,1,2,2,0,3,4},
+                                {0,1,1,2,2,3,3,4},
+                                {0,1,1,2,2,3,4,0},
+                                {0,1,1,2,3,0,2,4},
+                                {0,1,2,3,3,2,1,4}, //166
+                                // 6 colors
+                                {0,0,0,1,2,3,4,5}, //167
+                                {0,0,1,1,2,3,4,5},
+                                {0,0,1,2,1,3,4,5},
+                                {0,0,1,2,2,3,4,5},
+                                {0,0,1,2,3,4,0,5},
+                                {0,0,1,2,3,4,1,5},
+                                {0,0,1,2,3,4,2,5},
+                                {0,0,1,2,3,4,5,5},
+                                {0,1,1,0,2,3,4,5},
+                                {0,1,1,2,1,3,4,5},
+                                {0,1,1,2,2,3,4,5},
+                                {0,1,1,2,3,0,4,5},
+                                {0,1,1,2,3,4,4,5},
+                                {0,1,1,2,3,4,5,3},
+                                {0,1,2,3,3,2,4,5}, //181
+                                // 7 colors
+                                {0,0,1,2,3,4,5,6}, //182
+                                {0,1,1,2,3,4,5,6},
+                                {0,1,2,3,3,4,5,6}, //184
+                                // 8 colors
+                                {0,1,2,3,4,5,6,7} //185
+
+                                };
+
+}//end namespace internal
+}//end namespace Mesh_3
+}//end namespace CGAL
+
+
+#endif // CGAL_MESH_3_FEATURES_DETECTION_COMBINATIONS_H

Mesh_3/include/CGAL/Mesh_3/features_detection/coordinates.h

@@ -0,0 +1,62 @@
+// Copyright (c) 2022 GeometryFactory (France).
+// 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)     : Sebastien Loriot
+//
+//******************************************************************************
+//
+//******************************************************************************
+
+
+#ifndef CGAL_MESH_3_FEATURES_DETECTION_COORDINATES_H
+#define CGAL_MESH_3_FEATURES_DETECTION_COORDINATES_H
+
+#include <CGAL/license/Mesh_3.h>
+
+#include <array>
+
+namespace CGAL
+{
+namespace Mesh_3
+{
+namespace internal
+{
+  using Coordinates = std::array<int, 3>;
+  constexpr Coordinates coordinates[8] = { {0, 0, 0},
+                                           {1, 0, 0},
+                                           {0, 1, 0},
+                                           {1, 1, 0},
+                                           {0, 0, 1},
+                                           {1, 0, 1},
+                                           {0, 1, 1},
+                                           {1, 1, 1} };
+
+  inline Coordinates minus(const Coordinates& b, const Coordinates& a) {
+    return { b[0] - a[0], b[1] - a[1], b[2] - a[2] };
+  }
+
+  inline Coordinates cross(Coordinates a, Coordinates b) {
+    return { a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0] };
+  }
+
+  inline Coordinates square(Coordinates c) {
+    return { c[0] * c[0], c[1] * c[1], c[2] * c[2] };
+  }
+
+  inline int dist(Coordinates a, Coordinates b) {
+    auto s = square(minus(b, a));
+    return s[0] + s[1] + s[2];
+  }
+
+}//end namespace internal
+}//end namespace Mesh_3
+}//end namespace CGAL
+
+#endif // CGAL_MESH_3_FEATURES_DETECTION_COORDINATES_H

Mesh_3/include/CGAL/Mesh_3/features_detection/cube_isometries.h

@@ -0,0 +1,89 @@
+// Copyright (c) 2022 GeometryFactory (France).
+// 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)     : Sebastien Loriot
+//
+//******************************************************************************
+//
+//******************************************************************************
+
+#ifndef CGAL_MESH_3_FEATURES_DETECTION_CUBE_ISOMETRIES_H
+#define CGAL_MESH_3_FEATURES_DETECTION_CUBE_ISOMETRIES_H
+
+#include <CGAL/license/Mesh_3.h>
+
+#include <array>
+
+namespace CGAL
+{
+namespace Mesh_3
+{
+namespace internal
+{
+  using Permutation = std::array<std::uint8_t, 8>;
+
+  constexpr Permutation cube_isometries[] = {
+      {0,1,2,3,4,5,6,7},
+      {1,0,3,2,5,4,7,6},
+      {4,5,0,1,6,7,2,3},
+      {5,4,1,0,7,6,3,2},
+      {6,7,4,5,2,3,0,1},
+      {7,6,5,4,3,2,1,0},
+      {2,3,6,7,0,1,4,5},
+      {3,2,7,6,1,0,5,4},
+      {1,5,3,7,0,4,2,6},
+      {5,1,7,3,4,0,6,2},
+      {5,4,7,6,1,0,3,2},
+      {4,5,6,7,0,1,2,3},
+      {4,0,6,2,5,1,7,3},
+      {0,4,2,6,1,5,3,7},
+      {1,3,0,2,5,7,4,6},
+      {3,1,2,0,7,5,6,4},
+      {3,2,1,0,7,6,5,4},
+      {2,3,0,1,6,7,4,5},
+      {2,0,3,1,6,4,7,5},
+      {0,2,1,3,4,6,5,7},
+      {1,0,5,4,3,2,7,6},
+      {0,1,4,5,2,3,6,7},
+      {7,3,5,1,6,2,4,0},
+      {3,7,1,5,2,6,0,4},
+      {7,6,3,2,5,4,1,0},
+      {6,7,2,3,4,5,0,1},
+      {2,6,0,4,3,7,1,5},
+      {6,2,4,0,7,3,5,1},
+      {4,6,5,7,0,2,1,3},
+      {6,4,7,5,2,0,3,1},
+      {7,5,6,4,3,1,2,0},
+      {5,7,4,6,1,3,0,2},
+      {0,4,1,5,2,6,3,7},
+      {4,0,5,1,6,2,7,3},
+      {3,1,7,5,2,0,6,4},
+      {1,3,5,7,0,2,4,6},
+      {5,7,1,3,4,6,0,2},
+      {7,5,3,1,6,4,2,0},
+      {3,7,2,6,1,5,0,4},
+      {7,3,6,2,5,1,4,0},
+      {0,2,4,6,1,3,5,7},
+      {2,0,6,4,3,1,7,5},
+      {5,1,4,0,7,3,6,2},
+      {1,5,0,4,3,7,2,6},
+      {6,2,7,3,4,0,5,1},
+      {2,6,3,7,0,4,1,5},
+      {6,4,2,0,7,5,3,1},
+      {4,6,0,2,5,7,1,3}
+      };
+
+  constexpr int num_isometries = 48;
+
+}//end namespace internal
+}//end namespace Mesh_3
+}//end namespace CGAL
+
+#endif // CGAL_MESH_3_FEATURES_DETECTION_CUBE_ISOMETRIES_H

Mesh_3/include/CGAL/Mesh_3/features_detection/features_detection_helpers.h

@@ -0,0 +1,72 @@
+// Copyright (c) 2022 GeometryFactory (France).
+// 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)     : Sebastien Loriot, Jane Tournois
+//
+//******************************************************************************
+//
+//******************************************************************************
+
+
+#ifndef CGAL_MESH_3_FEATURES_DETECTION_HELPERS_H
+#define CGAL_MESH_3_FEATURES_DETECTION_HELPERS_H
+
+#include <CGAL/license/Mesh_3.h>
+
+#include <string>
+
+namespace CGAL
+{
+namespace Mesh_3
+{
+namespace internal
+{
+
+  template<typename Word_type, bool b = (sizeof(Word_type) > 1)>
+  struct Color_transformation_helper
+  {
+    using type = std::unordered_map<Word_type, std::uint8_t>;
+    static void reset(type& t)
+    {
+      t.clear();
+    }
+    static bool is_valid(const type& t, const Word_type& w)
+    {
+      return t.find(w) != t.end();
+    }
+  };
+  template<typename Word_type>
+  struct Color_transformation_helper<Word_type, false>
+  {
+    using type = std::array<std::uint8_t, (1 << (sizeof(Word_type) * 8))>;
+    static void reset(type& t)
+    {
+      std::fill(t.begin(), t.end(), 8/*invalid_word*/);
+    }
+    static bool is_valid(const type& t, const Word_type& w)
+    {
+      return t[w] != 8;/*invalid_word*/
+    }
+  };
+
+  template<typename Array>
+  void debug_cerr(const char* title, const Array& tab)
+  {
+    std::cerr << title << " [";
+    for (const auto t : tab)
+      std::cout << std::to_string(t) << " ";
+    std::cout << "]" << std::endl;
+  }
+
+}//end namespace internal
+}//end namespace Mesh_3
+}//end namespace CGAL
+
+#endif // CGAL_MESH_3_FEATURES_DETECTION_HELPERS_H

Mesh_3/include/CGAL/Mesh_3/generate_label_weights.h

@@ -309,9 +309,9 @@ CGAL::Image_3 generate_label_weights_with_known_word_type(const CGAL::Image_3& i
 * @returns a `CGAL::Image_3` of weights used to build a quality `Labeled_mesh_domain_3`,
 * with the same dimensions as `image`
 */
-
+inline
 CGAL::Image_3 generate_label_weights(const CGAL::Image_3& image,
-                               const float& sigma)
+                                     const float& sigma)
 {
   CGAL_IMAGE_IO_CASE(image.image(),
     return generate_label_weights_with_known_word_type<Word>(image, sigma);

Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h

@@ -15,7 +15,6 @@
 
 #include <CGAL/license/Mesh_3.h>
 
-
 #include <vector>
 #include <map>
 #include <utility> // std::swap
@@ -23,17 +22,23 @@
 
 #include <CGAL/tuple.h>
 #include <CGAL/Image_3.h>
+#include <CGAL/number_utils.h>
+#include <CGAL/squared_distance_3.h>
+
 #include <CGAL/boost/graph/split_graph_into_polylines.h>
 #include <CGAL/Mesh_3/internal/Graph_manipulations.h>
 #include <boost/graph/adjacency_list.hpp>
-#include <CGAL/Labeled_mesh_domain_3.h> // for CGAL::Null_subdomain_index
+
-#include <CGAL/number_utils.h>
 #include <boost/utility.hpp> // for boost::prior
 #include <boost/optional.hpp>
 
 #include <CGAL/Search_traits_3.h>
 #include <CGAL/Orthogonal_incremental_neighbor_search.h>
 
+#include <CGAL/Mesh_3/Null_subdomain_index.h>
+
+#include <type_traits>
+
 namespace CGAL {
 namespace Mesh_3 {
 namespace internal {
@@ -1128,6 +1133,63 @@ polylines_to_protect(const CGAL::Image_3& cgal_image,
      existing_polylines_end);
 }
 
+
+template <typename PolylineRange1, typename PolylineRange2>
+void
+merge_and_snap_polylines(const CGAL::Image_3& image,
+                         PolylineRange1& polylines_to_snap,
+                         const PolylineRange2& existing_polylines)
+{
+  static_assert(std::is_same<typename PolylineRange1::value_type::value_type,
+                             typename PolylineRange2::value_type::value_type>::value,
+                "Polyline ranges should have same point type");
+  using P = typename PolylineRange1::value_type::value_type;
+  using K = typename Kernel_traits<P>::Kernel;
+
+  using CGAL::internal::polylines_to_protect_namespace::Vertex_info;
+  using Graph = boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS,
+                                      Vertex_info<P> >;
+  using vertex_descriptor = typename boost::graph_traits<Graph>::vertex_descriptor;
+
+  // build graph of polylines_to_snap
+  Graph graph;
+  typedef Mesh_3::internal::Returns_midpoint<K, int> Midpoint_fct;
+  Mesh_3::internal::Graph_manipulations<Graph,
+    P,
+    int,
+    Midpoint_fct> g_manip(graph);
+
+  for (const auto& polyline : polylines_to_snap)
+  {
+    if (polyline.size() < 2)
+      continue;
+
+    auto pit = polyline.begin();
+    while (boost::next(pit) != polyline.end())
+    {
+      vertex_descriptor v = g_manip.get_vertex(*pit, false);
+      vertex_descriptor w = g_manip.get_vertex(*boost::next(pit), false);
+      g_manip.try_add_edge(v, w);
+      ++pit;
+    }
+  }
+
+  // snap graph to existing_polylines
+  snap_graph_vertices(graph,
+    image.vx(), image.vy(), image.vz(),
+    boost::begin(existing_polylines), boost::end(existing_polylines),
+    K());
+
+  // rebuild polylines_to_snap
+  polylines_to_snap.clear();
+  Mesh_3::Polyline_visitor<P, Graph> visitor(polylines_to_snap, graph);
+  Less_for_Graph_vertex_descriptors<Graph> less(graph);
+  const Graph& const_graph = graph;
+  Mesh_3::Angle_tester<K> angle_tester(90.);
+  split_graph_into_polylines(const_graph, visitor, angle_tester, less);
+}
+
+
 } // namespace CGAL
 
 #endif // CGAL_MESH_3_POLYLINES_TO_PROTECT_H

Mesh_3/test/Mesh_3/CMakeLists.txt

@@ -31,6 +31,7 @@ create_single_source_cgal_program( "test_without_detect_features.cpp" )
 if(CGAL_ImageIO_USE_ZLIB)
   create_single_source_cgal_program( "test_meshing_3D_image.cpp" )
   create_single_source_cgal_program( "test_meshing_3D_image_deprecated.cpp" )
+  create_single_source_cgal_program( "test_meshing_3D_image_with_features.cpp" )
   create_single_source_cgal_program( "test_meshing_3D_gray_image.cpp" )
   create_single_source_cgal_program( "test_meshing_3D_gray_image_deprecated.cpp" )
 else()
@@ -64,6 +65,7 @@ foreach(target
     test_without_detect_features
     test_meshing_3D_image
     test_meshing_3D_image_deprecated
+    test_meshing_3D_image_with_features
     test_meshing_3D_gray_image
     test_meshing_3D_gray_image_deprecated
     test_meshing_implicit_function

Mesh_3/test/Mesh_3/test_meshing_3D_image_with_features.cpp

@@ -0,0 +1,199 @@
+// Copyright (c) 2009 INRIA Sophia-Antipolis (France).
+// 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)     : Stephane Tayeb, Jane Tournois
+//
+//******************************************************************************
+// File Description :
+//******************************************************************************
+
+#include "test_meshing_utilities.h"
+#include <CGAL/Image_3.h>
+#include <CGAL/Labeled_mesh_domain_3.h>
+#include <CGAL/Mesh_domain_with_polyline_features_3.h>
+#include <CGAL/Mesh_3/Detect_features_in_image.h>
+#include <CGAL/Mesh_3/Detect_features_on_image_bbox.h>
+#include <CGAL/use.h>
+#include <vector>
+
+#include <cstddef>
+#include <fstream>
+
+template <typename Point_3>
+bool read_polylines(const std::string fname,
+                    std::vector<std::vector<Point_3> >& polylines)
+{
+  std::ifstream ifs(fname);
+  if(ifs.bad()) return false;
+  std::size_t n;
+  while(ifs >> n) {
+    polylines.resize(polylines.size()+1);
+    std::vector<Point_3>& polyline = polylines.back();
+    while(n-- != 0) {
+      Point_3 p;
+      ifs >> p;
+      if(ifs.fail()) return false;
+      polyline.push_back(p);
+    }
+  }
+  if(ifs.bad()) return false;
+  else return ifs.eof();
+}
+
+template <typename Concurrency_tag = CGAL::Sequential_tag>
+struct Image_tester : public Tester<K_e_i>
+{
+  typedef CGAL::Image_3 Image;
+  typedef CGAL::Labeled_mesh_domain_3<K_e_i> Domain;
+  typedef CGAL::Mesh_domain_with_polyline_features_3<Domain> Mesh_domain;
+
+  typedef typename CGAL::Mesh_triangulation_3<
+    Mesh_domain,
+    CGAL::Kernel_traits<Mesh_domain>::Kernel,
+    Concurrency_tag>::type Tr;
+  typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+  typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+  void mesh_and_verify(Mesh_domain& domain, const Image& image, const double volume) const
+  {
+    namespace p = CGAL::parameters;
+    // Set mesh criteria
+    Mesh_criteria criteria(p::edge_size = 2 * image.vx(),
+                           p::facet_angle = 30,
+                           p::facet_size = 20 * image.vx(),
+                           p::facet_distance = 5 * image.vx(),
+                           p::cell_radius_edge_ratio = 3.,
+                           p::cell_size = 25 * image.vx());
+
+    // Mesh generation
+    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
+                                        CGAL::parameters::no_exude(),
+                                        CGAL::parameters::no_perturb());
+
+    c3t3.remove_isolated_vertices();
+
+    // Verify
+    this->verify_c3t3_volume(c3t3, volume * 0.95, volume * 1.05);
+    this->verify(c3t3, domain, criteria, Bissection_tag());
+
+    typedef typename Mesh_domain::Surface_patch_index Patch_id;
+    CGAL_static_assertion(CGAL::Output_rep<Patch_id>::is_specialized);
+    CGAL_USE_TYPE(Patch_id);
+  }
+
+public:
+  void image() const
+  {
+    namespace p = CGAL::parameters;
+    std::cout << "\tSeed is\t"
+      << CGAL::get_default_random().get_seed() << std::endl;
+
+    //-------------------------------------------------------
+    // Data generation
+    //-------------------------------------------------------
+    Image image;
+    image.read(CGAL::data_file_path("images/liver.inr.gz"));
+
+    Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain
+     (p::image = image,
+      p::relative_error_bound = 1e-9,
+      CGAL::parameters::p_rng = &CGAL::get_default_random(),
+      CGAL::parameters::features_detector = CGAL::Mesh_3::Detect_features_in_image());
+
+    mesh_and_verify(domain, image, 1772330.);
+  }
+
+  void image_in_bbox() const
+  {
+    namespace p = CGAL::parameters;
+    std::cout << "\tSeed is\t"
+      << CGAL::get_default_random().get_seed() << std::endl;
+
+    //-------------------------------------------------------
+    // Data generation
+    //-------------------------------------------------------
+    Image image;
+    image.read(CGAL::data_file_path("images/40420.inr"));
+
+    Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain
+     (p::image = image,
+      p::relative_error_bound = 1e-9,
+      CGAL::parameters::p_rng = &CGAL::get_default_random(),
+      CGAL::parameters::features_detector = CGAL::Mesh_3::Detect_features_on_image_bbox());
+
+    mesh_and_verify(domain, image, 625044.);
+  }
+
+  void image_with_input_features() const
+  {
+    namespace p = CGAL::parameters;
+    std::cout << "\tSeed is\t"
+      << CGAL::get_default_random().get_seed() << std::endl;
+
+    //-------------------------------------------------------
+    // Data generation
+    //-------------------------------------------------------
+    Image image;
+    image.read(CGAL::data_file_path("images/40420.inr"));
+
+    const std::string lines_fname = CGAL::data_file_path("images/420.polylines.txt");
+    using Point_3 = Domain::Point_3;
+    std::vector<std::vector<Point_3> > features_input;
+    if (!read_polylines(lines_fname, features_input)) // see file "read_polylines.h"
+      assert(false);
+
+    Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain
+     (p::image = image,
+      p::relative_error_bound = 1e-9,
+      CGAL::parameters::p_rng = &CGAL::get_default_random(),
+      CGAL::parameters::features_detector = CGAL::Mesh_3::Detect_features_on_image_bbox(),
+      CGAL::parameters::input_features = std::cref(features_input));
+
+    mesh_and_verify(domain, image, 632091.);
+  }
+};
+
+
+
+int main()
+{
+  Image_tester<> test_epic;
+  std::cerr << "Mesh generation from a 3D image"
+            << " with detection of triple lines:\n";
+  test_epic.image();
+
+  std::cerr << "Mesh generation from a 3D image"
+            << " with detection of triple lines on bbox only:\n";
+  test_epic.image_in_bbox();
+
+  std::cerr << "Mesh generation from a 3D image"
+            << " with detection of triple lines on bbox"
+            << " and input feature polylines:\n";
+  test_epic.image_with_input_features();
+
+#ifdef CGAL_LINKED_WITH_TBB
+  Image_tester<CGAL::Parallel_tag> test_epic_p;
+  std::cerr << "Parallel mesh generation from a 3D image"
+            << " with detection of triple lines:\n";
+  test_epic_p.image();
+
+  std::cerr << "Parallel mesh generation from a 3D image"
+            << " with detection of triple lines on bbox only:\n";
+  test_epic_p.image_in_bbox();
+
+  std::cerr << "Parallel mesh generation from a 3D image"
+            << " with detection of triple lines on bbox"
+            << " and input feature polylines:\n";
+  test_epic_p.image_with_input_features();
+#endif
+
+  return EXIT_SUCCESS;
+}

Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp

@@ -277,12 +277,16 @@ void Mesh_3_plugin::mesh_3_volume()
   mesh_3(Mesh_type::VOLUME);
 }
 
-boost::optional<QString> Mesh_3_plugin::get_items_or_return_error_string() const {
+boost::optional<QString> Mesh_3_plugin::get_items_or_return_error_string() const
+{
   using boost::get;
   items = {};
   features_protection_available = false;
   item = nullptr;
-  for (int ind : scene->selectionIndices()) {
+  Scene_polylines_item* polylines_item = nullptr;
+
+  for (int ind : scene->selectionIndices())
+  {
     try {
       if (auto sm_item =
               qobject_cast<Scene_surface_mesh_item*>(scene->item(ind))) {
@@ -313,9 +317,11 @@ boost::optional<QString> Mesh_3_plugin::get_items_or_return_error_string() const
           return tr("An image items cannot be mixed with other items type");
       }
 #  endif
-      else if (auto polylines_item =
+      else if ((polylines_item =
-                   qobject_cast<Scene_polylines_item*>(scene->item(ind))) {
+                qobject_cast<Scene_polylines_item*>(scene->item(ind))))
-        if (!items) items = Polyhedral_mesh_items{};
+      {
+        if (!items)
+          continue;
         auto poly_items_ptr = get<Polyhedral_mesh_items>(&*items);
         if(poly_items_ptr) {
           if (poly_items_ptr->polylines_item) {
@@ -323,12 +329,16 @@ boost::optional<QString> Mesh_3_plugin::get_items_or_return_error_string() const
           } else {
             poly_items_ptr->polylines_item = polylines_item;
           }
-        } else {
+        }
-          auto image_items = get<Image_mesh_items>(*items);
+        else {
-          if (image_items.polylines_item) {
+          if(auto image_items_ptr = get<Image_mesh_items>(&*items))
-            return tr("Only one polyline item is accepted");
+          {
-          } else {
+            if (image_items_ptr->polylines_item) {
-            image_items.polylines_item = polylines_item;
+              return tr("Only one polyline item is accepted");
+            }
+            else {
+              image_items_ptr->polylines_item = polylines_item;
+            }
           }
         }
       }
@@ -341,6 +351,20 @@ boost::optional<QString> Mesh_3_plugin::get_items_or_return_error_string() const
       }
     } catch (const boost::bad_get&) { return tr("Wrong selection of items"); }
   } // end for loop on selected items
+
+  //attach polylines_item to one or the other item
+  //if it could not be done in the for loop
+  //because of selection order
+  if (polylines_item != nullptr && items != boost::none)
+  {
+    auto poly_items_ptr = get<Polyhedral_mesh_items>(&*items);
+    auto image_items_ptr = get<Image_mesh_items>(&*items);
+    if(poly_items_ptr && poly_items_ptr == nullptr)
+      poly_items_ptr->polylines_item = polylines_item;
+    else if(image_items_ptr && image_items_ptr == nullptr)
+      image_items_ptr->polylines_item = polylines_item;
+  }
+
   if (!items) { return tr("Selected objects can't be meshed"); }
   item = nullptr;
   features_protection_available = false;
@@ -428,6 +452,8 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type,
       get<Polyhedral_mesh_items>(&*items)
           ? get<Polyhedral_mesh_items>(&*items)->polylines_item
           : nullptr;
+  if (polylines_item == nullptr && get<Image_mesh_items>(&*items) != nullptr)
+    polylines_item = get<Image_mesh_items>(&*items)->polylines_item;
   Scene_implicit_function_item* function_item =
       get<Implicit_mesh_items>(&*items)
           ? get<Implicit_mesh_items>(&*items)->function_item.get()
@@ -571,18 +597,28 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type,
   ui.edgeLabel->setEnabled(ui.noEdgeSizing->isChecked());
   ui.edgeSizing->setEnabled(ui.noEdgeSizing->isChecked());
 
+  const QString sharp_and_boundary("Sharp and Boundary edges");
+  const QString boundary_only("Boundary edges only");
+  const QString sharp_edges("Sharp edges");
+  const QString input_polylines("Input polylines");
+  const QString on_cube("Polylines on cube");
+  const QString triple_lines("Triple+ lines");
   if (features_protection_available) {
     if (items->which() == POLYHEDRAL_MESH_ITEMS) {
       if (mesh_type == Mesh_type::SURFACE_ONLY) {
-        ui.protectEdges->addItem(QString("Sharp and Boundary edges"));
+        ui.protectEdges->addItem(sharp_and_boundary);
-        ui.protectEdges->addItem(QString("Boundary edges only"));
+        ui.protectEdges->addItem(boundary_only);
       } else
-        ui.protectEdges->addItem(QString("Sharp edges"));
+        ui.protectEdges->addItem(sharp_edges);
     } else if (items->which() == IMAGE_MESH_ITEMS) {
-      if (polylines_item != nullptr)
+      if (polylines_item != nullptr) {
-        ui.protectEdges->addItem(QString("Input polylines"));
+        ui.protectEdges->addItem(QString(input_polylines).append(" only"));
+        ui.protectEdges->addItem(QString(on_cube).append(" and input polylines"));
+        ui.protectEdges->addItem(QString(triple_lines).append(" and input polylines"));
+      }
       else {
-        ui.protectEdges->addItem(QString("Polylines on cube"));
+        ui.protectEdges->addItem(on_cube);
+        ui.protectEdges->addItem(triple_lines);
       }
     }
   }
@@ -629,10 +665,15 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type,
   approx = !ui.noApprox->isChecked() ? 0 : ui.approx->value();
   tets_shape = !ui.noTetShape->isChecked() ? 0 : ui.tetShape->value();
   tets_sizing = !ui.noTetSizing->isChecked() ? 0 : ui.tetSizing->value();
-  protect_features =
+
-      ui.protect->isChecked() && (ui.protectEdges->currentIndex() == 0);
+  const int pe_ci = ui.protectEdges->currentIndex();
-  protect_borders =
+  protect_borders = ui.protect->isChecked()
-      ui.protect->isChecked() && (ui.protectEdges->currentIndex() == 1);
+    && (  pe_ci == ui.protectEdges->findText(on_cube, Qt::MatchContains)
+       || pe_ci == ui.protectEdges->findText(boundary_only, Qt::MatchContains));
+  protect_features = ui.protect->isChecked()
+    && (  pe_ci == ui.protectEdges->findText(triple_lines, Qt::MatchContains)
+       || pe_ci == ui.protectEdges->findText(sharp_and_boundary, Qt::MatchContains));
+
   const bool detect_connected_components = ui.detectComponents->isChecked();
   const int manifold = (ui.manifoldCheckBox->isChecked() ? 1 : 0) +
                        (ui.facetTopology->isChecked() ? 2 : 0);
@@ -728,7 +769,7 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type,
      }
     break;
   }//end case POLYHEDRAL_MESH_ITEMS
-  // Image
+  // Implicit functions
 #  ifdef CGAL_MESH_3_DEMO_ACTIVATE_IMPLICIT_FUNCTIONS
   case IMPLICIT_MESH_ITEMS: {
     const Implicit_function_interface* pFunction = function_item->function();
@@ -747,11 +788,15 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type,
                               manifold,
                               mesh_type == Mesh_type::SURFACE_ONLY);
     break;
-  }
+  }//end case IMPLICIT_MESH_ITEMS
 #  endif
+  // Images
 #  ifdef CGAL_MESH_3_DEMO_ACTIVATE_SEGMENTED_IMAGES
   case IMAGE_MESH_ITEMS: {
     const Image* pImage = image_item->image();
+    auto& image_items = get<Image_mesh_items>(*items);
+    const auto img_polylines_item = image_items.polylines_item;
+
     if (nullptr == pImage) {
       QMessageBox::critical(mw, tr(""), tr("ERROR: no data in selected item"));
       return;
@@ -773,7 +818,7 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type,
 
     thread = cgal_code_mesh_3(
         pImage,
-        (polylines_item == nullptr) ? plc : polylines_item->polylines,
+        (img_polylines_item == nullptr) ? plc : img_polylines_item->polylines,
         angle,
         facets_sizing,
         approx,
@@ -781,6 +826,7 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type,
         edges_sizing,
         tets_shape,
         protect_features,
+        protect_borders,
         manifold,
         mesh_type == Mesh_type::SURFACE_ONLY,
         detect_connected_components,

Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp

@@ -13,6 +13,8 @@
 #include "Mesh_function.h"
 #include "Facet_extra_criterion.h"
 
+#include <CGAL/Mesh_3/Detect_features_in_image.h>
+#include <CGAL/Mesh_3/Detect_features_on_image_bbox.h>
 
 using namespace CGAL::Three;
 
@@ -291,7 +293,8 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage,
                                  const double tet_sizing,
                                  const double edge_size,
                                  const double tet_shape,
-                                 bool protect_features,
+                                 bool protect_features, //detect_polylines
+                                 const bool protect_borders,//polylines on bbox
                                  const int manifold,
                                  const bool surface_only,
                                  bool detect_connected_components,
@@ -303,11 +306,9 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage,
 {
   if (nullptr == pImage) { return nullptr; }
 
-  if(! polylines.empty()){
-    protect_features = true; // so that it will be passed in make_mesh_3
-  }
   Mesh_parameters param;
-  param.protect_features = protect_features;
+  param.protect_features
+    = protect_features || protect_borders || !polylines.empty();
   param.detect_connected_components = detect_connected_components;
   param.facet_angle = facet_angle;
   param.facet_sizing = facet_sizing;
@@ -343,7 +344,7 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage,
   {
     namespace p = CGAL::parameters;
 
-    Image_mesh_domain* p_domain;
+    Image_mesh_domain* p_domain = nullptr;
 #ifdef CGAL_USE_ITK
     if(nullptr != pWeights)
     {
@@ -359,35 +360,57 @@ Meshing_thread* cgal_code_mesh_3(const Image* pImage,
     }
     else
 #endif
-    {
+      if (protect_features)
-      p_domain = new Image_mesh_domain
+      {
-      (Image_mesh_domain::create_labeled_image_mesh_domain
+        p_domain = new Image_mesh_domain
-       (p::image = *pImage,
+          (Image_mesh_domain::create_labeled_image_mesh_domain
-        p::relative_error_bound = 1e-6,
+          (p::image = *pImage,
-        p::construct_surface_patch_index =
+            p::relative_error_bound = 1e-6,
-          [](int i, int j) { return (i * 1000 + j); }
+            p::construct_surface_patch_index =
-        )
+            [](int i, int j) { return (i * 1000 + j); },
-       );
+            p::features_detector = CGAL::Mesh_3::Detect_features_in_image(),
-    }
+            p::input_features = std::cref(polylines)
+          )
+        );
+      }
+      else if (protect_borders)//protect polylines on image Bbox
+      {
+        p_domain = new Image_mesh_domain
+          (Image_mesh_domain::create_labeled_image_mesh_domain
+          (p::image = *pImage,
+            p::relative_error_bound = 1e-6,
+            p::construct_surface_patch_index =
+            [](int i, int j) { return (i * 1000 + j); },
+            p::features_detector = CGAL::Mesh_3::Detect_features_on_image_bbox(),
+            p::input_features = std::cref(polylines)
+          )
+        );
+      }
+      else if (!polylines.empty())
+      {
+        p_domain = new Image_mesh_domain
+          (Image_mesh_domain::create_labeled_image_mesh_domain
+          (p::image = *pImage,
+            p::relative_error_bound = 1e-6,
+            p::construct_surface_patch_index =
+            [](int i, int j) { return (i * 1000 + j); },
+            p::input_features = std::cref(polylines)
+          )
+        );
+      }
 
-    if(protect_features && polylines.empty()){
+      if (p_domain == nullptr)
-      std::vector<std::vector<Bare_point> > polylines_on_bbox;
+      {
+        p_domain = new Image_mesh_domain
+          (Image_mesh_domain::create_labeled_image_mesh_domain
+          (p::image = *pImage,
+            p::relative_error_bound = 1e-6,
+            p::construct_surface_patch_index =
+            [](int i, int j) { return (i * 1000 + j); }
+          )
+        );
+      }
 
-      CGAL_IMAGE_IO_CASE(pImage->image(),
-                         {
-                           typedef Word Image_word_type;
-                           (CGAL::polylines_to_protect<
-                              Bare_point,
-                              Image_word_type>(*pImage, polylines_on_bbox));
-                           p_domain->add_features(polylines_on_bbox.begin(),
-                                                  polylines_on_bbox.end());
-                         }
-                         );
-    }
-    if(! polylines.empty()){
-      // Insert edge in domain
-      p_domain->add_features(polylines.begin(), polylines.end());
-    }
     typedef ::Mesh_function<Image_mesh_domain,
                             Mesh_fnt::Labeled_image_domain_tag> Mesh_function;
     Mesh_function* p_mesh_function = new Mesh_function(p_new_item->c3t3(),

Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.h

@@ -73,6 +73,7 @@ Meshing_thread* cgal_code_mesh_3(const CGAL::Image_3* pImage,
                                  const double edge_size,
                                  const double tet_shape,
                                  bool protect_features,
+                                 const bool protect_borders,
                                  const int manifold,
                                  const bool surface_only,
                                  bool detect_connected_components,

SMDS_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h

@@ -640,7 +640,7 @@ public:
     but are isolated from the complex at the end of the meshing process.
 
     This function removes these so-called \em isolated vertices, that belong to the
-    triangulation but not to any cell of the `C3T3`, from the triangulation.
+    triangulation but not to any simplex of the `C3T3`, from the triangulation.
   */
   void remove_isolated_vertices()
   {
@@ -669,7 +669,8 @@ public:
     std::vector<Vertex_handle> isolated;
     for (Vertex_handle v : tr.finite_vertex_handles())
     {
-      if (v->meshing_info() == 0.)
+      if (v->meshing_info() == 0.
+        && (v->in_dimension() > 1 || v->in_dimension() < 0))
         isolated.push_back(v);
     }
 

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -325,6 +325,8 @@ CGAL_add_named_parameter_with_compatibility(rng_t, rng, p_rng)
 CGAL_add_named_parameter_with_compatibility(null_subdomain_index_param_t,null_subdomain_index_param, null_subdomain_index)
 CGAL_add_named_parameter_with_compatibility(surface_patch_index_t, surface_patch_index, construct_surface_patch_index)
 CGAL_add_named_parameter_with_compatibility_ref_only(weights_param_t, weights_param, weights)
+CGAL_add_named_parameter_with_compatibility(features_detector_param_t, features_detector_param, features_detector)
+CGAL_add_named_parameter_with_compatibility(input_features_param_t, input_features_param, input_features)
 
 CGAL_add_named_parameter_with_compatibility(edge_size_param_t, edge_size_param, edge_size)
 CGAL_add_named_parameter_with_compatibility_ref_only(edge_sizing_field_param_t, edge_sizing_field_param, edge_sizing_field)