Merge remote-tracking branch 'cgal/master' into CGAL-Prepare_CHANGES.md_6.1-GF

.github/workflows/build_doc.yml

@@ -29,6 +29,7 @@ permissions:
 
 jobs:
   pre_build_checks:
+    if: ${{ inputs.pr_number || startsWith(github.event.comment.body, '/build:') || startsWith(github.event.comment.body, '/force-build:') }}
     runs-on: ubuntu-latest
     name: Trigger the build?
     outputs:

.gitignore

@@ -1148,3 +1148,27 @@ Polygonal_surface_reconstruction/examples/build*
 Polygonal_surface_reconstruction/test/build*
 Solver_interface/examples/build*
 /Mesh_3/examples/Mesh_3/indicator_0.inr.gz
+/*.off
+/*.xyz
+/r0*
+all_segments.polylines.txt
+Data/data/meshes/*.*.edge
+Data/data/meshes/*.*.ele
+Data/data/meshes/*.*.face
+Data/data/meshes/*.*.mesh
+Data/data/meshes/*.*.node
+Data/data/meshes/*.*.smesh
+Data/data/meshes/*.*.vtk
+Data/data/meshes/*.log
+Data/data/meshes/*.off-cdt-output.off
+dump_*.off
+dump_*.txt
+dump-*.binary.cgal
+dump-*.polylines.txt
+dump-*.xyz
+dump.off.mesh
+log.txt
+patches_after_merge.ply
+CMakeUserPresets.json
+/.cache
+compile_commands.json

AABB_tree/doc/AABB_tree/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgAABBTreeRef AABB Tree Reference
+/// \defgroup PkgAABBTreeRef Reference Manual
 
 /// \defgroup PkgAABBTreeConcepts Concepts
 /// \ingroup PkgAABBTreeRef

AABB_tree/include/CGAL/AABB_traits_2.h

@@ -395,15 +395,15 @@ public:
   public:
       CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_true) const
       {
-          return GeomTraits().do_intersect_2_object()
+          return do_intersect_circle_iso_rectangle_2
           (GeomTraits().construct_circle_2_object()
-           (p, GeomTraits().compute_squared_distance_2_object()(p, bound)), bb,true)?
+           (p, GeomTraits().compute_squared_distance_2_object()(p, bound)), bb)?
           CGAL::SMALLER : CGAL::LARGER;
       }
 
       CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_false) const
       {
-          return GeomTraits().do_intersect_2_object()
+          return do_intersect_circle_iso_rectangle_2
           (GeomTraits().construct_circle_2_object()
            (p, GeomTraits().compute_squared_distance_2_object()(p, bound)), bb)?
           CGAL::SMALLER : CGAL::LARGER;
@@ -433,6 +433,45 @@ public:
           CGAL::SMALLER :
           CGAL::LARGER;
       }
+
+      typename GeomTraits::Boolean do_intersect_circle_iso_rectangle_2(const typename GeomTraits::Circle_2& circle,
+        const typename GeomTraits::Iso_rectangle_2& rec) const
+      {
+        typedef typename GeomTraits::FT       FT;
+        typedef typename GeomTraits::Point_2  Point;
+
+        Point center = circle.center();
+
+        // Check that the minimum distance to the box is smaller than the radius, otherwise there is
+        // no intersection. `distance` stays at 0 if the center is inside or on `rec`.
+        FT distance = FT(0);
+        if (center.x() < rec.xmin())
+        {
+          FT d = rec.xmin() - center.x();
+          distance += d * d;
+        }
+        else if (center.x() > rec.xmax())
+        {
+          FT d = center.x() - rec.xmax();
+          distance += d * d;
+        }
+
+        if (center.y() < rec.ymin())
+        {
+          FT d = rec.ymin() - center.y();
+          distance += d * d;
+        }
+        else if (center.y() > rec.ymax())
+        {
+          FT d = center.y() - rec.ymax();
+          distance += d * d;
+        }
+
+        if (distance <= circle.squared_radius())
+          return true;
+
+        return false;
+      }
   };
 
   Closest_point closest_point_object() const {return Closest_point(*this);}

AABB_tree/include/CGAL/AABB_traits_3.h

@@ -302,12 +302,11 @@ public:
     typename AT::Bounding_box operator()(ConstPrimitiveIterator first,
                                          ConstPrimitiveIterator beyond) const
     {
-      typename AT::Bounding_box bbox = m_traits.compute_bbox(*first,m_traits.bbm);
-      for(++first; first != beyond; ++first)
-        {
-          bbox = bbox + m_traits.compute_bbox(*first,m_traits.bbm);
-        }
-      return bbox;
+      return std::accumulate(first, beyond,
+        typename AT::Bounding_box{} /* empty bbox */,
+        [this](const typename AT::Bounding_box& bbox, const Primitive& pr) {
+          return bbox + m_traits.compute_bbox(pr, m_traits.bbm);
+        });
     }
 
   };

AABB_tree/include/CGAL/AABB_tree.h

@@ -199,9 +199,7 @@ namespace CGAL {
     }
 
     /// returns the axis-aligned bounding box of the whole tree.
