diff --git a/.github/workflows/build_doc.yml b/.github/workflows/build_doc.yml
index 3484e24bfa5..a2ee453953d 100644
--- a/.github/workflows/build_doc.yml
+++ b/.github/workflows/build_doc.yml
@@ -3,7 +3,6 @@ name: Documentation
 on:
   issue_comment:
       types: [created]
-  workflow_dispatch:
 
 permissions:
   contents: read # to fetch code (actions/checkout)
@@ -67,7 +66,6 @@ jobs:
         with:
            repository: ${{ github.repository }}
            ref: refs/pull/${{ steps.get_pr_number.outputs.result }}/merge
-           token: ${{ secrets.PUSH_TO_CGAL_GITHUB_IO_TOKEN }}
            fetch-depth: 2
 
       - name: install dependencies
diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/overlay_unbounded.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/overlay_unbounded.cpp
index be62a491733..866c66016d1 100644
--- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/overlay_unbounded.cpp
+++ b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/overlay_unbounded.cpp
@@ -2,7 +2,6 @@
 // A face overlay of two arrangements with unbounded faces.
 
 #include <string>
-#include <boost/lexical_cast.hpp>
 
 #include <CGAL/basic.h>
 #include <CGAL/Arr_extended_dcel.h>
@@ -14,7 +13,7 @@
 // Define a functor for creating a label from a character and an integer.
 struct Overlay_label {
   std::string operator()(char c, unsigned int i) const
-  { return c + boost::lexical_cast<std::string>(i); }
+  { return c + std::to_string(i); }
 };
 
 typedef CGAL::Arr_face_extended_dcel<Traits, char>         Dcel_dlue;
diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_topology_traits/Arr_spherical_topology_traits_2_impl.h b/Arrangement_on_surface_2/include/CGAL/Arr_topology_traits/Arr_spherical_topology_traits_2_impl.h
index fbe5e059eef..9175587a317 100644
--- a/Arrangement_on_surface_2/include/CGAL/Arr_topology_traits/Arr_spherical_topology_traits_2_impl.h
+++ b/Arrangement_on_surface_2/include/CGAL/Arr_topology_traits/Arr_spherical_topology_traits_2_impl.h
@@ -215,16 +215,16 @@ is_in_face(const Face* f, const Point_2& p, const Vertex* v) const
 
     /*! We identify 2 main cases:
      * 1. The vertical ray intersects the boundary at a halfedge. In this
-     * case the x-possition of p is strictly larger than the x-possition of
-     * the current-curve source, and strictly smaller than x-possition of
+     * case the x-position of p is strictly larger than the x-position of
+     * the current-curve source, and strictly smaller than x-position of
      * the current-curve target, or vice versa.
      * 2. The vertical ray intersects the boundary at a vertex. In this case:
-     * a. the x-possition of p is strictly smaller than the x-position of the
+     * a. the x-position of p is strictly smaller than the x-position of the
      * current-curve source, and equal to the x-position of the current-curve
      * target, and
-     * b. the x-possition of p is equal to the x-position of the next-curve
+     * b. the x-position of p is equal to the x-position of the next-curve
      * source (not counting vertical curves in between), and strictly larger
-     * than the x-possition of the next-curve target, or vice verase (that is,
+     * than the x-position of the next-curve target, or vice verase (that is,
      * the "smaller" and "larger" interchanged).
      */
 
diff --git a/BGL/include/CGAL/draw_face_graph.h b/BGL/include/CGAL/draw_face_graph.h
index 02a96433146..7e1de8031ec 100644
--- a/BGL/include/CGAL/draw_face_graph.h
+++ b/BGL/include/CGAL/draw_face_graph.h
@@ -19,21 +19,32 @@
 #include <CGAL/Random.h>
 #include <CGAL/boost/graph/helpers.h>
 
-namespace CGAL
-{
+namespace CGAL {
 
 // Default color functor; user can change it to have its own face color
 struct DefaultColorFunctorFaceGraph
 {
   template<typename Graph>
-  CGAL::IO::Color operator()(const Graph&,
-                         typename boost::graph_traits<Graph>::face_descriptor fh) const
+  CGAL::IO::Color operator()(const Graph& /*g*/,
+                             typename boost::graph_traits<Graph>::face_descriptor /*f*/) const
   {
-    if (fh==boost::graph_traits<Graph>::null_face()) // use to get the mono color
-      return CGAL::IO::Color(100, 125, 200); // R G B between 0-255
-
     return get_random_color(CGAL::get_default_random());
   }
+
+  // edges and vertices are black by default
+  template<typename Graph>
+  CGAL::IO::Color operator()(const Graph& /*g*/,
+                             typename boost::graph_traits<Graph>::edge_descriptor /*e*/) const
+  {
+    return IO::black();
+  }
+
+  template<typename Graph>
+  CGAL::IO::Color operator()(const Graph& /*g*/,
+                             typename boost::graph_traits<Graph>::vertex_descriptor /*v*/) const
+  {
+    return IO::black();
+  }
 };
 
 class SimpleFaceGraphViewerQt : public Basic_viewer_qt
@@ -48,28 +59,28 @@ public:
   }
 
   /// Construct the viewer.
-  /// @param amesh the surface mesh to view
+  /// @param g the face graph to view
   /// @param title the title of the window
   /// @param anofaces if true, do not draw faces (faces are not computed; this can be
-  ///        useful for very big object where this time could be long)
-  template <typename SM>
+  ///        useful for very big objects where this time could be long)
+  template <typename Graph>
   SimpleFaceGraphViewerQt(QWidget* parent,
-                          const SM& amesh,
-                          const char* title="Basic Surface_mesh Viewer",
+                          const Graph& g,
+                          const char* title="Basic Face Graph Viewer",
                           bool anofaces=false) :
-    SimpleFaceGraphViewerQt(parent, amesh, title, anofaces, DefaultColorFunctorFaceGraph())
+    SimpleFaceGraphViewerQt(parent, g, title, anofaces, DefaultColorFunctorFaceGraph())
   {
   }
 
-  template <typename SM, typename ColorFunctor>
+  template <typename Graph, typename ColorFunctor>
   SimpleFaceGraphViewerQt(QWidget* parent,
-                          const SM& amesh,
+                          const Graph& g,
                           const char* title,
                           bool anofaces,
                           ColorFunctor fcolor) :
     // First draw: no vertex; edges, faces; mono-color; inverse normal
     Base(parent, title, false, true, true, true, false),
-    m_compute_elements_impl(compute_elements_functor(amesh, anofaces, fcolor))
+    m_compute_elements_impl(compute_elements_functor(g, anofaces, fcolor))
   {
   }
 
@@ -82,43 +93,42 @@ public:
     m_compute_elements_impl();
   }
 
