diff --git a/Lab/demo/Lab/Plugins/Mesh_3/Mesh_function.h b/Lab/demo/Lab/Plugins/Mesh_3/Mesh_function.h
index cf2ffaf9caa..99552b05b7f 100644
--- a/Lab/demo/Lab/Plugins/Mesh_3/Mesh_function.h
+++ b/Lab/demo/Lab/Plugins/Mesh_3/Mesh_function.h
@@ -257,9 +257,7 @@ initialize(const Mesh_criteria& criteria, Mesh_fnt::Labeled_image_domain_tag)
        p_.protect_features,
        p::mesh_3_options(p::pointer_to_stop_atomic_boolean = &stop_,
                          p::nonlinear_growth_of_balls = true).v,
-       p::internal::Initial_points_generator_generator<D_, C3t3>()
-          (p::initial_points_generator(
-            CGAL::Construct_initial_points_labeled_image(*p_.image_3_ptr)).v));
+       CGAL::Construct_initial_points_labeled_image<C3t3, Domain>(*p_.image_3_ptr, *domain_));
   }
   else
   {
@@ -278,6 +276,11 @@ initialize(const Mesh_criteria& criteria, Mesh_fnt::Gray_image_domain_tag)
   // features, or with the initial points (or both).
   if (p_.detect_connected_components)
   {
+    CGAL::Construct_initial_points_gray_image<C3t3, Domain, Compare_to_isovalue> generator
+       (*p_.image_3_ptr,
+        *domain_,
+        p_.iso_value,
+        Compare_to_isovalue(p_.iso_value, p_.inside_is_less));
     CGAL::Mesh_3::internal::C3t3_initializer<
         C3t3,
         Domain,
@@ -289,9 +292,7 @@ initialize(const Mesh_criteria& criteria, Mesh_fnt::Gray_image_domain_tag)
          p_.protect_features,
          p::mesh_3_options(p::pointer_to_stop_atomic_boolean = &stop_,
                            p::nonlinear_growth_of_balls = true).v,
-         p::internal::Initial_points_generator_generator<D_, C3t3>()
-           (p::initial_points_generator(
-             CGAL::Construct_initial_points_gray_image(*p_.image_3_ptr, p_.iso_value, Compare_to_isovalue(p_.iso_value, p_.inside_is_less))).v));
+         generator);
   }
   else
   {
diff --git a/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_gray_image.h b/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_gray_image.h
index 19ae3698780..9264d44d844 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_gray_image.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Construct_initial_points_gray_image.h
@@ -39,19 +39,21 @@ namespace CGAL
  * \sa `CGAL::make_mesh_3()`
  * \sa `CGAL::Construct_initial_points_labeled_image`
  */
-template <typename Functor = CGAL::Null_functor>
+template <typename C3t3, typename MeshDomain, typename Functor = CGAL::Null_functor>
 struct Construct_initial_points_gray_image
 {
   const CGAL::Image_3 & image_;
-  const double iso_value_;
+  const MeshDomain& domain_;
+  const typename MeshDomain::R::FT iso_value_;
   Functor image_values_to_subdomain_indices_;
 
-  template <typename FT>
   Construct_initial_points_gray_image(const CGAL::Image_3 & image,
-         const FT& iso_value,
+         const MeshDomain& domain,
+         const double iso_value,
          const Functor image_values_to_subdomain_indices = CGAL::Null_functor())
       : image_(image)
-      , iso_value_(iso_value)
+      , domain_(domain)
+      , iso_value_(static_cast<typename MeshDomain::R::FT>(iso_value))
       , image_values_to_subdomain_indices_(image_values_to_subdomain_indices)
   { }
 
@@ -67,15 +69,16 @@ struct Construct_initial_points_gray_image
   * \tparam C3t3 model of `MeshComplex_3InTriangulation_3`
   *
   */
-  template <typename OutputIterator, typename MeshDomain, typename C3t3>
-  OutputIterator operator()(OutputIterator pts, const MeshDomain& domain, const C3t3& c3t3, int n = 20) const
+  template <typename OutputIterator>
+  OutputIterator operator()(OutputIterator pts, int n = 20) const
   {
     using CGAL::Mesh_3::internal::Create_gray_image_values_to_subdomain_indices;
     typedef Create_gray_image_values_to_subdomain_indices<Functor> C_i_v_t_s_i;
     typedef typename C_i_v_t_s_i::type Image_values_to_subdomain_indices;
     Image_values_to_subdomain_indices transform_fct =
       C_i_v_t_s_i()(image_values_to_subdomain_indices_, iso_value_);
-    Construct_initial_points_labeled_image(image_).operator()(pts, domain, transform_fct, c3t3, n);
+    Construct_initial_points_labeled_image<C3t3, MeshDomain> init_pts{ image_, domain_ };
+    init_pts(pts, transform_fct, n);
     return pts;
   }
 };
diff --git a/Mesh_3/include/CGAL/make_mesh_3.h b/Mesh_3/include/CGAL/make_mesh_3.h
index 4cc6827311e..6e38a36bc7b 100644
--- a/Mesh_3/include/CGAL/make_mesh_3.h
+++ b/Mesh_3/include/CGAL/make_mesh_3.h
@@ -50,7 +50,12 @@ add_points_from_generator(C3T3& c3t3,
   typedef typename C3T3::Vertex_handle Vertex_handle;
   typedef CGAL::Mesh_3::Triangulation_helpers<Tr> Th;
 
-  using PointDimIndex = parameters::internal::Initial_point_type<MeshDomain, C3T3>;
+  struct PointDimIndex
+  {
+    typename Tr::Point m_wpt;
+    int m_dim;
+    typename MeshDomain::Index m_index;
+  };
 
   const auto& cwp = c3t3.triangulation().geom_traits().construct_weighted_point_3_object();
 
@@ -99,11 +104,18 @@ init_c3t3(C3T3& c3t3, const MeshDomain& domain, const MeshCriteria&,
           const int nb_initial_points,
           const InitializationOptions& init_options)
 {
-  add_points_from_generator(c3t3, domain, nb_initial_points, init_options.generator());
+  add_points_from_generator(c3t3, domain, nb_initial_points, init_options);
+
+  typedef CGAL::parameters::internal::Initialization_options<MeshDomain, C3T3> Default_init_options;
+  bool is_default_init = false;
+  if constexpr (std::is_same_v<InitializationOptions, Default_init_options>)
+  {
+    is_default_init = init_options.is_default();
+  }
 
   // If c3t3 initialization is not sufficient (may happen if
   // the user has not specified enough points ), add some surface points
-  bool need_more_init = c3t3.triangulation().dimension() != 3 || !init_options.is_default();
+  bool need_more_init = c3t3.triangulation().dimension() != 3 || !is_default_init;
   if(!need_more_init)
   {
     CGAL::Mesh_3::C3T3_helpers<C3T3, MeshDomain> helper(c3t3, domain);
@@ -287,8 +299,13 @@ struct C3t3_initializer < C3T3, MD, MC, true, CGAL::Tag_true>
 
       // If c3t3 initialization is not sufficient (may happen if there is only
       // a planar curve as feature for example), add some surface points
+      bool is_default_init = false;
+      if constexpr (std::is_same_v<InitOptions, Default_init_options>)
+      {
+        is_default_init = init_options.is_default();
+      }
 
-      bool need_more_init = c3t3.triangulation().dimension() != 3 || !init_options.is_default();
+      bool need_more_init = c3t3.triangulation().dimension() != 3 || !is_default_init;
       if(!need_more_init) {
         CGAL::Mesh_3::C3T3_helpers<C3T3, MD> helper(c3t3, domain);
         helper.update_restricted_facets();
@@ -537,7 +554,7 @@ C3T3 make_mesh_3(const MeshDomain& domain, const MeshCriteria& criteria, const C
     parameters::internal::Manifold_options manifold_options_param = choose_parameter(get_parameter(np, internal_np::manifold_param), parameters::internal::Manifold_options());
 
     // range of initial points
-    using Initial_point = parameters::internal::Initial_point_type<MeshDomain, C3T3>;
+    using Initial_point = std::pair<typename MeshDomain::Point_3, typename MeshDomain::Index>;
     using Initial_points_range_ref = typename internal_np::Lookup_named_param_def<internal_np::initial_points_param_t,
                                                                               CGAL_NP_CLASS,
                                                                               std::vector<Initial_point>>::reference;
@@ -609,8 +626,6 @@ void make_mesh_3_impl(C3T3& c3t3,
   CGAL::get_default_random() = CGAL::Random(0);
 #endif
 
-  using Init_options = parameters::internal::Initialization_options<MeshDomain, C3T3, InitPtsGenerator, InitPtsVec>;
-
   // Initialize c3t3
   Mesh_3::internal::C3t3_initializer<
     C3T3,
diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/mesh_option_classes.h b/STL_Extension/include/CGAL/STL_Extension/internal/mesh_option_classes.h
index 180dd59b5f6..8505194049c 100644
--- a/STL_Extension/include/CGAL/STL_Extension/internal/mesh_option_classes.h
+++ b/STL_Extension/include/CGAL/STL_Extension/internal/mesh_option_classes.h
@@ -14,6 +14,7 @@
 #include <functional>
 
 #include <CGAL/STL_Extension/internal/Has_features.h>
+#include <CGAL/Default.h>
 #include <boost/iterator/function_output_iterator.hpp>
 #include <type_traits>
 #include <iterator>
@@ -172,15 +173,6 @@ private:
 };
 
 // Mesh initialization
-
-template<typename MeshDomain, typename C3t3>
-struct Initial_point_type
-{
-  typename C3t3::Triangulation::Point m_weighted_point;
-  int m_dimension;
-  typename MeshDomain::Index m_index;
-};
-
 struct Dummy_initial_points_generator
 {
   template<typename OutputIterator>
@@ -191,20 +183,27 @@ struct Dummy_initial_points_generator
 // without knowing their types, into a single generator.
 template <typename MeshDomain,
           typename C3t3,
-          typename InitialPointsGenerator = Dummy_initial_points_generator,
-          typename InitialPointsRange = std::vector<Initial_point_type<MeshDomain, C3t3> > >
+          typename InitialPointsGenerator = CGAL::Default,
+          typename InitialPointsRange = CGAL::Default
+>
 struct Initialization_options
 {
-  using DefaultGenerator = Dummy_initial_points_generator;
-  using Initial_points_const_iterator = typename InitialPointsRange::const_iterator;
-  using Initial_point = typename std::iterator_traits<Initial_points_const_iterator>::value_type;
+  using Default_generator = Dummy_initial_points_generator;
+  using Initial_points_generator
+          = typename CGAL::Default::Get<InitialPointsGenerator, Default_generator>::type;
+  using Default_initial_point_type
+          = std::tuple<typename C3t3::Triangulation::Point, int, typename MeshDomain::Index>;
+  using Initial_points_range
+          = typename CGAL::Default::Get<InitialPointsRange, std::vector<Default_initial_point_type>>::type;
 
+  using Initial_points_const_iterator = typename Initial_points_range::const_iterator;
+  using Initial_point = typename std::iterator_traits<Initial_points_const_iterator>::value_type;
 
   Initialization_options()
   {}
 
-  Initialization_options(const InitialPointsGenerator& generator,
-                         const InitialPointsRange& initial_points = InitialPointsRange())
+  Initialization_options(const Initial_points_generator& generator,
+                         const Initial_points_range& initial_points = Initial_points_range())
     : initial_points_generator_(generator)
     , begin_it(initial_points.begin())
     , end_it(initial_points.end())
@@ -214,13 +213,13 @@ struct Initialization_options
   OutputIterator operator()(OutputIterator pts, int n = 0) const
   {
     // add initial_points
-    for (typename InitialPointsRange::const_iterator it = begin_it; it != end_it; ++it)
+    for (Initial_points_const_iterator it = begin_it; it != end_it; ++it)
       *pts++ = *it;
 
     return initial_points_generator_(pts, n);
   }
 
-  const InitialPointsGenerator& generator() const
+  const Initial_points_generator& generator() const
   {
     return initial_points_generator_;
   }
@@ -228,11 +227,11 @@ struct Initialization_options
   bool is_default() const
   {
     return begin_it == end_it
-        && std::is_same_v<InitialPointsGenerator, DefaultGenerator>;
+        && std::is_same_v<InitialPointsGenerator, Default_generator>;
   }
 
 private:
-  InitialPointsGenerator initial_points_generator_;
+  Initial_points_generator initial_points_generator_;
 
   // The two iterators point to the `initial_points` container
   Initial_points_const_iterator begin_it;
@@ -349,7 +348,7 @@ struct Initial_points_generator_generator
     if (input_points.empty()) {
       return Initialization_options(Initial_points_generator_domain_traductor(), input_points, true);
     }
-    return Initialization_options(Initial_points_generator_empty(), input_features, true);
+    return Initialization_options(Initial_points_generator_empty(), input_points, true);
   }
 
   // Default construction