-    /// \pre `!empty()`
     const Bounding_box bbox() const {
-      CGAL_precondition(!empty());
       if(size() > 1)
         return root_node()->bbox();
       else

AABB_tree/test/AABB_tree/aabb_test_closest_point_2.cpp

@@ -0,0 +1,80 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_rational.h>
+
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits_2.h>
+#include <CGAL/AABB_triangle_primitive_2.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+#include <array>
+#include <iostream>
+#include <string>
+#include <vector>
+
+template<typename Kernel>
+void test(const std::vector<CGAL::Simple_cartesian<double>::Point_2> &points, const std::vector<std::array<std::size_t, 3> > &faces) {
+  using Point_2 = typename Kernel::Point_2;
+  using Triangle_2 = typename Kernel::Triangle_2;
+  using Iterator = typename std::vector<Triangle_2>::const_iterator;
+  using Primitive = CGAL::AABB_triangle_primitive_2<Kernel, Iterator>;
+  using Tree_traits = CGAL::AABB_traits_2<Kernel, Primitive>;
+  using Tree = CGAL::AABB_tree<Tree_traits>;
+
+  std::vector<Triangle_2> triangles(faces.size());
+  for (std::size_t i = 0; i < faces.size(); ++i) {
+    const auto& f = faces[i];
+    triangles[i] = Triangle_2(Point_2(points[f[0]].x(), points[f[0]].y()), Point_2(points[f[1]].x(), points[f[1]].y()), Point_2(points[f[2]].x(), points[f[2]].y()));
+  }
+
+  Tree tree(triangles.begin(), triangles.end());
+
+  // Without hint
+  Point_2 query(-0.092372499264859229, -0.5067061545706153);
+  Point_2 closest_point = tree.closest_point(query);
+  std::cout << "Closest point to " << query << " is " << closest_point << std::endl;
+
+  // With hint
+  Point_2 hint(-0.077185400000000001, -0.42269299999999999);
+  Point_2 closest_point_hint = tree.closest_point(query, hint);
+  std::cout << "Closest point to " << query << " with hint " << hint << " is " << closest_point_hint << std::endl << std::endl;
+
+  assert(closest_point == closest_point_hint);
+}
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+
+  // Read the input
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/camel.off");
+  std::cout << "Reading " << filename << "..." << std::endl;
+
+  std::vector<CGAL::Simple_cartesian<double>::Point_3> points;
+  std::vector<std::array<std::size_t, 3> > faces;
+  if (!CGAL::IO::read_polygon_soup(filename, points, faces) || faces.empty())
+  {
+    std::cerr << "Invalid input:" << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "Input: " << points.size() << " points, " << faces.size() << " faces" << std::endl;
+
+  // Project onto the XY plane
+  std::vector<CGAL::Simple_cartesian<double>::Point_2> points_2(points.size());
+  for (std::size_t i = 0; i < points.size(); ++i)
+    points_2[i] = CGAL::Simple_cartesian<double>::Point_2(points[i].x(), points[i].y());
+
+  std::cout << "Testing closest point with Simple_cartesian<double>:" << std::endl;
+  test<CGAL::Simple_cartesian<double> >(points_2, faces);
+  std::cout << "Testing closest point with Epick:" << std::endl;
+  test<CGAL::Exact_predicates_inexact_constructions_kernel>(points_2, faces);
+  std::cout << "Testing closest point with Epeck:" << std::endl;
+  test<CGAL::Exact_predicates_exact_constructions_kernel>(points_2, faces);
+  std::cout << "Testing closest point with Simple_cartesian<Exact_rational>:" << std::endl;
+  test<CGAL::Simple_cartesian<CGAL::Exact_rational>>(points_2, faces);
+
+  std::cout << "Done." << std::endl;
+
+  return EXIT_SUCCESS;
+}

Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/Advancing_front_surface_reconstruction.txt

@@ -49,7 +49,7 @@ We describe next the algorithm and provide examples.
 
 \note A \ref tuto_reconstruction "detailed tutorial on surface reconstruction"
 is provided with a guide to choose the most appropriate method along
-with pre- and post-processing.
+with pre- and postprocessing.
 
 \section AFSR_Definitions Definitions and the Algorithm
 

Advancing_front_surface_reconstruction/doc/Advancing_front_surface_reconstruction/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgAdvancingFrontSurfaceReconstructionRef Advancing Front Surface Reconstruction Reference
+/// \defgroup PkgAdvancingFrontSurfaceReconstructionRef Reference Manual
 
 /// \defgroup PkgAdvancingFrontSurfaceReconstructionRefConcepts Concepts
 /// \ingroup PkgAdvancingFrontSurfaceReconstructionRef

Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h

@@ -382,7 +382,7 @@ namespace CGAL {
     int _facet_number;
 
     //---------------------------------------------------------------------
-    // For post-processing
+    // For postprocessing
     mutable int _postprocessing_counter;
     int _size_before_postprocessing;
 
@@ -2432,7 +2432,7 @@ namespace CGAL {
 
       std::size_t itmp, L_v_size_mem;
       L_v_size_mem = L_v.size();
-      if ((vh_on_border_inserted != 0)&& // to post-process only the borders
+      if ((vh_on_border_inserted != 0)&& // to postprocess only the borders
           (L_v.size() < .1 * _size_before_postprocessing))
         {
           {

Algebraic_foundations/doc/Algebraic_foundations/Algebraic_foundations.txt

@@ -80,12 +80,12 @@ using overloaded functions. However, for ease of use and backward
 compatibility all functionality is also
 accessible through global functions defined within namespace `CGAL`,
 e.g., \link sqrt `CGAL::sqrt(x)` \endlink. This is realized via function templates using
-the according functor of the traits class. For an overview see
-Section \ref PkgAlgebraicFoundationsRef in the reference manual.
+the according functor of the traits class. For an overview see the section "Global Functions" in the
+\ref PkgAlgebraicFoundationsRef.
 
 \subsection Algebraic_foundationsTagsinAlgebraicStructure Tags in Algebraic Structure Traits
 
-\subsection Algebraic_foundationsAlgebraicCategory Algebraic Category
+\subsubsection Algebraic_foundationsAlgebraicCategory Algebraic Category
 
 For a type `AS`, `Algebraic_structure_traits<AS>`
 provides several tags. The most important tag is the `Algebraic_category`
@@ -100,7 +100,7 @@ The tags are derived from each other such that they reflect the
 hierarchy of the algebraic structure concept, e.g.,
 `Field_with_sqrt_tag` is derived from `Field_tag`.
 
-\subsection Algebraic_foundationsExactandNumericalSensitive Exact and Numerical Sensitive
+\subsubsection Algebraic_foundationsExactandNumericalSensitive Exact and Numerical Sensitive
 
 Moreover, `Algebraic_structure_traits<AS>` provides the tags `Is_exact`
 and `Is_numerical_sensitive`, which are both `Boolean_tag`s.

Algebraic_foundations/doc/Algebraic_foundations/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgAlgebraicFoundationsRef Algebraic Foundations Reference
+/// \defgroup PkgAlgebraicFoundationsRef Reference Manual
 
 /// \defgroup PkgAlgebraicFoundationsAlgebraicStructuresConcepts Concepts
 /// \ingroup PkgAlgebraicFoundationsRef

Algebraic_kernel_d/doc/Algebraic_kernel_d/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgAlgebraicKernelDRef Algebraic Kernel Reference
+/// \defgroup PkgAlgebraicKernelDRef Reference Manual
 
 /// \defgroup PkgAlgebraicKernelDConcepts Concepts
 /// \ingroup PkgAlgebraicKernelDRef

Alpha_shapes_2/doc/Alpha_shapes_2/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgAlphaShapes2Ref 2D Alpha Shapes Reference
+/// \defgroup PkgAlphaShapes2Ref Reference Manual
 /// \defgroup PkgAlphaShapes2Concepts Concepts
 /// \ingroup PkgAlphaShapes2Ref
 

Alpha_shapes_3/doc/Alpha_shapes_3/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgAlphaShapes3Ref 3D Alpha Shapes Reference
+/// \defgroup PkgAlphaShapes3Ref Reference Manual
 /// \defgroup PkgAlphaShapes3Concepts Concepts
 /// \ingroup PkgAlphaShapes3Ref
 /*!

Alpha_wrap_3/doc/Alpha_wrap_3/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgAlphaWrap3Ref 3D Alpha Wrapping
+/// \defgroup PkgAlphaWrap3Ref Reference Manual
 
 /// \defgroup AW3_free_functions_grp Free Functions
 /// Functions to create a wrap from point clouds, triangle soups, and triangle meshes.

Alpha_wrap_3/examples/Alpha_wrap_3/volumetric_wrap.cpp

@@ -159,7 +159,9 @@ int main(int argc, char** argv)
   // edge length of regular tetrahedron with circumradius alpha
   const double l = 1.6329931618554521 * alpha; // sqrt(8/3)
 
-  CGAL::tetrahedral_isotropic_remeshing(tr, l, CGAL::parameters::remesh_boundaries(false));
+  CGAL::tetrahedral_isotropic_remeshing(tr, l,
+                                        CGAL::parameters::remesh_boundaries(false)
+                                                         .number_of_iterations(5));
 
   std::cout << "AFTER: " << tr.number_of_vertices() << " vertices, " << tr.number_of_cells() << " cells" << std::endl;
 

Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h

@@ -348,7 +348,7 @@ public:
 
 #ifdef CGAL_AW3_TIMER
       t.stop();
-      std::cout << "Manifoldness post-processing took: " << t.time() << " s." << std::endl;
+      std::cout << "Manifoldness postprocessing took: " << t.time() << " s." << std::endl;
       t.reset();
       t.start();
 #endif

Apollonius_graph_2/doc/Apollonius_graph_2/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgApolloniusGraph2Ref 2D Apollonius Graphs (Delaunay Graphs of Disks) Reference
+/// \defgroup PkgApolloniusGraph2Ref Reference Manual
 /// \defgroup PkgApolloniusGraph2Concepts Concepts
 /// \ingroup PkgApolloniusGraph2Ref
 /*!

Arrangement_on_surface_2/demo/Arrangement_on_surface_2_earth/Kml_reader.cpp

@@ -269,7 +269,7 @@ Kml::Nodes Kml::generate_ids_approx(Placemarks& placemarks, const double eps) {
 
         for (const auto& node : lring->nodes) {
           // check if there is a node sufficiently close to the current one
-          auto node_index = std::numeric_limits<std::size_t>::max();
+          auto node_index = (std::numeric_limits<std::size_t>::max)();
           for (std::size_t i = 0; i < nodes.size(); ++i) {
             const auto dist = node.distance_to(nodes[i]);
             if (dist < eps) {
@@ -278,7 +278,7 @@ Kml::Nodes Kml::generate_ids_approx(Placemarks& placemarks, const double eps) {
             }
           }
 
-          if (node_index == std::numeric_limits<std::size_t>::max()) {
+          if (node_index == (std::numeric_limits<std::size_t>::max)()) {
             // insert new node
             nodes.push_back(node);
             const auto node_id = nodes.size() - 1;
@@ -301,7 +301,7 @@ Kml::Nodes Kml::generate_ids_approx(Placemarks& placemarks, const double eps) {
   }
 
   // find the pair of closest nodes
-  double min_dist = std::numeric_limits<double>::max();
+  double min_dist = (std::numeric_limits<double>::max)();
   std::size_t ni1 = 0;
   std::size_t ni2 = 0;
   std::size_t num_nodes = nodes.size();

Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt

@@ -4196,7 +4196,7 @@ are both parameterized by a geometric kernel and model the concepts
 `AosTraits_2` and `AosLandmarkTraits_2`.
 \cgalFootnote{They also model the refined concept
 \cgalFootnoteCode{AosDirectionalXMonotoneTraits_2}, which enables Boolean set
-operations; see Package \ref PkgBooleanSetOperations2Ref.} The class
+operations; see Package \ref PkgBooleanSetOperations2.} The class
 template `Arr_non_caching_segment_traits_2<Kernel>` derives from the
 instance `Arr_non_caching_segment_basic_traits_2<Kernel>`, which
 models the `AosLandmarkTraits_2` traits concept but not the
@@ -4937,7 +4937,7 @@ template models the concepts `AosTraits_2` and
 `AosOpenBoundaryTraits_2`, but it does not model the
 `AosLandmarkTraits_2` concept. It also models the refined
 concept `AosDirectionalXMonotoneTraits_2`, which enables
-Boolean set operations; see Package \ref PkgBooleanSetOperations2Ref.
+Boolean set operations; see Package \ref PkgBooleanSetOperations2.
 Note that it is not a model of `AosLandmarkTraits_2` concept,
 so it is impossible to use the landmark point-location strategy with
 this traits class.
@@ -5167,7 +5167,7 @@ Every instance of the `Arr_Bezier_curve_traits_2` class templates
 models the concept `AosTraits_2` (but it does not model the
 `AosLandmarkTraits_2` concept). It also models the refined
 concept `AosDirectionalXMonotoneTraits_2`, which enables
-Boolean set operations; see Package \ref PkgBooleanSetOperations2Ref.
+Boolean set operations; see Package \ref PkgBooleanSetOperations2.
 
 <!-- ----------------------------------------------------------------------- -->
 \cgalFigureBegin{aos_fig-bezier_curves,bezier_curves.png}

Arrangement_on_surface_2/doc/Arrangement_on_surface_2/CGAL/Arr_observer.h

@@ -3,7 +3,7 @@ namespace CGAL {
 /*! \ingroup PkgArrangementOnSurface2Ref
  *
  * `Arr_observer<Arrangement_2>` is an alias for
- * Aos_observer<Arrangement_on_surface_2>`,
+ * `Aos_observer<Arrangement_on_surface_2>`,
  * where `Arrangement_2` derives from `Arrangement_on_surface_2` and the latter
  * is an instance of the template
  * `CGAL::Arrangement_on_surface_2<GeometryTraits, TopologyTraits>`.

Arrangement_on_surface_2/doc/Arrangement_on_surface_2/CGAL/draw_arrangement_2.h

@@ -31,7 +31,7 @@ namespace CGAL {
  *
  * A call to this function blocks the execution of the program until the drawing
  * window is closed. This function requires `CGAL_Qt6`, and is only available if
- * the macro `CGAL_USE_BASIC_VIEWER` is defined.  Linking with the cmake target
+ * the macro `CGAL_USE_BASIC_VIEWER` is defined. Linking with the cmake target
  * `CGAL::CGAL_Basic_viewer` will link with `CGAL_Qt6` and add the definition
  * `CGAL_USE_BASIC_VIEWER`.
  *

Arrangement_on_surface_2/doc/Arrangement_on_surface_2/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgArrangementOnSurface2Ref 2D Arrangement Reference
+/// \defgroup PkgArrangementOnSurface2Ref Reference Manual
 
 /// \defgroup PkgArrangementOnSurface2Concepts Concepts
 /// \ingroup PkgArrangementOnSurface2Ref

Arrangement_on_surface_2/include/CGAL/draw_arrangement_2.h

@@ -25,7 +25,7 @@
 #include <cstdlib>
 #include <random>
 
-#include <CGAL/Qt/Basic_viewer.h>
+#include <CGAL/Basic_viewer.h>
 #include <CGAL/Graphics_scene.h>
 #include <CGAL/Graphics_scene_options.h>
 #include <CGAL/Random.h>
@@ -587,8 +587,6 @@ void add_to_graphics_scene(const CGAL_ARR_TYPE& aos,
   add_to_graphics_scene(aos, graphics_scene, gso);
 }
 
-#ifdef CGAL_USE_BASIC_VIEWER
-
 /// Draw an arrangement on surface.
 template <typename GeometryTraits_2, typename TopologyTraits, class GSOptions>
 void draw(const CGAL_ARR_TYPE& aos, const GSOptions& gso,
@@ -609,8 +607,6 @@ void draw(const CGAL_ARR_TYPE& aos,
   draw_graphics_scene(graphics_scene, title);
 }
 
-#endif // CGAL_USE_BASIC_VIEWER
-
 #undef CGAL_ARR_TYPE
 
 } // namespace CGAL

BGL/doc/BGL/BGL.txt

@@ -698,7 +698,7 @@ using one of the following tags:
 - `CGAL::Alpha_expansion_boost_adjacency_list_tag` (default)
 - `CGAL::Alpha_expansion_boost_compressed_sparse_raw_tag`
 - `CGAL::Alpha_expansion_MaxFlow_tag`, released under GPL
-  license and provided by the \ref PkgSurfaceMeshSegmentationRef
+  license and provided by the \ref PkgSurfaceMeshSegmentation
   package
 
 All these implementations produce the exact same result but behave

BGL/doc/BGL/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgBGLRef CGAL and the Boost Graph Library Reference
+/// \defgroup PkgBGLRef Reference Manual
 
 /*! \defgroup PkgBGLConcepts Concepts
  \ingroup PkgBGLRef
@@ -31,7 +31,7 @@ digraph example {
 }
 \enddot
 
- \cgalHeading{Notations}
+\cgalHeading{Notations}
 
  <dl>
  <dt>`G`</dt>         <dd>A type that is a model of a graph concept.</dd>
@@ -51,8 +51,8 @@ and adds the requirement for traversal of all vertices in a graph.
 <center>
 <table>
  <tr>
-  <th width="800">Associated Type</th>
-  <th width="800">Description</th>
+  <th class="associated-type">Associated Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`boost::graph_traits<G>::%vertex_iterator`</td>
@@ -66,9 +66,9 @@ and adds the requirement for traversal of all vertices in a graph.
 
 <table>
  <tr>
-  <th width="392">Valid Expression</th>
-  <th width="392">Return Type</th>
-  <th width="800">Description</th>
+  <th class="valid-expression">Valid Expression</th>
+  <th class="return-type">Return Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`vertices(g)`</td>
@@ -92,8 +92,8 @@ and adds the requirement for traversal of all edges in a graph.
 <center>
 <table>
  <tr>
-  <th width="800">Associated Type</th>
-  <th width="800">Description</th>
+  <th class="associated-type">Associated Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`boost::graph_traits<G>::%edge_iterator`</td>
@@ -107,9 +107,9 @@ and adds the requirement for traversal of all edges in a graph.
 
 <table>
  <tr>
-  <th width="392">Valid Expression</th>
-  <th width="392">Return Type</th>
-  <th width="800">Description</th>
+  <th class="valid-expression">Valid Expression</th>
+  <th class="return-type">Return Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`edges(g)`</td>
@@ -143,8 +143,8 @@ and adds the notion of halfedges, where each edge corresponds to two opposite ha
 <center>
 <table>
  <tr>
-  <th width="800">Associated Type</th>
-  <th width="800">Description</th>
+  <th class="associated-type">Associated Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`boost::graph_traits<G>::%halfedge_descriptor`</td>
@@ -154,9 +154,9 @@ and adds the notion of halfedges, where each edge corresponds to two opposite ha
 
 <table>
  <tr>
-  <th width="392">Valid Expression</th>
-  <th width="392">Return Type</th>
-  <th width="800">Description</th>
+  <th class="valid-expression">Valid Expression</th>
+  <th class="return-type">Return Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`edge(h, g)`</td>
@@ -222,9 +222,9 @@ update the incidence information between vertices and halfedges.
 <center>
 <table>
  <tr>
-  <th width="392">Valid Expression</th>
-  <th width="392">Return Type</th>
-  <th width="800">Description</th>
+  <th class="valid-expression">Valid Expression</th>
+  <th class="return-type">Return Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`add_vertex(g)`</td>
@@ -272,8 +272,8 @@ and adds the requirements for traversal of all halfedges in the graph.
 <center>
 <table>
  <tr>
-  <th width="800">Associated Type</th>
-  <th width="800">Description</th>
+  <th class="associated-type">Associated Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`boost::graph_traits<G>::%halfedge_iterator`</td>
@@ -287,9 +287,9 @@ and adds the requirements for traversal of all halfedges in the graph.
 
 <table>
  <tr>
-  <th width="392">Valid Expression</th>
-  <th width="392">Return Type</th>
-  <th width="800">Description</th>
+  <th class="valid-expression">Valid Expression</th>
+  <th class="return-type">Return Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`halfedges(g)`</td>
@@ -315,8 +315,8 @@ face.
 <center>
 <table>
  <tr>
-  <th width="800">Associated Type</th>
-  <th width="800">Description</th>
+  <th class="associated-type">Associated Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`boost::graph_traits<G>::%face_descriptor`</td>
@@ -326,9 +326,9 @@ face.
 
 <table>
  <tr>
-  <th width="392">Valid Expression</th>
-  <th width="392">Return Type</th>
-  <th width="800">Description</th>
+  <th class="valid-expression">Valid Expression</th>
+  <th class="return-type">Return Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`face(h, g)`</td>
@@ -361,9 +361,9 @@ the requirement for operations to add faces and to modify face-halfedge relation
 <center>
 <table>
  <tr>
-  <th width="392">Valid Expression</th>
-  <th width="392">Return Type</th>
-  <th width="800">Description</th>
+  <th class="valid-expression">Valid Expression</th>
+  <th class="return-type">Return Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`add_face(g)`</td>
@@ -401,8 +401,8 @@ the requirement for traversal of all faces in a graph.
 <center>
 <table>
  <tr>
-  <th width="800">Associated Type</th>
-  <th width="800">Description</th>
+  <th class="associated-type">Associated Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`boost::graph_traits<G>::%face_iterator`</td>
@@ -416,9 +416,9 @@ the requirement for traversal of all faces in a graph.
 
 <table>
  <tr>
-  <th width="392">Valid Expression</th>
-  <th width="392">Return Type</th>
-  <th width="800">Description</th>
+  <th class="valid-expression">Valid Expression</th>
+  <th class="return-type">Return Type</th>
+  <th class="description">Description</th>
  </tr>
  <tr>
   <td>`faces(g)`</td>

BGL/include/CGAL/boost/graph/alpha_expansion_graphcut.h

@@ -509,7 +509,8 @@ class Alpha_expansion_MaxFlow_impl;
      \cgalParamNEnd
    \cgalNamedParamsEnd
 
-   \note The `MaxFlow` implementation is provided by the \ref PkgSurfaceMeshSegmentationRef
+
+   \note The `MaxFlow` implementation is provided by the \ref PkgSurfaceMeshSegmentation package
    under GPL license. The header `<CGAL/boost/graph/Alpha_expansion_MaxFlow_tag.h>`
    must be included if users want to use this implementation.
 */

BGL/include/CGAL/boost/graph/named_params_helper.h

@@ -227,6 +227,32 @@ struct GetGeomTraits
                                           NamedParametersVPM>::type type;
 };
 
+namespace internal {
+// Similar helper for polygon soups
+template <typename PointRange, typename PolygonRange>
+struct Polygon_types
+{
+  typedef typename boost::range_value<PointRange>::type                             Point_3;
+  typedef typename boost::range_value<PolygonRange>::type                           Polygon_3;
+
+  typedef typename boost::range_iterator<Polygon_3>::type                           V_ID_iterator;
+  typedef typename std::iterator_traits<V_ID_iterator>::value_type                  V_ID;
+  typedef typename std::vector<Polygon_3>::size_type                                P_ID;
+};
+}
+
+template <typename PointRange, typename PolygonRange, typename NamedParameters>
+struct GetPolygonGeomTraits
+{
+  typedef typename internal_np::Lookup_named_param_def <
+                     internal_np::geom_traits_t,
+                     NamedParameters,
+                     typename CGAL::Kernel_traits<
+                       typename internal::Polygon_types<
+                         PointRange, PolygonRange>::Point_3 >::type
+                   > ::type                                                         type;
+};
+
 // Define the following structs:
 //
 // GetInitializedVertexIndexMap

BGL/include/CGAL/boost/graph/split_graph_into_polylines.h

@@ -16,7 +16,6 @@
 
 #include <map>
 #include <vector>
-#include <utility>
 #include <boost/graph/adjacency_list.hpp>
 #include <CGAL/assertions.h>
 #include <CGAL/tags.h>
@@ -143,34 +142,29 @@ void duplicate_terminal_vertices(Graph& graph,
 {
   typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
   typedef typename boost::graph_traits<Graph>::edge_descriptor edge_descriptor;
-  typedef typename boost::graph_traits<Graph>::vertex_iterator vertex_iterator;
-  typedef typename boost::graph_traits<Graph>::out_edge_iterator out_edge_iterator;
 
-  vertex_iterator b,e;
-  std::tie(b,e) = vertices(graph);
-  std::vector<vertex_descriptor> V(b,e);
+  auto [b, e] = vertices(graph);
+  std::vector<vertex_descriptor> V(b, e); // copy vertices, because the graph may change
+  std::vector<edge_descriptor> out_edges_of_v; // used to store the out edges of a vertex
+                                               // created here to avoid allocating it in the loop
   for(vertex_descriptor v : V)
   {
-    typename boost::graph_traits<OrigGraph>::vertex_descriptor orig_v = graph[v];
-    typename boost::graph_traits<Graph>::degree_size_type deg = degree(v, graph);
+    auto orig_v = graph[v];
+    auto deg = degree(v, graph);
     if (deg != 2 || is_terminal(orig_v, orig))
+    {
+      auto [b, e] = out_edges(v, graph);
+      out_edges_of_v.assign(b, e); // same as creating a new vector from the range [b,e)
+      for (unsigned int i = 1; i < out_edges_of_v.size(); ++i)
       {
-        out_edge_iterator b, e;
-        std::tie(b, e) = out_edges(v, graph);
-        std::vector<edge_descriptor> out_edges_of_v(b, e);
-        for (unsigned int i = 1; i < out_edges_of_v.size(); ++i)
-          {
-            edge_descriptor e = out_edges_of_v[i];
-            typename boost::graph_traits<OrigGraph>::edge_descriptor orig_e =
-              graph[e];
-            vertex_descriptor w = target(e, graph);
-            remove_edge(e, graph);
-            vertex_descriptor vc = add_vertex(graph);
-            graph[vc] = orig_v;
-            const std::pair<edge_descriptor, bool> pair = add_edge(vc, w, graph);
-            graph[pair.first] = orig_e;
-          }
+        edge_descriptor e = out_edges_of_v[i];
+        auto orig_e = graph[e];
+        vertex_descriptor w = target(e, graph);
+        remove_edge(e, graph);
+        vertex_descriptor vc = add_vertex(orig_v, graph);
+        add_edge(vc, w, orig_e, graph);
       }
+    }
   }
 
   // check all vertices are of degree 1 or 2 and that the source
@@ -271,8 +265,7 @@ split_graph_into_polylines(const Graph& graph,
     V2vmap v2vmap;
 
     for(Graph_vertex_descriptor v : make_range(vertices(graph))){
-      vertex_descriptor vc = add_vertex(g_copy);
-      g_copy[vc] = v;
+      vertex_descriptor vc = add_vertex(v, g_copy);
       v2vmap[v] = vc;
     }
 
@@ -281,9 +274,7 @@ split_graph_into_polylines(const Graph& graph,
       Graph_vertex_descriptor vt = target(e,graph);
       CGAL_warning_msg(vs != vt, "ignore self loops");
       if(vs != vt){
-        const std::pair<edge_descriptor, bool> pair =
-          add_edge(v2vmap[vs],v2vmap[vt],g_copy);
-        g_copy[pair.first] = e;
+        add_edge(v2vmap[vs], v2vmap[vt], e, g_copy);
       }
     }
   }

BGL/test/BGL/CMakeLists.txt

@@ -63,6 +63,7 @@ else()
   message(STATUS "NOTICE: Tests that use OpenMesh will not be compiled.")
 endif()
 
+set(CMAKE_POLICY_DEFAULT_CMP0167 NEW)
 find_package(VTK 9.0 QUIET COMPONENTS CommonCore IOCore IOLegacy IOXML FiltersCore FiltersSources)
 if (VTK_FOUND AND VTK_LIBRARIES)
   message(STATUS "VTK ${VTK_VERSION} found ${VTK_LIBRARIES}")

Barycentric_coordinates_2/doc/Barycentric_coordinates_2/Barycentric_coordinates_2.txt

@@ -53,7 +53,7 @@ All analytic barycentric coordinates for polygons can be computed either by inst
 or through one of the free functions. Harmonic coordinates can be computed only by
 instantiating a class that must be parameterized by a model of the concept `DiscretizedDomain_2`.
 Segment and triangle coordinates can be computed only through the free functions.
-For more information see the \ref PkgBarycentricCoordinates2Ref "Reference Manual".
+For more information see the \ref PkgBarycentricCoordinates2Ref.
 
 Any point in the plane may be taken as a query point. However, we do not recommend using
 Wachspress and discrete harmonic coordinates with query points outside the closure

Barycentric_coordinates_2/doc/Barycentric_coordinates_2/PackageDescription.txt

@@ -1,8 +1,5 @@
-namespace CGAL {
-namespace Barycentric_coordinates {
-
 /*!
-\defgroup PkgBarycentricCoordinates2Ref 2D Generalized Barycentric Coordinates Reference
+\defgroup PkgBarycentricCoordinates2Ref Reference Manual
 
 \defgroup PkgBarycentricCoordinates2RefConcepts Concepts
 \ingroup PkgBarycentricCoordinates2Ref
@@ -77,6 +74,3 @@ coordinates from the Package \ref PkgInterpolation2.}
 - `discrete_harmonic_coordinates_2()`
 - `boundary_coordinates_2()`
 */
-
-} /* namespace Barycentric_coordinates */
-} /* namespace CGAL */

Basic_viewer/doc/Basic_viewer/Basic_viewer.txt

@@ -131,7 +131,7 @@ Example of drawing of a point cloud and a polyhedron in a same viewer.
 \cgalFigureEnd
 
 <!-- /////////////////////////////////////////////////////////////////////////////// -->
-\section BV_BasicViewer The Basic Viewer Class
+\section BV_BasicViewer The Qt Basic Viewer Class
 
 The class `CGAL::Qt::Basic_viewer` is a \qt widget that inherits from `QGLViewer` and mainly stores a `Graphics_scene` and allows to visualize it and interact with the scene. Since this class is a \qt widget, it can be used into more complex \qt code to create more advanced demos.
 

Basic_viewer/doc/Basic_viewer/CGAL/Graphics_scene.h

@@ -102,6 +102,36 @@ public:
 
   /// returns `true` if the scene is in 2D, i.e., lies on the XY or XZ or YZ plane.
   bool is_two_dimensional() const;
+
+  /// set the default color of faces
+  void set_default_color_face(const CGAL::IO::Color& c);
+
+  /// set the default color of points
+  void set_default_color_point(const CGAL::IO::Color& c);
+
+  /// set the default color of segments
+  void set_default_color_segment(const CGAL::IO::Color& c);
+
+  /// set the default color of rays
+  void set_default_color_ray(const CGAL::IO::Color& c);
+
+  /// set the default color of lines
+  void set_default_color_line(const CGAL::IO::Color& c);
+
+  /// returns the default color of faces
+  const CGAL::IO::Color &get_default_color_face() const;
+
+  /// returns the default color of points
+  const CGAL::IO::Color &get_default_color_point() const;
+
+  /// returns the default color of segments
+  const CGAL::IO::Color &get_default_color_segment() const;
+
+  /// returns the default color of rays
+  const CGAL::IO::Color &get_default_color_ray() const;
+
+  /// returns the default color of lines
+  const CGAL::IO::Color &get_default_color_line() const;
 };
 
 } // namespace CGAL

Basic_viewer/doc/Basic_viewer/CGAL/Qt/Basic_viewer.h

@@ -12,119 +12,180 @@ CGAL::QGLViewer is our internal fork of <a href="https://github.com/GillesDebunn
 class Basic_viewer : public CGAL::QGLViewer
 {
 public:
+  /// \name Constructors
+  /// @{
+
   /// Constructor given a pointer on a `QWidget` (can be a `nullptr`) and a `Graphics_scene`.
   /// `title` will be the title of the window.
   Basic_viewer(QWidget* parent,
                const Graphics_scene& scene,
                const char* title="");
 
-  /// enables or disables the drawing of vertices.
+  /// @}
+
+  /// \name Setters
+  /// @{
+
+  /// \brief Set size of vertices.
+  /// \param s The size of vertices.
+  void size_vertices(float s);
+
+  /// \brief Set size of edges.
+  /// \param s The size of edges.
+  void size_edges(float s);
+
+  /// \brief Set size of rays.
+  /// \param s The size of rays.
+  void size_rays(float s);
+
+  /// \brief Set size of lines.
+  /// \param s The size of lines.
+  void size_lines(float s);
+
+  /// \brief Enables or disables the drawing of vertices.
+  /// \param b Set to `true` to enable, `false` to disable.
   void draw_vertices(bool b);
 
-  /// enables or disables the drawing of edges.
+  /// \brief Enables or disables the drawing of edges.
+  /// \param b Set to `true` to enable, `false` to disable.
   void draw_edges(bool b);
 
-  /// enables or disables the drawing of rays.
+  /// \brief Enables or disables the drawing of rays.
+  /// \param b Set to `true` to enable, `false` to disable.
   void draw_rays(bool b);
 
-  /// enables or disables the drawing of lines.
+  /// \brief Enables or disables the drawing of lines.
+  /// \param b Set to `true` to enable, `false` to disable.
   void draw_lines(bool b);
 
-  /// enables or disables the drawing of faces.
+  /// \brief enables or disables the drawing of faces.
+  /// \param b Set to `true` to enable, `false` to disable.
   void draw_faces(bool b);
 
-  /// enables or disables the use of only one color (if `b` is `true`) or the use of multiple colors (if `b` is `false`).
-  void use_mono_color(bool b);
-
   /// enables or disables the drawing of texts.
+  /// \brief enables or disables the drawing of texts.
+  /// \param b Set to `true` to enable, `false` to disable.
   void draw_text(bool b);
 
-  /// sets the color used for vertices in mono color mode.
-  void vertices_mono_color(const CGAL::IO::Color& c);
+  /// \brief Enables or disables the drawing of mesh triangles.
+  /// \param b Set to `true` to enable, `false` to disable.
+  void draw_mesh_triangles(bool b);
 
-  /// sets the color used for edges in mono color mode.
-  void edges_mono_color(const CGAL::IO::Color& c);
+  /// \brief Enables or disables the use of only one color or the use of multiple colors.
+  /// \param b Set to `true` to use only one color, `false` to use multiple colors.
+  void use_default_color(bool b);
 
-  /// sets the color used for rays in mono color mode.
-  void rays_mono_color(const CGAL::IO::Color& c);
+  /// \brief Enables or disables the use of a single color for all normals.
+  /// \param b Set to `true` to enable, `false` to disable.
+  void use_default_color_normals(bool b);
 
-  /// sets the color used for lines in mono color mode.
-  void lines_mono_color(const CGAL::IO::Color& c);
+  /// \brief enables or disables the use of flat shading or the use of smooth shading.
+  /// \param b Set to `true` for flat shading, `false` for smooth shading.
+  void flat_shading(bool b);
 
-  /// sets the color used for faces in mono color mode.
-  void faces_mono_color(const CGAL::IO::Color& c);
+  /// \brief Enables or disables the reversal of normals.
+  /// \param b Set to `true` to reverse normals, `false` to keep normals as is.
+  void reverse_normal(bool b);
 
-  /// toggles the drawing of vertices.
+  /// \brief Sets the default color of the normals.
+  /// \param c The default color of the normals.
+  void default_color_normals(const CGAL::IO::Color& c);
+
+  /// \brief Sets the height factor value of the normals.
+  /// \param h The height factor value of the normals.
+  void normal_height_factor(float h);
+
+  /// \brief Toggles the drawing of vertices.
   void toggle_draw_vertices();
 
-  /// toggles the drawing of edges.
+  /// \brief Toggles the drawing of edges.
   void toggle_draw_edges();
 
-  /// toggles the drawing of rays.
+  /// \brief Toggles the drawing of rays.
   void toggle_draw_rays();
 
-  /// toggles the drawing of lines.
+  /// \brief Toggles the drawing of lines.
   void toggle_draw_lines();
 
-  /// toggles the drawing of faces.
+  /// \brief Toggles the drawing of faces.
   void toggle_draw_faces();
 
-  /// toggles the use of mono color mode.
-  void toggle_use_mono_color();
+  /// \brief Toggles the use of mono color mode.
+  void toggle_use_default_color();
 
-  /// toggles the drawing of text.
+  /// \brief Toggles the use of the default color mode for normals.
+  void toggle_use_default_color_normal();
+
+  /// \brief Toggles the use of flat shading.
+  void toggle_flat_shading();
+
+  /// \brief Toggles the drawing of text.
   void toggle_draw_text();
 
-  /// returns `true` if vertices are drawn.
+  /// \brief Reverses all normals of vertices and faces.
+  void reverse_all_normals();
+
+  /// @}
+
+  /// \name Getters
+  /// @{
+
+  /// \brief Checks if vertices are drawn.
+  /// \return `true` if vertices are drawn, `false` otherwise.
   bool draw_vertices() const;
 
-  /// returns `true` if edges are drawn.
+  /// \brief Checks if edges are drawn.
+  /// \return `true` if edges are drawn, `false` otherwise.
   bool draw_edges() const;
 
-  /// returns `true` if rays are drawn.
+  /// \brief Checks if rays are drawn.
+  /// \return `true` if rays are drawn, `false` otherwise.
   bool draw_rays() const;
 
-  /// returns `true` if lines are drawn.
+  /// \brief Checks if lines are drawn.
+  /// \return `true` if lines are drawn, `false` otherwise.
   bool draw_lines() const;
 
-  /// returns `true` if faces are drawn.
+  /// \brief Checks if faces are drawn.
+  /// \return `true` if faces are drawn, `false` otherwise.
   bool draw_faces() const;
 
-  /// returns `true` if mono color mode is used.
-  bool use_mono_color() const;
-
-  /// returns `true` if normals are reversed.
-  bool reverse_normal() const;
-
-  /// returns `true` if text are drawn.
+  /// \brief Checks if text is drawn.
+  /// \return `true` if text is drawn, `false` otherwise.
   bool draw_text() const;
 
-  /// returns the mono color used for vertices.
-  const CGAL::IO::Color& vertices_mono_color() const;
+  /// \brief Checks if the default color mode is used.
+  /// \return `true` if mono color mode is used, `false` otherwise.
+  bool use_default_color() const;
 
-  /// returns the mono color used for edges.
-  const CGAL::IO::Color& edges_mono_color() const;
+  /// \brief Checks if the default color mode for normals is used.
+  /// \return `true` if default color mode for normals is used, `false` otherwise.
+  bool use_default_color_normal() const;
 
-  /// returns the mono color used for rays.
-  const CGAL::IO::Color& rays_mono_color() const;
+  /// \brief Checks if normals are reversed.
+  /// \return `true` if normals are reversed, `false` otherwise.
+  bool reverse_normal() const;
 
-  /// returns the mono color used for lines.
-  const CGAL::IO::Color& lines_mono_color() const;
-
-  /// returns the mono color used for faces.
-  const CGAL::IO::Color& faces_mono_color() const;
-
-  /// returns `true` if the clipping plane is enabled.
+  /// \brief Checks if the clipping plane is enabled.
+  /// \return `true` if the clipping plane is enabled, `false` otherwise.
   bool clipping_plane_enabled() const;
 
-  /// returns the clipping plane when it is enabled.
+  /// \brief Checks if m_no_2D_mode is false and the graphics scene is two-dimensional.
+  /// \return `true` if m_no_2D_mode is false and the scene is 2D, `false` otherwise.
+  bool is_two_dimensional() const;
+
+  /// \brief Gets the clipping plane when enabled.
+  /// \return The clipping plane as a `CGAL::Exact_predicates_inexact_constructions_kernel::Plane_3` object.
   CGAL::Exact_predicates_inexact_constructions_kernel::Plane_3 clipping_plane() const;
 
-  /// returns the graphics scene of the viewer.
+  /// \brief Gets the graphics scene of the viewer.
+  /// \return A reference to the `Graphics_scene` object.
   const Graphics_scene& graphics_scene() const;
 
-  /// reverses all normals of vertices and faces.
-  void reverse_all_normals();
+  /// @}
+
+  /// \name Draw
+  /// @{
 
   /// draws the viewer without recomputing all internal buffers.
   virtual void draw();
@@ -132,6 +193,8 @@ public:
   /// redraws the viewer, i.e., recompute all internal buffers and update the window.
   virtual void redraw();
 
+  /// @}
+
   /// Function called when a key is pressed. Users can define their own function in order
   /// to add specific behavior.
   std::function<bool(QKeyEvent *, CGAL::Qt::Basic_viewer *)> on_key_pressed;

Basic_viewer/doc/Basic_viewer/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgBasicViewerRef Basic Viewer Reference
+/// \defgroup PkgBasicViewerRef Reference Manual
 
 /// \defgroup PkgBasicViewerConcepts Concepts
 /// \ingroup PkgBasicViewerRef

Basic_viewer/examples/Basic_viewer/draw_mesh_and_points.cpp

@@ -5,7 +5,7 @@
 #include <CGAL/draw_polyhedron.h>
 #include <CGAL/draw_point_set_3.h>
 #include <CGAL/Graphics_scene_options.h>
-#include <CGAL/Qt/Basic_viewer.h>
+#include <CGAL/Basic_viewer.h>
 
 #include <vector>
 #include <iostream>

Basic_viewer/examples/Basic_viewer/draw_several_windows.cpp

@@ -5,7 +5,7 @@
 #include <CGAL/draw_polyhedron.h>
 #include <CGAL/draw_point_set_3.h>
 #include <CGAL/Graphics_scene_options.h>
-#include <CGAL/Qt/Basic_viewer.h>
+#include <CGAL/Basic_viewer.h>
 
 #ifdef CGAL_USE_BASIC_VIEWER
 #include <QMainWindow>

Basic_viewer/include/CGAL/Basic_shaders.h

@@ -17,131 +17,299 @@
 
 namespace CGAL
 {
-
 //------------------------------------------------------------------------------
-const char vertex_source_color[]=R"DELIM(
+const char VERTEX_SOURCE_COLOR[]=R"DELIM(
 #version 150
-in highp vec4 vertex;
-in highp vec3 normal;
-in highp vec3 color;
+in highp   vec3 a_Pos;
+in highp   vec3 a_Normal;
+in mediump vec3 a_Color;
 
-uniform highp mat4 mvp_matrix;
-uniform highp mat4 mv_matrix;
-uniform highp float point_size;
+out highp   vec4 vs_fP; // view space position
+out highp   vec4 ls_fP; // local space position
+out highp   vec3 fN;
+out mediump vec4 fColor;
 
-out highp vec4 fP;
-out highp vec3 fN;
-out highp vec4 fColor;
-out highp vec4 m_vertex;
+uniform highp   mat4  u_Mvp;
+uniform highp   mat4  u_Mv;
+uniform mediump float u_PointSize;
+uniform mediump vec3  u_DefaultColor;
+uniform         bool  u_UseDefaultColor;
 
 void main(void)
 {
-  fP = mv_matrix * vertex;
-  fN = mat3(mv_matrix)* normal;
-  fColor = vec4(color, 1.0);
-  gl_PointSize = point_size;
+  fColor = vec4(a_Color, 1.0);
+  if (u_UseDefaultColor)
+  {
+    fColor = vec4(u_DefaultColor, 1.0);
+  }
 
-  m_vertex = vertex;
+  vec4 pos = vec4(a_Pos, 1.0);
 
-  gl_Position = mvp_matrix * vertex;
+  ls_fP = pos;
+  vs_fP = u_Mv * pos;
+
+  fN = mat3(u_Mv)* a_Normal;
+
+  gl_Position = u_Mvp * pos;
+  gl_PointSize = u_PointSize;
 }
 )DELIM";
 
-const char fragment_source_color[]=R"DELIM(
+const char FRAGMENT_SOURCE_COLOR[]=R"DELIM(
 #version 150
-in highp vec4 fP;
-in highp vec3 fN;
-in highp vec4 fColor;
-in highp vec4 m_vertex;
+in highp   vec4 vs_fP;
+in highp   vec4 ls_fP;
+in highp   vec3 fN;
+in mediump vec4 fColor;
 
-uniform highp vec4 light_pos;
-uniform highp vec4 light_diff;
-uniform highp vec4 light_spec;
-uniform highp vec4 light_amb;
-uniform highp float spec_power;
+out mediump vec4 out_color;
 
-uniform highp vec4 clipPlane;
-uniform highp vec4 pointPlane;
-uniform highp float rendering_mode;
-uniform highp float rendering_transparency;
+uniform highp   vec4  u_LightPos;
+uniform mediump vec4  u_LightDiff;
+uniform mediump vec4  u_LightSpec;
+uniform mediump vec4  u_LightAmb;
+uniform mediump float u_SpecPower;
 
-out highp vec4 out_color;
+uniform highp   vec4  u_ClipPlane;
+uniform highp   vec4  u_PointPlane;
+uniform mediump float u_RenderingMode;
+uniform mediump float u_RenderingTransparency;
 
 void main(void)
 {
-  highp vec3 L = light_pos.xyz - fP.xyz;
-  highp vec3 V = -fP.xyz;
+  highp vec3 L = u_LightPos.xyz - vs_fP.xyz;
+  highp vec3 V = -vs_fP.xyz;
 
-  highp vec3 N = normalize(fN);
+  highp vec3 a_Normal = normalize(fN);
   L = normalize(L);
   V = normalize(V);
 
-  highp vec3 R = reflect(-L, N);
-  highp vec4 diffuse = vec4(max(dot(N,L), 0.0) * light_diff.rgb * fColor.rgb, 1.0);
-  highp vec4 ambient = vec4(light_amb.rgb * fColor.rgb, 1.0);
-  highp vec4 specular = pow(max(dot(R,V), 0.0), spec_power) * light_spec;
+  highp vec3 R = reflect(-L, a_Normal);
+  highp vec4 diffuse = vec4(max(dot(a_Normal,L), 0.0) * u_LightDiff.rgb * fColor.rgb, 1.0);
+  highp vec4 ambient = vec4(u_LightAmb.rgb * fColor.rgb, 1.0);
+  highp vec4 specular = pow(max(dot(R,V), 0.0), u_SpecPower) * u_LightSpec;
 
   // onPlane == 1: inside clipping plane, should be solid;
   // onPlane == -1: outside clipping plane, should be transparent;
   // onPlane == 0: on clipping plane, whatever;
-  float onPlane = sign(dot((m_vertex.xyz-pointPlane.xyz), clipPlane.xyz));
+  float onPlane = sign(dot((ls_fP.xyz-u_PointPlane.xyz), u_ClipPlane.xyz));
 
   // rendering_mode == -1: draw all solid;
   // rendering_mode == 0: draw solid only;
   // rendering_mode == 1: draw transparent only;
-  if (rendering_mode == (onPlane+1)/2) {
+  if (u_RenderingMode == (onPlane+1)/2) {
     // discard other than the corresponding half when rendering
     discard;
   }
 
   // draw corresponding part
-  out_color = rendering_mode < 1 ? (diffuse + ambient) :
-                      vec4(diffuse.rgb + ambient.rgb, rendering_transparency);
+  out_color = u_RenderingMode < 1 ? (diffuse + ambient) :
+                      vec4(diffuse.rgb + ambient.rgb, u_RenderingTransparency);
 }
 )DELIM";
 
-const char vertex_source_p_l[]=R"DELIM(
+const char VERTEX_SOURCE_P_L[]=R"DELIM(
 #version 150
-in highp vec4 vertex;
-in highp vec3 color;
+in highp   vec3 a_Pos;
+in mediump vec3 a_Color;
 
-uniform highp mat4 mvp_matrix;
-uniform highp float point_size;
+out mediump vec4 fColor;
+out highp   vec4 ls_fP; // local space
 
-out highp vec4 fColor;
-out highp vec4 m_vertex;
+uniform highp   mat4  u_Mvp;
+uniform mediump float u_PointSize;
+uniform         bool  u_IsOrthographic;
+uniform mediump vec3  u_DefaultColor;
+uniform         bool  u_UseDefaultColor;
+
+bool EqualZero(float value)
+{
+  return abs(value) < 0.00001;
+}
 
 void main(void)
 {
-  gl_PointSize = point_size;
-  fColor = vec4(color, 1.0);
-  m_vertex = vertex;
-  gl_Position = mvp_matrix * vertex;
+  fColor = vec4(a_Color, 1.0);
+  if (u_UseDefaultColor)
+  {
+    fColor = vec4(u_DefaultColor, 1.0);
+  }
+
+  vec4 pos = vec4(a_Pos, 1.0);
+
+  ls_fP  = pos;
+
+  gl_Position  = u_Mvp * pos;
+
+  float distance = gl_Position.w;
+  if (u_IsOrthographic)
+  {
+    distance = u_PointSize;
+  }
+
+  float effectiveDistance = EqualZero(distance) ? 0.00001 : distance;
+  gl_PointSize = u_PointSize / effectiveDistance * 5.0;
 }
 )DELIM";
 
-const char fragment_source_p_l[]=R"DELIM(
+const char VERTEX_SOURCE_SHAPE[]=R"DELIM(
 #version 150
-in highp vec4 fColor;
-in highp vec4 m_vertex;
+in highp   vec3 a_Pos;
+in mediump vec3 a_Color;
 
-uniform highp vec4 clipPlane;
-uniform highp vec4 pointPlane;
-uniform highp float rendering_mode;
+out mediump vec4 gColor;
+out highp   vec4 ls_fP;
 
-out highp vec4 out_color;
+uniform highp   mat4 u_Mvp;
+uniform mediump vec3 u_DefaultColor;
+uniform         bool u_UseDefaultColor;
+
+void main(void)
+{
+  gColor = vec4(a_Color, 1.0);
+  if (u_UseDefaultColor)
+  {
+    gColor = vec4(u_DefaultColor, 1.0);
+  }
+
+  gl_Position = vec4(a_Pos, 1.0);
+}
+)DELIM";
+
+const char GEOMETRY_SOURCE_SPHERE[]=R"DELIM(
+#version 150
+layout(points) in;
+layout(triangle_strip, max_vertices = 72) out; // max_vertices = (resolution+1) * 2 * latResolution
+
+#define PI 3.14159265358979323846
+
+in mediump vec4 gColor[];
+
+out mediump vec4 fColor;
+out highp   vec4 ls_fP;
+
+uniform highp   mat4 u_Mvp;
+uniform mediump float u_Radius;
+
+void drawSphere(in vec4 center, in float radius, in float resolution)
+{
+  float latResolution = resolution*0.5;
+  float stepTheta = PI/latResolution;
+  float stepPhi = 2*PI/resolution;
+  for(int i=0; i<latResolution; ++i)
+  {
+    float theta1 = stepTheta*i;
+    float theta2 = stepTheta*(i+1);
+    for(int j=0; j<=resolution; ++j)
+    {
+      float phi = stepPhi*j;
+      float x1 = center.x + radius * sin(theta1) * cos(phi);
+      float y1 = center.y + radius * sin(theta1) * sin(phi);
+      float z1 = center.z + radius * cos(theta1);
+      ls_fP = vec4(x1, y1, z1, 1.0);
+      gl_Position = u_Mvp * ls_fP;
+      EmitVertex();
+
+      float x2 = center.x + radius * sin(theta2) * cos(phi);
+      float y2 = center.y + radius * sin(theta2) * sin(phi);
+      float z2 = center.z + radius * cos(theta2);
+      ls_fP = vec4(x2, y2, z2, 1.0);
+      gl_Position = u_Mvp * ls_fP;
+      EmitVertex();
+    }
+    EndPrimitive();
+  }
+}
+
+void main(void)
+{
+  fColor = gColor[0];
+
+  int resolution = 8;
+  vec4 center = gl_in[0].gl_Position;
+
+  drawSphere(center, u_Radius, resolution);
+}
+)DELIM";
+
+const char GEOMETRY_SOURCE_CYLINDER[]=R"DELIM(
+#version 150
+layout(lines) in;
+layout(triangle_strip, max_vertices = 22) out;
+
+#define PI 3.14159265358979323846
+
+in mediump vec4 gColor[];
+
+out mediump vec4 fColor;
+out highp   vec4 ls_fP;
+
+uniform highp   mat4 u_Mvp;
+uniform mediump float u_Radius;
+
+void drawCylinder(in vec3 u, in vec3 v, in vec4 bot, in vec4 top, in float radius, in float resolution)
+{
+  float step = 2*PI/resolution;
+  for(int i=0; i<=resolution; ++i)
+  {
+    float theta = step*i;
+    float cosf = radius*cos(theta);
+    float sinf = radius*sin(theta);
+    vec3 xAxis = cosf*u.xyz;
+    vec3 yAxis = sinf*v.xyz;
+    ls_fP = vec4(top.xyz+xAxis.xyz+yAxis.xyz, 1.0);
+    gl_Position = u_Mvp * ls_fP;
+    EmitVertex();
+    ls_fP = vec4(bot.xyz+xAxis.xyz+yAxis.xyz, 1.0);
+    gl_Position = u_Mvp * ls_fP;
+    EmitVertex();
+  }
+  EndPrimitive();
+}
+
+void main(void)
+{
+  fColor = gColor[0];
+
+  vec4 a = gl_in[0].gl_Position;
+  vec4 b = gl_in[1].gl_Position;
+
+  vec3 n = normalize(vec3(b.x-a.x, b.y-a.y, b.z-a.z)); // compute top normal
+
+  vec3 w = normalize(vec3(-n.z, n.x, n.y));
+
+  // Axis vectors
+  vec3 u = normalize(cross(n, w));
+  vec3 v = normalize(cross(n, u));
+
+  int resolution = 10;
+
+  drawCylinder(u, v, a, b, u_Radius, resolution);
+}
+)DELIM";
+
+const char FRAGMENT_SOURCE_P_L[]=R"DELIM(
+#version 150
+in mediump vec4 fColor;
+in highp   vec4 ls_fP;
+
+out mediump vec4 out_color;
+
+uniform highp   vec4  u_ClipPlane;
+uniform highp   vec4  u_PointPlane;
+uniform mediump float u_RenderingMode;
 
 void main(void)
 {
   // onPlane == 1: inside clipping plane, should be solid;
   // onPlane == -1: outside clipping plane, should be transparent;
   // onPlane == 0: on clipping plane, whatever;
-  float onPlane = sign(dot((m_vertex.xyz-pointPlane.xyz), clipPlane.xyz));
+  float onPlane = sign(dot((ls_fP.xyz-u_PointPlane.xyz), u_ClipPlane.xyz));
 
   // rendering_mode == -1: draw both inside and outside;
   // rendering_mode == 0: draw inside only;
   // rendering_mode == 1: draw outside only;
-  if (rendering_mode == (onPlane+1)/2) {
+  if (u_RenderingMode == (onPlane+1)/2) {
     // discard other than the corresponding half when rendering
     discard;
   }
@@ -150,106 +318,554 @@ void main(void)
 }
 )DELIM";
 
-const char vertex_source_clipping_plane[]=R"DELIM(
+const char VERTEX_SOURCE_CLIPPING_PLANE[]=R"DELIM(
 #version 150
-in highp vec4 vertex;
+in highp vec3 a_Pos;
 
-uniform highp mat4 vp_matrix;
-uniform highp mat4 m_matrix;
+uniform highp mat4 u_Vp;
+uniform highp mat4 u_M;
 
 void main(void)
 {
-  gl_Position = vp_matrix * m_matrix * vertex;
+  gl_Position = u_Vp * u_M * vec4(a_Pos, 1.0);
 }
 )DELIM";
 
-const char fragment_source_clipping_plane[]=R"DELIM(
+const char FRAGMENT_SOURCE_CLIPPING_PLANE[]=R"DELIM(
 #version 150
-out highp vec4 out_color;
+
+out mediump vec4 out_color;
+
 void main(void)
 {
   out_color = vec4(0.0, 0.0, 0.0, 1.0);
 }
 )DELIM";
 
