Merge branch 'master' into Aos_2-fixes-efif

AABB_tree/include/CGAL/AABB_face_graph_triangle_primitive.h

@@ -21,7 +21,6 @@
 #include <CGAL/AABB_primitive.h>
 #include <CGAL/boost/graph/property_maps.h>
 #include <CGAL/Default.h>
-#include <boost/mpl/if.hpp>
 
 namespace CGAL {
 
@@ -57,9 +56,9 @@ template < class FaceGraph,
            class CacheDatum=Tag_false >
 class AABB_face_graph_triangle_primitive
 #ifndef DOXYGEN_RUNNING
-  : public AABB_primitive<typename boost::mpl::if_<OneFaceGraphPerTree,
-                                                   typename boost::graph_traits<FaceGraph>::face_descriptor,
-                                                   std::pair<typename boost::graph_traits<FaceGraph>::face_descriptor, const FaceGraph*> >::type,
+  : public AABB_primitive<std::conditional_t<OneFaceGraphPerTree::value,
+                                             typename boost::graph_traits<FaceGraph>::face_descriptor,
+                                             std::pair<typename boost::graph_traits<FaceGraph>::face_descriptor, const FaceGraph*> >,
                         Triangle_from_face_descriptor_map<
                           FaceGraph,
                           typename Default::Get<VertexPointPMap,
@@ -76,7 +75,7 @@ class AABB_face_graph_triangle_primitive
 {
   typedef typename Default::Get<VertexPointPMap, typename boost::property_map< FaceGraph, vertex_point_t>::const_type >::type VertexPointPMap_;
   typedef typename boost::graph_traits<FaceGraph>::face_descriptor FD;
-  typedef typename boost::mpl::if_<OneFaceGraphPerTree, FD, std::pair<FD, const FaceGraph*> >::type Id_;
+  typedef std::conditional_t<OneFaceGraphPerTree::value, FD, std::pair<FD, const FaceGraph*> > Id_;
 
   typedef Triangle_from_face_descriptor_map<FaceGraph,VertexPointPMap_>  Triangle_property_map;
   typedef One_point_from_face_descriptor_map<FaceGraph,VertexPointPMap_> Point_property_map;

AABB_tree/include/CGAL/AABB_halfedge_graph_segment_primitive.h

@@ -24,7 +24,6 @@
 #include <iterator>
 #include <boost/mpl/and.hpp>
 #include <CGAL/type_traits/is_iterator.h>
-#include <boost/mpl/if.hpp>
 
 #include <CGAL/Default.h>
 
@@ -70,9 +69,9 @@ template < class HalfedgeGraph,
            class CacheDatum = Tag_false >
 class AABB_halfedge_graph_segment_primitive
 #ifndef DOXYGEN_RUNNING
-  : public AABB_primitive<  typename boost::mpl::if_<OneHalfedgeGraphPerTree,
-                                                     typename boost::graph_traits<HalfedgeGraph>::edge_descriptor,
-                                                     std::pair<typename boost::graph_traits<HalfedgeGraph>::edge_descriptor, const HalfedgeGraph*> >::type,
+  : public AABB_primitive< std::conditional_t<OneHalfedgeGraphPerTree::value,
+                                              typename boost::graph_traits<HalfedgeGraph>::edge_descriptor,
+                                              std::pair<typename boost::graph_traits<HalfedgeGraph>::edge_descriptor, const HalfedgeGraph*> >,
                             Segment_from_edge_descriptor_map<
                               HalfedgeGraph,
                               typename Default::Get<VertexPointPMap,
@@ -89,7 +88,7 @@ class AABB_halfedge_graph_segment_primitive
 {
   typedef typename Default::Get<VertexPointPMap,typename boost::property_map< HalfedgeGraph,vertex_point_t>::const_type >::type  VertexPointPMap_;
   typedef typename boost::graph_traits<HalfedgeGraph>::edge_descriptor ED;
-  typedef typename boost::mpl::if_<OneHalfedgeGraphPerTree, ED, std::pair<ED, const HalfedgeGraph*> >::type Id_;
+  typedef std::conditional_t<OneHalfedgeGraphPerTree::value, ED, std::pair<ED, const HalfedgeGraph*> > Id_;
 
   typedef Segment_from_edge_descriptor_map<HalfedgeGraph,VertexPointPMap_>  Segment_property_map;
   typedef Source_point_from_edge_descriptor_map<HalfedgeGraph,VertexPointPMap_> Point_property_map;

AABB_tree/include/CGAL/AABB_tree/internal/Is_ray_intersection_geomtraits.h

@@ -32,11 +32,11 @@ BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(Has_cartesian_const_iterator_3,Cartesian_const
 
 template<typename GeomTraits>
 struct Is_ray_intersection_geomtraits
-: boost::mpl::and_< Has_ray_3<GeomTraits>,
-                    Has_construct_source_3<GeomTraits>,
-                    Has_vector_3<GeomTraits>,
-                    Has_construct_cartesian_const_iterator_3<GeomTraits>,
-                    Has_cartesian_const_iterator_3<GeomTraits> >::type
+: std::bool_constant< Has_ray_3<GeomTraits>::value &&
+                      Has_construct_source_3<GeomTraits>::value &&
+                      Has_vector_3<GeomTraits>::value &&
+                      Has_construct_cartesian_const_iterator_3<GeomTraits>::value &&
+                      Has_cartesian_const_iterator_3<GeomTraits>::value >
 {};
 
 

Algebraic_foundations/include/CGAL/number_utils.h

@@ -194,21 +194,21 @@ root_of( int k, Input_iterator begin, Input_iterator end ) {
 template< class Number_type >
 inline
 // select a Is_zero functor
-typename boost::mpl::if_c<
- ::std::is_same< typename Algebraic_structure_traits< Number_type >::Is_zero,
- Null_functor  >::value ,
+typename std::conditional_t<
+ std::is_same_v< typename Algebraic_structure_traits< Number_type >::Is_zero,
+ Null_functor  >,
   typename Real_embeddable_traits< Number_type >::Is_zero,
   typename Algebraic_structure_traits< Number_type >::Is_zero
->::type::result_type
+>::result_type
 is_zero( const Number_type& x ) {
     // We take the Algebraic_structure_traits<>::Is_zero functor by default. If it
     //  is not available, we take the Real_embeddable_traits functor
-    typename ::boost::mpl::if_c<
-        ::std::is_same<
+    std::conditional_t<
+        std::is_same_v<
              typename Algebraic_structure_traits< Number_type >::Is_zero,
-             Null_functor >::value ,
+             Null_functor > ,
        typename Real_embeddable_traits< Number_type >::Is_zero,
-       typename Algebraic_structure_traits< Number_type >::Is_zero >::type
+       typename Algebraic_structure_traits< Number_type >::Is_zero >
        is_zero;
 return is_zero( x );
 }

Alpha_wrap_3/benchmark/Alpha_wrap_3/CMakeLists.txt

@@ -0,0 +1,15 @@
+# Created by the script cgal_create_cmake_script
+# This is the CMake script for compiling a CGAL application.
+
+cmake_minimum_required(VERSION 3.1...3.20)
+project(Alpha_wrap_3_Benchmark)
+
+find_package(CGAL REQUIRED)
+
+include_directories (BEFORE ../../include ./Quality ./Robustness) # AW3 includes
+include_directories (BEFORE ../../../CGAL-Patches/include)
+
+# create a target per cppfile
+create_single_source_cgal_program("Performance/performance_benchmark.cpp")
+create_single_source_cgal_program("Quality/quality_benchmark.cpp")
+create_single_source_cgal_program("Robustness/robustness_benchmark.cpp")

Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/compute_performance_benchmark_data.py

@@ -0,0 +1,61 @@
+# Copyright (c) 2019-2023 Google LLC (USA).
+# All rights reserved.
+#
+# This file is part of CGAL (www.cgal.org).
+#
+# $URL$
+# $Id$
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+#
+# Author(s)     : Pierre Alliez
+#                 Michael Hemmer
+#                 Cedric Portaneri
+#
+#!/usr/bin/python
+
+import os, sys, subprocess, datetime, time, getopt
+
+def compute_performance_benchmark_data(execname, filename, alpha):
+
+ output = ""
+ cmd = ("/usr/bin/time", "-v",
+        execname, "-i",
+        filename, "-a", alpha)
+ proc = subprocess.Popen(
+     cmd,
+     stdout=subprocess.PIPE,
+     stderr=subprocess.PIPE,
+     start_new_session=True)
+
+ outs, errs = proc.communicate()
+ output = outs.decode("utf-8") + errs.decode("utf-8")
+
+ for output_line in output.split("\n"):
+   if "User time (seconds): " in output_line:
+     print(output_line[len("User time (seconds):  "):])
+     continue
+   if "Maximum resident set size (kbytes): " in output_line:
+     print(output_line[len("Maximum resident set size (kbytes):  "):])
+     continue
+
+def main(argv):
+  execname=""
+  filename=""
+  alpha=""
+  try:
+    opts, args = getopt.getopt(sys.argv[1:], 'e:i:a:')
+  except getopt.GetoptError:
+    sys.exit(2)
+  for opt, arg in opts:
+    if opt == "-e":
+      execname = arg
+    elif opt == "-i":
+      filename = arg
+    elif opt == "-a":
+      alpha = arg
+
+  compute_performance_benchmark_data(execname, filename, alpha)
+
+if __name__ == "__main__":
+  main(sys.argv[1:])

Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/generate_performance_benchmark_charts.py

@@ -0,0 +1,156 @@
+# Copyright (c) 2019-2023 Google LLC (USA).
+# All rights reserved.
+#
+# This file is part of CGAL (www.cgal.org).
+#
+# $URL$
+# $Id$
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+#
+# Author(s)     : Pierre Alliez
+#                 Michael Hemmer
+#                 Cedric Portaneri
+#
+#!/usr/bin/python
+
+import os, sys, subprocess, datetime, time, signal, getopt
+import numpy as np
+import matplotlib.pyplot as plt
+
+def main(argv):
+
+  inputdir=""
+  outputdir=""
+  commit_hash=""
+  alpha=""
+  do_diff=False
+  diffdir=""
+  diff_hash=""
+  try:
+    opts, args = getopt.getopt(sys.argv[1:], 'i:a:o:c:d:p:')
+  except getopt.GetoptError:
+    sys.exit(2)
+  for opt, arg in opts:
+    if opt == "-i":
+      inputdir = arg
+    elif opt == "-a":
+      alpha = arg
+    elif opt == "-o":
+      outputdir = arg
+    elif opt == "-c":
+      commit_hash = arg
+    elif opt == "-d":
+      diff_hash = arg
+      do_diff = True
+    elif opt == "-p":
+      diffdir = arg
+
+  all_metric = {
+      "Time_(second)" : {},
+      "Memory_Peak_(kbytes)" : {}}
+  num_input = 0
+  for filename in os.listdir(inputdir) :
+    new_path = os.path.join(inputdir,filename)
+    new_file = open(new_path)
+    is_empty_new = os.path.getsize(new_path) <= 1
+    if do_diff :
+      old_path = os.path.join(diffdir,filename)
+      old_file = open(old_path)
+      is_empty_old = os.path.getsize(old_path) <= 1
+      for key in all_metric:
+        if is_empty_new or is_empty_old :
+          new_val = 0.
+          old_val = 0.
+        else :
+          new_val = float(new_file.readline().rstrip('\n'))
+          old_val = float(old_file.readline().rstrip('\n'))
+        mesh_id = str(filename.split('.')[0])
+        all_metric[key][mesh_id] = [new_val, old_val]
+    else :
+      for key in all_metric:
+         if is_empty_new :
+           new_val = 0.
+         else :
+           new_val = float(new_file.readline().rstrip('\n'))
+         mesh_id = str(filename.split('.')[0])
+         all_metric[key][mesh_id] = [new_val, new_val]
+    num_input = num_input+1
+
+  # update .pdf chart
+  date_now = datetime.datetime.now()
+  date_for_filename = str(date_now.year) +"_"+ str(date_now.month) +"_"+ str(date_now.day) +"_"+ str(date_now.hour) +"h"+ str(date_now.minute) +"mn"
+  for key in all_metric:
+    goal = 0
+    num_el = range(len(all_metric[key]))
+    avg_diff_to_goal = 0.
+    avg = 0.
+    x1 = []
+    x2 = []
+    for value in all_metric[key].values() :
+      avg += value[0]
+      diff_to_goal = abs(value[1]-goal) - abs(value[0]-goal)
+      avg_diff_to_goal += diff_to_goal
+      x1.append(value[0])
+      x2.append(value[1])
+    avg_diff_to_goal /= float(len(all_metric[key]))
+    avg /= float(len(all_metric[key]))
+
+    plt.figure(figsize=(8,8))
+    if do_diff :
+      plt.hist(x2, bins=100, color='tab:green', alpha=0.5)
+    plt.hist(x1, bins=100, color='tab:blue', alpha=0.5)
+    plt.vlines(x = goal, ymin=plt.ylim()[0], ymax=plt.ylim()[1], linestyles='dashed')
+
+    title = ""
+    if do_diff :
+      title += "Diff between " + commit_hash + " and " + diff_hash + " on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha
+    else :
+      title += "Benchmarking on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha
+
+    avg_str = str(format(abs(avg), '.2f'))
+    if key == "Time_(second)" :
+      title += "\nIn average we spend " + avg_str + " seconds"
+    else :
+      title += "\nIn average we use up to " +  avg_str + " kbytes"
+
+    if do_diff and avg_diff_to_goal == 0. :
+      title += "\nNo change between the two commits"
+    elif do_diff :
+      avg_diff_str = str(format(abs(avg_diff_to_goal), '.2f'))
+      if key == "Time_(second)" :
+        if avg_diff_to_goal < 0 :
+          title += "\nIn average we get slower by "
+        else :
+          title += "\nIn average we get faster "
+        title += avg_diff_str + " seconds"
+      else :
+        if avg_diff_to_goal < 0 :
+          title += "\nIn average we use " +  avg_diff_str + " more"
+        else :
+          title += "\nIn average we use " +  avg_diff_str + " less"
+        title += " kbytes"
+
+    plt.title(title, fontsize=15)
+    plt.xlabel(key.replace("_"," "), fontsize=14)
+    plt.ylabel("# of meshes", fontsize=14)
+    plt.tick_params(axis="x", labelsize=9)
+    plt.tick_params(axis="y", labelsize=9)
+
+    chart_filename = ""
+    if do_diff :
+      chart_filename += "diff_"+commit_hash+"_"+diff_hash+"_"+key+"_"+date_for_filename+".pdf"
+    else :
+      chart_filename += "results_"+commit_hash+"_"+key+"_"+date_for_filename+".pdf"
+    chart_path = os.path.join(outputdir+"/charts",chart_filename)
+    if os.path.isfile(chart_path) :
+      os.remove(chart_path)
+    plt.savefig(chart_path, bbox_inches="tight")
+    plt.close()
+
+  print("pdf updated")
+
+  sys.exit()
+
+if __name__ == "__main__":
+  main(sys.argv[1:])

Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/performance_benchmark.cpp

@@ -0,0 +1,65 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/alpha_wrap_3.h>
+
+#include <CGAL/IO/polygon_soup_io.h>
+
+#include <array>
+#include <iostream>
+#include <string>
+#include <vector>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point_3 = K::Point_3;
+using Vector_3 = K::Vector_3;
+
+using Mesh = CGAL::Surface_mesh<Point_3>;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int argc, char** argv)
+{
+  const int argc_check = argc - 1;
+  const char* entry_name_ptr = nullptr;
+  double relative_alpha_ratio = 20., relative_offset_ratio = 600.;
+
+  for(int i=1; i<argc; ++i)
+  {
+    if(!strcmp("-i", argv[i]) && i < argc_check)
+      entry_name_ptr = argv[++i];
+    else if(!strcmp("-a", argv[i]) && i < argc_check)
+      relative_alpha_ratio = std::stod(argv[++i]);
+    else if(!strcmp("-d", argv[i]) && i < argc_check)
+      relative_offset_ratio = std::stod(argv[++i]);
+  }
+
+  if(argc < 3 || relative_alpha_ratio <= 0.)
+  {
+    std::cerr << "Error: bad input parameters." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::vector<Point_3> points;
+  std::vector<std::array<std::size_t, 3> > faces;
+  if(!CGAL::IO::read_polygon_soup(entry_name_ptr, points, faces) || faces.empty())
+  {
+    std::cerr << "Error: Invalid input data." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  CGAL::Bbox_3 bbox;
+  for(const Point_3& p : points)
+    bbox += p.bbox();
+
+  const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                                       CGAL::square(bbox.ymax() - bbox.ymin()) +
+                                       CGAL::square(bbox.zmax() - bbox.zmin()));
+  const double alpha = diag_length / relative_alpha_ratio;
+  const double offset = diag_length / relative_offset_ratio;
+
+  Mesh wrap;
+  CGAL::alpha_wrap_3(points, faces, alpha, offset, wrap);
+
+  return EXIT_SUCCESS;
+}

Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/compute_quality_benchmark_data.py

@@ -0,0 +1,54 @@
+# Copyright (c) 2019-2023 Google LLC (USA).
+# All rights reserved.
+#
+# This file is part of CGAL (www.cgal.org).
+#
+# $URL$
+# $Id$
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+#
+# Author(s)     : Pierre Alliez
+#                 Michael Hemmer
+#                 Cedric Portaneri
+#
+#!/usr/bin/python
+
+import os, sys, subprocess, datetime, time, getopt
+
+def compute_quality_benchmark_data(execname, filename, alpha):
+
+ output = ""
+ cmd = (execname, "-i",
+        filename, "-a", alpha)
+ proc = subprocess.Popen(
+     cmd,
+     stdout=subprocess.PIPE,
+     stderr=subprocess.PIPE,
+     start_new_session=True)
+
+ outs, errs = proc.communicate()
+ output = outs.decode("utf-8") + errs.decode("utf-8")
+
+ print(output)
+
+def main(argv):
+  execname=""
+  filename=""
+  alpha=""
+  try:
+    opts, args = getopt.getopt(sys.argv[1:], 'e:i:a:')
+  except getopt.GetoptError:
+    sys.exit(2)
+  for opt, arg in opts:
+    if opt == "-e":
+      execname = arg
+    elif opt == "-i":
+      filename = arg
+    elif opt == "-a":
+      alpha = arg
+
+  compute_quality_benchmark_data(execname, filename, alpha)
+
+if __name__ == "__main__":
+  main(sys.argv[1:])

Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/distance_utils.h

@@ -0,0 +1,151 @@
+// Copyright (c) 2019-2022 Google LLC (USA).
+// 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)     : Pierre Alliez
+//                 Michael Hemmer
+//                 Cedric Portaneri
+
+#ifndef CGAL_ALPHA_WRAP_3_BENCHMARK_ALPHA_WRAP_3_QUALITY_DISTANCE_H_
+#define CGAL_ALPHA_WRAP_3_BENCHMARK_ALPHA_WRAP_3_QUALITY_DISTANCE_H_
+
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/bounding_box.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+#include <CGAL/point_generators_3.h>
+#include <CGAL/Search_traits_3.h>
+
+namespace Aw3i {
+
+enum Distance_metric { HAUSDORFF = 0, MEAN = 1, RMS = 2 };
+
+template <typename Point, typename AABBTree>
+inline double approximate_hausdorff_distance(const std::vector<Point>& sample_points,
+                                             const AABBTree& tree,
+                                             Point& hint)
+{
+  double hdist = 0;
+  for(const Point& pt : sample_points)
+  {
+    hint = tree.closest_point(pt, hint);
+    auto dist = CGAL::squared_distance(hint, pt);
+    double d = CGAL::to_double(CGAL::approximate_sqrt(dist));
+    if(d > hdist)
+      hdist = d;
+  }
+
+  return hdist;
+}
+
+template <typename Point, typename AABBTree>
+inline double approximate_mean_distance(const std::vector<Point>& sample_points,
+                                        const AABBTree& tree,
+                                        Point& hint)
+{
+  double mdist = 0;
+  for(const Point& pt : sample_points)
+  {
+    hint = tree.closest_point(pt, hint);
+    auto dist = CGAL::squared_distance(hint, pt);
+    double d = CGAL::to_double(CGAL::approximate_sqrt(dist));
+    mdist += d;
+  }
+
+  return mdist / sample_points.size();
+}
+
+template <typename Point, typename AABBTree>
+inline double approximate_rms_distance(const std::vector<Point>& sample_points,
+                                       const AABBTree& tree,
+                                       Point& hint)
+{
+  double rmsdist = 0;
+  for(const Point& pt : sample_points)
+  {
+    hint = tree.closest_point(pt, hint);
+    auto dist = CGAL::squared_distance(hint, pt);
+    rmsdist += CGAL::to_double(dist);
+  }
+
+  return CGAL::to_double(CGAL::approximate_sqrt(rmsdist / sample_points.size()));
+}
+
+template <typename TriangleMesh>
+inline double approximate_distance(const TriangleMesh& tm1,
+                                   const TriangleMesh& tm2,
+                                   const Distance_metric& metric)
+{
+  using GT = typename CGAL::GetGeomTraits<TriangleMesh>::type;
+  using Point_3 = typename GT::Point_3;
+
+  using Primitive = CGAL::AABB_face_graph_triangle_primitive<TriangleMesh>;
+  using AABB_traits = CGAL::AABB_traits<GT, Primitive>;
+  using AABB_tree = CGAL::AABB_tree<AABB_traits>;
+
+  using CGAL::parameters::choose_parameter;
+  using CGAL::parameters::get_parameter;
+
+  std::vector<Point_3> original_sample_points;
+  CGAL::Polygon_mesh_processing::sample_triangle_mesh(tm1, std::back_inserter(original_sample_points),
+                                                      CGAL::parameters::all_default());
+
+  std::vector<Point_3> sample_points(std::begin(original_sample_points),
+                                     std::end(original_sample_points));
+  CGAL::spatial_sort(sample_points.begin(), sample_points.end());
+
+  AABB_tree tree(faces(tm2).first, faces(tm2).second, tm2);
+  tree.build();
+
+  auto vpm_2 = get(CGAL::vertex_point, tm2);
+  Point_3 hint = get(vpm_2, *vertices(tm2).first);
+
+  if(metric == HAUSDORFF)
+    return approximate_hausdorff_distance(sample_points, tree, hint);
+  else if(metric == MEAN)
+    return approximate_mean_distance(sample_points, tree, hint);
+  else if(metric == RMS)
+    return approximate_rms_distance(sample_points, tree, hint);
+  else
+    std::cerr << "Metric unknown\n" << std::endl;
+
+  return -1.0;
+}
+
+template <typename TriangleMesh>
+double get_longest_diag_bbox(const TriangleMesh& tm)
+{
+  CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(tm);
+  return std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                   CGAL::square(bbox.ymax() - bbox.ymin()) +
+                   CGAL::square(bbox.zmax() - bbox.zmin()));
+}
+
+template <typename TriangleMesh>
+inline double approximate_distance_relative_to_bbox(const TriangleMesh& tm1,
+                                                    const TriangleMesh& tm2,
+                                                    const Distance_metric& metric)
+{
+  double longest_diag_length = get_longest_diag_bbox(tm1);
+  return approximate_distance(tm1, tm2, metric) / longest_diag_length;
+}
+
+template <typename TriangleMesh, typename FT>
+inline double approximate_distance_relative_to_bbox(const TriangleMesh& tm1,
+                                                    const TriangleMesh& tm2,
+                                                    const Distance_metric& metric,
+                                                    const FT& longest_diag_length)
+{
+  return approximate_distance(tm1, tm2, metric) / CGAL::to_double(longest_diag_length);
+}
+
+} // namespace Aw3i
+
+#endif // CGAL_CGAL_ALPHA_WRAP_3_BENCHMARK_ALPHA_WRAP_3_QUALITY_DISTANCE_H_

Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/generate_quality_benchmark_charts.py

@@ -0,0 +1,182 @@
+# Copyright (c) 2019-2023 Google LLC (USA).
+# All rights reserved.
+#
+# This file is part of CGAL (www.cgal.org).
+#
+# $URL$
+# $Id$
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+#
+# Author(s)     : Pierre Alliez
+#                 Michael Hemmer
+#                 Cedric Portaneri
+#
+#!/usr/bin/python
+
+import os, sys, subprocess, datetime, time, signal, getopt
+import numpy as np
+import matplotlib.pyplot as plt
+
+def main(argv):
+
+  inputdir=""
+  outputdir=""
+  commit_hash=""
+  alpha=""
+  do_diff=False
+  diffdir=""
+  diff_hash=""
+  try:
+    opts, args = getopt.getopt(sys.argv[1:], 'i:a:o:c:d:p:')
+  except getopt.GetoptError:
+    sys.exit(2)
+  for opt, arg in opts:
+    if opt == "-i":
+      inputdir = arg
+    elif opt == "-a":
+      alpha = arg
+    elif opt == "-o":
+      outputdir = arg
+    elif opt == "-c":
+      commit_hash = arg
+    elif opt == "-d":
+      diff_hash = arg
+      do_diff = True
+    elif opt == "-p":
+      diffdir = arg
+
+  all_metric = {
+      "Mean_Min_Angle_(degree)" : {},
+      "Mean_Max_Angle_(degree)" : {},
+      "Mean_Radius_Ratio" : {},
+      "Mean_Edge_Ratio" : {},
+      "Mean_Aspect_Ratio" : {},
+      "Complexity_(#_of_triangle)" : {},
+      "#_of_almost_degenerate_triangle" : {},
+      "Hausdorff_distance_output_to_input_(%_of_bbox_diag)" : {}}
+  num_input = 0
+  print("inputdir = ", inputdir)
+  for filename in os.listdir(inputdir) :
+    new_path = os.path.join(inputdir,filename)
+    new_file = open(new_path)
+    if do_diff :
+      old_path = os.path.join(diffdir,filename)
+      old_file = open(old_path)
+      is_empty_old = os.path.getsize(old_path) <= 1
+      for key in all_metric :
+        try :
+          new_val = float(new_file.readline().rstrip('\n'))
+          old_val = float(old_file.readline().rstrip('\n'))
+          mesh_id = str(filename.split('.')[0])
+          all_metric[key][mesh_id] = [new_val, old_val]
+        except ValueError:
+          pass
+    else :
+      for key in all_metric :
+        try :
+          new_val = float(new_file.readline().rstrip('\n'))
+          mesh_id = str(filename.split('.')[0])
+          all_metric[key][mesh_id] = [new_val, new_val]
+        except ValueError:
+          pass
+    num_input = num_input+1
+
+  # update .pdf chart
+  date_now = datetime.datetime.now()
+  date_for_filename = str(date_now.year) +"_"+ str(date_now.month) +"_"+ str(date_now.day) +"_"+ str(date_now.hour) +"h"+ str(date_now.minute) +"mn"
+  for key in all_metric:
+    goal = 0
+    if key == "Mean_Min_Angle_(degree)" or key == "Mean_Max_Angle_(degree)":
+      goal = 60
+    elif key == "Mean_Radius_Ratio" or key == "Mean_Edge_Ratio" or key == "Mean_Aspect_Ratio" :
+      goal = 1
+
+    num_el = range(len(all_metric[key]))
+    avg_diff_to_goal = 0.
+    avg = 0.
+    x1 = []
+    x2 = []
+    for value in all_metric[key].values() :
+      avg += value[0]
+      diff_to_goal = abs(value[1]-goal) - abs(value[0]-goal)
+      avg_diff_to_goal += diff_to_goal
+      x1.append(value[0])
+      x2.append(value[1])
+    avg_diff_to_goal /= float(len(all_metric[key]))
+    avg /= float(len(all_metric[key]))
+
+    plt.figure(figsize=(8,8))
+    if do_diff :
+      plt.hist(x2, bins=100, color='tab:green', alpha=0.5)
+    plt.hist(x1, bins=100, color='tab:blue', alpha=0.5)
+    plt.vlines(x = goal, ymin=plt.ylim()[0], ymax=plt.ylim()[1], linestyles='dashed')
+
+    title = ""
+    if do_diff :
+      title += "Diff between " + commit_hash + " and " + diff_hash + " on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha
+    else :
+      title += "Benchmarking on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha
+
+    avg_str = str(format(abs(avg), '.2f'))
+    if key == "Mean_Min_Angle_(degree)" or key == "Mean_Max_Angle_(degree)":
+      title += "\nIn average we have " + avg_str + "°"
+    elif key == "Mean_Radius_Ratio" or key == "Mean_Edge_Ratio" or key == "Mean_Aspect_Ratio" :
+      title += "\nIn average we have a ratio of " +  avg_str
+    elif key == "Hausdorff_distance_output_to_input_(%_of_bbox_diag)" :
+      title += "\nIn average we have a distance of " +  avg_str + "% of bbox diag"
+    elif key == "Complexity_(#_of_triangle)" or key == "#_of_almost_degenerate_triangle" :
+      title += "\nIn average we have " +  avg_str + " triangles"
+
+    if do_diff and avg_diff_to_goal == 0. :
+      title += "\nNo change between the two commits"
+    elif do_diff :
+      avg_diff_str = str(format(abs(avg_diff_to_goal), '.2f'))
+      if key == "Mean_Min_Angle_(degree)" or key == "Mean_Max_Angle_(degree)":
+        if avg_diff_to_goal < 0 :
+          title += "\nIn average we loose "
+        else :
+          title += "\nIn average we gain "
+        title += avg_diff_str + "° toward 60°"
+      elif key == "Mean_Radius_Ratio" or key == "Mean_Edge_Ratio" or key == "Mean_Aspect_Ratio" :
+        if avg_diff_to_goal < 0 :
+          title += "\nIn average we loose "
+        else :
+          title += "\nIn average we gain "
+        title += avg_diff_str + " of ratio toward 1"
+      elif key == "Hausdorff_distance_output_to_input_(%_of_bbox_diag)" :
+        if avg_diff_to_goal < 0 :
+          title += "\nIn average we increase by "
+        else :
+          title += "\nIn average we reduce by "
+        title += avg_diff_str + " the bbox ratio"
+      elif key == "Complexity_(#_of_triangle)" or key == "#_of_almost_degenerate_triangle" :
+        if avg_diff_to_goal < 0 :
+          title += "\nIn average we get " + avg_diff_str + " more"
+        else :
+          title += "\nIn average we get " +  avg_diff_str + " less"
+        title += " triangles"
+
+    plt.title(title, fontsize=15)
+    plt.xlabel(key.replace("_"," "), fontsize=14)
+    plt.ylabel("# of meshes", fontsize=14)
+    plt.tick_params(axis="x", labelsize=9)
+    plt.tick_params(axis="y", labelsize=9)
+
+    chart_filename = ""
+    if do_diff :
+      chart_filename += "diff_"+commit_hash+"_"+diff_hash+"_"+key+"_"+date_for_filename+".pdf"
+    else :
+      chart_filename += "results_"+commit_hash+"_"+key+"_"+date_for_filename+".pdf"
+    chart_path = os.path.join(outputdir+"/charts",chart_filename)
+    if os.path.isfile(chart_path) :
+      os.remove(chart_path)
+    plt.savefig(chart_path, bbox_inches="tight")
+    plt.close()
+
+  print("pdf updated")
+
+  sys.exit()
+
+if __name__ == "__main__":
+  main(sys.argv[1:])

Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/quality_benchmark.cpp

@@ -0,0 +1,271 @@
+#include <distance_utils.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/alpha_wrap_3.h>
+
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+
+#include <array>
+#include <cmath>
+#include <iostream>
+
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point_3 = Kernel::Point_3;
+using Vector_3 = Kernel::Vector_3;
+using Triangle_3 = Kernel::Triangle_3;
+using FT = Kernel::FT;
+
+using Mesh = CGAL::Surface_mesh<Point_3>;
+using face_descriptor = boost::graph_traits<Mesh>::face_descriptor;
+
+using Oracle = CGAL::Alpha_wraps_3::internal::Triangle_mesh_oracle<Kernel>;
+using Dt = CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle>::Triangulation;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+std::array<FT, 3> triangle_angles(const Triangle_3& tr)
+{
+  FT sq_a = CGAL::squared_distance(tr[0], tr[1]);
+  FT sq_b = CGAL::squared_distance(tr[1], tr[2]);
+  FT sq_c = CGAL::squared_distance(tr[2], tr[0]);
+
+  FT two_ab = 2. * CGAL::sqrt(sq_a) * CGAL::sqrt(sq_b);
+  FT two_bc = 2. * CGAL::sqrt(sq_b) * CGAL::sqrt(sq_c);
+  FT two_ca = 2. * CGAL::sqrt(sq_c) * CGAL::sqrt(sq_a);
+
+  FT angle_a = (sq_b + sq_c - sq_a) / two_bc;
+  FT angle_b = (sq_c + sq_a - sq_b) / two_ca;
+  FT angle_c = (sq_a + sq_b - sq_c) / two_ab;
+  if(angle_a < -1.) angle_a = -1.;
+  if(angle_b < -1.) angle_b = -1.;
+  if(angle_c < -1.) angle_c = -1.;
+  if(angle_a > 1.) angle_a = 1.;
+  if(angle_b > 1.) angle_b = 1.;
+  if(angle_c > 1.) angle_c = 1.;
+  angle_a = std::acos(angle_a);
+  angle_b = std::acos(angle_b);
+  angle_c = std::acos(angle_c);
+
+  return {angle_a, angle_b, angle_c};
+}
+
+bool is_almost_degenerate(const Triangle_3& tr,
+                          double threshold)
+{
+  FT sq_area = tr.squared_area();
+  return (CGAL::sqrt(CGAL::to_double(sq_area)) < threshold);
+}
+
+auto surface_mesh_face_to_triangle(const face_descriptor fd,
+                                   const Mesh& sm)
+{
+  typename boost::graph_traits<Mesh>::halfedge_descriptor hd = halfedge(fd,sm);
+   return Triangle_3(sm.point(target(hd,sm)),
+                     sm.point(target(next(hd,sm),sm)),
+                     sm.point(target(next(next(hd,sm),sm),sm)));
+}
+
+double mean_min_angle(const Mesh& mesh)
+{
+  double mean_min_angle = 0.;
+  for(const face_descriptor f : faces(mesh))
+  {
+    const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh);
+    std::array<FT, 3> angles = triangle_angles(tr);
+
+    FT min_angle = std::min({angles[0], angles[1], angles[2]});
+
+    min_angle = min_angle * (180.0 / CGAL_PI);
+    mean_min_angle += min_angle;
+  }
+
+  mean_min_angle /= static_cast<double>(mesh.number_of_faces());
+  return mean_min_angle;
+}
+
+double mean_max_angle(const Mesh& mesh)
+{
+  double mean_max_angle = 0.;
+  for(const face_descriptor f : faces(mesh))
+  {
+    const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh);
+    std::array<FT, 3> angles = triangle_angles(tr);
+
+    FT max_angle = std::max({angles[0], angles[1], angles[2]});
+
+    max_angle = max_angle * (180.0 / CGAL_PI);
+    mean_max_angle += max_angle;
+  }
+
+  mean_max_angle /= static_cast<double>(mesh.number_of_faces());
+  return mean_max_angle;
+}
+
+double mean_radius_ratio(const Mesh& mesh,
+                         double degenerate_threshold)
+{
+  double mean_radius_ratio = 0.;
+  size_t num_almost_degenerate_tri = 0;
+  for(const face_descriptor f : faces(mesh))
+  {
+    const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh);
+    if(is_almost_degenerate(tr, degenerate_threshold))
+    {
+      ++num_almost_degenerate_tri;
+      continue;
+    }
+
+    FT circumsphere_radius = std::sqrt(CGAL::squared_radius(tr[0], tr[1], tr[2]));
+
+    FT a = std::sqrt(CGAL::squared_distance(tr[0], tr[1]));
+    FT b = std::sqrt(CGAL::squared_distance(tr[1], tr[2]));
+    FT c = std::sqrt(CGAL::squared_distance(tr[2], tr[0]));
+    FT s = 0.5 * (a + b + c);
+    FT inscribed_radius = std::sqrt((s * (s - a) * (s - b) * (s - c)) / s);
+    FT radius_ratio = circumsphere_radius / inscribed_radius;
+    radius_ratio /= 2.;  // normalized
+    mean_radius_ratio += radius_ratio;
+  }
+
+  mean_radius_ratio /= static_cast<double>(mesh.number_of_faces() - num_almost_degenerate_tri);
+  return mean_radius_ratio;
+}
+
+double mean_edge_ratio(const Mesh& mesh,
+                       double degenerate_threshold)
+{
+  double mean_edge_ratio = 0.;
+  size_t num_almost_degenerate_tri = 0;
+
+  for(const face_descriptor f : faces(mesh))
+  {
+    const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh);
+    if(is_almost_degenerate(tr, degenerate_threshold))
+    {
+      ++num_almost_degenerate_tri;
+      continue;
+    }
+
+    FT a = std::sqrt(CGAL::squared_distance(tr[0], tr[1]));
+    FT b = std::sqrt(CGAL::squared_distance(tr[1], tr[2]));
+    FT c = std::sqrt(CGAL::squared_distance(tr[2], tr[0]));
+    FT min_edge = std::min({a, b, c});
+    FT max_edge = std::max({a, b, c});
+    FT edge_ratio = max_edge / min_edge;
+
+    mean_edge_ratio += edge_ratio;
+  }
+
+  mean_edge_ratio /= static_cast<double>(mesh.number_of_faces() - num_almost_degenerate_tri);
+  return mean_edge_ratio;
+}
+
+double mean_aspect_ratio(const Mesh& mesh,
+                         double degenerate_threshold)
+{
+  double mean_aspect_ratio = 0.;
+  size_t num_almost_degenerate_tri = 0;
+  for(const face_descriptor f : faces(mesh))
+  {
+    const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh);
+    if(is_almost_degenerate(tr, degenerate_threshold))
+    {
+      ++num_almost_degenerate_tri;
+      continue;
+    }
+
+    FT a = std::sqrt(CGAL::squared_distance(tr[0], tr[1]));
+    FT b = std::sqrt(CGAL::squared_distance(tr[1], tr[2]));
+    FT c = std::sqrt(CGAL::squared_distance(tr[2], tr[0]));
+    FT s = 0.5 * (a + b + c);
+    FT inscribed_radius = std::sqrt((s * (s - a) * (s - b) * (s - c)) / s);
+    FT max_edge = std::max({a, b, c});
+    FT aspect_ratio = max_edge / inscribed_radius;
+    aspect_ratio /= (2. * std::sqrt(3.));  // normalized
+    mean_aspect_ratio += aspect_ratio;
+  }
+
+  mean_aspect_ratio /= static_cast<double>(mesh.number_of_faces() - num_almost_degenerate_tri);
+  return mean_aspect_ratio;
+}
+
+size_t num_almost_degenerate_tri(const Mesh& mesh,
+                                 double degenerate_threshold)
+{
+  size_t num_almost_degenerate_tri = 0;
+  for(const face_descriptor f : faces(mesh))
+  {
+    const Triangle_3 tr = surface_mesh_face_to_triangle(f, mesh);
+    if(is_almost_degenerate(tr, degenerate_threshold))
+    {
+      ++num_almost_degenerate_tri;
+    }
+  }
+  return num_almost_degenerate_tri;
+}
+
+int main(int argc, char** argv)
+{
+  const int argc_check = argc - 1;
+  char *entry_name_ptr = nullptr;
+  double relative_alpha_ratio = 20.;
+  double relative_offset_ratio = 600.;
+
+  for(int i=1; i<argc; ++i)
+  {
+    if(!strcmp("-i", argv[i]) && i < argc_check) {
+      entry_name_ptr = argv[++i];
+    } else if(!strcmp("-a", argv[i]) && i < argc_check) {
+      relative_alpha_ratio = std::stod(argv[++i]);
+    } else if(!strcmp("-d", argv[i]) && i < argc_check) {
+      relative_offset_ratio = std::stod(argv[++i]);
+    }
+  }
+
+  if(argc < 3 || relative_alpha_ratio <= 0.)
+  {
+    std::cerr << "Error: bad input parameters." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  Mesh input_mesh;
+  if(!PMP::IO::read_polygon_mesh(entry_name_ptr, input_mesh) ||
+     is_empty(input_mesh) ||
+     !is_triangle_mesh(input_mesh))
+  {
+    return EXIT_FAILURE;
+  }
+
+  CGAL::Bbox_3 bbox = PMP::bbox(input_mesh);
+  const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                                       CGAL::square(bbox.ymax() - bbox.ymin()) +
+                                       CGAL::square(bbox.zmax() - bbox.zmin()));
+  const double alpha = diag_length / relative_alpha_ratio;
+  const double offset = diag_length / relative_offset_ratio;
+
+  Mesh wrap;
+  CGAL::alpha_wrap_3(input_mesh, alpha, offset, wrap);
+
+  double degenerate_threshold = 0.;
+  for(const face_descriptor f : faces(wrap))
+  {
+    const Triangle_3 tr = surface_mesh_face_to_triangle(f, wrap);
+    degenerate_threshold += CGAL::sqrt(CGAL::to_double(tr.squared_area()));
+  }
+
+  degenerate_threshold /= wrap.number_of_faces();
+  degenerate_threshold /= 1000;
+
+  std::cout << mean_min_angle(wrap) << "\n";
+  std::cout << mean_max_angle(wrap) << "\n";
+  std::cout << mean_radius_ratio(wrap, degenerate_threshold) << "\n";
+  std::cout << mean_edge_ratio(wrap, degenerate_threshold) << "\n";
+  std::cout << mean_aspect_ratio(wrap, degenerate_threshold) << "\n";
+  std::cout << wrap.number_of_faces() << "\n";
+  std::cout << num_almost_degenerate_tri(wrap, degenerate_threshold) << "\n";
+  std::cout << 100. * approximate_distance_relative_to_bbox(wrap, input_mesh, Aw3i::HAUSDORFF) << "\n";
+
+  return 0;
+}

