diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h b/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h
index 0a3d73f0341e..3b41f515bb5b 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h
@@ -32,7 +32,6 @@
 #include <CGAL/Time_stamper.h>
 
 #include <CGAL/iterator.h>
-#include <CGAL/tuple.h>
 
 #include <CGAL/Mesh_3/Concurrent_mesher_config.h>
 
@@ -74,8 +73,13 @@ class Mesh_global_optimizer_base
   // The sizing field info is stored inside the move vector because it is computed
   // when the move is computed. This is because the parallel version uses the threadsafe
   // version of incident_cells (which thus requires points to not be moving yet)
-  typedef std::vector<std::tuple<typename Tr::Vertex_handle, Vector_3, FT> >
-                                                            Moves_vector;
+  struct Move
+  {
+    typename Tr::Vertex_handle vertex_;
+    Vector_3 move_;
+    FT size_;
+  };
+  typedef std::vector<Move>                                 Moves_vector;
   typedef unsigned int                                      Nb_frozen_points_type;
 
   Mesh_global_optimizer_base(const Bbox_3 &, int)
@@ -129,8 +133,14 @@ class Mesh_global_optimizer_base<Tr, Parallel_tag>
   typedef typename GT::FT                                   FT;
   typedef typename GT::Vector_3                             Vector_3;
   typedef typename Tr::Lock_data_structure                  Lock_data_structure;
-  typedef tbb::concurrent_vector<std::tuple<typename Tr::Vertex_handle, Vector_3, FT> >
-                                                            Moves_vector;
+
+  struct Move
+  {
+    typename Tr::Vertex_handle vertex_;
+    Vector_3 move_;
+    FT size_;
+  };
+  typedef tbb::concurrent_vector<Move>                      Moves_vector;
   typedef std::atomic<unsigned int>                         Nb_frozen_points_type ;
 
   Mesh_global_optimizer_base(const Bbox_3 &bbox, int num_grid_cells_per_axis)