-//------------------------------------------------------------------------------
-//  compatibility shaders
+const char VERTEX_SOURCE_LINE[]=R"DELIM(
+#version 150
+in highp   vec3 a_Pos;
+in mediump vec3 a_Color;
 
-const char vertex_source_color_comp[]=R"DELIM(
-varying highp vec4 vertex;
-varying highp vec3 normal;
-varying highp vec3 color;
+out VS_OUT {
+  mediump vec4 color;
+} vs_out; // vertex shader output
 
-uniform highp mat4 mvp_matrix;
-uniform highp mat4 mv_matrix;
-uniform highp float point_size;
-
-varying highp vec4 fP;
-varying highp vec3 fN;
-varying highp vec4 fColor;
+uniform mediump vec3 u_DefaultColor;
+uniform         bool u_UseDefaultColor;
 
 void main(void)
 {
-  fP = mv_matrix * vertex;
+  vs_out.color = vec4(a_Color, 1.0);
+  if (u_UseDefaultColor)
+  {
+    vs_out.color = vec4(u_DefaultColor, 1.0);
+  }
+
+  gl_Position = vec4(a_Pos, 1.0);
+}
+)DELIM";
+
+const char GEOMETRY_SOURCE_ARROW[]=R"DELIM(
+#version 150
+layout(lines) in;
+layout(triangle_strip, max_vertices = 82) out; // max_vertices = resolution * 2 + 2 (cylinder) + resolution * 3 (disc) + resolution * 3 (cone)
+
+#define PI 3.14159265358979323846
+
+in VS_OUT {
+  mediump vec4 color;
+} gs_in[]; // geometry shader input
+
+out mediump vec4 fColor;
+
+uniform highp   mat4  u_Mvp;
+uniform mediump float u_SceneRadius;
+
+void drawTriangle(in vec4 v1, in vec4 v2, in vec4 v3)
+{
+  gl_Position = u_Mvp*v1;
+  EmitVertex();
+  gl_Position = u_Mvp*v2;
+  EmitVertex();
+  gl_Position = u_Mvp*v3;
+  EmitVertex();
+  EndPrimitive();
+}
+
+void drawTriangleFan(in vec3 u, in vec3 v, in vec4 center, in vec4 edge0, in float radius, in int resolution)
+{
+  float step = 2*PI/resolution;
+  for(int i=0; i<resolution; ++i)
+  {
+    float theta = step*i;
+    float cosf = radius*cos(theta);
+    float sinf = radius*sin(theta);
+    vec3 xAxis = cosf*u.xyz;
+    vec3 yAxis = sinf*v.xyz;
+    vec4 edge1 = vec4(edge0.xyz+xAxis.xyz+yAxis.xyz, 1.0);
+    theta = step*(i+1);
+    cosf = radius*cos(theta);
+    sinf = radius*sin(theta);
+    xAxis = cosf*u.xyz;
+    yAxis = sinf*v.xyz;
+    vec4 edge2 = vec4(edge0.xyz+xAxis.xyz+yAxis.xyz, 1.0);
+    drawTriangle(center, edge1, edge2);
+  }
+}
+
+void drawDisc(in vec3 u, in vec3 v, in vec4 center, in float radius, in int resolution)
+{
+  drawTriangleFan(u, v, center, center, radius, resolution);
+}
+
+void drawCone(in vec3 u, in vec3 v, in vec3 n, in vec4 center, in float radius, in float height, in int resolution)
+{
+  drawTriangleFan(u, v, center, vec4(center.xyz-height*n.xyz, 1.0), radius, resolution);
+}
+
+void drawCylinder(in vec3 u, in vec3 v, in vec4 bot, in vec4 top, in float radius, in float resolution)
+{
+  float step = 2*PI/resolution;
+  for(int i=0; i<=resolution; ++i)
+  {
+    float theta = step*i;
+    float cosf = radius*cos(theta);
+    float sinf = radius*sin(theta);
+    vec3 xAxis = cosf*u.xyz;
+    vec3 yAxis = sinf*v.xyz;
+    gl_Position = u_Mvp * vec4(top.xyz+xAxis.xyz+yAxis.xyz, 1.0);
+    EmitVertex();
+    gl_Position = u_Mvp * vec4(bot.xyz+xAxis.xyz+yAxis.xyz, 1.0);
+    EmitVertex();
+  }
+  EndPrimitive();
+}
+
+void main(void)
+{
+  fColor = gs_in[0].color;
+
+  vec4 a = gl_in[0].gl_Position;
+  vec4 b = gl_in[1].gl_Position;
+
+  vec3 n = normalize(vec3(b.x-a.x, b.y-a.y, b.z-a.z)); // compute top normal
+  vec3 w = normalize(vec3(-n.z, n.x, n.y));
+
+  // Axis vectors
+  vec3 u = normalize(cross(n, w));
+  vec3 v = normalize(cross(n, u));
+
+  float radius = 0.013 * u_SceneRadius;
+  float height = 0.035 * u_SceneRadius;
+  int resolution = 10;
+
+  vec4 c = vec4(b.xyz-height*n.xyz, 1.0);
+  drawDisc(u, v, c, radius, resolution);
+  drawCone(u, v, n, b, radius, height, resolution);
+  drawCylinder(u, v, a, c, radius*0.5, resolution);
+}
+)DELIM";
+
+const char GEOMETRY_SOURCE_LINE[]=R"DELIM(
+#version 150
+layout(lines) in;
+layout(line_strip, max_vertices = 2) out;
+
+in VS_OUT {
+  mediump vec4 color;
+} gs_in[]; // geometry shader input
+
+out mediump vec4 fColor;
+
+uniform highp mat4 u_Mvp;
+
+void main(void)
+{
+  fColor = gs_in[0].color;
+
+  gl_Position = u_Mvp * gl_in[0].gl_Position;
+  EmitVertex();
+  gl_Position = u_Mvp * gl_in[1].gl_Position;
+  EmitVertex();
+
+  EndPrimitive();
+}
+)DELIM";
+
+const char FRAGMENT_SOURCE_LINE[]=R"DELIM(
+#version 150
+
+in mediump vec4 fColor;
+
+out mediump vec4 out_color;
+
+void main(void)
+{
+  out_color = fColor;
+}
+)DELIM";
+
+const char VERTEX_SOURCE_NORMAL[]=R"DELIM(
+#version 150
+in highp vec3 a_Pos;
+in highp vec3 a_Normal;
+
+out VS_OUT {
+  mediump vec4  color;
+  highp   vec3  normal;
+} vs_out; // vertex shader output
+
+uniform highp   mat4 u_Mv;
+uniform mediump vec3 u_DefaultColor;
+uniform         bool u_UseDefaultColor;
+
+void main(void)
+{
+  vs_out.color = vec4(abs(normalize(a_Normal)), 1.0);
+  if (u_UseDefaultColor)
+  {
+    vs_out.color = vec4(u_DefaultColor, 1.0);
+  }
+
+  mat3 normalMatrix = mat3(transpose(inverse(u_Mv)));
+  vs_out.normal = normalize(vec3(vec4(normalMatrix * a_Normal, 0.0)));
+
+  gl_Position = vec4(a_Pos, 1.0);
+}
+)DELIM";
+
+const char GEOMETRY_SOURCE_NORMAL[]=R"DELIM(
+#version 150
+layout (triangles) in;
+layout (line_strip, max_vertices = 6) out;
+
+in VS_OUT {
+  mediump vec4  color;
+  highp   vec3  normal;
+} gs_in[]; // geometry shader input
+
+out mediump vec4 fColor;
+out highp   vec4 ls_fP;
+
+uniform highp   mat4  u_Projection;
+uniform mediump float u_Factor;
+uniform mediump float u_SceneRadius;
+uniform highp   mat4  u_Mv;
+uniform bool          u_DisplayFaceNormal;
+
+void GenerateLine(int index)
+{
+  fColor = gs_in[index].color;
+
+  ls_fP = gl_in[index].gl_Position;
+  gl_Position = u_Projection * u_Mv * gl_in[index].gl_Position;
+  EmitVertex();
+
+  vec4 newPosition = u_Mv * gl_in[index].gl_Position + vec4(gs_in[index].normal, 0.0) * u_SceneRadius * u_Factor;
+  ls_fP = inverse(u_Mv) * newPosition;
+  gl_Position = u_Projection * newPosition;
+  EmitVertex();
+
+  EndPrimitive();
+}
+
+void DrawVerticesNormal()
+{
+  GenerateLine(0); // first vertex normal
+  GenerateLine(1); // second vertex normal
+  GenerateLine(2); // third vertex normal
+}
+
+void DrawFaceNormal()
+{
+  fColor = (gs_in[0].color + gs_in[1].color + gs_in[2].color) / 3;
+  vec4 center = (gl_in[0].gl_Position + gl_in[1].gl_Position + gl_in[2].gl_Position) / 3;
+  ls_fP = center;
+  gl_Position = u_Projection * u_Mv * center;
+  EmitVertex();
+
+  vec3 n = normalize((gs_in[0].normal.xyz + gs_in[1].normal.xyz + gs_in[2].normal.xyz) / 3);
+
+  vec4 newPosition = u_Mv * center + vec4(n, 0.0) * u_SceneRadius * u_Factor;
+  ls_fP = inverse(u_Mv) * newPosition;
+  gl_Position = u_Projection * newPosition;
+  EmitVertex();
+}
+
+void main()
+{
+  if (u_DisplayFaceNormal)
+  {
+    DrawFaceNormal();
+  }
+  else
+  {
+    DrawVerticesNormal();
+  }
+}
+)DELIM";
+
+const char VERTEX_SOURCE_TRIANGLE[]=R"DELIM(
+#version 150
+in highp   vec3 a_Pos;
+
+out VS_OUT {
+  mediump vec4 color;
+  highp   vec4 ls_fP;
+} vs_out;
+
+uniform highp mat4 u_Mvp;
+
+void main(void)
+{
+  vec4 pos = vec4(a_Pos, 1.0);
+
+  vs_out.color = vec4(0.85, 0.85, 0.85, 1.0);
+  vs_out.ls_fP = pos;
+
+  gl_Position = u_Mvp * pos;
+}
+)DELIM";
+
+const char GEOMETRY_SOURCE_TRIANGLE[]=R"DELIM(
+#version 150
+layout (triangles) in;
+layout (line_strip, max_vertices=4) out;
+
+in VS_OUT {
+  mediump vec4 color;
+  highp   vec4 ls_fP;
+} gs_in[];
+
+out mediump vec4 fColor;
+out highp   vec4 ls_fP;
+
+void main(void)
+{
+  fColor = gs_in[0].color;
+  ls_fP = gs_in[0].ls_fP;
+  gl_Position = gl_in[0].gl_Position;
+  EmitVertex();
+
+  ls_fP = gs_in[1].ls_fP;
+  gl_Position = gl_in[1].gl_Position;
+  EmitVertex();
+
+  ls_fP = gs_in[2].ls_fP;
+  gl_Position = gl_in[2].gl_Position;
+  EmitVertex();
+
+  ls_fP = gs_in[0].ls_fP;
+  gl_Position = gl_in[0].gl_Position;
+  EmitVertex();
+}
+)DELIM";
+
+const char VERTEX_SOURCE_LINE_WIDTH[]=R"DELIM(
+#version 150
+
+in highp   vec3 a_Pos;
+in mediump vec3 a_Color;
+
+out VS_OUT {
+  mediump float pointSize;
+  mediump vec4 color;
+  highp   vec4 ls_fP;
+} vs_out;
+
+uniform highp   mat4  u_Mvp;
+uniform mediump float u_PointSize;
+uniform         bool  u_IsOrthographic;
+uniform mediump vec3  u_DefaultColor;
+uniform         bool  u_UseDefaultColor;
+
+bool EqualZero(float value)
+{
+  return abs(value) < 0.00001;
+}
+
+void main(void)
+{
+  vec4 pos = vec4(a_Pos, 1.0);
+
+  vs_out.ls_fP = pos;
+  vs_out.color = vec4(a_Color, 1.0);
+  if (u_UseDefaultColor)
+  {
+    vs_out.color = vec4(u_DefaultColor, 1.0);
+  }
+
+  gl_Position  = u_Mvp * pos;
+
+  float distance = gl_Position.w;
+  if (u_IsOrthographic)
+  {
+    distance = u_PointSize;
+  }
+
+  float effectiveDistance = EqualZero(distance) ? 0.00001 : distance;
+  vs_out.pointSize = u_PointSize / effectiveDistance;
+}
+)DELIM";
+
+const char GEOMETRY_SOURCE_LINE_WIDTH[]=R"DELIM(
+#version 150
+layout (lines) in;
+layout (triangle_strip, max_vertices = 4) out;
+
+in mediump vec4 g_Color[];
+
+in VS_OUT {
+  mediump float pointSize;
+  mediump vec4 color;
+  highp   vec4 ls_fP;
+} gs_in[];
+
+out mediump vec4 fColor;
+out highp   vec4 ls_fP;
+
+uniform mediump float u_PointSize;
+uniform mediump vec2  u_Viewport;
+uniform highp   mat4  u_Mvp;
+
+vec2 ToScreenSpace(vec4 vertex)
+{
+  return vec2(vertex.xy / vertex.w) * u_Viewport;
+}
+
+vec4 ToWorldSpace(vec4 vertex)
+{
+  return vec4((vertex.xy * vertex.w) / u_Viewport, vertex.zw);
+}
+
+void main(void)
+{
+  vec2 p0 = ToScreenSpace(gl_in[0].gl_Position);
+  vec2 p1 = ToScreenSpace(gl_in[1].gl_Position);
+  vec2 v0 = normalize(p1 - p0);
+  vec2 n0 = vec2(-v0.y, v0.x) * u_PointSize * 0.5;
+
+  // line start
+  gl_Position = ToWorldSpace(vec4(p0 - n0 * gs_in[0].pointSize, gl_in[0].gl_Position.zw));
+  fColor = gs_in[0].color;
+  ls_fP = inverse(u_Mvp) * gl_Position;
+  EmitVertex();
+
+  gl_Position = ToWorldSpace(vec4(p0 + n0 * gs_in[0].pointSize, gl_in[0].gl_Position.zw));
+  fColor = gs_in[0].color;
+  ls_fP = inverse(u_Mvp) * gl_Position;
+  EmitVertex();
+
+  // line end
+  gl_Position = ToWorldSpace(vec4(p1 - n0 * gs_in[1].pointSize, gl_in[1].gl_Position.zw));
+  fColor = gs_in[1].color;
+  ls_fP = inverse(u_Mvp) * gl_Position;
+  EmitVertex();
+
+  gl_Position = ToWorldSpace(vec4(p1 + n0 * gs_in[1].pointSize, gl_in[1].gl_Position.zw));
+  fColor = gs_in[1].color;
+  ls_fP = inverse(u_Mvp) * gl_Position;
+  EmitVertex();
+}
+)DELIM";
+
+//------------------------------------------------------------------------------
+//  compatibility shaders
+
+const char VERTEX_SOURCE_COLOR_COMP[]=R"DELIM(
+varying highp   vec3 a_Pos;
+varying highp   vec3 a_Normal;
+varying mediump vec3 a_Color;
+
+varying highp   vec4 vs_fP; // view space position
+varying highp   vec3 fN;
+varying mediump vec4 fColor;
+
+uniform highp   mat4  u_Mvp;
+uniform highp   mat4  u_Mv;
+uniform mediump float u_PointSize;
+uniform mediump vec3  u_DefaultColor;
+uniform         bool  u_UseDefaultColor;
+
+void main(void)
+{
+  vec4 pos = vec4(a_Pos, 1.0);
+
+  vs_fP = u_Mv * pos;
   highp mat3 mv_matrix_3;
   mv_matrix_3[0] = mv_matrix[0].xyz;
   mv_matrix_3[1] = mv_matrix[1].xyz;
   mv_matrix_3[2] = mv_matrix[2].xyz;
-  fN = mv_matrix_3* normal;
-  fColor = vec4(color, 1.0);
-  gl_PointSize = point_size;
+  fN = mv_matrix_3* a_Normal;
 
-  gl_Position = mvp_matrix * vertex;
+  fColor = vec4(a_Color, 1.0);
+  if (u_UseDefaultColor)
+  {
+    fColor = vec4(u_DefaultColor, 1.0);
+  }
+
+  gl_PointSize = u_PointSize;
+
+  gl_Position = u_Mvp * pos;
 }
 )DELIM";
 
-const char fragment_source_color_comp[]=R"DELIM(
-varying highp vec4 fP;
-varying highp vec3 fN;
-varying highp vec4 fColor;
+const char FRAGMENT_SOURCE_COLOR_COMP[]=R"DELIM(
+varying highp   vec4 vs_fP;
+varying highp   vec3 fN;
+varying mediump vec4 fColor;
 
-uniform highp vec4 light_pos;
-uniform highp vec4 light_diff;
-uniform highp vec4 light_spec;
-uniform highp vec4 light_amb;
-uniform highp float spec_power ;
+uniform highp   vec4 u_LightPos;
+uniform mediump vec4 u_LightDiff;
+uniform mediump vec4 u_LightSpec;
+uniform mediump vec4 u_LightAmb;
+uniform mediump float u_SpecPower ;
 
 void main(void)
 {
-  highp vec3 L = light_pos.xyz - fP.xyz;
-  highp vec3 V = -fP.xyz;
+  highp vec3 L = u_LightPos.xyz - vs_fP.xyz;
+  highp vec3 V = -vs_fP.xyz;
 
-  highp vec3 N = normalize(fN);
+  highp vec3 a_Normal = normalize(fN);
   L = normalize(L);
   V = normalize(V);
 
-  highp vec3 R = reflect(-L, N);
-  highp vec4 diffuse = max(dot(N,L), 0.0) * light_diff * fColor;
-  highp vec4 specular = pow(max(dot(R,V), 0.0), spec_power) * light_spec;
+  highp vec3 R = reflect(-L, a_Normal);
+  highp vec4 diffuse = max(dot(a_Normal,L), 0.0) * u_LightDiff * fColor;
+  highp vec4 specular = pow(max(dot(R,V), 0.0), u_SpecPower) * u_LightSpec;
 
-  gl_FragColor = light_amb*fColor + diffuse;
+  gl_FragColor = u_LightAmb*fColor + diffuse;
 }
 )DELIM";
 
-const char vertex_source_p_l_comp[]=R"DELIM(
-varying highp vec4 vertex;
-varying highp vec3 color;
+const char VERTEX_SOURCE_P_L_COMP[]=R"DELIM(
+varying highp   vec3 a_Pos;
+varying mediump vec3 a_Color;
 
-uniform highp mat4 mvp_matrix;
-uniform highp float point_size;
+varying mediump vec4 fColor;
 
-varying highp vec4 fColor;
+uniform highp   mat4  u_Mvp;
+uniform mediump float u_PointSize;
+uniform mediump vec3  u_DefaultColor;
+uniform         bool  u_UseDefaultColor;
 
 void main(void)
 {
-  gl_PointSize = point_size;
-  fColor = vec4(color, 1.0);
-  gl_Position = mvp_matrix * vertex;
+  fColor = vec4(a_Color, 1.0);
+  if (u_UseDefaultColor)
+  {
+    fColor = vec4(u_DefaultColor, 1.0);
+  }
+
+  gl_PointSize = u_PointSize;
+  gl_Position = u_Mvp * vec4(a_Pos, 1.0);
 }
 )DELIM";
 