Alpha_wrap_3/benchmark/Alpha_wrap_3/Robustness/compute_robustness_benchmark_data.py

@@ -0,0 +1,97 @@
+# Copyright (c) 2019-2023 Google LLC (USA).
+# All rights reserved.
+#
+# This file is part of CGAL (www.cgal.org).
+#
+# $URL$
+# $Id$
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+#
+# Author(s)     : Pierre Alliez
+#                 Michael Hemmer
+#                 Cedric Portaneri
+#
+#!/usr/bin/python
+
+import os, sys, subprocess, datetime, time, signal, getopt
+
+def signal_handler(signum, frame):
+  raise Exception("Timed out!")
+
+def compute_robustness_benchmark_data(execname, filename, alpha, max_time):
+
+  exit_codes = {
+   0 : "VALID_SOLID_OUTPUT",
+   1 : "INPUT_IS_INVALID",
+   2 : "OUTPUT_IS_NOT_TRIANGLE_MESH",
+   3 : "OUTPUT_IS_COMBINATORIAL_NON_MANIFOLD",
+   4 : "OUTPUT_HAS_BORDERS",
+   5 : "OUTPUT_HAS_DEGENERATED_FACES",
+   6 : "OUTPUT_HAS_GEOMETRIC_SELF_INTERSECTIONS",
+   7 : "OUTPUT_DOES_NOT_BOUND_VOLUME",
+   8 : "OUTPUT_DOES_NOT_CONTAIN_INPUT",
+   9 : "OUTPUT_DISTANCE_IS_TOO_LARGE",
+   10 : "SIGSEGV",
+   11 : "SIGABRT",
+   12 : "SIGFPE",
+   13 : "TIMEOUT"
+  }
+
+  exit_code = 0
+  output = ""
+  cmd = ("/usr/bin/time", "-v",
+         execname, "-i",
+         filename, "-a", alpha)
+  proc = subprocess.Popen(
+      cmd,
+      stdout=subprocess.PIPE,
+      stderr=subprocess.PIPE,
+      start_new_session=True)
+
+  try:
+    outs, errs = proc.communicate(timeout=int(max_time))
+    exit_code = proc.returncode
+    output = outs.decode("utf-8") + errs.decode("utf-8")
+
+    for output_line in output.split("\n"):
+      if output_line == "Command terminated by signal 11":
+        exit_code = 10
+        continue
+      elif output_line == "Command terminated by signal 6":
+        exit_code = 11
+        continue
+      elif output_line == "Command terminated by signal 8":
+        exit_code = 12
+        continue
+
+  except subprocess.TimeoutExpired:
+    os.killpg(os.getpgid(proc.pid), signal.SIGTERM)
+    exit_code = 13
+    output = "process ran too long"
+
+  print(exit_codes[exit_code])
+
+def main(argv):
+  execname=""
+  filename=""
+  alpha=""
+  max_time=""
+  try:
+    opts, args = getopt.getopt(sys.argv[1:], 'e:i:a:t:')
+  except getopt.GetoptError:
+    sys.exit(2)
+  for opt, arg in opts:
+    if opt == "-e":
+      execname = arg
+    elif opt == "-i":
+      filename = arg
+    elif opt == "-a":
+      alpha = arg
+    elif opt == "-t":
+      max_time = arg
+
+  compute_robustness_benchmark_data(execname, filename, alpha, max_time)
+
+if __name__ == "__main__":
+  main(sys.argv[1:])

Alpha_wrap_3/benchmark/Alpha_wrap_3/Robustness/generate_robustness_benchmark_charts.py

@@ -0,0 +1,156 @@
+# Copyright (c) 2019-2023 Google LLC (USA).
+# All rights reserved.
+#
+# This file is part of CGAL (www.cgal.org).
+#
+# $URL$
+# $Id$
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+#
+# Author(s)     : Pierre Alliez
+#                 Michael Hemmer
+#                 Cedric Portaneri
+#
+#!/usr/bin/python
+
+import os, sys, subprocess, datetime, time, signal, getopt
+import numpy as np
+import matplotlib.pyplot as plt
+
+def main(argv):
+
+  inputdir=""
+  outputdir=""
+  commit_hash=""
+  alpha=""
+  try:
+    opts, args = getopt.getopt(sys.argv[1:], 'i:a:o:c:')
+  except getopt.GetoptError:
+    sys.exit(2)
+  for opt, arg in opts:
+    if opt == "-i":
+      inputdir = arg
+    elif opt == "-a":
+      alpha = arg
+    elif opt == "-o":
+      outputdir = arg
+    elif opt == "-c":
+      commit_hash = arg
+
+  exit_codes = {
+   0 : "VALID_SOLID_OUTPUT",
+   1 : "INPUT_IS_INVALID",
+   2 : "OUTPUT_IS_NOT_TRIANGLE_MESH",
+   3 : "OUTPUT_IS_COMBINATORIAL_NON_MANIFOLD",
+   4 : "OUTPUT_HAS_BORDERS",
+   5 : "OUTPUT_HAS_DEGENERATED_FACES",
+   6 : "OUTPUT_HAS_GEOMETRIC_SELF_INTERSECTIONS",
+   7 : "OUTPUT_DOES_NOT_BOUND_VOLUME",
+   8 : "OUTPUT_DOES_NOT_CONTAIN_INPUT",
+   9 : "OUTPUT_DISTANCE_IS_TOO_LARGE",
+   10 : "SIGSEGV",
+   11 : "SIGABRT",
+   12 : "SIGFPE",
+   13 : "TIMEOUT"
+  }
+
+  current_run_data = {
+   "VALID_SOLID_OUTPUT" : 0,
+   "INPUT_IS_INVALID" : 0,
+   "OUTPUT_IS_NOT_TRIANGLE_MESH" : 0,
+   "OUTPUT_IS_COMBINATORIAL_NON_MANIFOLD" : 0,
+   "OUTPUT_HAS_BORDERS" : 0,
+   "OUTPUT_HAS_DEGENERATED_FACES" : 0,
+   "OUTPUT_HAS_GEOMETRIC_SELF_INTERSECTIONS" : 0,
+   "OUTPUT_DOES_NOT_BOUND_VOLUME" : 0,
+   "OUTPUT_DOES_NOT_CONTAIN_INPUT" : 0,
+   "OUTPUT_DISTANCE_IS_TOO_LARGE" : 0,
+   "SIGSEGV" : 0,
+   "SIGABRT" : 0,
+   "SIGFPE" : 0,
+   "TIMEOUT" : 0
+  }
+
+  filenames_per_codes = {}
+  for key in current_run_data :
+   filenames_per_codes[key] = []
+
+  print("inputdir = ", inputdir)
+
+  num_input = 0
+  for filename in os.listdir(inputdir) :
+    print("filename = ", filename)
+
+    f = open(os.path.join(inputdir,filename))
+    status = f.readline().rstrip('\n');
+    current_run_data[status] += 1
+    filenames_per_codes[status].append(filename.rstrip('.log'))
+    num_input = num_input+1
+
+  # sort current_run_data by value
+  current_run_data = {k: v for k, v in sorted(current_run_data.items(), key=lambda item: item[1], reverse=True)}
+
+  # update chart data files
+  date_now = datetime.datetime.now()
+  date = str(date_now.year) +"-"+ str(date_now.month) +"-"+ str(date_now.day) +" "+ str(date_now.hour) +"h"+ str(date_now.minute) +"mn"
+
+  for key_filename in current_run_data:
+    f = open(os.path.join(outputdir+"/charts_data", key_filename+".txt"), "a+")
+    f.write(str(current_run_data[key_filename]) + " " + commit_hash + " " + date + "\n")
+
+  print("chart data updated")
+
+  # update .pdf chart
+  chart = plt.figure(figsize=(10, 7))
+  colormap = ["tab:blue","tab:orange","tab:green","tab:red","tab:purple","tab:brown","tab:pink","tab:gray","tab:olive","tab:cyan","b","palegreen", "peachpuff"]
+  plt.gca().set_prop_cycle('color', colormap)
+  plt.style.use('tableau-colorblind10')
+  for key_filename in current_run_data:
+
+    f = open(os.path.join(outputdir+"/charts_data", key_filename+".txt"), "r")
+    lines = f.readlines()
+    x_number_values = []
+    y_number_values = []
+    i = 0
+    for line in lines :
+      if i < (len(lines) - 10) :
+        i=i+1
+        continue
+
+      i=i+1
+      words = line.strip().split()
+      x_number_values.append(words[1]+"\n"+words[2]+"\n"+words[3])
+      y_number_values.append(int(words[0]))
+    plt.plot(x_number_values, y_number_values, marker='o', label=key_filename+": "+str(current_run_data[key_filename]))
+
+  plt.xlabel("Version", fontsize=14)
+  plt.ylabel("# of mesh", fontsize=14)
+  plt.tick_params(axis="both", labelsize=9)
+  plt.title("Benchmarking on " + str(num_input) + " meshes from Thingi10K\nAlpha = Bbox diag length / " + alpha, fontsize=15)
+  plt.legend(loc='lower left', bbox_to_anchor= (1.01, 0.58), ncol=1,
+  borderaxespad=0, frameon=False)
+
+  date_for_filename = str(date_now.year) +"-"+ str(date_now.month) +"-"+ str(date_now.day) +"-"+ str(date_now.hour) +"h"+ str(date_now.minute) +"mn"
+  chart_filename = os.path.join(outputdir+"/charts","benchmarking_version_"+commit_hash+"-"+date_for_filename+".pdf")
+  if os.path.isfile(chart_filename) :
+     os.remove(chart_filename)
+  chart.savefig(chart_filename, bbox_inches="tight")
+  plt.close(chart)
+
+  print("pdf updated")
+
+  # dump filenames per codes
+  log_dirname = os.path.join(outputdir, "log/"+commit_hash+"-"+date_for_filename)
+  if not os.path.exists(log_dirname):
+    os.mkdir(log_dirname)
+  for key in filenames_per_codes :
+   file = open(os.path.join(log_dirname, key+".txt"), "w+")
+   for filename in filenames_per_codes[key] :
+     file.write(filename + "\n")
+   file.close()
+
+  sys.exit()
+
+if __name__ == "__main__":
+  main(sys.argv[1:])

Alpha_wrap_3/benchmark/Alpha_wrap_3/Robustness/robustness_benchmark.cpp

@@ -0,0 +1,111 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/alpha_wrap_3.h>
+#include <CGAL/Alpha_wrap_3/internal/validation.h>
+
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point_3 = Kernel::Point_3;
+using Mesh = CGAL::Surface_mesh<Point_3>;
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+namespace {
+
+enum Robustness_benchmark_exit_code
+{
+  // Success
+  VALID_SOLID_OUTPUT = 0,
+
+  // Failure
+  INPUT_IS_INVALID = 1,
+  OUTPUT_IS_NOT_TRIANGLE_MESH = 2,
+  OUTPUT_IS_COMBINATORIAL_NON_MANIFOLD = 3,
+  OUTPUT_HAS_BORDERS = 4,
+  OUTPUT_HAS_DEGENERATED_FACES = 5,
+  OUTPUT_HAS_GEOMETRIC_SELF_INTERSECTIONS = 6,
+  OUTPUT_DOES_NOT_BOUND_VOLUME = 7,
+  OUTPUT_DOES_NOT_CONTAIN_INPUT = 8,
+  OUTPUT_DISTANCE_IS_TOO_LARGE = 9,
+};
+
+} // namespace
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+namespace AW3i = CGAL::Alpha_wraps_3::internal;
+
+int main(int argc, char** argv)
+{
+  const int argc_check = argc - 1;
+  char* entry_name_ptr = nullptr;
+  double relative_alpha_ratio = 20.;
+  double relative_offset_ratio = 600.;
+
+  for(int i=1; i<argc; ++i)
+  {
+    if(!strcmp("-i", argv[i]) && i < argc_check) {
+      entry_name_ptr = argv[++i];
+    } else if(!strcmp("-a", argv[i]) && i < argc_check) {
+      relative_alpha_ratio = std::stod(argv[++i]);
+    } else if(!strcmp("-d", argv[i]) && i < argc_check) {
+      relative_offset_ratio = std::stod(argv[++i]);
+    }
+  }
+
+  if(argc < 3 || relative_alpha_ratio <= 0.)
+    return AW3i::INPUT_IS_INVALID;
+
+  Mesh input_mesh;
+  if(!PMP::IO::read_polygon_mesh(entry_name_ptr, input_mesh) ||
+     is_empty(input_mesh) ||
+     !is_triangle_mesh(input_mesh)
+#ifndef CGAL_ALPHA_WRAP_3_TOLERATE_DEGENERACIES
+     || AW3i::has_degenerated_faces(input_mesh)
+#endif
+     )
+  {
+    return AW3i::INPUT_IS_INVALID;
+  }
+
+  const CGAL::Bbox_3 bbox = PMP::bbox(input_mesh);
+  const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                                       CGAL::square(bbox.ymax() - bbox.ymin()) +
+                                       CGAL::square(bbox.zmax() - bbox.zmin()));
+  const double alpha = diag_length / relative_alpha_ratio;
+  const double offset = diag_length / relative_offset_ratio;
+
+  Mesh wrap;
+  alpha_wrap_3(input_mesh, alpha, offset, wrap);
+
+  if(!is_triangle_mesh(wrap))
+    return AW3i::OUTPUT_IS_NOT_TRIANGLE_MESH;
+
+  if(!is_closed(wrap))
+    return AW3i::OUTPUT_HAS_BORDERS;
+
+  if(AW3i::has_degenerated_faces(wrap))
+    return AW3i::OUTPUT_HAS_DEGENERATED_FACES;
+
+  if(AW3i::is_combinatorially_non_manifold(wrap))
+    return AW3i::OUTPUT_IS_COMBINATORIAL_NON_MANIFOLD;
+
+  if(PMP::does_self_intersect(wrap))
+    return AW3i::OUTPUT_HAS_GEOMETRIC_SELF_INTERSECTIONS;
+
+  if(!PMP::does_bound_a_volume(wrap))
+    return AW3i::OUTPUT_DOES_NOT_BOUND_VOLUME;
+
+  if(!AW3i::is_outer_wrap_of_triangle_mesh(wrap, input_mesh))
+    return AW3i::OUTPUT_DOES_NOT_CONTAIN_INPUT;
+
+  if(!AW3i::has_expected_Hausdorff_distance(wrap, input_mesh, alpha, offset))
+     return AW3i::OUTPUT_DISTANCE_IS_TOO_LARGE;
+
+  return AW3i::VALID_SOLID_OUTPUT;
+}

Alpha_wrap_3/benchmark/Alpha_wrap_3/benchmarking.sh

@@ -0,0 +1,86 @@
+# Copyright (c) 2019-2023 Google LLC (USA).
+# All rights reserved.
+#
+# This file is part of the 3D Alpha Wrapping package, which is being prepared for
+# submission to CGAL (www.cgal.org).
+#
+# $URL$
+# $Id$
+# SPDX-License-Identifier: GPL-3.0-or-later
+#
+#
+# Author(s)     : Pierre Alliez
+#                 Michael Hemmer
+#                 Cedric Portaneri
+#                 Mael Rouxel-Labbé
+#
+
+#!/bin/bash
+
+# $1: directory containing the alpha wrap project
+# $2: directory containing the output results
+# $3: alpha value
+# $4: timeout value in seconds
+# $5: hash of the latest commit
+# $6: the input file path
+function compute_benchmark_data() {
+  filename=$(basename -- "$6")
+  filename="${filename%.*}"
+
+  python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Robustness/compute_robustness_benchmark_data.py \
+   -e $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/build-release/robustness_benchmark -i $6 -a $3 -t $4 \
+   > $2/Robustness/results/$5/$filename.log
+
+  python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/compute_performance_benchmark_data.py \
+   -e $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/build-release/performance_benchmark -i $6 -a $3 \
+   > $2/Performance/results/$5/$filename.log
+
+  python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/compute_quality_benchmark_data.py \
+   -e $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/build-release/quality_benchmark -i $6 -a $3 \
+   > $2/Quality/results/$5/$filename.log
+}
+export -f compute_benchmark_data
+
+# $1: directory containing the alpha wrap project
+# $2: directory containing the input data folder
+# $3: directory containing the output results
+# $4: alpha value
+# $5: timeout value for robustness benchmark in seconds
+# $6: number of virtual thread used
+# $7: hash of the latest commit
+# $8: hash of a commit to perform the diff with latest
+cd $1
+
+mkdir -p $3/Robustness/results/$7
+mkdir -p $3/Performance/results/$7
+mkdir -p $3/Quality/results/$7
+mkdir -p $3/Robustness/charts_data
+mkdir -p $3/Performance/charts_data
+mkdir -p $3/Quality/charts_data
+mkdir -p $3/Robustness/charts
+mkdir -p $3/Performance/charts
+mkdir -p $3/Quality/charts
+mkdir -p $3/Robustness/log
+mkdir -p $3/Performance/log
+mkdir -p $3/Quality/log
+mkdir -p $3/charts
+
+find $2 -mindepth 1 | parallel -j$6 compute_benchmark_data $1 $3 $4 $5 $7 :::
+
+python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Robustness/generate_robustness_benchmark_charts.py -i $3/Robustness/results/$7 -o $3/Robustness -a $4 -c $7
+
+if [ -z "$8" ]; then
+  python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/generate_performance_benchmark_charts.py -i $3/Performance/results/$7 -o $3/Performance -a $4 -c $7;
+else
+  python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Performance/generate_performance_benchmark_charts.py -i $3/Performance/results/$7 -o $3/Performance -a $4 -c $7 -p $3/Performance/results/$8 -d $8;
+fi
+
+if [ -z "$8" ]; then
+  python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/generate_quality_benchmark_charts.py -i $3/Quality/results/$7 -o $3/Quality -a $4 -c $7;
+else
+  python3 $1/Alpha_wrap_3/benchmark/Alpha_wrap_3/Quality/generate_quality_benchmark_charts.py -i $3/Quality/results/$7 -o $3/Quality -a $4 -c $7 -p $3/Quality/results/$8 -d $8;
+fi
+
+charts_path="$(ls "$3/Robustness/charts"/* -dArt | tail -n 1) $(ls "$3/Performance/charts"/* -dArt | tail -n 2) $(ls "$3/Quality/charts"/* -dArt | tail -n 9)"
+
+pdfjam --nup 2x6 $charts_path --outfile $3/charts/results_$7_$8_alpha_$4_$(date '+%Y-%m-%d_%H:%M:%S').pdf