@@ -411,6 +421,7 @@ class Mesh_global_optimizer
     {
       typename GT::Construct_point_3 cp = m_gt.construct_point_3_object();
       typename GT::Construct_translated_point_3 translate = m_gt.construct_translated_point_3_object();
+      typedef typename Moves_vector_::value_type Move;
 
       Vector_3 move = m_mgo.compute_move(oldv);
       if ( CGAL::NULL_VECTOR != move )
@@ -429,7 +440,7 @@ class Mesh_global_optimizer
         //note : this is not happening for Lloyd and ODT so it's commented
         //       maybe for a new global optimizer it should be de-commented
 
-        m_moves.push_back(std::make_tuple(oldv, move, size));
+        m_moves.push_back(Move{oldv, move, size});
       }
       else // CGAL::NULL_VECTOR == move
       {
@@ -517,13 +528,13 @@ class Mesh_global_optimizer
     {
       for( size_t i = r.begin() ; i != r.end() ; ++i)
       {
-        const Vertex_handle& v = std::get<0>(m_moves[i]);
-        const Vector_3& move = std::get<1>(m_moves[i]);
+        const Vertex_handle& v = m_moves[i].vertex_;
+        const Vector_3& move = m_moves[i].move_;
 
         // Get size at new position
         if ( MGO::Sizing_field::is_vertex_update_needed )
         {
-          FT size = std::get<2>(m_moves[i]);
+          const FT size = m_moves[i].size_;
 
           // Move point
           bool could_lock_zone;
@@ -846,7 +857,7 @@ compute_moves(Moving_vertices_set& moving_vertices)
           size = sizing_field_(new_position, oldv);
         }
 
-        moves.push_back(std::make_tuple(oldv, move, size));
+        moves.push_back({oldv, move, size});
       }
       else // CGAL::NULL_VECTOR == move
       {
@@ -965,12 +976,12 @@ update_mesh(const Moves_vector& moves,
          it != moves.end() ;
          ++it )
     {
-      const Vertex_handle& v = std::get<0>(*it);
-      const Vector_3& move = std::get<1>(*it);
+      const Vertex_handle& v = it->vertex_;
+      const Vector_3& move = it->move_;
       // Get size at new position
       if ( Sizing_field::is_vertex_update_needed )
       {
-        FT size = std::get<2>(*it);
+        const FT size = it->size_;
 
 #ifdef CGAL_MESH_3_OPTIMIZER_VERY_VERBOSE
         std::cerr << "Moving #" << it - moves.begin()
diff --git a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
index e84e7dd933fe..aa7d3dfd6af9 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
@@ -54,7 +54,6 @@
 #  include <boost/math/special_functions/next.hpp> // for float_prior
 #endif
 #include <optional>
-#include <boost/tuple/tuple.hpp>
 #include <boost/unordered_map.hpp>
 #include <boost/unordered_set.hpp>
 
@@ -1019,29 +1018,30 @@ Protect_edges_sizing_field<C3T3, MD, Sf, Df>::
 insert_balls_on_edges()
 {
   // Get features
-  typedef std::tuple<Curve_index,
-                             std::pair<Bare_point,Index>,
-                             std::pair<Bare_point,Index> >    Feature_tuple;
-  typedef std::vector<Feature_tuple>                          Input_features;
-
-  Input_features input_features;
+  struct Feature_tuple
+  {
+    Curve_index curve_index_;
+    std::pair<Bare_point, Index> point_s_;
+    std::pair<Bare_point, Index> point_t_;
+  };
+  std::vector<Feature_tuple> input_features;
   domain_.get_curves(std::back_inserter(input_features));
 
   // Iterate on edges
   for (const Feature_tuple& ft : input_features)
   {
     if(forced_stop()) break;
-    const Curve_index& curve_index = std::get<0>(ft);
+    const Curve_index& curve_index = ft.curve_index_;
     if ( ! is_treated(curve_index) )
     {
 #if CGAL_MESH_3_PROTECTION_DEBUG & 1
       std::cerr << "\n** treat curve #" << curve_index << std::endl;
 #endif
-      const Bare_point& p = std::get<1>(ft).first;
-      const Bare_point& q = std::get<2>(ft).first;
+      const Bare_point& p = ft.point_s_.first;
+      const Bare_point& q = ft.point_t_.first;
 
-      const Index& p_index = std::get<1>(ft).second;
-      const Index& q_index = std::get<2>(ft).second;
+      const Index& p_index = ft.point_s_.second;
+      const Index& q_index = ft.point_t_.second;
 
       Vertex_handle vp,vq;
       if ( ! domain_.is_loop(curve_index) )
diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h
index 8eed7475b808..72bc54181d15 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h
@@ -46,7 +46,7 @@
 #include <optional>
 #include <boost/optional/optional_io.hpp>
 #include <boost/mpl/has_xxx.hpp>
-#include <CGAL/tuple.h>
+
 #include <sstream>
 #include <atomic>
 
@@ -54,30 +54,40 @@ namespace CGAL {
 
 namespace Mesh_3 {
 
+  template<typename Facet>
+  struct Facets_erase_counters
+  {
+    const Facet f1_;
+    const unsigned int f1_erase_counter_;
+    const Facet f2_;
+    const unsigned int f2_erase_counter_;
+  };
+
   // Predicate to know if a facet in a refinement queue is a zombie
   // A facet is a pair <cell, index of the opposite vertex>.
   // A facet is a "zombie" if at least one of its two adjacent cells
   // has been erased. We test it thanks to the "erase counter" which
   // is inside each cell (incremented by the compact container).
-  // In the refinement queue, we store a tuple
+  // In the refinement queue, we store a struct
   // <facet1, facet1_erase counter, facet2, facet2_erase_counter>
   // where facet2 = mirror_facet(facet1) and facetx_erase_counter is
   // the erase_counter of facetx's cell when added to the queue>
   template<typename Facet>
   class Facet_to_refine_is_not_zombie
   {
+    typedef Facets_erase_counters<Facet> Facets_ec;
+
   public:
     Facet_to_refine_is_not_zombie() {}
 
-    bool operator()(const std::tuple<
-      Facet, unsigned int, Facet, unsigned int> &f) const
+    bool operator()(const Facets_ec& f) const
     {
 #ifdef _DEBUG
       /*
-      int f1_current_erase_counter = std::get<0>(f).first->erase_counter();
-      int f1_saved_erase_counter = std::get<1>(f);
-      int f2_current_erase_counter = std::get<2>(f).first->erase_counter();
-      int f2_saved_erase_counter = std::get<3>(f);
+      int f1_current_erase_counter = f.f1_.first->erase_counter();
+      int f1_saved_erase_counter = f.f1_erase_counter_;
+      int f2_current_erase_counter = f.f2_.first->erase_counter();
+      int f2_saved_erase_counter = f.f2_erase_counter_;
       */
       //f1_current_erase_counter - f1_saved_erase_counter + f2_current_erase_counter - f2_saved_erase_counter == 1
 
@@ -89,7 +99,7 @@ namespace Mesh_3 {
 #endif
 
         std::stringstream sstr;
-        Facet facet = std::get<0>(f);
+        Facet facet = f.f1_;
         sstr << "Facet 1 { " << std::endl
         << "  - " << *facet.first->vertex((facet.second+1)%4)  << std::endl
         << "  - " << *facet.first->vertex((facet.second+2)%4)  << std::endl
@@ -97,7 +107,7 @@ namespace Mesh_3 {
         << "  - 4th vertex in cell: " << *facet.first->vertex(facet.second)  << std::endl
         << "}" << std::endl;
 
-        facet = std::get<2>(f);
+        facet = f.f2_;
         sstr << "Facet 2 { " << std::endl
         << "  - " << *facet.first->vertex((facet.second+1)%4)  << std::endl
         << "  - " << *facet.first->vertex((facet.second+2)%4)  << std::endl
@@ -109,8 +119,8 @@ namespace Mesh_3 {
         std::cerr << s << std::endl;
       }*/
 #endif
-      return (std::get<0>(f).first->erase_counter() == std::get<1>(f)
-        && std::get<2>(f).first->erase_counter() == std::get<3>(f) );
+      return (f.f1_.first->erase_counter() == f.f1_erase_counter_
+           && f.f2_.first->erase_counter() == f.f2_erase_counter_ );
     }
   };
 
@@ -123,6 +133,8 @@ namespace Mesh_3 {
 template <typename Index, typename Facet, typename Concurrency_tag>
 class Refine_facets_3_handle_queue_base
 {
+  typedef Facets_erase_counters<Facet> Facets_counters;
+
 protected:
   Refine_facets_3_handle_queue_base() : m_last_vertex_index() {}
 
@@ -139,21 +151,20 @@ class Refine_facets_3_handle_queue_base
 #if defined(CGAL_MESH_3_USE_LAZY_SORTED_REFINEMENT_QUEUE) \
  || defined(CGAL_MESH_3_USE_LAZY_UNSORTED_REFINEMENT_QUEUE)
 
-  std::tuple<Facet, unsigned int, Facet, unsigned int>
+  Facets_counters
   from_facet_to_refinement_queue_element(const Facet &facet,
                                          const Facet &mirror) const
   {
-    return std::make_tuple(
-      facet, facet.first->erase_counter(),
-      mirror, mirror.first->erase_counter());
+    return Facets_counters{facet, facet.first->erase_counter(),
+                           mirror, mirror.first->erase_counter()};
   }
 
 public:
   template<typename Container_element>
   Facet extract_element_from_container_value(const Container_element &e) const
   {
-    // We get the first Facet inside the tuple
-    return std::get<0>(e);
+    // We get the first Facet inside the struct
+    return e.f1_;
   }
 
 #else
@@ -198,21 +209,20 @@ class Refine_facets_3_handle_queue_base<Index, Facet, Parallel_tag>
     m_last_vertex_index.local() = i;
   }
 
-  std::tuple<Facet, unsigned int, Facet, unsigned int>
+  Facets_erase_counters<Facet>
   from_facet_to_refinement_queue_element(const Facet &facet,
                                          const Facet &mirror) const
   {
-    return std::make_tuple(
-      facet, facet.first->erase_counter(),
-      mirror, mirror.first->erase_counter());
+    return Facets_erase_counters<Facet>{facet, facet.first->erase_counter(),
+                                        mirror, mirror.first->erase_counter()};
   }
 
 public:
   template<typename Container_element>
   Facet extract_element_from_container_value(const Container_element &e) const
   {
-    // We get the first Facet inside the tuple
-    return std::get<0>(e);
+    // We get the first Facet inside the struct
+    return e.f1_;
   }
 
 protected:
@@ -247,6 +257,8 @@ class Refine_facets_3_base
   typedef typename GT::Ray_3 Ray_3;
   typedef typename GT::Line_3 Line_3;
 
+  typedef typename MeshDomain::Intersection Intersection;
+
 public:
   Refine_facets_3_base(Tr& tr, Complex3InTriangulation3& c3t3,
                        const MeshDomain& oracle,
@@ -356,6 +368,19 @@ class Refine_facets_3_base
     return sstr.str();
   }
 
+  int dimension(const Intersection& intersection) const
+  {
+    return std::get<2>(intersection);
+  }
+  const typename MeshDomain::Index& index(const Intersection& intersection) const
+  {
+    return std::get<1>(intersection);
+  }
+  const typename MeshDomain::Point_3& point(const Intersection& intersection) const
+  {
+    return std::get<0>(intersection);
+  }
+
 protected:
 
   // Functor for scan_triangulation_impl function
@@ -395,9 +420,13 @@ class Refine_facets_3_base
   typedef typename MeshDomain::Surface_patch_index Surface_patch_index;
   typedef typename MeshDomain::Index Index;
 
-  typedef typename std::optional<
-    std::tuple<Surface_patch_index, Index, Bare_point> >
-                                                      Facet_properties;
+  struct Facet_prop
+  {
+    Surface_patch_index surface_patch_index;
+    Index index;
+    Bare_point point;
+  };
+  typedef typename std::optional<Facet_prop> Facet_properties;
 
 
   /// Returns canonical facet of facet
@@ -648,8 +677,7 @@ template<class Tr,
           Meshes::Filtered_multimap_container
 # endif
           <
-            std::tuple<typename Tr::Facet, unsigned int,
-                               typename Tr::Facet, unsigned int>,
+            Facets_erase_counters<typename Tr::Facet>,
             typename Criteria::Facet_quality,
             Facet_to_refine_is_not_zombie<typename Tr::Facet>,
             Concurrency_tag
@@ -658,8 +686,7 @@ template<class Tr,
 # ifdef CGAL_MESH_3_USE_LAZY_UNSORTED_REFINEMENT_QUEUE
           Meshes::Filtered_deque_container
           <
-            std::tuple<typename Tr::Facet, unsigned int,
-                               typename Tr::Facet, unsigned int>,
+            Facets_erase_counters<typename Tr::Facet>,
             typename Criteria::Facet_quality,
             Facet_to_refine_is_not_zombie<typename Tr::Facet>,
             Concurrency_tag
@@ -667,8 +694,7 @@ template<class Tr,
 # elif defined(CGAL_MESH_3_USE_LAZY_SORTED_REFINEMENT_QUEUE)
           Meshes::Filtered_multimap_container
           <
-            std::tuple<typename Tr::Facet, unsigned int,
-                               typename Tr::Facet, unsigned int>,
+            Facets_erase_counters<typename Tr::Facet>,
             typename Criteria::Facet_quality,
             Facet_to_refine_is_not_zombie<typename Tr::Facet>,
             Concurrency_tag
@@ -686,8 +712,7 @@ template<class Tr,
 # ifdef CGAL_MESH_3_USE_LAZY_UNSORTED_REFINEMENT_QUEUE
           Meshes::Filtered_deque_container
           <
-            std::tuple<typename Tr::Facet, unsigned int,
-                               typename Tr::Facet, unsigned int>,
+            Facets_erase_counters<typename Tr::Facet>,
             typename Criteria::Facet_quality,
             Facet_to_refine_is_not_zombie<typename Tr::Facet>,
             Concurrency_tag
@@ -695,8 +720,7 @@ template<class Tr,
 # elif defined(CGAL_MESH_3_USE_LAZY_SORTED_REFINEMENT_QUEUE)
           Meshes::Filtered_multimap_container
           <
-            std::tuple<typename Tr::Facet, unsigned int,
-                               typename Tr::Facet, unsigned int>,
+            Facets_erase_counters<typename Tr::Facet>,
             typename Criteria::Facet_quality,
             Facet_to_refine_is_not_zombie<typename Tr::Facet>,
             Concurrency_tag
@@ -1351,9 +1375,7 @@ conflicts_zone_impl(const Weighted_point& point
       this->compute_facet_properties(facet, properties, /*force_exact=*/true);
       if ( properties )
       {
-        const Surface_patch_index& surface_index = std::get<0>(*properties);
-        const Index& surface_center_index = std::get<1>(*properties);
-        const Bare_point& surface_center = std::get<2>(*properties);
+        const auto& [surface_index, surface_center_index, surface_center] = *properties;
 
         // Facet is on surface: set facet properties
         this->set_facet_surface_center(facet, surface_center, surface_center_index);
@@ -1413,9 +1435,7 @@ conflicts_zone_impl(const Weighted_point& point
       this->compute_facet_properties(facet, properties, /*force_exact=*/true);
       if ( properties )
       {
-        const Surface_patch_index& surface_index = std::get<0>(*properties);
-        const Index& surface_center_index = std::get<1>(*properties);
-        const Bare_point& surface_center = std::get<2>(*properties);
+        const auto& [surface_index, surface_center_index, surface_center] = *properties;
 
         // Facet is on surface: set facet properties
         this->set_facet_surface_center(facet, surface_center, surface_center_index);
@@ -1581,9 +1601,7 @@ treat_new_facet(Facet& facet)
   compute_facet_properties(facet, properties);
   if ( properties )
   {
-    const Surface_patch_index& surface_index = std::get<0>(*properties);
-    const Index& surface_center_index = std::get<1>(*properties);
-    const Bare_point& surface_center = std::get<2>(*properties);
+    const auto& [surface_index, surface_center_index, surface_center] = *properties;
 
     // Facet is on surface: set facet properties
     set_facet_surface_center(facet, surface_center, surface_center_index);
@@ -1681,16 +1699,12 @@ compute_facet_properties(const Facet& facet,
       // test "intersect == Intersection()" (aka empty intersection), but
       // the later does not work.
       Surface_patch surface =
-        (std::get<2>(intersect) == 0) ? Surface_patch() :
+        (dimension(intersect) == 0) ? Surface_patch() :
         Surface_patch(
-          r_oracle_.surface_patch_index(std::get<1>(intersect)));
+          r_oracle_.surface_patch_index(index(intersect)));
       if(surface)
 #endif // CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
-      {
-        fp =  Facet_properties(std::make_tuple(*surface,
-                                      std::get<1>(intersect),
-                                      std::get<0>(intersect)));
-      }
+      fp =  Facet_prop{*surface, index(intersect), point(intersect)};
     }
   }
   // If the dual is a ray
@@ -1721,15 +1735,13 @@ compute_facet_properties(const Facet& facet,
       Intersection intersect = construct_intersection(ray);
 #ifdef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
       Surface_patch surface =
-        (std::get<2>(intersect) == 0) ? Surface_patch() :
+        (dimension(intersect) == 0) ? Surface_patch() :
         Surface_patch(
-          r_oracle_.surface_patch_index(std::get<1>(intersect)));
+          r_oracle_.surface_patch_index(index(intersect)));
       if(surface)
 #endif // CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
       {
-        fp = Facet_properties(std::make_tuple(*surface,
-                                      std::get<1>(intersect),
-                                      std::get<0>(intersect)));
+        fp = Facet_prop{*surface, index(intersect), point(intersect)};
       }
     }
   }
diff --git a/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h b/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h
index 787b443489c4..7b3dde34edcd 100644
--- a/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h
+++ b/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h
@@ -20,7 +20,6 @@
 #include <utility> // std::swap
 #include <algorithm> // std::min
 
-#include <CGAL/tuple.h>
 #include <CGAL/Image_3.h>
 #include <CGAL/number_utils.h>
 #include <CGAL/squared_distance_3.h>
@@ -515,11 +514,7 @@ polylines_to_protect
       for(int j = 0; j < ydim; j+= (axis == 1 ? (std::max)(1, ydim-1) : 1 ) )
         for(int k = 0; k < zdim; k+= (axis == 2 ? (std::max)(1, zdim-1) : 1 ) )
         {
-
           using std::array;
-          using std::tuple;
-          using std::get;
-
           typedef array<int, 3> Pixel;
 
 #ifdef CGAL_MESH_3_DEBUG_POLYLINES_TO_PROTECT
diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h
index ddb1ba1109f6..d31ac841756f 100644
--- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h
+++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h
@@ -961,9 +961,9 @@ get_curves(OutputIterator out) const
       q_index = p_index;
     }
 
-    *out++ = std::make_tuple(eit->first,
-                                     std::make_pair(p,p_index),
-                                     std::make_pair(q,q_index));
+    *out++ = {eit->first,
+              std::make_pair(p,p_index),
+              std::make_pair(q,q_index)};
   }
 
   return out;