-const char fragment_source_p_l_comp[]=R"DELIM(
-varying highp vec4 fColor;
+const char FRAGMENT_SOURCE_P_L_COMP[]=R"DELIM(
+varying mediump vec4 fColor;
+
 void main(void)
 {
   gl_FragColor = fColor;

Basic_viewer/include/CGAL/Basic_viewer.h

@@ -0,0 +1,44 @@
+// Copyright (c) 2018  GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Guillaume Damiand <guillaume.damiand@liris.cnrs.fr>
+
+#ifndef CGAL_BASIC_VIEWER_H
+#define CGAL_BASIC_VIEWER_H
+
+#include <CGAL/license/GraphicsView.h>
+#include <CGAL/Graphics_scene.h>
+
+// compatibility
+#if defined(CGAL_USE_BASIC_VIEWER) && !defined(CGAL_USE_BASIC_VIEWER_QT)
+#define CGAL_USE_BASIC_VIEWER_QT 1
+#endif
+
+#if defined(CGAL_USE_BASIC_VIEWER_QT) && !defined(CGAL_USE_BASIC_VIEWER)
+#define CGAL_USE_BASIC_VIEWER 1
+#endif
+
+#if defined(CGAL_USE_BASIC_VIEWER_QT)
+#include <CGAL/Qt/Basic_viewer.h>
+// #elif defined(CGAL_USE_BASIC_VIEWER_GLFW)
+// #include <CGAL/GLFW/Basic_viewer.h>
+#else
+namespace CGAL
+{
+  inline
+  void draw_graphics_scene(const Graphics_scene&,
+                           const char* ="CGAL Basic Viewer")
+  {
+    std::cerr<<"Impossible to draw, CGAL_USE_BASIC_VIEWER is not defined."<<std::endl;
+  }
+} // End namespace CGAL
+#endif
+
+#endif // CGAL_BASIC_VIEWER_H

Basic_viewer/include/CGAL/Buffer_for_vao.h

@@ -146,7 +146,7 @@ namespace internal
         {
           typename CDT::Edge e(fh,i);
           auto n = fh->neighbor(i);
-          if (n->info().m_nesting_level==-1)
+          if (n!=nullptr && n->info().m_nesting_level==-1)
           {
             if (tri.is_constrained(e)) { border.push(e); }
             else { queue.push(n); }
@@ -884,13 +884,13 @@ protected:
   typedef CGAL::Constrained_Delaunay_triangulation_2<P_traits, TDS, Itag>    CDT;
 
 protected:
-  std::vector<BufferType>* m_pos_buffer;
-  std::vector<IndexType>* m_index_buffer;
-  std::vector<BufferType>* m_color_buffer;
-  mutable std::vector<BufferType>* m_flat_normal_buffer;
-  mutable std::vector<BufferType>* m_gouraud_normal_buffer;
+  std::vector<BufferType>* m_pos_buffer=nullptr;
+  std::vector<IndexType>* m_index_buffer=nullptr;
+  std::vector<BufferType>* m_color_buffer=nullptr;
+  mutable std::vector<BufferType>* m_flat_normal_buffer=nullptr;
+  mutable std::vector<BufferType>* m_gouraud_normal_buffer=nullptr;
 
-  CGAL::Bbox_3* m_bb;
+  CGAL::Bbox_3* m_bb=nullptr;
 
   bool m_zero_x; /// True iff all points have x==0
   bool m_zero_y; /// True iff all points have y==0

Basic_viewer/include/CGAL/Graphics_scene.h

@@ -64,64 +64,82 @@ public:
   typedef Local_kernel::Vector_3 Local_vector;
 
   Graphics_scene()
-      : m_buffer_for_mono_points(&arrays[POS_MONO_POINTS], nullptr,
-                                 &m_bounding_box, nullptr, nullptr, nullptr),
-        m_buffer_for_colored_points(&arrays[POS_COLORED_POINTS], nullptr,
-                                    &m_bounding_box, &arrays[COLOR_POINTS],
-                                    nullptr, nullptr),
-        m_buffer_for_mono_segments(&arrays[POS_MONO_SEGMENTS], nullptr,
-                                   &m_bounding_box, nullptr, nullptr, nullptr),
-        m_buffer_for_colored_segments(&arrays[POS_COLORED_SEGMENTS], nullptr,
-                                      &m_bounding_box, &arrays[COLOR_SEGMENTS],
-                                      nullptr, nullptr),
-        m_buffer_for_mono_rays(&arrays[POS_MONO_RAYS], nullptr, &m_bounding_box,
-                               nullptr, nullptr),
-        m_buffer_for_colored_rays(&arrays[POS_COLORED_RAYS], nullptr,
-                                  &m_bounding_box, &arrays[COLOR_RAYS], nullptr,
-                                  nullptr),
-        m_buffer_for_mono_lines(&arrays[POS_MONO_RAYS], nullptr,
-                                &m_bounding_box, nullptr, nullptr),
-        m_buffer_for_colored_lines(&arrays[POS_COLORED_LINES], nullptr,
-                                   &m_bounding_box, &arrays[COLOR_LINES],
-                                   nullptr, nullptr),
-        m_buffer_for_mono_faces(
-            &arrays[POS_MONO_FACES], nullptr, &m_bounding_box, nullptr,
-            &arrays[FLAT_NORMAL_MONO_FACES], &arrays[SMOOTH_NORMAL_MONO_FACES]),
-        m_buffer_for_colored_faces(&arrays[POS_COLORED_FACES], nullptr,
-                                   &m_bounding_box, &arrays[COLOR_FACES],
-                                   &arrays[FLAT_NORMAL_COLORED_FACES],
-                                   &arrays[SMOOTH_NORMAL_COLORED_FACES])
+      : m_buffer_for_points(&arrays[POS_POINTS], nullptr,
+                            &m_bounding_box, &arrays[COLOR_POINTS]),
+        m_buffer_for_segments(&arrays[POS_SEGMENTS], nullptr,
+                              &m_bounding_box, &arrays[COLOR_SEGMENTS]),
+        m_buffer_for_rays(&arrays[POS_RAYS], nullptr, &m_bounding_box,
+                          &arrays[COLOR_RAYS]),
+        m_buffer_for_lines(&arrays[POS_RAYS], nullptr,
+                           &m_bounding_box, &arrays[COLOR_LINES]),
+        m_buffer_for_faces(&arrays[POS_FACES], nullptr, &m_bounding_box, &arrays[COLOR_FACES],
+                           &arrays[FLAT_NORMAL_FACES], &arrays[SMOOTH_NORMAL_FACES]),
+        m_default_color_face(60, 60, 200),
+        m_default_color_point(200, 60, 60),
+        m_default_color_segment(0, 0, 0),
+        m_default_color_ray(0, 0, 0),
+        m_default_color_line(0, 0, 0)
   {}
 
-  const Buffer_for_vao &get_buffer_for_mono_points() const
-  { return m_buffer_for_mono_points; }
+  inline
+  const CGAL::IO::Color &get_default_color_face() const
+  { return m_default_color_face; }
 
-  const Buffer_for_vao &get_buffer_for_colored_points() const
-  { return m_buffer_for_colored_points; }
+  inline
+  const CGAL::IO::Color &get_default_color_point() const
+  { return m_default_color_point; }
 
-  const Buffer_for_vao &get_buffer_for_mono_segments() const
-  { return m_buffer_for_mono_segments; }
+  inline
+  const CGAL::IO::Color &get_default_color_segment() const
+  { return m_default_color_segment; }
 
-  const Buffer_for_vao &get_buffer_for_colored_segments() const
-  { return m_buffer_for_colored_segments; }
+  inline
+  const CGAL::IO::Color &get_default_color_ray() const
+  { return m_default_color_ray; }
 
-  const Buffer_for_vao &get_buffer_for_mono_rays() const
-  { return m_buffer_for_mono_rays; }
+  inline
+  const CGAL::IO::Color &get_default_color_line() const
+  { return m_default_color_line; }
 
-  const Buffer_for_vao &get_buffer_for_colored_rays() const
-  { return m_buffer_for_colored_rays; }
+  inline
+  const Buffer_for_vao &get_buffer_for_points() const
+  { return m_buffer_for_points; }
 
-  const Buffer_for_vao &get_buffer_for_mono_lines() const
-  { return m_buffer_for_mono_lines; }
+  inline
+  void set_default_color_face(const CGAL::IO::Color& c)
+  { m_default_color_face = c; }
 
-  const Buffer_for_vao &get_buffer_for_colored_lines() const
-  { return m_buffer_for_colored_lines; }
+  inline
+  void set_default_color_point(const CGAL::IO::Color& c)
+  { m_default_color_point = c; }
 
-  const Buffer_for_vao &get_buffer_for_mono_faces() const
-  { return m_buffer_for_mono_faces; }
+  inline
+  void set_default_color_segment(const CGAL::IO::Color& c)
+  { m_default_color_segment = c; }
 
-  const Buffer_for_vao &get_buffer_for_colored_faces() const
-  { return m_buffer_for_colored_faces; }
+  inline
+  void set_default_color_ray(const CGAL::IO::Color& c)
+  { m_default_color_ray = c; }
+
+  inline
+  void set_default_color_line(const CGAL::IO::Color& c)
+  { m_default_color_line = c; }
+
+  inline
+  const Buffer_for_vao &get_buffer_for_segments() const
+  { return m_buffer_for_segments; }
+
+  inline
+  const Buffer_for_vao &get_buffer_for_rays() const
+  { return m_buffer_for_rays; }
+
+  inline
+  const Buffer_for_vao &get_buffer_for_lines() const
+  { return m_buffer_for_lines; }
+
+  inline
+  const Buffer_for_vao &get_buffer_for_faces() const
+  { return m_buffer_for_faces; }
 
   const CGAL::Bbox_3 &bounding_box() const { return m_bounding_box; }
 
@@ -159,31 +177,30 @@ public:
 
   void reverse_all_normals() const
   {
-    m_buffer_for_mono_faces.negate_normals();
-    m_buffer_for_colored_faces.negate_normals();
+    m_buffer_for_faces.negate_normals();
   }
 
   template <typename KPoint> void add_point(const KPoint &p)
-  { m_buffer_for_mono_points.add_point(p); }
+  { m_buffer_for_points.add_point(p, m_default_color_point); }
 
   template <typename KPoint>
   void add_point(const KPoint &p, const CGAL::IO::Color &acolor)
-  { m_buffer_for_colored_points.add_point(p, acolor); }
+  { m_buffer_for_points.add_point(p, acolor); }
 
   template <typename KPoint>
   void add_segment(const KPoint &p1, const KPoint &p2)
-  { m_buffer_for_mono_segments.add_segment(p1, p2); }
+  { m_buffer_for_segments.add_segment(p1, p2, m_default_color_segment); }
 
   template <typename KPoint>
   void add_segment(const KPoint &p1, const KPoint &p2,
                    const CGAL::IO::Color &acolor)
-  { m_buffer_for_colored_segments.add_segment(p1, p2, acolor); }
+  { m_buffer_for_segments.add_segment(p1, p2, acolor); }
 
   template <typename KPoint, typename KVector>
   void add_ray(const KPoint &p, const KVector &v)
   {
     double bigNumber = 1e30;
-    m_buffer_for_mono_rays.add_ray_segment(p, (p + (bigNumber)*v));
+    m_buffer_for_rays.add_ray_segment(p, (p + (bigNumber)*v), m_default_color_ray);
   }
 
   template <typename KPoint, typename KVector>
@@ -191,15 +208,15 @@ public:
                const CGAL::IO::Color &acolor)
   {
     double bigNumber = 1e30;
-    m_buffer_for_colored_rays.add_ray_segment(p, (p + (bigNumber)*v), acolor);
+    m_buffer_for_rays.add_ray_segment(p, (p + (bigNumber)*v), acolor);
   }
 
   template <typename KPoint, typename KVector>
   void add_line(const KPoint &p, const KVector &v)
   {
     double bigNumber = 1e30;
-    m_buffer_for_mono_lines.add_line_segment((p - (bigNumber)*v),
-                                             (p + (bigNumber)*v));
+    m_buffer_for_lines.add_line_segment((p - (bigNumber)*v),
+                                             (p + (bigNumber)*v), m_default_color_line);
   }
 
   template <typename KPoint, typename KVector>
@@ -207,33 +224,28 @@ public:
                 const CGAL::IO::Color &acolor)
   {
     double bigNumber = 1e30;
-    m_buffer_for_colored_lines.add_line_segment((p - (bigNumber)*v),
+    m_buffer_for_lines.add_line_segment((p - (bigNumber)*v),
                                                 (p + (bigNumber)*v), acolor);
   }
 
   template <typename KPoint> bool add_point_in_face(const KPoint &kp)
   {
-    if (m_buffer_for_mono_faces.is_a_face_started())
-    { return m_buffer_for_mono_faces.add_point_in_face(kp); }
-    else if (m_buffer_for_colored_faces.is_a_face_started())
-    { return m_buffer_for_colored_faces.add_point_in_face(kp); }
+    if (m_buffer_for_faces.is_a_face_started())
+    { return m_buffer_for_faces.add_point_in_face(kp); }
     return false;
   }
 
   template <typename KPoint, typename KVector>
   bool add_point_in_face(const KPoint &kp, const KVector &p_normal)
   {
-    if (m_buffer_for_mono_faces.is_a_face_started())
-    { return m_buffer_for_mono_faces.add_point_in_face(kp, p_normal); }
-    else if (m_buffer_for_colored_faces.is_a_face_started())
-    { return m_buffer_for_colored_faces.add_point_in_face(kp, p_normal); }
+    if (m_buffer_for_faces.is_a_face_started())
+    { return m_buffer_for_faces.add_point_in_face(kp, p_normal); }
     return false;
   }
 
   bool a_face_started() const
   {
-    return m_buffer_for_mono_faces.is_a_face_started() ||
-           m_buffer_for_colored_faces.is_a_face_started();
+    return m_buffer_for_faces.is_a_face_started();
   }
 
   void face_begin()
@@ -245,7 +257,7 @@ public:
           << std::endl;
     }
     else
-    { m_buffer_for_mono_faces.face_begin(); }
+    { m_buffer_for_faces.face_begin(m_default_color_face); }
   }
 
   void face_begin(const CGAL::IO::Color &acolor)
@@ -257,15 +269,13 @@ public:
           << std::endl;
     }
     else
-    { m_buffer_for_colored_faces.face_begin(acolor); }
+    { m_buffer_for_faces.face_begin(acolor); }
   }
 
   void face_end()
   {
-    if (m_buffer_for_mono_faces.is_a_face_started())
-    { m_buffer_for_mono_faces.face_end(); }
-    else if (m_buffer_for_colored_faces.is_a_face_started())
-    { m_buffer_for_colored_faces.face_end(); }
+    if (m_buffer_for_faces.is_a_face_started())
+    { m_buffer_for_faces.face_end(); }
   }
 
   template <typename KPoint>
@@ -281,58 +291,38 @@ public:
 
   bool empty() const
   {
-    return (m_buffer_for_mono_points.is_empty() &&
-            m_buffer_for_colored_points.is_empty() &&
-            m_buffer_for_mono_segments.is_empty() &&
-            m_buffer_for_colored_segments.is_empty() &&
-            m_buffer_for_mono_rays.is_empty() &&
-            m_buffer_for_colored_rays.is_empty() &&
-            m_buffer_for_mono_lines.is_empty() &&
-            m_buffer_for_colored_lines.is_empty() &&
-            m_buffer_for_mono_faces.is_empty() &&
-            m_buffer_for_colored_faces.is_empty());
+    return (m_buffer_for_points.is_empty() &&
+            m_buffer_for_segments.is_empty() &&
+            m_buffer_for_rays.is_empty() &&
+            m_buffer_for_lines.is_empty() &&
+            m_buffer_for_faces.is_empty());
   }
 
   bool has_zero_x() const
   {
-    return m_buffer_for_mono_points.has_zero_x() &&
-           m_buffer_for_colored_points.has_zero_x() &&
-           m_buffer_for_mono_segments.has_zero_x() &&
-           m_buffer_for_colored_segments.has_zero_x() &&
-           m_buffer_for_mono_faces.has_zero_x() &&
-           m_buffer_for_colored_faces.has_zero_x() &&
-           m_buffer_for_mono_rays.has_zero_x() &&
-           m_buffer_for_colored_rays.has_zero_x() &&
-           m_buffer_for_mono_lines.has_zero_x() &&
-           m_buffer_for_colored_lines.has_zero_x();
+    return m_buffer_for_points.has_zero_x() &&
+           m_buffer_for_segments.has_zero_x() &&
+           m_buffer_for_faces.has_zero_x() &&
+           m_buffer_for_rays.has_zero_x() &&
+           m_buffer_for_lines.has_zero_x();
   }
 
   bool has_zero_y() const
   {
-    return m_buffer_for_mono_points.has_zero_y() &&
-           m_buffer_for_colored_points.has_zero_y() &&
-           m_buffer_for_mono_segments.has_zero_y() &&
-           m_buffer_for_colored_segments.has_zero_y() &&
-           m_buffer_for_mono_faces.has_zero_y() &&
-           m_buffer_for_colored_faces.has_zero_y() &&
-           m_buffer_for_mono_rays.has_zero_y() &&
-           m_buffer_for_colored_rays.has_zero_y() &&
-           m_buffer_for_mono_lines.has_zero_y() &&
-           m_buffer_for_colored_lines.has_zero_y();
+    return m_buffer_for_points.has_zero_y() &&
+           m_buffer_for_segments.has_zero_y() &&
+           m_buffer_for_faces.has_zero_y() &&
+           m_buffer_for_rays.has_zero_y() &&
+           m_buffer_for_lines.has_zero_y();
   }
 
   bool has_zero_z() const
   {
-    return m_buffer_for_mono_points.has_zero_z() &&
-           m_buffer_for_colored_points.has_zero_z() &&
-           m_buffer_for_mono_segments.has_zero_z() &&
-           m_buffer_for_colored_segments.has_zero_z() &&
-           m_buffer_for_mono_faces.has_zero_z() &&
-           m_buffer_for_colored_faces.has_zero_z() &&
-           m_buffer_for_mono_rays.has_zero_z() &&
-           m_buffer_for_colored_rays.has_zero_z() &&
-           m_buffer_for_mono_lines.has_zero_z() &&
-           m_buffer_for_colored_lines.has_zero_z();
+    return m_buffer_for_points.has_zero_z() &&
+           m_buffer_for_segments.has_zero_z() &&
+           m_buffer_for_faces.has_zero_z() &&
+           m_buffer_for_rays.has_zero_z() &&
+           m_buffer_for_lines.has_zero_z();
   }
 
   // Returns true if the data structure lies on a XY or XZ or YZ plane
@@ -343,16 +333,11 @@ public:
 
   void clear()
   {
-    m_buffer_for_mono_points.clear();
-    m_buffer_for_colored_points.clear();
-    m_buffer_for_mono_segments.clear();
-    m_buffer_for_colored_segments.clear();
-    m_buffer_for_mono_rays.clear();
-    m_buffer_for_colored_rays.clear();
-    m_buffer_for_mono_lines.clear();
-    m_buffer_for_colored_lines.clear();
-    m_buffer_for_mono_faces.clear();
-    m_buffer_for_colored_faces.clear();
+    m_buffer_for_points.clear();
+    m_buffer_for_segments.clear();
+    m_buffer_for_rays.clear();
+    m_buffer_for_lines.clear();
+    m_buffer_for_faces.clear();
     m_texts.clear();
     m_bounding_box=CGAL::Bbox_3();
   }
@@ -385,16 +370,11 @@ public:
   // vectors.
   enum Buffers {
     BEGIN_POS = 0,
-    POS_MONO_POINTS = BEGIN_POS,
-    POS_COLORED_POINTS,
-    POS_MONO_SEGMENTS,
-    POS_COLORED_SEGMENTS,
-    POS_MONO_RAYS,
-    POS_COLORED_RAYS,
-    POS_MONO_LINES,
-    POS_COLORED_LINES,
-    POS_MONO_FACES,
-    POS_COLORED_FACES,
+    POS_POINTS = BEGIN_POS,
+    POS_SEGMENTS,
+    POS_RAYS,
+    POS_LINES,
+    POS_FACES,
     END_POS,
     BEGIN_COLOR = END_POS,
     COLOR_POINTS = BEGIN_COLOR,
@@ -404,25 +384,24 @@ public:
     COLOR_FACES,
     END_COLOR,
     BEGIN_NORMAL = END_COLOR,
-    SMOOTH_NORMAL_MONO_FACES = BEGIN_NORMAL,
-    FLAT_NORMAL_MONO_FACES,
-    SMOOTH_NORMAL_COLORED_FACES,
-    FLAT_NORMAL_COLORED_FACES,
+    SMOOTH_NORMAL_FACES = BEGIN_NORMAL,
+    FLAT_NORMAL_FACES,
     END_NORMAL,
     LAST_INDEX = END_NORMAL
   };
 
 protected:
-  Buffer_for_vao m_buffer_for_mono_points;
-  Buffer_for_vao m_buffer_for_colored_points;
-  Buffer_for_vao m_buffer_for_mono_segments;
-  Buffer_for_vao m_buffer_for_colored_segments;
-  Buffer_for_vao m_buffer_for_mono_rays;
-  Buffer_for_vao m_buffer_for_colored_rays;
-  Buffer_for_vao m_buffer_for_mono_lines;
-  Buffer_for_vao m_buffer_for_colored_lines;
-  Buffer_for_vao m_buffer_for_mono_faces;
-  Buffer_for_vao m_buffer_for_colored_faces;
+  Buffer_for_vao m_buffer_for_points;
+  Buffer_for_vao m_buffer_for_segments;
+  Buffer_for_vao m_buffer_for_rays;
+  Buffer_for_vao m_buffer_for_lines;
+  Buffer_for_vao m_buffer_for_faces;
+
+  CGAL::IO::Color m_default_color_face;
+  CGAL::IO::Color m_default_color_point;
+  CGAL::IO::Color m_default_color_segment;
+  CGAL::IO::Color m_default_color_ray;
+  CGAL::IO::Color m_default_color_line;
 
   std::vector<std::tuple<Local_point, std::string>> m_texts;
 

Boolean_set_operations_2/doc/Boolean_set_operations_2/Boolean_set_operations_2.txt

@@ -116,7 +116,7 @@ In our context, a polygon must uphold the following conditions:
 <OL>
 <LI><I>Closed Boundary</I> - the polygon's outer boundary must be a connected sequence of curves, that start and end at the same vertex.
 <LI><I>Simplicity</I> - the polygon must be simple.
-<LI><I>Orientation</I> - the polygon's outer boundary must be <I>counter-clockwise oriented</I>.
+<LI><I>Orientation</I> - the polygon's outer boundary must be <I>counterclockwise oriented</I>.
 </OL>
 
 \subsection bso_ssecpolygon_with_holes_validation Conditions for Valid Polygons with Holes

Boolean_set_operations_2/doc/Boolean_set_operations_2/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgBooleanSetOperations2Ref 2D Regularized Boolean Set-Operations Reference
+/// \defgroup PkgBooleanSetOperations2Ref Reference Manual
 /// \defgroup PkgBooleanSetOperations2Concepts Concepts
 /// \ingroup PkgBooleanSetOperations2Ref
 

Boolean_set_operations_2/include/CGAL/Boolean_set_operations_2/Gps_bfs_xor_visitor.h

@@ -52,7 +52,7 @@ public:
     return (ic % 2) == 1;
   }
 
-  //! after_scan post-processing after bfs scan.
+  //! after_scan postprocessing after bfs scan.
 /*! The function fixes some of the curves, to be in the same direction as the
     half-edges.
 

Boolean_set_operations_2/include/CGAL/draw_polygon_set_2.h

@@ -19,10 +19,10 @@
 
 #include <CGAL/license/Boolean_set_operations_2.h>
 
-#include <CGAL/Qt/Basic_viewer.h>
 #include <CGAL/Graphics_scene.h>
 #include <CGAL/Graphics_scene_options.h>
 #include <CGAL/Polygon_set_2.h>
+#include <CGAL/Basic_viewer.h>
 
 namespace CGAL {
 
@@ -166,7 +166,7 @@ void compute_elements(const PWH& pwh,
 
 } // End namespace draw_function_for_boolean_set_2
 
-#ifdef CGAL_USE_BASIC_VIEWER
+#if defined(CGAL_USE_BASIC_VIEWER)
 
 template <typename PolygonSet_2, typename GSOptions>
 class Polygon_set_2_basic_viewer_qt : public Basic_viewer
@@ -196,6 +196,8 @@ public:
     CGAL::qglviewer::Vec maxv(bbox.xmax(), bbox.ymax(), 0);
     auto diameter = (maxv - minv).norm();
     m_pixel_ratio = diameter / gso.height();
+
+    add_elements();
   }
 
   /*! Intercept the resizing of the window.
@@ -281,13 +283,16 @@ void add_to_graphics_scene(const CGAL_PS2_TYPE& ap2,
   draw_function_for_boolean_set_2::compute_elements(ap2, graphics_scene, gso);
 }
 
-#ifdef CGAL_USE_BASIC_VIEWER
-
 // Specialization of draw function.
 template<class T, class C, class D, class GSOptions>
 void draw(const CGAL_PS2_TYPE& ps, GSOptions& gso,
           const char* title = "Polygon_set_2 Basic Viewer")
 {
+  CGAL_USE(ps);
+  CGAL_USE(gso);
+  CGAL_USE(title);
+
+#ifdef CGAL_USE_BASIC_VIEWER
 #if defined(CGAL_TEST_SUITE)
   bool cgal_test_suite = true;
 #else
@@ -306,6 +311,7 @@ void draw(const CGAL_PS2_TYPE& ps, GSOptions& gso,
     basic_viewer.show();
     app.exec();
   }
+#endif // CGAL_USE_BASIC_VIEWER
 }
 
 template<class T, class C, class D>
@@ -319,7 +325,6 @@ void draw(const CGAL_PS2_TYPE& ps,
   draw(ps, gso, title);
 }
 
-#endif // CGAL_USE_BASIC_VIEWER
 
 #undef CGAL_PS2_TYPE
 

Bounding_volumes/doc/Bounding_volumes/CGAL/Min_sphere_annulus_d_traits_2.h

@@ -16,7 +16,6 @@ optimization algorithms using the two-dimensional \cgal kernel.
 \sa `CGAL::Min_annulus_d<Traits>`
 \sa `CGAL::Min_sphere_annulus_d_traits_3<K,ET,NT>`
 \sa `CGAL::Min_sphere_annulus_d_traits_d<K,ET,NT>`
-\sa `MinSphereAnnulusDTraits`
 
 */
 template< typename K, typename ET, typename NT >

Bounding_volumes/doc/Bounding_volumes/CGAL/Min_sphere_annulus_d_traits_3.h

@@ -16,8 +16,6 @@ optimization algorithms using the three-dimensional \cgal kernel.
 \sa `CGAL::Min_annulus_d<Traits>`
 \sa `CGAL::Min_sphere_annulus_d_traits_2<K,ET,NT>`
 \sa `CGAL::Min_sphere_annulus_d_traits_d<K,ET,NT>`