Alpha_wrap_3/examples/Alpha_wrap_3/CMakeLists.txt

@@ -13,3 +13,5 @@ create_single_source_cgal_program("point_set_wrap.cpp")
 create_single_source_cgal_program("wrap_from_cavity.cpp")
 create_single_source_cgal_program("mixed_inputs_wrap.cpp")
 create_single_source_cgal_program("volumetric_wrap.cpp")
+create_single_source_cgal_program("successive_wraps.cpp")
+create_single_source_cgal_program("pause_and_resume_wrapping.cpp")

Alpha_wrap_3/examples/Alpha_wrap_3/mixed_inputs_wrap.cpp

@@ -103,10 +103,10 @@ int main(int argc, char** argv)
   oracle.add_segment_soup(segments, CGAL::parameters::default_values());
   oracle.add_point_set(ps_points, CGAL::parameters::default_values());
 
-  CGAL::Alpha_wraps_3::internal::Alpha_wrap_3<Oracle> aw3(oracle);
+  CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle> aw3(oracle);
 
-  Mesh output_mesh;
-  aw3(alpha, offset, output_mesh);
+  Mesh wrap;
+  aw3(alpha, offset, wrap);
 
   t.stop();
   std::cout << "Took " << t.time() << std::endl;
@@ -120,10 +120,11 @@ int main(int argc, char** argv)
   std::string ps_name = std::string(ps_filename);
   ps_name = ps_name.substr(ps_name.find_last_of("/") + 1, ps_name.length() - 1);
   ps_name = ps_name.substr(0, ps_name.find_last_of("."));
-  std::string output_name = ts_name + "_" + ss_name + "_"  + ps_name + "_" + std::to_string(static_cast<int>(relative_alpha))
-                            + "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
+  std::string output_name = ts_name + "_" + ss_name + "_"  + ps_name + "_"
+                                  + std::to_string(static_cast<int>(relative_alpha)) + "_"
+                                  + std::to_string(static_cast<int>(relative_offset)) + ".off";
   std::cout << "Writing to " << output_name << std::endl;
-  CGAL::IO::write_polygon_mesh(output_name, output_mesh, CGAL::parameters::stream_precision(17));
+  CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
 
   return EXIT_SUCCESS;
 }

Alpha_wrap_3/examples/Alpha_wrap_3/output_helper.h

@@ -0,0 +1,19 @@
+#ifndef CGAL_ALPHA_WRAP_3_EXAMPLES_OUTPUT_HELPER_H
+#define CGAL_ALPHA_WRAP_3_EXAMPLES_OUTPUT_HELPER_H
+
+#include <string>
+
+std::string generate_output_name(std::string input_name,
+                                 const double alpha,
+                                 const double offset)
+{
+  input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);
+  input_name = input_name.substr(0, input_name.find_last_of("."));
+  std::string output_name = input_name
+                            + "_" + std::to_string(static_cast<int>(alpha))
+                            + "_" + std::to_string(static_cast<int>(offset)) + ".off";
+
+  return output_name;
+}
+
+#endif // CGAL_ALPHA_WRAP_3_EXAMPLES_OUTPUT_HELPER_H

Alpha_wrap_3/examples/Alpha_wrap_3/pause_and_resume_wrapping.cpp

@@ -0,0 +1,164 @@
+// This example demonstrates how to interrupt the wrapping process before it has terminated,
+// and how to resume afterwards.
+//
+// -------------------------------- !! Warning !! --------------------------------------------------
+// By default, the wrapper uses an unsorted LIFO queue of faces to refine. This means that
+// the intermediate result is not very useful because the algorithm carves deep and not wide
+// (somewhat like a DFS vs a BFS).
+//
+// The sorted queue option is enabled with the macro below to make the refinement algorithm
+// more uniform. The downside is that it is slower.
+// -------------------------------------------------------------------------------------------------
+#define CGAL_AW3_USE_SORTED_PRIORITY_QUEUE
+
+#include "output_helper.h"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/alpha_wrap_3.h>
+#include <CGAL/IO/polygon_soup_io.h>
+#include <CGAL/Polygon_mesh_processing/measure.h>
+#include <CGAL/Random.h>
+#include <CGAL/Timer.h>
+
+#include <iostream>
+#include <string>
+
+namespace AW3 = CGAL::Alpha_wraps_3;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point_3 = K::Point_3;
+
+using Points = std::vector<Point_3>;
+using Face = std::array<std::size_t, 3>;
+using Faces = std::vector<Face>;
+
+using Mesh = CGAL::Surface_mesh<Point_3>;
+using face_descriptor = boost::graph_traits<Mesh>::face_descriptor;
+
+struct Interrupter_visitor
+  : public AW3::internal::Wrapping_default_visitor
+{
+  using Base = AW3::internal::Wrapping_default_visitor;
+
+  CGAL::Real_timer timer;
+  double max_time = -1; // in seconds
+
+public:
+  void set_max_time(double t) { max_time = t; }
+
+public:
+  template <typename AlphaWrapper>
+  void on_flood_fill_begin(const AlphaWrapper&)
+  {
+    std::cout << "Starting timer..." << std::endl;
+    timer.start();
+  }
+
+  template <typename Wrapper>
+  bool go_further(const Wrapper&)
+  {
+    if(timer.time() > max_time)
+    {
+      timer.stop();
+      std::cout << "Paused after " << timer.time() << " s." << std::endl;
+      return false;
+    }
+
+    return true;
+  }
+};
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  CGAL::Random rng;
+  std::cout << "Random seed = " << rng.get_seed() << std::endl;
+
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off");
+
+  // = read the soup
+  Points points;
+  Faces faces;
+  if(!CGAL::IO::read_polygon_soup(filename, points, faces) || faces.empty())
+  {
+    std::cerr << "Invalid soup input: " << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "Input: " << points.size() << " points, " << faces.size() << " faces" << std::endl;
+
+  // Compute the alpha and offset values
+  const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : rng.get_double(150., 200.);
+  const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;
+  std::cout << "relative_alpha = " << relative_alpha << std::endl;
+
+  CGAL::Bbox_3 bbox;
+  for(const Point_3& p : points)
+    bbox += p.bbox();
+
+  const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                                       CGAL::square(bbox.ymax() - bbox.ymin()) +
+                                       CGAL::square(bbox.zmax() - bbox.zmin()));
+
+  const double alpha = diag_length / relative_alpha;
+  const double offset = diag_length / relative_offset;
+
+  // Build the wrapper
+  using Oracle = CGAL::Alpha_wraps_3::internal::Triangle_soup_oracle<K>;
+  Oracle oracle(alpha);
+  oracle.add_triangle_soup(points, faces, CGAL::parameters::default_values());
+  CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle> aw3(oracle);
+
+  // --- Launch the wrapping, and pause when the algorithm has spent 1s flooding
+  Interrupter_visitor interrupter;
+  interrupter.set_max_time(1.);
+
+  Mesh wrap;
+  aw3(alpha, offset, wrap, CGAL::parameters::visitor(interrupter));
+  std::cout << ">>> The current wrap has " << num_vertices(wrap) << " vertices" << std::endl;
+  CGAL::IO::write_polygon_mesh("stopped_1.off", wrap, CGAL::parameters::stream_precision(17));
+
+  // --- Restart from the previous state, and pause a bit further
+  interrupter.set_max_time(2.);
+  aw3(alpha, offset, wrap, CGAL::parameters::visitor(interrupter)
+                                            .refine_triangulation(true));
+  std::cout << ">>> The current wrap has " << num_vertices(wrap) << " vertices" << std::endl;
+  CGAL::IO::write_polygon_mesh("stopped_2.off", wrap, CGAL::parameters::stream_precision(17));
+
+  // --- Restart from the previous state, and let it finish
+  aw3(alpha, offset, wrap, CGAL::parameters::refine_triangulation(true));
+  std::cout << ">>> The final (resumed) wrap has " << num_vertices(wrap) << " vertices" << std::endl;
+  std::string output_name = generate_output_name(filename, relative_alpha, relative_offset);
+  std::cout << "Writing to " << "resumed_" + output_name << std::endl;
+  CGAL::IO::write_polygon_mesh("resumed_" + output_name, wrap, CGAL::parameters::stream_precision(17));
+
+  // --- Get the final wrap, in one go:
+  Mesh single_pass_wrap;
+  CGAL::alpha_wrap_3(points, faces, alpha, offset, single_pass_wrap);
+  std::cout << ">>> The final (from scratch) wrap has " << num_vertices(single_pass_wrap) << " vertices" << std::endl;
+
+  output_name = generate_output_name(filename, relative_alpha, relative_offset);
+  std::cout << "Writing to " << output_name << std::endl;
+  CGAL::IO::write_polygon_mesh(output_name, single_pass_wrap, CGAL::parameters::stream_precision(17));
+
+  // --- Compare the results to ensure both approaches yield identical meshes
+  std::vector<std::pair<face_descriptor, face_descriptor> > common;
+  std::vector<face_descriptor> m1_only;
+  std::vector<face_descriptor> m2_only;
+  PMP::match_faces(wrap, single_pass_wrap,
+                   std::back_inserter(common),
+                   std::back_inserter(m1_only),
+                   std::back_inserter(m2_only));
+  if(!m1_only.empty() || !m2_only.empty())
+  {
+    std::cerr << "Error: The two wraps should have been identical!" << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  return EXIT_SUCCESS;
+}

Alpha_wrap_3/examples/Alpha_wrap_3/point_set_wrap.cpp

@@ -1,3 +1,5 @@
+#include "output_helper.h"
+
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
 
@@ -23,7 +25,7 @@ int main(int argc, char** argv)
   Point_container points;
   if(!CGAL::IO::read_points(filename, std::back_inserter(points)) || points.empty())
   {
-    std::cerr << "Invalid input." << std::endl;
+    std::cerr << "Invalid input:" << filename << std::endl;
     return EXIT_FAILURE;
   }
 
@@ -53,11 +55,7 @@ int main(int argc, char** argv)
   std::cout << "Took " << t.time() << " s." << std::endl;
 
   // Save the result
-  std::string input_name = std::string(filename);
-  input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);
-  input_name = input_name.substr(0, input_name.find_last_of("."));
-  std::string output_name = input_name + "_" + std::to_string(static_cast<int>(relative_alpha))
-                            + "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
+  const std::string output_name = generate_output_name(filename, relative_alpha, relative_offset);
   std::cout << "Writing to " << output_name << std::endl;
   CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
 

Alpha_wrap_3/examples/Alpha_wrap_3/successive_wraps.cpp

@@ -0,0 +1,135 @@
+// In this example, we reuse the underlying triangulation of the previous state, and carve using
+// a new (smaller) alpha value. This enables considerable speed-up: the cumulated time taken
+// to run `n` successive instances of `{alpha_wrap(alpha_i)}_(i=1...n)` will be roughly equal
+// to the time taken to the single instance of alpha_wrap(alpha_n) from scratch.
+//
+// The speed-up increases with the number of intermediate results, and on the gap between
+// alpha values: if alpha_2 is close to alpha_1, practically no new computation are required,
+// and the speed-up is almost 100%.
+//
+// -------------------------------- !! Warning !! --------------------------------------------------
+// The result of:
+//   > alpha_wrap(alpha_1, ...)
+//   > alpha_wrap(alpha_2, ..., reuse)
+// is not exactly identical to calling directly:
+//   > alpha_wrap(alpha_2, ..., do_not_reuse)
+// because the queues are sorted slightly differently and the AABB tree is rebuilt differently
+// to optimize the runtime.
+// -------------------------------------------------------------------------------------------------
+
+#include "output_helper.h"
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/alpha_wrap_3.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+#include <CGAL/Real_timer.h>
+
+#include <iostream>
+#include <string>
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using FT = K::FT;
+using Point_3 = K::Point_3;
+
+using Mesh = CGAL::Surface_mesh<Point_3>;
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  // Read the input
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cube.off");
+  std::cout << "Reading " << filename << "..." << std::endl;
+
+  Mesh mesh;
+  if(!PMP::IO::read_polygon_mesh(filename, mesh) || is_empty(mesh) || !is_triangle_mesh(mesh))
+  {
+    std::cerr << "Invalid input:" << filename << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "Input: " << num_vertices(mesh) << " vertices, " << num_faces(mesh) << " faces" << std::endl;
+
+  const CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(mesh);
+  const double diag_length = std::sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
+                                       CGAL::square(bbox.ymax() - bbox.ymin()) +
+                                       CGAL::square(bbox.zmax() - bbox.zmin()));
+
+  // We want decreasing alphas, and these are relative ratios, so they need to be increasing
+  const std::vector<FT> relative_alphas = { 1, 50, 100, 150, 200, 250 };
+  const FT relative_offset = 600;
+
+  // ===============================================================================================
+  // Naive approach:
+
+  CGAL::Real_timer t;
+  double total_time = 0.;
+
+  for(std::size_t i=0; i<relative_alphas.size(); ++i)
+  {
+    t.reset();
+    t.start();
+
+    const double alpha = diag_length / relative_alphas[i];
+    const double offset = diag_length / relative_offset;
+    std::cout << ">>> [" << i << "] alpha: " << alpha << " offset: " << offset << std::endl;
+
+    Mesh wrap;
+    CGAL::alpha_wrap_3(mesh, alpha, offset, wrap);
+
+    t.stop();
+    std::cout << "  Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;
+    std::cout << "  Elapsed time: " << t.time() << " s." << std::endl;
+
+    total_time += t.time();
+  }
+
+  std::cout << "Total elapsed time (naive): " << total_time << " s.\n" << std::endl;
+
+  // ===============================================================================================
+  // Re-use approach
+
+  total_time = 0.;
+  t.reset();
+
+  using Oracle = CGAL::Alpha_wraps_3::internal::Triangle_mesh_oracle<K>;
+  using Wrapper = CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle>;
+  Wrapper wrapper; // contains the triangulation that is being refined iteratively
+
+  for(std::size_t i=0; i<relative_alphas.size(); ++i)
+  {
+    t.reset();
+    t.start();
+
+    const double alpha = diag_length / relative_alphas[i];
+    const double offset = diag_length / relative_offset;
+    std::cout << ">>> [" << i << "] alpha: " << alpha << " offset: " << offset << std::endl;
+
+    // The triangle mesh oracle should be initialized with alpha to internally perform a split
+    // of too-big facets while building the AABB Tree. This split in fact yields a significant
+    // speed-up for meshes with elements that are large compared to alpha. This speed-up makes it
+    // faster to re-build the AABB tree for every value of alpha than to use a non-optimized tree.
+    Oracle oracle(alpha);
+    oracle.add_triangle_mesh(mesh, CGAL::parameters::default_values());
+    wrapper.oracle() = oracle;
+
+    Mesh wrap;
+    wrapper(alpha, offset, wrap, CGAL::parameters::refine_triangulation((i != 0)));
+
+    t.stop();
+    std::cout << "  Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;
+    std::cout << "  Elapsed time: " << t.time() << " s." << std::endl;
+
+    total_time += t.time();
+  }
+
+  std::cout << "Total elapsed time (successive): " << total_time << " s." << std::endl;
+
+  return EXIT_SUCCESS;
+}

Alpha_wrap_3/examples/Alpha_wrap_3/triangle_mesh_wrap.cpp

@@ -1,3 +1,5 @@
+#include "output_helper.h"
+
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
 
@@ -25,7 +27,7 @@ int main(int argc, char** argv)
   Mesh mesh;
   if(!PMP::IO::read_polygon_mesh(filename, mesh) || is_empty(mesh) || !is_triangle_mesh(mesh))
   {
-    std::cerr << "Invalid input." << std::endl;
+    std::cerr << "Invalid input:" << filename << std::endl;
     return EXIT_FAILURE;
   }
 
@@ -56,12 +58,7 @@ int main(int argc, char** argv)
   std::cout << "Took " << t.time() << " s." << std::endl;
 
   // Save the result
-  std::string input_name = std::string(filename);
-  input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);
-  input_name = input_name.substr(0, input_name.find_last_of("."));
-  std::string output_name = input_name
-                            + "_" + std::to_string(static_cast<int>(relative_alpha))
-                            + "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
+  const std::string output_name = generate_output_name(filename, relative_alpha, relative_offset);
   std::cout << "Writing to " << output_name << std::endl;
   CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
 

Alpha_wrap_3/examples/Alpha_wrap_3/triangle_soup_wrap.cpp

@@ -1,3 +1,5 @@
+#include "output_helper.h"
+
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
 
@@ -30,7 +32,7 @@ int main(int argc, char** argv)
   std::vector<std::array<std::size_t, 3> > faces;
   if(!CGAL::IO::read_polygon_soup(filename, points, faces) || faces.empty())
   {
-    std::cerr << "Invalid input." << std::endl;
+    std::cerr << "Invalid input:" << filename << std::endl;
     return EXIT_FAILURE;
   }
 
@@ -63,12 +65,7 @@ int main(int argc, char** argv)
   std::cout << "Took " << t.time() << " s." << std::endl;
 
   // Save the result
-  std::string input_name = std::string(filename);
-  input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);
-  input_name = input_name.substr(0, input_name.find_last_of("."));
-  std::string output_name = input_name
-                            + "_" + std::to_string(static_cast<int>(relative_alpha))
-                            + "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
+  const std::string output_name = generate_output_name(filename, relative_alpha, relative_offset);
   std::cout << "Writing to " << output_name << std::endl;
   CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
 

Alpha_wrap_3/examples/Alpha_wrap_3/volumetric_wrap.cpp

@@ -1,3 +1,5 @@
+#include "output_helper.h"
+
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
 
@@ -70,7 +72,7 @@ int main(int argc, char** argv)
   Faces faces;
   if(!CGAL::IO::read_polygon_soup(filename, points, faces) || faces.empty())
   {
-    std::cerr << "Invalid input." << std::endl;
+    std::cerr << "Invalid input:" << filename << std::endl;
     return EXIT_FAILURE;
   }
 
@@ -101,7 +103,7 @@ int main(int argc, char** argv)
   Oracle oracle(K{});
   oracle.add_triangle_soup(points, faces, CGAL::parameters::default_values());
 
-  CGAL::Alpha_wraps_3::internal::Alpha_wrap_3<Oracle, Delaunay_triangulation> aw3(oracle);
+  CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle, Delaunay_triangulation> aw3(oracle);
   Mesh wrap;
   aw3(alpha, offset, wrap);
 
@@ -113,12 +115,7 @@ int main(int argc, char** argv)
   auto dt = aw3.triangulation();
 
   // Save the result
-  std::string input_name = std::string(filename);
-  input_name = input_name.substr(input_name.find_last_of("/") + 1, input_name.length() - 1);
-  input_name = input_name.substr(0, input_name.find_last_of("."));
-  std::string output_name = input_name
-                            + "_" + std::to_string(static_cast<int>(relative_alpha))
-                            + "_" + std::to_string(static_cast<int>(relative_offset)) + ".off";
+  const std::string output_name = generate_output_name(filename, relative_alpha, relative_offset);
   std::cout << "Writing to " << output_name << std::endl;
   CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
 

Alpha_wrap_3/examples/Alpha_wrap_3/wrap_from_cavity.cpp

@@ -1,3 +1,5 @@
+#include "output_helper.h"
+
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
 
@@ -25,13 +27,13 @@ int main(int argc, char** argv)
   if(!PMP::IO::read_polygon_mesh(filename, input) ||
      is_empty(input) || !is_triangle_mesh(input))
   {
-    std::cerr << "Invalid input." << std::endl;
+    std::cerr << "Invalid input:" << filename << std::endl;
     return EXIT_FAILURE;
   }
 
   std::cout << "Input: " << num_vertices(input) << " vertices, " << num_faces(input) << " faces" << std::endl;
 
-  const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 30.;
+  const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 40.;
   const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;
 
   // Compute the alpha and offset values

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

@@ -20,18 +20,26 @@ namespace CGAL {
 namespace Alpha_wraps_3 {
 namespace internal {
 
+enum class Cell_label
+{
+  // Cells that have been carved
+  OUTSIDE,
+  // Cells that have not yet been carved
+  INSIDE,
+  // OUTSIDE cells that have been labeled "inside" again as to make the result manifold
+  MANIFOLD
+};
+
 template < typename GT,
            typename Cb = CGAL::Delaunay_triangulation_cell_base_with_circumcenter_3<GT> >
 class Alpha_wrap_triangulation_cell_base_3
   : public Cb
 {
-private:
-  bool outside = false;
-
 public:
   typedef typename Cb::Vertex_handle                   Vertex_handle;
   typedef typename Cb::Cell_handle                     Cell_handle;
 
+public:
   template < typename TDS2 >
   struct Rebind_TDS
   {
@@ -39,6 +47,14 @@ public:
     using Other = Alpha_wrap_triangulation_cell_base_3<GT, Cb2>;
   };
 
+private:
+  Cell_label m_label = Cell_label::INSIDE;
+
+#ifndef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE
+  unsigned int m_erase_counter;
+#endif
+
+public:
   Alpha_wrap_triangulation_cell_base_3()
     : Cb()
   {}
@@ -55,8 +71,26 @@ public:
     : Cb(v0, v1, v2, v3, n0, n1, n2, n3)
   {}
 
-  bool is_outside() const { return outside; }
-  bool& is_outside() { return outside; }
+public:
+  Cell_label label() const { return m_label; }
+  void set_label(const Cell_label label) { m_label = label; }
+  bool is_inside() const { return m_label == Cell_label::INSIDE; }
+  bool is_outside() const { return m_label == Cell_label::OUTSIDE; }
+
+#ifndef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE
+  unsigned int erase_counter() const
+  {
+    return m_erase_counter;
+  }
+  void set_erase_counter(unsigned int c)
+  {
+    m_erase_counter = c;
+  }
+  void increment_erase_counter()
+  {
+    ++m_erase_counter;
+  }
+#endif
 };
 
 template <typename Cb>

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

@@ -321,7 +321,6 @@ public:
   typename AABB_tree::Bounding_box bbox() const
   {
     CGAL_precondition(!empty());
-
     return tree().bbox();
   }
 

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

@@ -28,6 +28,7 @@
 #include <iostream>
 #include <iterator>
 #include <functional>
+#include <memory>
 #include <vector>
 
 namespace CGAL {
@@ -41,6 +42,7 @@ struct PS_oracle_traits
   using Geom_traits = Alpha_wrap_AABB_geom_traits<GT_>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
 
   using Points = std::vector<typename GT_::Point_3>;
+  using Points_ptr = std::shared_ptr<Points>;
   using PR_iterator = typename Points::const_iterator;
 
   using Primitive = AABB_primitive<PR_iterator,
@@ -69,26 +71,29 @@ public:
 
 private:
   using Points = typename PSOT::Points;
+  using Points_ptr = typename PSOT::Points_ptr;
   using AABB_tree = typename PSOT::AABB_tree;
   using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, CGAL::Default, BaseOracle>;
 
 private:
-  Points m_points;
+  Points_ptr m_points_ptr;
 
 public:
   // Constructors
-  Point_set_oracle()
-    : Oracle_base(BaseOracle(), Base_GT())
-  { }
-
   Point_set_oracle(const BaseOracle& base_oracle,
                    const Base_GT& gt = Base_GT())
     : Oracle_base(base_oracle, gt)
-  { }
+  {
+    m_points_ptr = std::make_shared<Points>();
+  }
 
   Point_set_oracle(const Base_GT& gt,
                    const BaseOracle& base_oracle = BaseOracle())
-    : Oracle_base(base_oracle, gt)
+    : Point_set_oracle(base_oracle, gt)
+  { }
+
+  Point_set_oracle()
+    : Point_set_oracle(BaseOracle(), Base_GT())
   { }
 
 public:
@@ -101,21 +106,27 @@ public:
     if(points.empty())
     {
 #ifdef CGAL_AW3_DEBUG
-      std::cout << "Warning: Input is empty " << std::endl;
+      std::cout << "Warning: Input is empty (PS)" << std::endl;
 #endif
       return;
     }
 
-    const std::size_t old_size = m_points.size();
-    m_points.insert(std::cend(m_points), std::cbegin(points), std::cend(points));
+    const std::size_t old_size = m_points_ptr->size();
+    m_points_ptr->insert(std::cend(*m_points_ptr), std::cbegin(points), std::cend(points));
 
 #ifdef CGAL_AW3_DEBUG
     std::cout << "Insert into AABB tree (points)..." << std::endl;
 #endif
 
-    this->tree().insert(std::next(std::cbegin(m_points), old_size), std::cend(m_points));
+    this->tree().insert(std::next(std::cbegin(*m_points_ptr), old_size), std::cend(*m_points_ptr));
 
-    CGAL_postcondition(this->tree().size() == m_points.size());
+    // Manually constructing it here purely for profiling reasons: if we keep the lazy approach,
+    // it will be done at the first treatment of a facet that needs a Steiner point.
+    // So if one wanted to bench the flood fill runtime, it would be skewed by the time it takes
+    // to accelerate the tree.
+    this->tree().accelerate_distance_queries();
+
+    CGAL_postcondition(this->tree().size() == m_points_ptr->size());
   }
 };
 

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

@@ -28,6 +28,7 @@
 #include <iostream>
 #include <iterator>
 #include <functional>
+#include <memory>
 #include <vector>
 
 namespace CGAL {
@@ -40,7 +41,9 @@ struct SS_oracle_traits
 {
   using Geom_traits = Alpha_wrap_AABB_geom_traits<GT_>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
 
-  using Segments = std::vector<typename GT_::Segment_3>;
+  using Segment = typename GT_::Segment_3;
+  using Segments = std::vector<Segment>;
+  using Segments_ptr = std::shared_ptr<Segments>;
   using SR_iterator = typename Segments::const_iterator;
 
   using Primitive = AABB_primitive<SR_iterator,
@@ -68,27 +71,31 @@ public:
   using Geom_traits = typename SSOT::Geom_traits;
 
 private:
+  using Segment = typename SSOT::Segment;
   using Segments = typename SSOT::Segments;
+  using Segments_ptr = typename SSOT::Segments_ptr;
   using AABB_tree = typename SSOT::AABB_tree;
   using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, CGAL::Default, BaseOracle>;
 
 private:
-  Segments m_segments;
+  Segments_ptr m_segments_ptr;
 
 public:
   // Constructors
-  Segment_soup_oracle()
-    : Oracle_base(BaseOracle(), Base_GT())
-  { }
-
   Segment_soup_oracle(const BaseOracle& base_oracle,
                       const Base_GT& gt = Base_GT())
     : Oracle_base(base_oracle, gt)
-  { }
+  {
+    m_segments_ptr = std::make_shared<Segments>();
+  }
 
   Segment_soup_oracle(const Base_GT& gt,
                       const BaseOracle& base_oracle = BaseOracle())
-    : Oracle_base(base_oracle, gt)
+    : Segment_soup_oracle(base_oracle, gt)
+  { }
+
+  Segment_soup_oracle()
+    : Segment_soup_oracle(BaseOracle(), Base_GT())
   { }
 
 public:
@@ -100,20 +107,40 @@ public:
     if(segments.empty())
     {
 #ifdef CGAL_AW3_DEBUG
-      std::cout << "Warning: Input is empty " << std::endl;
+      std::cout << "Warning: Input is empty (SS)" << std::endl;
 #endif
       return;
     }
 
-    const std::size_t old_size = m_segments.size();
-    m_segments.insert(std::cend(m_segments), std::cbegin(segments), std::cend(segments));
+    typename Geom_traits::Is_degenerate_3 is_degenerate = this->geom_traits().is_degenerate_3_object();
+
+    const std::size_t old_size = m_segments_ptr->size();
+
+    for(const Segment& s : segments)
+    {
+      if(is_degenerate(s))
+      {
+#ifdef CGAL_AW3_DEBUG
+        std::cerr << "Warning: ignoring degenerate segment " << s << std::endl;
+#endif
+        continue;
+      }
+
+      m_segments_ptr->push_back(s);
+    }
 
 #ifdef CGAL_AW3_DEBUG
     std::cout << "Insert into AABB tree (segments)..." << std::endl;
 #endif
-    this->tree().insert(std::next(std::cbegin(m_segments), old_size), std::cend(m_segments));
+    this->tree().insert(std::next(std::cbegin(*m_segments_ptr), old_size), std::cend(*m_segments_ptr));
 
-    CGAL_postcondition(this->tree().size() == m_segments.size());
+    // Manually constructing it here purely for profiling reasons: if we keep the lazy approach,
+    // it will be done at the first treatment of a facet that needs a Steiner point.
+    // So if one wanted to bench the flood fill runtime, it would be skewed by the time it takes
+    // to accelerate the tree.
+    this->tree().accelerate_distance_queries();
+
+    CGAL_postcondition(this->tree().size() == m_segments_ptr->size());
   }
 };
 

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

@@ -135,7 +135,7 @@ public:
     if(is_empty(tmesh))
     {
 #ifdef CGAL_AW3_DEBUG
-      std::cout << "Warning: Input is empty " << std::endl;
+      std::cout << "Warning: Input is empty (TM)" << std::endl;
 #endif
       return;
     }
@@ -153,7 +153,12 @@ public:
     for(face_descriptor f : faces(tmesh))
     {
       if(Polygon_mesh_processing::is_degenerate_triangle_face(f, tmesh, np))
+      {
+#ifdef CGAL_AW3_DEBUG
+        std::cerr << "Warning: ignoring degenerate face " << f << std::endl;
+#endif
         continue;
+      }
 
       const Point_ref p0 = get(vpm, source(halfedge(f, tmesh), tmesh));
       const Point_ref p1 = get(vpm, target(halfedge(f, tmesh), tmesh));
@@ -164,6 +169,12 @@ public:
       Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits());
     }
 
+    // Manually constructing it here purely for profiling reasons: if we keep the lazy approach,
+    // it will be done at the first treatment of a facet that needs a Steiner point.
+    // So if one wanted to bench the flood fill runtime, it would be skewed by the time it takes
+    // to accelerate the tree.
+    this->tree().accelerate_distance_queries();
+
 #ifdef CGAL_AW3_DEBUG
     std::cout << "Tree: " << this->tree().size() << " primitives (" << num_faces(tmesh) << " faces in input)" << std::endl;
 #endif

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

@@ -133,7 +133,7 @@ public:
     if(points.empty() || faces.empty())
     {
 #ifdef CGAL_AW3_DEBUG
-      std::cout << "Warning: Input is empty " << std::endl;
+      std::cout << "Warning: Input is empty (TS)" << std::endl;
 #endif
       return;
     }
@@ -164,11 +164,22 @@ public:
 
       const Triangle_3 tr = triangle(p0, p1, p2);
       if(is_degenerate(tr))
+      {
+#ifdef CGAL_AW3_DEBUG
+        std::cerr << "Warning: ignoring degenerate face " << tr << std::endl;
+#endif
         continue;
+      }
 
       Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits());
     }
 