-  template <typename SM, typename ColorFunctor>
-  void set_face_graph(const SM& amesh,
+  template <typename Graph, typename ColorFunctor>
+  void set_face_graph(const Graph& g,
                       bool anofaces,
                       ColorFunctor fcolor) {
-    m_compute_elements_impl = compute_elements_functor(amesh, anofaces, fcolor);
+    m_compute_elements_impl = compute_elements_functor(g, anofaces, fcolor);
   }
 
-  template <typename SM, typename ColorFunctor>
-  void set_face_graph(const SM& amesh,
+  template <typename Graph, typename ColorFunctor>
+  void set_face_graph(const Graph& g,
                       bool anofaces=false) {
-    set_mesh(amesh, anofaces, DefaultColorFunctorFaceGraph());
+    set_mesh(g, anofaces, DefaultColorFunctorFaceGraph());
   }
 protected:
-  template <typename SM, typename ColorFunctor>
+  template <typename Graph, typename ColorFunctor>
   std::function<void()>
-  compute_elements_functor(const SM& sm,
+  compute_elements_functor(const Graph& g,
                            bool anofaces,
                            ColorFunctor fcolor)
   {
-    using Point =
-        typename boost::property_map_value<SM, CGAL::vertex_point_t>::type;
+    using Point = typename boost::property_map_value<Graph, CGAL::vertex_point_t>::type;
     using Kernel = typename CGAL::Kernel_traits<Point>::Kernel;
     using Vector = typename Kernel::Vector_3;
 
-    auto vnormals = get(CGAL::dynamic_vertex_property_t<Vector>(), sm);
-    auto point_pmap = get(CGAL::vertex_point, sm);
-    for (auto v : vertices(sm))
+    auto vnormals = get(CGAL::dynamic_vertex_property_t<Vector>(), g);
+    auto point_pmap = get(CGAL::vertex_point, g);
+    for (auto v : vertices(g))
     {
       Vector n(NULL_VECTOR);
       int i=0;
-      for (auto h : halfedges_around_target(halfedge(v, sm), sm))
+      for (auto h : halfedges_around_target(halfedge(v, g), g))
       {
-        if (!is_border(h, sm))
+        if (!is_border(h, g))
         {
           Vector ni = CGAL::cross_product(
-                        Vector(get(point_pmap, source(h, sm)), get(point_pmap, target(h, sm))),
-                        Vector(get(point_pmap, target(h, sm)), get(point_pmap, target(next(h, sm), sm))));
+                        Vector(get(point_pmap, source(h, g)), get(point_pmap, target(h, g))),
+                        Vector(get(point_pmap, target(h, g)), get(point_pmap, target(next(h, g), g))));
           if (ni != NULL_VECTOR)
           {
             n+=ni;
@@ -131,41 +141,41 @@ protected:
 
     // This function return a lambda expression, type-erased in a
     // `std::function<void()>` object.
-    return [this, &sm, vnormals, anofaces, fcolor, point_pmap]()
+    return [this, &g, vnormals, anofaces, fcolor, point_pmap]()
     {
       this->clear();
 
       if (!anofaces)
       {
-        for (auto fh: faces(sm))
+        for (auto fh: faces(g))
         {
-          if (fh!=boost::graph_traits<SM>::null_face())
-          {
-            CGAL::IO::Color c=fcolor(sm, fh);
-            face_begin(c);
-            auto hd=halfedge(fh, sm);
-            const auto first_hd = hd;
-            do
-              {
-                auto v = source(hd, sm);
-                add_point_in_face(get(point_pmap, v), get(vnormals, v));
-                hd=next(hd, sm);
-              }
-            while(hd!=first_hd);
-            face_end();
-          }
+          const CGAL::IO::Color& c = fcolor(g, fh);
+          face_begin(c);
+          auto hd=halfedge(fh, g);
+          const auto first_hd = hd;
+          do
+            {
+              auto v = source(hd, g);
+              add_point_in_face(get(point_pmap, v), get(vnormals, v));
+              hd=next(hd, g);
+            }
+          while(hd!=first_hd);
+          face_end();
         }
       }
 
-      for (auto e: edges(sm))
+      for (auto e: edges(g))
       {
-        add_segment(get(point_pmap, source(halfedge(e, sm), sm)),
-                    get(point_pmap, target(halfedge(e, sm), sm)));
+        const CGAL::IO::Color& c = fcolor(g, e);
+        add_segment(get(point_pmap, source(halfedge(e, g), g)),
+                    get(point_pmap, target(halfedge(e, g), g)),
+                    c);
       }
 
-      for (auto v: vertices(sm))
+      for (auto v: vertices(g))
       {
-        this->add_point(get(point_pmap, v));
+        const CGAL::IO::Color& c = fcolor(g, v);
+        this->add_point(get(point_pmap, v), c);
       }
     };
   }
diff --git a/Cone_spanners_2/examples/Cone_spanners_2/theta_io.cpp b/Cone_spanners_2/examples/Cone_spanners_2/theta_io.cpp
index 4096cf1c71f..3b49478f47f 100644
--- a/Cone_spanners_2/examples/Cone_spanners_2/theta_io.cpp
+++ b/Cone_spanners_2/examples/Cone_spanners_2/theta_io.cpp
@@ -4,7 +4,6 @@
 #include <iterator>
 #include <string>
 #include <vector>
-#include <boost/lexical_cast.hpp>
 #include <CGAL/Exact_predicates_exact_constructions_kernel_with_root_of.h>
 #include <CGAL/Construct_theta_graph_2.h>
 #include <CGAL/gnuplot_output_2.h>
@@ -79,7 +78,7 @@ int main(int argc, char ** argv)
   // obtain the number of vertices in the constructed graph
   boost::graph_traits<Graph>::vertices_size_type n = boost::num_vertices(g);
   // generate gnuplot files for plotting this graph
-  std::string file_prefix = "t" + boost::lexical_cast<std::string>(k) + "n" + boost::lexical_cast<std::string>(n);
+  std::string file_prefix = "t" + std::to_string(k) + "n" + std::to_string(n);
   CGAL::gnuplot_output_2(g, file_prefix);
 
   return 0;
diff --git a/Cone_spanners_2/test/Cone_spanners_2/theta_exact.cpp b/Cone_spanners_2/test/Cone_spanners_2/theta_exact.cpp
index 0173c9798ac..6bdb2e0eb8d 100644
--- a/Cone_spanners_2/test/Cone_spanners_2/theta_exact.cpp
+++ b/Cone_spanners_2/test/Cone_spanners_2/theta_exact.cpp
@@ -10,7 +10,6 @@
 #include <iterator>
 #include <string>
 #include <vector>
-#include <boost/lexical_cast.hpp>
 #include <CGAL/Exact_predicates_exact_constructions_kernel_with_root_of.h>
 #include <CGAL/Construct_theta_graph_2.h>
 #include <CGAL/gnuplot_output_2.h>
@@ -70,7 +69,7 @@ int main(int argc, char ** argv)
   // obtain the number of vertices in the constructed graph
   boost::graph_traits<Graph>::vertices_size_type n = boost::num_vertices(g);
   // generate gnuplot files for plotting this graph
-  std::string file_prefix = "t" + boost::lexical_cast<std::string>(k) + "n" + boost::lexical_cast<std::string>(n);
+  std::string file_prefix = "t" + std::to_string(k) + "n" + std::to_string(n);
   CGAL::gnuplot_output_2(g, file_prefix);
 
   return 0;
diff --git a/Cone_spanners_2/test/Cone_spanners_2/theta_inexact.cpp b/Cone_spanners_2/test/Cone_spanners_2/theta_inexact.cpp
index 4d305d0fb30..f62d9d34dc9 100644
--- a/Cone_spanners_2/test/Cone_spanners_2/theta_inexact.cpp
+++ b/Cone_spanners_2/test/Cone_spanners_2/theta_inexact.cpp
@@ -11,7 +11,6 @@
 #include <iterator>
 #include <string>
 #include <vector>
-#include <boost/lexical_cast.hpp>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Construct_theta_graph_2.h>
 #include <CGAL/gnuplot_output_2.h>
@@ -71,7 +70,7 @@ int main(int argc, char ** argv)
   // obtain the number of vertices in the constructed graph
   boost::graph_traits<Graph>::vertices_size_type n = boost::num_vertices(g);
   // generate gnuplot files for plotting this graph
-  std::string file_prefix = "t" + boost::lexical_cast<std::string>(k) + "n" + boost::lexical_cast<std::string>(n);
+  std::string file_prefix = "t" + std::to_string(k) + "n" + std::to_string(n);
   CGAL::gnuplot_output_2(g, file_prefix);
 
   return 0;
diff --git a/Cone_spanners_2/test/Cone_spanners_2/yao_exact.cpp b/Cone_spanners_2/test/Cone_spanners_2/yao_exact.cpp
index 9b6c3668777..4e19687de49 100644
--- a/Cone_spanners_2/test/Cone_spanners_2/yao_exact.cpp
+++ b/Cone_spanners_2/test/Cone_spanners_2/yao_exact.cpp
@@ -12,7 +12,6 @@
 #include <iterator>
 #include <string>
 #include <vector>
-#include <boost/lexical_cast.hpp>
 #include <CGAL/Exact_predicates_exact_constructions_kernel_with_root_of.h>
 #include <CGAL/Construct_yao_graph_2.h>
 #include <CGAL/gnuplot_output_2.h>
@@ -72,7 +71,7 @@ int main(int argc, char ** argv)
   boost::graph_traits<Graph>::vertices_size_type n = boost::num_vertices(g);
 
   // generate gnuplot files for plotting this graph
-  std::string fileprefix = "y" + boost::lexical_cast<std::string>(k) + "n" + boost::lexical_cast<std::string>(n);
+  std::string fileprefix = "y" + std::to_string(k) + "n" + std::to_string(n);
   CGAL::gnuplot_output_2(g, fileprefix);
 
   return 0;
diff --git a/Cone_spanners_2/test/Cone_spanners_2/yao_inexact.cpp b/Cone_spanners_2/test/Cone_spanners_2/yao_inexact.cpp
index 8f6221935ed..5ddf1a2c3f8 100644
--- a/Cone_spanners_2/test/Cone_spanners_2/yao_inexact.cpp
+++ b/Cone_spanners_2/test/Cone_spanners_2/yao_inexact.cpp
@@ -12,7 +12,6 @@
 #include <iterator>
 #include <string>
 #include <vector>
-#include <boost/lexical_cast.hpp>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Construct_yao_graph_2.h>
 #include <CGAL/gnuplot_output_2.h>
@@ -72,7 +71,7 @@ int main(int argc, char ** argv)
   boost::graph_traits<Graph>::vertices_size_type n = boost::num_vertices(g);
 
   // generate gnuplot files for plotting this graph
-  std::string fileprefix = "y" + boost::lexical_cast<std::string>(k) + "n" + boost::lexical_cast<std::string>(n);
+  std::string fileprefix = "y" + std::to_string(k) + "n" + std::to_string(n);
   CGAL::gnuplot_output_2(g, fileprefix);
 
   return 0;
diff --git a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/External_structure_builder.h b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/External_structure_builder.h
index 823a9b7c4ce..8c4c546590d 100644
--- a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/External_structure_builder.h
+++ b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/External_structure_builder.h
@@ -17,6 +17,7 @@
 
 
 #include <CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Nef_3/SNC_external_structure.h>
 #include <CGAL/Nef_3/SNC_intersection.h>
 
 #undef CGAL_NEF_DEBUG
diff --git a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Ray_hit_generator.h b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Ray_hit_generator.h
index ffc0868e952..2768f43ce72 100644
--- a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Ray_hit_generator.h
+++ b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Ray_hit_generator.h
@@ -17,6 +17,7 @@
 
 #include <CGAL/Modifier_base.h>
 #include <CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Nef_3/SNC_point_locator.h>
 #include <CGAL/Nef_3/SNC_intersection.h>
 #include <CGAL/Convex_decomposition_3/SM_walls.h>
 
diff --git a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/SFace_separator.h b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/SFace_separator.h
index 2c42e64fe72..cd9ed4851d0 100644
--- a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/SFace_separator.h
+++ b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/SFace_separator.h
@@ -14,8 +14,9 @@
 
 #include <CGAL/license/Convex_decomposition_3.h>
 
-
+#include <CGAL/Modifier_base.h>
 #include <CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Nef_3/SNC_point_locator.h>
 #include <CGAL/Convex_decomposition_3/SM_walls.h>
 
 namespace CGAL {
diff --git a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/SM_walls.h b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/SM_walls.h
index 27a1c4f697a..e1c5b099414 100644
--- a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/SM_walls.h
+++ b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/SM_walls.h
@@ -14,6 +14,8 @@
 
 #include <CGAL/license/Convex_decomposition_3.h>
 
+#include <CGAL/Nef_S2/SM_decorator.h>
+#include <CGAL/Nef_S2/SM_point_locator.h>
 
 #undef CGAL_NEF_DEBUG
 #define CGAL_NEF_DEBUG 227
diff --git a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Single_wall_creator2.h b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Single_wall_creator2.h
index bd4d22bac38..954138100e5 100644
--- a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Single_wall_creator2.h
+++ b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Single_wall_creator2.h
@@ -14,10 +14,12 @@
 
 #include <CGAL/license/Convex_decomposition_3.h>
 
-
+#include <CGAL/Modifier_base.h>
 #include <CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Nef_3/SNC_point_locator.h>
 #include <CGAL/Nef_3/SNC_intersection.h>
 #include <CGAL/Nef_S2/Normalizing.h>
+#include <CGAL/Convex_decomposition_3/Ray_hit_generator.h>
 #include <CGAL/Convex_decomposition_3/SM_walls.h>
 
 #undef CGAL_NEF_DEBUG
diff --git a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Single_wall_creator3.h b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Single_wall_creator3.h
index fc26c54ceb0..2277f58d5cf 100644
--- a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Single_wall_creator3.h
+++ b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/Single_wall_creator3.h
@@ -14,10 +14,12 @@
 
 #include <CGAL/license/Convex_decomposition_3.h>
 
-
+#include <CGAL/Modifier_base.h>
 #include <CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Nef_3/SNC_point_locator.h>
 #include <CGAL/Nef_3/SNC_intersection.h>
 #include <CGAL/Nef_S2/Normalizing.h>
+#include <CGAL/Convex_decomposition_3/Ray_hit_generator.h>
 #include <CGAL/Convex_decomposition_3/SM_walls.h>
 
 #undef CGAL_NEF_DEBUG
diff --git a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/YVertical_wall_builder.h b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/YVertical_wall_builder.h
index 0db21512f4c..e43ef7d92d3 100644
--- a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/YVertical_wall_builder.h
+++ b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/YVertical_wall_builder.h
@@ -14,8 +14,9 @@
 
 #include <CGAL/license/Convex_decomposition_3.h>
 
-
-#include<CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Modifier_base.h>
+#include <CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Convex_decomposition_3/is_reflex_sedge.h>
 #include <CGAL/Convex_decomposition_3/Single_wall_creator3.h>
 #include <CGAL/Convex_decomposition_3/External_structure_builder.h>
 
diff --git a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/is_reflex_sedge.h b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/is_reflex_sedge.h
index 5a881e30331..9409c04846c 100644
--- a/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/is_reflex_sedge.h
+++ b/Convex_decomposition_3/include/CGAL/Convex_decomposition_3/is_reflex_sedge.h
@@ -14,6 +14,8 @@
 
 #include <CGAL/license/Convex_decomposition_3.h>
 
+#include <CGAL/enum.h>
+#include <CGAL/Origin.h>
 
 #undef CGAL_NEF_DEBUG
 #define CGAL_NEF_DEBUG 239
diff --git a/Documentation/doc/biblio/geom.bib b/Documentation/doc/biblio/geom.bib
index e9a313dde70..930ef665ca5 100644
--- a/Documentation/doc/biblio/geom.bib
+++ b/Documentation/doc/biblio/geom.bib
@@ -137313,7 +137313,7 @@ Contains C code."
 
 @inproceedings{ss-kaud-88
 , author =	"Th. Strothotte and J.-R. Sack"
-, title =	"Knowledge Aquisition using Diagrams"
+, title =	"Knowledge Acquisition using Diagrams"
 , booktitle =	"Proc. 3rd IFIP Conference on Man-Machine Systems"
 , site =	"Oulo, Finland"
 , year =	1988
diff --git a/GraphicsView/demo/Largest_empty_rect_2/Largest_empty_rectangle_2.cpp b/GraphicsView/demo/Largest_empty_rect_2/Largest_empty_rectangle_2.cpp
index 012ef71d3e5..67e4f336075 100644
--- a/GraphicsView/demo/Largest_empty_rect_2/Largest_empty_rectangle_2.cpp
+++ b/GraphicsView/demo/Largest_empty_rect_2/Largest_empty_rectangle_2.cpp
@@ -98,6 +98,8 @@ public Q_SLOTS:
 
   void on_actionClear_triggered();
 
+  void on_actionOpen_triggered();
+
   void processInput(CGAL::Object);
 
   void on_actionRecenter_triggered();
@@ -105,7 +107,7 @@ public Q_SLOTS:
   void on_actionGeneratePointsInSquare_triggered();
   void on_actionGeneratePointsInDisc_triggered();
   void clear();
-
+  void open(QString fileName);
   void update_largest_empty_rectangle();
 
 Q_SIGNALS:
@@ -229,6 +231,50 @@ MainWindow::on_actionClear_triggered()
   Q_EMIT( changed());
 }
 
+void
+MainWindow::on_actionOpen_triggered()
+{
+  QString fileName = QFileDialog::getOpenFileName(this,
+                                                  tr("Open points file"),
+                                                  "."
+                                                  ,tr("xy files (*.xy)")
+                                                  );
+  if(! fileName.isEmpty()){
+    open(fileName);
+  }
+
+}
+
+void
+MainWindow::open(QString fileName)
+{
+  // wait cursor
+  QApplication::setOverrideCursor(Qt::WaitCursor);
+  std::ifstream ifs(qPrintable(fileName));
+
+  clear();
+
+  Point_2 p;
+  while(ifs >> p){
+    points.push_back(p);
+  }
+
+  CGAL::Bbox_2 bbox = CGAL::bbox_2(points.begin(), points.end());
+  square = Iso_rectangle_2(bbox);
+
+  ler = Largest_empty_iso_rectangle_2(square);
+  ler.insert(points.begin(), points.end());
+
+  frame[0]->setLine(convert(Segment_2(square.vertex(0),square.vertex(1))));
+  frame[1]->setLine(convert(Segment_2(square.vertex(1), square.vertex(2))));
+  frame[2]->setLine(convert(Segment_2(square.vertex(2), square.vertex(3))));
+  frame[3]->setLine(convert(Segment_2(square.vertex(3), square.vertex(0))));
+
+  QApplication::restoreOverrideCursor();
+  on_actionRecenter_triggered();
+  Q_EMIT( changed());
+}
+
 void
 MainWindow::on_actionRecenter_triggered()
 {
diff --git a/GraphicsView/demo/Largest_empty_rect_2/Largest_empty_rectangle_2.ui b/GraphicsView/demo/Largest_empty_rect_2/Largest_empty_rectangle_2.ui
index 19c46f07374..a6ddf170ca2 100644
--- a/GraphicsView/demo/Largest_empty_rect_2/Largest_empty_rectangle_2.ui
+++ b/GraphicsView/demo/Largest_empty_rect_2/Largest_empty_rectangle_2.ui
@@ -80,7 +80,7 @@
      <x>0</x>
      <y>0</y>
      <width>500</width>
-     <height>26</height>
+     <height>22</height>
     </rect>
    </property>
    <widget class="QMenu" name="menuFile">
@@ -90,6 +90,7 @@
     <addaction name="separator"/>
     <addaction name="actionClear"/>
     <addaction name="separator"/>
+    <addaction name="actionOpen"/>
     <addaction name="actionQuit"/>
    </widget>
    <widget class="QMenu" name="menuTools">
@@ -199,6 +200,11 @@
     <string>Generate Segment Fans</string>
    </property>
   </action>
+  <action name="actionOpen">
+   <property name="text">
+    <string>Open</string>
+   </property>
+  </action>
  </widget>
  <resources>
   <include location="Largest_empty_rectangle_2.qrc"/>
diff --git a/Inscribed_areas/doc/Inscribed_areas/CGAL/Largest_empty_iso_rectangle_2.h b/Inscribed_areas/doc/Inscribed_areas/CGAL/Largest_empty_iso_rectangle_2.h
index f11af4ccdf1..080fc203bfb 100644
--- a/Inscribed_areas/doc/Inscribed_areas/CGAL/Largest_empty_iso_rectangle_2.h
+++ b/Inscribed_areas/doc/Inscribed_areas/CGAL/Largest_empty_iso_rectangle_2.h
@@ -142,9 +142,11 @@ Iso_rectangle_2 get_bounding_box();
 /// @{
 
 /*!
-Inserts point `p` in the point set, if it is not already in the set.
+Inserts point `p` in the point set, if it is not already in the set
+and on the bounded side of the bounding rectangle.
+\note Points on the boundary can be ignored as they lead to the same result.
 */
-void
+bool
 insert(const Point_2& p);
 
 /*!
diff --git a/Inscribed_areas/examples/Inscribed_areas/largest_empty_rectangle.cpp b/Inscribed_areas/examples/Inscribed_areas/largest_empty_rectangle.cpp
index 5128e958d84..bbf8989b7b8 100644
--- a/Inscribed_areas/examples/Inscribed_areas/largest_empty_rectangle.cpp
+++ b/Inscribed_areas/examples/Inscribed_areas/largest_empty_rectangle.cpp
@@ -1,8 +1,7 @@
 #include <CGAL/Simple_cartesian.h>
-#include <CGAL/Iso_rectangle_2.h>
 #include <CGAL/Largest_empty_iso_rectangle_2.h>
 
-#include <fstream>
+#include <iostream>
 
 typedef double                                 Number_Type;
 typedef CGAL::Simple_cartesian<Number_Type>    K;
diff --git a/Inscribed_areas/include/CGAL/Largest_empty_iso_rectangle_2.h b/Inscribed_areas/include/CGAL/Largest_empty_iso_rectangle_2.h
index 95391cd85a5..4ef236892ab 100644
--- a/Inscribed_areas/include/CGAL/Largest_empty_iso_rectangle_2.h
+++ b/Inscribed_areas/include/CGAL/Largest_empty_iso_rectangle_2.h
@@ -762,7 +762,7 @@ bool
 Largest_empty_iso_rectangle_2<T>::insert(const Point_2& _p)
 {
   // check that the point is inside the bounding box
-  if(bbox_p.has_on_unbounded_side(_p)) {
+  if(! bbox_p.has_on_bounded_side(_p)) {
     return(false);
   }
 
diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md
index 30a49be1ae5..c5106c6f430 100644
--- a/Installation/CHANGES.md
+++ b/Installation/CHANGES.md
@@ -20,8 +20,8 @@ Release date: June 2023
   `CGAL::Polygon_mesh_processing::triangulate_and_refine_hole()`, and `CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole()`
   which have output iterators for vertices and faces as parameter. They are replaced by overloads with two additional named parameters.
 
--   Added the function `CGAL::Polygon_mesh_processing::surface_Delaunay_remeshing()`, that remeshes a surface triangle mesh following the
-CGAL tetrahedral Delaunay refinement algorithm.
+-   Added the function `CGAL::Polygon_mesh_processing::surface_Delaunay_remeshing()`, that remeshes a surface triangle mesh using
+  the Delaunay refinement algorithm from the 3D Mesh Generation package.
 
 -   Added the function `CGAL::Polygon_mesh_processing::remove_almost_degenerate_faces()` to remove badly shaped triangles faces in a mesh.
 
diff --git a/Intersections_2/include/CGAL/Intersections_2/Segment_2_Segment_2.h b/Intersections_2/include/CGAL/Intersections_2/Segment_2_Segment_2.h
index 57f2e31b5bb..28d7920bd2a 100644
--- a/Intersections_2/include/CGAL/Intersections_2/Segment_2_Segment_2.h
+++ b/Intersections_2/include/CGAL/Intersections_2/Segment_2_Segment_2.h
@@ -435,6 +435,22 @@ Segment_2_Segment_2_pair<K>::intersection_type() const
                                            : CGAL::make_array( _seg2->point(s2s2_id[c][2]), _seg2->point(s2s2_id[c][3]),
                                                                _seg1->point(s2s2_id[c][0]), _seg1->point(s2s2_id[c][1]) );
 
+    // special case for vertical and horizontal segments
+    if (std::is_floating_point<typename K::FT>::value &&
+        std::is_same<typename K::Kernel_tag, Cartesian_tag>::value)
+    {
+      if (pts[0].x()==pts[1].x() && pts[2].y()==pts[3].y())
+      {
+        _intersection_point = K().construct_point_2_object()(pts[0].x(), pts[2].y());
+        return _result;
+      }
+      if (pts[0].y()==pts[1].y() && pts[2].x()==pts[3].x())
+      {
+        _intersection_point = K().construct_point_2_object()(pts[2].x(), pts[0].y());
+        return _result;
+      }
+    }
+
     typename K::FT alpha =  s2s2_alpha(pts[0].x(), pts[0].y(), pts[1].x(), pts[1].y(), pts[2].x(), pts[2].y(), pts[3].x(), pts[3].y());
 
     _intersection_point = K().construct_barycenter_2_object()(pts[0], alpha, pts[1]);
diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Triangle_3_intersection.h
index 7f76c9d6000..d16b33d6570 100644
--- a/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Triangle_3_intersection.h
+++ b/Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Triangle_3_intersection.h
@@ -18,6 +18,7 @@
 
 #include <CGAL/enum.h>
 #include <CGAL/kernel_assertions.h>
+#include <iterator>
 
 namespace CGAL {
 namespace Intersections {
@@ -141,7 +142,7 @@ intersection(const typename K::Plane_3& plane,
   CGAL_kernel_assertion(pts.size() == 2);
 
   return intersection_return<typename K::Intersect_3, typename K::Plane_3, typename K::Triangle_3>(
-           k.construct_segment_3_object()(*pts.begin(), *boost::prior(pts.end())));
+           k.construct_segment_3_object()(*pts.begin(), *std::prev(pts.end())));
 }
 
 template <class K>
diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h
index fa148adcba2..424e1fc9361 100644
--- a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h
+++ b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h
@@ -51,7 +51,7 @@ void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p,
   for (Iterator it=inter_pts.begin();it!=inter_pts.end();++it)
     orientations[ &(*it) ]=orient(p,q,r,*it);
 
-  int pt_added = 0;
+  CGAL_kernel_assertion_code(int pt_added = 0;)
 
   const typename Kernel::Point_3* prev = &(*boost::prior(inter_pts.end()));
   Iterator stop = inter_pts.size() > 2 ? inter_pts.end() : boost::prior(inter_pts.end());
@@ -75,7 +75,7 @@ void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p,
 
       prev = &(*inter_pts.insert(it,*inter));
       orientations[prev] = COLLINEAR;
-      ++pt_added;
+      CGAL_kernel_assertion_code(++pt_added;)
     }
 
     prev = &(*it);
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Aff_transformation_2.h b/Kernel_23/doc/Kernel_23/CGAL/Aff_transformation_2.h
index 0300279b5ee..4bd99fa3771 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Aff_transformation_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Aff_transformation_2.h
@@ -107,7 +107,7 @@ approximates the rotation over the angle indicated by direction
 `d`, such that the differences between the sines and cosines
 of the rotation given by d and the approximating rotation
 are at most \f$ num/den\f$ each.
-\pre \f$ num/den>0\f$ and \f$ d != 0\f$.
+\pre `num/den > 0` and  `d != 0`.
 */
 Aff_transformation_2(const Rotation,
 const Direction_2<Kernel> &d,
@@ -116,7 +116,7 @@ const Kernel::RT &den = RT(1));
 
 /*!
 introduces a rotation by the angle `rho`.
-\pre \f$ sine\_rho^2 + cosine\_rho^2 == hw^2\f$.
+\pre <tt>sine\_rho<sup>2</sup> + cosine\_rho<sup>2</sup> == hw<sup>2</sup></tt>.
 */
 Aff_transformation_2(const Rotation,
 const Kernel::RT &sine_rho,
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Bbox_2.h b/Kernel_23/doc/Kernel_23/CGAL/Bbox_2.h
index 61267300f39..ec861b9cb72 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Bbox_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Bbox_2.h
@@ -77,13 +77,13 @@ double ymax() const;
 
 /*!
 Returns `xmin()` if `i==0` or `ymin()` if `i==1`.
-\pre i==0 or i==1
+\pre `i==0` or `i==1`
 */
 double min(int i) const;
 
 /*!
 Returns `xmax()` if `i==0` or `ymax()` if `i==1`.
-\pre i==0 or i==1
+\pre `i==0` or `i==1`
 */
 double max(int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h b/Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h
index ad69f1d00aa..a61d5e339af 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Bbox_3.h
@@ -90,14 +90,14 @@ double zmax() const;
 /*!
 Returns `xmin()` if `i==0` or `ymin()` if `i==1`
 or `zmin()` if `i==2`.
-\pre i>=0 and i<=2
+\pre `i>=0` and `i<=2`
 */
 double min(int i) const;
 
 /*!
 Returns `xmax()` if `i==0` or `ymax()` if `i==1`
 or `zmax()` if `i==2`.
-\pre i>=0 and i<=2
+\pre `i>=0` and `i<=2`
 */
 double max(int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Circle_2.h b/Kernel_23/doc/Kernel_23/CGAL/Circle_2.h
index 01238a5f044..1e9eacdeb0e 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Circle_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Circle_2.h
@@ -28,7 +28,7 @@ introduces a variable `c` of type `Circle_2`.
 It is initialized to the circle with center `center`,
 squared radius `squared_radius` and orientation
 `ori`.
-\pre `ori` \f$ \neq\f$ `COLLINEAR`, and further, `squared_radius` \f$ \geq\f$ 0.
+\pre `ori != COLLINEAR` and `squared_radius >= 0`.
 */
 Circle_2(const Point_2<Kernel> &center,
 const Kernel::FT &squared_radius,
@@ -52,7 +52,7 @@ const Point_2<Kernel> &r);
 introduces a variable `c` of type `Circle_2`.
 It is initialized to the circle with diameter \f$ \overline{pq}\f$
 and orientation `ori`.
-\pre `ori` \f$ \neq\f$ `COLLINEAR`.
+\pre `ori != COLLINEAR`.
 */
 Circle_2( const Point_2<Kernel> &p,
 const Point_2<Kernel> &q,
@@ -63,7 +63,7 @@ const Orientation &ori = COUNTERCLOCKWISE);
 introduces a variable `c` of type `Circle_2`.
 It is initialized to the circle with center `center`, squared
 radius zero and orientation `ori`.
-\pre `ori` \f$ \neq\f$ `COLLINEAR`.
+\pre `ori != COLLINEAR`.
 \post `c.is_degenerate()` = `true`.
 */
 Circle_2( const Point_2<Kernel> &center,
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Circle_3.h b/Kernel_23/doc/Kernel_23/CGAL/Circle_3.h
index 96063004383..1fcdcab6e68 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Circle_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Circle_3.h
@@ -21,7 +21,7 @@ public:
 introduces a variable `c` of type `Circle_3`.
 It is initialized to the circle of center `center` and
 squared radius `sq_r` in plane `plane`.
-\pre `center` lies in `plane` and                 `sq_r` \f$ \geq\f$ 0.
+\pre `center` lies in `plane` and `sq_r >= 0`.
 */
 Circle_3(const Point_3<Kernel> &center,
          const Kernel::FT &sq_r,
@@ -32,7 +32,7 @@ introduces a variable `c` of type `Circle_3`.
 It is initialized to the circle of center `center` and
 squared radius `sq_r` in a plane normal to
 the vector `n`.
-\pre `sq_r` \f$ \geq\f$ 0.
+\pre `sq_r >= 0`.
 */
 Circle_3(const Point_3<Kernel> & center,
          const Kernel::FT & sq_r,
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Direction_2.h b/Kernel_23/doc/Kernel_23/CGAL/Direction_2.h
index 5b2843b994e..00a2e8d74a5 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Direction_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Direction_2.h
@@ -60,7 +60,7 @@ Direction_2(const Kernel::RT &x, const Kernel::RT &y);
 
 /*!
 returns values, such that `d``== Direction_2<Kernel>(delta(0),delta(1))`.
-\pre \f$ 0 \leq i \leq1\f$.
+\pre `0 <= i <= 1`.
 */
 Kernel::RT delta(int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Direction_3.h b/Kernel_23/doc/Kernel_23/CGAL/Direction_3.h
index f11bf023cb3..81b096ac359 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Direction_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Direction_3.h
@@ -57,7 +57,7 @@ Direction_3(const Kernel::RT &x, const Kernel::RT &y, const Kernel::RT &z);
 
 /*!
 returns values, such that `d``== Direction_3<Kernel>(delta(0),delta(1),delta(2))`.
-\pre \f$ 0 \leq i \leq2\f$.
+\pre `0 <= i <= 2`.
 */
 Kernel::RT delta(int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Iso_cuboid_3.h b/Kernel_23/doc/Kernel_23/CGAL/Iso_cuboid_3.h
index 8cbbd395085..a32da1d2490 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Iso_cuboid_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Iso_cuboid_3.h
@@ -39,7 +39,7 @@ const Point_3<Kernel> &q);
 introduces an iso-oriented cuboid `c` with diagonal
 opposite vertices `p` and `q`. The `int` argument value
 is only used to distinguish the two overloaded functions.
-\pre `p.x()<=q.x()`, `p.y()<=q.y()`and `p.z()<=q.z()`.
+\pre `p.x()<=q.x()`, `p.y()<=q.y()` and `p.z()<=q.z()`.
 */
 Iso_cuboid_3(const Point_3<Kernel> &p,
 const Point_3<Kernel> &q, int);
@@ -65,7 +65,7 @@ introduces an iso-oriented cuboid `c` with diagonal
 opposite vertices
 (`min_hx/hw`, `min_hy/hw`, `min_hz/hw`) and
 (`max_hx/hw`, `max_hy/hw`, `max_hz/hw`).
-\pre `hw` \f$ \neq\f$ 0.
+\pre `hw != 0`.
 */
 Iso_cuboid_3(
 const Kernel::RT& min_hx, const Kernel::RT& min_hy, const Kernel::RT& min_hz,
@@ -156,14 +156,14 @@ Kernel::FT zmax() const;
 /*!
 returns `i`-th %Cartesian coordinate of
 the smallest vertex of `c`.
-\pre \f$ 0 \leq i \leq2\f$.
+\pre `0 <= i <= 2`.
 */
 Kernel::FT min_coord(int i) const;
 
 /*!
 returns `i`-th %Cartesian coordinate of
 the largest vertex of `c`.
-\pre \f$ 0 \leq i \leq2\f$.
+\pre `0 <= i <= 2`.
 */
 Kernel::FT max_coord(int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Iso_rectangle_2.h b/Kernel_23/doc/Kernel_23/CGAL/Iso_rectangle_2.h
index 175d1261bea..d7228803e0e 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Iso_rectangle_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Iso_rectangle_2.h
@@ -62,7 +62,7 @@ const Point_2<Kernel> &top);
 introduces an iso-oriented rectangle `r` with diagonal
 opposite vertices (`min_hx/hw`, `min_hy/hw`) and
 (`max_hx/hw`, `max_hy/hw`).
-\pre `hw` \f$ \neq\f$ 0.
+\pre `hw != 0`.
 */
 Iso_rectangle_2(const Kernel::RT& min_hx, const Kernel::RT& min_hy,
 const Kernel::RT& max_hx, const Kernel::RT& max_hy,
@@ -134,14 +134,14 @@ Kernel::FT ymax() const;
 /*!
 returns the `i`'th %Cartesian coordinate of the
 lower left vertex of `r`.
-\pre \f$ 0 \leq i \leq1\f$.
+\pre `0 <= i <= 1`.
 */
 Kernel::FT min_coord(int i) const;
 
 /*!
 returns the `i`'th %Cartesian coordinate of the
 upper right vertex of `r`.
-\pre \f$ 0 \leq i \leq1\f$.
+\pre `0 <= i <= 1`.
 */
 Kernel::FT max_coord(int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Kernel/global_functions.h b/Kernel_23/doc/Kernel_23/CGAL/Kernel/global_functions.h
index 3e4bba13a43..45c56fe7710 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Kernel/global_functions.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Kernel/global_functions.h
@@ -89,7 +89,7 @@ Angle angle(const CGAL::Point_3<Kernel>&p,
 /*!
 returns an approximation of the angle between `p-q` and `r-q`.
 The angle is given in degrees.
-\pre `p` and `r` are not equal to `q`.
+\pre `p != q` and `r != q`.
 */
 template <typename Kernel>
 Kernel::FT approximate_angle(const CGAL::Point_3<Kernel>& p,
@@ -341,7 +341,7 @@ const CGAL::Point_3<Kernel>& p4, const Kernel::FT&w4);
 /*!
 constructs the bisector line of the two points `p` and `q`.
 The bisector is oriented in such a way that `p` lies on its
-positive side. \pre `p` and `q` are not equal.
+positive side. \pre `p != q`.
 */
 template <typename Kernel>
 CGAL::Line_2<Kernel> bisector(const CGAL::Point_2<Kernel> &p,
@@ -367,7 +367,7 @@ const CGAL::Line_2<Kernel> &l2);
 /*!
 constructs the bisector plane of the two points `p` and `q`.
 The bisector is oriented in such a way that `p` lies on its
-positive side. \pre `p` and `q` are not equal.
+positive side. \pre `p != q'.
 */
 template <typename Kernel>
 CGAL::Plane_3<Kernel> bisector(const CGAL::Point_3<Kernel> &p,
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Point_2.h b/Kernel_23/doc/Kernel_23/CGAL/Point_2.h
index 8ac3eff82d6..00ea07ef2fb 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Point_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Point_2.h
@@ -71,7 +71,7 @@ Point_2(double x, double y);
 
 /*!
 introduces a point `p` initialized to `(hx/hw,hy/hw)`.
-\pre `hw` \f$ \neq\f$ `Kernel::RT(0)`.
+\pre `hw != Kernel::RT(0)`.
 */
 Point_2(const Kernel::RT &hx, const Kernel::RT &hy, const Kernel::RT &hw = RT(1));
 
@@ -159,19 +159,19 @@ Kernel::FT y() const;
 
 /*!
 returns the i'th homogeneous coordinate of `p`.
-\pre \f$ 0\leq i \leq2\f$.
+\pre `0 <= i <= 2`.
 */
 Kernel::RT homogeneous(int i) const;
 
 /*!
 returns the i'th %Cartesian coordinate of `p`.
-\pre \f$ 0\leq i \leq1\f$.
+\pre `0 <= i <= 1`.
 */
 Kernel::FT cartesian(int i) const;
 
 /*!
 returns `cartesian(i)`.
-\pre \f$ 0\leq i \leq1\f$.
+\pre `0 <= i <= 1`.
 */
 Kernel::FT operator[](int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Point_3.h b/Kernel_23/doc/Kernel_23/CGAL/Point_3.h
index deed3e522d7..0babd70a4c8 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Point_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Point_3.h
@@ -56,7 +56,7 @@ Point_3(double x, double y, double z);
 
 /*!
 introduces a point `p` initialized to `(hx/hw,hy/hw, hz/hw)`.
-\pre `hw` \f$ \neq\f$ 0.
+\pre `hw != 0`.
 */
 Point_3(const Kernel::RT &hx, const Kernel::RT &hy, const Kernel::RT &hz, const Kernel::RT &hw = RT(1));
 
@@ -154,19 +154,19 @@ Kernel::FT z() const;
 
 /*!
 returns the i'th homogeneous coordinate of `p`.
-\pre \f$ 0\leq i \leq3\f$.
+\pre `0 <= i <= 3`.
 */
 Kernel::RT homogeneous(int i) const;
 
 /*!
 returns the i'th %Cartesian coordinate of `p`.
-\pre \f$ 0\leq i \leq2\f$.
+\pre `0 <= i <= 2`.
 */
 Kernel::FT cartesian(int i) const;
 
 /*!
 returns `cartesian(i)`.
-\pre \f$ 0\leq i \leq2\f$.
+\pre `0 <= i <= 2`.
 */
 Kernel::FT operator[](int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Ray_2.h b/Kernel_23/doc/Kernel_23/CGAL/Ray_2.h
index 87cd676cde9..c153b0f9ee3 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Ray_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Ray_2.h
@@ -65,7 +65,7 @@ Point_2<Kernel> source() const;
 /*!
 returns a point on `r`. `point(0)` is the source,
 `point(i)`, with `i>0`, is different from the
-source. \pre \f$ i \geq0\f$.
+source. \pre `i >= 0`.
 */
 Point_2<Kernel> point(const Kernel::FT i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Ray_3.h b/Kernel_23/doc/Kernel_23/CGAL/Ray_3.h
index dd6926c3565..86e6efce248 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Ray_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Ray_3.h
@@ -65,7 +65,7 @@ Point_3<Kernel> source() const;
 /*!
 returns a point on `r`. `point(0)` is the source.
 `point(i)`, with `i>0`, is different from the
-source. \pre \f$ i \geq0\f$.
+source. \pre `i >= 0`.
 */
 Point_3<Kernel> point(const Kernel::FT i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Sphere_3.h b/Kernel_23/doc/Kernel_23/CGAL/Sphere_3.h
index 053bac81119..f8541de6dd3 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Sphere_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Sphere_3.h
@@ -28,7 +28,7 @@ introduces a variable `c` of type `Sphere_3`.
 It is initialized to the sphere with center `center`,
 squared radius `squared_radius` and orientation
 `orientation`.
-\pre `orientation` \f$ \neq\f$ \ref COPLANAR, and furthermore, `squared_radius` \f$ \geq\f$ 0.
+\pre `orientation != COPLANAR` and `squared_radius >= 0`.
 */
 Sphere_3( const Point_3<Kernel> & center,
           const Kernel::FT & squared_radius,
@@ -53,7 +53,7 @@ const Point_3<Kernel> & s);
 introduces a variable `c` of type `Sphere_3`.
 It is initialized to the smallest sphere which passes through
 the points `p`, `q`, and `r`. The orientation of
-the sphere is `o`. \pre `o` is not \ref COPLANAR.
+the sphere is `o`. \pre `o != COPLANAR`.
 */
 Sphere_3( const Point_3<Kernel> & p,
 const Point_3<Kernel> & q,
@@ -65,7 +65,7 @@ const Orientation& o = COUNTERCLOCKWISE);
 introduces a variable `c` of type `Sphere_3`.
 It is initialized to the smallest sphere which passes through
 the points `p` and `q`. The orientation of
-the sphere is `o`. \pre `o` is not \ref COPLANAR.
+the sphere is `o`. \pre `o != COPLANAR`.
 */
 Sphere_3( const Point_3<Kernel> & p,
 const Point_3<Kernel> & q,
@@ -76,7 +76,7 @@ const Orientation& o = COUNTERCLOCKWISE);
 introduces a variable `c` of type `Sphere_3`.
 It is initialized to the sphere with center `center`, squared
 radius zero and orientation `orientation`.
-\pre `orientation` \f$ \neq\f$ \ref COPLANAR.
+\pre `orientation != COPLANAR`.
 \post `c.is_degenerate()` = `true`.
 */
 Sphere_3( const Point_3<Kernel> & center,
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Vector_2.h b/Kernel_23/doc/Kernel_23/CGAL/Vector_2.h
index 42176a6daed..eaeab03a882 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Vector_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Vector_2.h
@@ -71,7 +71,7 @@ Vector_2(double x, double y);
 
 /*!
 introduces a vector `v` initialized to `(hx/hw,hy/hw)`.
-\pre \f$ hw\neq0\f$.
+\pre `hw != 0`.
 */
 Vector_2(const Kernel::RT &hx, const Kernel::RT &hy, const Kernel::RT &hw = RT(1));
 
@@ -126,19 +126,19 @@ Kernel::FT y() const;
 
 /*!
 returns the i'th homogeneous coordinate of `v`.
-\pre \f$ 0\leq i \leq2\f$.
+\pre `0 <= i <= 2`.
 */
 Kernel::RT homogeneous(int i) const;
 
 /*!
 returns the i'th Cartesian coordinate of `v`.
-\pre \f$ 0\leq i \leq1\f$.
+\pre `0 <= i <= 1`.
 */
 Kernel::FT cartesian(int i) const;
 
 /*!
 returns `cartesian(i)`.
-\pre \f$ 0\leq i \leq1\f$.
+\pre `0 <= i <= 1`.
 */
 Kernel::FT operator[](int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Vector_3.h b/Kernel_23/doc/Kernel_23/CGAL/Vector_3.h
index 036dafdc8f2..d126bd97740 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Vector_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Vector_3.h
@@ -137,19 +137,19 @@ Kernel::FT z() const;
 
 /*!
 returns the i'th homogeneous coordinate of `v`.
-\pre \f$ 0\leq i \leq3\f$.
+\pre `0 <= i <= 3`.
 */
 Kernel::RT homogeneous(int i) const;
 
 /*!
 returns the i'th %Cartesian coordinate of `v`.
-\pre \f$ 0\leq i \leq2\f$.
+\pre  `0 <= i <= 2`
 */
 Kernel::FT cartesian(int i) const;
 
 /*!
 returns `cartesian(i)`.
-\pre \f$ 0\leq i \leq2\f$.
+\pre  `0 <= i <= 2`
 */
 Kernel::FT operator[](int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Weighted_point_2.h b/Kernel_23/doc/Kernel_23/CGAL/Weighted_point_2.h
index 4c2d0a15cbc..01a147e6d64 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Weighted_point_2.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Weighted_point_2.h
@@ -147,19 +147,19 @@ public:
 
   /*!
   returns the i'th homogeneous coordinate of `p`.
-  \pre \f$ 0\leq i \leq2\f$.
+  \pre  `0 <= i <= 2`
   */
   Kernel::RT homogeneous(int i) const;
 
   /*!
   returns the i'th %Cartesian coordinate of `p`.
-  \pre \f$ 0\leq i \leq1\f$.
+  \pre  `0 <= i <= 1`
   */
   Kernel::FT cartesian(int i) const;
 
   /*!
   returns `cartesian(i)`.
-  \pre \f$ 0\leq i \leq1\f$.
+  \pre  `0 <= i <= 1`
   */
   Kernel::FT operator[](int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/Weighted_point_3.h b/Kernel_23/doc/Kernel_23/CGAL/Weighted_point_3.h
index e652a47049d..3dcb39e1786 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/Weighted_point_3.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/Weighted_point_3.h
@@ -157,19 +157,19 @@ public:
 
   /*!
   returns the i'th homogeneous coordinate of `p`.
-  \pre \f$ 0\leq i \leq3\f$.
+  \pre  `0 <= i <= 3`
   */
   Kernel::RT homogeneous(int i) const;
 
   /*!
   returns the i'th %Cartesian coordinate of `p`.
-  \pre \f$ 0\leq i \leq2\f$.
+  \pre  `0 <= i <= 2`
   */
   Kernel::FT cartesian(int i) const;
 
   /*!
   returns `cartesian(i)`.
-  \pre \f$ 0\leq i \leq2\f$.
+  \pre  `0 <= i <= 2`
   */
   Kernel::FT operator[](int i) const;
 
diff --git a/Kernel_23/doc/Kernel_23/CGAL/rational_rotation.h b/Kernel_23/doc/Kernel_23/CGAL/rational_rotation.h
index 38dab692034..c56d44465fe 100644
--- a/Kernel_23/doc/Kernel_23/CGAL/rational_rotation.h
+++ b/Kernel_23/doc/Kernel_23/CGAL/rational_rotation.h
@@ -7,7 +7,7 @@ computes integers `sin_num`, `cos_num` and `denom`, such
 that `sin_num`/`denom` approximates the sine of direction
 \f$ (\f$`dirx`,`diry`\f$ )\f$. The difference between the sine and
 the approximating rational is bounded by `eps_num`/`eps_den`.
-\pre `eps_num` \f$ \neq0\f$.
+\pre `eps_num != 0`.
 
 \cgalHeading{Implementation}
 
diff --git a/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h b/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
index c770d0b7b40..ac8e47c1e54 100644
--- a/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
+++ b/Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h
@@ -3879,7 +3879,7 @@ public:
   /*!
     constructs the bisector of `p` and `q`.
     The bisector is oriented in such a way that `p` lies on its
-    positive side. \pre `p` and `q` are not equal.
+    positive side. \pre `p != q`.
   */
   Kernel::Line_2 operator()(const Kernel::Point_2&p,
                             const Kernel::Point_2&q );
@@ -3920,7 +3920,7 @@ public:
   /*!
     constructs the bisector plane of `p` and `q`.
     The bisector is oriented in such a way that `p` lies on its
-    positive side. \pre `p` and `q` are not equal.
+    positive side. \pre `p != q`.
   */
   Kernel::Plane_3 operator()(const Kernel::Point_3&p,
                              const Kernel::Point_3&q );
@@ -4200,7 +4200,7 @@ public:
     It is initialized to the circle with center `center`,
     squared radius `squared_radius` and orientation
     `orientation`.
-    \pre `orientation` \f$ \neq\f$ \ref CGAL::COLLINEAR, and further, `squared_radius` \f$ \geq\f$ 0.
+    \pre `orientation != CGAL::COLLINEAR` and `squared_radius >= 0`.
   */
   Kernel::Circle_2 operator()( Kernel::Point_2 const& center,
                                Kernel::FT const& squared_radius,
@@ -4225,7 +4225,7 @@ public:
     introduces a variable of type `Kernel::Circle_2`.
     It is initialized to the circle with diameter `pq`
     and orientation `orientation`.
-    \pre `orientation` \f$ \neq\f$ \ref CGAL::COLLINEAR.
+    \pre `orientation != CGAL::COLLINEAR`.
   */
   Kernel::Circle_2 operator()( Kernel::Point_2 const& p,
                                Kernel::Point_2 const& q,
@@ -4237,7 +4237,7 @@ public:
     introduces a variable of type `Kernel::Circle_2`.
     It is initialized to the circle with center `center`, squared
     radius zero and orientation `orientation`.
-    \pre `orientation` \f$ \neq\f$ \ref CGAL::COLLINEAR.
+    \pre `orientation != CGAL::COLLINEAR`.
     \post .`is_degenerate()` = `true`.
   */
   Kernel::Circle_2 operator()( Kernel::Point_2 const& center,
@@ -4269,7 +4269,7 @@ public:
     introduces a variable of type `Kernel::Circle_3`.
     It is initialized to the circle with center `center`,
     and squared radius `sq_r` in the plane `plane`.
-    \pre `center` lies in `plane` and                 `sq_r` \f$ \geq\f$ 0.
+    \pre `center` lies in `plane` and  `sq_r >= 0`.
   */
   Kernel::Circle_3 operator()
   ( Kernel::Point_3 const& center,
@@ -4281,7 +4281,7 @@ public:
     It is initialized to the circle with center `center`,
     and squared radius `sq_r` in the plane
     containing `center` and normal to `n`.
-    \pre `sq_r` \f$ \geq\f$ 0.
+    \pre `sq_r >= 0`.
   */
   Kernel::Circle_3 operator()
   ( Kernel::Point_3 const& center,
@@ -5637,7 +5637,7 @@ public:
     introduces a direction orthogonal to `d`. If `o` is
     \ref CGAL::CLOCKWISE, `d` is rotated clockwise; if `o` is
     \ref CGAL::COUNTERCLOCKWISE, `d` is rotated counterclockwise.
-    \pre `o` is not \ref CGAL::COLLINEAR.
+    \pre `o != CGAL::COLLINEAR.`
   */
   Kernel::Direction_2 operator()(const Kernel::Direction_2& d,
                                  Orientation o);
@@ -5753,8 +5753,7 @@ public:
   /*!
     returns `v` rotated clockwise by 90 degrees, if `o` is
     \ref CGAL::CLOCKWISE, and rotated counterclockwise otherwise.
-    \pre `o` is not \ref CGAL::COLLINEAR.
-
+    \pre `o != CGAL::COLLINEAR`.
   */
   Kernel::Vector_2 operator()(const Kernel::Vector_2& v,
                               Orientation o);
@@ -6580,7 +6579,7 @@ public:
     introduces a sphere initialized to the sphere with center `center`,
     squared radius `squared_radius` and orientation
     `orientation`.
-    \pre `orientation` \f$ \neq\f$ \ref CGAL::COPLANAR, and furthermore, `squared_radius` \f$ \geq\f$ 0.
+    \pre `orientation != CGAL::COPLANAR` and `squared_radius >= 0`.
   */
   Kernel::Sphere_3 operator()(const Kernel::Point_3 & center,
                               const Kernel::FT & squared_radius,
@@ -6601,7 +6600,7 @@ public:
   /*!
     introduces a sphere initialized to the smallest sphere which passes
     through the points `p`, `q`, and `r`. The orientation of
-    the sphere is `o`. \pre `o` is not \ref CGAL::COPLANAR.
+    the sphere is `o`. \pre `o != CGAL::COPLANAR`.
   */
   Kernel::Sphere_3 operator()(const Kernel::Point_3 & p,
                               const Kernel::Point_3 & q,
@@ -6611,7 +6610,7 @@ public:
   /*!
     introduces a sphere initialized to the smallest sphere which passes
     through the points `p` and `q`. The orientation of
-    the sphere is `o`. \pre `o` is not \ref CGAL::COPLANAR.
+    the sphere is `o`. \pre `o != CGAL::COPLANAR`.
   */
   Kernel::Sphere_3 operator()(const Kernel::Point_3 & p,
                               const Kernel::Point_3 & q,
@@ -6620,7 +6619,7 @@ public:
   /*!
     introduces a sphere `s` initialized to the sphere with center
     `center`, squared radius zero and orientation `orientation`.
-    \pre `orientation` \f$ \neq\f$ \ref CGAL::COPLANAR.
+    \pre `orientation != CGAL::COPLANAR`.
     \post `s.is_degenerate()` = `true`.
   */
   Kernel::Sphere_3 operator()( const Kernel::Point_3 & center,
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 714430116f6..7b83d8b05bc 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h
@@ -974,7 +974,7 @@ update_mesh(const Moves_vector& moves,
       {
         FT size = std::get<2>(*it);
 
-#ifdef CGAL_MESH_3_OPTIMIZER_VERBOSE
+#ifdef CGAL_MESH_3_OPTIMIZER_VERY_VERBOSE
         std::cerr << "Moving #" << it - moves.begin()
                   << " addr: " << &*v
                   << " pt: " << tr_.point(v)
diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h b/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h
index 4f5115fb13c..185efe6229c 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h
@@ -575,7 +575,7 @@ refine_mesh(std::string dump_after_refine_surface_prefix)
   nbsteps = 0;
 
   facets_visitor_.activate();
-  dump_c3t3(r_c3t3_, dump_after_refine_surface_prefix);
+
   std::cerr << "Start volume scan...";
   CGAL_MESH_3_TASK_BEGIN(scan_cells_task_handle);
   cells_mesher_.scan_triangulation();
@@ -584,6 +584,7 @@ refine_mesh(std::string dump_after_refine_surface_prefix)
   std::cerr << "end scan. [Bad tets:" << cells_mesher_.size() << "]";
   std::cerr << std::endl << std::endl;
   elapsed_time += timer.time();
+  dump_c3t3(r_c3t3_, dump_after_refine_surface_prefix);
   timer.stop(); timer.reset(); timer.start();
 
   std::cerr << "Refining...\n";
diff --git a/Mesh_3/include/CGAL/Mesh_3/config.h b/Mesh_3/include/CGAL/Mesh_3/config.h
index 6fc06788902..cc2e5187f9b 100644
--- a/Mesh_3/include/CGAL/Mesh_3/config.h
+++ b/Mesh_3/include/CGAL/Mesh_3/config.h
@@ -57,4 +57,10 @@
 #  endif
 #endif
 
+#ifdef CGAL_MESH_3_VERY_VERBOSE
+#  ifndef CGAL_MESH_3_OPTIMIZER_VERY_VERBOSE
+#    define CGAL_MESH_3_OPTIMIZER_VERY_VERBOSE 1
+#  endif
+#endif
+
 #endif // CGAL_MESH_3_CONFIG_H
diff --git a/Nef_3/include/CGAL/Nef_3/Binary_operation.h b/Nef_3/include/CGAL/Nef_3/Binary_operation.h
index 558e3e09f1c..fa26c870f77 100644
--- a/Nef_3/include/CGAL/Nef_3/Binary_operation.h
+++ b/Nef_3/include/CGAL/Nef_3/Binary_operation.h
@@ -23,6 +23,7 @@
 #include <CGAL/Nef_S2/Normalizing.h>
 #include <CGAL/Unique_hash_map.h>
 #include <CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Nef_3/SNC_const_decorator.h>
 #include <CGAL/Nef_S2/SM_decorator.h>
 #include <CGAL/Nef_S2/SM_point_locator.h>
 #include <CGAL/Nef_3/SNC_SM_overlayer.h>
@@ -74,6 +75,7 @@ class Binary_operation : public CGAL::SNC_decorator<Map> {
   typedef typename SNC_structure::Items                Items;
   typedef typename Map::Sphere_map                     Sphere_map;
   typedef CGAL::SNC_decorator<SNC_structure>           SNC_decorator;
+  typedef CGAL::SNC_const_decorator<SNC_structure>     SNC_const_decorator;
   typedef SNC_decorator                                Base;
   typedef CGAL::SNC_constructor<Items, SNC_structure>  SNC_constructor;
   typedef CGAL::SNC_external_structure<Items, SNC_structure>
@@ -82,10 +84,13 @@ class Binary_operation : public CGAL::SNC_decorator<Map> {
   typedef CGAL::SNC_SM_overlayer<Items, SM_decorator>  SM_overlayer;
   typedef CGAL::SM_point_locator<SM_decorator>         SM_point_locator;
   typedef CGAL::SNC_point_locator<SNC_decorator>       SNC_point_locator;
+  typedef CGAL::SNC_point_locator<SNC_const_decorator> SNC_const_point_locator;
 
   typedef typename SNC_structure::Vertex_handle Vertex_handle;
   typedef typename SNC_structure::Halfedge_handle Halfedge_handle;
+  typedef typename SNC_structure::Halfedge_const_handle Halfedge_const_handle;
   typedef typename SNC_structure::Halffacet_handle Halffacet_handle;
+  typedef typename SNC_structure::Halffacet_const_handle Halffacet_const_handle;
   typedef typename SNC_structure::Volume_handle Volume_handle;
   typedef typename SNC_structure::SVertex_handle SVertex_handle;
   typedef typename SNC_structure::SHalfedge_handle SHalfedge_handle;
@@ -144,12 +149,12 @@ class Binary_operation : public CGAL::SNC_decorator<Map> {
     return v01;
   }
 
-  Vertex_handle create_local_view_on( const Point_3& p, Halfedge_handle e) {
+  Vertex_handle create_local_view_on( const Point_3& p, Halfedge_const_handle e) {
     SNC_constructor C(*this->sncp());
     return C.create_from_edge( e, p);
   }
 
-  Vertex_handle create_local_view_on( const Point_3& p, Halffacet_handle f) {
+  Vertex_handle create_local_view_on( const Point_3& p, Halffacet_const_handle f) {
     SNC_constructor C(*this->sncp());
     return C.create_from_facet( f, p);
   }
@@ -167,14 +172,14 @@ class Binary_operation : public CGAL::SNC_decorator<Map> {
             typename Selection,
             typename Association>
   class Intersection_call_back :
-    public SNC_point_locator::Intersection_call_back
+    public SNC_const_point_locator::Intersection_call_back
   {
     typedef typename SNC_decorator::Decorator_traits Decorator_traits;
     typedef typename Decorator_traits::Halfedge_handle Halfedge_handle;
     typedef typename Decorator_traits::Halffacet_handle Halffacet_handle;
 
   public:
-    Intersection_call_back( SNC_structure& s0, SNC_structure& s1,
+    Intersection_call_back( const SNC_structure& s0, const SNC_structure& s1,
                             const Selection& _bop, SNC_structure& r,
                             bool invert_order, Association& Ain) :
       snc0(s0), snc1(s1), bop(_bop), result(r),
@@ -456,12 +461,10 @@ class Binary_operation : public CGAL::SNC_decorator<Map> {
 
     // CGAL_NEF_SETDTHREAD(19*509*43*131);
 
-    Intersection_call_back<SNC_decorator, Selection, Association> call_back0
-      ( const_cast<SNC_structure&>(snc1), const_cast<SNC_structure&>(snc2),
-        BOP, *this->sncp(), false, A);
-    Intersection_call_back<SNC_decorator, Selection, Association> call_back1
-      ( const_cast<SNC_structure&>(snc2), const_cast<SNC_structure&>(snc2),
-        BOP, *this->sncp(), true, A);
+    Intersection_call_back<SNC_const_decorator, Selection, Association> call_back0
+      ( snc1, snc2, BOP, *this->sncp(), false, A);
+    Intersection_call_back<SNC_const_decorator, Selection, Association> call_back1
+      ( snc2, snc2, BOP, *this->sncp(), true, A);
 
 #ifdef CGAL_NEF3_TIMER_INTERSECTION
     double split_intersection = timer_overlay.time();
@@ -503,10 +506,8 @@ class Binary_operation : public CGAL::SNC_decorator<Map> {
                     << this->sncp()->number_of_vertices());
 #else
     CGAL_NEF_TRACEN("intersection by fast box intersection");
-        binop_intersection_test_segment_tree<SNC_decorator> binop_box_intersection;
-        binop_box_intersection(call_back0, call_back1,
-                               const_cast<SNC_structure&>(snc1),
-                               const_cast<SNC_structure&>(snc2));
+        binop_intersection_test_segment_tree<SNC_const_decorator> binop_box_intersection;
+        binop_box_intersection(call_back0, call_back1, snc1, snc2);
 #endif
 
 #ifdef CGAL_NEF3_TIMER_INTERSECTION
diff --git a/Nef_3/include/CGAL/Nef_3/Halfedge.h b/Nef_3/include/CGAL/Nef_3/Halfedge.h
index b0c21962585..b4fc7065740 100644
--- a/Nef_3/include/CGAL/Nef_3/Halfedge.h
+++ b/Nef_3/include/CGAL/Nef_3/Halfedge.h
@@ -21,6 +21,7 @@
 
 #include <string>
 #include <sstream>
+#include <CGAL/IO/io.h>
 #include <CGAL/IO/Verbose_ostream.h>
 #include <CGAL/Nef_3/SNC_iteration.h>
 #include <CGAL/Origin.h>
diff --git a/Nef_3/include/CGAL/Nef_3/Infimaximal_box.h b/Nef_3/include/CGAL/Nef_3/Infimaximal_box.h
index 05aabb0edd2..b1eb8d6b9ee 100644
--- a/Nef_3/include/CGAL/Nef_3/Infimaximal_box.h
+++ b/Nef_3/include/CGAL/Nef_3/Infimaximal_box.h
@@ -24,6 +24,7 @@
 #include <CGAL/Extended_homogeneous.h>
 
 #include <CGAL/Nef_3/SNC_intersection.h>
+#include <CGAL/Nef_3/SNC_iteration.h>
 
 namespace CGAL {
 
diff --git a/Nef_3/include/CGAL/Nef_3/K3_tree.h b/Nef_3/include/CGAL/Nef_3/K3_tree.h
index e147621b163..81eb55e4c4a 100644
--- a/Nef_3/include/CGAL/Nef_3/K3_tree.h
+++ b/Nef_3/include/CGAL/Nef_3/K3_tree.h
@@ -20,6 +20,7 @@
 #include <CGAL/basic.h>
 #include <CGAL/Unique_hash_map.h>
 #include <CGAL/Nef_3/quotient_coordinates_to_homogeneous_point.h>
+#include <CGAL/Nef_3/SNC_iteration.h>
 #include <CGAL/Lazy_kernel.h>
 #include <CGAL/Cartesian.h>
 #include <boost/container/deque.hpp>
diff --git a/Nef_3/include/CGAL/Nef_3/Mark_bounded_volumes.h b/Nef_3/include/CGAL/Nef_3/Mark_bounded_volumes.h
index f5af77cf1d5..0ce560af241 100644
--- a/Nef_3/include/CGAL/Nef_3/Mark_bounded_volumes.h
+++ b/Nef_3/include/CGAL/Nef_3/Mark_bounded_volumes.h
@@ -15,8 +15,8 @@
 
 #include <CGAL/license/Nef_3.h>
 
-
-#include <CGAL/Nef_polyhedron_3.h>
+#include <CGAL/Modifier_base.h>
+#include <CGAL/Nef_3/SNC_iteration.h>
 
 namespace CGAL {
 
diff --git a/Nef_3/include/CGAL/Nef_3/Nef_box.h b/Nef_3/include/CGAL/Nef_3/Nef_box.h
index 03a3a883ba8..c05d1e6c9ea 100644
--- a/Nef_3/include/CGAL/Nef_3/Nef_box.h
+++ b/Nef_3/include/CGAL/Nef_3/Nef_box.h
@@ -16,7 +16,11 @@
 
 #include <CGAL/box_intersection_d.h>
 #include <CGAL/Box_intersection_d/box_limits.h>
+#include <CGAL/Extended_cartesian.h>
+#include <CGAL/Extended_homogeneous.h>
 #include <CGAL/Filtered_kernel/internal/Static_filters/tools.h>
+#include <CGAL/Nef_3/SNC_intersection.h>
+#include <CGAL/Nef_3/SNC_const_decorator.h>
 
 namespace CGAL {
 
diff --git a/Nef_3/include/CGAL/Nef_3/SHalfedge.h b/Nef_3/include/CGAL/Nef_3/SHalfedge.h
index 31694604135..19339312f20 100644
--- a/Nef_3/include/CGAL/Nef_3/SHalfedge.h
+++ b/Nef_3/include/CGAL/Nef_3/SHalfedge.h
@@ -21,6 +21,7 @@
 
 #include <string>
 #include <sstream>
+#include <CGAL/IO/io.h>
 #include <CGAL/IO/Verbose_ostream.h>
 #include <CGAL/Nef_3/SNC_iteration.h>
 
diff --git a/Nef_3/include/CGAL/Nef_3/SHalfloop.h b/Nef_3/include/CGAL/Nef_3/SHalfloop.h
index 59c449c150f..c611bfeb85d 100644
--- a/Nef_3/include/CGAL/Nef_3/SHalfloop.h
+++ b/Nef_3/include/CGAL/Nef_3/SHalfloop.h
@@ -21,6 +21,7 @@
 
 #include <string>
 #include <sstream>
+#include <CGAL/IO/io.h>
 #include <CGAL/IO/Verbose_ostream.h>
 #include <CGAL/Nef_3/SNC_iteration.h>
 
diff --git a/Nef_3/include/CGAL/Nef_3/SNC_constructor.h b/Nef_3/include/CGAL/Nef_3/SNC_constructor.h
index c79d68de6d5..f0c67661879 100644
--- a/Nef_3/include/CGAL/Nef_3/SNC_constructor.h
+++ b/Nef_3/include/CGAL/Nef_3/SNC_constructor.h
@@ -24,12 +24,13 @@
 #include <CGAL/Nef_S2/Normalizing.h>
 #include <CGAL/Nef_3/bounded_side_3.h>
 #include <CGAL/Nef_3/Pluecker_line_3.h>
-#include <CGAL/Nef_3/SNC_decorator.h>
+#include <CGAL/Nef_3/SNC_const_decorator.h>
 #include <CGAL/Nef_3/SNC_SM_overlayer.h>
 #include <CGAL/Nef_S2/SM_point_locator.h>
+#include <CGAL/Nef_3/SNC_halfedge_key.h>
 #include <CGAL/Nef_3/SNC_sphere_map.h>
+#include <CGAL/Nef_3/SNC_structure.h>
 #include <CGAL/Nef_3/SNC_intersection.h>
-#include <CGAL/Nef_3/SNC_external_structure.h>
 #ifdef SM_VISUALIZOR
 #include <CGAL/Nef_3/SNC_SM_visualizor.h>
 #endif // SM_VISUALIZOR
diff --git a/Nef_3/include/CGAL/Nef_3/SNC_external_structure.h b/Nef_3/include/CGAL/Nef_3/SNC_external_structure.h
index 038843e9fa1..336981c8385 100644
--- a/Nef_3/include/CGAL/Nef_3/SNC_external_structure.h
+++ b/Nef_3/include/CGAL/Nef_3/SNC_external_structure.h
@@ -25,7 +25,7 @@
 #include <CGAL/Nef_3/SNC_point_locator.h>
 #include <CGAL/Nef_S2/SM_point_locator.h>
 #include <CGAL/Nef_3/SNC_FM_decorator.h>
-#include <CGAL/Nef_3/SNC_io_parser.h>
+#include <CGAL/Nef_3/SNC_halfedge_key.h>
 #include <CGAL/Nef_3/SNC_indexed_items.h>
 #include <CGAL/Nef_3/SNC_simplify.h>
 #include <map>
@@ -40,77 +40,6 @@
 
 namespace CGAL {
 
-struct int_lt {
-  bool operator()(const int& i1, const int& i2) const { return i1<i2; }
-};
-template <typename Edge_handle>
-struct Halfedge_key_lt4 {
-
-  bool operator()(const Edge_handle& e1, const Edge_handle& e2) const {
-    if(CGAL::sign(e1->point().x()) != 0) {
-      if(e1->source() != e2->source())
-        return CGAL::compare_x(e1->source()->point(), e2->source()->point()) < 0;
-      else
-        return e1->point().x() < 0;
-    }
-    if(CGAL::sign(e1->point().y()) != 0) {
-      if(e1->source() != e2->source())
-        return CGAL::compare_y(e1->source()->point(), e2->source()->point()) < 0;
-      else
-        return e1->point().y() < 0;
-    }
-    if(e1->source() != e2->source())
-      return CGAL::compare_z(e1->source()->point(), e2->source()->point()) < 0;
-    return e1->point().z() < 0;
-  }
-};
-
-template <typename Edge_handle>
-struct Halfedge_key_lt3 {
-
-  bool operator()(const Edge_handle& e1, const Edge_handle& e2) const {
-    if(e1->source() != e2->source())
-      return CGAL::lexicographically_xyz_smaller(e1->source()->point(), e2->source()->point());
-    if(CGAL::sign(e1->point().x()) != 0)
-      return e1->point().x() < 0;
-    if(CGAL::sign(e1->point().y()) != 0)
-      return e1->point().y() < 0;
-    return e1->point().z() < 0;
-  }
-};
-
-template <typename Point, typename Edge>
-struct Halfedge_key {
-  typedef Halfedge_key<Point,Edge> Self;
-  Point p; int i; Edge e;
-  Halfedge_key(Point pi, int ii, Edge ei) :
-    p(pi), i(ii), e(ei) {}
-  Halfedge_key(const Self& k) : p(k.p), i(k.i), e(k.e) {}
-  Self& operator=(const Self& k) { p=k.p; i=k.i; e=k.e; return *this; }
-  bool operator==(const Self& k) const { return p==k.p && i==k.i; }
-  bool operator!=(const Self& k) const { return !operator==(k); }
-};
-
-template <typename Point, typename Edge, class Decorator>
-struct Halfedge_key_lt {
-  typedef Halfedge_key<Point,Edge> Key;
-  typedef typename Point::R R;
-  typedef typename R::Vector_3 Vector;
-  typedef typename R::Direction_3 Direction;
-  bool operator()( const Key& k1, const Key& k2) const {
-    if( k1.e->source() == k2.e->source())
-      return (k1.i < k2.i);
-    Direction l(k1.e->vector());
-    if( k1.i < 0) l = -l;
-    return (Direction( k2.p - k1.p) == l);
-  }
-};
-
-template <typename Point, typename Edge>
-std::ostream& operator<<(std::ostream& os,
-                         const Halfedge_key<Point,Edge>& k )
-{ os << k.p << " " << k.i; return os; }
-
 template <typename R>
 int sign_of(const CGAL::Plane_3<R>& h)
 { if ( h.c() != 0 ) return CGAL_NTS sign(h.c());
diff --git a/Nef_3/include/CGAL/Nef_3/SNC_halfedge_key.h b/Nef_3/include/CGAL/Nef_3/SNC_halfedge_key.h
new file mode 100644
index 00000000000..4c2738e63f8
--- /dev/null
+++ b/Nef_3/include/CGAL/Nef_3/SNC_halfedge_key.h
@@ -0,0 +1,95 @@
+// Copyright (c) 1997-2002  Max-Planck-Institute Saarbruecken (Germany).
+// 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)     : Peter Hachenberger <hachenberger@mpi-sb.mpg.de>
+
+#ifndef CGAL_SNC_HALFEDGE_KEY_H
+#define CGAL_SNC_HALFEDGE_KEY_H
+
+#include <CGAL/license/Nef_3.h>
+
+#include <CGAL/Kernel/global_functions.h>
+
+namespace CGAL {
+
+struct int_lt {
+  bool operator()(const int& i1, const int& i2) const { return i1<i2; }
+};
+
+template <typename Edge_handle>
+struct Halfedge_key_lt4 {
+
+  bool operator()(const Edge_handle& e1, const Edge_handle& e2) const {
+    if(CGAL::sign(e1->point().x()) != 0) {
+      if(e1->source() != e2->source())
+        return CGAL::compare_x(e1->source()->point(), e2->source()->point()) < 0;
+      else
+        return e1->point().x() < 0;
+    }
+    if(CGAL::sign(e1->point().y()) != 0) {
+      if(e1->source() != e2->source())
+        return CGAL::compare_y(e1->source()->point(), e2->source()->point()) < 0;
+      else
+        return e1->point().y() < 0;
+    }
+    if(e1->source() != e2->source())
+      return CGAL::compare_z(e1->source()->point(), e2->source()->point()) < 0;
+    return e1->point().z() < 0;
+  }
+};
+
+template <typename Edge_handle>
+struct Halfedge_key_lt3 {
+
+  bool operator()(const Edge_handle& e1, const Edge_handle& e2) const {
+    if(e1->source() != e2->source())
+      return CGAL::lexicographically_xyz_smaller(e1->source()->point(), e2->source()->point());
+    if(CGAL::sign(e1->point().x()) != 0)
+      return e1->point().x() < 0;
+    if(CGAL::sign(e1->point().y()) != 0)
+      return e1->point().y() < 0;
+    return e1->point().z() < 0;
+  }
+};
+
+template <typename Point, typename Edge>
+struct Halfedge_key {
+  typedef Halfedge_key<Point,Edge> Self;
+  Point p; int i; Edge e;
+  Halfedge_key(Point pi, int ii, Edge ei) :
+    p(pi), i(ii), e(ei) {}
+  Halfedge_key(const Self& k) : p(k.p), i(k.i), e(k.e) {}
+  Self& operator=(const Self& k) { p=k.p; i=k.i; e=k.e; return *this; }
+  bool operator==(const Self& k) const { return p==k.p && i==k.i; }
+  bool operator!=(const Self& k) const { return !operator==(k); }
+};
+
+template <typename Point, typename Edge, class Decorator>
+struct Halfedge_key_lt {
+  typedef Halfedge_key<Point,Edge> Key;
+  typedef typename Point::R R;
+  typedef typename R::Vector_3 Vector;
+  typedef typename R::Direction_3 Direction;
+  bool operator()( const Key& k1, const Key& k2) const {
+    if( k1.e->source() == k2.e->source())
+      return (k1.i < k2.i);
+    Direction l(k1.e->vector());
+    if( k1.i < 0) l = -l;
+    return (Direction( k2.p - k1.p) == l);
+  }
+};
+
+template <typename Point, typename Edge>
+std::ostream& operator<<(std::ostream& os,
+                         const Halfedge_key<Point,Edge>& k )
+{ os << k.p << " " << k.i; return os; }
+
+}
+#endif //CGAL_SNC_HALFEDGE_KEY_H
diff --git a/Nef_3/include/CGAL/Nef_3/SNC_io_parser.h b/Nef_3/include/CGAL/Nef_3/SNC_io_parser.h
index aab28de6795..9aa445f38a1 100644
--- a/Nef_3/include/CGAL/Nef_3/SNC_io_parser.h
+++ b/Nef_3/include/CGAL/Nef_3/SNC_io_parser.h
@@ -24,7 +24,6 @@
 #include <CGAL/Nef_S2/SM_decorator.h>
 #include <CGAL/Nef_3/SNC_structure.h>
 #include <CGAL/Nef_3/SNC_decorator.h>
-#include <CGAL/Nef_3/SNC_constructor.h>
 #include <CGAL/Nef_2/Object_index.h>
 #include <CGAL/Nef_S2/Normalizing.h>
 #include <vector>
diff --git a/Nef_3/include/CGAL/Nef_3/SNC_sphere_map.h b/Nef_3/include/CGAL/Nef_3/SNC_sphere_map.h
index dbdb461db60..0b84b289122 100644
--- a/Nef_3/include/CGAL/Nef_3/SNC_sphere_map.h
+++ b/Nef_3/include/CGAL/Nef_3/SNC_sphere_map.h
@@ -24,6 +24,7 @@
 #include <CGAL/Nef_S2/Generic_handle_map.h>
 #include <CGAL/Nef_2/iterator_tools.h>
 #include <CGAL/Nef_3/Infimaximal_box.h>
+#include <CGAL/Nef_3/SNC_list.h>
 #include <CGAL/Nef_S2/Sphere_geometry.h>
 #include <list>
 #undef CGAL_NEF_DEBUG
diff --git a/Nef_3/include/CGAL/Nef_3/SNC_structure.h b/Nef_3/include/CGAL/Nef_3/SNC_structure.h
index c6dccf530a3..27a3c5c414d 100644
--- a/Nef_3/include/CGAL/Nef_3/SNC_structure.h
+++ b/Nef_3/include/CGAL/Nef_3/SNC_structure.h
@@ -56,6 +56,7 @@ void merge_sets( Object o1, Object o2, Hash_map& hash, Union_find& uf) {
 template <typename K, typename I, typename M> class SNC_sphere_map;
 template <typename S> class SM_decorator;
 template <typename S> class SNC_decorator;
+template <typename S> class SNC_io_parser;
 
 /*{\Manpage {SNC_structure}{Items}{Selective Nef Complex}{C}}*/
 
diff --git a/Nef_3/include/CGAL/Nef_3/binop_intersection_tests.h b/Nef_3/include/CGAL/Nef_3/binop_intersection_tests.h
index 03773258a68..b64fe1e99d2 100644
--- a/Nef_3/include/CGAL/Nef_3/binop_intersection_tests.h
+++ b/Nef_3/include/CGAL/Nef_3/binop_intersection_tests.h
@@ -17,9 +17,8 @@
 
 #include <CGAL/Nef_3/Nef_box.h>
 #include <CGAL/Nef_3/Infimaximal_box.h>
+#include <CGAL/Nef_3/SNC_const_decorator.h>
 #include <vector>
-#include <iostream>
-#include <CGAL/Timer.h>
 
 namespace CGAL {
 
@@ -41,21 +40,6 @@ struct binop_intersection_test_segment_tree {
   struct Bop_edge0_face1_callback {
     Callback           &cb;
 
-    struct Pair_hash_function {
-      typedef std::size_t result_type;
-
-      template <class H>
-      std::size_t
-      operator() (const H& h) const {
-        return
-          std::size_t(&*(h.first)) / sizeof
-          (typename std::iterator_traits<typename H::first_type>::value_type)
-              +
-          std::size_t(&*(h.second)) / sizeof
-          (typename std::iterator_traits<typename H::second_type>::value_type);
-      }
-    };
-
     Bop_edge0_face1_callback(Callback &cb)
     : cb(cb)
     {}
@@ -76,7 +60,6 @@ struct binop_intersection_test_segment_tree {
     }
   };
 
-
   template<class Callback>
   struct Bop_edge1_face0_callback {
     Callback         &cb;
@@ -126,55 +109,40 @@ struct binop_intersection_test_segment_tree {
   template<class Callback>
   void operator()(Callback& cb0,
                   Callback& cb1,
-                  SNC_structure& sncp,
-                  SNC_structure& snc1i)
+                  const SNC_structure& snc0,
+                  const SNC_structure& snc1)
   {
     Halfedge_iterator e0, e1;
     Halffacet_iterator f0, f1;
-    std::vector<Nef_box> a, b;
+    std::vector<Nef_box> e0boxes, e1boxes, f0boxes, f1boxes;
+
+    e0boxes.reserve(snc0.number_of_halfedges());
+    e1boxes.reserve(snc1.number_of_halfedges());
+    f0boxes.reserve(snc0.number_of_halffacets());
+    f1boxes.reserve(snc1.number_of_halffacets());
+
+    CGAL_forall_edges( e0, snc0) e0boxes.push_back( Nef_box( e0 ) );
+    CGAL_forall_edges( e1, snc1) e1boxes.push_back( Nef_box( e1 ) );
+    CGAL_forall_facets( f0, snc0) f0boxes.push_back( Nef_box( f0 ) );
+    CGAL_forall_facets( f1, snc1) f1boxes.push_back( Nef_box( f1 ) );
 
     CGAL_NEF_TRACEN("start edge0 edge1");
     Bop_edge0_edge1_callback<Callback> callback_edge0_edge1( cb0 );
-    CGAL_forall_edges( e0, sncp)  a.push_back( Nef_box( e0 ) );
-    CGAL_forall_edges( e1, snc1i) b.push_back( Nef_box( e1 ) );
-#ifdef CGAL_NEF3_BOX_INTERSECTION_CUTOFF
-    box_intersection_d( a.begin(), a.end(), b.begin(), b.end(),
-                        callback_edge0_edge1,
-                        CGAL_NEF3_BOX_INTERSECTION_CUTOFF,);
-#else
-    box_intersection_d( a.begin(), a.end(), b.begin(), b.end(),
+    box_intersection_d( e0boxes.begin(), e0boxes.end(),
+                        e1boxes.begin(), e1boxes.end(),
                         callback_edge0_edge1);
-#endif
-    a.clear();
-    b.clear();
 
     CGAL_NEF_TRACEN("start edge0 face1");
     Bop_edge0_face1_callback<Callback> callback_edge0_face1( cb0 );
-    CGAL_forall_edges( e0, sncp ) a.push_back( Nef_box( e0 ) );
-    CGAL_forall_facets( f1, snc1i)    b.push_back( Nef_box( f1 ) );
-#ifdef CGAL_NEF3_BOX_INTERSECTION_CUTOFF
-    box_intersection_d( a.begin(), a.end(), b.begin(), b.end(),
-                        callback_edge0_face1,
-                        CGAL_NEF3_BOX_INTERSECTION_CUTOFF);
-#else
-    box_intersection_d( a.begin(), a.end(), b.begin(), b.end(),
+    box_intersection_d( e0boxes.begin(), e0boxes.end(),
+                        f1boxes.begin(), f1boxes.end(),
                         callback_edge0_face1);
-#endif
-    a.clear();
-    b.clear();
 
     CGAL_NEF_TRACEN("start edge1 face0");
     Bop_edge1_face0_callback<Callback> callback_edge1_face0( cb1 );
-    CGAL_forall_edges( e1, snc1i)  a.push_back( Nef_box( e1 ) );
-    CGAL_forall_facets( f0, sncp ) b.push_back( Nef_box( f0 ) );
-#ifdef CGAL_NEF3_BOX_INTERSECTION_CUTOFF
-    box_intersection_d( a.begin(), a.end(), b.begin(), b.end(),
-                        callback_edge1_face0,
-                        CGAL_NEF3_BOX_INTERSECTION_CUTOFF);
-#else
-    box_intersection_d( a.begin(), a.end(), b.begin(), b.end(),
+    box_intersection_d( e1boxes.begin(), e1boxes.end(),
+                        f0boxes.begin(), f0boxes.end(),
                         callback_edge1_face0);
-#endif
   }
 };
 
diff --git a/Nef_3/include/CGAL/Nef_3/polygon_mesh_to_nef_3.h b/Nef_3/include/CGAL/Nef_3/polygon_mesh_to_nef_3.h
index c0076cfbca5..a20a6643f50 100644
--- a/Nef_3/include/CGAL/Nef_3/polygon_mesh_to_nef_3.h
+++ b/Nef_3/include/CGAL/Nef_3/polygon_mesh_to_nef_3.h
@@ -19,10 +19,13 @@
 
 #include <CGAL/license/Nef_3.h>
 
-
+#include <CGAL/Circulator_project.h>
 #include <CGAL/normal_vector_newell_3.h>
 #include <CGAL/Nef_S2/SM_point_locator.h>
 #include <CGAL/Nef_3/SNC_indexed_items.h>
+#include <CGAL/Plane_3.h>
+#include <CGAL/Point_3.h>
+#include <CGAL/Vector_3.h>
 #include <CGAL/boost/graph/helpers.h>
 
 #undef CGAL_NEF_DEBUG
diff --git a/Nef_3/include/CGAL/Nef_3/vertex_cycle_to_nef_3.h b/Nef_3/include/CGAL/Nef_3/vertex_cycle_to_nef_3.h
index 1b7426fee5a..750fcf86ffe 100644
--- a/Nef_3/include/CGAL/Nef_3/vertex_cycle_to_nef_3.h
+++ b/Nef_3/include/CGAL/Nef_3/vertex_cycle_to_nef_3.h
@@ -28,10 +28,7 @@
 #include <CGAL/Constrained_triangulation_plus_2.h>
 
 // Nef polyhedra
-#include <CGAL/Nef_polyhedron_3.h>
-#include <CGAL/Nef_3/SNC_structure.h>
-#include <CGAL/Nef_3/SNC_constructor.h>
-#include <CGAL/Nef_3/SNC_point_locator.h>
+#include <CGAL/Nef_3/SNC_indexed_items.h>
 
 namespace CGAL {
 
diff --git a/Polygon/include/CGAL/Polygon_2.h b/Polygon/include/CGAL/Polygon_2.h
index e4305cccfed..b2b7ab7340a 100644
--- a/Polygon/include/CGAL/Polygon_2.h
+++ b/Polygon/include/CGAL/Polygon_2.h
@@ -242,8 +242,11 @@ class Polygon_2 {
     /// Erases the vertex pointed to by `i`.
     Vertex_circulator erase(Vertex_circulator i)
       {
-        return Vertex_circulator(&d_container,
-                                 d_container.erase(i.mod_iterator()));
+        auto it = d_container.erase(i.mod_iterator());
+        if(it == d_container.end()){
+          it = d_container.begin();
+        }
+        return Vertex_circulator(&d_container, it);
       }
 
     /// Erases the vertices in the range `[first, last)`.
diff --git a/Polygon/test/Polygon/issue7228.cpp b/Polygon/test/Polygon/issue7228.cpp
new file mode 100644
index 00000000000..623a8633130
--- /dev/null
+++ b/Polygon/test/Polygon/issue7228.cpp
@@ -0,0 +1,29 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Polygon_2.h>
+
+#include <vector>
+#include <array>
+#include <cassert>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_2 Point;
+typedef CGAL::Polygon_2<K> Polygon_2;
+typedef Polygon_2::Vertex_circulator Vertex_circulator;
+
+int main()
+{
+  std::array<Point,4> points =  { Point(0,0), Point(1,0), Point(1,1), Point(0,1) };
+  Polygon_2 poly(points.begin(), points.end());
+
+  Vertex_circulator vc = poly.vertices_circulator();
+
+  ++vc;
+  ++vc;
+  ++vc;
+
+  vc = poly.erase(vc);
+
+  assert(*vc == Point(0,0));
+
+  return 0;
+}
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp
index 048739c136f..159ff7cd708 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp
@@ -43,7 +43,7 @@ int main(int argc, char* argv[])
   }
 
   double target_edge_length = (argc > 2) ? std::stod(std::string(argv[2])) : 0.04;
-  unsigned int nb_iter = 3;
+  unsigned int nb_iter = (argc > 3) ? std::stoi(std::string(argv[3])) : 10;
 
   std::cout << "Split border...";
 
@@ -59,6 +59,8 @@ int main(int argc, char* argv[])
                            CGAL::parameters::number_of_iterations(nb_iter)
                                             .protect_constraints(true)); //i.e. protect border, here
 
+  CGAL::IO::write_polygon_mesh("out.off", mesh, CGAL::parameters::stream_precision(17));
+
   std::cout << "Remeshing done." << std::endl;
 
   return 0;
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/refine_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/refine_impl.h
index 54369aee910..06ddc868f0d 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/refine_impl.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/refine_impl.h
@@ -15,12 +15,6 @@
 
 #include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
 
-
-#include <cmath>
-#include <map>
-#include <set>
-#include <list>
-
 #include <CGAL/assertions.h>
 #ifdef CGAL_PMP_FAIR_DEBUG
 #include <CGAL/Timer.h>
@@ -30,8 +24,14 @@
 #include <CGAL/boost/graph/iterator.h>
 #include <CGAL/boost/graph/Euler_operations.h>
 #include <CGAL/boost/graph/properties.h>
+#include <CGAL/Polygon_mesh_processing/repair_degeneracies.h>
 #include <CGAL/property_map.h>
 
+#include <cmath>
+#include <map>
+#include <set>
+#include <list>
+
 namespace CGAL {
 
 namespace Polygon_mesh_processing {
@@ -49,15 +49,28 @@ class Refine_Polyhedron_3 {
   typedef Halfedge_around_face_circulator<PolygonMesh>    Halfedge_around_facet_circulator;
   typedef Halfedge_around_target_circulator<PolygonMesh>  Halfedge_around_vertex_circulator;
 
+  typedef typename CGAL::Kernel_traits<Point_3>::type                     Traits;
+
 private:
   PolygonMesh& pmesh;
   VertexPointMap vpmap;
+  Traits traits = {};
 
-  bool flippable(halfedge_descriptor h) {
+  bool flippable(halfedge_descriptor h)
+  {
     // this check is added so that edge flip does not break manifoldness
     // it might happen when there is an edge where flip_edge(h) will be placed (i.e. two edges collide after flip)
     vertex_descriptor v_tip_0 = target(next(h,pmesh),pmesh);
     vertex_descriptor v_tip_1 = target(next(opposite(h,pmesh),pmesh),pmesh);
+
+#ifdef CGAL_PMP_REFINE_DEBUG_PP
+    std::cout << "flippable() " << h << std::endl;
+    std::cout << "\t" << source(h, pmesh) << ": " << pmesh.point(source(h, pmesh)) << std::endl;
+    std::cout << "\t" << target(h, pmesh) << ": " << pmesh.point(target(h, pmesh)) << std::endl;
+    std::cout << "\t" << v_tip_0 << ": " << pmesh.point(v_tip_0) << std::endl;
+    std::cout << "\t" << v_tip_1 << ": " << pmesh.point(v_tip_1) << std::endl;
+#endif
+
     Halfedge_around_vertex_circulator v_cir(next(h,pmesh), pmesh), v_end(v_cir);
     do {
       if(target(opposite(*v_cir, pmesh),pmesh) == v_tip_1) { return false; }
@@ -74,13 +87,21 @@ private:
 
   bool relax(halfedge_descriptor h)
   {
+#ifdef CGAL_PMP_REFINE_DEBUG_PP
     typedef typename boost::property_traits<VertexPointMap>::reference Point_3_ref;
-    Point_3_ref p = get(vpmap, target(h,pmesh));
-    Point_3_ref q = get(vpmap, target(opposite(h,pmesh),pmesh));
+    Point_3_ref p = get(vpmap, source(h,pmesh));
+    Point_3_ref q = get(vpmap, target(h,pmesh));
     Point_3_ref r = get(vpmap, target(next(h,pmesh),pmesh));
     Point_3_ref s = get(vpmap, target(next(opposite(h,pmesh),pmesh),pmesh));
-    if( (CGAL::ON_UNBOUNDED_SIDE  != CGAL::side_of_bounded_sphere(p,q,r,s)) ||
-        (CGAL::ON_UNBOUNDED_SIDE  != CGAL::side_of_bounded_sphere(p,q,s,r)) )
+
+    std::cout << "relax() " << h << std::endl;
+    std::cout << "\t" << source(h, pmesh) << ": " << p << std::endl;
+    std::cout << "\t" << target(h, pmesh) << ": " << q << std::endl;
+    std::cout << "\t" << target(next(h,pmesh),pmesh) << ": " << r << std::endl;
+    std::cout << "\t" << target(next(opposite(h,pmesh),pmesh),pmesh) << ": " << s << std::endl;
+#endif
+
+    if(internal::should_flip(edge(h, pmesh), pmesh, vpmap, traits))
     {
       if(flippable(h)) {
         Euler::flip_edge(h,pmesh);
@@ -238,8 +259,7 @@ private:
       Halfedge_around_face_circulator<PolygonMesh> circ(halfedge(fd,pmesh),pmesh), done(circ);
       do {
         vertex_descriptor v = target(*circ,pmesh);
-        std::pair<typename std::map<vertex_descriptor, double>::iterator, bool> v_insert
-          = scale_attribute.insert(std::make_pair(v, 0));
+        auto v_insert = scale_attribute.emplace(v, 0);
         if(!v_insert.second) { continue; } // already calculated
         v_insert.first->second = average_length(v, interior_map, accept_internal_facets);
       } while(++circ != done);
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h
index 2340993b99d..ac5a340d816 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h
@@ -24,6 +24,7 @@
 #include <CGAL/Named_function_parameters.h>
 #include <CGAL/boost/graph/named_params_helper.h>
 
+#include <CGAL/Polygon_mesh_processing/border.h>
 
 #include <CGAL/Lazy.h> // needed for CGAL::exact(FT)/CGAL::exact(Lazy_exact_nt<T>)
 
@@ -31,6 +32,7 @@
 #include <boost/graph/graph_traits.hpp>
 #include <boost/dynamic_bitset.hpp>
 
+#include <vector>
 #include <utility>
 #include <algorithm>
 #include <unordered_set>
@@ -308,6 +310,7 @@ face_border_length(typename boost::graph_traits<PolygonMesh>::halfedge_descripto
   *   - `first`: a halfedge on the longest border.
   *     The return type `halfedge_descriptor` is a halfedge descriptor. It is
   *     deduced from the graph traits corresponding to the type `PolygonMesh`.
+  *     `first` is among the halfedges reported by `extract_boundary_cycles()`.
   *   - `second`: the length of the longest border
   *     The return type `FT` is a number type either deduced from the `geom_traits`
   *     \ref bgl_namedparameters "Named Parameters" if provided,
@@ -318,6 +321,7 @@ face_border_length(typename boost::graph_traits<PolygonMesh>::halfedge_descripto
   * will be performed approximately.
   *
   * @see `face_border_length()`
+  * @see `extract_boundary_cycles()`
   */
 template<typename PolygonMesh,
          typename NamedParameters = parameters::Default_named_parameters>
@@ -334,28 +338,18 @@ longest_border(const PolygonMesh& pmesh,
             typename property_map_value<PolygonMesh, CGAL::vertex_point_t>::type>::Kernel::FT  FT;
   typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor                       halfedge_descriptor;
 
-  std::unordered_set<halfedge_descriptor> visited;
+  std::vector<halfedge_descriptor> boundary_cycles;
+  extract_boundary_cycles(pmesh, std::back_inserter(boundary_cycles));
   halfedge_descriptor result_halfedge = boost::graph_traits<PolygonMesh>::null_halfedge();
   FT result_len = 0;
-  for(halfedge_descriptor h : halfedges(pmesh))
+  for(halfedge_descriptor h : boundary_cycles)
   {
-    if(visited.find(h)== visited.end())
-    {
-      if(is_border(h, pmesh))
-      {
-        FT len = 0;
-        for(halfedge_descriptor haf : halfedges_around_face(h, pmesh))
-        {
-          len += edge_length(haf, pmesh, np);
-          visited.insert(haf);
-        }
+    FT len = face_border_length(h, pmesh, np);
 
-        if(result_len < len)
-        {
-          result_len = len;
-          result_halfedge = h;
-        }
-      }
+    if(result_len < len)
+    {
+      result_len = len;
+      result_halfedge = h;
     }
   }
   return std::make_pair(result_halfedge, result_len);
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_degeneracies.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_degeneracies.h
index bb7e48292a6..798499f430b 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_degeneracies.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_degeneracies.h
@@ -294,7 +294,6 @@ get_best_edge_orientation(typename boost::graph_traits<TriangleMesh>::edge_descr
   return boost::graph_traits<TriangleMesh>::null_halfedge();
 }
 
-// adapted from triangulate_faces
 template <typename TriangleMesh, typename VPM, typename Traits>
 bool should_flip(typename boost::graph_traits<TriangleMesh>::edge_descriptor e,
                  const TriangleMesh& tmesh,
@@ -303,49 +302,23 @@ bool should_flip(typename boost::graph_traits<TriangleMesh>::edge_descriptor e,
 {
   typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
 
-  typedef typename Traits::FT                                             FT;
+  typedef typename Traits:: FT                                            FT;
   typedef typename boost::property_traits<VPM>::reference                 Point_ref;
-  typedef typename Traits::Vector_3                                       Vector_3;
 
   CGAL_precondition(!is_border(e, tmesh));
 
-  halfedge_descriptor h = halfedge(e, tmesh);
+  typename Traits::Compute_approximate_angle_3 angle = gt.compute_approximate_angle_3_object();
 
-  Point_ref p0 = get(vpm, target(h, tmesh));
-  Point_ref p1 = get(vpm, target(next(h, tmesh), tmesh));
-  Point_ref p2 = get(vpm, source(h, tmesh));
-  Point_ref p3 = get(vpm, target(next(opposite(h, tmesh), tmesh), tmesh));
+  const halfedge_descriptor h = halfedge(e, tmesh);
 
-  /* Chooses the diagonal that will split the quad in two triangles that maximize
-   * the scalar product of of the un-normalized normals of the two triangles.
-   * The lengths of the un-normalized normals (computed using cross-products of two vectors)
-   *  are proportional to the area of the triangles.
-   * Maximize the scalar product of the two normals will avoid skinny triangles,
-   * and will also taken into account the cosine of the angle between the two normals.
-   * In particular, if the two triangles are oriented in different directions,
-   * the scalar product will be negative.
-   */
+  const Point_ref p0 = get(vpm, target(h, tmesh));
+  const Point_ref p1 = get(vpm, target(next(h, tmesh), tmesh));
+  const Point_ref p2 = get(vpm, source(h, tmesh));
+  const Point_ref p3 = get(vpm, target(next(opposite(h, tmesh), tmesh), tmesh));
 
-//  CGAL::cross_product(p2-p1, p3-p2) * CGAL::cross_product(p0-p3, p1-p0);
-//  CGAL::cross_product(p1-p0, p1-p2) * CGAL::cross_product(p3-p2, p3-p0);
-
-  const Vector_3 v01 = gt.construct_vector_3_object()(p0, p1);
-  const Vector_3 v12 = gt.construct_vector_3_object()(p1, p2);
-  const Vector_3 v23 = gt.construct_vector_3_object()(p2, p3);
-  const Vector_3 v30 = gt.construct_vector_3_object()(p3, p0);
-
-  const FT p1p3 = gt.compute_scalar_product_3_object()(
-                    gt.construct_cross_product_vector_3_object()(v12, v23),
-                    gt.construct_cross_product_vector_3_object()(v30, v01));
-
-  const Vector_3 v21 = gt.construct_opposite_vector_3_object()(v12);
-  const Vector_3 v03 = gt.construct_opposite_vector_3_object()(v30);
-
-  const FT p0p2 = gt.compute_scalar_product_3_object()(
-                    gt.construct_cross_product_vector_3_object()(v01, v21),
-                    gt.construct_cross_product_vector_3_object()(v23, v03));
-
-  return p0p2 <= p1p3;
+  const FT ap1 = angle(p0,p1,p2);
+  const FT ap3 = angle(p2,p3,p0);
+  return (ap1 + ap3 > FT(180));
 }
 
 template <class TriangleMesh, class VPM, class Traits, class Functor>
diff --git a/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_model_complexity_control.cpp b/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_model_complexity_control.cpp
index e132d62975d..cf1ce77b8b1 100644
--- a/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_model_complexity_control.cpp
+++ b/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_model_complexity_control.cpp
@@ -39,9 +39,9 @@ typedef CGAL::Nth_of_tuple_property_map<2, PNI>
 * candidate generation are cached and reused.
 */
 
-int main()
+int main(int argc, char* argv[])
 {
-    const std::string& input_file(CGAL::data_file_path("points_3/building.ply"));
+  const std::string input_file = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/building.ply");
   std::ifstream input_stream(input_file.c_str());
 
   std::vector<PNI> points; // store points
diff --git a/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_user_provided_planes.cpp b/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_user_provided_planes.cpp
index 2f2ed2df2c9..5db8d3f3700 100644
--- a/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_user_provided_planes.cpp
+++ b/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_user_provided_planes.cpp
@@ -39,9 +39,9 @@ typedef CGAL::Nth_of_tuple_property_map<2, PNI>
 * the point is not assigned to a plane).
 */
 
-int main()
+int main(int argc, char* argv[])
 {
-    const std::string& input_file(CGAL::data_file_path("points_3/ball.ply"));
+  const std::string input_file = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/ball.ply");
   std::ifstream input_stream(input_file.c_str());
 
   std::vector<PNI> points; // store points
diff --git a/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_with_region_growing.cpp b/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_with_region_growing.cpp
index 493c8c138aa..530d1422de4 100644
--- a/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_with_region_growing.cpp
+++ b/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_with_region_growing.cpp
@@ -83,13 +83,13 @@ private:
 * the surface model from the planes.
 */
 
-int main()
+int main(int argc, char* argv[])
 {
   Point_vector points;
 
   // Load point set from a file.
-  const std::string input_file(CGAL::data_file_path("points_3/cube.pwn"));
-    std::ifstream input_stream(input_file.c_str());
+  const std::string input_file = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/cube.pwn");
+  std::ifstream input_stream(input_file.c_str());
   if (input_stream.fail()) {
     std::cerr << "Failed open file \'" << input_file << "\'" << std::endl;
     return EXIT_FAILURE;
diff --git a/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_without_input_planes.cpp b/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_without_input_planes.cpp
index 4b2ed8f91a4..d14411babd1 100644
--- a/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_without_input_planes.cpp
+++ b/Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_without_input_planes.cpp
@@ -52,13 +52,13 @@ typedef CGAL::Surface_mesh<Point>
 * the surface model from the planes.
 */
 
-int main()
+int main(int argc, char* argv[])
 {
   Point_vector points;
 
   // Loads point set from a file.
-  const std::string input_file(CGAL::data_file_path("points_3/cube.pwn"));
-    std::ifstream input_stream(input_file.c_str());
+  const std::string input_file = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/cube.pwn");
+  std::ifstream input_stream(input_file.c_str());
   if (input_stream.fail()) {
     std::cerr << "failed open file \'" <<input_file << "\'" << std::endl;
     return EXIT_FAILURE;
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp
index ee9c3fa0e69..a15e1f14dfd 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp
@@ -347,8 +347,7 @@ public Q_SLOTS:
       }
       // Create dialog box
       QDialog dialog(mw);
-      Ui::Isotropic_remeshing_dialog ui
-        = remeshing_dialog(&dialog, poly_item, selection_item);
+      initialize_remeshing_dialog(&dialog, poly_item, selection_item);
 
       // Get values
       int i = dialog.exec();
@@ -357,6 +356,7 @@ public Q_SLOTS:
         std::cout << "Remeshing aborted" << std::endl;
         return;
       }
+
       bool edges_only = ui.splitEdgesOnly_checkbox->isChecked();
       bool preserve_duplicates = ui.preserveDuplicates_checkbox->isChecked();
       double target_length = ui.edgeLength_dspinbox->value();
@@ -710,7 +710,7 @@ public Q_SLOTS:
         if (target_length == 0.)//parameters have not been set yet
         {
         QDialog dialog(mw);
-        Ui::Isotropic_remeshing_dialog ui = remeshing_dialog(&dialog, poly_item);
+        initialize_remeshing_dialog(&dialog, poly_item);
         ui.objectName->setText(QString::number(scene->selectionIndices().size())
           .append(QString(" items to be remeshed")));
         int i = dialog.exec();
@@ -937,32 +937,73 @@ private:
   };
 #endif
 
-  Ui::Isotropic_remeshing_dialog
-  remeshing_dialog(QDialog* dialog,
-                   Scene_facegraph_item* poly_item,
-                   Scene_polyhedron_selection_item* selection_item = nullptr)
+public Q_SLOTS:
+  void update_after_protect_checkbox_click()
+  {
+    if(ui.protect_checkbox->isChecked())
+    {
+      ui.smooth1D_label->setEnabled(false);
+      ui.smooth1D_checkbox->setEnabled(false);
+      ui.smooth1D_checkbox->setChecked(false);
+    }
+    else
+    {
+      ui.smooth1D_label->setEnabled(true);
+      ui.smooth1D_checkbox->setEnabled(true);
+    }
+  }
+
+  void update_after_splitEdgesOnly_click()
+  {
+    if(ui.splitEdgesOnly_checkbox->isChecked())
+    {
+      ui.nbIterations_label->setEnabled(false);
+      ui.nbIterations_spinbox->setEnabled(false);
+      ui.nbSmoothing_label->setEnabled(false);
+      ui.nbSmoothing_spinbox->setEnabled(false);
+
+      ui.protect_label->setEnabled(false);
+      ui.protect_checkbox->setEnabled(false);
+      ui.protect_checkbox->setChecked(false);
+
+      ui.smooth1D_label->setEnabled(false);
+      ui.smooth1D_checkbox->setEnabled(false);
+      ui.smooth1D_checkbox->setChecked(false);
+    }
+    else
+    {
+      ui.nbIterations_label->setEnabled(true);
+      ui.nbIterations_spinbox->setEnabled(true);
+      ui.nbSmoothing_label->setEnabled(true);
+      ui.nbSmoothing_spinbox->setEnabled(true);
+
+      ui.protect_label->setEnabled(true);
+      ui.protect_checkbox->setEnabled(true);
+
+      ui.smooth1D_label->setEnabled(true);
+      ui.smooth1D_checkbox->setEnabled(true);
+    }
+  }
+
+public:
+  void
+  initialize_remeshing_dialog(QDialog* dialog,
+                              Scene_facegraph_item* poly_item,
+                              Scene_polyhedron_selection_item* selection_item = nullptr)
   {
-    Ui::Isotropic_remeshing_dialog ui;
     ui.setupUi(dialog);
     connect(ui.buttonBox, SIGNAL(accepted()), dialog, SLOT(accept()));
     connect(ui.buttonBox, SIGNAL(rejected()), dialog, SLOT(reject()));
 
     //connect checkbox to spinbox
-    connect(ui.splitEdgesOnly_checkbox, SIGNAL(toggled(bool)),
-            ui.nbIterations_spinbox, SLOT(setDisabled(bool)));
-    connect(ui.splitEdgesOnly_checkbox, SIGNAL(toggled(bool)),
-            ui.protect_checkbox, SLOT(setDisabled(bool)));
-    connect(ui.splitEdgesOnly_checkbox, SIGNAL(toggled(bool)),
-            ui.smooth1D_checkbox, SLOT(setDisabled(bool)));
-    connect(ui.splitEdgesOnly_checkbox, SIGNAL(toggled(bool)),
-            ui.nbSmoothing_spinbox, SLOT(setDisabled(bool)));
-    connect(ui.protect_checkbox, SIGNAL(toggled(bool)),
-            ui.smooth1D_checkbox, SLOT(setDisabled(bool)));
     connect(ui.preserveDuplicates_checkbox, SIGNAL(toggled(bool)),
             ui.protect_checkbox, SLOT(setChecked(bool)));
     connect(ui.preserveDuplicates_checkbox, SIGNAL(toggled(bool)),
             ui.protect_checkbox, SLOT(setDisabled(bool)));
 
+    connect(ui.protect_checkbox, SIGNAL(clicked(bool)), this, SLOT(update_after_protect_checkbox_click()));
+    connect(ui.splitEdgesOnly_checkbox, SIGNAL(clicked(bool)), this, SLOT(update_after_splitEdgesOnly_click()));
+
     //Set default parameters
     Scene_interface::Bbox bbox = poly_item != nullptr ? poly_item->bbox()
       : (selection_item != nullptr ? selection_item->bbox()
@@ -1003,14 +1044,12 @@ private:
       ui.preserveDuplicates_checkbox->setDisabled(true);
       ui.preserveDuplicates_checkbox->setChecked(false);
     }
-
-    return ui;
   }
 
 
 private:
   QAction* actionIsotropicRemeshing_;
-
+  Ui::Isotropic_remeshing_dialog ui;
 }; // end Polyhedron_demo_isotropic_remeshing_plugin
 
 #include "Isotropic_remeshing_plugin.moc"
diff --git a/Polyhedron/demo/Polyhedron/Scene.cpp b/Polyhedron/demo/Polyhedron/Scene.cpp
index b148e8df46b..c70929ab4a5 100644
--- a/Polyhedron/demo/Polyhedron/Scene.cpp
+++ b/Polyhedron/demo/Polyhedron/Scene.cpp
@@ -1928,6 +1928,7 @@ void Scene::removeViewer(Viewer_interface *viewer)
   if(viewer->property("is_destroyed").toBool())
     return;
 
+  viewer->makeCurrent();
   vaos[viewer]->destroy();
   vaos[viewer]->deleteLater();
   vaos.remove(viewer);
diff --git a/STL_Extension/include/CGAL/exceptions.h b/STL_Extension/include/CGAL/exceptions.h
index edb926383f9..3b269aea302 100644
--- a/STL_Extension/include/CGAL/exceptions.h
+++ b/STL_Extension/include/CGAL/exceptions.h
@@ -14,22 +14,9 @@
 #ifndef CGAL_EXCEPTIONS_H
 #define CGAL_EXCEPTIONS_H
 
-#include <CGAL/config.h>
 #include <stdexcept>
 #include <string>
 
-// Address the warning C4003: not enough actual parameters for macro 'BOOST_PP_SEQ_DETAIL_IS_NOT_EMPTY'
-// lexical_cast.hpp includes files from boost/preprocessor
-// This concerns boost 1_67_0
-#if defined(BOOST_MSVC)
-#  pragma warning(push)
-#  pragma warning(disable: 4003)
-#endif
-#include <boost/lexical_cast.hpp>
-#if defined(BOOST_MSVC)
-#  pragma warning(pop)
-#endif
-
 
 namespace CGAL {
 
@@ -87,7 +74,7 @@ public:
     std::logic_error( lib + std::string( " ERROR: ") + kind + std::string( "!")
         + ((expr.empty()) ? (std::string("")) : (std::string("\nExpr: ")+expr))
         + std::string( "\nFile: ") + file
-        + std::string( "\nLine: ") + boost::lexical_cast<std::string>(line)
+        + std::string( "\nLine: ") + std::to_string(line)
         + ((msg.empty()) ? (std::string(""))
                          : (std::string("\nExplanation: ") + msg))),
     m_lib( lib),
diff --git a/Surface_mesh/examples/Surface_mesh/draw_surface_mesh.cpp b/Surface_mesh/examples/Surface_mesh/draw_surface_mesh.cpp
index cbd42788fcc..488282915ba 100644
--- a/Surface_mesh/examples/Surface_mesh/draw_surface_mesh.cpp
+++ b/Surface_mesh/examples/Surface_mesh/draw_surface_mesh.cpp
@@ -20,6 +20,25 @@ int main(int argc, char* argv[])
     return EXIT_FAILURE;
   }
 
+  // Internal color property maps are used if they exist and are called "v:color", "e:color" and "f:color".
+  auto vcm = sm.add_property_map<Mesh::Vertex_index, CGAL::IO::Color>("v:color").first;
+  auto ecm = sm.add_property_map<Mesh::Edge_index, CGAL::IO::Color>("e:color").first;
+  auto fcm = sm.add_property_map<Mesh::Face_index>("f:color", CGAL::IO::white() /*default*/).first;
+
+  for(auto v : vertices(sm))
+  {
+    if(v.idx()%2)
+      put(vcm, v, CGAL::IO::black());
+    else
+      put(vcm, v, CGAL::IO::blue());
+  }
+
+  for(auto e : edges(sm))
+    put(ecm, e, CGAL::IO::gray());
+
+  CGAL_USE(fcm);
+
+  // Draw!
   CGAL::draw(sm);
 
   return EXIT_SUCCESS;
diff --git a/Surface_mesh/include/CGAL/draw_surface_mesh.h b/Surface_mesh/include/CGAL/draw_surface_mesh.h
index 3459392aec1..8dd25acd5b2 100644
--- a/Surface_mesh/include/CGAL/draw_surface_mesh.h
+++ b/Surface_mesh/include/CGAL/draw_surface_mesh.h
@@ -36,8 +36,58 @@ void draw(const SM& asm);
 #include <CGAL/Qt/init_ogl_context.h>
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/draw_face_graph.h>
-namespace CGAL
+#include <CGAL/use.h>
+
+namespace CGAL {
+
+// Check if there are any color maps that could be used, random otherwise
+template <typename K>
+struct Surface_mesh_basic_viewer_color_map
+  : DefaultColorFunctorFaceGraph
 {
+  using Base = DefaultColorFunctorFaceGraph;
+
+  using SM = ::CGAL::Surface_mesh<K>;
+  using vertex_descriptor = typename boost::graph_traits<SM>::vertex_descriptor;
+  using edge_descriptor = typename boost::graph_traits<SM>::edge_descriptor;
+  using face_descriptor = typename boost::graph_traits<SM>::face_descriptor;
+
+  using vertex_colors = typename SM::template Property_map<vertex_descriptor, CGAL::IO::Color>;
+  using edge_colors = typename SM::template Property_map<edge_descriptor, CGAL::IO::Color>;
+  using face_colors = typename SM::template Property_map<face_descriptor, CGAL::IO::Color>;
+
+  Surface_mesh_basic_viewer_color_map(const SM& amesh)
+  {
+    bool found = false;
+    std::tie(vcolors, found) = amesh.template property_map<vertex_descriptor, CGAL::IO::Color>("v:color");
+    std::tie(ecolors, found) = amesh.template property_map<edge_descriptor, CGAL::IO::Color>("e:color");
+    std::tie(fcolors, found) = amesh.template property_map<face_descriptor, CGAL::IO::Color>("f:color");
+    CGAL_USE(found);
+  }
+
+  CGAL::IO::Color operator()(const Surface_mesh<K>& amesh,
+                             const vertex_descriptor v) const
+  {
+    return vcolors ? get(vcolors, v) : Base::operator()(amesh, v);
+  }
+
+  CGAL::IO::Color operator()(const Surface_mesh<K>& amesh,
+                             const edge_descriptor e) const
+  {
+    return ecolors ? get(ecolors, e) : Base::operator()(amesh, e);
+  }
+
+  CGAL::IO::Color operator()(const Surface_mesh<K>& amesh,
+                             const face_descriptor f) const
+  {
+    return fcolors ? get(fcolors, f) : Base::operator()(amesh, f);
+  }
+
+private:
+  vertex_colors vcolors;
+  edge_colors ecolors;
+  face_colors fcolors;
+};
 
 // Specialization of draw function.
 template<class K>
@@ -57,8 +107,8 @@ void draw(const Surface_mesh<K>& amesh,
     int argc=1;
     const char* argv[2]={"surface_mesh_viewer", nullptr};
     QApplication app(argc,const_cast<char**>(argv));
-    SimpleFaceGraphViewerQt mainwindow(app.activeWindow(), amesh, title,
-                                       nofill);
+    SimpleFaceGraphViewerQt mainwindow(app.activeWindow(), amesh, title, nofill,
+                                       Surface_mesh_basic_viewer_color_map<K>(amesh));
     mainwindow.show();
     app.exec();
   }