-\sa `Min_sphere_annulusDTraits`
-
 */
 template< typename K, typename ET, typename NT >
 class Min_sphere_annulus_d_traits_3 {

Bounding_volumes/doc/Bounding_volumes/CGAL/Min_sphere_annulus_d_traits_d.h

@@ -17,7 +17,6 @@ optimization algorithms using the \f$ d\f$-dimensional \cgal kernel.
 \sa `CGAL::Polytope_distance_d<Traits>`
 \sa `CGAL::Min_sphere_annulus_d_traits_2<K,ET,NT>`
 \sa `CGAL::Min_sphere_annulus_d_traits_3<K,ET,NT>`
-\sa `MinSphereAnnulusDTraits`
 
 */
 template< typename K, typename ET, typename NT >

Bounding_volumes/doc/Bounding_volumes/Concepts/MinQuadrilateralTraits_2.h

@@ -74,25 +74,25 @@ typedef unspecified_type Equal_2;
 
 /*!
 a model for
-`Kernel::Less_xy_2`.
+`Kernel::LessXY_2`.
 */
 typedef unspecified_type Less_xy_2;
 
 /*!
 a model for
-`Kernel::Less_yx_2`.
+`Kernel::LessYX_2`.
 */
 typedef unspecified_type Less_yx_2;
 
 /*!
 a model for
-`Kernel::Has_on_negative_side_2`.
+`Kernel::HasOnNegativeSide_2`.
 */
 typedef unspecified_type Has_on_negative_side_2;
 
 /*!
 a model for
-`Kernel::Compare_angle_with_x_axis_2`.
+`Kernel::CompareAngleWithXAxis_2`.
 */
 typedef unspecified_type Compare_angle_with_x_axis_2;
 
@@ -101,7 +101,7 @@ AdaptableFunctor
 `op`:
 `Rectangle_2` \f$ \times\f$ `Rectangle_2` \f$ \rightarrow\f$
 `bool`.
-`op(r1,r2)` returns true, iff the area of \f$ r1\f$ is
+`op(r1,r2)` returns `true`, iff the area of \f$ r1\f$ is
 strictly less than the area of \f$ r2\f$.
 */
 typedef unspecified_type Area_less_rectangle_2;
@@ -112,7 +112,7 @@ AdaptableFunctor
 `Parallelogram_2` \f$ \times\f$
 `Parallelogram_2` \f$ \rightarrow\f$ `bool`.
 
-`op(p1,p2)` returns true, iff the area of \f$ p1\f$ is strictly less
+`op(p1,p2)` returns `true`, iff the area of \f$ p1\f$ is strictly less
 than the area of \f$ p2\f$.
 */
 typedef unspecified_type Area_less_parallelogram_2;
@@ -122,7 +122,7 @@ AdaptableFunctor
 `op`:
 `Strip_2` \f$ \times\f$ `Strip_2` \f$ \rightarrow\f$ `bool`.
 
-`op(s1,s2)` returns true, iff the width of \f$ s1\f$ is strictly less
+`op(s1,s2)` returns `true`, iff the width of \f$ s1\f$ is strictly less
 than the width of \f$ s2\f$.
 */
 typedef unspecified_type Width_less_strip_2;
@@ -143,7 +143,7 @@ typedef unspecified_type Orientation_2;
 
 /*!
 a model for
-`Kernel::Construct_vector_2`.
+`Kernel::ConstructVector_2`.
 */
 typedef unspecified_type Construct_vector_2;
 
@@ -158,25 +158,25 @@ typedef unspecified_type Construct_vector_from_direction_2;
 
 /*!
 a model for
-`Kernel::Construct_perpendicular_vector_2`.
+`Kernel::ConstructPerpendicularVector_2`.
 */
 typedef unspecified_type Construct_perpendicular_vector_2;
 
 /*!
 a model for
-`Kernel::Construct_direction_2`.
+`Kernel::ConstructDirection_2`.
 */
 typedef unspecified_type Construct_direction_2;
 
 /*!
 a model for
-`Kernel::Construct_opposite_direction_2`.
+`Kernel::ConstructOppositeDirection_2`.
 */
 typedef unspecified_type Construct_opposite_direction_2;
 
 /*!
 a model for
-`Kernel::Construct_line_2`.
+`Kernel::ConstructLine_2`.
 */
 typedef unspecified_type Construct_line_2;
 

Bounding_volumes/doc/Bounding_volumes/Concepts/RectangularPCenterTraits_2.h

@@ -34,31 +34,31 @@ typedef unspecified_type Point_2;
 
 /*!
 model for
-`Kernel::Iso_rectangle_2`.
+`Kernel::IsoRectangle_2`.
 */
 typedef unspecified_type Iso_rectangle_2;
 
 /*!
 model for
-`Kernel::Less_x_2`.
+`Kernel::LessX_2`.
 */
 typedef unspecified_type Less_x_2;
 
 /*!
 model for
-`Kernel::Less_y_2`.
+`Kernel::LessY_2`.
 */
 typedef unspecified_type Less_y_2;
 
 /*!
 model for
-`Kernel::Construct_vertex_2`.
+`Kernel::ConstructVertex_2`.
 */
 typedef unspecified_type Construct_vertex_2;
 
 /*!
 model for
-`Kernel::Construct_iso_rectangle_2`.
+`Kernel::ConstructIsoRectangle_2`.
 */
 typedef unspecified_type Construct_iso_rectangle_2;
 

Bounding_volumes/doc/Bounding_volumes/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgBoundingVolumesRef Bounding Volumes Reference
+/// \defgroup PkgBoundingVolumesRef Reference Manual
 /// \defgroup PkgBoundingVolumesConcepts Concepts
 /// \ingroup PkgBoundingVolumesRef
 /*!

Bounding_volumes/doc/Bounding_volumes/dependencies

@@ -6,3 +6,4 @@ Circulator
 Stream_support
 Matrix_search
 Polytope_distance_d
+Number_types

Box_intersection_d/doc/Box_intersection_d/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgBoxIntersectionDRef Intersecting Sequences of dD Iso-oriented Boxes Reference
+/// \defgroup PkgBoxIntersectionDRef Reference Manual
 /// \defgroup PkgBoxIntersectionDConcepts Concepts
 /// \ingroup PkgBoxIntersectionDRef
 

CGAL_ipelets/doc/CGAL_ipelets/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgCGALIpeletsRef CGAL Ipelets Reference
+/// \defgroup PkgCGALIpeletsRef Reference Manual
 /*!
 \addtogroup PkgCGALIpeletsRef
 \cgalPkgDescriptionBegin{CGAL Ipelets,PkgCGALIpelets}

Circular_kernel_2/doc/Circular_kernel_2/Concepts/AlgebraicKernelForCircles.h

@@ -29,7 +29,7 @@ A model of `RingNumberType`.
 typedef unspecified_type RT;
 
 /*!
-A model of `FieldNumberType``<RT>`.
+A model of `FieldNumberType<RT>`.
 */
 typedef unspecified_type FT;
 

Circular_kernel_2/doc/Circular_kernel_2/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgCircularKernel2Ref 2D Circular Geometry Kernel Reference
+/// \defgroup PkgCircularKernel2Ref Reference Manual
 
 /// \defgroup PkgCircularKernel2GeometricConcepts Geometric Concepts
 /// \ingroup PkgCircularKernel2Ref

Circular_kernel_3/doc/Circular_kernel_3/Concepts/AlgebraicKernelForSpheres--PolynomialForSpheres_2_3.h

@@ -12,7 +12,7 @@ are of a type that is a model of the concept
 \cgalRefines{CopyConstructible,Assignable,DefaultConstructible}
 
 \cgalHasModelsBegin
-\cgalHasModels{GAL::Polynomial_for_spheres_2_}
+\cgalHasModels{CGAL::Polynomial_for_spheres_2_}
 \cgalHasModelsEnd
 
 \sa `AlgebraicKernelForSpheres`

Circular_kernel_3/doc/Circular_kernel_3/Concepts/AlgebraicKernelForSpheres.h

@@ -29,7 +29,7 @@ A model of `RingNumberType`.
 typedef unspecified_type RT;
 
 /*!
-A model of `FieldNumberType``<RT>`.
+A model of `FieldNumberType<RT>`.
 */
 typedef unspecified_type FT;
 

Circular_kernel_3/doc/Circular_kernel_3/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgCircularKernel3Ref 3D Spherical Geometry Kernel Reference
+/// \defgroup PkgCircularKernel3Ref Reference Manual
 
 /// \defgroup PkgCircularKernel3GeometricConcepts Geometric Concepts
 /// \ingroup PkgCircularKernel3Ref

Circulator/doc/Circulator/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgHandlesAndCirculatorsRef Handles and Circulators Reference
+/// \defgroup PkgHandlesAndCirculatorsRef Reference Manual
 /// \defgroup PkgHandlesAndCirculatorsConcepts Concepts
 /// \ingroup PkgHandlesAndCirculatorsRef
 

Circulator/include/CGAL/circulator.h

@@ -22,10 +22,11 @@
 #include <CGAL/assertions.h>
 #include <CGAL/use.h>
 #include <CGAL/tags.h>
+#include <CGAL/type_traits.h>
 
 #include <cstddef>
-#include <functional>
 #include <iterator>
+#include <utility>
 
 // These are name redefinitions for backwards compatibility
 // with the pre iterator-traits style adaptors.
@@ -640,6 +641,34 @@ operator+( Dist n, const Iterator_from_circulator<C,Ref,Ptr>& circ) {
     return tmp += n;
 }
 
+template <class Circ>
+class Range_from_circulator {
+private:
+  Circ anchor;
+public:
+  using pointer = CGAL::cpp20::remove_cvref_t<decltype(anchor.operator->())>;
+  using const_pointer = CGAL::cpp20::remove_cvref_t<decltype(std::as_const(anchor).operator->())>;
+  using reference = decltype(*anchor);
+  using const_reference = decltype(*std::as_const(anchor));
+  using iterator = Iterator_from_circulator<Circ, reference, pointer>;
+  using const_iterator = Iterator_from_circulator<Circ, const_reference, const_pointer>;
+
+  iterator begin() {
+    return iterator(&anchor, 0);
+  }
+  const_iterator begin() const {
+    return const_iterator(&anchor, 0);
+  }
+  iterator end() {
+    return anchor == nullptr ? iterator(&anchor, 0) : iterator(&anchor, 1);
+  }
+  const_iterator end() const {
+    return anchor == nullptr ? const_iterator(&anchor, 0) : const_iterator(&anchor, 1);
+  }
+  Range_from_circulator() = default;
+  Range_from_circulator(const Circ& c) : anchor(get_min_circulator(c)) {}
+};
+
 template < class  C >
 class Container_from_circulator {
 private:

Classification/doc/Classification/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgClassificationRef Classification Reference
+/// \defgroup PkgClassificationRef Reference Manual
 
 /*!
 \defgroup PkgClassificationConcepts Concepts

Combinatorial_map/doc/Combinatorial_map/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgCombinatorialMapsRef Combinatorial Maps Reference
+/// \defgroup PkgCombinatorialMapsRef Reference Manual
 
 /// \defgroup PkgCombinatorialMapsConcepts Concepts
 /// \ingroup PkgCombinatorialMapsRef

Cone_spanners_2/doc/Cone_spanners_2/PackageDescription.txt

@@ -1,4 +1,4 @@
-/// \defgroup PkgConeSpanners2Ref Cone-Based Spanners Reference
+/// \defgroup PkgConeSpanners2Ref Reference Manual
 
 /*!
 \addtogroup PkgConeSpanners2Ref

Cone_spanners_2/include/CGAL/Compute_cone_boundaries_2.h

@@ -137,7 +137,7 @@ public:
 
       The direction of the first ray can be specified by the parameter
       initial_direction, which allows the first ray to start at any direction.
-      The remaining directions are calculated in counter-clockwise order.
+      The remaining directions are calculated in counterclockwise order.
 
       \param cone_number The number of cones
       \param initial_direction The direction of the first ray

Cone_spanners_2/include/CGAL/Construct_theta_graph_2.h

@@ -169,7 +169,7 @@ protected:
     /* Construct edges in one cone bounded by two directions.
 
      \param cwBound      The direction of the clockwise boundary of the cone.
-     \param ccwBound     The direction of the counter-clockwise boundary.
+     \param ccwBound     The direction of the counterclockwise boundary.
      \param g            The Theta graph to be built.
     */
     void add_edges_in_cone(const Direction_2& cwBound, const Direction_2& ccwBound, Graph_& g) {

Cone_spanners_2/include/CGAL/Construct_yao_graph_2.h

@@ -164,7 +164,7 @@ protected:
     /* Construct edges in one cone bounded by two directions.
 
      \param cwBound      The direction of the clockwise boundary of the cone.
-     \param ccwBound     The direction of the counter-clockwise boundary.
+     \param ccwBound     The direction of the counterclockwise boundary.
      \param g            The Yao graph to be built.
     */
     void add_edges_in_cone(const Direction_2& cwBound, const Direction_2& ccwBound, Graph_& g) {

Constrained_triangulation_3/doc/Constrained_triangulation_3/Concepts/ConformingConstrainedDelaunayTriangulationCellBase_3.h

@@ -0,0 +1,30 @@
+/*!
+\ingroup PkgConstrainedTriangulation3Concepts
+\cgalConcept
+
+The concept `ConformingConstrainedDelaunayTriangulationCellBase_3` refines the concept
+`TriangulationCellBase_3` and and describes the requirements for a base cell class of
+the `CGAL::Conforming_constrained_Delaunay_triangulation_3` class.
+
+\cgalRefines{TriangulationCellBase_3, BaseWithTimeStamp}
+
+\cgalHasModelsBegin
+\cgalHasModels{CGAL::Conforming_constrained_Delaunay_triangulation_cell_base_3}
+\cgalHasModelsEnd
+
+\sa `ConformingConstrainedDelaunayTriangulationVertexBase_3`
+*/
+class ConformingConstrainedDelaunayTriangulationCellBase_3 {
+public:
+
+  /// @name Access Functions
+  ///
+  /// The following functions return a reference to an object of type
+  /// `CGAL::Conforming_constrained_Delaunay_triangulation_cell_data_3`, that contains
+  /// the per-cell data required by the implementation of the
+  /// `CGAL::Conforming_constrained_Delaunay_triangulation_3` class.
+  /// @{
+  CGAL::Conforming_constrained_Delaunay_triangulation_cell_data_3& ccdt_3_data();
+  const CGAL::Conforming_constrained_Delaunay_triangulation_cell_data_3& ccdt_3_data() const;
+  /// @}
+};

Constrained_triangulation_3/doc/Constrained_triangulation_3/Concepts/ConformingConstrainedDelaunayTriangulationTraits_3.h

@@ -0,0 +1,87 @@
+/*!
+\ingroup PkgConstrainedTriangulation3Concepts
+\cgalConcept
+
+The concept `ConformingConstrainedDelaunayTriangulationTraits_3` specifies the requirements
+for the geometric traits class of the triangulation used as the first template
+parameter `%Triangulation` in the function template
+`CGAL::make_conforming_constrained_Delaunay_triangulation_3()`.
+
+\cgalRefines{DelaunayTriangulationTraits_3, ProjectionTraitsGeometricTraits_3}
+
+\cgalHasModelsBegin
+\cgalHasModels{CGAL::Exact_predicates_inexact_constructions_kernel (recommended)}
+\cgalHasModelsBare{All models of the concept `Kernel`}
+\cgalHasModelsEnd
+
+*/
+class ConformingConstrainedDelaunayTriangulationTraits_3 {
+public:
+
+  /// \name Operations
+  /// The following functions give access to the predicate and construction objects:
+  /// @{
+
+  /*!
+   * returns a function object model of `Kernel::Angle_3`
+   */
+  unspecified_type angle_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::CompareAngle_3`
+   */
+  unspecified_type compare_angle_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ComputeScalarProduct_3`
+   */
+  unspecified_type compute_scalar_product_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ComputeSquaredLength_3`
+   */
+  unspecified_type compute_squared_length_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ConstructCrossProductVector_3`
+   */
+  unspecified_type construct_cross_product_vector_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ConstructMidpoint_3`
+   */
+  unspecified_type construct_midpoint_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ConstructVector_3`
+   */
+  unspecified_type construct_vector_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ConstructScaledVector_3`
+   */
+  unspecified_type construct_scaled_vector_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ConstructSumOfVectors_3`
+   */
+  unspecified_type construct_sum_of_vectors_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ConstructTranslatedPoint_3`
+   */
+  unspecified_type construct_translated_point_3_object();
+
+  /*!
+   * returns a function object model of `Kernel::ConstructVertex_3`
+   */
+  unspecified_type construct_vertex_3_object();
+
+  /*!
+   * returns a predicate object model of `Kernel::IsDegenerate_3`
+   */
+  unspecified_type is_degenerate_3_object();
+
+  /// @}
+
+};

Constrained_triangulation_3/doc/Constrained_triangulation_3/Concepts/ConformingConstrainedDelaunayTriangulationVertexBase_3.h

@@ -0,0 +1,31 @@
+/*!
+\ingroup PkgConstrainedTriangulation3Concepts
+\cgalConcept
+
+The concept `ConformingConstrainedDelaunayTriangulationVertexBase_3` refines the concept
+`TriangulationVertexBase_3` and describes the requirements for a base vertex class of
+the `CGAL::Conforming_constrained_Delaunay_triangulation_3` class.
+
+\cgalRefines{TriangulationVertexBase_3, BaseWithTimeStamp}
+
+\cgalHasModelsBegin
+\cgalHasModels{CGAL::Conforming_constrained_Delaunay_triangulation_vertex_base_3}
+\cgalHasModelsEnd
+
+\sa `ConformingConstrainedDelaunayTriangulationCellBase_3`
+
+*/
+class ConformingConstrainedDelaunayTriangulationVertexBase_3 {
+public:
+
+  /// @name Access Functions
+  ///
+  /// The following functions return a reference to an object of type
+  /// `CGAL::Conforming_constrained_Delaunay_triangulation_vertex_data_3`, that contains
+  /// the per-vertex data required by the implementation of the
+  /// `CGAL::Conforming_constrained_Delaunay_triangulation_3` class.
+  /// @{
+  CGAL::Conforming_constrained_Delaunay_triangulation_vertex_data_3& ccdt_3_data();
+  const CGAL::Conforming_constrained_Delaunay_triangulation_vertex_data_3& ccdt_3_data() const;
+  /// @}
+}

Constrained_triangulation_3/doc/Constrained_triangulation_3/Constrained_triangulation_3.txt

@@ -0,0 +1,345 @@
+namespace CGAL {
+/*!
+
+\mainpage User Manual
+\anchor Chapter_CT_3
+\anchor userchapterct3
+
+\cgalAutoToc
+\author Laurent Rineau and Jane Tournois
+
+\cgalFigureAnchor{CT_3_pyramid_fig}
+<img src="cdt_title_pyramid.png" style="max-width:60%;min-width=20%"/>
+<BR>
+
+\section CT_3_CCDT_3 Constrained Triangulations in 3D
+
+3D triangulations partition space and are useful in many applications. In some cases, it is
+important to ensure that specific faces, such as those representing the sharp features of an object, appear in the output.
+When a triangulation exactly respects these constraints, it is called a _constrained_ triangulation.
+However, it is sometimes only possible to preserve the geometry of the constraints, but not their exact
+combinatorics. In such cases, additional points, called _Steiner_ _points_, must be inserted. This process
+results in a _conforming_ triangulation.
+
+This package implements an algorithm for constructing conforming triangulations of 3D polygonal
+constraints. Specifically, it requires that these piecewise linear constraints are provided as a
+_piecewise linear complex_&nbsp;(PLC). The resulting triangulations are of type `Triangulation_3`,
+as described in the chapter \ref PkgTriangulation3.
+
+The article by Cohen-Steiner et al. \cgalCite{cgal:cohen2002conforming} discusses the problem of
+constructing conforming Delaunay triangulations and proposes an algorithm to address it.
+Si et al.'s work \cgalCite{si2005meshing}, \cgalCite{cgal:si2008cdt3}, \cgalCite{si2015tetgen},
+ presents an algorithm for computing conforming constrained
+Delaunay triangulations in 3D.
+
+\section CT_3_definitions Definitions
+
+This section introduces the key concepts necessary to understand and use this package effectively.
+
+\subsection CT_3_PLC Piecewise Linear Complex
+
+A _piecewise linear complex_&nbsp;(PLC) is the three-dimensional generalization of a
+planar straight-line graph. It consists of a finite set of vertices, edges, and polygonal faces
+that satisfy the following properties:
+
+- The vertices and edges of the PLC form a simplicial complex: two edges may intersect only at a
+  shared vertex.
+- The boundary of each polygonal face in the PLC is an ordered list of vertices from the PLC, forming
+  one closed loop.
+- Each polygonal face must be a simple polygon, i.e., its edges don't intersect,
+  except consecutive edges, which intersect at their common vertex.
+- Each polygonal face must be planar, meaning all its vertices lie on the same plane.
+- Each polygonal face may have one or more holes, each of them also represented by an ordered list of vertices
+  from the PLC, forming a closed loop.
+- If two polygons in the PLC intersect, their intersection is a union of edges and vertices from the
+  PLC. In particular, the interiors of two polygons cannot overlap.
+
+Polygons in a PLC may be non-convex and may have holes.
+
+\cgalFigureAnchor{CT_3_plc_fig}
+<center>
+<img src="plc.png" style="max-width:60%;"/>
+</center>
+\cgalFigureCaptionBegin{CT_3_plc_fig}
+A piecewise linear complex, composed of planar faces connected by edges and vertices.
+\cgalFigureCaptionEnd
+
+
+\subsection CT_3_CDT Conforming Constrained Delaunay Triangulation
+
+The algorithms developed in this package are designed to compute a constrained Delaunay
+triangulation that contains a given set of polygonal constraints in 3D as a subcomplex.
+
+A triangulation is a _Delaunay triangulation_ if the circumscribing sphere of any simplex
+in the triangulation contains no vertex in its interior (see chapter \ref PkgTriangulation3
+for more details on Delaunay triangulations).
+
+A _constrained Delaunay triangulation_ of a PLC is a constrained triangulation that is as close as
+possible to being Delaunay, given that some faces are marked as _constrained_. More precisely, a
+triangulation is _constrained Delaunay_ if, for any simplex \f$s\f$ of the triangulation, the
+interior of its circumscribing sphere contains no vertex of the triangulation that is _visible_ from
+any point in the interior of the simplex \f$s\f$. Two points are _visible_ if the open line segment
+joining them does not intersect any polygonal face of the PLC, except for polygons that are coplanar with
+the segment.
+
+In 3D, constrained triangulations do not always exist. This can be demonstrated using the example of
+Sch&ouml;nhardt polyhedra \cgalCite{s-udzvd-28} (see \cgalFigureRef{CT_3_schonhardt_fig}),
+\cgalCite{b-ip-48a}. Shewchuk \cgalCite{cgal:shewchuk1998condition} demonstrated that for any PLC,
+there exists a refined
+version of the original PLC that admits a constrained Delaunay triangulation. This refinement is
+achieved by adding Steiner vertices to the input edges and polygons. The constrained triangulation
+built on this refined PLC is known as a _conforming constrained Delaunay triangulation_ (CCDT for
+short). \cgalFigureRef{CT_3_plc2cdt_fig} illustrates an example of a conforming constrained
+Delaunay triangulation constructed from a PLC.
+
+
+\cgalFigureAnchor{CT_3_schonhardt_fig}
+<center>
+<img src="schonhardt.png" style="max-width:25%;"/>
+</center>
+\cgalFigureCaptionBegin{CT_3_schonhardt_fig}
+A Sch&ouml;nhardt polyhedron.
+\cgalFigureCaptionEnd
+
+
+\cgalFigureAnchor{CT_3_plc2cdt_fig}
+<center>
+<img src="plc_to_cdt.png" style="max-width:70%;"/>
+</center>
+\cgalFigureCaptionBegin{CT_3_plc2cdt_fig}
+Left: PLC (360 vertices);
+Right: CCDT (2452 vertices).
+\cgalFigureCaptionEnd
+
+The algorithm implemented in this package is based on the work of Hang Si et al., who developed particular
+algorithms for constructing conforming constrained Delaunay triangulations from PLCs.
+The corresponding implementation takes with floating point numbers as coordinates
+\cgalCite{si2005meshing}, \cgalCite{cgal:si2008cdt3}, \cgalCite{si2015tetgen}.
+
+
+\section CT_3_design Software Design
+
+\subsection CT_3_representation_of_PLCs Representation of Piecewise Linear Complexes
+
+There is no universal or canonical way to represent all possible PLCs in \cgal.
+
+Any polyhedral surface is a PLC, so any model of `FaceListGraph`, such as `CGAL::Surface_mesh`, can be
+used to represent such a PLC.
+In this representation, the geometric structure of the PLC is directly mapped to the elements
+of the `CGAL::Surface_mesh`:
+- vertices of the PLC geometrically correspond to vertices of the surface mesh,
+- edges of the PLC correspond to edges of the surface mesh,
+- and polygonal faces of the PLC correspond to faces of the surface mesh,
+covering the surface of a geometric object.
+However, PLCs represented in this way are limited to being
+manifold (that is, each edge belongs to exactly two faces), and their faces cannot have holes.
+
+A PLC can also be represented as a polygon soup: a collection of vertices and a set of polygons, where
+each polygon is defined by an ordered list of vertices, without explicit connectivity information between
+polygons. For a polygon soup to represent a valid PLC, its polygons must satisfy the properties described
+in the previous section. This approach allows for the representation of non-manifold geometries; however,
+polygons in a polygon soup cannot have holes.
+
+This package also provides a way to group polygons into distinct surface patches using a property map,
+named `plc_face_id`.
+Each polygon can be assigned a _patch_ identifier, allowing multiple polygons to form a continuous surface patch,
+which may include holes. Some necessary geometric conditions must be satisfied for these patches to be
+used in the conforming constrained Delaunay triangulation construction:
+- Each patch must be planar, meaning all polygonal faces in the patch lie on the same plane;
+- The polygonal faces of the patch must not intersect except at their shared edges.
+
+When this property map is provided, the input PLC is interpreted in terms of its polygonal faces,
+edges and vertices as follows:
+- Each polygonal face of the PLC is defined as the union of input polygons sharing the same patch identifier;
+- The edges of the PLC are those from the surface mesh or polygon soup that satisfy one of the following conditions:
+-- they are adjacent to only one polygonal face;
+-- they are adjacent to two polygonal faces with different patch identifiers;
+-- they are adjacent to more than two polygonal faces with differing patch identifiers, indicating non-manifold features of the PLC.
+- The vertices of the PLC are the ones lying on the boundaries of surface patches in the original surface mesh or polygon soup.
+
+\subsection CT_3_api API
+
+This package provides a primary class, `CGAL::Conforming_constrained_Delaunay_triangulation_3`.
+This class is templated by a
+geometric traits class and an underlying triangulation class, allowing for flexibility and
+customization.
+
+In addition to the main class, the package includes several auxiliary classes that define the types
+of vertices, cells, and associated metadata used within the triangulation. These supporting classes
+enable users to extend or adapt the triangulation data structure to their specific needs.
+
+Two overloads of the constructor function \link
+PkgConstrainedTriangulation3FunctionsPolygonSoupOrMesh
+`CGAL::make_conforming_constrained_Delaunay_triangulation_3()`\endlink are provided to facilitate the creation
+of a `CGAL::Conforming_constrained_Delaunay_triangulation_3` object from either a surface mesh or a polygon soup.
+
+\subsection CT_3_geomtraits Traits and Kernel Choice
+
+The requirements for geometric objects and operations are specified by the traits class concept
+`ConformingConstrainedDelaunayTriangulationTraits_3`. Any CGAL kernel is a model of this concept.
+However, because this package builds upon the 3D Triangulation package, it inherits the requirement
+that the traits class must provide exact predicates.
+
+A key aspect of this algorithm is the creation of new points, known as Steiner points, which are
+inserted on the segments and polygons of the input PLC. If a traits class with inexact constructions
+is used, it cannot be guaranteed that these points will lie exactly on the intended segments or polygons.
+As a result, the output will only approximate the input, with the accuracy limited by the rounding
+of the computed Steiner points.
+
+Furthermore, when using inexact constructions, the algorithm may fail if the input PLC contains
+non-adjacent simplices that are too close to each other. In such cases, the triangulation process
+will emit an error if the distance between simplices falls below an internally computed threshold.
+An error message describing the involved simplices will be displayed on the standard output.
+If the issue is caused by poorly shaped triangles, functions such as
+`CGAL::Polygon_mesh_processing::remove_almost_degenerate_faces()` may help resolve the problem.
+
+\section CT_3_examples Examples
+
+\subsection CT_3_example_ccdt Build a Conforming Constrained Delaunay Triangulation
+
+The following example illustrates how to use the helper function
+`CGAL::make_conforming_constrained_Delaunay_triangulation_3()` to construct a conforming constrained
+Delaunay triangulation from a given PLC.
+
+The triangulation is saved in the MEDIT file format, using the
+\link PkgCDT3IOFunctions `CGAL::IO::write_MEDIT()` \endlink function.
+
+\cgalExample{Constrained_triangulation_3/conforming_constrained_Delaunay_triangulation_3.cpp }
+
+
+\subsection CT_3_example_ccdt_soup Build a Conforming Constrained Delaunay Triangulation from a Polygon Soup
+
+You can also construct a conforming constrained Delaunay triangulation from a polygon soup.
+The following example demonstrates how to create such a triangulation from a collection of polygons
+without explicit connectivity information.
+
+\cgalExample{Constrained_triangulation_3/ccdt_3_from_soup.cpp }
+
+
+\subsection CT_3_example_ccdt_fimap Build a Conforming Constrained Delaunay Triangulation with Known Polygon Identifiers
+
+If the user already knows the set of polygonal face identifiers to associate with each PLC face, this information can be
+provided and preserved throughout the construction of the conforming constrained Delaunay
+triangulation.
+
+The following example demonstrates how to detect planar surface patches and remesh them as coarsely
+as possible using
+\link CGAL::Polygon_mesh_processing::remesh_planar_patches(const TriangleMeshIn&,PolygonMeshOut&,const NamedParametersIn&,const NamedParametersOut&) `CGAL::Polygon_mesh_processing::remesh_planar_patches()` \endlink
+from the \ref PkgPolygonMeshProcessing package.
+The resulting patches and segmentation are then used to build a conforming constrained Delaunay triangulation.
+
+When the named parameter `plc_face_id` is specified, each constrained facet in the 3D triangulation
+is assigned to the corresponding input PLC face, as identified by the provided property map.
+If this parameter is not specified, each input polygonal face is assigned a unique face index.
+
+Figure \cgalFigureRef{CT_3_ccdt_examples_fig} shows the benefit of using the `plc_face_id` property map.
+On the last line of the figure, the input PLC is enriched with a segmentation of the planar faces,
+provided via the `plc_face_id` property map. In the resulting conforming constrained Delaunay triangulation,
+only the boundary edges of the PLC faces are constrained, while the other edges never get inserted as
+edges of the 3D triangulation.
+
+Without the `plc_face_id` property map, all edges of the PLC faces are constrained,
+each PLC face is considered as a constraint,
+possibly resulting in a 3D triangulation with surfaces that are more refined than necessary.
+
+\cgalExample{Constrained_triangulation_3/conforming_constrained_Delaunay_triangulation_3_fimap.cpp}
+
+\cgalFigureRef{CT_3_ccdt_examples_fig} shows the input and output of this triangulation construction example.
+
+
+\subsection CT_3_example_ccdt_region_growing_fimap Build a Conforming Constrained Delaunay Triangulation with Detected Polygon Identifiers
+
+If the user does not know the set of polygonal face identifiers to associate with each PLC face, this information can be
+automatically detected using the
+\link CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces(const PolygonMesh& mesh,RegionMap region_map,const NamedParameters& np) `CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces()`\endlink
+function from the \ref PkgPolygonMeshProcessing package.
+
+The following example demonstrates how to detect planar surface patches and build a conforming
+constrained Delaunay triangulation using the detected segmentation. The named parameter `plc_face_id`
+is used to associate each facet of the triangulation with the corresponding input PLC face.
+
+\cgalExample{Constrained_triangulation_3/ccdt_3_fimap_region_growing.cpp}
+
+
+\subsection CT_3_examples_preprocessing Preprocessing the Input for Conforming Constrained Delaunay Triangulations
+
+Given a PLC, the algorithms in this package can construct a conforming constrained Delaunay triangulation, provided
+the input surface can be represented as a valid surface mesh or a collection of surface meshes,
+and does not contain self-intersections.
+Several preprocessing functions are available in the \ref PkgPolygonMeshProcessing package to help ensure these
+preconditions are met.
+
+The following example demonstrates how to construct a conforming constrained Delaunay triangulation from
+an input mesh that is not triangulated and may contain self-intersections,
+using autorefinement.
+
+\cgalExample{Constrained_triangulation_3/ccdt_3_after_autorefinement.cpp }
+
+The function
+\link CGAL::Polygon_mesh_processing::does_self_intersect(const FaceRange&, const TriangleMesh&, const NamedParameters&) `CGAL::Polygon_mesh_processing::does_self_intersect()` \endlink
+is used to detect self-intersections,
+but it requires the input mesh to be triangulated. Therefore, the input mesh must first be triangulated
+using
+\link CGAL::Polygon_mesh_processing::triangulate_faces(FaceRange,PolygonMesh&,const NamedParameters&) `CGAL::Polygon_mesh_processing::triangulate_faces()`\endlink
+before performing the self-intersection check.
+
+If self-intersections are found, the triangulated mesh is converted into a triangle soup, which is then
+processed with
+\link CGAL::Polygon_mesh_processing::autorefine_triangle_soup(PointRange&,TriangleRange&,const NamedParameters&) `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()`\endlink
+to resolve the self-intersections.
+
+
+\subsection CT_3_example_remesh Remeshing a Conforming Constrained Delaunay Triangulation
+
+After constructing the triangulation, you can improve its quality or adapt it to a specific sizing
+field by applying the
+\link CGAL::tetrahedral_isotropic_remeshing(CGAL::Triangulation_3<Traits, TDS, SLDS>&, const SizingFunction&, const NamedParameters&) `CGAL::tetrahedral_isotropic_remeshing()`\endlink
+ function from the \ref PkgTetrahedralRemeshing package.
+
+The following example demonstrates how to remesh a conforming constrained Delaunay triangulation.
+
+\cgalExample{Constrained_triangulation_3/remesh_constrained_Delaunay_triangulation_3.cpp }
+
+
+\subsection CT_3_examples_figure Figures
+
+The following table of figures (\cgalFigureRef{CT_3_ccdt_examples_fig}) illustrates some results of the examples provided in this package.
+The left column shows the input PLC, while the right column displays the resulting conforming
+constrained Delaunay triangulation.
+
+From top to bottom, the lines show different input PLC, from the same input triangulated surface and,
+for each of them, the resulting conforming constrained Delaunay triangulation.
+The input data are:
+<ul>
+<li> the input PLC with no preprocessing;
+<li> the input PLC with a segmentation of the surface patches done with `CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces()` as done in example \ref CT_3_example_ccdt_region_growing_fimap;</li>
+<li> the input PLC, remeshed and not triangulated, with a segmentation of the surface patches, done with `CGAL::Polygon_mesh_processing::remesh_planar_patches()` as in example \ref CT_3_example_ccdt_fimap;</li>
+<li> the input PLC, isotropically remeshed, with a segmentation of the surface patches, done with `CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces()` as in example \ref CT_3_example_ccdt_region_growing_fimap.</li>
+</ul>
+
+On the fourth line, the input PLC is remeshed using `CGAL::Polygon_mesh_processing::isotropic_remeshing()`.
+The resulting conforming constrained Delaunay triangulation contains fewer vertices than the input
+remeshed and segmented input PLC. This reduction occurs because only the boundary edges of the PLC faces
+are marked as constraints in the triangulation;
+interior edges that do not lie on the boundaries of surface patches (as defined by `plc_face_id`) are ignored.
+As a result, these non-boundary edges are omitted from the triangulation, leading to a coarser triangulation.
+
+\cgalFigureAnchor{CT_3_ccdt_examples_fig}
+<center>
+<img src="ccdt_examples.png" style="max-width:50%;"/>
+</center>
+\cgalFigureCaptionBegin{CT_3_ccdt_examples_fig}
+A collection of conforming constrained Delaunay triangulations built from different inputs.
+The left column shows the input PLC, while the right column displays the resulting 3D triangulation.<br>
+\cgalFigureCaptionEnd
+
+
+\section CT_3_history Implementation History
+
+The initial version of this package was implemented by Laurent Rineau and released in
+\cgal&nbsp;6.1&nbsp;(2025). Jane Tournois contributed to the documentation and helped improve the API.
+The package design and algorithms are grounded in the theoretical work of
+Hang&nbsp;Si et al. on meshing algorithms&nbsp;\cgalCite{si2005meshing},&nbsp;\cgalCite{si2015tetgen}.
+
+*/
+} /* namespace CGAL */

Constrained_triangulation_3/doc/Constrained_triangulation_3/Doxyfile.in

@@ -0,0 +1,9 @@
+@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
+PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Constrained Triangulations"
+HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/plc.png \
+                   ${CGAL_PACKAGE_DOC_DIR}/fig/schonhardt.png \
+                   ${CGAL_PACKAGE_DOC_DIR}/fig/plc_to_cdt.png \
+                   ${CGAL_PACKAGE_DOC_DIR}/fig/ccdt_fpmap.png \
+                   ${CGAL_PACKAGE_DOC_DIR}/fig/cdt_title_pyramid.png \
+                   ${CGAL_PACKAGE_DOC_DIR}/fig/ccdt_examples.png
+QUIET = YES

Constrained_triangulation_3/doc/Constrained_triangulation_3/PackageDescription.txt

@@ -0,0 +1,77 @@
+/*!
+\defgroup PkgConstrainedTriangulation3Ref Reference Manual
+
+\defgroup PkgConstrainedTriangulation3Concepts Concepts
+\ingroup PkgConstrainedTriangulation3Ref
+
+\defgroup PkgConstrainedTriangulation3Classes Classes and Class Templates
+\ingroup PkgConstrainedTriangulation3Ref
+
+\defgroup PkgConstrainedTriangulation3FunctionsPolygonSoupOrMesh  Free Functions for Creating Conforming Constrained Delaunay Triangulations
+\ingroup PkgConstrainedTriangulation3Ref
+
+\defgroup PkgCDT3IOFunctions Output Functions
+\ingroup PkgConstrainedTriangulation3Ref
+
+\defgroup PkgDrawCDT_3 Draw a 3D Constrained Triangulation
+\ingroup PkgConstrainedTriangulation3Ref
+
+\defgroup PkgConstrainedTriangulation3Functions Other Functions
+\ingroup PkgConstrainedTriangulation3Ref
+
+\addtogroup PkgConstrainedTriangulation3Ref
+\cgalPkgDescriptionBegin{3D Constrained Triangulations,PkgConstrainedTriangulation3}
+\cgalPkgPicture{small-pyramid.png}
+
+\cgalPkgSummaryBegin
+\cgalPkgAuthors{Laurent Rineau and Jane Tournois}
+\cgalPkgDesc{This package implements the construction of a 3D Constrained Delaunay triangulation.
+This triangulation is a generalization of a 3D Delaunay Triangulation which conforms to the set of faces
+of a 3D _piecewise linear complex_ (PLC), ensuring that these faces are part of the triangulation.
+As not all PLCs are tetrahedralizable,
+the algorithm may insert Steiner points to construct the constrained triangulation.
+}
+\cgalPkgManuals{Chapter_CT_3,PkgConstrainedTriangulation3Ref}
+\cgalPkgSummaryEnd
+
+\cgalPkgShortInfoBegin
+\cgalPkgSince{6.1}
+\cgalPkgDependsOn{\ref PkgTriangulation3 "3D Triangulations"}
+\cgalPkgBib{cgal:rt-cdt3}
+\cgalPkgLicense{\ref licensesGPL "GPL"}
+\cgalPkgDemo{CGAL Lab,CGALlab.zip}
+\cgalPkgShortInfoEnd
+
+
+\cgalPkgDescriptionEnd
+
+\cgalClassifedRefPages
+
+\cgalCRPSection{Concepts}
+
+- `ConformingConstrainedDelaunayTriangulationTraits_3`
+- `ConformingConstrainedDelaunayTriangulationVertexBase_3`
+- `ConformingConstrainedDelaunayTriangulationCellBase_3`
+
+\cgalCRPSection{Functions}
+
+- \link PkgConstrainedTriangulation3FunctionsPolygonSoupOrMesh `CGAL::make_conforming_constrained_Delaunay_triangulation_3()` \endlink: the main function to create
+a conforming constrained Delaunay triangulation in 3D, an instance of the class template
+`CGAL::Conforming_constrained_Delaunay_triangulation_3`.
+- `CGAL::Tetrahedral_remeshing::get_remeshing_triangulation()`
+
+\cgalCRPSubsection{Classes}
+
+- `CGAL::Conforming_constrained_Delaunay_triangulation_3<Traits, Triangulation>`
+- `CGAL::Conforming_constrained_Delaunay_triangulation_vertex_base_3<Traits, Vertex_base>`
+- `CGAL::Conforming_constrained_Delaunay_triangulation_cell_base_3<Traits, Cell_base>`
+
+\cgalCRPSection{Draw a 3D Constrained Triangulation}
+
+- \link PkgDrawCDT_3 `CGAL::draw()` \endlink
+
+\cgalCRPSection{Output Functions}
+
+- \link PkgCDT3IOFunctions `CGAL::IO::write_MEDIT()` \endlink
+
+*/

Constrained_triangulation_3/doc/Constrained_triangulation_3/dependencies

@@ -0,0 +1,13 @@
+Basic_viewer
+BGL
+Kernel_23
+Manual
+Polygon_mesh_processing
+Property_map
+SMDS_3
+STL_Extension
+Stream_support
+Surface_mesh
+TDS_3
+Tetrahedral_remeshing
+Triangulation_3

Constrained_triangulation_3/doc/Constrained_triangulation_3/examples.txt

@@ -0,0 +1,48 @@
+/*!
+\example Constrained_triangulation_3/conforming_constrained_Delaunay_triangulation_3.cpp
+@brief
+
+Simple example demonstrating the usage of the Constrained_triangulation_3 package.<br>
+
+It constructs a 3D constrained Delaunay triangulation from a polygon mesh.
+*/
+
+/*!
+\example Constrained_triangulation_3/ccdt_3_from_soup.cpp
+@brief
+
+Simple example demonstrating the usage of the Constrained_triangulation_3 package.<br>
+
+It constructs a constrained Delaunay triangulation from a polygon soup.
+*/
+
+/*!
+\example Constrained_triangulation_3/remesh_constrained_Delaunay_triangulation_3.cpp
+@brief
+
+How to use `CGAL::tetrahedral_isotropic_remeshing` with a constrained Delaunay triangulation.
+*/
+
+
+/*!
+\example Constrained_triangulation_3/conforming_constrained_Delaunay_triangulation_3_fimap.cpp
+@brief
+
+From a non-manifold OFF file, construct the constrained Delaunay triangulation.
+*/
+
+
+/*!
+\example Constrained_triangulation_3/ccdt_3_after_autorefinement.cpp
+@brief
+
+From a self-intersecting and non-triangulated surface in an OFF file,
+construct the constrained Delaunay triangulation after preprocessing by autorefinement.
+*/
+
+/*!
+\example Constrained_triangulation_3/ccdt_3_check_preconditions.cpp
+\example Constrained_triangulation_3/ccdt_3_fimap_region_growing.cpp
+\example Constrained_triangulation_3/ccdt_3_from_soup_fimap.cpp
+\example Constrained_triangulation_3/ccdt_3_preprocessing.cpp
+*/

Constrained_triangulation_3/examples/Constrained_triangulation_3/CMakeLists.txt

@@ -0,0 +1,40 @@
+cmake_minimum_required(VERSION 3.12...3.31)
+project(Constrained_triangulation_3_Examples)
+
+find_package(CGAL REQUIRED COMPONENTS Qt6)
+
+find_package(Eigen3 QUIET)
+include(CGAL_Eigen3_support)
+
+add_compile_definitions(QT_NO_KEYWORDS)
+
+create_single_source_cgal_program(conforming_constrained_Delaunay_triangulation_3.cpp)
+create_single_source_cgal_program(conforming_constrained_Delaunay_triangulation_3_fimap.cpp)
+create_single_source_cgal_program(ccdt_3_from_soup.cpp)
+create_single_source_cgal_program(ccdt_3_from_soup_fimap.cpp)
+create_single_source_cgal_program(ccdt_3_after_autorefinement.cpp)
+create_single_source_cgal_program(ccdt_3_preprocessing.cpp)
+create_single_source_cgal_program(ccdt_3_check_preconditions.cpp)
+create_single_source_cgal_program(remesh_constrained_Delaunay_triangulation_3.cpp)
+
+if(TARGET CGAL::Eigen3_support)
+  create_single_source_cgal_program(ccdt_3_fimap_region_growing.cpp)
+  target_link_libraries(ccdt_3_fimap_region_growing PUBLIC CGAL::Eigen3_support)
+endif()
+
+if(CGAL_Qt6_FOUND)
+  target_link_libraries(conforming_constrained_Delaunay_triangulation_3 PUBLIC CGAL::CGAL_Basic_viewer)
+  target_link_libraries(ccdt_3_from_soup PUBLIC CGAL::CGAL_Basic_viewer)
+  target_link_libraries(ccdt_3_from_soup_fimap PUBLIC CGAL::CGAL_Basic_viewer)
+  target_link_libraries(ccdt_3_after_autorefinement PUBLIC CGAL::CGAL_Basic_viewer)
+  target_link_libraries(ccdt_3_preprocessing PUBLIC CGAL::CGAL_Basic_viewer)
+  target_link_libraries(ccdt_3_check_preconditions PUBLIC CGAL::CGAL_Basic_viewer)
+else()
+  message(STATUS "NOTICE: The example 'conforming_constrained_Delaunay_triangulation_3' cannot draw the result without Qt6.")
+endif()
+
+find_package(TBB QUIET)
+include(CGAL_TBB_support)
+if(TARGET CGAL::TBB_support)
+  target_link_libraries(ccdt_3_from_soup_fimap PUBLIC CGAL::TBB_support)
+endif()

Constrained_triangulation_3/examples/Constrained_triangulation_3/ccdt_3_after_autorefinement.cpp

@@ -0,0 +1,71 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/IO/polygon_mesh_io.h>
+#include <CGAL/IO/write_MEDIT.h>
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_self_intersections.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/Surface_mesh/Surface_mesh.h>
+#include <CGAL/draw_constrained_triangulation_3.h>
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point = K::Point_3;
+using Surface_mesh = CGAL::Surface_mesh<Point>;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char* argv[])
+{
+  const auto filename = (argc > 1) ? argv[1]
+                       : CGAL::data_file_path("meshes/spheres_intersecting.off");
+
+  CGAL::Surface_mesh<K::Point_3> mesh;
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh)) {
+    std::cerr << "Error: cannot read file " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "Number of facets in " << filename << ": "
+    << mesh.number_of_faces() << "\n";
+
+  CGAL::Conforming_constrained_Delaunay_triangulation_3<K> ccdt;
+
+  if(PMP::does_self_intersect(mesh))
+  {
+    std::cout << "Mesh self-intersects, performing autorefine...\n";
+
+    // use a polygon soup as container as the output will most likely be non-manifold
+    std::vector<Point> points;
+    std::vector<std::vector<std::size_t>> polygons;
+    PMP::polygon_mesh_to_polygon_soup(mesh, points, polygons);
+    PMP::autorefine_triangle_soup(points, polygons);
+    std::cout << "Number of input triangles after autorefine: "
+              << polygons.size() << "\n";
+
+    if(PMP::does_polygon_soup_self_intersect(points, polygons))
+    {
+      std::cerr << "Error: mesh still self-intersects after autorefine\n";
+      std::cerr << "Run autorefinement again with an exact kernel ";
+      std::cerr << "such as CGAL::Exact_predicates_exact_constructions_kernel \n";
+      return EXIT_FAILURE;
+    }
+
+    ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(points, polygons);
+
+  }
+  else
+  {
+    ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(mesh);
+  }
+  std::cout << "Number of constrained facets in the CDT: "
+            << ccdt.number_of_constrained_facets() << '\n';
+
+  std::ofstream ofs(argc > 2 ? argv[2] : "out.mesh");
+  ofs.precision(17);
+  CGAL::IO::write_MEDIT(ofs, ccdt);
+
+  CGAL::draw(ccdt);
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/examples/Constrained_triangulation_3/ccdt_3_check_preconditions.cpp

@@ -0,0 +1,108 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Surface_mesh/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/boost/graph/copy_face_graph.h>
+#include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
+
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+
+#include <CGAL/draw_constrained_triangulation_3.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/IO/write_MEDIT.h>
+
+#include <vector>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point = K::Point_3;
+using Surface_mesh = CGAL::Surface_mesh<Point>;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+// Function to verify preconditions for a mesh input
+bool verify_preconditions_mesh(const Surface_mesh& mesh)
+{
+  if(CGAL::is_triangle_mesh(mesh))
+    return !CGAL::Polygon_mesh_processing::does_self_intersect(mesh);
+
+  Surface_mesh triangle_mesh;
+  CGAL::copy_face_graph(mesh, triangle_mesh);
+  bool tri_ok = CGAL::Polygon_mesh_processing::triangulate_faces(triangle_mesh);
+  if(!tri_ok)
+    return false;
+
+  return !CGAL::Polygon_mesh_processing::does_self_intersect(triangle_mesh);
+}
+
+// Function to verify preconditions for a polygon soup input
+template <typename PointRange, typename PolygonRange>
+bool verify_preconditions_soup(const PointRange& points, const PolygonRange& polygons)
+{
+  return !CGAL::Polygon_mesh_processing::does_polygon_soup_self_intersect(points, polygons);
+}
+
+int main(int argc, char* argv[])
+{
+  const auto filename = (argc > 1) ? argv[1]
+                       : CGAL::data_file_path("meshes/cubes.off");
+
+  // Read polygon soup
+  using Points = std::vector<Point>;
+  using PLC_face = std::vector<std::size_t>;
+  using PLC_faces = std::vector<PLC_face>;
+  Points points;
+  PLC_faces faces;
+  if(!CGAL::IO::read_polygon_soup(filename, points, faces))
+  {
+    std::cerr << "Error: cannot read file " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  // When possible, convert polygon soup to polygon mesh
+  CGAL::Surface_mesh<K::Point_3> mesh;
+  if(CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(faces))
+  {
+    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, faces, mesh);
+    std::cout << "Number of facets in " << filename << ": " << mesh.number_of_faces() << "\n";
+  }
+
+  // Verify preconditions for the mesh input
+  const bool is_polygon_mesh = !mesh.is_empty();
+  if(is_polygon_mesh && !verify_preconditions_mesh(mesh))
+  {
+    std::cerr << "Error: input mesh is not a valid input for CCDT_3\n";
+    return EXIT_FAILURE;
+  }
+  else if(!verify_preconditions_soup(points, faces))
+  {
+    std::cerr << "Error: input polygon soup is not a valid input for CCDT_3\n";
+    return EXIT_FAILURE;
+  }
+
+  // Build conforming constrained Delaunay triangulation
+  auto ccdt = is_polygon_mesh
+    ? CGAL::make_conforming_constrained_Delaunay_triangulation_3(mesh)
+    : CGAL::make_conforming_constrained_Delaunay_triangulation_3(points, faces);
+
+  // Print mesh details
+  if(ccdt.number_of_constrained_facets() == 0)
+  {
+    std::cerr << "Error: no constrained facets in the CDT.\n";
+    std::cerr << "Invalid input.\n";
+    return EXIT_SUCCESS;
+  }
+
+  std::cout << "Number of constrained facets in the CDT: "
+            << ccdt.number_of_constrained_facets() << '\n';
+
+  // Output
+  std::ofstream ofs(argc > 2 ? argv[2] : "out.mesh");
+  ofs.precision(17);
+  CGAL::IO::write_MEDIT(ofs, ccdt);
+
+  CGAL::draw(ccdt);
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/examples/Constrained_triangulation_3/ccdt_3_fimap_region_growing.cpp

@@ -0,0 +1,76 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/write_MEDIT.h>
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/region_growing.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+
+using Mesh = CGAL::Surface_mesh<K::Point_3>;
+using face_descriptor = boost::graph_traits<Mesh>::face_descriptor;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char* argv[])
+{
+  std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross_quad.off");
+
+  Mesh mesh;
+  if(!PMP::IO::read_polygon_mesh(filename, mesh)) {
+    std::cerr << "Invalid input." << std::endl;
+    return 1;
+  }
+
+  std::cout << "Read " << mesh.number_of_vertices() << " vertices and "
+                       << mesh.number_of_faces() << " facets\n";
+
+  auto [face_patch_map, _] =mesh.add_property_map<face_descriptor, std::size_t>("f:patch_id");
+
+  const auto bbox = CGAL::Polygon_mesh_processing::bbox(mesh);
+  const double bbox_max_span = (std::max)({bbox.x_span(), bbox.y_span(), bbox.z_span()});
+
+  std::cout << "Merging facets into coplanar patches..." << std::endl;
+
+  auto number_of_patches = CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces(
+      mesh,
+      face_patch_map,
+      CGAL::parameters::maximum_distance(bbox_max_span * 1.e-6)
+                       .maximum_angle(5.));
+
+  for(auto f: faces(mesh))
+  {
+    // if region growing did not assign a patch id, assign one
+    if(get(face_patch_map, f) == static_cast<std::size_t>(-1)) {
+      put(face_patch_map, f, number_of_patches++);
+    }
+  }
+
+  std::cout << "Number of patches: " << number_of_patches << std::endl;
+
+  filename = argc > 2 ? argv[2] : "mesh.ply";
+  CGAL::IO::write_polygon_mesh(filename, mesh,
+                               CGAL::parameters::stream_precision(17)
+                                                .use_binary_mode(false)
+                                                .face_patch_map(face_patch_map));
+  std::cout << "-- Wrote segmented mesh to \"" << filename << "\"\n";
+
+  std::cout << "Creating a conforming constrained Delaunay triangulation...\n";
+  auto ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(mesh,
+                CGAL::parameters::plc_face_id(face_patch_map));
+
+  std::cout << "Number of vertices in the CDT: "
+            << ccdt.triangulation().number_of_vertices() << '\n'
+            << "Number of constrained facets in the CDT: "
+            << ccdt.number_of_constrained_facets() << '\n';
+
+  // Write the CDT to a file, with the PLC face ids
+  filename = argc > 3 ? argv[3] : "out.mesh";
+  std::ofstream out(filename);
+  out.precision(17);
+  CGAL::IO::write_MEDIT(out, ccdt, CGAL::parameters::with_plc_face_id(true));
+  std::cout << "-- Wrote CDT to \"" << filename << "\"\n";
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/examples/Constrained_triangulation_3/ccdt_3_from_soup.cpp

@@ -0,0 +1,39 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/IO/write_MEDIT.h>
+#include <CGAL/draw_constrained_triangulation_3.h>
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+
+#include <vector>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+
+int main(int argc, char* argv[])
+{
+  auto filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cubes.off");
+
+  std::vector<K::Point_3> points;
+  std::vector<std::vector<std::size_t>> polygons;
+  if(!CGAL::IO::read_polygon_soup(filename, points, polygons)) {
+    std::cerr << "Error: cannot read file " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "Read " << points.size() << " vertices and "
+            << polygons.size() << " polygons" << std::endl;
+
+  auto ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(points, polygons);
+
+  std::cout << "Number of vertices in the CDT: "
+            << ccdt.triangulation().number_of_vertices() << '\n'
+            << "Number of constrained facets in the CDT: "
+            << ccdt.number_of_constrained_facets() << '\n';
+
+  std::ofstream ofs(argc > 2 ? argv[2] : "out.mesh");
+  ofs.precision(17);
+  CGAL::IO::write_MEDIT(ofs, ccdt);
+
+  CGAL::draw(ccdt);
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/examples/Constrained_triangulation_3/ccdt_3_from_soup_fimap.cpp

@@ -0,0 +1,68 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+#include <CGAL/IO/write_MEDIT.h>
+#include <CGAL/Polygon_mesh_processing/connected_components.h>
+#include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
+#include <CGAL/Surface_mesh/Surface_mesh.h>
+#include <CGAL/draw_constrained_triangulation_3.h>
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+
+#include <algorithm>
+#include <vector>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point = K::Point_3;
+using Surface_mesh = CGAL::Surface_mesh<Point>;
+using face_descriptor = Surface_mesh::Face_index;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char* argv[])
+{
+  auto filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cubes_one_patch_id_per_cc.off");
+
+  CGAL::Surface_mesh<K::Point_3> mesh;
+    if(!CGAL::IO::read_polygon_mesh(filename, mesh)) {
+    std::cerr << "Error: cannot read file " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  auto fpmap = mesh.add_property_map<face_descriptor, std::size_t>("f:cc_id").first;
+
+  std::cout << "Read " << mesh.number_of_vertices() << " vertices and "
+            << mesh.number_of_faces() << " polygons" << std::endl;
+
+  PMP::connected_components(mesh, fpmap);
+
+  std::cout << "Number of connected components: "
+            << fpmap[*std::max_element(mesh.faces_begin(), mesh.faces_end(),
+                                [&](face_descriptor f1, face_descriptor f2) {
+                                  return fpmap[f1] < fpmap[f2];
+                                })] + 1
+            << std::endl;
+
+  std::vector<K::Point_3> points;
+  std::vector<std::vector<std::size_t>> polygons;
+  PMP::polygon_mesh_to_polygon_soup(mesh, points, polygons);
+
+  auto polygon_to_patch_id = [&](std::size_t i) {
+    return fpmap[*std::next(mesh.faces_begin(), i)];
+  };
+
+  auto soup_fpmap = boost::make_function_property_map<std::size_t>(polygon_to_patch_id);
+  auto ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(
+      points, polygons, CGAL::parameters::plc_face_id(soup_fpmap));
+
+  std::cout << "Number of vertices in the CDT: "
+            << ccdt.triangulation().number_of_vertices() << '\n'
+            << "Number of constrained facets in the CDT: "
+            << ccdt.number_of_constrained_facets() << '\n';
+
+  std::ofstream ofs(argc > 2 ? argv[2] : "out.mesh");
+  ofs.precision(17);
+  CGAL::IO::write_MEDIT(ofs, ccdt);
+
+  CGAL::draw(ccdt);
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/examples/Constrained_triangulation_3/ccdt_3_preprocessing.cpp

@@ -0,0 +1,95 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/IO/write_MEDIT.h>
+#include <CGAL/Surface_mesh/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/Polygon_mesh_processing/autorefinement.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/boost/graph/copy_face_graph.h>
+
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+
+#include <CGAL/draw_constrained_triangulation_3.h>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point = K::Point_3;
+using Surface_mesh = CGAL::Surface_mesh<Point>;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char* argv[])
+{
+  const auto filename = (argc > 1) ? argv[1]
+                       : CGAL::data_file_path("meshes/mpi_and_sphere.off");
+
+  CGAL::Surface_mesh<K::Point_3> mesh;
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh)) {
+    std::cerr << "Error: cannot read file " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "Number of facets in " << filename << ": "
+            << mesh.number_of_faces() << "\n";
+
+  // Check if the mesh is a triangle mesh
+  bool triangle_mesh = CGAL::is_triangle_mesh(mesh);
+  if(!triangle_mesh)
+  {
+    std::cout << "Mesh is not a triangle mesh, triangulate faces"
+              << " to check self-intersections...\n";
+
+    CGAL::Surface_mesh<K::Point_3> trimesh;
+    CGAL::copy_face_graph(mesh, trimesh);
+    PMP::triangulate_faces(trimesh);
+
+    if(PMP::does_self_intersect(trimesh))
+    {
+      std::cout << "Mesh self-intersects, let's keep the triangulated version"
+                << " for future autorefinement\n";
+      CGAL::copy_face_graph(trimesh, mesh);
+      mesh = std::move(trimesh);
+      triangle_mesh = true;
+    }
+  }
+
+  CGAL::Conforming_constrained_Delaunay_triangulation_3<K> ccdt;
+  if(triangle_mesh && PMP::does_self_intersect(mesh))
+  {
+    std::cout << "Mesh is a self-intersecting triangle mesh, perform autorefinement...\n";
+
+    // use a polygon soup as container as the output will most likely be non-manifold
+    std::vector<K::Point_3> points;
+    std::vector<std::vector<std::size_t>> polygons;
+    PMP::polygon_mesh_to_polygon_soup(mesh, points, polygons);
+    PMP::autorefine_triangle_soup(points, polygons);
+    std::cout << "Number of facets after preprocessing: "
+              << polygons.size() << "\n";
+
+    if(PMP::does_polygon_soup_self_intersect(points, polygons))
+    {
+      std::cerr << "Error: mesh still self-intersects after autorefine\n";
+      return EXIT_FAILURE;
+    }
+
+    ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(points, polygons);
+
+  }
+  else
+  {
+    std::cout << "Number of facets after preprocessing: "
+              << mesh.number_of_faces() << "\n";
+
+    ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(mesh);
+  }
+
+  std::cout << "Number of constrained facets in the CDT: "
+            << ccdt.number_of_constrained_facets() << '\n';
+
+  std::ofstream ofs(argc > 2 ? argv[2] : "out.mesh");
+  ofs.precision(17);
+  CGAL::IO::write_MEDIT(ofs, ccdt);
+
+  CGAL::draw(ccdt);
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/examples/Constrained_triangulation_3/conforming_constrained_Delaunay_triangulation_3.cpp

@@ -0,0 +1,49 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+#include <CGAL/IO/write_MEDIT.h>
+
+#include <cassert>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+
+int main(int argc, char* argv[])
+{
+  auto filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mpi.off");
+
+  CGAL::Surface_mesh<K::Point_3> mesh;
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh)) {
+    std::cerr << "Error: cannot read file " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "Read " << mesh.number_of_vertices() << " vertices and "
+            << mesh.number_of_faces() << " faces" << std::endl;
+
+  auto ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(mesh);
+
+  //! [use of ccdt.triangulation()]
+  std::cout << "Number of vertices in the CDT: "
+            << ccdt.triangulation().number_of_vertices() << '\n';
+  //! [use of ccdt.triangulation()]
+  std::cout << "Number of constrained facets in the CDT: "
+            << ccdt.number_of_constrained_facets() << '\n';
+
+  std::ofstream ofs(argc > 2 ? argv[2] : "out.mesh");
+  ofs.precision(17);
+  CGAL::IO::write_MEDIT(ofs, ccdt);
+
+  //! [move ccdt to tr]
+  auto tr = std::move(ccdt).triangulation();
+  // Now `tr` is a valid `CGAL::Triangulation_3` object that can be used for further processing.
+  // and the triangulation of `ccdt` is empty.
+  std::cout << "Number of vertices in the triangulation `tr`: "
+            << tr.number_of_vertices() << '\n';
+  std::cout << "Number of vertices in `ccdt`: "
+            << ccdt.triangulation().number_of_vertices() << '\n';
+  assert(ccdt.triangulation().number_of_vertices() == 0);
+  //! [move ccdt to tr]
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/examples/Constrained_triangulation_3/conforming_constrained_Delaunay_triangulation_3_fimap.cpp

@@ -0,0 +1,60 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/remesh_planar_patches.h>
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+
+#include <CGAL/IO/write_MEDIT.h>
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+
+using Mesh = CGAL::Surface_mesh<K::Point_3>;
+using face_descriptor = boost::graph_traits<Mesh>::face_descriptor;
+
+
+int main(int argc, char* argv[])
+{
+  std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross.off");
+
+  Mesh input;
+  if(!CGAL::Polygon_mesh_processing::IO::read_polygon_mesh(filename, input)) {
+    std::cerr << "Invalid input." << std::endl;
+    return 1;
+  }
+
+  std::cout << "Read " << input.number_of_vertices() << " vertices and "
+                       << input.number_of_faces() << " facets\n";
+
+  Mesh mesh;
+  auto plc_facet_map = get(CGAL::face_patch_id_t<int>(), mesh);
+
+  // Remesh planar patches and segment the mesh into planar patches
+  CGAL::Polygon_mesh_processing::remesh_planar_patches(input, mesh,
+                             CGAL::parameters::default_values(),
+                             CGAL::parameters::face_patch_map(plc_facet_map)
+                                              .do_not_triangulate_faces(true));
+
+
+  filename = argc > 2 ? argv[2] : "mesh.ply";
+  CGAL::IO::write_polygon_mesh(filename, mesh,
+      CGAL::parameters::stream_precision(17));
+  std::cout << "Wrote segmented mesh to " << filename << "\n";
+
+  // Create a conforming constrained Delaunay triangulation from the mesh
+  auto ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3(mesh,
+                CGAL::parameters::plc_face_id(plc_facet_map));
+
+  std::cout << "Number of vertices in the CDT: "
+            << ccdt.triangulation().number_of_vertices() << '\n'
+            << "Number of constrained facets in the CDT: "
+            << ccdt.number_of_constrained_facets() << '\n';
+
+  filename = argc > 3 ? argv[3] : "out.mesh";
+  std::ofstream out(filename);
+  out.precision(17);
+  CGAL::IO::write_MEDIT(out, ccdt, CGAL::parameters::with_plc_face_id(true));
+  std::cout << "Wrote CDT to " << filename << "\n";
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/examples/Constrained_triangulation_3/remesh_constrained_Delaunay_triangulation_3.cpp

@@ -0,0 +1,79 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Conforming_constrained_Delaunay_triangulation_cell_base_3.h>
+#include <CGAL/Conforming_constrained_Delaunay_triangulation_vertex_base_3.h>
+#include <CGAL/Tetrahedral_remeshing/Remeshing_cell_base_3.h>
+#include <CGAL/Tetrahedral_remeshing/Remeshing_vertex_base_3.h>
+
+#include <CGAL/make_conforming_constrained_Delaunay_triangulation_3.h>
+#include <CGAL/tetrahedral_remeshing.h>
+
+#include <CGAL/IO/File_medit.h>
+#include <CGAL/IO/polygon_mesh_io.h>
+#include <CGAL/IO/write_MEDIT.h>
+#include <CGAL/property_map.h>
+
+#include <unordered_set>
+#include <fstream>
+#include <string>
+
+using K    = CGAL::Exact_predicates_inexact_constructions_kernel;
+
+using Vbb  = CGAL::Tetrahedral_remeshing::Remeshing_vertex_base_3<K>;
+using Vb   = CGAL::Conforming_constrained_Delaunay_triangulation_vertex_base_3<K, Vbb>;
+
+using Cbb  = CGAL::Tetrahedral_remeshing::Remeshing_cell_base_3<K>;
+using Cb   = CGAL::Conforming_constrained_Delaunay_triangulation_cell_base_3<K, Cbb>;
+
+using Tds  = CGAL::Triangulation_data_structure_3<Vb, Cb>;
+using Tr   = CGAL::Triangulation_3<K, Tds>;
+using CCDT = CGAL::Conforming_constrained_Delaunay_triangulation_3<K, Tr>;
+
+// Triangulation for Remeshing
+using CCDT_Tr = CCDT::Triangulation;
+using Triangulation_3 = CGAL::Triangulation_3<K, CCDT_Tr::Triangulation_data_structure>;
+
+using Vertex_handle = Triangulation_3::Vertex_handle;
+using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
+using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
+using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
+
+
+int main(int argc, char* argv[])
+{
+  std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mpi.off");
+
+  CGAL::Surface_mesh<K::Point_3> mesh;
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh))
+  {
+    std::cerr << "Error: cannot read file " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+  CCDT ccdt = CGAL::make_conforming_constrained_Delaunay_triangulation_3<CCDT>(mesh);
+
+  std::ofstream out("ccdt.mesh");
+  CGAL::IO::write_MEDIT(out, ccdt);
+  out.close();
+
+  Constraints_set constraints;
+  Constraints_pmap constraints_pmap(constraints);
+
+  namespace np = CGAL::parameters;
+  namespace Tet_remesh = CGAL::Tetrahedral_remeshing;
+  Tr tr = Tet_remesh::get_remeshing_triangulation(std::move(ccdt),
+                                                  np::edge_is_constrained_map(constraints_pmap));
+  std::cout << "Number of vertices in tr: " << tr.number_of_vertices() << std::endl;
+
+  CGAL::tetrahedral_isotropic_remeshing(tr,
+                                        1., // target edge length
+                                        np::number_of_iterations(3)
+                                        .edge_is_constrained_map(constraints_pmap));
+
+  std::cout << "Number of vertices in tr: "
+            << tr.number_of_vertices() << std::endl;
+
+  std::ofstream ofs("tr.mesh");
+  CGAL::IO::write_MEDIT(ofs, tr);
+
+  return EXIT_SUCCESS;
+}

Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h

@@ -0,0 +1,993 @@
+// Copyright (c) 2019-2024  GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Laurent Rineau
+
+#ifndef CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H
+#define CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H
+
+#include <CGAL/license/Constrained_triangulation_3.h>
+
+#include <CGAL/Constrained_triangulation_3/internal/config.h>
+
+#include <CGAL/Conforming_constrained_Delaunay_triangulation_vertex_data_3.h>
+#include <CGAL/Triangulation_2/internal/Polyline_constraint_hierarchy_2.h>
+#include <CGAL/Triangulation_segment_traverser_3.h>
+#include <CGAL/unordered_flat_set.h>
+
+#include <CGAL/Mesh_3/io_signature.h>
+#include <CGAL/IO/File_binary_mesh_3.h>
+
+#include <boost/container/flat_set.hpp>
+#include <boost/container/map.hpp>
+#include <boost/container/small_vector.hpp>
+#include <boost/iterator/function_output_iterator.hpp>
+
+#include <fstream>
+#include <bitset>
+
+#ifndef DOXYGEN_RUNNING
+
+namespace CGAL {
+
+template <typename T_3>
+class Conforming_Delaunay_triangulation_3 : public T_3 {
+public:
+  static constexpr bool t_3_is_not_movable =
+      CGAL::cdt_3_msvc_2019_or_older() || false == CGAL::is_nothrow_movable_v<T_3>;
+  using Geom_traits = typename T_3::Geom_traits;
+  using Vertex_handle = typename T_3::Vertex_handle;
+  using Edge = typename T_3::Edge;
+  using Facet = typename T_3::Facet;
+  using Cell_handle = typename T_3::Cell_handle;
+  using Point = typename T_3::Point;
+  using Line = typename T_3::Geom_traits::Line_3;
+  using Locate_type = typename T_3::Locate_type;
+
+  inline static With_offset_tag with_offset{};
+  inline static With_point_tag with_point{};
+  inline static With_point_and_info_tag with_point_and_info{};
+
+  Conforming_Delaunay_triangulation_3(const Geom_traits& gt = Geom_traits())
+    : T_3(gt)
+  {
+  }
+
+protected:
+  struct Compare_vertex_handle {
+    const T_3* tr;
+    Compare_vertex_handle(const T_3* tr) : tr(tr) {}
+    bool operator()(const Vertex_handle va, const Vertex_handle vb) const {
+      return va < vb;
+    }
+  };
+
+  using Constraint_hierarchy =
+      Polyline_constraint_hierarchy_2<Vertex_handle, Compare_vertex_handle, Point>;
+public:
+  using Constrained_polyline_id = typename Constraint_hierarchy::Constraint_id;
+
+protected:
+  using Subconstraint = typename Constraint_hierarchy::Subconstraint;
+
+  auto display_vert(Vertex_handle v) const{
+    std::stringstream os;
+    os.precision(17);
+    os << IO::oformat(v, with_point);
+    return os.str();
+  }
+
+  auto display_subcstr(Subconstraint subconstraint) const {
+    auto [va, vb] = subconstraint;
+    std::stringstream os;
+    os << "(" << IO::oformat(va, with_offset) << ", " << IO::oformat(vb, with_offset) << ")"
+       << ": [ " << display_vert(va) << " - " << display_vert(vb) << " ]";
+    return os.str();
+  }
+
+  class Insert_in_conflict_visitor
+  {
+    Conforming_Delaunay_triangulation_3<T_3>* self;
+  public:
+    Insert_in_conflict_visitor(Conforming_Delaunay_triangulation_3* self) : self(self) {}
+
+    template <class InputIterator>
+    void process_cells_in_conflict(InputIterator cell_it, InputIterator end) {
+      auto d = self->tr().dimension();
+      for( ; cell_it != end; ++cell_it )
+        for( int i = 0; i < d; ++i )
+          for( int j = i+1; j <= d; ++j ) {
+            auto v1 = (*cell_it)->vertex(i);
+            auto v2 = (*cell_it)->vertex(j);
+            if(self->tr().is_infinite(v1) || self->tr().is_infinite(v2)) continue;
+            if(self->use_finite_edges_map()) {
+              if(v1 > v2) std::swap(v1, v2);
+              auto v1_index = v1->time_stamp();
+              [[maybe_unused]] auto nb_erased = self->all_finite_edges[v1_index].erase(v2);
+              if constexpr (cdt_3_can_use_cxx20_format()) if(self->debug_finite_edges_map() && nb_erased > 0) {
+                std::cerr << cdt_3_format("erasing edge {} {}\n", self->display_vert((std::min)(v1, v2)),
+                                        self->display_vert((std::max)(v1, v2)));
+              }
+            }
+            auto [contexts_begin, contexts_end] =
+              self->constraint_hierarchy.contexts(v1, v2);
+            if(contexts_begin != contexts_end) {
+              self->add_to_subconstraints_to_conform(v1, v2, contexts_begin->id());
+            }
+          }
+    }
+
+    void after_insertion(Vertex_handle v) const {
+      if(!self->use_finite_edges_map()) return;
+      CGAL_assertion(self->dimension() > 1);
+      self->incident_edges(v, boost::make_function_output_iterator([this](Edge e) { self->new_edge(e); }));
+      self->incident_cells(v, boost::make_function_output_iterator([this, v](Cell_handle c) {
+        auto v_index = c->index(v);
+        if(self->dimension() == 2) {
+          auto j = self->cw(v_index);
+          auto k = self->ccw(v_index);
+          self->new_edge(Edge(c, j, k));
+        } else {
+          for(int i = 0; i < 3; ++i) {
+            auto j = self->vertex_triple_index(v_index, i);
+            auto k = self->vertex_triple_index(v_index, self->cw(i));
+            self->new_edge(Edge(c, j, k));
+          }
+        }
+      }));
+    }
+
+    void reinsert_vertices(Vertex_handle v) const { after_insertion(v); }
+    Vertex_handle replace_vertex(Cell_handle c, int index, const Point& ) const
+    {
+      return c->vertex(index);
+    }
+    void hide_point(Cell_handle, const Point& ) const {}
+
+    void insert_Steiner_point_on_constraint([[maybe_unused]] Constrained_polyline_id constraint,
+                                            [[maybe_unused]] Vertex_handle va,
+                                            [[maybe_unused]] Vertex_handle vb,
+                                            [[maybe_unused]] Vertex_handle v_Steiner) const
+    {
+    }
+
+    Vertex_handle insert_in_triangulation(const Point& p, Locate_type lt, Cell_handle c, int li, int lj) {
+      return self->insert_impl_do_not_split(p, lt, c, li, lj, *this);
+    }
+  };
+
+  auto make_subconstraint(Vertex_handle va, Vertex_handle vb) {
+    return make_sorted_pair(va, vb);
+  }
+
+  void add_to_subconstraints_to_conform(Vertex_handle va, Vertex_handle vb,
+                                        Constrained_polyline_id id) {
+    const auto pair = make_subconstraint(va, vb);
+#if CGAL_DEBUG_CDT_3 & 32
+    std::cerr << "tr().subconstraints_to_conform.push("
+              << display_subcstr(pair) << ")\n";
+#endif // CGAL_DEBUG_CDT_3
+    subconstraints_to_conform.push({pair, id});
+  }
+
+  template <typename Visitor>
+  Constrained_polyline_id insert_constrained_edge_impl(Vertex_handle va, Vertex_handle vb,
+                                             Visitor&) {
+    if(va != vb) {
+      if(segment_vertex_epsilon != 0.) {
+        auto [min_dist, min_vertex] = min_distance_and_vertex_between_constraint_and_encroaching_vertex(va, vb);
+        check_segment_vertex_distance_or_throw(va, vb, min_vertex, CGAL::to_double(min_dist),
+                                               Check_distance::NON_SQUARED_DISTANCE);
+      }
+      const Constrained_polyline_id c_id = constraint_hierarchy.insert_constraint(va, vb);
+      pair_of_vertices_to_cid.emplace(make_sorted_pair(va, vb), c_id);
+      // traverse all the vertices along [va, vb] and add pairs of consecutive
+      // vertices as sub-constraints.
+      std::for_each(tr().segment_traverser_simplices_begin(va, vb), tr().segment_traverser_simplices_end(),
+                    [&, prev = Vertex_handle{}](auto simplex) mutable {
+                      // std::cerr << "- " << oformat(simplex, With_point_tag{}) << '\n';
+                      if(simplex.dimension() == 0) {
+                        const auto v = static_cast<Vertex_handle>(simplex);
+                        v->ccdt_3_data().set_on_constraint(c_id);
+                        if(prev != Vertex_handle{}) {
+                          if(v != vb) {
+                            v->ccdt_3_data().set_vertex_type(CDT_3_vertex_type::STEINER_ON_EDGE);
+                            constraint_hierarchy.add_Steiner(prev, vb, v);
+                          }
+                          add_to_subconstraints_to_conform(prev, v, c_id);
+                        }
+                        prev = v;
+                      }
+                    });
+      return c_id;
+    }
+    return {};
+  }
+
+  void new_edge(Edge e)
+  {
+    if(!update_all_finite_edges_) return;
+    auto [v1, v2] = make_sorted_pair(tr().vertices(e));
+    if(tr().is_infinite(v1) || tr().is_infinite(v2))
+      return;
+    [[maybe_unused]] auto [_, inserted] = all_finite_edges[v1->time_stamp()].insert(v2);
+    if constexpr (cdt_3_can_use_cxx20_format()) if (debug_finite_edges_map() && inserted) {
+      if(v2 < v1) std::swap(v1, v2);
+      std::cerr << cdt_3_format("new_edge({}, {})\n", display_vert(v1), display_vert(v2));
+    }
+  }
+
+  void new_cell(Cell_handle c) {
+    auto d = tr().dimension();
+    for(int i = 0; i < d; ++i) {
+      for(int j = i+1; j <= d; ++j) {
+        new_edge(Edge(c, i, j));
+      }
+    }
+  }
+
+  auto new_cells_output_iterator()
+  {
+    return boost::function_output_iterator([this](Cell_handle c) {
+      if(use_finite_edges_map())
+        this->new_cell(c);
+    });
+  }
+
+  template <typename Visitor>
+  Vertex_handle insert_impl_do_not_split(const Point &p, Locate_type lt, Cell_handle c,
+                                         int li, int lj, Visitor& visitor)
+  {
+    switch (tr().dimension()) {
+    case 3: {
+      typename T_3::Conflict_tester_3 tester(p, this);
+      auto v = tr().insert_in_conflict(p, lt, c, li, lj, tester,
+                                     visitor);
+      new_vertex(v);
+      return v;
+    } // dim 3
+    case 2: {
+      typename T_3::Conflict_tester_2 tester(p, this);
+      auto v = tr().insert_in_conflict(p, lt, c, li, lj, tester, visitor);
+      if(use_finite_edges_map()) {
+        new_vertex(v);
+        tr().incident_edges(v, boost::make_function_output_iterator([&](Edge e) { this->new_edge(e); }));
+      }
+      return v;
+    } // dim 2
+    default:
+      // dimension <= 1
+      // Do not use the generic insert.
+      auto v = tr().insert(p, c);
+      if(use_finite_edges_map()) {
+        new_vertex(v);
+        all_finite_edges.clear();
+        if (debug_finite_edges_map()) std::cerr << "all_finite_edges.clear()\n";
+        for(auto e: tr().all_edges()) {
+          new_edge(e);
+        }
+      }
+      return v;
+    }
+  }
+
+  template <typename Visitor>
+  Vertex_handle insert_impl(const Point &p, Locate_type lt, Cell_handle c,
+                            int li, int lj, Visitor& visitor)
+  {
+    Vertex_handle v1, v2;
+    bool split_constrained_edge = false;
+    if(lt == T_3::EDGE) {
+      v1 = c->vertex(li);
+      v2 = c->vertex(lj);
+      if(constraint_hierarchy.is_subconstraint(v1, v2)) {
+        split_constrained_edge = true;
+      }
+    }
+    Vertex_handle new_vertex = insert_impl_do_not_split(p, lt, c, li, lj, visitor);
+    if(split_constrained_edge) {
+      constraint_hierarchy.split_constraint(v1, v2, new_vertex);
+    }
+    return new_vertex;
+  }
+
+  void update_bbox(const Point& p) {
+    bbox = bbox + p.bbox();
+    if(max_bbox_edge_length) {
+      update_max_bbox_edge_length();
+    }
+  }
+
+  void update_max_bbox_edge_length() {
+    double d_x = bbox.xmax() - bbox.xmin();
+    double d_y = bbox.ymax() - bbox.ymin();
+    double d_z = bbox.zmax() - bbox.zmin();
+
+    max_bbox_edge_length = (std::max)(d_x, (std::max)(d_y, d_z));
+  }
+
+public:
+  void set_segment_vertex_epsilon(double epsilon) {
+    segment_vertex_epsilon = epsilon;
+  }
+
+  bool debug_Steiner_points() const {
+    return debug_flags[static_cast<int>(Debug_flags::Steiner_points)];
+  }
+
+  void debug_Steiner_points(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::Steiner_points), b);
+  }
+
+  bool debug_input_faces() const {
+    return debug_flags[static_cast<int>(Debug_flags::input_faces)];
+  }
+
+  void debug_input_faces(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::input_faces), b);
+  }
+
+  bool debug_missing_region() const {
+    return debug_flags[static_cast<int>(Debug_flags::missing_region)];
+  }
+
+  void debug_missing_region(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::missing_region), b);
+  }
+
+  bool debug_regions() const {
+    return debug_flags[static_cast<int>(Debug_flags::regions)];
+  }
+
+  void debug_regions(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::regions), b);
+  }
+
+  bool debug_copy_triangulation_into_hole() const {
+    return debug_flags[static_cast<int>(Debug_flags::copy_triangulation_into_hole)];
+  }
+
+  void debug_copy_triangulation_into_hole(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::copy_triangulation_into_hole), b);
+  }
+
+  bool debug_validity() const {
+    return debug_flags[static_cast<int>(Debug_flags::validity)];
+  }
+
+  void debug_validity(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::validity), b);
+  }
+
+  bool use_older_cavity_algorithm() const {
+    return debug_flags[static_cast<int>(Debug_flags::use_older_cavity_algorithm)];
+  }
+
+  void use_older_cavity_algorithm(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::use_older_cavity_algorithm), b);
+  }
+
+  bool debug_finite_edges_map() const {
+    return debug_flags[static_cast<int>(Debug_flags::debug_finite_edges_map)];
+  }
+
+  void debug_finite_edges_map(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::debug_finite_edges_map), b);
+  }
+
+  bool use_finite_edges_map() const {
+    return update_all_finite_edges_ && debug_flags[static_cast<int>(Debug_flags::use_finite_edges_map)];
+  }
+
+  void use_finite_edges_map(bool b) {
+    debug_flags.set(static_cast<int>(Debug_flags::use_finite_edges_map), b);
+  }
+
+  Vertex_handle insert(const Point &p, Locate_type lt, Cell_handle c,
+                       int li, int lj)
+  {
+    update_bbox(p);
+    auto v = insert_impl(p, lt, c, li, lj, insert_in_conflict_visitor);
+    restore_Delaunay(insert_in_conflict_visitor);
+    return v;
+  }
+
+  Vertex_handle insert(const Point &p, Cell_handle start = {}) {
+    Locate_type lt;
+    int li, lj;
+
+    Cell_handle c = tr().locate(p, lt, li, lj, start);
+    return insert(p, lt, c, li, lj);
+  }
+
+  Constrained_polyline_id insert_constrained_edge(Vertex_handle va, Vertex_handle vb)
+  {
+    const auto id = insert_constrained_edge_impl(va, vb, insert_in_conflict_visitor);
+    restore_Delaunay(insert_in_conflict_visitor);
+    return id;
+  }
+
+  bool is_edge(Vertex_handle va, Vertex_handle vb) const {
+    const bool is_edge_v1 =
+        ((debug_finite_edges_map() && use_finite_edges_map()) || !use_finite_edges_map()) && tr().tds().is_edge(va, vb);
+
+    if(use_finite_edges_map() && va > vb) std::swap(va, vb);
+    const auto va_index = va->time_stamp();
+    const bool is_edge_v2 =
+        use_finite_edges_map() && all_finite_edges[va_index].find(vb) != all_finite_edges[va_index].end();
+
+    if(debug_finite_edges_map() && use_finite_edges_map() && is_edge_v1 != is_edge_v2) {
+      std::cerr << "!! Inconsistent edge status\n";
+      std::cerr << "  -> constraint " << display_vert(va) << "     " << display_vert(vb) << '\n';
+      std::cerr << "  ->     edge " << (is_edge_v1 ? "is" : "is not") << " in the triangulation\n";
+      std::cerr << "  ->     edge " << (is_edge_v2 ? "is" : "is not") << " in all_finite_edges\n";
+      debug_dump("bug-inconsistent-edge-status");
+      CGAL_error();
+    }
+    const bool is_edge = use_finite_edges_map() ? is_edge_v2 : is_edge_v1;
+    return is_edge;
+  }
+
+  using T_3::is_edge;
+
+  bool is_conforming() const {
+    return std::all_of(constraint_hierarchy.subconstraints_begin(),
+                       constraint_hierarchy.subconstraints_end(),
+                       [this](const auto &sc) {
+                         const auto [va, vb] = sc;
+                         const auto is_edge = this->is_edge(va, vb);
+#if CGAL_DEBUG_CDT_3 & 128 && CGAL_CAN_USE_CXX20_FORMAT
+                         std::cerr << cdt_3_format("is_conforming>> Edge is 3D: {}  ({} , {})\n",
+                                                  is_edge,
+                                                  CGAL::IO::oformat(va, with_point_and_info),
+                                                  CGAL::IO::oformat(vb, with_point_and_info));
+#endif // CGAL_DEBUG_CDT_3
+                         return is_edge;
+                       });
+  }
+
+  enum class Check_distance { SQUARED_DISTANCE, NON_SQUARED_DISTANCE };
+
+  void check_segment_vertex_distance_or_throw(Vertex_handle va,
+                                              Vertex_handle vb,
+                                              Vertex_handle min_vertex,
+                                              double min_dist,
+                                              Check_distance option)
+  {
+    if(!max_bbox_edge_length) {
+      update_max_bbox_edge_length();
+    }
+    if((option == Check_distance::NON_SQUARED_DISTANCE && min_dist < segment_vertex_epsilon * *max_bbox_edge_length) ||
+       (option == Check_distance::SQUARED_DISTANCE &&
+        min_dist < CGAL::square(segment_vertex_epsilon * *max_bbox_edge_length)))
+    {
+      std::stringstream ss;
+      ss.precision(std::cerr.precision());
+      ss << "A constrained segment is too close to a vertex.\n";
+      ss << "  -> vertex " << display_vert(min_vertex) << '\n';
+      ss << "  -> constrained segment " << display_vert(va) << "  -  " << display_vert(vb) << '\n';
+      ss << "  -> distance = " << min_dist << '\n';
+      ss << "  -> max_bbox_edge_length = " << *max_bbox_edge_length << '\n';
+      CGAL_error_msg(ss.str().c_str());
+    }
+  }
+
+  auto ancestors_of_Steiner_vertex_on_edge(Vertex_handle v) const {
+    std::pair<Vertex_handle, Vertex_handle> result;
+    CGAL_precondition(v->ccdt_3_data().is_Steiner_vertex_on_edge());
+    CGAL_assertion(v->ccdt_3_data().number_of_incident_constraints() == 1);
+    const auto v_time_stamp = v->time_stamp();
+    const auto constrained_polyline_id = v->ccdt_3_data().constrained_polyline_id(*this);
+    const auto first = this->constraint_hierarchy.vertices_in_constraint_begin(constrained_polyline_id);
+    const auto end =   this->constraint_hierarchy.  vertices_in_constraint_end(constrained_polyline_id);
+    std::cerr << "ancestors_of_Steiner_vertex_on_edge " << display_vert(v) << '\n';
+    for(auto it = first; it != end; ++it) {
+      std::cerr << "  - " << display_vert(*it) << '\n';
+    }
+    CGAL_assertion(first != end);
+    const auto last = std::prev(this->constraint_hierarchy.vertices_in_constraint_end(constrained_polyline_id));
+    CGAL_assertion(first != last);
+    CGAL_assertion((*first)->time_stamp() < v_time_stamp);
+    CGAL_assertion((*last)->time_stamp() < v_time_stamp);
+    for(auto it = first; (*it) != v; ++it) {
+      if((*it)->time_stamp() < v_time_stamp) {
+        result.first = *it;
+      }
+      CGAL_assertion(it != last);
+    }
+    for(auto it = last; (*it) != v; --it) {
+      if((*it)->time_stamp() < v_time_stamp) {
+        result.second = *it;
+      }
+      CGAL_assertion(it != first);
+    }
+    return result;
+  }
+
+  auto min_distance_and_vertex_between_constraint_and_encroaching_vertex(Vertex_handle va, Vertex_handle vb) const {
+    struct Result {
+      typename T_3::Geom_traits::FT min_dist = (std::numeric_limits<double>::max)();
+      Vertex_handle v = {};
+    } result;
+    const auto vector_of_encroaching_vertices = encroaching_vertices(va, vb);
+    for(auto v: vector_of_encroaching_vertices) {
+      const auto dist = CGAL::approximate_sqrt(squared_distance(tr().point(v), Line{tr().point(va), tr().point(vb)}));
+      if(dist < result.min_dist) {
+        result = Result{dist, v};
+      }
+    }
+    return result;
+  }
+
+  bool write_missing_segments_file(std::ostream &out) {
+    bool any_missing_segment = false;
+    std::for_each(
+        constraint_hierarchy.subconstraints_begin(), constraint_hierarchy.subconstraints_end(),
+        [this, &out, &any_missing_segment](const auto &sc) {
+          const auto [v0, v1] = sc;
+          if (!this->is_edge(v0, v1)) {
+            out << "2 " << this->tr().point(v0) << " " << this->tr().point(v1)
+                << '\n';
+            any_missing_segment = true;
+          }
+        });
+    return any_missing_segment;
+  }
+
+  void write_all_segments_file(std::ostream &out) {
+    std::for_each(
+        constraint_hierarchy.subconstraints_begin(), constraint_hierarchy.subconstraints_end(),
+        [this, &out](const auto &sc) {
+          const auto [v0, v1] = sc;
+          out << "2 " << this->tr().point(v0) << " " << this->tr().point(v1) << '\n';
+        });
+  }
+
+  /// @{
+  /// remove functions cannot be called
+  void remove(Vertex_handle) = delete;
+  void remove_cluster() = delete;
+  /// @}
+
+  static std::string io_signature() {
+    return Get_io_signature<T_3>()();
+  }
+protected:
+  void debug_dump(std::string filename) const {
+    {
+      std::ofstream dump(filename + ".binary.cgal", std::ios::binary);
+      CGAL::IO::save_binary_file(dump, *this);
+    }
+    {
+      std::ofstream dump(filename + "_point.xyz");
+      dump.precision(17);
+      for(auto vh: this->finite_vertex_handles()){
+        dump << this->point(vh) << '\n';
+      }
+    }
+  }
+
+  template <typename Visitor>
+  void restore_Delaunay(Visitor& visitor) {
+    update_all_finite_edges();
+    while(!subconstraints_to_conform.empty()) {
+      const auto [subconstraint, constrained_polyline_id] = subconstraints_to_conform.top();
+      subconstraints_to_conform.pop();
+      const auto [va, vb] = subconstraint;
+      if(!constraint_hierarchy.is_subconstraint(va, vb)) {
+        continue;
+      }
+#if CGAL_DEBUG_CDT_3 & 32
+      std::cerr << "tr().subconstraints_to_conform.pop()="
+                << display_subcstr(subconstraint) << "\n";
+#endif // CGAL_DEBUG_CDT_3
+      conform_subconstraint(subconstraint, constrained_polyline_id, visitor);
+    }
+  }
+
+  template <typename Visitor>
+  Vertex_handle insert_Steiner_point_on_subconstraint(
+      Point steiner_pt, Cell_handle hint,
+      Subconstraint subconstraint, Constrained_polyline_id constraint, Visitor& visitor)
+  {
+    const Vertex_handle va = subconstraint.first;
+    const Vertex_handle vb = subconstraint.second;
+    Locate_type lt;
+    int li, lj;
+    const Cell_handle c = tr().locate(steiner_pt, lt, li, lj, hint);
+    const Vertex_handle v = visitor.insert_in_triangulation(steiner_pt, lt, c, li, lj);
+    v->ccdt_3_data().set_vertex_type(CDT_3_vertex_type::STEINER_ON_EDGE);
+    if(lt != T_3::VERTEX) {
+      v->ccdt_3_data().set_on_constraint(constraint);
+    }
+    constraint_hierarchy.add_Steiner(va, vb, v);
+    visitor.insert_Steiner_point_on_constraint(constraint, va, vb, v);
+    add_to_subconstraints_to_conform(va, v, constraint);
+    add_to_subconstraints_to_conform(v, vb, constraint);
+
+    new_vertex(v);
+
+    return v;
+  }
+
+  /// Return `true` if a Steiner point was inserted
+  template <typename Visitor>
+  bool conform_subconstraint(Subconstraint subconstraint,
+                             Constrained_polyline_id constraint,
+                             Visitor& visitor)
+  {
+    const Vertex_handle va = subconstraint.first;
+    const Vertex_handle vb = subconstraint.second;
+    CGAL_assertion(va != vb);
+    if(!this->is_edge(va, vb)) {
+      const auto& [steiner_pt, hint, ref_vertex] = construct_Steiner_point(constraint, subconstraint);
+      [[maybe_unused]] const auto v =
+          insert_Steiner_point_on_subconstraint(steiner_pt, hint, subconstraint, constraint, visitor);
+      if(debug_Steiner_points()) {
+        const auto [c_start, c_end] = constraint_extremities(constraint);
+        std::cerr << "(" << IO::oformat(va, with_offset) << ", " << IO::oformat(vb, with_offset) << ")";
+        std::cerr << ": [ " << display_vert(c_start) << " - " << display_vert(c_end) << " ] ";
+        std::cerr << "  new vertex " << display_vert(v) << '\n';
+      }
+      return true;
+    }
+
+    return false;
+  }
+
+  Constrained_polyline_id constraint_from_extremities(Vertex_handle va, Vertex_handle vb) const {
+    if(va->ccdt_3_data().number_of_incident_constraints() == 0 || vb->ccdt_3_data().number_of_incident_constraints() == 0)
+    {
+      return {};
+    }
+    auto it = pair_of_vertices_to_cid.find(make_sorted_pair(va, vb));
+    if(it != pair_of_vertices_to_cid.end()) {
+      return it->second;
+    }
+    return {};
+    // @TODO: cleanup the rest of the function, and `constraint_around`
+    Constrained_polyline_id c_id = constraint_around(va, vb, false);
+    if(c_id != Constrained_polyline_id{}) return c_id;
+    c_id = constraint_around(vb, va, false);
+    if(c_id != Constrained_polyline_id{}) return c_id;
+    c_id = constraint_around(va, vb, true);
+    return c_id;
+  }
+
+  auto constraint_extremities(Constrained_polyline_id c_id) const {
+      CGAL_assertion(std::find(this->constraint_hierarchy.constraints_begin(),
+                               this->constraint_hierarchy.constraints_end(), c_id) != this->constraint_hierarchy.constraints_end());
+      CGAL_assertion(this->constraint_hierarchy.vertices_in_constraint_begin(c_id) !=
+                     this->constraint_hierarchy.vertices_in_constraint_end(c_id));
+#if CGAL_DEBUG_CDT_3 & 8
+      std::cerr << "constraint " << (void*) c_id.vl_ptr() << " has "
+                << c_id.vl_ptr()->skip_size() << " vertices\n";
+#endif // CGAL_DEBUG_CDT_3
+      const auto begin = this->constraint_hierarchy.vertices_in_constraint_begin(c_id);
+      const auto end = this->constraint_hierarchy.vertices_in_constraint_end(c_id);
+      const auto c_va = *begin;
+      const auto c_vb = *std::prev(end);
+    return std::make_pair(c_va, c_vb);
+  }
+
+  Constrained_polyline_id constraint_around(Vertex_handle va, Vertex_handle vb, bool expensive = true) const {
+    auto constraint_id_goes_to_vb = [this, va, vb](Constrained_polyline_id c_id) {
+      const auto [c_va, c_vb] = constraint_extremities(c_id);
+      if (va == c_va && vb == c_vb)
+        return true;
+      if (vb == c_va && va == c_vb)
+        return true;
+      return false;
+    }; // end lambda constraint_id_goes_to_vb
+    if (va->ccdt_3_data().number_of_incident_constraints() == 1)
+    {
+      const Constrained_polyline_id c_id = va->ccdt_3_data().constrained_polyline_id(*this);
+      CGAL_assertion(c_id != Constrained_polyline_id{});
+      if(constraint_id_goes_to_vb(c_id)) return c_id;
+    } else if (expensive == true && va->ccdt_3_data().number_of_incident_constraints() > 1) {
+      boost::container::small_vector<Vertex_handle, 64> adj_vertices;
+      this->finite_adjacent_vertices(va, std::back_inserter(adj_vertices));
+      for(auto other_v: adj_vertices) {
+        for(auto context: this->constraint_hierarchy.contexts(va, other_v)) {
+          const Constrained_polyline_id c_id = context.id();
+          if(constraint_id_goes_to_vb(c_id)) return c_id;
+        }
+      }
+    }
+    return Constrained_polyline_id{};
+  }
+
+  struct Construct_Steiner_point_return_type {
+    typename T_3::Geom_traits::Point_3 point;
+    Cell_handle hint;
+    Vertex_handle reference_vertex;
+  };
+
+  auto encroaching_vertices(Vertex_handle va, Vertex_handle vb) const {
+    auto& gt = tr().geom_traits();
+    auto angle_functor = gt.angle_3_object();
+
+    const auto& pa = tr().point(va);
+    const auto& pb = tr().point(vb);
+
+    namespace bc = boost::container;
+    bc::flat_set<Vertex_handle, std::less<Vertex_handle>,
+                 bc::small_vector<Vertex_handle, 256>>
+        encroaching_vertices;
+    auto register_vertex = [this,&encroaching_vertices](Vertex_handle v) {
+      if(tr().is_infinite(v)) return;
+      // std::cerr << "register_vertex " << display_vert(v) << '\n';
+      encroaching_vertices.insert(v);
+    };
+    auto fill_encroaching_vertices = [&](const auto simplex) {
+#if CGAL_DEBUG_CDT_3 & 0x10
+      std::cerr << " - " << IO::oformat(simplex, With_point_tag{}) << '\n';
+#endif // CGAL_DEBUG_CDT_3
+      auto visit_cell = [&](Cell_handle cell) {
+        for(int i = 0, end = this->tr().dimension() + 1; i < end; ++i) {
+          const auto v = cell->vertex(i);
+          register_vertex(v);
+        }
+      };
+      switch(simplex.dimension()) {
+      case 3: {
+        const auto cell = static_cast<Cell_handle>(simplex);
+        visit_cell(cell);
+      } break;
+      case 2: {
+        const auto [cell, facet_index] = static_cast<Facet>(simplex);
+        visit_cell(cell);
+        if(tr().dimension() > 2) {
+          const auto [other_cell, other_index] = tr().mirror_facet({cell, facet_index});
+          register_vertex(other_cell->vertex(other_index));
+        }
+        break;
+      }
+      case 1: {
+        auto edge = static_cast<Edge>(simplex);
+        if(tr().dimension() < 3) {
+          auto [cell, i, j] = edge;
+          visit_cell(cell);
+          if(tr().dimension() < 2) break;
+          auto neighbor_cell = cell->neighbor(3 - i - j);
+          visit_cell(neighbor_cell);
+          break;
+        }
+        auto circ = tr().incident_cells(edge);
+        CGAL_assertion(circ != nullptr);
+        const auto end = circ;
+        do {
+          visit_cell(circ);
+        } while(++circ != end);
+      } break;
+      case 0: {
+        const auto v = static_cast<Vertex_handle>(simplex);
+        if(v != va && v != vb) {
+          std::cerr << "!! The constraint passes through a vertex!\n";
+          std::cerr << "  -> constraint " << display_vert(va) << "     " << display_vert(vb) << '\n';
+          std::cerr << "  ->     vertex " << display_vert(v) << '\n';
+          debug_dump("bug-through-vertex");
+          CGAL_error();
+        }
+      } break;
+      default: CGAL_unreachable();
+      } // end switch
+    };
+    std::for_each(tr().segment_traverser_simplices_begin(va, vb), tr().segment_traverser_simplices_end(),
+                  fill_encroaching_vertices);
+    auto vector_of_encroaching_vertices = encroaching_vertices.extract_sequence();
+#if CGAL_DEBUG_CDT_3 & 0x10
+    std::cerr << "  -> vector_of_encroaching_vertices (before filter):\n";
+    std::for_each(vector_of_encroaching_vertices.begin(),
+                  vector_of_encroaching_vertices.end(),
+                  [this](Vertex_handle v){
+                    std::cerr << "    " << this->display_vert(v) << '\n';
+                  });
+#endif // CGAL_DEBUG_CDT_3
+    auto end = std::remove_if(vector_of_encroaching_vertices.begin(),
+                              vector_of_encroaching_vertices.end(),
+                              [va, vb, pa, pb, &angle_functor, this](Vertex_handle v) {
+                                if(va == v || vb == v) return true;
+                                return angle_functor(pa,
+                                                    this->tr().point(v),
+                                                    pb) == ACUTE;
+                              });
+#if CGAL_DEBUG_CDT_3 & 0x10
+    std::cerr << "  -> vector_of_encroaching_vertices (after filter):\n";
+    std::for_each(vector_of_encroaching_vertices.begin(), end, [&](Vertex_handle v) {
+      std::cerr << "    " << this->display_vert(v) << "  angle " << approximate_angle(pa, this->tr().point(v), pb)
+                << '\n';
+    });
+#endif // CGAL_DEBUG_CDT_3
+    vector_of_encroaching_vertices.erase(end, vector_of_encroaching_vertices.end());
+    return vector_of_encroaching_vertices;
+  }
+
+  Construct_Steiner_point_return_type
+  construct_Steiner_point(Constrained_polyline_id constrained_polyline_id, Subconstraint subconstraint)
+  {
+    auto& gt = tr().geom_traits();
+    auto compare_angle_functor = gt.compare_angle_3_object();
+    auto vector_functor = gt.construct_vector_3_object();
+    auto midpoint_functor = gt.construct_midpoint_3_object();
+    auto scaled_vector_functor = gt.construct_scaled_vector_3_object();
+    auto sq_length_functor = gt.compute_squared_length_3_object();
+    auto sc_product_functor = gt.compute_scalar_product_3_object();
+    auto translate_functor = gt.construct_translated_point_3_object();
+
+    const Vertex_handle va = subconstraint.first;
+    const Vertex_handle vb = subconstraint.second;
+    const auto& pa = tr().point(va);
+    const auto& pb = tr().point(vb);
+    const auto [orig_va, orig_vb] = constraint_extremities(constrained_polyline_id);
+    const auto& orig_pa = tr().point(orig_va);
+    const auto& orig_pb = tr().point(orig_vb);
+
+    if(this->dimension() < 2) {
+      std::cerr << "dim < 2: midpoint\n";
+      return {midpoint_functor(pa, pb), va->cell(), va};
+    }
+
+#if CGAL_DEBUG_CDT_3 & 0x10
+    std::cerr << "construct_Steiner_point( " << display_vert(va) << " , "
+              << display_vert(vb) << " )\n";
+#endif // CGAL_DEBUG_CDT_3
+
+    const auto vector_of_encroaching_vertices = encroaching_vertices(va, vb);
+    CGAL_assertion(vector_of_encroaching_vertices.size() > 0);
+
+    const auto reference_vertex_it = std::max_element(
+        vector_of_encroaching_vertices.begin(), vector_of_encroaching_vertices.end(),
+        [pa, pb, &compare_angle_functor, this](Vertex_handle v1,
+                                               Vertex_handle v2) {
+          return compare_angle_functor(pa, this->tr().point(v1), pb,
+                                       pa, this->tr().point(v2), pb) == SMALLER;
+        });
+    CGAL_assertion(reference_vertex_it != vector_of_encroaching_vertices.end());
+#if CGAL_CDT_3_DEBUG_CONFORMING
+    std::cerr << "  -> reference point: " << display_vert(*reference_vertex_it)
+              << '\n';
+#endif // CGAL_CDT_3_DEBUG_CONFORMING
+    const auto reference_vertex = *reference_vertex_it;
+    const auto& reference_point = tr().point(reference_vertex);
+
+    const auto vector_ab = vector_functor(pa, pb);
+
+    if(reference_vertex->ccdt_3_data().is_Steiner_vertex_on_edge()) {
+      CGAL_assertion(reference_vertex->ccdt_3_data().number_of_incident_constraints() == 1);
+      const auto ref_constrained_polyline_id = reference_vertex->ccdt_3_data().constrained_polyline_id(*this);
+      const auto [ref_va, ref_vb] = constraint_extremities(ref_constrained_polyline_id);
+#if CGAL_CDT_3_DEBUG_CONFORMING
+      std::cerr << "  reference point is on constraint: " << display_vert(ref_va)
+                << "    " << display_vert(ref_vb) << '\n'
+                << "  original constraint:              " << display_vert(orig_va)
+                << "    " << display_vert(orig_vb) << '\n';
+#endif // CGAL_CDT_3_DEBUG_CONFORMING
+      const auto vector_orig_ab = vector_functor(orig_pa, orig_pb);
+      const auto length_ab = CGAL::approximate_sqrt(sq_length_functor(vector_ab));
+      auto return_orig_result_point =
+          [&](auto lambda, Point orig_pa, Point orig_pb)
+              -> Construct_Steiner_point_return_type
+          {
+            const auto vector_orig_ab = vector_functor(orig_pa, orig_pb);
+            const auto inter_point = translate_functor(orig_pa, scaled_vector_functor(vector_orig_ab, lambda));
+            const auto dist_a_result = CGAL::approximate_sqrt(sq_length_functor(vector_functor(pa, inter_point)));
+            const auto ratio = dist_a_result / length_ab;
+            const auto result_point = (ratio < 0.2 || ratio > 0.8)
+                                          ? midpoint_functor(pa, pb)
+                                          : inter_point;
+
+#if CGAL_CDT_3_DEBUG_CONFORMING
+            std::cerr << "  ref ratio = " << ratio << '\n';
+            std::cerr << "  -> Steiner point: " << result_point << '\n';
+#endif // CGAL_CDT_3_DEBUG_CONFORMING
+            return {result_point, reference_vertex->cell(), reference_vertex};
+          };
+
+      const auto length_orig_ab = CGAL::approximate_sqrt(sq_length_functor(vector_orig_ab));
+      if(ref_va == orig_va || ref_vb == orig_va) {
+        const auto vector_orig_a_ref = vector_functor(orig_pa, reference_point);
+        const auto length_orig_a_ref = CGAL::approximate_sqrt(sq_length_functor(vector_orig_a_ref));
+        const auto lambda = length_orig_a_ref / length_orig_ab;
+        return return_orig_result_point(lambda, orig_pa, orig_pb);
+      } else if(ref_va == orig_vb || ref_vb == orig_vb) {
+        const auto vector_orig_b_ref = vector_functor(orig_pb, reference_point);
+        const auto length_orig_b_ref = CGAL::approximate_sqrt(sq_length_functor(vector_orig_b_ref));
+        const auto lambda = length_orig_b_ref / length_orig_ab;
+        return return_orig_result_point(lambda, orig_pb, orig_pa);
+      }
+    } else {
+      if(segment_vertex_epsilon > 0) {
+        if(!max_bbox_edge_length) {
+          update_max_bbox_edge_length();
+        }
+        auto sq_dist = squared_distance(reference_point, Line{orig_pa, orig_pb});
+        check_segment_vertex_distance_or_throw(orig_va, orig_vb, reference_vertex, CGAL::to_double(sq_dist),
+                                               Check_distance::SQUARED_DISTANCE);
+      }
+    }
+    // compute the projection of the reference point
+    const auto vector_a_ref = vector_functor(pa, reference_point);
+    const auto lambda = sc_product_functor(vector_a_ref, vector_ab) / sq_length_functor(vector_ab);
+    const auto result_point = (lambda < 0.2 || lambda > 0.8)
+                                  ? midpoint_functor(pa, pb)
+                                  : translate_functor(pa, scaled_vector_functor(vector_ab, lambda));
+
+#if CGAL_CDT_3_DEBUG_CONFORMING
+    std::cerr << "  lambda = " << lambda << '\n';
+    std::cerr << "  -> Steiner point: " << result_point << '\n';
+#endif // CGAL_CDT_3_DEBUG_CONFORMING
+    return {result_point, reference_vertex->cell(), reference_vertex};
+  }
+
+protected:
+  T_3& tr() { return *this; };
+  const T_3& tr() const { return *this; };
+
+  Compare_vertex_handle comp = {this};
+  Constraint_hierarchy constraint_hierarchy = {comp};
+  static_assert(CGAL::cdt_3_msvc_2019_or_older() || CGAL::is_nothrow_movable_v<Constraint_hierarchy>);
+  Bbox_3 bbox{};
+  double segment_vertex_epsilon = 1e-8;
+  std::optional<double> max_bbox_edge_length;
+  using Pair_of_vertex_handles = std::pair<Vertex_handle, Vertex_handle>;
+  boost::container::map<Pair_of_vertex_handles, Constrained_polyline_id> pair_of_vertices_to_cid;
+  Insert_in_conflict_visitor insert_in_conflict_visitor = {this};
+
+  using Stack_info = std::pair<Subconstraint, Constrained_polyline_id>;
+  using Subconstraints_to_conform = std::stack<Stack_info, std::vector<Stack_info>>;
+  Subconstraints_to_conform subconstraints_to_conform;
+
+  std::vector<CGAL::unordered_flat_set<Vertex_handle>> all_finite_edges;
+  bool update_all_finite_edges_ = false;
+
+  void update_all_finite_edges() {
+    if(!update_all_finite_edges_) {
+      update_all_finite_edges_ = true;
+      if(use_finite_edges_map()) {
+        all_finite_edges.clear();
+        all_finite_edges.resize(tr().number_of_vertices()+1);
+        for(auto e: tr().all_edges()) {
+          new_edge(e);
+        }
+      }
+    }
+  }
+
+  void new_vertex(Vertex_handle v) {
+    if(use_finite_edges_map() && v->time_stamp() >= all_finite_edges.size()) {
+      all_finite_edges.emplace_back();
+      CGAL_assertion(v->time_stamp() == all_finite_edges.size() - 1);
+    }
+  }
+
+  enum class Debug_flags {
+    Steiner_points = 0,
+    conforming,
+    input_faces,
+    missing_region,
+    regions,
+    copy_triangulation_into_hole,
+    validity,
+    use_older_cavity_algorithm,
+    debug_finite_edges_map,
+    use_finite_edges_map,
+    nb_of_flags
+  };
+  std::bitset<static_cast<int>(Debug_flags::nb_of_flags)> debug_flags{};
+  bool is_Delaunay = true;
+};
+
+} // end CGAL
+
+#endif // not DOXYGEN_RUNNING
+
+#
+
+#endif // CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H

Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3_fwd.h

@@ -0,0 +1,24 @@
+// Copyright (c) 2025  GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Laurent Rineau
+
+#ifndef CGAL_CONFORMING_CONSTRAINED_DELAUNAY_TRIANGULATION_3_FWD_H
+#define CGAL_CONFORMING_CONSTRAINED_DELAUNAY_TRIANGULATION_3_FWD_H
+
+#include <CGAL/license/Constrained_triangulation_3.h>
+
+namespace CGAL {
+
+template <typename Traits, typename Tr>
+class Conforming_constrained_Delaunay_triangulation_3;
+
+} // namespace CGAL
+
+#endif // CGAL_CONFORMING_CONSTRAINED_DELAUNAY_TRIANGULATION_3_FWD_H

Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_cell_base_3.h

@@ -0,0 +1,109 @@
+// Copyright (c) 2019-2024  GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Laurent Rineau
+
+#ifndef CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_CELL_BASE_3_H
+#define CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_CELL_BASE_3_H
+
+#include <CGAL/license/Constrained_triangulation_3.h>
+
+#include <CGAL/Triangulation_simplex_base_with_time_stamp.h>
+#include <CGAL/Conforming_constrained_Delaunay_triangulation_cell_data_3.h>
+#include <CGAL/Triangulation_cell_base_3.h>
+#include <CGAL/SMDS_3/io_signature.h>
+
+namespace CGAL {
+
+/**
+ * @ingroup PkgConstrainedTriangulation3Classes
+ * @brief Cell base class for the 3D conforming constrained Delaunay triangulation.
+ *
+ * This class is derived from its template parameter `CellBase` and provides additional functionality
+ * required by `Conforming_constrained_Delaunay_triangulation_3`.
+ *
+ * @tparam Traits The geometric traits class, which must be a model of `ConformingConstrainedDelaunayTriangulationTraits_3`.
+ *         It should be the same as the geometric traits class of the triangulation.
+ * @tparam CellBase The base class for the cell, which must be a model of `TriangulationCellBase_3`.
+ *
+ * @cgalModels{ConformingConstrainedDelaunayTriangulationCellBase_3}
+ *
+ * \sa `CGAL::Conforming_constrained_Delaunay_triangulation_vertex_base_3`
+ */
+template <typename Traits, typename CellBase = Triangulation_cell_base_3<Traits> >
+class Conforming_constrained_Delaunay_triangulation_cell_base_3
+    : public Triangulation_simplex_base_with_time_stamp<CellBase>
+{
+  using Base = Triangulation_simplex_base_with_time_stamp<CellBase>;
+  Conforming_constrained_Delaunay_triangulation_cell_data_3 ccdt_3_data_;
+
+public:
+  // To get correct cell type in TDS
+  template < class TDS3 >
+  struct Rebind_TDS {
+    typedef typename CellBase::template Rebind_TDS<TDS3>::Other Cb3;
+    typedef Conforming_constrained_Delaunay_triangulation_cell_base_3 <Traits, Cb3> Other;
+  };
+
+  // Constructors inherited from the base class
+  using Base::Base;
+
+  Conforming_constrained_Delaunay_triangulation_cell_data_3& ccdt_3_data() {
+    return ccdt_3_data_;
+  }
+
+  const Conforming_constrained_Delaunay_triangulation_cell_data_3& ccdt_3_data() const {
+    return ccdt_3_data_;
+  }
+
+  static std::string io_signature() {
+    static_assert(
+        std::is_same_v<
+            decltype(std::declval<Conforming_constrained_Delaunay_triangulation_cell_data_3>().face_constraint_index(0)), int>);
+
+    return Get_io_signature<Base>()() + "+(" + Get_io_signature<int>()() + ")[4]";
+  }
+
+  friend std::ostream&
+  operator<<(std::ostream& os,
+             const Conforming_constrained_Delaunay_triangulation_cell_base_3& c)
+  {
+    os << static_cast<const Base&>(c);
+    for( unsigned li = 0; li < 4; ++li ) {
+      if(IO::is_ascii(os)) {
+        os << " " << c.ccdt_3_data().face_constraint_index(li);
+      } else {
+        CGAL::write(os, c.ccdt_3_data().face_constraint_index(li));
+      }
+    }
+    return os;
+  }
+  friend std::istream&
+  operator>>(std::istream& is,
+             Conforming_constrained_Delaunay_triangulation_cell_base_3& c)
+  {
+    is >> static_cast<Base&>(c);
+    if(!is) return is;
+    for( int li = 0; li < 4; ++li ) {
+      int i;
+      if(IO::is_ascii(is)) {
+        is >> i;
+      } else {
+        CGAL::read(is, i);
+      }
+      if(!is) return is;
+      c.face_id[li] = i;
+    }
+    return is;
+  }
+};
+
+} // namespace CGAL
+
+#endif // CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_CELL_BASE_3_H

Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_cell_data_3.h

@@ -0,0 +1,86 @@
+// Copyright (c) 2019-2024  GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Laurent Rineau
+
+#ifndef CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_CELL_DATA_3_H
+#define CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_CELL_DATA_3_H
+
+#include <CGAL/license/Constrained_triangulation_3.h>
+
+#include <CGAL/Constrained_triangulation_3/internal/config.h>
+
+#include <array>
+#include <bitset>
+
+namespace CGAL {
+enum class CDT_3_cell_marker {
+  CLEAR = 0,
+  IN_REGION = 1,
+  VISITED = 1,
+  ON_REGION_BOUNDARY = 2,
+  nb_of_markers
+};
+
+/*!
+ * @ingroup PkgConstrainedTriangulation3Classes
+ * @brief Internal per-cell data for \cgal 3D conforming constrained Delaunay triangulations
+ *
+ * This class is an internal detail of the implementation of \cgal 3D conforming constrained Delaunay triangulations.
+ *
+ * Any model of the `ConformingConstrainedDelaunayTriangulationCellBase_3` concept must include one object of this type
+ * as a non-static data member.
+ */
+class Conforming_constrained_Delaunay_triangulation_cell_data_3 {
+  /// @cond SKIP_IN_MANUAL
+  template <typename Tr> friend class Conforming_constrained_Delaunay_triangulation_3_impl;
+  /// @endcond
+
+  std::array<CDT_3_signed_index, 4> face_id = { -1, -1, -1, -1 };
+  std::array<void*, 4> facet_2d = {nullptr, nullptr, nullptr, nullptr};
+  std::bitset<static_cast<unsigned>(CDT_3_cell_marker::nb_of_markers)> markers;
+
+  bool is_marked() const { return markers.any(); }
+  bool is_marked(CDT_3_cell_marker m) const { return markers.test(static_cast<unsigned>(m)); }
+  void set_mark(CDT_3_cell_marker m) { markers.set(static_cast<unsigned>(m)); }
+  void clear_mark(CDT_3_cell_marker m) { markers.reset(static_cast<unsigned>(m)); }
+  void clear_marks() { markers.reset(); }
+
+  template <typename Facet_handle>
+  void set_facet_constraint(int i, CDT_3_signed_index face_id,
+                            Facet_handle facet_2d)
+  {
+    this->face_id[unsigned(i)] = face_id;
+    this->facet_2d[unsigned(i)] = static_cast<void*>(facet_2d == Facet_handle{} ?  nullptr : std::addressof(*facet_2d));
+  }
+
+  template <typename CDT_2>
+  auto face_2 (const CDT_2& cdt, int i) const {
+    using Face = typename CDT_2::Face;
+    auto ptr = static_cast<Face*>(facet_2d[unsigned(i)]);
+    return cdt.tds().faces().iterator_to(*ptr);
+  }
+public:
+  /// @{
+  // @cond SKIP_IN_MANUAL
+  bool is_facet_constrained(int i) const { return face_id[unsigned(i)] >= 0; }
+
+  CDT_3_signed_index face_constraint_index(int i) const {
+    return face_id[unsigned(i)];
+  }
+
+  void set_face_constraint_index(int i, CDT_3_signed_index index) {
+    face_id[unsigned(i)] = index;
+  }
+  /// @endcond
+};
+
+} // namespace CGAL
+
+#endif // CGAL_CONSTRAINED_DELAUNAY_TRIANGULATION_CELL_DATA_3_H