+    // Manually constructing it here purely for profiling reasons: if we keep the lazy approach,
+    // it will be done at the first treatment of a facet that needs a Steiner point.
+    // So if one wanted to bench the flood fill runtime, it would be skewed by the time it takes
+    // to accelerate the tree.
+    this->tree().accelerate_distance_queries();
+
 #ifdef CGAL_AW3_DEBUG
     std::cout << "Tree: " << this->tree().size() << " primitives (" << faces.size() << " faces in input)" << std::endl;
 #endif
@@ -179,12 +190,31 @@ public:
   void add_triangle_soup(const TriangleRange& triangles,
                          const CGAL_NP_CLASS& /*np*/ = CGAL::parameters::default_values())
   {
+    if(triangles.empty())
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cout << "Warning: Input is empty (TS)" << std::endl;
+#endif
+      return;
+    }
+
+#ifdef CGAL_AW3_DEBUG
+    std::cout << "Insert into AABB Tree (triangles)..." << std::endl;
+#endif
+
     typename Geom_traits::Is_degenerate_3 is_degenerate = this->geom_traits().is_degenerate_3_object();
 
+    Splitter_base::reserve(triangles.size());
+
     for(const Triangle_3& tr : triangles)
     {
       if(is_degenerate(tr))
+      {
+#ifdef CGAL_AW3_DEBUG
+        std::cerr << "Warning: ignoring degenerate triangle " << tr << std::endl;
+#endif
         continue;
+      }
 
       Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits());
     }

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

@@ -27,27 +27,29 @@ namespace CGAL {
 namespace Alpha_wraps_3 {
 namespace internal {
 
+#ifdef CGAL_AW3_USE_SORTED_PRIORITY_QUEUE
+
 // Represents an alpha-traversable facet in the mutable priority queue
-template <typename DT3>
+template <typename Tr>
 class Gate
 {
-  using Facet = typename DT3::Facet;
-  using FT = typename DT3::Geom_traits::FT;
+  using Facet = typename Tr::Facet;
+  using FT = typename Tr::Geom_traits::FT;
 
 private:
   Facet m_facet;
   FT m_priority; // circumsphere sq_radius
-  bool m_is_artificial_facet;
+  bool m_is_permissive_facet;
 
 public:
   // Constructors
   Gate(const Facet& facet,
        const FT& priority,
-       const bool is_artificial_facet)
+       const bool is_permissive_facet)
     :
       m_facet(facet),
       m_priority(priority),
-      m_is_artificial_facet(is_artificial_facet)
+      m_is_permissive_facet(is_permissive_facet)
   {
     CGAL_assertion(priority >= 0);
   }
@@ -60,34 +62,85 @@ public:
 public:
   const Facet& facet() const { return m_facet; }
   const FT& priority() const { return m_priority; }
-  bool is_artificial_facet() const { return m_is_artificial_facet; }
+  bool is_permissive_facet() const { return m_is_permissive_facet; }
 };
 
 struct Less_gate
 {
-  template <typename DT3>
-  bool operator()(const Gate<DT3>& a, const Gate<DT3>& b) const
+  template <typename Tr>
+  bool operator()(const Gate<Tr>& a, const Gate<Tr>& b) const
   {
-    // @fixme? make it a total order by comparing addresses if both gates are bbox facets
-    if(a.is_artificial_facet())
-      return true;
-    else if(b.is_artificial_facet())
-      return false;
+    // If one is permissive and the other is not, give priority to the permissive facet.
+    //
+    // The permissive facet are given highest priority because they need to be treated
+    // regardless of their circumradius. Treating them first allow the part that depends
+    // on alpha to be treated uniformly in a way: whatever the alpha, all permissive faces
+    // will first be treated.
+    if(a.is_permissive_facet() != b.is_permissive_facet())
+      return a.is_permissive_facet();
+
+    if(a.priority() == b.priority())
+    {
+      // arbitrary, the sole purpose is to make it a total order for determinism
+      if(a.facet().first->time_stamp() == b.facet().first->time_stamp())
+        return a.facet().second < b.facet().second;
+
+      return a.facet().first->time_stamp() < b.facet().first->time_stamp();
+    }
+
     return a.priority() > b.priority();
   }
 };
 
-template <typename DT3>
+#else // CGAL_AW3_USE_SORTED_PRIORITY_QUEUE
+
+// Represents an alpha-traversable facet in the mutable priority queue
+template <typename Tr>
+class Gate
+{
+  using Facet = typename Tr::Facet;
+  using FT = typename Tr::Geom_traits::FT;
+
+private:
+  Facet m_facet, m_mirror_facet;
+  const unsigned int m_erase_counter_mem;
+  const unsigned int m_mirror_erase_counter_mem;
+
+public:
+  // Constructors
+  Gate(const Facet& facet,
+       const Tr& tr)
+    :
+      m_facet(facet),
+      m_mirror_facet(tr.mirror_facet(facet)),
+      m_erase_counter_mem(m_facet.first->erase_counter()),
+      m_mirror_erase_counter_mem(m_mirror_facet.first->erase_counter())
+  {
+  }
+
+public:
+  const Facet& facet() const { return m_facet; }
+
+  bool is_zombie() const
+  {
+    return (m_facet.first->erase_counter() != m_erase_counter_mem) ||
+           (m_mirror_facet.first->erase_counter() != m_mirror_erase_counter_mem);
+  }
+};
+
+#endif // CGAL_AW3_USE_SORTED_PRIORITY_QUEUE
+
+template <typename Tr>
 struct Gate_ID_PM
 {
-  using key_type = Gate<DT3>;
+  using key_type = Gate<Tr>;
   using value_type = std::size_t;
   using reference = std::size_t;
   using category = boost::readable_property_map_tag;
 
   inline friend value_type get(Gate_ID_PM, const key_type& k)
   {
-    using Facet = typename DT3::Facet;
+    using Facet = typename Tr::Facet;
 
     const Facet& f = k.facet();
     return (4 * f.first->time_stamp() + f.second);

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

@@ -40,16 +40,16 @@ struct Orientation_of_circumcenter
   }
 };
 
-template <typename Dt>
+template <typename Tr>
 bool
-less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha,
-                                        const typename Dt::Facet& fh,
-                                        const Dt& dt)
+less_squared_radius_of_min_empty_sphere(typename Tr::Geom_traits::FT sq_alpha,
+                                        const typename Tr::Facet& fh,
+                                        const Tr& tr)
 {
-  using Cell_handle = typename Dt::Cell_handle;
-  using Point = typename Dt::Point;
+  using Cell_handle = typename Tr::Cell_handle;
+  using Point = typename Tr::Point;
 
-  using CK = typename Dt::Geom_traits;
+  using CK = typename Tr::Geom_traits;
   using Exact_kernel = typename Exact_kernel_selector<CK>::Exact_kernel;
   using Approximate_kernel = Simple_cartesian<Interval_nt_advanced>;
   using C2A = Cartesian_converter<CK, Approximate_kernel>;
@@ -61,21 +61,30 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha,
 
   Orientation_of_circumcenter orientation_of_circumcenter;
 
+#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
+  std::cout << "Checking for traversability of facet" << std::endl;
+#endif
+
   const Cell_handle c = fh.first;
   const int ic = fh.second;
   const Cell_handle n = c->neighbor(ic);
 
-  const Point& p1 = dt.point(c, Dt::vertex_triple_index(ic,0));
-  const Point& p2 = dt.point(c, Dt::vertex_triple_index(ic,1));
-  const Point& p3 = dt.point(c, Dt::vertex_triple_index(ic,2));
+  const Point& p1 = tr.point(c, Tr::vertex_triple_index(ic,0));
+  const Point& p2 = tr.point(c, Tr::vertex_triple_index(ic,1));
+  const Point& p3 = tr.point(c, Tr::vertex_triple_index(ic,2));
 
   // This is not actually possible in the context of alpha wrapping, but keeping it for genericity
   // and because it does not cost anything.
-  if(dt.is_infinite(n))
+  if(tr.is_infinite(n))
   {
+#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
+    std::cerr << "Warning: computing less_squared_radius_of_min_empty_sphere() with an infinite neighbor?" << std::endl;
+#endif
+    CGAL_assertion(!tr.is_infinite(c));
+
     Orientation ori = orientation_of_circumcenter(p1, p2, p3,
-                                                  dt.point(c, 0), dt.point(c, 1),
-                                                  dt.point(c, 2), dt.point(c, 3));
+                                                  tr.point(c, 0), tr.point(c, 1),
+                                                  tr.point(c, 2), tr.point(c, 3));
 
     if(ori == POSITIVE)
     {
@@ -84,18 +93,22 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha,
     }
     else
     {
-      Comparison_result cr = compare_squared_radius(dt.point(c, 0), dt.point(c, 1),
-                                                    dt.point(c, 2), dt.point(c, 3),
+      Comparison_result cr = compare_squared_radius(tr.point(c, 0), tr.point(c, 1),
+                                                    tr.point(c, 2), tr.point(c, 3),
                                                     sq_alpha);
       return cr == LARGER;
     }
   }
 
-  if(dt.is_infinite(c))
+  if(tr.is_infinite(c))
   {
     Orientation ori = orientation_of_circumcenter(p1, p2, p3,
-                                                  dt.point(n, 0), dt.point(n, 1),
-                                                  dt.point(n, 2), dt.point(n, 3));
+                                                  tr.point(n, 0), tr.point(n, 1),
+                                                  tr.point(n, 2), tr.point(n, 3));
+
+#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
+    std::cout << "Cell 'c' is infinite; Orientation: " << ori << std::endl;
+#endif
 
     if(ori == NEGATIVE)
     {
@@ -104,8 +117,8 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha,
     }
     else
     {
-      Comparison_result cr = compare_squared_radius(dt.point(n, 0), dt.point(n, 1),
-                                                    dt.point(n, 2), dt.point(n, 3),
+      Comparison_result cr = compare_squared_radius(tr.point(n, 0), tr.point(n, 1),
+                                                    tr.point(n, 2), tr.point(n, 3),
                                                     sq_alpha);
       return cr == LARGER;
     }
@@ -113,40 +126,40 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha,
 
   // both c and n are finite
   if(orientation_of_circumcenter(p1, p2, p3,
-                                 dt.point(c, 0), dt.point(c, 1), dt.point(c, 2), dt.point(c, 3)) !=
+                                 tr.point(c, 0), tr.point(c, 1), tr.point(c, 2), tr.point(c, 3)) !=
      orientation_of_circumcenter(p1, p2, p3,
-                                 dt.point(n, 0), dt.point(n, 1), dt.point(n, 2), dt.point(n, 3)))
+                                 tr.point(n, 0), tr.point(n, 1), tr.point(n, 2), tr.point(n, 3)))
   {
     Comparison_result cr = compare_squared_radius(p1, p2, p3, sq_alpha);
 #ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
     std::cout << "dual crosses the face; CR: "
-              << typename Dt::Geom_traits().compute_squared_radius_3_object()(p1, p2, p3)
+              << typename Tr::Geom_traits().compute_squared_radius_3_object()(p1, p2, p3)
               << " sq alpha " << sq_alpha << std::endl;
 #endif
     return cr == LARGER;
   }
   else
   {
-    Comparison_result cr = compare_squared_radius(dt.point(c, 0), dt.point(c, 1),
-                                                  dt.point(c, 2), dt.point(c, 3),
+    Comparison_result cr = compare_squared_radius(tr.point(c, 0), tr.point(c, 1),
+                                                  tr.point(c, 2), tr.point(c, 3),
                                                   sq_alpha);
 #ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
     std::cout << "dual does not cross the face; CR(c): "
-              << typename Dt::Geom_traits().compute_squared_radius_3_object()(dt.point(c, 0), dt.point(c, 1),
-                                                                              dt.point(c, 2), dt.point(c, 3))
+              << typename Tr::Geom_traits().compute_squared_radius_3_object()(tr.point(c, 0), tr.point(c, 1),
+                                                                              tr.point(c, 2), tr.point(c, 3))
               << " sq alpha " << sq_alpha << std::endl;
 #endif
 
     if(cr != LARGER)
       return false;
 
-    cr = compare_squared_radius(dt.point(n, 0), dt.point(n, 1),
-                                dt.point(n, 2), dt.point(n, 3),
+    cr = compare_squared_radius(tr.point(n, 0), tr.point(n, 1),
+                                tr.point(n, 2), tr.point(n, 3),
                                 sq_alpha);
 #ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
     std::cout << "dual does not cross the face; CR(n): "
-              << typename Dt::Geom_traits().compute_squared_radius_3_object()(dt.point(n, 0), dt.point(n, 1),
-                                                                              dt.point(n, 2), dt.point(n, 3))
+              << typename Tr::Geom_traits().compute_squared_radius_3_object()(tr.point(n, 0), tr.point(n, 1),
+                                                                              tr.point(n, 2), tr.point(n, 3))
               << " sq alpha " << sq_alpha << std::endl;
 #endif
 
@@ -154,6 +167,100 @@ less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha,
   }
 }
 
+template <typename Tr>
+typename Tr::Geom_traits::FT
+smallest_squared_radius_3(const typename Tr::Facet& fh,
+                          const Tr& tr)
+{
+  using Cell_handle = typename Tr::Cell_handle;
+  using Point = typename Tr::Point;
+  using FT = typename Tr::Geom_traits::FT;
+
+  using CK = typename Tr::Geom_traits;
+  using Exact_kernel = typename Exact_kernel_selector<CK>::Exact_kernel;
+  using Approximate_kernel = Simple_cartesian<Interval_nt_advanced>;
+  using C2A = Cartesian_converter<CK, Approximate_kernel>;
+  using C2E = typename Exact_kernel_selector<CK>::C2E;
+
+  using Orientation_of_circumcenter = Filtered_predicate<Orientation_of_circumcenter<Exact_kernel>,
+                                                         Orientation_of_circumcenter<Approximate_kernel>,
+                                                         C2E, C2A>;
+
+  Orientation_of_circumcenter orientation_of_circumcenter;
+
+  auto squared_radius = tr.geom_traits().compute_squared_radius_3_object();
+
+#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
+  std::cout << "Computing circumradius of facet" << std::endl;
+#endif
+
+  CGAL_precondition(!tr.is_infinite(fh));
+
+  const Cell_handle c = fh.first;
+  const int ic = fh.second;
+  const Cell_handle n = c->neighbor(ic);
+
+  const Point& p1 = tr.point(c, Tr::vertex_triple_index(ic,0));
+  const Point& p2 = tr.point(c, Tr::vertex_triple_index(ic,1));
+  const Point& p3 = tr.point(c, Tr::vertex_triple_index(ic,2));
+
+  // This is not actually possible in the context of alpha wrapping, but keeping it for genericity
+  // and because it does not cost anything.
+  if(tr.is_infinite(n))
+  {
+    CGAL_assertion(!tr.is_infinite(c));
+
+    Orientation ori = orientation_of_circumcenter(p1, p2, p3,
+                                                  tr.point(c, 0), tr.point(c, 1),
+                                                  tr.point(c, 2), tr.point(c, 3));
+    if(ori == POSITIVE)
+      return squared_radius(p1, p2, p3);
+    else
+      return squared_radius(tr.point(c, 0), tr.point(c, 1), tr.point(c, 2), tr.point(c, 3));
+  }
+
+  if(tr.is_infinite(c))
+  {
+    Orientation ori = orientation_of_circumcenter(p1, p2, p3,
+                                                  tr.point(n, 0), tr.point(n, 1),
+                                                  tr.point(n, 2), tr.point(n, 3));
+
+#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
+    std::cout << "Cell 'c' is infinite; Orientation: " << ori << std::endl;
+#endif
+
+    if(ori == NEGATIVE)
+      return squared_radius(p1, p2, p3);
+    else
+      return squared_radius(tr.point(n, 0), tr.point(n, 1), tr.point(n, 2), tr.point(n, 3));
+  }
+
+  // both c and n are finite
+  if(orientation_of_circumcenter(p1, p2, p3,
+                                 tr.point(c, 0), tr.point(c, 1), tr.point(c, 2), tr.point(c, 3)) !=
+     orientation_of_circumcenter(p1, p2, p3,
+                                 tr.point(n, 0), tr.point(n, 1), tr.point(n, 2), tr.point(n, 3)))
+  {
+#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
+    std::cout << "dual crosses the face; CR: " << squared_radius(p1, p2, p3) << std::endl;
+#endif
+
+    return squared_radius(p1, p2, p3);
+  }
+  else
+  {
+    const FT cr = squared_radius(tr.point(c, 0), tr.point(c, 1), tr.point(c, 2), tr.point(c, 3));
+    const FT cnr = squared_radius(tr.point(n, 0), tr.point(n, 1), tr.point(n, 2), tr.point(n, 3));
+
+#ifdef CGAL_AW3_DEBUG_TRAVERSABILITY
+    std::cout << "dual does not cross the face; CR(c): " << cr << " CRn: " << cnr << std::endl;
+#endif
+
+    return (CGAL::min)(cr, cnr);
+  }
+}
+
+
 } // namespace internal
 } // namespace Alpha_wraps_3
 } // namespace CGAL

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

@@ -117,7 +117,7 @@ bool has_expected_Hausdorff_distance(const TriangleMesh& wrap,
 
 template <typename TriangleMesh, typename NamedParameters = parameters::Default_named_parameters>
 bool is_valid_wrap(const TriangleMesh& wrap,
-                   const bool check_manifoldness = true,
+                   const bool check_manifoldness,
                    const NamedParameters& np = parameters::default_values())
 {
   namespace PMP = CGAL::Polygon_mesh_processing;
@@ -203,6 +203,13 @@ bool is_valid_wrap(const TriangleMesh& wrap,
   return true;
 }
 
+template <typename TriangleMesh, typename NamedParameters = parameters::Default_named_parameters>
+bool is_valid_wrap(const TriangleMesh& wrap,
+                   const NamedParameters& np = parameters::default_values())
+{
+  return is_valid_wrap(wrap, true /*consider manifoldness*/, np);
+}
+
 template <typename InputTriangleMesh, typename OutputTriangleMesh,
           typename InputNamedParameters = parameters::Default_named_parameters,
           typename OutputNamedParameters = parameters::Default_named_parameters>

Alpha_wrap_3/include/CGAL/alpha_wrap_3.h

@@ -105,7 +105,7 @@ void alpha_wrap_3(const PointRange& points,
   using NP_helper = Point_set_processing_3_np_helper<PointRange, InputNamedParameters>;
   using Geom_traits = typename NP_helper::Geom_traits;
   using Oracle = Alpha_wraps_3::internal::Triangle_soup_oracle<Geom_traits>;
-  using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3<Oracle>;
+  using AW3 = Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle>;
 
   Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
 
@@ -254,7 +254,7 @@ void alpha_wrap_3(const TriangleMesh& tmesh,
 
   using Geom_traits = typename GetGeomTraits<TriangleMesh, InputNamedParameters>::type;
   using Oracle = Alpha_wraps_3::internal::Triangle_mesh_oracle<Geom_traits>;
-  using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3<Oracle>;
+  using AW3 = Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle>;
 
   Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
 
@@ -350,7 +350,7 @@ void alpha_wrap_3(const PointRange& points,
   using NP_helper = Point_set_processing_3_np_helper<PointRange, InputNamedParameters>;
   using Geom_traits = typename NP_helper::Geom_traits;
   using Oracle = Alpha_wraps_3::internal::Point_set_oracle<Geom_traits>;
-  using AW3 = Alpha_wraps_3::internal::Alpha_wrap_3<Oracle>;
+  using AW3 = Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle>;
 
   Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
 

Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp

@@ -1,12 +1,13 @@
 #define CGAL_AW3_TIMER
 #define CGAL_AW3_DEBUG
+#define CGAL_AW3_DEBUG_MANIFOLDNESS
 // #define CGAL_AW3_DEBUG_INITIALIZATION
 
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
 #include <CGAL/alpha_wrap_3.h>
-#include "alpha_wrap_validation.h"
+#include <CGAL/Alpha_wrap_3/internal/validation.h>
 
 #include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
 
@@ -27,7 +28,7 @@ void generate_random_seeds(const Oracle& oracle,
                            Seeds& seeds,
                            CGAL::Random& r)
 {
-  const auto bbox = CGAL::Alpha_wraps_3::internal::Alpha_wrap_3<Oracle>(oracle).construct_bbox(offset);
+  const auto bbox = CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle>(oracle).construct_bbox(offset);
   const double sq_offset = CGAL::square(offset);
 
   while(seeds.size() < 3)
@@ -69,7 +70,7 @@ void alpha_wrap_triangle_mesh(Mesh& input_mesh,
 
   Oracle oracle;
   oracle.add_triangle_mesh(input_mesh);
-  AW3::internal::Alpha_wrap_3<Oracle> aw3(oracle);
+  AW3::internal::Alpha_wrapper_3<Oracle> aw3(oracle);
 
   if(seeds.empty())
     generate_random_seeds(oracle, offset, seeds, r);

Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp

@@ -1,11 +1,12 @@
 #define CGAL_AW3_TIMER
 #define CGAL_AW3_DEBUG
+#define CGAL_AW3_DEBUG_MANIFOLDNESS
 //#define CGAL_AW3_DEBUG_STEINER_COMPUTATION
 //#define CGAL_AW3_DEBUG_INITIALIZATION
 //#define CGAL_AW3_DEBUG_QUEUE
 
 #include <CGAL/alpha_wrap_3.h>
-#include "alpha_wrap_validation.h"
+#include <CGAL/Alpha_wrap_3/internal/validation.h>
 
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp

@@ -1,10 +1,11 @@
 #define CGAL_AW3_TIMER
 #define CGAL_AW3_DEBUG
+#define CGAL_AW3_DEBUG_MANIFOLDNESS
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
 #include <CGAL/alpha_wrap_3.h>
-#include "alpha_wrap_validation.h"
+#include <CGAL/Alpha_wrap_3/internal/validation.h>
 
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/IO/polygon_soup_io.h>
@@ -52,7 +53,7 @@ void alpha_wrap_triangle_soup(Points& pr,
   // AW3
   Oracle oracle;
   oracle.add_triangle_soup(pr, fr);
-  AW3::internal::Alpha_wrap_3<Oracle> aw3(oracle);
+  AW3::internal::Alpha_wrapper_3<Oracle> aw3(oracle);
 
   Mesh wrap;
   aw3(alpha, offset, wrap, CGAL::parameters::do_enforce_manifoldness(false));

Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp

@@ -1,12 +1,12 @@
 #define CGAL_AW3_TIMER
 //#define CGAL_AW3_DEBUG
-//#define CGAL_AW3_DEBUG_MANIFOLDNESS
+#define CGAL_AW3_DEBUG_MANIFOLDNESS
 //#define CGAL_AW3_DEBUG_STEINER_COMPUTATION
 //#define CGAL_AW3_DEBUG_INITIALIZATION
 //#define CGAL_AW3_DEBUG_QUEUE
 
 #include <CGAL/alpha_wrap_3.h>
-#include "alpha_wrap_validation.h"
+#include <CGAL/Alpha_wrap_3/internal/validation.h>
 
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

Arrangement_on_surface_2/include/CGAL/Arr_batched_point_location.h

@@ -26,7 +26,6 @@
 #include <CGAL/Surface_sweep_2/Arr_batched_pl_ss_visitor.h>
 
 #include <vector>
-#include <boost/mpl/if.hpp>
 #include <boost/type_traits.hpp>
 
 namespace CGAL {
@@ -120,7 +119,7 @@ locate(const Arrangement_on_surface_2<GeometryTraits_2, TopologyTraits>& arr,
    * Use the form 'A a(*b);' and not ''A a = b;' to handle the case where A has
    * only an implicit constructor, (which takes *b as a parameter).
    */
-  typename boost::mpl::if_<std::is_same<Gt2, Bgt2>, const Bgt2&, Bgt2>::type
+  std::conditional_t<std::is_same_v<Gt2, Bgt2>, const Bgt2&, Bgt2>
     ex_traits(*geom_traits);
 
   // Define the sweep-line visitor and perform the sweep.

Arrangement_on_surface_2/include/CGAL/Arr_overlay_2.h

@@ -25,7 +25,6 @@
 
 #include <vector>
 #include <optional>
-#include <boost/mpl/if.hpp>
 #include <boost/mpl/or.hpp>
 #include <boost/type_traits.hpp>
 
@@ -247,8 +246,8 @@ overlay(const Arrangement_on_surface_2<GeometryTraitsA_2, TopologyTraitsA>& arr1
    * Use the form 'A a(*b);' and not ''A a = b;' to handle the case where A has
    * only an implicit constructor, (which takes *b as a parameter).
    */
-  typename boost::mpl::if_<std::is_same<Gt_adaptor_2, Ovl_gt2>,
-                           const Ovl_gt2&, Ovl_gt2>::type
+  std::conditional_t<std::is_same_v<Gt_adaptor_2, Ovl_gt2>,
+                     const Ovl_gt2&, Ovl_gt2>
     ex_traits(*traits_adaptor);
 
   Ovl_visitor visitor(&arr1, &arr2, &arr, &ovl_tr);

Arrangement_on_surface_2/include/CGAL/Arr_tags.h

@@ -21,7 +21,6 @@
 #include <CGAL/config.h>
 #include <boost/type_traits.hpp>
 #include <boost/mpl/bool.hpp>
-#include <boost/mpl/if.hpp>
 #include <boost/mpl/and.hpp>
 #include <boost/mpl/or.hpp>
 #include <boost/mpl/not.hpp>
@@ -246,14 +245,14 @@ struct Arr_all_sides_not_finite_tag :
 struct Arr_not_all_sides_not_finite_tag :
   public virtual Arr_not_all_sides_oblivious_tag {};
 
-typedef boost::mpl::bool_<true>                                 Arr_true;
-typedef boost::mpl::bool_<false>                                Arr_false;
+typedef std::true_type                                          Arr_true;
+typedef std::false_type                                         Arr_false;
 
 template <typename ArrSideCategory>
 struct Arr_is_side_oblivious {
   typedef ArrSideCategory                                       Side_cat;
   typedef std::is_same<Side_cat, Arr_oblivious_side_tag>        Is_same;
-  typedef boost::mpl::if_<Is_same, Arr_true, Arr_false>         result;
+  typedef std::bool_constant<Is_same::value>                    result;
   typedef typename result::type                                 type;
 };
 
@@ -261,7 +260,7 @@ template <typename ArrSideCategory>
 struct Arr_is_side_open {
   typedef ArrSideCategory                                       Side_cat;
   typedef std::is_same<Side_cat, Arr_open_side_tag>             Is_same;
-  typedef boost::mpl::if_<Is_same, Arr_true, Arr_false>         result;
+  typedef std::bool_constant<Is_same::value>                    result;
   typedef typename result::type                                 type;
 };
 
@@ -269,15 +268,19 @@ template <typename ArrSideCategory>
 struct Arr_is_side_identified {
   typedef ArrSideCategory                                       Side_cat;
   typedef std::is_same<Side_cat, Arr_identified_side_tag>       Is_same;
-  typedef boost::mpl::if_<Is_same, Arr_true, Arr_false>         result;
+  typedef std::bool_constant<Is_same::value>                    result;
   typedef typename result::type                                 type;
 };
 
+template <typename ArrSideCategory>
+inline constexpr bool Arr_is_side_identified_v =
+  Arr_is_side_identified<ArrSideCategory>::type::value;
+
 template <typename ArrSideCategory>
 struct Arr_is_side_contracted {
   typedef ArrSideCategory                                       Side_cat;
   typedef std::is_same<Side_cat, Arr_contracted_side_tag>       Is_same;
-  typedef boost::mpl::if_<Is_same, Arr_true, Arr_false>         result;
+  typedef std::bool_constant<Is_same::value>                    result;
   typedef typename result::type                                 type;
 };
 
@@ -285,7 +288,7 @@ template <typename ArrSideCategory>
 struct Arr_is_side_closed {
   typedef ArrSideCategory                                       Side_cat;
   typedef std::is_same<Side_cat, Arr_closed_side_tag>           Is_same;
-  typedef boost::mpl::if_<Is_same, Arr_true, Arr_false>         result;
+  typedef std::bool_constant<Is_same::value>                    result;
   typedef typename result::type                                 type;
 };
 
@@ -307,10 +310,10 @@ struct Arr_all_sides_oblivious_category {
   /*! Boolean tag that is Arr_all_sides_oblivious_tag if all sides are
    * oblivious, otherwise Arr_not_all_sides_oblivious_tag
    */
-  typedef typename boost::mpl::if_<boost::mpl::and_<Lef_obl, Rig_obl,
-                                                    Bot_obl, Top_obl>,
-                                   Arr_all_sides_oblivious_tag,
-                                   Arr_not_all_sides_oblivious_tag>::type
+  typedef std::conditional_t<Lef_obl::value && Rig_obl::value &&
+                             Bot_obl::value && Top_obl::value,
+                             Arr_all_sides_oblivious_tag,
+                             Arr_not_all_sides_oblivious_tag>
     result;
 };
 
@@ -331,19 +334,19 @@ private:
   typedef typename Arr_is_side_open<Bot_side_cat>::result       Bot_ope;
   typedef typename Arr_is_side_open<Top_side_cat>::result       Top_ope;
 
-  typedef boost::mpl::not_<Lef_ope>                             Lef_not_ope;
-  typedef boost::mpl::not_<Rig_ope>                             Rig_not_ope;
-  typedef boost::mpl::not_<Bot_ope>                             Bot_not_ope;
-  typedef boost::mpl::not_<Top_ope>                             Top_not_ope;
+  static inline constexpr bool lef_not_ope = !Lef_ope::value;
+  static inline constexpr bool rig_not_ope = !Rig_ope::value;
+  static inline constexpr bool bot_not_ope = !Bot_ope::value;
+  static inline constexpr bool top_not_ope = !Top_ope::value;
 
 public:
   /*! Boolean tag that is Arr_all_sides_not_open_tag if all sides are not-open,
    * otherwise Arr_not_all_sides_not_open_tag
    */
-  typedef typename boost::mpl::if_<boost::mpl::and_<Lef_not_ope, Rig_not_ope,
-                                                    Bot_not_ope, Top_not_ope>,
-                                   Arr_all_sides_not_open_tag,
-                                   Arr_not_all_sides_not_open_tag>::type
+  typedef std::conditional_t<lef_not_ope && rig_not_ope &&
+                             bot_not_ope && top_not_ope,
+                             Arr_all_sides_not_open_tag,
+                             Arr_not_all_sides_not_open_tag>
     result;
 };
 
@@ -374,23 +377,23 @@ private:
   typedef typename Arr_is_side_open<Bot_side_cat>::result       Bot_ope;
   typedef typename Arr_is_side_open<Top_side_cat>::result       Top_ope;
 
-  typedef boost::mpl::or_<Lef_obl, Lef_ope>                     Lef_obl_or_ope;
-  typedef boost::mpl::or_<Rig_obl, Rig_ope>                     Rig_obl_or_ope;
-  typedef boost::mpl::or_<Bot_obl, Bot_ope>                     Bot_obl_or_ope;
-  typedef boost::mpl::or_<Top_obl, Top_ope>                     Top_obl_or_ope;
+  static inline constexpr bool lef_obl_or_ope = Lef_obl::value || Lef_ope::value;
+  static inline constexpr bool rig_obl_or_ope = Rig_obl::value || Rig_ope::value;
+  static inline constexpr bool bot_obl_or_ope = Bot_obl::value || Bot_ope::value;
+  static inline constexpr bool top_obl_or_ope = Top_obl::value || Top_ope::value;
 
-  typedef typename boost::mpl::if_<boost::mpl::and_<Lef_obl_or_ope,
-                                                    Rig_obl_or_ope,
-                                                    Bot_obl_or_ope,
-                                                    Top_obl_or_ope>,
-                                   Arr_all_sides_not_finite_tag,
-                                   Arr_not_all_sides_not_finite_tag>::type
+  typedef std::conditional_t<lef_obl_or_ope &&
+                             rig_obl_or_ope &&
+                             bot_obl_or_ope &&
+                             top_obl_or_ope,
+                             Arr_all_sides_not_finite_tag,
+                             Arr_not_all_sides_not_finite_tag>
     tmp;
 
 public:
-  typedef typename boost::mpl::if_<boost::mpl::and_<Lef_obl, Rig_obl,
-                                                    Bot_obl, Top_obl>,
-                                   Arr_all_sides_oblivious_tag, tmp>::type
+  typedef std::conditional_t<Lef_obl::value && Rig_obl::value &&
+                             Bot_obl::value && Top_obl::value,
+                             Arr_all_sides_oblivious_tag, tmp>
     result;
 };
 
@@ -404,32 +407,27 @@ struct Arr_sane_identified_tagging {
   typedef ArrBottomSideCategory                         Bot_side_cat;
   typedef ArrTopSideCategory                            Top_side_cat;
 
-  typedef typename Arr_is_side_identified<Lef_side_cat>::result Lef_ide;
-  typedef typename Arr_is_side_identified<Rig_side_cat>::result Rig_ide;
-  typedef typename Arr_is_side_identified<Bot_side_cat>::result Bot_ide;
-  typedef typename Arr_is_side_identified<Top_side_cat>::result Top_ide;
+  static inline constexpr bool lef_ide = Arr_is_side_identified_v<Lef_side_cat>;
+  static inline constexpr bool rig_ide = Arr_is_side_identified_v<Rig_side_cat>;
+  static inline constexpr bool bot_ide = Arr_is_side_identified_v<Bot_side_cat>;
+  static inline constexpr bool top_ide = Arr_is_side_identified_v<Top_side_cat>;
 
-  typedef boost::mpl::and_<Lef_ide, Rig_ide>            LR_ide;
-  typedef boost::mpl::and_<Bot_ide, Top_ide>            BT_ide;
+  static inline constexpr bool lr_ide = lef_ide && rig_ide;
+  static inline constexpr bool bt_ide = bot_ide && top_ide;
 
-  typedef boost::mpl::not_<Lef_ide>                     Lef_not_ide;
-  typedef boost::mpl::not_<Rig_ide>                     Rig_not_ide;
-  typedef boost::mpl::not_<Bot_ide>                     Bot_not_ide;
-  typedef boost::mpl::not_<Top_ide>                     Top_not_ide;
+  static inline constexpr bool lr_not_ide = !lef_ide && !rig_ide;
 
-  typedef boost::mpl::and_<Lef_not_ide, Rig_not_ide>    LR_not_ide;
+  static inline constexpr bool bt_not_ide = !bot_ide && !top_ide;
 
-  typedef boost::mpl::and_<Bot_not_ide, Top_not_ide>    BT_not_ide;
+  static inline constexpr bool lr_ok = lr_ide || lr_not_ide;
+  static inline constexpr bool bt_ok = bt_ide || bt_not_ide;
 
-  typedef boost::mpl::or_<LR_ide, LR_not_ide>           LR_ok;
-  typedef boost::mpl::or_<BT_ide, BT_not_ide>           BT_ok;
-
-  /*! Boolean tag that is bool_<true> if opposite sides are either
+  /*! Boolean tag that is bool_constant<true> if opposite sides are either
    * both identified or both not-identified,
-   * otherwise bool_<false>
+   * otherwise bool_constant<false>
    */
-  typedef boost::mpl::and_<LR_ok, BT_ok>                result;
-  static constexpr bool value = result::value;
+  typedef std::bool_constant<lr_ok && bt_ok>              result;
+  static inline constexpr bool value = result::value;
 };
 
 /*! Checks whether one of two boundary sides are identified
@@ -448,10 +446,10 @@ struct Arr_has_identified_sides {
   typedef typename Arr_is_side_identified<Side_one_cat>::result Side_one_ide;
   typedef typename Arr_is_side_identified<Side_two_cat>::result Side_two_ide;
 
-  /*! Boolean tag that is bool_<true> if one side is identified,
-   * otherwise bool_<false>
+  /*! Boolean tag that is bool_constant<true> if one side is identified,
+   * otherwise bool_constant<false>
    */
-  typedef boost::mpl::or_<Side_one_ide, Side_two_ide>   result;
+  typedef std::bool_constant<Side_one_ide::value || Side_two_ide::value> result;
 };
 
 /*! Checks whether one of two boundary sides are contracted
@@ -464,10 +462,11 @@ struct Arr_has_contracted_sides_two {
   typedef typename Arr_is_side_contracted<Side_one_cat>::result Side_one_con;
   typedef typename Arr_is_side_contracted<Side_two_cat>::result Side_two_con;
 
-  /*!\ Boolean tag that is bool_<true> if one side is identified,
-   * otherwise bool_<false>
+  /*!\ Boolean tag that is bool_constant<true> if one side is identified,
+   * otherwise bool_constant<false>
    */
-  typedef boost::mpl::or_<Side_one_con, Side_two_con>           result;
+  typedef std::bool_constant<Side_one_con::value ||
+                             Side_two_con::value>               result;
 };
 
 /*! Checks whether one of two boundary sides are closed
@@ -480,10 +479,11 @@ struct Arr_has_closed_sides_two {
   typedef typename Arr_is_side_closed<Side_one_cat>::result     Side_one_clo;
   typedef typename Arr_is_side_closed<Side_two_cat>::result     Side_two_clo;
 
-  /*! Boolean tag that is bool_<true> if one side is identified,
-   * otherwise bool_<false>
+  /*! Boolean tag that is bool_constant<true> if one side is identified,
+   * otherwise bool_constant<false>
    */
-  typedef boost::mpl::or_<Side_one_clo, Side_two_clo>           result;
+  typedef std::bool_constant<Side_one_clo::value ||
+                             Side_two_clo::value>               result;
 };
 
 /*! Checks whether one of two boundary sides are open
@@ -496,10 +496,11 @@ struct Arr_has_open_sides_two {
   typedef typename Arr_is_side_open<Side_one_cat>::result       Side_one_ope;
   typedef typename Arr_is_side_open<Side_two_cat>::result       Side_two_ope;
 
-  /*! Boolean tag that is bool_<true> if one side is identified,
-   * otherwise bool_<false>
+  /*! Boolean tag that is bool_constant<true> if one side is identified,
+   * otherwise bool_constant<false>
    */
-  typedef boost::mpl::or_<Side_one_ope, Side_two_ope>           result;
+  typedef std::bool_constant<Side_one_ope::value ||
+                             Side_two_ope::value>               result;
 };
 
 /*! Categorizes two boundary sides:
@@ -532,11 +533,11 @@ struct Arr_two_sides_category {
     Is_open;
 
 public:
-  typedef typename boost::mpl::if_<Is_identified, Arr_has_identified_side_tag,
-    typename boost::mpl::if_<Is_contracted, Arr_has_contracted_side_tag,
-      typename boost::mpl::if_<Is_closed, Arr_has_closed_side_tag,
-        typename boost::mpl::if_<Is_open, Arr_has_open_side_tag,
-          Arr_all_sides_oblivious_tag>::type>::type>::type>::type
+  typedef std::conditional_t<Is_identified::value, Arr_has_identified_side_tag,
+    std::conditional_t<Is_contracted::value, Arr_has_contracted_side_tag,
+      std::conditional_t<Is_closed::value, Arr_has_closed_side_tag,
+        std::conditional_t<Is_open::value, Arr_has_open_side_tag,
+          Arr_all_sides_oblivious_tag>>>>
     result;
 };
 

Arrangement_on_surface_2/include/CGAL/Arr_vertical_decomposition_2.h

@@ -24,7 +24,6 @@
 #include <CGAL/Surface_sweep_2/Arr_vert_decomp_ss_visitor.h>
 
 #include <vector>
-#include <boost/mpl/if.hpp>
 #include <boost/type_traits.hpp>
 
 namespace CGAL {
@@ -124,7 +123,7 @@ decompose(const Arrangement_on_surface_2<GeometryTraits_2, TopologyTraits>& arr,
    * Use the form 'A a(*b);' and not ''A a = b;' to handle the case where A has
    * only an implicit constructor, (which takes *b as a parameter).
    */
-  typename boost::mpl::if_<std::is_same<Gt2, Vgt2>, const Vgt2&, Vgt2>::type
+  std::conditional_t<std::is_same_v<Gt2, Vgt2>, const Vgt2&, Vgt2>
     ex_traits(*geom_traits);
 
   // Define the sweep-line visitor and perform the sweep.

Arrangement_on_surface_2/include/CGAL/Arrangement_2/Arr_traits_adaptor_2_dispatching.h

@@ -17,7 +17,6 @@
 
 
 #include <boost/mpl/bool.hpp>
-#include <boost/mpl/if.hpp>
 #include <boost/mpl/and.hpp>
 #include <boost/mpl/or.hpp>
 #include <boost/mpl/not.hpp>
@@ -59,25 +58,22 @@ public:
 
 private:
 
-  typedef boost::mpl::bool_< true > true_;
-  typedef boost::mpl::bool_< false > false_;
+  typedef std::conditional_t<
+       std::is_same_v< Arr_smaller_implementation_tag, Arr_use_traits_tag >,
+       std::true_type, std::false_type > Smaller_traits;
 
-  typedef boost::mpl::if_<
-       std::is_same< Arr_smaller_implementation_tag, Arr_use_traits_tag >,
-       true_, false_ > Smaller_traits;
-
-  typedef boost::mpl::if_<
-       std::is_same< Arr_larger_implementation_tag, Arr_use_traits_tag >,
-       true_, false_ > Larger_traits;
+  typedef std::conditional_t<
+       std::is_same_v< Arr_larger_implementation_tag, Arr_use_traits_tag >,
+       std::true_type, std::false_type > Larger_traits;
 
 public:
 
   //! the result type (if one side asks for traits, then ask traits!
   //! Or vice versa: If both ask for dummy, then dummy!)
-  typedef typename boost::mpl::if_<
-              boost::mpl::or_< Smaller_traits, Larger_traits >,
+  typedef std::conditional_t<
+              Smaller_traits::value || Larger_traits::value,
               Arr_use_traits_tag,
-              Arr_use_dummy_tag >::type type;
+              Arr_use_dummy_tag > type;
 
 };
 

Arrangement_on_surface_2/include/CGAL/Arrangement_2/Arrangement_on_surface_2_global.h

@@ -23,8 +23,6 @@
 
 #include <list>
 #include <boost/type_traits.hpp>
-#include <boost/mpl/if.hpp>
-#include <boost/type_traits.hpp>
 #include <list>
 
 #include <CGAL/Arr_accessor.h>
@@ -66,7 +64,7 @@ namespace Ss2 = Surface_sweep_2;
 //
 // error: no matching function for call to `do_intersect(Arrangement_2<>&,
 // const Arr_segment_2&, const Arr_walk_along_line_point_location<>&,
-// mpl_::bool_< true>)'
+// std::bool_constant< true>)'
 //
 template <typename GeometryTraits_2, typename TopologyTraits,
           typename PointLocation, typename ZoneVisitor>
@@ -137,7 +135,7 @@ void insert(Arrangement_on_surface_2<GeometryTraits_2, TopologyTraits>& arr,
 //
 // error: no matching function for call to `do_intersect(Arrangement_2<>&,
 // const Arr_segment_2&, const Arr_walk_along_line_point_location<>&,
-// mpl_::bool_< true>)'
+// std::bool_constant< true>)'
 //
 //
 template <typename GeometryTraits_2, typename TopologyTraits,
@@ -271,7 +269,7 @@ insert_empty(Arrangement_on_surface_2<GeometryTraits_2, TopologyTraits>& arr,
    * Use the form 'A a(*b);' and not ''A a = b;' to handle the case where A has
    * only an implicit constructor, (which takes *b as a parameter).
    */
-  typename boost::mpl::if_<std::is_same<Gt2, Cgt2>, const Cgt2&, Cgt2>::type
+  std::conditional_t<std::is_same_v<Gt2, Cgt2>, const Cgt2&, Cgt2>
     traits(*geom_traits);
 
   // Define a surface-sweep instance and perform the sweep:
@@ -326,7 +324,7 @@ void insert_empty(Arrangement_on_surface_2<GeometryTraits_2, TopologyTraits>&
    * Use the form 'A a(*b);' and not ''A a = b;' to handle the case where A has
    * only an implicit constructor, (which takes *b as a parameter).
    */
-  typename boost::mpl::if_<std::is_same<Gt2, Cgt2>, const Cgt2&, Cgt2>::type
+  std::conditional_t<std::is_same_v<Gt2, Cgt2>, const Cgt2&, Cgt2>
     traits(*geom_traits);
 
   // Define a surface-sweep instance and perform the sweep.
@@ -379,7 +377,7 @@ void insert_non_empty(Arrangement_on_surface_2<GeometryTraits_2,
    * Use the form 'A a(*b);' and not ''A a = b;' to handle the case where A has
    * only an implicit constructor, (which takes *b as a parameter).
    */
-  typename boost::mpl::if_<std::is_same<Gt2, Igt2>, const Igt2&, Igt2>::type
+  std::conditional_t<std::is_same_v<Gt2, Igt2>, const Igt2&, Igt2>
     traits(*geom_traits);
 
   // Create a set of existing as well as new curves and points.
@@ -411,7 +409,7 @@ void insert_non_empty(Arrangement_on_surface_2<GeometryTraits_2,
 //
 // error: no matching function for call to `do_intersect(Arrangement_2<>&,
 // const Arr_segment_2&, const Arr_walk_along_line_point_location<>&,
-// mpl_::bool_< true>)'
+// std::bool_constant< true>)'
 //
 template <typename GeometryTraits_2, typename TopologyTraits,
           typename InputIterator>
@@ -465,7 +463,7 @@ void insert(Arrangement_on_surface_2<GeometryTraits_2, TopologyTraits>& arr,
 //
 // error: no matching function for call to `do_intersect(Arrangement_2<>&,
 // const Arr_segment_2&, const Arr_walk_along_line_point_location<>&,
-// mpl_::bool_< true>)'
+// std::bool_constant< true>)'
 //
 template <typename GeometryTraits_2, typename TopologyTraits,
           typename InputIterator>
@@ -979,7 +977,7 @@ non_intersecting_insert_non_empty(Arrangement_on_surface_2<GeometryTraits_2,
    * Use the form 'A a(*b);' and not ''A a = b;' to handle the case where A has
    * only an implicit constructor, (which takes *b as a parameter).
    */
-  typename boost::mpl::if_<std::is_same<Gt2, Igt2>, const Igt2&, Igt2>::type
+  std::conditional_t<std::is_same_v<Gt2, Igt2>, const Igt2&, Igt2>
     traits(*geom_traits);
 
   // Create a set of existing as well as new curves and points.
@@ -1526,7 +1524,7 @@ zone(Arrangement_on_surface_2<GeometryTraits_2, TopologyTraits>& arr,
 // workaround since it didn't compile in FC3_g++-3.4.4 with the error of:
 //
 // error: no matching function for call to `do_intersect(Arrangement_on_surface_2<>&,
-// const Arr_segment_2&, const Arr_walk_along_line_point_location<>&, mpl_::bool_< true>)'
+// const Arr_segment_2&, const Arr_walk_along_line_point_location<>&, std::bool_constant< true>)'
 //
 template <typename GeometryTraits_2, typename TopologyTraits,
           typename PointLocation>
@@ -1564,7 +1562,7 @@ do_intersect(Arrangement_on_surface_2<GeometryTraits_2, TopologyTraits>& arr,
 //
 // error: no matching function for call to
 // `do_intersect(Arrangement_on_surface_2<>&,
-// const Arr_segment_2&, const Arr_walk_along_line_point_location<>&, mpl_::bool_< true>)'
+// const Arr_segment_2&, const Arr_walk_along_line_point_location<>&, std::bool_constant< true>)'
 //
 template <typename GeometryTraits_2, typename TopologyTraits,
           typename PointLocation>

Arrangement_on_surface_2/include/CGAL/Surface_sweep_2/Arr_overlay_ss_visitor.h

@@ -295,7 +295,7 @@ protected:
    * \param tag The tag used for dispatching.
    */
   void _map_boundary_vertices(Event* event, Vertex_handle v,
-                              boost::mpl::bool_<true> /* tag */);
+                              std::bool_constant<true> /* tag */);
 
   /*!
    * Update the boundary vertices map.
@@ -306,7 +306,7 @@ protected:
    * \param tag The tag used for dispatching.
    */
   void _map_boundary_vertices(Event* event, Vertex_handle v,
-                              boost::mpl::bool_<false> /* tag */);
+                              std::bool_constant<false> /* tag */);
 
   /*!
    * Update a newly created vertex using the overlay traits.
@@ -319,7 +319,7 @@ protected:
    * \param tag The tag used for dispatching.
    */
   void _create_vertex(Event* event, Vertex_handle res_v, Subcurve* sc,
-                      boost::mpl::bool_<true> /* tag */);
+                      std::bool_constant<true> /* tag */);
 
   /*!
    * Update a newly created vertex using the overlay traits.
@@ -331,7 +331,7 @@ protected:
    * \param tag The tag used for dispatching.
    */
   void _create_vertex(Event* event, Vertex_handle res_v, Subcurve* sc,
-                      boost::mpl::bool_<false> /* tag */);
+                      std::bool_constant<false> /* tag */);
 
   /*!
    * Update a newly created edge using the overlay traits.
@@ -922,7 +922,7 @@ _map_halfedge_and_twin(Halfedge_handle he,
 //
 template <typename OvlHlpr, typename OvlTr, typename Vis>
 void Arr_overlay_ss_visitor<OvlHlpr, OvlTr, Vis>::
-_map_boundary_vertices(Event* event, Vertex_handle v, boost::mpl::bool_<true>)
+_map_boundary_vertices(Event* event, Vertex_handle v, std::bool_constant<true>)
 {
   // Update the red and blue object if the last event on sc is on the boundary.
   if ((event->parameter_space_in_x() != ARR_INTERIOR) ||
@@ -960,7 +960,7 @@ _map_boundary_vertices(Event* event, Vertex_handle v, boost::mpl::bool_<true>)
 template <typename OvlHlpr, typename OvlTr, typename Vis>
 void Arr_overlay_ss_visitor<OvlHlpr, OvlTr, Vis>::
 _map_boundary_vertices(Event* /* event */, Vertex_handle /* v */,
-                       boost::mpl::bool_<false>)
+                       std::bool_constant<false>)
 {}
 
 /* Notify the overlay traits about a newly created vertex.
@@ -974,7 +974,7 @@ void Arr_overlay_ss_visitor<OvlHlpr, OvlTr, Vis>::
 _create_vertex(Event* event,
                Vertex_handle new_v,
                Subcurve* sc,
-               boost::mpl::bool_<true>)
+               std::bool_constant<true>)
 {
   const Point_2& pt = event->point();
   const Cell_handle_red* red_handle = pt.red_cell_handle();
@@ -1011,7 +1011,7 @@ _create_vertex(Event* event,
     return;
   }
 
-  _create_vertex(event, new_v, sc, boost::mpl::bool_<false>());
+  _create_vertex(event, new_v, sc, std::bool_constant<false>());
 }
 
 /* Notify the overlay traits about a newly created vertex. */
@@ -1020,7 +1020,7 @@ void Arr_overlay_ss_visitor<OvlHlpr, OvlTr, Vis>::
 _create_vertex(Event* event,
                Vertex_handle new_v,
                Subcurve* sc,
-               boost::mpl::bool_<false>)
+               std::bool_constant<false>)
 {
   const Point_2& pt = event->point();
   const Cell_handle_red* red_handle = pt.red_cell_handle();

Arrangement_on_surface_2/test/Arrangement_on_surface_2/test_tags.cpp

@@ -6,7 +6,6 @@
 
 #include <boost/mpl/assert.hpp>
 #include <boost/mpl/bool.hpp>
-#include <boost/mpl/if.hpp>
 
 struct Traits1 {
   typedef CGAL::Arr_open_side_tag Left_side_category;

Arrangement_on_surface_2/test/Arrangement_on_surface_2/test_traits_dispatching.cpp

@@ -4,7 +4,6 @@
 #include <CGAL/Arrangement_2/Arr_traits_adaptor_2_dispatching.h>
 
 #include <boost/mpl/bool.hpp>
-#include <boost/mpl/if.hpp>
 
 int dispatch(CGAL::Arr_use_dummy_tag) {
   return 0;

Arrangement_on_surface_2/test/Arrangement_on_surface_2/utils.h

@@ -55,13 +55,13 @@ private:
   // This function is invoked for traits classes where at least one
   // boundary is not oblivious and all boundaries are not identified.
   bool operator()(const Point_2& p1, const Point_2& p2,
-                  boost::mpl::bool_<false>) const
+                  std::bool_constant<false>) const
   { return (m_traits.compare_xy_2_object()(p1, p2) == CGAL::SMALLER); }
 
   // This function should be invoked for traits classes where at least one
   // boundary is identified.
   bool operator()(const Point_2& p1, const Point_2& p2,
-                  boost::mpl::bool_<true>) const
+                  std::bool_constant<true>) const
   {
     // Compare in y boundaries:
     CGAL::Arr_parameter_space ps_y1 =

BGL/doc/BGL/Doxyfile.in

@@ -1,7 +1,8 @@
 @INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
 
 PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - CGAL and the Boost Graph Library"
-INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/Euler_operations.h \
+INPUT += ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/IO/polygon_mesh_io.h \
+         ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/Euler_operations.h \
          ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/iterator.h \
          ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/helpers.h \
          ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/boost/graph/generators.h \

BGL/examples/BGL_LCC/distance_lcc.cpp

@@ -1,7 +1,7 @@
 #include <CGAL/Simple_cartesian.h>
 
 #include <CGAL/boost/graph/graph_traits_Linear_cell_complex_for_combinatorial_map.h>
-#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 
 #include <CGAL/boost/graph/breadth_first_search.h> // wrapper to suppress a warning
 

BGL/examples/BGL_LCC/incident_vertices_lcc.cpp

@@ -2,7 +2,7 @@
 
 #include <CGAL/boost/graph/graph_traits_Linear_cell_complex_for_combinatorial_map.h>
 #include <CGAL/boost/graph/iterator.h>
-#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 
 #include <iostream>
 #include <fstream>

BGL/examples/BGL_LCC/normals_lcc.cpp

@@ -1,7 +1,7 @@
 #include <CGAL/Simple_cartesian.h>
 
 #include <CGAL/boost/graph/graph_traits_Linear_cell_complex_for_combinatorial_map.h>
-#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 #include <CGAL/property_map.h>
 
 #include <boost/graph/graph_traits.hpp>

BGL/examples/BGL_LCC/range_lcc.cpp

@@ -2,7 +2,7 @@
 
 #include <CGAL/Iterator_range.h>
 #include <CGAL/boost/graph/graph_traits_Linear_cell_complex_for_combinatorial_map.h>
-#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 
 #include <iostream>
 #include <fstream>

BGL/examples/BGL_LCC/transform_iterator_lcc.cpp

@@ -3,7 +3,7 @@
 #include <CGAL/boost/graph/graph_traits_Linear_cell_complex_for_combinatorial_map.h>
 #include <CGAL/boost/graph/iterator.h>
 #include <CGAL/boost/iterator/transform_iterator.hpp>
-#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 
 #include <algorithm>
 #include <fstream>

BGL/examples/BGL_OpenMesh/TriMesh.cpp

@@ -6,7 +6,7 @@
 #include <CGAL/boost/graph/graph_traits_TriMesh_ArrayKernelT.h>
 #include <CGAL/boost/graph/iterator.h>
 #include <CGAL/boost/graph/Euler_operations.h>
-#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 #include <CGAL/mesh_segmentation.h>
 #include <CGAL/property_map.h>
 

BGL/examples/BGL_surface_mesh/surface_mesh_partition.cpp

@@ -3,7 +3,7 @@
 
 #include <CGAL/boost/graph/Face_filtered_graph.h>
 #include <CGAL/boost/graph/partition.h>
-#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+#include <CGAL/IO/polygon_mesh_io.h>
 
 #include <fstream>
 #include <iostream>

BGL/include/CGAL/IO/polygon_mesh_io.h

@@ -0,0 +1,264 @@
+// Copyright (c) 2020  GeometryFactory (France).  All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Maxime Gimeno
+//                 Mael Rouxel-Labbé
+
+#ifndef CGAL_IO_POLYGON_MESH_IO_H
+#define CGAL_IO_POLYGON_MESH_IO_H
+
+#include <CGAL/boost/graph/IO/3MF.h>
+#include <CGAL/boost/graph/IO/GOCAD.h>
+#include <CGAL/boost/graph/IO/INP.h>
+#include <CGAL/boost/graph/IO/OBJ.h>
+#include <CGAL/boost/graph/IO/OFF.h>
+#include <CGAL/boost/graph/IO/PLY.h>
+#include <CGAL/boost/graph/IO/STL.h>
+#include <CGAL/boost/graph/IO/VTK.h>
+#include <CGAL/boost/graph/IO/WRL.h>
+#include <CGAL/IO/helpers.h>
+
+#include <fstream>
+#include <string>
+
+namespace CGAL {
+
+namespace IO {
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// Read
+
+//not for now : some readers will return "ok" despite not managing to read anything
+/*
+template <class Graph, typename NamedParameters = parameters::Default_named_parameters>
+bool read_polygon_mesh(std::istream& is,
+                       Graph& g,
+                       const NamedParameters& np = parameters::default_values())
+{
+  bool ok = false;
+  ok = read_OFF(is, g, np, false);
+  if(ok)
+    return true;
+  g.clear();
+  is.clear();//reset the error state
+  is.seekg (0, is.beg);
+  ok = read_OBJ(is, g, np, false);
+  if(ok)
+    return true;
+  g.clear();
+  is.clear();
+  is.seekg (0, is.beg);
+  ok = read_PLY(is, g, np, false);
+  if(ok)
+    return true;
+  g.clear();
+  is.clear();
+  is.seekg (0, is.beg);
+  ok = read_STL(is, g, np, false);
+  if(ok)
+    return true;
+  g.clear();
+  is.clear();
+  is.seekg (0, is.beg);
+  ok = read_GOCAD(is, g, np, false);
+  return ok;
+}
+
+*/
+
+/*!
+ * \ingroup PkgBGLIOFct
+ *
+ * \brief reads a polygon mesh from a file.
+ *
+ * Supported file formats are the following:
+ * - \ref IOStreamOFF (`.off`)
+ * - \ref IOStreamOBJ (`.obj`)
+ * - \ref IOStreamSTL (`.stl`)
+ * - \ref IOStreamPLY (`.ply`)
+ * - \ref IOStreamGocad (`.ts`)
+ * - \ref IOStreamVTK (`.vtp`)
+ *
+ * The format is detected from the filename extension (letter case is not important).
+ *
+ * The data is expected to represent a 2-manifold (possibly with borders).
+ *
+ * \tparam Graph a model of `MutableFaceGraph`
+ * \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+ *
+ * \param fname the name of the file
+ * \param g the mesh
+ * \param np optional \ref bgl_namedparameters "Named Parameters" described below
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{vertex_point_map}
+ *     \cgalParamDescription{a property map associating points to the vertices of `g`}
+ *     \cgalParamType{a class model of `WritablePropertyMap` with `boost::graph_traits<Graph>::%vertex_descriptor`
+ *                    as key type and `%Point_3` as value type}
+ *     \cgalParamDefault{`boost::get(CGAL::vertex_point, g)`}
+ *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+ *                     must be available in `Graph`.}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{verbose}
+ *     \cgalParamDescription{whether extra information is printed when an incident occurs during reading}
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`false`}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ * Other named parameters may be used according to the file extension, see \ref PkgBGLIOFct for an exhaustive list.
+ *
+ * \return `true` if reading was successful, `false` otherwise.
+ *
+ * \sa \link PMP_IO_grp `CGAL::Polygon_mesh_processing::IO::read_polygon_mesh()`\endlink if the data is not 2-manifold
+*/
+template <class Graph, typename NamedParameters = parameters::Default_named_parameters>
+bool read_polygon_mesh(const std::string& fname,
+                       Graph& g,
+                       const NamedParameters& np = parameters::default_values())
+{
+  const bool verbose = parameters::choose_parameter(parameters::get_parameter(np, internal_np::verbose), false);
+
+  const std::string ext = internal::get_file_extension(fname);
+  if(ext == std::string())
+  {
+    if(verbose)
+      std::cerr << "Error: cannot read from file without extension" << std::endl;
+    return false;
+  }
+
+  if(ext == "obj")
+    return read_OBJ(fname, g, np);
+  else if(ext == "off")
+    return read_OFF(fname, g, np);
+  else if(ext == "ply")
+    return read_PLY(fname, g, np);
+  else if(ext == "stl")
+    return read_STL(fname, g, np);
+  else if(ext == "ts")
+    return read_GOCAD(fname, g, np);
+#ifdef CGAL_USE_VTK
+  else if(ext == "vtp")
+    return read_VTP(fname, g, np);
+#endif
+
+  if(verbose)
+  {
+    std::cerr << "Error: unknown input file extension: " << ext << "\n"
+              << "Please refer to the documentation for the list of supported file formats" << std::endl;
+  }
+
+  return false;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// Write
+
+/*!
+ * \ingroup PkgBGLIOFct
+ *
+ * \brief writes a polygon mesh in a file.
+ *
+ * Supported file formats are the following:
+ * - \ref IOStreamOFF (`.off`)
+ * - \ref IOStreamOBJ (`.obj`)
+ * - \ref IOStreamSTL (`.stl`)
+ * - \ref IOStreamPLY (`.ply`)
+ * - \ref IOStreamGocad (`.ts`)
+ * - \ref IOStreamVTK (`.vtp`)
+ *
+ * The format is detected from the filename extension (letter case is not important).
+ *
+ * \tparam Graph a model of `FaceListGraph` and `HalfedgeListGraph`
+ * \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+ *
+ * \param fname the name of the file
+ * \param g the mesh to be output
+ * \param np optional \ref bgl_namedparameters "Named Parameters" described below
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{vertex_point_map}
+ *     \cgalParamDescription{a property map associating points to the vertices of `g`}
+ *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<Graph>::%vertex_descriptor`
+ *                    as key type and `%Point_3` as value type}
+ *     \cgalParamDefault{`boost::get(CGAL::vertex_point, g)`}
+ *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+ *                     must be available in `Graph`.}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{stream_precision}
+ *     \cgalParamDescription{a parameter used to set the precision (i.e. how many digits are generated) of the output stream}
+ *     \cgalParamType{int}
+ *     \cgalParamDefault{`6`}
+ *     \cgalParamExtra{This parameter is only meaningful while using \ascii encoding.}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{use_binary_mode}
+ *     \cgalParamDescription{indicates whether data should be written in binary (`true`) or in \ascii (`false`)}
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`true`}
+ *     \cgalParamExtra{This parameter is only meaningful for formats that support binary encoding.}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{verbose}
+ *     \cgalParamDescription{whether extra information is printed when an incident occurs during reading}
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`false`}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ * Other named parameters may be used according to the file extension, see \ref PkgBGLIOFct for an exhaustive list.
+ *
+ * \return `true` if writing was successful, `false` otherwise.
+ */
+template <class Graph, typename NamedParameters = parameters::Default_named_parameters>
+bool write_polygon_mesh(const std::string& fname,
+                        Graph& g,
+                        const NamedParameters& np = parameters::default_values())
+{
+  const bool verbose = parameters::choose_parameter(parameters::get_parameter(np, internal_np::verbose), false);
+
+  const std::string ext = internal::get_file_extension(fname);
+  if(ext == std::string())
+  {
+    if(verbose)
+      std::cerr << "Error: trying to output to file without extension" << std::endl;
+    return false;
+  }
+
+  if(ext == "obj")
+    return write_OBJ(fname, g, np);
+  else if(ext == "off")
+    return write_OFF(fname, g, np);
+  else if(ext == "ply")
+    return write_PLY(fname, g, np);
+  else if(ext == "stl")
+    return write_STL(fname, g, np);
+  else if(ext == "ts")
+    return write_GOCAD(fname, g, np);
+#ifdef CGAL_USE_VTK
+  else if(ext == "vtp")
+    return write_VTP(fname, g, np);
+#endif
+
+  if(verbose)
+  {
+    std::cerr << "Error: unknown output file extension: " << ext << "\n"
+              << "Please refer to the documentation for the list of supported file formats" << std::endl;
+  }
+
+  return false;
+}
+
+} // namespace IO
+} // namespace CGAL
+
+#endif // CGAL_IO_POLYGON_MESH_IO_H

BGL/include/CGAL/boost/graph/IO/polygon_mesh_io.h

@@ -1,4 +1,4 @@
-// Copyright (c) 2020  GeometryFactory (France).  All rights reserved.
+// Copyright (c) 2023  GeometryFactory (France).  All rights reserved.
 //
 // This file is part of CGAL (www.cgal.org)
 //
@@ -6,258 +6,11 @@
 // $Id$
 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
 //
-// Author(s)     : Maxime Gimeno
-//                 Mael Rouxel-Labbé
+// Author(s)     : Mael Rouxel-Labbé
 
 #ifndef CGAL_BOOST_GRAPH_POLYGON_MESH_IO_H
 #define CGAL_BOOST_GRAPH_POLYGON_MESH_IO_H
 
-#include <CGAL/boost/graph/IO/3MF.h>
-#include <CGAL/boost/graph/IO/GOCAD.h>
-#include <CGAL/boost/graph/IO/INP.h>
-#include <CGAL/boost/graph/IO/OBJ.h>
-#include <CGAL/boost/graph/IO/OFF.h>
-#include <CGAL/boost/graph/IO/PLY.h>
-#include <CGAL/boost/graph/IO/STL.h>
-#include <CGAL/boost/graph/IO/VTK.h>
-#include <CGAL/boost/graph/IO/WRL.h>
-#include <CGAL/IO/helpers.h>
-
-#include <fstream>
-#include <string>
-
-namespace CGAL {
-
-namespace IO {
-
-////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////
-// Read
-
-//not for now : some readers will return "ok" despite not managing to read anything
-/*
-template <class Graph, typename NamedParameters = parameters::Default_named_parameters>
-bool read_polygon_mesh(std::istream& is,
-                       Graph& g,
-                       const NamedParameters& np = parameters::default_values())
-{
-  bool ok = false;
-  ok = read_OFF(is, g, np, false);
-  if(ok)
-    return true;
-  g.clear();
-  is.clear();//reset the error state
-  is.seekg (0, is.beg);
-  ok = read_OBJ(is, g, np, false);
-  if(ok)
-    return true;
-  g.clear();
-  is.clear();
-  is.seekg (0, is.beg);
-  ok = read_PLY(is, g, np, false);
-  if(ok)
-    return true;
-  g.clear();
-  is.clear();
-  is.seekg (0, is.beg);
-  ok = read_STL(is, g, np, false);
-  if(ok)
-    return true;
-  g.clear();
-  is.clear();
-  is.seekg (0, is.beg);
-  ok = read_GOCAD(is, g, np, false);
-  return ok;
-}
-
-*/
-
-/*!
- * \ingroup PkgBGLIOFct
- *
- * \brief reads a polygon mesh from a file.
- *
- * Supported file formats are the following:
- * - \ref IOStreamOFF (`.off`)
- * - \ref IOStreamOBJ (`.obj`)
- * - \ref IOStreamSTL (`.stl`)
- * - \ref IOStreamPLY (`.ply`)
- * - \ref IOStreamGocad (`.ts`)
- * - \ref IOStreamVTK (`.vtp`)
- *
- * The format is detected from the filename extension (letter case is not important).
- *
- * The data is expected to represent a 2-manifold (possibly with borders).
- *
- * \tparam Graph a model of `MutableFaceGraph`
- * \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
- *
- * \param fname the name of the file
- * \param g the mesh
- * \param np optional \ref bgl_namedparameters "Named Parameters" described below
- *
- * \cgalNamedParamsBegin
- *   \cgalParamNBegin{vertex_point_map}
- *     \cgalParamDescription{a property map associating points to the vertices of `g`}
- *     \cgalParamType{a class model of `WritablePropertyMap` with `boost::graph_traits<Graph>::%vertex_descriptor`
- *                    as key type and `%Point_3` as value type}
- *     \cgalParamDefault{`boost::get(CGAL::vertex_point, g)`}
- *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
- *                     must be available in `Graph`.}
- *   \cgalParamNEnd
- *
- *   \cgalParamNBegin{verbose}
- *     \cgalParamDescription{whether extra information is printed when an incident occurs during reading}
- *     \cgalParamType{Boolean}
- *     \cgalParamDefault{`false`}
- *   \cgalParamNEnd
- * \cgalNamedParamsEnd
- *
- * Other named parameters may be used according to the file extension, see \ref PkgBGLIOFct for an exhaustive list.
- *
- * \return `true` if reading was successful, `false` otherwise.
- *
- * \sa \link PMP_IO_grp `CGAL::Polygon_mesh_processing::IO::read_polygon_mesh()`\endlink if the data is not 2-manifold
-*/
-template <class Graph, typename NamedParameters = parameters::Default_named_parameters>
-bool read_polygon_mesh(const std::string& fname,
-                       Graph& g,
-                       const NamedParameters& np = parameters::default_values())
-{
-  const bool verbose = parameters::choose_parameter(parameters::get_parameter(np, internal_np::verbose), false);
-
-  const std::string ext = internal::get_file_extension(fname);
-  if(ext == std::string())
-  {
-    if(verbose)
-      std::cerr << "Error: cannot read from file without extension" << std::endl;
-    return false;
-  }
-
-  if(ext == "obj")
-    return read_OBJ(fname, g, np);
-  else if(ext == "off")
-    return read_OFF(fname, g, np);
-  else if(ext == "ply")
-    return read_PLY(fname, g, np);
-  else if(ext == "stl")
-    return read_STL(fname, g, np);
-  else if(ext == "ts")
-    return read_GOCAD(fname, g, np);
-#ifdef CGAL_USE_VTK
-  else if(ext == "vtp")
-    return read_VTP(fname, g, np);
-#endif
-
-  if(verbose)
-  {
-    std::cerr << "Error: unknown input file extension: " << ext << "\n"
-              << "Please refer to the documentation for the list of supported file formats" << std::endl;
-  }
-
-  return false;
-}
-
-////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////
-// Write
-
-/*!
- * \ingroup PkgBGLIOFct
- *
- * \brief writes a polygon mesh in a file.
- *
- * Supported file formats are the following:
- * - \ref IOStreamOFF (`.off`)
- * - \ref IOStreamOBJ (`.obj`)
- * - \ref IOStreamSTL (`.stl`)
- * - \ref IOStreamPLY (`.ply`)
- * - \ref IOStreamGocad (`.ts`)
- * - \ref IOStreamVTK (`.vtp`)
- *
- * The format is detected from the filename extension (letter case is not important).
- *
- * \tparam Graph a model of `FaceListGraph` and `HalfedgeListGraph`
- * \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
- *
- * \param fname the name of the file
- * \param g the mesh to be output
- * \param np optional \ref bgl_namedparameters "Named Parameters" described below
- *
- * \cgalNamedParamsBegin
- *   \cgalParamNBegin{vertex_point_map}
- *     \cgalParamDescription{a property map associating points to the vertices of `g`}
- *     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<Graph>::%vertex_descriptor`
- *                    as key type and `%Point_3` as value type}
- *     \cgalParamDefault{`boost::get(CGAL::vertex_point, g)`}
- *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
- *                     must be available in `Graph`.}
- *   \cgalParamNEnd
- *
- *   \cgalParamNBegin{stream_precision}
- *     \cgalParamDescription{a parameter used to set the precision (i.e. how many digits are generated) of the output stream}
- *     \cgalParamType{int}
- *     \cgalParamDefault{`6`}
- *     \cgalParamExtra{This parameter is only meaningful while using \ascii encoding.}
- *   \cgalParamNEnd
- *
- *   \cgalParamNBegin{use_binary_mode}
- *     \cgalParamDescription{indicates whether data should be written in binary (`true`) or in \ascii (`false`)}
- *     \cgalParamType{Boolean}
- *     \cgalParamDefault{`true`}
- *     \cgalParamExtra{This parameter is only meaningful for formats that support binary encoding.}
- *   \cgalParamNEnd
- *
- *   \cgalParamNBegin{verbose}
- *     \cgalParamDescription{whether extra information is printed when an incident occurs during reading}
- *     \cgalParamType{Boolean}
- *     \cgalParamDefault{`false`}
- *   \cgalParamNEnd
- * \cgalNamedParamsEnd
- *
- * Other named parameters may be used according to the file extension, see \ref PkgBGLIOFct for an exhaustive list.
- *
- * \return `true` if writing was successful, `false` otherwise.
- */
-template <class Graph, typename NamedParameters = parameters::Default_named_parameters>
-bool write_polygon_mesh(const std::string& fname,
-                        Graph& g,
-                        const NamedParameters& np = parameters::default_values())
-{
-  const bool verbose = parameters::choose_parameter(parameters::get_parameter(np, internal_np::verbose), false);
-
-  const std::string ext = internal::get_file_extension(fname);
-  if(ext == std::string())
-  {
-    if(verbose)
-      std::cerr << "Error: trying to output to file without extension" << std::endl;
-    return false;
-  }
-
-  if(ext == "obj")
-    return write_OBJ(fname, g, np);
-  else if(ext == "off")
-    return write_OFF(fname, g, np);
-  else if(ext == "ply")
-    return write_PLY(fname, g, np);
-  else if(ext == "stl")
-    return write_STL(fname, g, np);
-  else if(ext == "ts")
-    return write_GOCAD(fname, g, np);
-#ifdef CGAL_USE_VTK
-  else if(ext == "vtp")
-    return write_VTP(fname, g, np);
-#endif
-
-  if(verbose)
-  {
-    std::cerr << "Error: unknown output file extension: " << ext << "\n"
-              << "Please refer to the documentation for the list of supported file formats" << std::endl;
-  }
-
-  return false;
-}
-
-}} // namespace CGAL::IO
+#include <CGAL/IO/polygon_mesh_io.h>
 
 #endif // CGAL_BOOST_GRAPH_POLYGON_MESH_IO_H

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

@@ -953,6 +953,11 @@ void swap_edges(const typename boost::graph_traits<FaceGraph>::halfedge_descript
   if (fo2 != nf && halfedge(fo2, g)==oh2) set_halfedge(fo2, oh1, g);
 }
 
+template <typename Graph>
+void collect_garbage(Graph&)
+{
+  // nothing by default
+}
 
 } //end of internal namespace
 

BGL/include/CGAL/boost/graph/internal/OM_iterator_from_circulator.h

@@ -20,8 +20,6 @@
 #include <boost/operators.hpp>
 #include <boost/concept_check.hpp>
 
-#include <boost/mpl/if.hpp>
-
 namespace CGAL {
 
 // adapted from circulator.h, does not support
@@ -45,23 +43,20 @@ public:
 
   typedef typename I__traits::iterator_category iterator_category;
 
-  typedef typename
-  boost::mpl::if_c<  Prevent_deref
-                   , C
-                   , typename C::value_type
-                  >::type             value_type;
+  typedef std::conditional_t<  Prevent_deref
+                             , C
+                             , typename C::value_type>
+                                                         value_type;
 
   typedef typename C::difference_type difference_type;
-  typedef typename
-  boost::mpl::if_c<  Prevent_deref
-                   , C&
-                   , typename C::reference
-                  >::type             reference;
-  typedef typename
-  boost::mpl::if_c<  Prevent_deref
-                   , C*
-                   , typename C::reference
-                  >::type             pointer;
+  typedef std::conditional_t<  Prevent_deref
+                             , C&
+                             , typename C::reference
+                            >             reference;
+  typedef std::conditional_t<  Prevent_deref
+                             , C*
+                             , typename C::reference
+                            >             pointer;
 
   OM_iterator_from_circulator(){}
 

BGL/include/CGAL/boost/graph/internal/initialized_index_maps_helpers.h

@@ -276,8 +276,8 @@ class GetInitializedIndexMap
 {
 public:
   // Check if there is an internal property map; if not, we must a dynamic property map
-  typedef typename boost::mpl::if_c<
-      CGAL::graph_has_property<Graph, Tag>::value, Tag, DynamicTag>::type    Final_tag;
+  typedef std::conditional_t<
+      CGAL::graph_has_property<Graph, Tag>::value, Tag, DynamicTag>    Final_tag;
 
   typedef typename internal_np::Lookup_named_param_def<
       PropertyTag,

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

@@ -20,7 +20,6 @@
 #include <CGAL/Named_function_parameters.h>
 #include <CGAL/property_map.h>
 
-#include <boost/mpl/if.hpp>
 #include <boost/mpl/has_xxx.hpp>
 
 #include <fstream>
@@ -43,14 +42,13 @@ class property_map_selector
 {
 public:
   typedef typename graph_has_property<PolygonMesh, PropertyTag>::type Has_internal_pmap;
-  typedef typename boost::mpl::if_c<Has_internal_pmap::value,
-                                    typename boost::property_map<PolygonMesh, PropertyTag>::type,
-                                    typename boost::cgal_no_property::type
-                                    >::type type;
-  typedef typename boost::mpl::if_c<Has_internal_pmap::value,
-                                    typename boost::property_map<PolygonMesh, PropertyTag>::const_type,
-                                    typename boost::cgal_no_property::const_type
-                                    >::type const_type;
+  typedef std::conditional_t<Has_internal_pmap::value,
+                             typename boost::property_map<PolygonMesh, PropertyTag>::type,
+                             typename boost::cgal_no_property::type> type;
+  typedef std::conditional_t<Has_internal_pmap::value,
+                             typename boost::property_map<PolygonMesh, PropertyTag>::const_type,
+                             typename boost::cgal_no_property::const_type
+                             > const_type;
 
   type get_pmap(const PropertyTag& p, PolygonMesh& pmesh)
   {
@@ -209,10 +207,10 @@ struct GetGeomTraits_impl<PolygonMesh, internal_np::Param_not_found, NamedParame
 
   struct Fake_GT {}; // to be used if there is no internal vertex_point_map in PolygonMesh
 
-  typedef typename boost::mpl::if_c<Has_internal_pmap::value ||
-                                    !std::is_same<internal_np::Param_not_found, NP_vpm>::value,
-                                    typename GetK<PolygonMesh, NamedParametersVPM>::Kernel,
-                                    Fake_GT>::type type;
+  typedef std::conditional_t<Has_internal_pmap::value ||
+                             !std::is_same<internal_np::Param_not_found, NP_vpm>::value,
+                             typename GetK<PolygonMesh, NamedParametersVPM>::Kernel,
+                             Fake_GT> type;
 };
 
 template <typename PolygonMesh,

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

@@ -140,9 +140,9 @@ struct Point_accessor<Handle, ValueType, ConstReference, true>
   typedef ValueType                      value_type;
   typedef Handle                         key_type;
 
-  typedef typename boost::mpl::if_< std::is_reference<ConstReference>,
-                                    ValueType&,
-                                    ValueType >::type Reference;
+  typedef std::conditional_t< std::is_reference_v<ConstReference>,
+                              ValueType&,
+                              ValueType > Reference;
 
   Point_accessor() {}
   Point_accessor(Point_accessor<Handle, ValueType, Reference, false>) {}
@@ -172,9 +172,9 @@ struct Is_writable_property_map<PropertyMap, boost::read_write_property_map_tag>
 // property map must define.
 template <typename PropertyMap>
 struct Is_writable_property_map<PropertyMap, boost::lvalue_property_map_tag>
-  : boost::mpl::if_c<std::is_const<typename std::remove_reference<
+  : std::conditional_t<std::is_const<typename std::remove_reference<
                        typename boost::property_traits<PropertyMap>::reference>::type>::value,
-                     CGAL::Tag_false, CGAL::Tag_true>::type
+                       CGAL::Tag_false, CGAL::Tag_true>
 { };
 
 } // namespace internal

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

@@ -13,7 +13,6 @@
 #include <CGAL/assertions.h>
 #include <CGAL/boost/graph/properties.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <boost/mpl/if.hpp>
 
 #ifndef OPEN_MESH_CLASS
   #error OPEN_MESH_CLASS is not defined
@@ -29,13 +28,13 @@ namespace CGAL {
 template <typename Mesh, typename Descriptor, typename Value>
 class OM_pmap {
 public:
-  typedef typename boost::mpl::if_<std::is_same<Descriptor, typename boost::graph_traits<Mesh>::vertex_descriptor>,
-                                   OpenMesh::VPropHandleT<Value>,
-                                   typename boost::mpl::if_<std::is_same<Descriptor, typename boost::graph_traits<Mesh>::face_descriptor>,
-                                                            OpenMesh::FPropHandleT<Value>,
-                                                            typename boost::mpl::if_<std::is_same<Descriptor, typename boost::graph_traits<Mesh>::halfedge_descriptor>,
-                                                                                     OpenMesh::HPropHandleT<Value>,
-                                                                                     OpenMesh::EPropHandleT<Value> >::type>::type>::type H;
+  typedef std::conditional_t<std::is_same_v<Descriptor, typename boost::graph_traits<Mesh>::vertex_descriptor>,
+                             OpenMesh::VPropHandleT<Value>,
+                             std::conditional_t<std::is_same_v<Descriptor, typename boost::graph_traits<Mesh>::face_descriptor>,
+                                                OpenMesh::FPropHandleT<Value>,
+                                                std::conditional_t<std::is_same_v<Descriptor, typename boost::graph_traits<Mesh>::halfedge_descriptor>,
+                                                                   OpenMesh::HPropHandleT<Value>,
+                                                                   OpenMesh::EPropHandleT<Value> >>> H;
 
   typedef boost::lvalue_property_map_tag category;
 

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

@@ -15,7 +15,6 @@
 
 #include <CGAL/property_map.h>
 #include <CGAL/boost/graph/properties.h>
-#include <boost/mpl/if.hpp>
 
 namespace CGAL{
 

CGAL_ipelets/include/CGAL/CGAL_Ipelet_base_v6.h

@@ -636,12 +636,11 @@ public:
     template<class iterator>
     void
     draw_in_ipe(const iterator begin,const iterator end,const Iso_rectangle_2& bbox,bool make_grp=true,bool deselect_all=false,
-     std::enable_if_t<  boost::mpl::or_< std::is_same<typename std::iterator_traits<iterator>::value_type,Point_2> ,
-                        boost::mpl::or_< std::is_same<typename std::iterator_traits<iterator>::value_type,Segment_2> ,
-                        boost::mpl::or_< std::is_same<typename std::iterator_traits<iterator>::value_type,Circle_2> ,
-                        boost::mpl::or_< std::is_same<typename std::iterator_traits<iterator>::value_type,Circular_arc_2> ,
-                                         std::is_same<typename std::iterator_traits<iterator>::value_type,Polygon_2>
-                                      > > > >::value
+     std::enable_if_t<  std::is_same_v<typename std::iterator_traits<iterator>::value_type,Point_2> ||
+                        std::is_same_v<typename std::iterator_traits<iterator>::value_type,Segment_2> ||
+                        std::is_same_v<typename std::iterator_traits<iterator>::value_type,Circle_2> ||
+                        std::is_same_v<typename std::iterator_traits<iterator>::value_type,Circular_arc_2> ||
+                        std::is_same_v<typename std::iterator_traits<iterator>::value_type,Polygon_2>
                     >* = nullptr) const
     {
       for (iterator it=begin;it!=end;++it)

CGAL_ipelets/include/CGAL/CGAL_Ipelet_base_v7.h

@@ -645,12 +645,11 @@ public:
     template<class iterator>
     void
     draw_in_ipe(const iterator begin,const iterator end,const Iso_rectangle_2& bbox,bool make_grp=true,bool deselect_all=false,
-     std::enable_if_t<  boost::mpl::or_< std::is_same<typename std::iterator_traits<iterator>::value_type,Point_2> ,
-                        boost::mpl::or_< std::is_same<typename std::iterator_traits<iterator>::value_type,Segment_2> ,
-                        boost::mpl::or_< std::is_same<typename std::iterator_traits<iterator>::value_type,Circle_2> ,
-                        boost::mpl::or_< std::is_same<typename std::iterator_traits<iterator>::value_type,Circular_arc_2> ,
-                                         std::is_same<typename std::iterator_traits<iterator>::value_type,Polygon_2>
-                                                > > > >::value
+     std::enable_if_t<  std::is_same_v<typename std::iterator_traits<iterator>::value_type,Point_2> ||
+                        std::is_same_v<typename std::iterator_traits<iterator>::value_type,Segment_2> ||
+                        std::is_same_v<typename std::iterator_traits<iterator>::value_type,Circle_2> ||
+                        std::is_same_v<typename std::iterator_traits<iterator>::value_type,Circular_arc_2> ||
+                        std::is_same_v<typename std::iterator_traits<iterator>::value_type,Polygon_2>
                     >* = nullptr) const
     {
       for (iterator it=begin;it!=end;++it)

Combinatorial_map/include/CGAL/Combinatorial_map_iterators_base.h

@@ -18,7 +18,6 @@
 #include <CGAL/Compact_container_with_index.h>
 
 #include <queue>
-#include <boost/mpl/if.hpp>
 
 #include <type_traits>
 
@@ -68,24 +67,24 @@ namespace CGAL {
   class CMap_dart_iterator;
 
   template < typename Map_,bool Const>
-  class CMap_dart_iterator<Map_, Const, Tag_false>: public boost::mpl::if_c< Const,
+  class CMap_dart_iterator<Map_, Const, Tag_false>: public std::conditional_t< Const,
       typename Map_::Dart_container::const_iterator,
-      typename Map_::Dart_container::iterator>::type
+      typename Map_::Dart_container::iterator>
     //public internal::CC_iterator<typename Map_::Dart_container,Const>
   {
   public:
     typedef CMap_dart_iterator<Map_,Const> Self;
 
-    typedef typename boost::mpl::if_c< Const,
+    typedef std::conditional_t< Const,
           typename Map_::Dart_container::const_iterator,
-          typename Map_::Dart_container::iterator>::type Base;
+          typename Map_::Dart_container::iterator> Base;
     // typedef internal::CC_iterator<typename Map_::Dart_container,Const> Base;
 
-    typedef typename boost::mpl::if_c< Const,
-                                       typename Map_::Dart_const_descriptor,
-                                       typename Map_::Dart_descriptor>::type
+    typedef std::conditional_t< Const,
+                                typename Map_::Dart_const_descriptor,
+                                typename Map_::Dart_descriptor>
                                      Dart_descriptor;
-    typedef typename boost::mpl::if_c< Const, const Map_, Map_>::type Map;
+    typedef std::conditional_t< Const, const Map_, Map_> Map;
 
     typedef std::input_iterator_tag iterator_category;
     typedef typename Base::value_type value_type;
@@ -180,25 +179,24 @@ namespace CGAL {
 
   template < typename Map_,bool Const >
   class CMap_dart_iterator<Map_, Const, Tag_true>:
-      /*public boost::mpl::if_c< Const,
+      /*public std::conditional_t< Const,
       typename Map_::Dart_container::const_iterator,
-      typename Map_::Dart_container::iterator>::type*/
+      typename Map_::Dart_container::iterator>*/
     public internal::CC_iterator_with_index<typename Map_::Dart_container,Const>
   {
   public:
     typedef CMap_dart_iterator<Map_,Const> Self;
 
-    /*typedef typename boost::mpl::if_c< Const,
+    /*typedef std::conditional_t< Const,
           typename Map_::Dart_container::const_iterator,
-          typename Map_::Dart_container::iterator>::type Base;*/
+          typename Map_::Dart_container::iterator> Base;*/
     typedef internal::CC_iterator_with_index<typename Map_::Dart_container,Const> Base;
 
-    typedef typename boost::mpl::if_c< Const,
-                                       typename Map_::Dart_const_descriptor,
-                                       typename Map_::Dart_descriptor>::type
+    typedef std::conditional_t< Const,
+                                typename Map_::Dart_const_descriptor,
+                                typename Map_::Dart_descriptor>
                                      Dart_descriptor;
-    typedef typename boost::mpl::if_c< Const, const Map_,
-                                       Map_>::type Map;
+    typedef std::conditional_t< Const, const Map_, Map_> Map;
 
     typedef std::input_iterator_tag iterator_category;
     typedef typename Base::value_type value_type;

Combinatorial_map/include/CGAL/Compact_container_with_index.h

@@ -861,14 +861,13 @@ namespace internal {
     typedef typename DSC::value_type                  value_type;
     typedef typename DSC::size_type                   size_type;
     typedef typename DSC::difference_type             difference_type;
-    typedef typename boost::mpl::if_c< Const, const value_type*,
-                                       value_type*>::type pointer;
-    typedef typename boost::mpl::if_c< Const, const value_type&,
-                                       value_type&>::type reference;
+    typedef std::conditional_t< Const, const value_type*,
+                                       value_type*>   pointer;
+    typedef std::conditional_t< Const, const value_type&,
+                                       value_type&>  reference;
     typedef std::bidirectional_iterator_tag           iterator_category;
 
-    typedef typename boost::mpl::if_c< Const, const DSC*, DSC*>::type
-    cc_pointer;
+    typedef std::conditional_t< Const, const DSC*, DSC*> cc_pointer;
 
     CC_iterator_with_index(): m_ptr_to_cc(nullptr),
       m_index(0)

Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h

@@ -26,7 +26,9 @@
 #include <CGAL/assertions.h>
 #include <CGAL/boost/graph/Euler_operations.h>
 // For interior_polyhedron_3
+#ifndef CGAL_CH3_DUAL_WITHOUT_QP_SOLVER
 #include <CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h>
+#endif
 #include <CGAL/Number_types/internal/Exact_type_selector.h>
 
 #include <unordered_map>
@@ -227,7 +229,11 @@ namespace CGAL
     template <class PlaneIterator, class Polyhedron>
     void halfspace_intersection_3 (PlaneIterator begin, PlaneIterator end,
                                    Polyhedron &P,
-                                   std::optional<typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3> origin = std::nullopt) {
+                                   std::optional<typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3> origin
+#ifndef CGAL_CH3_DUAL_WITHOUT_QP_SOLVER
+                                   = std::nullopt
+#endif
+    ) {
         // Checks whether the intersection is a polyhedron
         CGAL_assertion_msg(Convex_hull_3::internal::is_intersection_dim_3(begin, end), "halfspace_intersection_3: intersection not a polyhedron");
 
@@ -238,8 +244,10 @@ namespace CGAL
 
         // if a point inside is not provided find one using linear programming
         if (!origin) {
+#ifndef CGAL_CH3_DUAL_WITHOUT_QP_SOLVER
           // find a point inside the intersection
           origin = halfspace_intersection_interior_point_3(begin, end);
+#endif
 
           CGAL_assertion_msg(origin!=std::nullopt, "halfspace_intersection_3: problem when determining a point inside the intersection");
           if (origin==std::nullopt)

Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h

@@ -23,7 +23,9 @@
 #include <CGAL/assertions.h>
 
 // For interior_polyhedron_3
+#ifndef CGAL_CH3_DUAL_WITHOUT_QP_SOLVER
 #include <CGAL/Convex_hull_3/dual/halfspace_intersection_interior_point_3.h>
+#endif
 #include <CGAL/Number_types/internal/Exact_type_selector.h>
 
 #include <unordered_map>
@@ -99,7 +101,9 @@ namespace CGAL
           // if a point inside is not provided find one using linear programming
           if (!origin) {
             // find a point inside the intersection
+#ifndef CGAL_CH3_DUAL_WITHOUT_QP_SOLVER
             origin = halfspace_intersection_interior_point_3(pbegin, pend);
+#endif
 
             CGAL_assertion_msg(origin!=std::nullopt, "halfspace_intersection_with_constructions_3: problem when determining a point inside the intersection");
             if (origin==std::nullopt)
@@ -134,7 +138,11 @@ namespace CGAL
           void halfspace_intersection_with_constructions_3 (PlaneIterator pbegin,
                                                             PlaneIterator pend,
                                                             Polyhedron &P,
-                                                            std::optional<typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3> const& origin = std::nullopt) {
+                                                            std::optional<typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel::Point_3> const& origin
+#ifndef CGAL_CH3_DUAL_WITHOUT_QP_SOLVER
+                                                            = std::nullopt
+#endif
+        ) {
           typedef typename Kernel_traits<typename std::iterator_traits<PlaneIterator>::value_type>::Kernel K;
           typedef typename K::Point_3 Point_3;
           typedef typename Convex_hull_3::internal::Default_traits_for_Chull_3<Point_3>::type Traits;

Data/data/meshes/regular_tetrahedron.off

@@ -0,0 +1,11 @@
+OFF
+4 4 0
+-1 0 -0.707107
+1 0 -0.707107
+0 -1 0.707107
+0 1 0.707107
+3 0 1 2
+3 0 3 1
+3 0 2 3
+3 1 3 2
+

Documentation/doc/biblio/cgal_manual.bib

@@ -778,7 +778,7 @@ Teillaud"
 , volume      = 44
 , year        = 2011
 , pages       = "160--168"
-, url         = "https://theses.hal.science/inria-00560388/"
+, url         = "https://hal.science/inria-00560388/"
 , doi         = "10.1016/j.comgeo.2010.09.010"
 }
 

Filtered_kernel/include/CGAL/Lazy_kernel.h

@@ -27,7 +27,6 @@
 #include <CGAL/Filtered_kernel/internal/Static_filters/tools.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <boost/none.hpp>
-#include <boost/mpl/if.hpp>
 #include <boost/mpl/or.hpp>
 #include <CGAL/Lazy_exact_nt.h>
 
@@ -199,9 +198,9 @@ private:
                                               boost::mpl::eval_if< std::is_same< typename internal::Lazy_result_type<Construction>::type,
                                                                                    CGAL::Object >,
                                                                    boost::mpl::int_<OBJECT>,
-                                                                   boost::mpl::eval_if< boost::mpl::or_<
-                                                                                          std::is_same< typename internal::Lazy_result_type<Construction>::type, CGAL::Bbox_2 >,
-                                                                                          std::is_same< typename internal::Lazy_result_type<Construction>::type, CGAL::Bbox_3 > >,
+                                                                   boost::mpl::eval_if< std::bool_constant<
+                                                                                          std::is_same_v< typename internal::Lazy_result_type<Construction>::type, CGAL::Bbox_2 > ||
+                                                                                          std::is_same_v< typename internal::Lazy_result_type<Construction>::type, CGAL::Bbox_3 > >,
                                                                                         boost::mpl::int_<BBOX>,
                                                                                         boost::mpl::int_<NONE> > > >,
                          boost::mpl::int_<NONE> >::type {};

Generalized_map/include/CGAL/Generalized_map_iterators_base.h

@@ -19,7 +19,6 @@
 #include <CGAL/Compact_container.h>
 #include <queue>
 #include <type_traits>
-#include <boost/mpl/if.hpp>
 
 namespace CGAL {
 

Homogeneous_kernel/include/CGAL/Homogeneous/PointH2.h

@@ -53,8 +53,8 @@ public:
 
     template < typename Tx, typename Ty >
     PointH2(const Tx & x, const Ty & y,
-            std::enable_if_t< boost::mpl::and_<std::is_convertible<Tx, RT>,
-                                               std::is_convertible<Ty, RT> >::value >* = 0)
+            std::enable_if_t< std::is_convertible_v<Tx, RT> &&
+                              std::is_convertible_v<Ty, RT> >* = 0)
       : base(x, y) {}
 
     PointH2(const FT& x, const FT& y)

Homogeneous_kernel/include/CGAL/Homogeneous/PointH3.h

@@ -52,9 +52,9 @@ public:
 
   template < typename Tx, typename Ty, typename Tz >
   PointH3(const Tx & x, const Ty & y, const Tz & z,
-          std::enable_if_t< boost::mpl::and_< boost::mpl::and_< std::is_convertible<Tx, RT>,
-                                                                         std::is_convertible<Ty, RT> >,
-                                                       std::is_convertible<Tz, RT> >::value >* = 0)
+          std::enable_if_t<std::is_convertible_v<Tx, RT> &&
+                           std::is_convertible_v<Ty, RT> &&
+                           std::is_convertible_v<Tz, RT>>* = 0)
     : base(x, y, z) {}
 
   PointH3(const FT& x, const FT& y, const FT& z)

Homogeneous_kernel/include/CGAL/Homogeneous/VectorH2.h

@@ -57,8 +57,8 @@ public:
 
    template < typename Tx, typename Ty >
    VectorH2(const Tx & x, const Ty & y,
-            std::enable_if_t< boost::mpl::and_<std::is_convertible<Tx, RT>,
-                                                        std::is_convertible<Ty, RT> >::value >* = 0)
+            std::enable_if_t<std::is_convertible_v<Tx, RT> &&
+                             std::is_convertible_v<Ty, RT>>* = 0)
       : base(CGAL::make_array<RT>(x, y, RT(1))) {}
 
    VectorH2(const FT& x, const FT& y)

Homogeneous_kernel/include/CGAL/Homogeneous/VectorH3.h

@@ -67,9 +67,9 @@ public:
 
   template < typename Tx, typename Ty, typename Tz >
   VectorH3(const Tx & x, const Ty & y, const Tz & z,
-           std::enable_if_t< boost::mpl::and_< boost::mpl::and_< std::is_convertible<Tx, RT>,
-                                                                          std::is_convertible<Ty, RT> >,
-                                                        std::is_convertible<Tz, RT> >::value >* = 0)
+           std::enable_if_t< std::is_convertible_v<Tx, RT> &&
+                             std::is_convertible_v<Ty, RT> &&
+                             std::is_convertible_v<Tz, RT>>* = 0)
     : base(CGAL::make_array<RT>(x, y, z, RT(1))) {}
 
   VectorH3(const FT& x, const FT& y, const FT& z)

Installation/CHANGES.md

@@ -14,7 +14,7 @@ Release date: October 2023
 - **Breaking change**: The usage of `boost::shared_ptr` has been replaced by `std::shared_ptr`. Packages affected are 2D Straight Line Skeleton and Shape Detection.
 - **Breaking change**: The usage of `boost::optional` has been replaced by `std::optional`. Packages affected are 2D Straight Line Skeleton, 3D Fast Intersection and Distance Computation (AABB Tree), and the Kernel intersection.
 - **Breaking change**: The usage of `boost::variant` has been replaced by `std::variant`. Packages affected are 2D Arrangements, and the Kernel intersection.
-- **Breaking chahge**: The file CMake file `UseCGAL.cmake` has been removed from CGAL. Usages of the CMake variables `${CGAL_USE_FILE}` and `${CGAL_LIBRARIES}` must be replaced by a link to the imported target `CGAL::CGAL`, for example: `target_link_library(the_target PRIVATE CGAL::CGAL)`.
+- **Breaking change**: The file CMake file `UseCGAL.cmake` has been removed from CGAL. Usages of the CMake variables `${CGAL_USE_FILE}` and `${CGAL_LIBRARIES}` must be replaced by a link to the imported target `CGAL::CGAL`, for example: `target_link_library(the_target PRIVATE CGAL::CGAL)`.
 
 
 #### 2D Arrangements
@@ -44,6 +44,11 @@ Release date: October 2023
     the mean and Gaussian curvatures, as well as the principal curvature and directions.
 -   Added the option to use a variable sizing field for `CGAL::Polygon_mesh_processing::isotropic_remeshing()`,
     and a sizing function based on a measure of local curvature for adaptive remeshing.
+-   Added the function `CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel()` that refines a polygon mesh along an isocurve.
+-   Added the function
+    `CGAL::Polygon_mesh_processing::autorefine_triangle_soup()` that refines a soup of triangles so that no pair of triangles intersects
+     in their interiors. Also added, the function `autorefine()` operating directly on a triangle mesh and updating it
+     using the aforementioned function on a triangle soup.
 
 ### [2D Arrangements](https://doc.cgal.org/6.0/Manual/packages.html#PkgArrangementOnSurface2)
 -   Fixed a bug in the zone construction code applied to arrangements of geodesic arcs on a sphere,

Installation/CMakeLists.txt

@@ -1253,6 +1253,6 @@ if(RUNNING_CGAL_AUTO_TEST OR CGAL_TEST_SUITE)
     "USING BOOST_VERSION = '${Boost_MAJOR_VERSION}.${Boost_MINOR_VERSION}.${Boost_SUBMINOR_VERSION}'"
     )
   if(Qt6_FOUND)
-    message(STATUS "USING Qt6_VERSION = '${Qt6Core_VERSION_STRING}'")
+    message(STATUS "USING Qt6_VERSION = '${Qt6Core_VERSION}'")
   endif()#Qt6_FOUND
 endif()#RUNNING_CGAL_AUTO_TEST OR CGAL_TEST_SUITE

Installation/include/CGAL/license/Polygon_mesh_processing/autorefinement.h

@@ -0,0 +1,54 @@
+// Copyright (c) 2016  GeometryFactory SARL (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s) : Andreas Fabri
+//
+// Warning: this file is generated, see include/CGAL/license/README.md
+
+#ifndef CGAL_LICENSE_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H
+#define CGAL_LICENSE_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H
+
+#include <CGAL/config.h>
+#include <CGAL/license.h>
+
+#ifdef CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_COMMERCIAL_LICENSE
+
+#  if CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE
+
+#    if defined(CGAL_LICENSE_WARNING)
+
+       CGAL_pragma_warning("Your commercial license for CGAL does not cover "
+                           "this release of the Polygon Mesh Processing - Autorefinement package.")
+#    endif
+
+#    ifdef CGAL_LICENSE_ERROR
+#      error "Your commercial license for CGAL does not cover this release \
+              of the Polygon Mesh Processing - Autorefinement package. \
+              You get this error, as you defined CGAL_LICENSE_ERROR."
+#    endif // CGAL_LICENSE_ERROR
+
+#  endif // CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE
+
+#else // no CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_COMMERCIAL_LICENSE
+
+#  if defined(CGAL_LICENSE_WARNING)
+     CGAL_pragma_warning("\nThe macro CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_COMMERCIAL_LICENSE is not defined."
+                          "\nYou use the CGAL Polygon Mesh Processing - Autorefinement package under "
+                          "the terms of the GPLv3+.")
+#  endif // CGAL_LICENSE_WARNING
+
+#  ifdef CGAL_LICENSE_ERROR
+#    error "The macro CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_COMMERCIAL_LICENSE is not defined.\
+            You use the CGAL Polygon Mesh Processing - Autorefinement package under the terms of \
+            the GPLv3+. You get this error, as you defined CGAL_LICENSE_ERROR."
+#  endif // CGAL_LICENSE_ERROR
+
+#endif // no CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_COMMERCIAL_LICENSE
+
+#endif // CGAL_LICENSE_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H

Installation/include/CGAL/license/Polygon_mesh_processing/interpolated_corrected_curvatures.h

@@ -9,10 +9,10 @@
 //
 // Author(s) : Andreas Fabri
 //
-// Warning: this file is generated, see include/CGAL/licence/README.md
+// Warning: this file is generated, see include/CGAL/license/README.md
 
-#ifndef CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H
-#define CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H
+#ifndef CGAL_LICENSE_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H
+#define CGAL_LICENSE_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H
 
 #include <CGAL/config.h>
 #include <CGAL/license.h>
@@ -51,4 +51,4 @@
 
 #endif // no CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_COMMERCIAL_LICENSE
 
-#endif // CGAL_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H
+#endif // CGAL_LICENSE_POLYGON_MESH_PROCESSING_INTERPOLATED_CORRECTED_CURVATURES_H

Installation/include/CGAL/license/SMDS_3.h

@@ -24,12 +24,12 @@
 #    if defined(CGAL_LICENSE_WARNING)
 
        CGAL_pragma_warning("Your commercial license for CGAL does not cover "
-                           "this release of the 3D Mesh Data Structure package.")
+                           "this release of the 3D Simplicial Mesh Data Structure package.")
 #    endif
 
 #    ifdef CGAL_LICENSE_ERROR
 #      error "Your commercial license for CGAL does not cover this release \
-              of the 3D Mesh Data Structure package. \
+              of the 3D Simplicial Mesh Data Structure package. \
               You get this error, as you defined CGAL_LICENSE_ERROR."
 #    endif // CGAL_LICENSE_ERROR
 

Installation/include/CGAL/license/gpl_package_list.txt

@@ -61,6 +61,7 @@ Polygon_mesh_processing/geometric_repair Polygon Mesh Processing - Geometric Rep
 Polygon_mesh_processing/miscellaneous Polygon Mesh Processing - Miscellaneous
 Polygon_mesh_processing/detect_features Polygon Mesh Processing - Feature Detection
 Polygon_mesh_processing/collision_detection Polygon Mesh Processing - Collision Detection
+Polygon_mesh_processing/autorefinement Polygon Mesh Processing - Autorefinement
 Polyhedron 3D Polyhedral Surface
 Polyline_simplification_2 2D Polyline Simplification
 Polytope_distance_d Optimal Distances

Interpolation/include/CGAL/sibson_gradient_fitting.h

@@ -23,7 +23,6 @@
 #include <functional>
 
 #include <any>
-#include <boost/mpl/if.hpp>
 #include <boost/utility/result_of.hpp>
 
 #include <iterator>

Intersections_3/include/CGAL/Intersections_3/internal/Plane_3_Plane_3_Plane_3_intersection.h

@@ -84,10 +84,13 @@ intersection(const typename K::Plane_3& plane1,
                                                   typename K::Line_3,
                                                   typename K::Plane_3> > result_type;
 
-  typedef typename K::Point_3      Point_3;
   typedef typename K::Line_3       Line_3;
   typedef typename K::Plane_3      Plane_3;
 
+  auto res = intersection_point(plane1,plane2,plane3, k);
+  if (res)
+    return result_type(*res);
+
   // Intersection between plane1 and plane2 can either be
   // a line, a plane, or empty.
   typename Intersection_traits<K, Plane_3, Plane_3>::result_type
@@ -97,26 +100,19 @@ intersection(const typename K::Plane_3& plane1,
   {
     if(const Line_3* l = intersect_get<Line_3>(o12))
     {
-      // either point or line
-      typename Intersection_traits<K, Plane_3, Line_3>::result_type
-          v = internal::intersection(plane3, *l, k);
-      if(v)
-      {
-        if(const Point_3* p = intersect_get<Point_3>(v))
-          return result_type(*p);
-        else if(const Line_3* l = intersect_get<Line_3>(v))
-          return result_type(*l);
-      }
+      if (internal::do_intersect(*l, plane3, k))
+        return result_type(*l);
     }
-    else if(const Plane_3 *pl = intersect_get<Plane_3>(o12))
+    else
     {
+      CGAL_assertion(intersect_get<Plane_3>(o12) != nullptr);
       // either line or plane
       typename Intersection_traits<K, Plane_3, Plane_3>::result_type
-          v = internal::intersection(plane3, *pl, k);
+          v = internal::intersection(plane3, plane1, k);
       if(v)
       {
-        if(const Plane_3* p = intersect_get<Plane_3>(v))
-          return result_type(*p);
+        if( intersect_get<Plane_3>(v)!=nullptr)
+          return result_type(plane1);
         else if(const Line_3* l = intersect_get<Line_3>(v))
           return result_type(*l);
       }

Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h

@@ -14,6 +14,8 @@
 #ifndef CGAL_INTERNAL_INTERSECTIONS_TRIANGLE_3_TRIANGLE_3_INTERSECTION_H
 #define CGAL_INTERNAL_INTERSECTIONS_TRIANGLE_3_TRIANGLE_3_INTERSECTION_H
 
+//#define CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+
 #include <CGAL/Intersection_traits_3.h>
 #include <CGAL/Intersections_3/internal/Line_3_Triangle_3_intersection.h>
 #include <CGAL/Intersections_3/internal/Line_3_Line_3_intersection.h>
@@ -21,6 +23,9 @@
 
 #include <CGAL/kernel_assertions.h>
 
+#include <boost/next_prior.hpp>
+#include <boost/container/flat_set.hpp>
+
 #include <list>
 #include <map>
 #include <vector>
@@ -30,53 +35,338 @@ namespace Intersections {
 namespace internal{
 
 template <class Kernel>
-void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p,
-                                            const typename Kernel::Point_3& q,
-                                            const typename Kernel::Point_3& r,
-                                            const Kernel& k,
-                                            std::list<typename Kernel::Point_3>& inter_pts)
+struct Point_on_triangle
 {
-  typedef typename std::list<typename Kernel::Point_3>::iterator Iterator;
+  // triangle points are not stored in this class but are expected
+  // to always be passed in the same order. For a triangle pqr,
+  // edge 0 is pq, edge 1 qr and edge 2 rp. Point 0 is p, 1 is q and 2 is r.
+  //
+  // (id, -1) point on t1
+  // (-1, id) point on t2
+  // (id1, id2) intersection of edges
+  std::pair<int, int> t1_t2_ids;
+  boost::container::flat_set<int> extra_t1; // store other ids of edges containing the point
+  typename Kernel::FT alpha; //
+
+//////
+
+  static
+  inline
+  const typename Kernel::Point_3&
+  point_from_id(const typename Kernel::Point_3& p,
+                const typename Kernel::Point_3& q,
+                const typename Kernel::Point_3& r,
+                int id)
+  {
+    switch(id)
+    {
+      case 0:
+        return p;
+      case 1:
+        return q;
+      default:
+        return r;
+    }
+  }
+
+  Point_on_triangle(int i1, int i2=-1, typename Kernel::FT alpha = 0.) // TODO add global zero()?
+    : t1_t2_ids(i1,i2)
+    , alpha(alpha)
+  {}
+
+  // orientation of the current point wrt to edge id1 (p1q1)
+  Orientation
+  orientation (const typename Kernel::Point_3& p1, // source of edge edge_id1
+               const typename Kernel::Point_3& q1, // target of edge edge_id1
+               const typename Kernel::Point_3& r1,
+               int edge_id1,
+               const typename Kernel::Point_3& p2,
+               const typename Kernel::Point_3& q2,
+               const typename Kernel::Point_3& r2,
+               const Kernel& k) const
+  {
+    if (t1_t2_ids.first!=-1)
+    {
+      if (t1_t2_ids.second==-1)
+        return (edge_id1==t1_t2_ids.first || (edge_id1+1)%3==t1_t2_ids.first) ? ZERO:POSITIVE; // it is a point on t1
+      // this is an intersection point
+
+      if (t1_t2_ids.first==edge_id1)
+        return ZERO;
+      if (t1_t2_ids.first==(edge_id1+1)%3)
+      {
+        if (alpha==0) return ZERO;
+        return alpha>=0 ? POSITIVE:NEGATIVE;
+      }
+      CGAL_assertion((t1_t2_ids.first+1)%3==edge_id1);
+      if (alpha==1) return ZERO;
+      return alpha<=1?POSITIVE:NEGATIVE;
+    }
+    else
+    {
+      //this is an input point of t2
+      typename Kernel::Coplanar_orientation_3 orient = k.coplanar_orientation_3_object();
+      const typename Kernel::Point_3& query = point_from_id(p2,q2,r2,t1_t2_ids.second);
+      return orient(p1,q1,r1,query);
+    }
+  }
+
+  int id1() const { return t1_t2_ids.first; }
+  int id2() const { return t1_t2_ids.second; }
+
+  // construct the intersection point from the info stored
+  typename Kernel::Point_3
+  point(const typename Kernel::Point_3& p1,
+        const typename Kernel::Point_3& q1,
+        const typename Kernel::Point_3& r1,
+        const typename Kernel::Point_3& p2,
+        const typename Kernel::Point_3& q2,
+        const typename Kernel::Point_3& r2,
+        const Kernel& k) const
+  {
+    if (t1_t2_ids.first==-1)
+      return point_from_id(p2,q2,r2,t1_t2_ids.second);
+    if (t1_t2_ids.second==-1)
+      return point_from_id(p1,q1,r1,t1_t2_ids.first);
+
+    return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,(t1_t2_ids.first+1)%3), alpha, point_from_id(p1,q1,r1,t1_t2_ids.first)) ;
+  }
+};
+
+// the intersection of two triangles is computed by interatively intersection t2
+// with halfspaces defined by edges of t1. The following function is called
+// for each each on t1 on edge of the current intersection.
+// pq is such an edge and p1q1 from t1 defines the halfspace intersection
+// we are currently interseted in. We return the intersection point of
+// pq with p1q1
+template <class Kernel>
+Point_on_triangle<Kernel>
+intersection(const Point_on_triangle<Kernel>& p,
+             const Point_on_triangle<Kernel>& q,
+             int edge_id_t1,
+             const typename Kernel::Point_3& p1,
+             const typename Kernel::Point_3& q1,
+//             const typename Kernel::Point_3& r1,
+             const typename Kernel::Point_3& p2,
+             const typename Kernel::Point_3& q2,
+             const typename Kernel::Point_3& r2,
+             const Kernel& k)
+{
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "    calling intersection: ";
+  std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) -";
+  std::cout << " (" << q.id1() << "," << q.id2() << ",[" << q.alpha << "]) || e" << edge_id_t1;
+#endif
+  typename Kernel::Compute_alpha_for_coplanar_triangle_intersection_3 compute_alpha
+    = k.compute_alpha_for_coplanar_triangle_intersection_3_object();
+  typedef Point_on_triangle<Kernel> Pot;
+  switch(p.id1())
+  {
+    case -1:
+    {
+      switch(q.id1())
+      {
+        case -1: // A: (-1, ip2) - (-1, iq2)
+        {
+          CGAL_assertion((p.id2()+1)%3 == q.id2() || (q.id2()+1)%3 == p.id2());
+//          CGAL_assertion(p.extra_t1.empty() && q.extra_t1.empty());  // TMP to see if it's worth implementing special case
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+          std::cout << " -- case 1\n";
+#endif
+          typename Kernel::FT alpha = compute_alpha(p1, q1,
+                                                    Pot::point_from_id(p2, q2, r2, p.id2()),
+                                                    Pot::point_from_id(p2, q2, r2, q.id2()));
+          int id2 = (p.id2()+1)%3 == q.id2() ? p.id2() : q.id2();
+          return  Point_on_triangle<Kernel>(edge_id_t1, id2, alpha); // intersection with an original edge of t2
+        }
+        default:
+          if (q.id2()!=-1) // B: (-1, ip2) - (iq1, iq2)
+          {
+            if (p.id2() == q.id2() || p.id2() == (q.id2()+1)%3)
+            {
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+              std::cout << " -- case 2\n";
+#endif
+              // points are on the same edge of t2 --> we shorten an already cut edge
+              typename Kernel::FT alpha = compute_alpha(p1, q1,
+                                                        Pot::point_from_id(p2, q2, r2,  q.id2()),
+                                                        Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3));
+
+              return  Point_on_triangle<Kernel>(edge_id_t1, q.id2(), alpha);
+            }
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+            std::cout << " -- case 3\n";
+#endif
+            // point of t1: look for an edge of t1 containing both points
+            CGAL_assertion( p.extra_t1.count(q.id1())!=0 || p.extra_t1.count(3-q.id1()-edge_id_t1)!=0 );
+            int eid1 = p.extra_t1.count(q.id1())!=0 ? q.id1() : 3-q.id1()-edge_id_t1;
+            return  Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
+          }
+          // C: (-1, ip2) - (iq1, -1)
+          //vertex of t1, special case t1 edge passed thru a vertex of t2
+          CGAL_assertion(edge_id_t1 == 2);
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+          std::cout << " -- case 4\n";
+#endif
+          CGAL_assertion(q.id1()==1);
+          CGAL_assertion(!p.extra_t1.empty());
+          return  Point_on_triangle<Kernel>(p.extra_t1.count(0)==1?0:2,-1);
+      }
+    }
+    default:
+    {
+      switch(p.id2())
+      {
+        case -1:
+        {
+          switch(q.id1())
+          {
+            case -1: // G: (ip1, -1) - (-1, iq2)
+              //vertex of t1, special case t1 edge passed thru a vertex of t2
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+              std::cout << " -- case 5\n";
+#endif
+              CGAL_assertion(edge_id_t1 == 2);
+              CGAL_assertion(p.id1()==1);
+              CGAL_assertion(!q.extra_t1.empty());
+              return  Point_on_triangle<Kernel>(q.extra_t1.count(0)==1?0:2,-1);
+            default:
+            {
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+              std::cout << " -- case 6\n";
+#endif
+              CGAL_assertion(q.id2()!=-1); // I: (ip1, -1) - (iq2, -1)
+              //H: (ip1,-1), (iq1, iq2)
+              CGAL_assertion(edge_id_t1==2);
+              // p and q are on the same edge of t1
+              CGAL_assertion(p.id1()==q.id1() || p.id1()==(q.id1()+1)%3);
+              return  Point_on_triangle<Kernel>((q.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1);
+            }
+          }
+        }
+        default:
+        {
+          switch(q.id1())
+          {
+            case -1: // D: (ip1, ip2) - (-1, iq2)
+            {
+              if (q.id2() == p.id2() || q.id2() == (p.id2()+1)%3)
+              {
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+                std::cout << " -- case 7\n";
+#endif
+                // points are on the same edge of t2 --> we shorten an already cut edge
+                typename Kernel::FT alpha = compute_alpha(p1, q1,
+                                                          Pot::point_from_id(p2, q2, r2,  p.id2()),
+                                                          Pot::point_from_id(p2, q2, r2, (p.id2()+1)%3));
+
+                return  Point_on_triangle<Kernel>(edge_id_t1, p.id2(), alpha);
+              }
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+              std::cout << " -- case 8\n";
+#endif
+              // point of t1
+              //std::cout << "q.extra_t1: "; for(int qet1 : q.extra_t1) std::cout << " " << qet1; std::cout << "\n";
+              CGAL_assertion( q.extra_t1.count(p.id1())!=0 || q.extra_t1.count(3-p.id1()-edge_id_t1)!=0 );
+              int eid1 = q.extra_t1.count(p.id1())!=0 ? p.id1() : 3-p.id1()-edge_id_t1;
+              return  Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
+            }
+            default:
+            {
+              switch(q.id2())
+              {
+                case -1: // F: (ip1, ip2) - (iq1, -1)
+                {
+                  // p and q are on the same edge of t1
+                  CGAL_assertion(q.id1()==p.id1() ||  q.id1()==(p.id1()+1)%3);
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+                  std::cout << " -- case 9\n";
+#endif
+                  return  Point_on_triangle<Kernel>((p.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1);
+                }
+                default: // E: (ip1, ip2) - (iq1, iq2)
+                {
+                  if (p.id2()==q.id2())
+                  {
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+                    std::cout << " -- case 10\n";
+#endif
+                    typename Kernel::FT alpha = compute_alpha(p1, q1,
+                                                              Pot::point_from_id(p2, q2, r2,  q.id2()),
+                                                              Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3));
+                    return  Point_on_triangle<Kernel>(edge_id_t1, q.id2(), alpha);
+                  }
+                  // we are intersecting an edge of t1
+                  CGAL_assertion(p.id1()==q.id1() || edge_id_t1==2);
+                  int eid1 = p.id1()==q.id1() ? p.id1() : 1;
+                  return  Point_on_triangle<Kernel>((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+}
+
+template <class Kernel>
+void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p1,
+                                            const typename Kernel::Point_3& q1,
+                                            const typename Kernel::Point_3& r1,
+                                            int edge_id,
+                                            const typename Kernel::Point_3& p2,
+                                            const typename Kernel::Point_3& q2,
+                                            const typename Kernel::Point_3& r2,
+                                            const Kernel& k,
+                                            std::list<Point_on_triangle<Kernel>>& inter_pts)
+{
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "  cutoff using e" << edge_id << ": "
+            << to_double(p1.x()) << " " << to_double(p1.y()) << " " << to_double(p1.z()) << " "
+            << to_double(q1.x()) << " " << to_double(q1.y()) << " " << to_double(q1.z()) << "\n";
+#endif
+  typedef typename std::list<Point_on_triangle<Kernel>>::iterator Iterator;
 
   if(inter_pts.empty())
     return;
 
-  typename Kernel::Coplanar_orientation_3 orient = k.coplanar_orientation_3_object();
-  typename Kernel::Construct_line_3 line = k.construct_line_3_object();
+  //orient(p1,q1,r1,r1) is POSITIVE
+  std::map<const Point_on_triangle<Kernel>*,Orientation> orientations; // TODO skip map
+  for (Point_on_triangle<Kernel>& pot : inter_pts)
+  {
+    orientations[ &pot ]=pot.orientation(p1,q1,r1,edge_id,p2,q2,r2,k);
+    if (pot.id1()==-1 && orientations[ &pot ]==COLLINEAR)
+      pot.extra_t1.insert(edge_id);
+  }
 
-  //orient(p,q,r,r) is POSITIVE
-  std::map<const typename Kernel::Point_3*,Orientation> orientations;
-  for (Iterator it=inter_pts.begin();it!=inter_pts.end();++it)
-    orientations[ &(*it) ]=orient(p,q,r,*it);
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "    Orientations:";
+  for (const Point_on_triangle<Kernel>& pot : inter_pts)
+    std::cout << " " << orientations[ &pot ];
+  std::cout << "\n";
+#endif
+  CGAL_kernel_assertion_code(int pt_added = 0);
 
-  CGAL_kernel_assertion_code(int pt_added = 0;)
+  Iterator prev = std::prev(inter_pts.end());
 
-  const typename Kernel::Point_3* prev = &(*std::prev(inter_pts.end()));
   Iterator stop = inter_pts.size() > 2 ? inter_pts.end() : std::prev(inter_pts.end());
   for(Iterator it=inter_pts.begin(); it!=stop; ++it)
   {
-    const typename Kernel::Point_3& curr = *it;
-    Orientation or_prev = orientations[prev],
-                or_curr = orientations[&curr];
+    Orientation or_prev = orientations[&(*prev)],
+                or_curr = orientations[&(*it)];
 
     if((or_prev == POSITIVE && or_curr == NEGATIVE) ||
        (or_prev == NEGATIVE && or_curr == POSITIVE))
     {
-      typename Intersection_traits<Kernel, typename Kernel::Line_3, typename Kernel::Line_3>::result_type
-        obj = intersection(line(p,q), line(*prev,curr), k);
+      Point_on_triangle<Kernel> new_pt = intersection(*prev, *it, edge_id, p1, q1, p2, q2, r2, k);
 
-      // assert "not empty"
-      CGAL_kernel_assertion(bool(obj));
-
-      const typename Kernel::Point_3* inter = intersect_get<typename Kernel::Point_3>(obj);
-      CGAL_kernel_assertion(inter != nullptr);
-
-      prev = &(*inter_pts.insert(it,*inter));
-      orientations[prev] = COLLINEAR;
-      CGAL_kernel_assertion_code(++pt_added;)
+      prev = inter_pts.insert(it,new_pt);
+      orientations[&(*prev)] = COLLINEAR;
+      CGAL_kernel_assertion_code(++pt_added);
     }
 
-    prev = &(*it);
+    prev = it;
   }
 
   CGAL_kernel_assertion(pt_added<3);
@@ -96,35 +386,77 @@ intersection_coplanar_triangles(const typename K::Triangle_3& t1,
                                 const typename K::Triangle_3& t2,
                                 const K& k)
 {
-  const typename K::Point_3& p = t1.vertex(0),
-                             q = t1.vertex(1),
-                             r = t1.vertex(2);
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  auto to_string = [](const typename K::Triangle_3& t)
+  {
+    std::stringstream sstr;
+    sstr << "4 "
+         << to_double(t[0].x()) << " "  << to_double(t[0].y()) << " "  << to_double(t[0].z()) << " "
+         << to_double(t[1].x()) << " "  << to_double(t[1].y()) << " "  << to_double(t[1].z()) << " "
+         << to_double(t[2].x()) << " "  << to_double(t[2].y()) << " "  << to_double(t[2].z()) << " "
+         << to_double(t[0].x()) << " "  << to_double(t[0].y()) << " "  << to_double(t[0].z()) << "\n";
+    return sstr.str();
+  };
 
-  std::list<typename K::Point_3> inter_pts;
-  inter_pts.push_back(t2.vertex(0));
-  inter_pts.push_back(t2.vertex(1));
-  inter_pts.push_back(t2.vertex(2));
+  std::cout << "intersection_coplanar_triangles\n";
+  std::ofstream("/tmp/t1.polylines.txt") << std::setprecision(17) << to_string(t1) << "\n";
+  std::ofstream("/tmp/t2.polylines.txt") << std::setprecision(17) << to_string(t2) << "\n";
+#endif
+  const typename K::Point_3& p1 = t1.vertex(0),
+                             q1 = t1.vertex(1),
+                             r1 = t1.vertex(2);
 
+  const typename K::Point_3& p2 = t2.vertex(0),
+                             q2 = t2.vertex(1),
+                             r2 = t2.vertex(2);
+
+  std::list<Point_on_triangle<K>> inter_pts;
+  inter_pts.push_back(Point_on_triangle<K>(-1,0));
+  inter_pts.push_back(Point_on_triangle<K>(-1,1));
+  inter_pts.push_back(Point_on_triangle<K>(-1,2));
+
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  auto print_points = [&]()
+  {
+    for(auto p : inter_pts) std::cout << "  (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) "; std::cout <<"\n";
+  };
+  std::cout << "  ipts size: " << inter_pts.size() << "\n";
+  print_points();
+#endif
   //intersect t2 with the three half planes which intersection defines t1
-  intersection_coplanar_triangles_cutoff(p,q,r,k,inter_pts); //line pq
-  intersection_coplanar_triangles_cutoff(q,r,p,k,inter_pts); //line qr
-  intersection_coplanar_triangles_cutoff(r,p,q,k,inter_pts); //line rp
+  intersection_coplanar_triangles_cutoff(p1,q1,r1,0,p2,q2,r2,k,inter_pts); //line pq
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "  ipts size: " << inter_pts.size() << "\n";
+  print_points();
+#endif
+  intersection_coplanar_triangles_cutoff(q1,r1,p1,1,p2,q2,r2,k,inter_pts); //line qr
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "  ipts size: " << inter_pts.size() << "\n";
+  print_points();
+#endif
+  intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,inter_pts); //line rp
+#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION
+  std::cout << "  ipts size: " << inter_pts.size() << "\n";
+  print_points();
+#endif
 
+  auto point = [&](const Point_on_triangle<K>& pot){ return pot.point(p1,q1,r1,p2,q2,r2,k); };
   switch(inter_pts.size())
   {
     case 0:
       return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>();
     case 1:
-      return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(*inter_pts.begin());
+      return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(point(*inter_pts.begin()));
     case 2:
       return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
-            k.construct_segment_3_object()(*inter_pts.begin(), *std::next(inter_pts.begin())) );
+            k.construct_segment_3_object()(point(*inter_pts.begin()), point(*std::next(inter_pts.begin()))) );
     case 3:
       return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
-            k.construct_triangle_3_object()(*inter_pts.begin(), *std::next(inter_pts.begin()), *std::prev(inter_pts.end())) );
+            k.construct_triangle_3_object()(point(*inter_pts.begin()), point(*std::next(inter_pts.begin())), point(*std::prev(inter_pts.end()))) );
     default:
       return intersection_return<typename K::Intersect_3, typename K::Triangle_3, typename K::Triangle_3>(
-            std::vector<typename K::Point_3>(inter_pts.begin(),inter_pts.end()));
+            std::vector<typename K::Point_3>(boost::make_transform_iterator(inter_pts.begin(), point),
+                                             boost::make_transform_iterator(inter_pts.end(), point)));
   }
 }
 

Intersections_3/test/Intersections_3/triangle_3_triangle_3_intersection.cpp

@@ -96,6 +96,13 @@ void test_coplanar_triangles(){
   assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
+  // TK10 case C'
+  t1=Triangle(Point(88.7921, 89.0007, 1.25), Point(88.1912, 88.3997, 1.25), Point(89.8224, 90.031, 1.25));
+  t2=Triangle(Point(88.0497, 88.2583, 1.25), Point(82.9292, 81.8747, 1.25), Point(91.1726, 91.3812, 1.25));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Triangle>(&obj)!=nullptr);
 //Intersection is a point
   //edges  are collinear, one vertex in common
   t1=Triangle( Point(0,0,0),Point(0,1,0),Point(1,0,0) );
@@ -153,6 +160,13 @@ void test_coplanar_triangles(){
   assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
+  // TK10 case D
+  t1=Triangle(Point(-34.893700000000003, -16.0351, 3.1334899999999998e-12), Point(-34.893700000000003, -18.5351, 3.1334899999999998e-12), Point(-42.393700000000003, -16.0351, 3.1334899999999998e-12));
+  t2=Triangle(Point(-34.893700000000003, -32.0351, 3.1334899999999998e-12), Point(-34.893700000000003, -9.7851400000000002, 3.1334899999999998e-12), Point(-31.643699999999999, -17.201799999999999, 3.1334899999999998e-12));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Segment>(&obj)!=nullptr);
 //Intersection is a polygon
   //David's star
   t1=Triangle( Point(0,0,0),Point(1,0,0),Point(0.5,1.5,0) );
@@ -181,6 +195,51 @@ void test_coplanar_triangles(){
   obj=CGAL::intersection(t2,t1);
   assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
   assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case A
+  t1=Triangle(Point(3.74861, 12.4822, 14.0112), Point(5.40582, 12.4822, 15.6895), Point(5.37748, 12.4822, 15.7206));
+  t2=Triangle(Point(5.49972, 12.4822, 13.491), Point(5.27627, 12.4822, 15.8106), Point(5.32119, 12.4822, 15.8126));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case C
+  t1=Triangle(Point(5, -94.6659, 3.85175), Point(5, -94.5682, 3.08638), Point(5, -94.8182, 3.08638));
+  t2=Triangle(Point(5, -94.4317, 3.76399), Point(5, -97.6182, 3.08638), Point(5, -94.5659, 2.99682));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case E
+  t1=Triangle(Point(-955.858, -45.032, -0.016), Point(-955.856, -45.032, -0.004), Point(-955.856, -45.032, -0.002));
+  t2=Triangle(Point(-955.856, -45.032, 0.006), Point(-955.854, -45.032, -0.002), Point(-955.876, -45.032, -0.034));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case F
+  t1=Triangle(Point(141.172, 20.576, 155.764), Point(141.172, 20.588, 155.766), Point(141.172, 20.59, 155.766));
+  t2=Triangle(Point(141.172, 20.602, 155.768), Point(141.172, 20.594, 155.766), Point(141.172, 20.574, 155.764));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  // TK10 case D
+  t1=Triangle(Point(152.864, 126.324, 0.950001), Point(152.77, 126.483, 0.950001), Point(153.072, 125.973, 0.950001));
+  t2=Triangle(Point(153.322, 125.551, 0.950001), Point(152.218, 127.415, 0.950001), Point(153.66, 124.768, 0.950001));
+  obj=CGAL::intersection(t1,t2);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
+  obj=CGAL::intersection(t2,t1);
+  assert(CGAL::object_cast<Polygon2>(&obj)!=nullptr);
+  assert(CGAL::object_cast<Polygon2>(&obj)->size()==4);
 //Intersection is empty
   t1=Triangle( Point(0,0,0),Point(0,1,0),Point(1,0,0) );
   t2=Triangle( Point(-0.1,-0.1,0),Point(-0.1,-0.9,0),Point(-1,-0.1,0) );

Kernel_23/include/CGAL/Kernel/function_objects.h

@@ -978,11 +978,42 @@ namespace CommonKernelFunctors {
      const Vector_3 ad = vector(a,d);
 
      const Vector_3 abad = cross_product(ab,ad);
-     const double x = CGAL::to_double(scalar_product(cross_product(ab,ac), abad));
-     const double l_ab = CGAL::sqrt(CGAL::to_double(sq_distance(a,b)));
-     const double y = l_ab * CGAL::to_double(scalar_product(ac,abad));
+     const Vector_3 abac = cross_product(ab,ac);
 
-     return FT(std::atan2(y, x) * 180 / CGAL_PI );
+     // The dihedral angle we are interested in is the angle around the oriented
+     // edge ab which is the same (in absolute value) as the angle between the
+     // vectors ab^ac and ab^ad (cross-products).
+     // (abac points inside the tetra abcd if its orientation is positive and outside otherwise)
+     //
+     // We consider the vector abad in the basis defined by the three vectors
+     //    (<ab>, <abac>, <ab^abac>)
+     // where <u> denote the normalized vector u/|u|.
+     //
+     // In this orthonormal basis, the vector adab has the coordinates
+     //    x = <ab>      * abad
+     //    y = <abac>    * abad
+     //    z = <ab^abac> * abad
+     // We have x == 0, because abad and ab are orthogonal, and thus abad is in
+     // the plane (yz) of the new basis.
+     //
+     // In that basis, the dihedral angle is the angle between the y axis and abad
+     // which is the arctan of y/z, or atan2(z, y).
+     //
+     // (Note that ab^abac is in the plane abc, pointing outside the tetra if
+     //  its orientation is positive and inside otherwise).
+     //
+     // For the normalization, abad appears in both scalar products
+     // in the quotient so we can ignore its norm. For the second
+     // terms of the scalar products, we are left with ab^abac and abac.
+     // Since ab and abac are orthogonal, the sinus of the angle between the
+     // two vectors is 1.
+     // So the norms are |ab|.|abac| vs |abac|, which is why we have a
+     // multiplication by |ab| in y below.
+     const double l_ab = CGAL::sqrt(CGAL::to_double(sq_distance(a,b)));
+     const double y = l_ab * CGAL::to_double(scalar_product(abac, abad));
+     const double z = CGAL::to_double(scalar_product(cross_product(ab,abac),abad));
+
+     return FT(std::atan2(z, y) * 180 / CGAL_PI );
    }
  };
 
@@ -2193,6 +2224,107 @@ namespace CommonKernelFunctors {
     }
   };
 
+  template <typename K>
+  class Construct_planes_intersection_point_3
+  {
+    typedef typename K::Plane_3 Plane;
+    typedef typename K::Point_3 Point;
+    typename K::Construct_plane_3 construct_plane;
+  public:
+    typedef Point result_type;
+
+    Point
+    operator()(const Point& p1, const Point& q1, const Point& r1,
+               const Point& p2, const Point& q2, const Point& r2,
+               const Point& p3, const Point& q3, const Point& r3) const
+    {
+      Plane plane1 = construct_plane(p1, q1, r1);
+      Plane plane2 = construct_plane(p2, q2, r2);
+      Plane plane3 = construct_plane(p3, q3, r3);
+
+      const auto res = typename K::Intersect_3()(plane1, plane2, plane3);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+
+    Point
+    operator()(const Plane& plane1, const Plane& plane2, const Plane& plane3) const
+    {
+      const auto res = typename K::Intersect_3()(plane1, plane2, plane3);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+  };
+
+  template <typename K>
+  class Construct_coplanar_segments_intersection_point_3
+  {
+    typedef typename K::Segment_3 Segment;
+    typedef typename K::Point_3 Point;
+    typename K::Construct_segment_3 construct_segment;
+  public:
+    typedef Point result_type;
+
+    Point
+    operator()(const Point& p1, const Point& q1,
+               const Point& p2, const Point& q2) const
+    {
+      Segment s1 = construct_segment(p1, q1);
+      Segment s2 = construct_segment(p2, q2);
+
+      const auto res = typename K::Intersect_3()(s1, s2);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+
+    Point
+    operator()(const Segment& s1, const Segment& s2) const
+    {
+      const auto res = typename K::Intersect_3()(s1, s2);
+      CGAL_assertion(res!=std::nullopt);
+      const Point* e_pt = std::get_if<Point>(&(*res));
+      CGAL_assertion(e_pt!=nullptr);
+      return *e_pt;
+    }
+  };
+
+  template <typename K>
+  class Compute_alpha_for_coplanar_triangle_intersection_3
+  {
+    typedef typename K::Point_3 Point_3;
+    typedef typename K::Vector_3 Vector_3;
+  public:
+    typedef typename K::FT  result_type;
+    result_type
+    operator()(const Point_3& p1, const Point_3& p2,       // segment 1
+               const Point_3& p3, const Point_3& p4) const // segment 2
+    {
+      typename K::Construct_vector_3 vector = K().construct_vector_3_object();
+      typename K::Construct_cross_product_vector_3 cross_product =
+        K().construct_cross_product_vector_3_object();
+
+      const Vector_3 v1 = vector(p1, p2);
+      const Vector_3 v2 = vector(p3, p4);
+
+      CGAL_assertion(K().coplanar_3_object()(p1,p2,p3,p4));
+
+      const Vector_3 v3 = vector(p1, p3);
+      const Vector_3 v3v2 = cross_product(v3,v2);
+      const Vector_3 v1v2 = cross_product(v1,v2);
+      const typename K::FT sl = K().compute_squared_length_3_object()(v1v2);
+      CGAL_assertion(!certainly(is_zero(sl)));
+
+      const typename K::FT t = ((v3v2.x()*v1v2.x()) + (v3v2.y()*v1v2.y()) + (v3v2.z()*v1v2.z())) / sl;
+      return t; // p1 + (p2-p1) * t
+    }
+  };
+
   template <typename K>
   class Construct_point_on_2
   {

Kernel_23/include/CGAL/Kernel/interface_macros.h

@@ -398,6 +398,12 @@ CGAL_Kernel_cons(Construct_plane_3,
                  construct_plane_3_object)
 CGAL_Kernel_cons(Construct_plane_line_intersection_point_3,
                  construct_plane_line_intersection_point_3_object)
+CGAL_Kernel_cons(Construct_planes_intersection_point_3,
+                 construct_planes_intersection_point_3_object)
+CGAL_Kernel_cons(Construct_coplanar_segments_intersection_point_3,
+                 construct_coplanar_segments_intersection_point_3_object)
+CGAL_Kernel_cons(Compute_alpha_for_coplanar_triangle_intersection_3,
+                 compute_alpha_for_coplanar_triangle_intersection_3_object)
 CGAL_Kernel_cons(Construct_point_on_2,
                  construct_point_on_2_object)
 CGAL_Kernel_cons(Construct_point_on_3,

Kernel_23/include/CGAL/Kernel_traits.h

@@ -19,7 +19,6 @@
 
 #include <CGAL/Kernel_traits_fwd.h>
 #include <boost/mpl/has_xxx.hpp>
-#include <boost/mpl/if.hpp>
 
 namespace CGAL {
 

Kernel_23/test/Kernel_23/include/CGAL/_test_io.h

@@ -17,26 +17,19 @@
 #ifndef CGAL__TEST_IO_H
 #define CGAL__TEST_IO_H
 
-#include <fstream>
+#include <sstream>
 #include <cassert>
 
-#ifndef TEST_FILENAME
-#  define TEST_FILENAME "Test_IO.out"
-#endif
-
 template <class T>
 void
 _test_io_for(const T& t)
 {
-    {
-        std::ofstream oFile(TEST_FILENAME, std::ios::out);
-        oFile << t << std::endl;
-    }
+    std::stringstream ss;
+    ss << t << std::endl;
 
-    std::ifstream iFile(TEST_FILENAME, std::ios::in);
     T u = t;
-    iFile >> u;
-    assert(!iFile.fail());
+    ss >> u;
+    assert(! ss.fail() );
     assert(u == t);
 }
 

Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp

@@ -10,10 +10,49 @@ struct query {
   double expected_angle;
 };
 
+void sign_test()
+{
+  K::Point_3 a(0,0,0), b(1,0,0), c(0,1, 0), d(0,0,1);
+
+  assert( CGAL::approximate_dihedral_angle(a, b, c, d) > 0);
+  assert( CGAL::approximate_dihedral_angle(c, a, b, d) > 0);
+  assert( CGAL::approximate_dihedral_angle(a, d, b, c) > 0);
+  assert( CGAL::approximate_dihedral_angle(c, b, d, a) > 0);
+  assert( CGAL::approximate_dihedral_angle(d, b, a, c) > 0);
+  assert( CGAL::approximate_dihedral_angle(d, c, b, a) > 0);
+
+  assert( CGAL::approximate_dihedral_angle(a, b, d, c) < 0);
+  assert( CGAL::approximate_dihedral_angle(c, a, d, b) < 0);
+  assert( CGAL::approximate_dihedral_angle(a, d, c, b) < 0);
+  assert( CGAL::approximate_dihedral_angle(c, b, a, d) < 0);
+  assert( CGAL::approximate_dihedral_angle(d, b, c, a) < 0);
+  assert( CGAL::approximate_dihedral_angle(d, c, a, b) < 0);
+}
+
+auto almost_equal_angle(double a, double b) {
+  return (std::min)(std::abs(a - b), std::abs(a + 360 - b)) < 0.1;
+}
+
+void test_regular_tetrahedron()
+{
+  auto half_root_of_2 = std::sqrt(2) / 2;
+
+  // Regular tetrahedron
+  Point_3 a{ -1,  0, -half_root_of_2};
+  Point_3 b{  1,  0, -half_root_of_2};
+  Point_3 c{  0,  1,  half_root_of_2};
+  Point_3 d{  0, -1,  half_root_of_2};
+  assert(orientation(a, b, c, d) == CGAL::POSITIVE);
+  assert(almost_equal_angle(CGAL::approximate_dihedral_angle(a, b, c, d), 70.5288));
+}
+
 int main() {
+  std::cout.precision(17);
+  sign_test();
+  test_regular_tetrahedron();
+
   Point_3 a = {0, 0, 0};
-  Point_3 b = {0, 1, 0};
-  Point_3 c = {1, 0, 0};
+  Point_3 b = {0, -1, 0}; // ab is oriented so that it sees the plan xz positively.
 
   const query queries[] = {
     { {  1, 0,  0},    0.},
@@ -26,11 +65,28 @@ int main() {
     { {  1, 0, -1},  -45.},
   };
 
-  for(auto query: queries) {
-    const auto& expected = query.expected_angle;
-    const auto& p = query.p;
-    auto approx = CGAL::approximate_dihedral_angle(a, b, c, p);
-    std::cout << approx << "  -- " << expected << '\n';
-    assert( std::abs(approx - expected) < 0.1 );
+  auto cnt = 0u;
+  for(double yc = -10; yc < 10; yc += 0.1) {
+    // c can be any point in the half-plane xy, with x>0
+    Point_3 c{1, yc, 0};
+    // std::cout << "c = " << c << '\n';
+    for(const auto& query : queries) {
+      for(double yp = -10; yp < 10; yp += 0.3) {
+        const auto& expected = query.expected_angle;
+        const Point_3 p{query.p.x(), yp, query.p.z()};
+        // std::cout << "p = " << p << '\n';
+        auto approx = CGAL::approximate_dihedral_angle(a, b, c, p);
+        // std::cout << approx << "  -- " << expected << '\n';
+        if(!almost_equal_angle(approx, expected)) {
+          std::cout << "ERROR:\n";
+          std::cout << "CGAL::approximate_dihedral_angle(" << a << ", " << b << ", " << c << ", " << p << ") = " << approx << '\n';
+          std::cout << "expected: " << expected << '\n';
+          return 1;
+        }
+        ++cnt;
+      }
+    }
   }
+  std::cout << "OK (" << cnt << " tests)\n";
+  assert(cnt > 10000);
 }

Maintenance/test_handling/create_testresult_page

@@ -474,7 +474,7 @@ sub print_platform_descriptions()
 <th>BOOST</th>
 <th>MPFR</th>
 <th>GMP</th>
-<th>QT5</th>
+<th>QT</th>
 <th>LEDA</th>
 <th>CXXFLAGS</th>
 <th>LDFLAGS</th>
@@ -544,6 +544,7 @@ EOF
             print OUTPUT ">$pf_short</a>";
 	    my $index = 12;
 	    while ($index) {
+        $index--;
 		print OUTPUT "<td>?</td>\n";
 	    }
 	}

Maintenance/test_handling/to_zipped_format

@@ -51,7 +51,7 @@ sub reformat_results($)
     $_ = $line;
     open (PLATFORM_INFO,">${platform}.info") or return;
     open (PLATFORM_NEW_RESULTS,">${platform}.new_results") or return;
-    my ($CGAL_VERSION,$LEDA_VERSION,$COMPILER,$TESTER_NAME,$TESTER_ADDRESS,$GMP,$MPFR,$ZLIB,$OPENGL,$BOOST,$QT,$QT4,$QT5,$CMAKE) = ("-","-","-","-","-","-","-","-","-","-","-","-","-","-","-","no");
+    my ($CGAL_VERSION,$LEDA_VERSION,$COMPILER,$TESTER_NAME,$TESTER_ADDRESS,$GMP,$MPFR,$ZLIB,$OPENGL,$BOOST,$QT,$CMAKE) = ("-","-","-","-","-","-","-","-","-","-","-","-","-","no");
     my ($LDFLAGS,$CXXFLAGS) = ("", "");
     while (! /^------/) {
         if(/^\s*$/) {
@@ -97,10 +97,13 @@ sub reformat_results($)
 	    $QT="$1";
 	}
 	if (/QT4_VERSION = '([^']+)'/) {
-	    $QT4="$1";
+	    $QT="$1";
 	}
 	if (/Qt5_VERSION = '([^']+)'/) {
-	    $QT5="$1";
+	    $QT="$1";
+	}
+	if (/Qt6_VERSION = '([^']+)'/) {
+	    $QT="$1";
 	}
 	if (/BOOST_VERSION = '([^']+)'/) {
 	    $BOOST="$1";
@@ -147,7 +150,7 @@ $CMAKE
 $BOOST
 $MPFR
 $GMP
-$QT5
+$QT
 $LEDA_VERSION
 $CXXFLAGS
 $LDFLAGS

Mesh_3/benchmark/Mesh_3/StdAfx.h

@@ -42,11 +42,9 @@
 #include <boost/mpl/clear.hpp>
 #include <boost/mpl/contains.hpp>
 #include <boost/mpl/end.hpp>
-#include <boost/mpl/eval_if.hpp>
 #include <boost/mpl/has_key.hpp>
 #include <boost/mpl/has_xxx.hpp>
 #include <boost/mpl/identity.hpp>
-#include <boost/mpl/if.hpp>
 #include <boost/mpl/insert.hpp>
 #include <boost/mpl/insert_fwd.hpp>
 #include <boost/mpl/iterator_range.hpp>