Merge remote-tracking branch 'cgal/master' into PMP-BF_autorefine

# Conflicts:
#	Installation/CHANGES.md
#	Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

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>

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
 

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"
 }
 

Installation/CHANGES.md

@@ -44,6 +44,7 @@ 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

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 );
    }
  };
 

Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp

@@ -10,10 +10,51 @@ 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.
+  [[maybe_unused]] Point_3 c = {1, 0, 0};
+  // c can be any point in the half-plane xy, with x>0
 
   const query queries[] = {
     { {  1, 0,  0},    0.},
@@ -26,11 +67,27 @@ 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) {
+    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);
 }

Mesh_3/include/CGAL/Mesh_criteria_3.h

@@ -70,22 +70,15 @@ public:
   // is not a problem
   template <typename CGAL_NP_TEMPLATE_PARAMETERS>
   Mesh_criteria_3_impl(const CGAL_NP_CLASS& np)
-    :edge_criteria_(parameters::choose_parameter(parameters::get_parameter(np, internal_np::edge_size_param),
-                                                 parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::edge_sizing_field_param),
-                                                                              parameters::choose_parameter(parameters::get_parameter(np, internal_np::sizing_field_param), FT(DBL_MAX)))),
+    :edge_criteria_(parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::edge_size_param), FT(DBL_MAX)),
                     parameters::choose_parameter(parameters::get_parameter(np, internal_np::edge_min_size_param), FT(0))),
     facet_criteria_(parameters::choose_parameter(parameters::get_parameter(np, internal_np::facet_angle_param), FT(0)),
-                    parameters::choose_parameter(parameters::get_parameter(np, internal_np::facet_size_param),
-                                                 parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::facet_sizing_field_param),
-                                                                              parameters::choose_parameter(parameters::get_parameter(np, internal_np::sizing_field_param), FT(0)))),
+                    parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::facet_size_param), FT(0)),
                     parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::facet_distance_param), FT(0)),
                     parameters::choose_parameter(parameters::get_parameter(np, internal_np::facet_topology_param), CGAL::FACET_VERTICES_ON_SURFACE),
                     parameters::choose_parameter(parameters::get_parameter(np, internal_np::facet_min_size_param), FT(0))),
-    cell_criteria_(parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_radius_edge_ratio_param),
-                                                parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_radius_edge_param), FT(0))),
-                   parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_size_param),
-                                                parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::cell_sizing_field_param),
-                                                                             parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::sizing_field_param), FT(0)))),
+    cell_criteria_(parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_radius_edge_ratio_param), FT(0)),
+                   parameters::choose_parameter(parameters::get_parameter_reference(np, internal_np::cell_size_param), FT(0)),
                    parameters::choose_parameter(parameters::get_parameter(np, internal_np::cell_min_size_param), FT(0)))
   { }
 

Mesh_3/test/Mesh_3/test_mesh_capsule_var_distance_bound.cpp

@@ -58,7 +58,9 @@ int main()
                          facet_distance=field);
 
   // Mesh generation
-  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
+                                      perturb(time_limit = 0, sliver_bound = 0),
+                                      exude(time_limit = 0, sliver_bound = 0));
 
 //  // Output
 //  std::ofstream medit_file("out.mesh");

Mesh_3/test/Mesh_3/test_mesh_criteria_creation.cpp

@@ -75,28 +75,12 @@ int main()
   Mc ec1(edge_size = 1);
   assert( ec1.edge_criteria_object().sizing_field(bp1,1,index) == 1 );
 
-  Mc ec2(edge_sizing_field = Esf(2));
+  Mc ec2(edge_size = Esf(2));
   assert( ec2.edge_criteria_object().sizing_field(bp1,1,index) == 2 );
 
-  Mc ec3(edge_sizing_field = 3.);
+  Mc ec3(edge_size = 3.);
   assert( ec3.edge_criteria_object().sizing_field(bp1,1,index) == 3 );
 
-  Mc ec4(edge_size = 4.,
-         edge_sizing_field = Esf(5.));
-  assert( ec4.edge_criteria_object().sizing_field(bp1,1,index) == 4. );
-
-  Mc ec5(sizing_field = 5.);
-  assert( ec5.edge_criteria_object().sizing_field(bp1,1,index) == 5 );
-
-  Mc ec6(sizing_field = 6.,
-         edge_sizing_field = 7.);
-  assert( ec6.edge_criteria_object().sizing_field(bp1,1,index) == 7. );
-
-  Mc ec7(sizing_field = 7.,
-         edge_size = 8.);
-  assert( ec7.edge_criteria_object().sizing_field(bp1,1,index) == 8. );
-
-
   // -----------------------------------
   // Test facet criteria
   // -----------------------------------
@@ -109,48 +93,29 @@ int main()
                                               cp(tr.point(ch, k+3))));
 
   FT facet_size_ok = radius_facet*FT(10);
-  FT facet_size_nok = radius_facet/FT(10);
 
   Mc fc1(facet_size = facet_size_ok);
   assert( ! fc1.facet_criteria_object()(tr, f) );
 
-  Mc fc2(facet_sizing_field = facet_size_ok);
-  assert( ! fc2.facet_criteria_object()(tr, f) );
-
-  Mc fc3(facet_sizing_field = Fsf(facet_size_ok));
+  Mc fc3(facet_size = Fsf(facet_size_ok));
   assert( ! fc3.facet_criteria_object()(tr, f) );
 
-  Mc fc4(facet_sizing_field = facet_size_nok,
-         facet_size = facet_size_ok);
-  assert( ! fc4.facet_criteria_object()(tr, f) );
-
-  Mc fc5(sizing_field = facet_size_ok);
-  assert( ! fc5.facet_criteria_object()(tr, f) );
-
-  Mc fc6(facet_size = facet_size_ok,
-         facet_sizing_field = facet_size_nok,
-         sizing_field = facet_size_nok);
-  assert( ! fc6.facet_criteria_object()(tr, f) );
-
-  Mc fc7(facet_sizing_field = Fsf(facet_size_ok),
-         sizing_field = facet_size_nok);
-  assert( ! fc7.facet_criteria_object()(tr, f) );
-
   Mc fc8(facet_distance = 8.);
   Mc fc9(facet_angle = 9.);
   Mc fc10(facet_angle = 10.1,
           facet_distance = 10.2,
           facet_size = 10.3,
-          facet_min_size = 0.2,
-          facet_sizing_field = Fsf(10.4),
-          sizing_field = 10.5);
+          facet_min_size = 0.2);
+  Mc fc1Ob(facet_angle = 10.11,
+          facet_distance = 10.12,
+          facet_min_size = 0.3,
+          facet_size = Fsf(10.14));
 
   // Test construction from int
   Mc fc11(facet_size = 11);
   Mc fc11b(facet_size = 11, facet_min_size = 1);
-  Mc fc12(facet_sizing_field = 12);
-  Mc fc12b(facet_sizing_field = 12, facet_min_size = 2);
-  Mc fc13(sizing_field = 13);
+  Mc fc12(facet_size = Fsf(12));
+  Mc fc12b(facet_size = Fsf(12), facet_min_size = 1);
 
   // Test topological criterion creation
   Mc fc14(facet_topology = CGAL::FACET_VERTICES_ON_SURFACE);
@@ -168,46 +133,27 @@ int main()
                                              cp(tr.point(ch, 3))));
 
   FT cell_size_ok = radius_cell*FT(10);
-  FT cell_size_nok = radius_cell/FT(10);
 
   Mc cc1(cell_size = cell_size_ok);
   assert( ! cc1.cell_criteria_object()(tr, ch) );
 
-  Mc cc2(cell_sizing_field = cell_size_ok);
-  assert( ! cc2.cell_criteria_object()(tr, ch) );
-
-  Mc cc3(cell_sizing_field = Fsf(cell_size_ok));
+  Mc cc3(cell_size = Fsf(cell_size_ok));
   assert( ! cc3.cell_criteria_object()(tr, ch) );
 
-  Mc cc4(cell_sizing_field = cell_size_nok,
-         cell_size = cell_size_ok);
-  assert( ! cc4.cell_criteria_object()(tr, ch) );
-
-  Mc cc5(sizing_field = cell_size_ok);
-  assert( ! cc5.cell_criteria_object()(tr, ch) );
-
-  Mc cc6(cell_size = cell_size_ok,
-         cell_sizing_field = cell_size_nok,
-         sizing_field = cell_size_nok);
-  assert( ! cc6.cell_criteria_object()(tr, ch) );
-
-  Mc cc7(cell_sizing_field = Csf(cell_size_ok),
-         sizing_field = cell_size_nok);
+  Mc cc7(cell_size = Csf(cell_size_ok));
   assert( ! cc7.cell_criteria_object()(tr, ch) );
 
   Mc cc8(cell_radius_edge_ratio = 8.);
   Mc cc9(cell_radius_edge_ratio = 9.1,
-         sizing_field = Csf(9.2) );
+         cell_size = Csf(9.2) );
   Mc cc10(cell_radius_edge_ratio = 10.1,
           cell_size = 10.2,
+          cell_min_size = 0.1);
+  Mc cc10b(cell_radius_edge_ratio = 10.1,
           cell_min_size = 0.1,
-          cell_sizing_field = Csf(10.3),
-          sizing_field = 10.4);
+          cell_size = Csf(10.3));
 
   // Test construction from int
   Mc cc11(cell_size = 11);
   Mc cc11b(cell_size = 11, cell_min_size = 1);
-  Mc cc12(cell_sizing_field = 12);
-  Mc cc12b(cell_sizing_field = 12, cell_min_size = 2);
-  Mc cc13(sizing_field = 13);
 }

Nef_S2/include/CGAL/Nef_S2/Sphere_segment.h

@@ -33,7 +33,10 @@ template <class R_> class Sphere_segment_rep
   typedef Sphere_segment_rep<R_> Rep;
   friend class Sphere_segment<R_>;
 public:
-Sphere_segment_rep() { ps_ = pt_ = Point(); c_ = Circle(); }
+
+Sphere_segment_rep() :
+  ps_(), pt_(), c_()
+{}
 
 Sphere_segment_rep(const Point& p1, const Point& p2,
                    bool shorter_arc=true) :

Periodic_3_mesh_3/test/Periodic_3_mesh_3/test_implicit_shapes_bunch.cpp

@@ -230,7 +230,7 @@ int main()
     Periodic_mesh_criteria criteria(facet_angle = 30,
                                     facet_size = 0.05,
                                     facet_distance = 0.025,
-                                    cell_radius_edge = 2,
+                                    cell_radius_edge_ratio = 2,
                                     cell_size = 0.05);
 
     // Mesh generation

Poisson_surface_reconstruction_3/include/CGAL/Reconstruction_triangulation_3.h

@@ -26,6 +26,7 @@
 #include <CGAL/Delaunay_triangulation_3.h>
 #include <CGAL/Delaunay_triangulation_cell_base_3.h>
 #include <CGAL/Triangulation_cell_base_with_info_3.h>
+#include <CGAL/Triangulation_structural_filtering_traits.h>
 
 #include <CGAL/algorithm.h>
 #include <CGAL/bounding_box.h>
@@ -148,6 +149,12 @@ struct Reconstruction_triangulation_default_geom_traits_3 : public BaseGt
 };
 
 
+template < class BaseGt >
+struct Triangulation_structural_filtering_traits<Reconstruction_triangulation_default_geom_traits_3<BaseGt> > {
+  typedef typename Triangulation_structural_filtering_traits<BaseGt>::Use_structural_filtering_tag  Use_structural_filtering_tag;
+};
+
+
 /// \internal
 /// The Reconstruction_triangulation_3 class
 /// provides the interface requested by the Poisson_reconstruction_function class:

Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSizingField.h

@@ -16,6 +16,7 @@ class PMPSizingField{
 public:
 
 /// @name Types
+///  These types are used for the documentation of the functions of the concept and not needed implementation wise.
 /// @{
 
 /// Vertex descriptor type
@@ -38,27 +39,30 @@ typedef unspecified_type FT;
 /// @name Functions
 /// @{
 
-/// a function that returns the sizing value at `v`.
-FT at(const vertex_descriptor v) const;
+/// returns the sizing value at `v` (used during tangential relaxation).
+FT at(const vertex_descriptor v, const PolygonMesh& pmesh) const;
 
-/// a function controlling edge split and edge collapse,
-/// returning the ratio of the current edge length and the local target edge length between
-/// the points of `va` and `vb` in case the current edge is too long, and `std::nullopt` otherwise.
+/// returns the ratio of the current edge squared length and the local target edge squared length between
+/// the points of `va` and `vb` in case the current edge is too long, and `std::nullopt` otherwise
+/// (used for triggering edge splits and preventing some edge collapses).
 std::optional<FT> is_too_long(const vertex_descriptor va,
-                              const vertex_descriptor vb) const;
+                              const vertex_descriptor vb,
+                              const PolygonMesh& pmesh) const;
 
-/// a function controlling edge collapse by returning the ratio of the squared length of `h` and the
-/// local target edge length if it is too short, and `std::nullopt` otherwise.
+/// returns the ratio of the squared length of `h` and the
+/// local target edge squared length if it is too short, and `std::nullopt` otherwise
+/// (used for triggering edge collapses).
 std::optional<FT> is_too_short(const halfedge_descriptor h,
                                const PolygonMesh& pmesh) const;
 
-/// a function returning the location of the split point of the edge of `h`.
+/// returns the position of the new vertex created when splitting the edge of `h`.
 Point_3 split_placement(const halfedge_descriptor h,
                         const PolygonMesh& pmesh) const;
 
-/// a function that updates the sizing field value at the vertex `v`.
-void update(const vertex_descriptor v,
-            const PolygonMesh& pmesh);
+/// function called after the addition of the split vertex `v` in `pmesh`.
+/// This function can be used for example to update a pre-computed sizing field.
+void register_split_vertex(const vertex_descriptor v,
+                           const PolygonMesh& pmesh);
 
 /// @}
 };

Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt

@@ -270,6 +270,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - `CGAL::Polygon_mesh_processing::triangle()`
 - `CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces()`
 - `CGAL::Polygon_mesh_processing::detect_corners_of_regions()`
+- `CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel()`
 
 \cgalCRPSection{I/O Functions}
 - \link PMP_IO_grp `CGAL::Polygon_mesh_processing::IO::read_polygon_mesh()`\endlink

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -51,6 +51,7 @@ create_single_source_cgal_program("match_faces.cpp")
 create_single_source_cgal_program("cc_compatible_orientations.cpp")
 create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp")
 create_single_source_cgal_program("hausdorff_bounded_error_distance_example.cpp")
+create_single_source_cgal_program("isotropic_remeshing_with_custom_sizing_example.cpp")
 create_single_source_cgal_program("triangle_mesh_autorefinement.cpp")
 create_single_source_cgal_program("soup_autorefinement.cpp")
 
@@ -77,6 +78,8 @@ if(TARGET CGAL::Eigen3_support)
   target_link_libraries(delaunay_remeshing_example PUBLIC CGAL::Eigen3_support)
   create_single_source_cgal_program("remesh_almost_planar_patches.cpp")
   target_link_libraries(remesh_almost_planar_patches PUBLIC CGAL::Eigen3_support)
+  create_single_source_cgal_program("geodesic_isolevel_refinement.cpp")
+  target_link_libraries(geodesic_isolevel_refinement PUBLIC CGAL::Eigen3_support)
   create_single_source_cgal_program("interpolated_corrected_curvatures_SM.cpp")
   target_link_libraries(interpolated_corrected_curvatures_SM PUBLIC CGAL::Eigen3_support)
   create_single_source_cgal_program("interpolated_corrected_curvatures_PH.cpp")

Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example.cpp

@@ -7,14 +7,14 @@
 #include <iostream>
 #include <string>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel       K;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 
-typedef K::Point_3                                                Point;
-typedef K::Vector_3                                               Vector;
+typedef K::Point_3                                          Point;
+typedef K::Vector_3                                         Vector;
 
-typedef CGAL::Surface_mesh<Point>                                 Surface_mesh;
-typedef boost::graph_traits<Surface_mesh>::vertex_descriptor      vertex_descriptor;
-typedef boost::graph_traits<Surface_mesh>::face_descriptor        face_descriptor;
+typedef CGAL::Surface_mesh<Point>                           Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor        vertex_descriptor;
+typedef boost::graph_traits<Mesh>::face_descriptor          face_descriptor;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -22,7 +22,7 @@ int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
 
-  Surface_mesh mesh;
+  Mesh mesh;
   if(!PMP::IO::read_polygon_mesh(filename, mesh))
   {
     std::cerr << "Invalid input." << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp

@@ -16,9 +16,9 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel   K;
 typedef K::Point_3                                            Point;
 typedef K::Vector_3                                           Vector;
 
-typedef CGAL::Polyhedron_3<K>                                 Polyhedron;
-typedef boost::graph_traits<Polyhedron>::vertex_descriptor    vertex_descriptor;
-typedef boost::graph_traits<Polyhedron>::face_descriptor      face_descriptor;
+typedef CGAL::Polyhedron_3<K>                                 Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor          vertex_descriptor;
+typedef boost::graph_traits<Mesh>::face_descriptor            face_descriptor;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -26,7 +26,7 @@ int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
 
-  Polyhedron mesh;
+  Mesh mesh;
   if(!PMP::IO::read_polygon_mesh(filename, mesh))
   {
     std::cerr << "Invalid input." << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_LCC.cpp

@@ -12,7 +12,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel   Kernel;
 typedef Kernel::Point_3                                       Point;
 
 typedef CGAL::Linear_cell_complex_traits<3, Kernel>           MyTraits;
-typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper<2, 3, MyTraits>::type LCC;
+typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper<2, 3, MyTraits>::type Mesh;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -21,7 +21,7 @@ int main(int argc, char* argv[])
   const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
   const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off");
 
-  LCC mesh1, mesh2;
+  Mesh mesh1, mesh2;
   if(!PMP::IO::read_polygon_mesh(filename1, mesh1) || !PMP::IO::read_polygon_mesh(filename2, mesh2))
   {
     std::cerr << "Invalid input." << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/extrude.cpp

@@ -8,13 +8,13 @@
 #include <fstream>
 #include <string>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel  Kernel;
-typedef CGAL::Surface_mesh<Kernel::Point_3> SM;
-typedef boost::graph_traits<SM>::vertex_descriptor           vertex_descriptor;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel    Kernel;
+typedef CGAL::Surface_mesh<Kernel::Point_3>                    Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor           vertex_descriptor;
 
-typedef Kernel::Vector_3                                     Vector;
-typedef boost::property_map<SM, CGAL::vertex_point_t>::type  VPMap;
-typedef SM::template Property_map<vertex_descriptor, Vector> VNMap;
+typedef Kernel::Vector_3                                       Vector;
+typedef boost::property_map<Mesh, CGAL::vertex_point_t>::type  VPMap;
+typedef Mesh::template Property_map<vertex_descriptor, Vector> VNMap;
 
 struct Bottom
 {
@@ -52,7 +52,7 @@ struct Top
 
 int main(int argc, char* argv[])
 {
-  SM in, out;
+  Mesh in, out;
 
   std::string filename = (argc > 1) ? std::string(argv[1]) : CGAL::data_file_path("meshes/cube-ouvert.off");
   double vlen = (argc > 2) ? std::stod(argv[2]) : 0.1;

Polygon_mesh_processing/examples/Polygon_mesh_processing/geodesic_isolevel_refinement.cpp

@@ -0,0 +1,72 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h>
+
+#include <CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h>
+#include <CGAL/Polygon_mesh_processing/connected_components.h>
+
+#include <fstream>
+#include <iostream>
+
+typedef CGAL::Simple_cartesian<double>                       Kernel;
+typedef Kernel::Point_3                                      Point_3;
+typedef CGAL::Surface_mesh<Point_3>                          Triangle_mesh;
+
+typedef boost::graph_traits<Triangle_mesh>::vertex_descriptor vertex_descriptor;
+typedef boost::graph_traits<Triangle_mesh>::edge_descriptor edge_descriptor;
+
+typedef Triangle_mesh::Property_map<vertex_descriptor,double> Vertex_distance_map;
+typedef CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<Triangle_mesh> Heat_method;
+
+int main(int argc, char* argv[])
+{
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off");
+
+  Triangle_mesh tm;
+  if(!CGAL::IO::read_polygon_mesh(filename, tm) ||
+     CGAL::is_empty(tm) || !CGAL::is_triangle_mesh(tm))
+  {
+    std::cerr << "Invalid input file." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  // default isovalues for cutting the mesh
+  std::vector<double> isovalues = {0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 10};
+
+  if (argc>2)
+  {
+    isovalues.clear();
+    for (int i=2; i<argc; ++i)
+      isovalues.push_back(atof(argv[i]));
+  }
+
+  //property map for the distance values to the source set
+  Vertex_distance_map vertex_distance = tm.add_property_map<vertex_descriptor, double>("v:distance", 0).first;
+
+  Heat_method hm(tm);
+
+  //use heat method to compute approximated geodesic distances to the source vertex `s`
+  vertex_descriptor s = *(vertices(tm).first);
+  hm.add_source(s);
+  hm.estimate_geodesic_distances(vertex_distance);
+
+  // property map to flag new cut edge added in the mesh
+  auto ecm = tm.add_property_map<edge_descriptor, bool>("e:is_constrained", 0).first;
+
+  // refine the mesh along isovalues
+  for (double isovalue : isovalues)
+    CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(tm, vertex_distance, isovalue, CGAL::parameters::edge_is_constrained_map(ecm));
+
+  // split the mesh in connected components bounded by the isocurves
+  std::vector<Triangle_mesh> edges_split;
+  CGAL::Polygon_mesh_processing::split_connected_components(tm, edges_split, CGAL::parameters::edge_is_constrained_map(ecm));
+
+  assert(argc!=1 || edges_split.size() == 22);
+
+  // export each submesh in a file
+  for(std::size_t i=0; i<edges_split.size(); ++i)
+    std::ofstream("out_"+std::to_string(i)+".off") << std::setprecision(17) << edges_split[i];
+
+  return 0;
+}

Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example.cpp

@@ -12,11 +12,11 @@
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel   Kernel;
 
-typedef CGAL::Polyhedron_3<Kernel>                            Polyhedron;
+typedef CGAL::Polyhedron_3<Kernel>                            Mesh;
 
-typedef Polyhedron::Vertex_handle                             Vertex_handle;
-typedef Polyhedron::Halfedge_handle                           Halfedge_handle;
-typedef Polyhedron::Facet_handle                              Facet_handle;
+typedef Mesh::Vertex_handle                                   Vertex_handle;
+typedef Mesh::Halfedge_handle                                 Halfedge_handle;
+typedef Mesh::Facet_handle                                    Facet_handle;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -24,7 +24,7 @@ int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mech-holes-shark.off");
 
-  Polyhedron poly;
+  Mesh poly;
   if(!PMP::IO::read_polygon_mesh(filename, poly))
   {
     std::cerr << "Invalid input." << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_LCC.cpp

@@ -14,11 +14,11 @@
 typedef CGAL::Exact_predicates_inexact_constructions_kernel                              Kernel;
 typedef Kernel::Point_3                                                                  Point;
 typedef CGAL::Linear_cell_complex_traits<3, Kernel>                                      MyTraits;
-typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper<2, 3, MyTraits>::type LCC;
+typedef CGAL::Linear_cell_complex_for_bgl_combinatorial_map_helper<2, 3, MyTraits>::type Mesh;
 
-typedef boost::graph_traits<LCC>::vertex_descriptor                                      vertex_descriptor;
-typedef boost::graph_traits<LCC>::halfedge_descriptor                                    halfedge_descriptor;
-typedef boost::graph_traits<LCC>::face_descriptor                                        face_descriptor;
+typedef boost::graph_traits<Mesh>::vertex_descriptor                                     vertex_descriptor;
+typedef boost::graph_traits<Mesh>::halfedge_descriptor                                   halfedge_descriptor;
+typedef boost::graph_traits<Mesh>::face_descriptor                                       face_descriptor;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -26,7 +26,7 @@ int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mech-holes-shark.off");
 
-  LCC mesh;
+  Mesh mesh;
   if(!PMP::IO::read_polygon_mesh(filename, mesh))
   {
     std::cerr << "Invalid input." << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_PH.cpp

@@ -7,17 +7,17 @@
 #include <boost/graph/graph_traits.hpp>
 
 #include <iostream>
-#include <unordered_map>
+#include <string>
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel;
-typedef CGAL::Polyhedron_3<Epic_kernel> Polyhedron;
-typedef boost::graph_traits<Polyhedron>::vertex_descriptor vertex_descriptor;
+typedef CGAL::Polyhedron_3<Epic_kernel>                     Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor        vertex_descriptor;
 
 int main(int argc, char* argv[])
 {
-  Polyhedron polyhedron;
+  Mesh polyhedron;
   const std::string filename = (argc > 1) ?
     argv[1] :
     CGAL::data_file_path("meshes/sphere.off");
@@ -29,11 +29,11 @@ int main(int argc, char* argv[])
   }
 
   // define property map to store curvature value and directions
-  boost::property_map<Polyhedron, CGAL::dynamic_vertex_property_t<Epic_kernel::FT>>::type
+  boost::property_map<Mesh, CGAL::dynamic_vertex_property_t<Epic_kernel::FT>>::type
     mean_curvature_map = get(CGAL::dynamic_vertex_property_t<Epic_kernel::FT>(), polyhedron),
     Gaussian_curvature_map = get(CGAL::dynamic_vertex_property_t<Epic_kernel::FT>(), polyhedron);
 
-  boost::property_map<Polyhedron, CGAL::dynamic_vertex_property_t<PMP::Principal_curvatures_and_directions<Epic_kernel>>>::type
+  boost::property_map<Mesh, CGAL::dynamic_vertex_property_t<PMP::Principal_curvatures_and_directions<Epic_kernel>>>::type
     principal_curvatures_and_directions_map =
     get(CGAL::dynamic_vertex_property_t<PMP::Principal_curvatures_and_directions<Epic_kernel>>(), polyhedron);
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_SM.cpp

@@ -11,12 +11,12 @@
 namespace PMP = CGAL::Polygon_mesh_processing;
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel;
-typedef CGAL::Surface_mesh<Epic_kernel::Point_3> Surface_Mesh;
-typedef boost::graph_traits<Surface_Mesh>::vertex_descriptor vertex_descriptor;
+typedef CGAL::Surface_mesh<Epic_kernel::Point_3> Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
 
 int main(int argc, char* argv[])
 {
-  Surface_Mesh smesh;
+  Mesh smesh;
   const std::string filename = (argc > 1) ?
     argv[1] :
     CGAL::data_file_path("meshes/sphere.off");
@@ -29,7 +29,7 @@ int main(int argc, char* argv[])
 
   // creating and tying surface mesh property maps for curvatures (with defaults = 0)
   bool created = false;
-  Surface_Mesh::Property_map<vertex_descriptor, Epic_kernel::FT>
+  Mesh::Property_map<vertex_descriptor, Epic_kernel::FT>
     mean_curvature_map, Gaussian_curvature_map;
 
   boost::tie(mean_curvature_map, created) =
@@ -41,7 +41,7 @@ int main(int argc, char* argv[])
   assert(created);
 
   // we use a tuple of 2 scalar values and 2 vectors for principal curvatures and directions
-  Surface_Mesh::Property_map<vertex_descriptor, PMP::Principal_curvatures_and_directions<Epic_kernel>>
+  Mesh::Property_map<vertex_descriptor, PMP::Principal_curvatures_and_directions<Epic_kernel>>
     principal_curvatures_and_directions_map;
 
   boost::tie(principal_curvatures_and_directions_map, created) =
@@ -67,4 +67,5 @@ int main(int argc, char* argv[])
       << ", GC = " << Gaussian_curvature_map[v] << "\n"
       << ", PC = [ " << PC.min_curvature << " , " << PC.max_curvature << " ]\n";
   }
+  return 0;
 }

Polygon_mesh_processing/examples/Polygon_mesh_processing/interpolated_corrected_curvatures_vertex.cpp

@@ -6,17 +6,18 @@
 #include <boost/graph/graph_traits.hpp>
 
 #include <iostream>
+#include <string>
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic_kernel;
-typedef CGAL::Surface_mesh<Epic_kernel::Point_3> Surface_Mesh;
-typedef boost::graph_traits<Surface_Mesh>::vertex_descriptor vertex_descriptor;
+typedef CGAL::Surface_mesh<Epic_kernel::Point_3>            Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor        vertex_descriptor;
 
 int main(int argc, char* argv[])
 {
   // instantiating and reading mesh
-  Surface_Mesh smesh;
+  Mesh smesh;
   const std::string filename = (argc > 1) ?
     argv[1] :
     CGAL::data_file_path("meshes/sphere.off");
@@ -47,4 +48,5 @@ int main(int argc, char* argv[])
       << ", GC = " << g << "\n"
       << ", PC = [ " << p.min_curvature << " , " << p.max_curvature << " ]\n";
   }
+  return 0;
 }

Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_with_custom_sizing_example.cpp

@@ -0,0 +1,101 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/remesh.h>
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Surface_mesh<K::Point_3> Mesh;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+// a sizing fied that is increasing the size of edge along the y-axis
+// starting at a minimum size at y-max and ending at a maximum size at
+// y-min, with a linear interpolation of sizes in between the two extreme
+// sizing values
+struct My_sizing_field
+{
+  double min_size, max_size;
+  double ymin, ymax;
+
+  My_sizing_field(double min_size, double max_size, double ymin, double ymax)
+    : min_size(min_size)
+    , max_size(max_size)
+    , ymin(ymin)
+    , ymax(ymax)
+  {}
+
+  double at(K::Point_3 p) const
+  {
+    double y=p.y();
+    return CGAL::square( (y-ymin)/(ymax-ymin) * (min_size - max_size) + max_size );
+  }
+  double at(const Mesh::Vertex_index v, const Mesh& mesh) const { return at(mesh.point(v)); }
+
+  std::optional<double> is_too_long(const Mesh::Vertex_index va,
+                                    const Mesh::Vertex_index vb,
+                                    const Mesh& mesh) const
+  {
+    // TODO: no mesh as parameters?
+    K::Point_3 mp = CGAL::midpoint(mesh.point(va), mesh.point(vb));
+    double sql_at = at(mp);
+    double sql = CGAL::squared_distance(mesh.point(va), mesh.point(vb));
+    if (sql > sql_at)
+      return sql / sql_at;
+    return std::nullopt;
+  }
+
+  std::optional<double> is_too_short(const Mesh::Halfedge_index h,
+                                     const Mesh& mesh) const
+  {
+    K::Point_3 mp = CGAL::midpoint(mesh.point(source(h, mesh)), mesh.point(target(h, mesh)));
+    double sql_at = at(mp);
+    double sql = CGAL::squared_distance(mesh.point(source(h, mesh)), mesh.point(target(h, mesh)));
+    if (sql < sql_at)
+      return sql / sql_at;
+    return std::nullopt;
+  }
+
+  K::Point_3 split_placement(const Mesh::Halfedge_index h,
+                             const Mesh& mesh) const
+  {
+    return CGAL::midpoint(mesh.point(source(h, mesh)), mesh.point(target(h, mesh)));
+  }
+
+  void register_split_vertex(const Mesh::Vertex_index, const Mesh&) {}
+};
+
+
+int main(int argc, char* argv[])
+{
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elk.off");
+
+  Mesh mesh;
+  if (!PMP::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) {
+    std::cerr << "Not a valid input file." << std::endl;
+    return 1;
+  }
+
+  std::cout << "Start remeshing of " << filename
+    << " (" << num_faces(mesh) << " faces)..." << std::endl;
+
+  CGAL::Bbox_3 bb = PMP::bbox(mesh);
+  My_sizing_field sizing_field(0.1, 30, bb.ymin(), bb.ymax());
+  unsigned int nb_iter = 5;
+
+  PMP::isotropic_remeshing(
+      faces(mesh),
+      sizing_field,
+      mesh,
+      CGAL::parameters::number_of_iterations(nb_iter)
+                       .number_of_relaxation_steps(3)
+      );
+
+  CGAL::IO::write_polygon_mesh("custom_remesh_out.off", mesh, CGAL::parameters::stream_precision(17));
+
+  std::cout << "Remeshing done." << std::endl;
+
+  return 0;
+}

Polygon_mesh_processing/examples/Polygon_mesh_processing/orient_polygon_soup_example.cpp

@@ -15,7 +15,7 @@
 #include <vector>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel          K;
-typedef CGAL::Polyhedron_3<K, CGAL::Polyhedron_items_with_id_3>      Polyhedron;
+typedef CGAL::Polyhedron_3<K, CGAL::Polyhedron_items_with_id_3>      Mesh;
 
 // Optional visitor for orientating a polygon soup to demonstrate usage for some functions.
 // inherits from the default class as some functions are not overloaded
@@ -59,12 +59,12 @@ int main(int argc, char* argv[])
   Visitor visitor;
   CGAL::Polygon_mesh_processing::orient_polygon_soup(points, polygons, CGAL::parameters::visitor(visitor));
 
-  Polyhedron mesh;
+  Mesh mesh;
   CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh);
 
   // Number the faces because 'orient_to_bound_a_volume' needs a face <--> index map
   int index = 0;
-  for(Polyhedron::Face_iterator fb=mesh.facets_begin(), fe=mesh.facets_end(); fb!=fe; ++fb)
+  for(Mesh::Face_iterator fb=mesh.facets_begin(), fe=mesh.facets_end(); fb!=fe; ++fb)
     fb->id() = index++;
 
   if(CGAL::is_closed(mesh))

Polygon_mesh_processing/examples/Polygon_mesh_processing/point_inside_example.cpp

@@ -12,14 +12,14 @@
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel    K;
 typedef K::Point_3                                             Point;
-typedef CGAL::Polyhedron_3<K>                                  Polyhedron;
+typedef CGAL::Polyhedron_3<K>                                  Mesh;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
-double max_coordinate(const Polyhedron& poly)
+double max_coordinate(const Mesh& poly)
 {
   double max_coord = -std::numeric_limits<double>::infinity();
-  for(Polyhedron::Vertex_handle v : vertices(poly))
+  for(Mesh::Vertex_handle v : vertices(poly))
   {
     Point p = v->point();
     max_coord = (std::max)(max_coord, CGAL::to_double(p.x()));
@@ -33,14 +33,14 @@ int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
 
-  Polyhedron poly;
+  Mesh poly;
   if(!PMP::IO::read_polygon_mesh(filename, poly) || CGAL::is_empty(poly) || !CGAL::is_triangle_mesh(poly))
   {
     std::cerr << "Invalid input." << std::endl;
     return 1;
   }
 
-  CGAL::Side_of_triangle_mesh<Polyhedron, K> inside(poly);
+  CGAL::Side_of_triangle_mesh<Mesh, K> inside(poly);
 
   double size = max_coordinate(poly);
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope.cpp

@@ -9,13 +9,13 @@ int main(int argc, char* argv[])
 {
   typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
   typedef Kernel::Point_3 Point_3;
-  typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
-  typedef boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
+  typedef CGAL::Surface_mesh<Point_3> Mesh;
+  typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
 
   typedef CGAL::Polyhedral_envelope<Kernel> Envelope;
 
   std::ifstream in((argc>1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"));
-  Surface_mesh tmesh;
+  Mesh tmesh;
 
   in >> tmesh;
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/polyhedral_envelope_mesh_containment.cpp

@@ -13,19 +13,19 @@ int main(int argc, char* argv[])
 {
   typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
   typedef Kernel::Point_3 Point_3;
-  typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
+  typedef CGAL::Surface_mesh<Point_3> Mesh;
 
   typedef CGAL::Polyhedral_envelope<Kernel> Envelope;
 
   std::ifstream in((argc>1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"));
-  Surface_mesh tmesh;
+  Mesh tmesh;
   in >> tmesh;
 
   // remesh the input using the longest edge size as target edge length
-  Surface_mesh query = tmesh;
-  Surface_mesh::Edge_iterator longest_edge_it =
+  Mesh query = tmesh;
+  Mesh::Edge_iterator longest_edge_it =
     std::max_element(edges(query).begin(), edges(query).end(),
-                     [&query](Surface_mesh::Edge_index e1, Surface_mesh::Edge_index e2)
+                     [&query](Mesh::Edge_index e1, Mesh::Edge_index e2)
                      {
                         return PMP::edge_length(halfedge(e1, query), query) <
                                PMP::edge_length(halfedge(e2, query), query);

Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp

@@ -10,9 +10,9 @@
 typedef CGAL::Exact_predicates_inexact_constructions_kernel   K;
 typedef K::Point_3                                            Point;
 
-typedef CGAL::Surface_mesh<Point>                             Surface_mesh;
-typedef boost::graph_traits<Surface_mesh>::vertex_descriptor  vertex_descriptor;
-typedef boost::graph_traits<Surface_mesh>::face_descriptor    face_descriptor;
+typedef CGAL::Surface_mesh<Point>                             Mesh;
+typedef boost::graph_traits<Mesh>::vertex_descriptor          vertex_descriptor;
+typedef boost::graph_traits<Mesh>::face_descriptor            face_descriptor;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -20,7 +20,7 @@ int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/eight.off");
 
-  Surface_mesh mesh;
+  Mesh mesh;
   if(!PMP::IO::read_polygon_mesh(filename, mesh))
   {
     std::cerr << "Invalid input." << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp

@@ -13,8 +13,8 @@
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 
-typedef CGAL::Polyhedron_3<Kernel>                          Polyhedron;
-typedef Polyhedron::Vertex_handle                           Vertex_handle;
+typedef CGAL::Polyhedron_3<Kernel>                          Mesh;
+typedef Mesh::Vertex_handle                                 Vertex_handle;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -34,7 +34,7 @@ void extract_k_ring(Vertex_handle v,
   {
     v = qv[current_index++];
 
-    Polyhedron::Halfedge_around_vertex_circulator e(v->vertex_begin()), e_end(e);
+    Mesh::Halfedge_around_vertex_circulator e(v->vertex_begin()), e_end(e);
     do {
       Vertex_handle new_v = e->opposite()->vertex();
       if (D.insert(std::make_pair(new_v, dist_v + 1)).second)
@@ -47,14 +47,14 @@ int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
 
-  Polyhedron poly;
+  Mesh poly;
   if(!PMP::IO::read_polygon_mesh(filename, poly) || !CGAL::is_triangle_mesh(poly))
   {
     std::cerr << "Invalid input." << std::endl;
     return 1;
   }
 
-  std::vector<Polyhedron::Facet_handle>  new_facets;
+  std::vector<Mesh::Facet_handle>  new_facets;
   std::vector<Vertex_handle> new_vertices;
 
   PMP::refine(poly, faces(poly),
@@ -68,7 +68,7 @@ int main(int argc, char* argv[])
   refined_off.close();
   std::cout << "Refinement added " << new_vertices.size() << " vertices." << std::endl;
 
-  Polyhedron::Vertex_iterator v = poly.vertices_begin();
+  Mesh::Vertex_iterator v = poly.vertices_begin();
   std::advance(v, 82/*e.g.*/);
   std::vector<Vertex_handle> region;
   extract_k_ring(v, 12/*e.g.*/, region);

Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_almost_planar_patches.cpp

@@ -14,13 +14,13 @@
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::Point_3 Point_3;
-typedef CGAL::Surface_mesh<Kernel::Point_3> Surface_mesh;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
 int main()
 {
-  Surface_mesh sm;
+  Mesh sm;
   CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/fandisk.off"), sm);
 
   //apply a perturbation to input vertices so that points are no longer coplanar
@@ -51,7 +51,7 @@ int main()
                                                      edge_is_constrained_map(CGAL::make_random_access_property_map(ecm)));
 
   // run the remeshing algorithm using filled properties
-  Surface_mesh out;
+  Mesh out;
   PMP::remesh_almost_planar_patches(sm,
                                     out,
                                     nb_regions, nb_corners,

Polygon_mesh_processing/examples/Polygon_mesh_processing/remesh_planar_patches.cpp

@@ -12,12 +12,12 @@
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::Point_3 Point_3;
-typedef CGAL::Surface_mesh<Kernel::Point_3> Surface_mesh;
+typedef CGAL::Surface_mesh<Kernel::Point_3> Mesh;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 int main()
 {
-  Surface_mesh sm;
+  Mesh sm;
   CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/cube_quad.off"), sm);
 
   // triangulate faces;
@@ -25,8 +25,8 @@ int main()
   std::cout << "Input mesh has " << faces(sm).size() << " faces" << std::endl;
   assert(faces(sm).size()==12);
 
-  Surface_mesh::Property_map<Surface_mesh::Edge_index, bool> ecm =
-    sm.add_property_map<Surface_mesh::Edge_index, bool>("ecm",false).first;
+  Mesh::Property_map<Mesh::Edge_index, bool> ecm =
+    sm.add_property_map<Mesh::Edge_index, bool>("ecm",false).first;
 
   // detect sharp edges of the cube
   PMP::detect_sharp_edges(sm, 60, ecm);
@@ -37,7 +37,7 @@ int main()
   assert(faces(sm).size()>100);
 
   // decimate the mesh
-  Surface_mesh out;
+  Mesh out;
   PMP::remesh_planar_patches(sm, out);
   CGAL::IO::write_polygon_mesh("cube_decimated.off", out, CGAL::parameters::stream_precision(17));
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp

@@ -9,7 +9,7 @@
 #include <string>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel   K;
-typedef CGAL::Polyhedron_3<K>                                 Polyhedron;
+typedef CGAL::Polyhedron_3<K>                                 Mesh;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -17,7 +17,7 @@ int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/quads_to_stitch.off");
 
-  Polyhedron mesh;
+  Mesh mesh;
   if(!PMP::IO::read_polygon_mesh(filename, mesh))
   {
     std::cerr << "Invalid input." << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example_OM.cpp

@@ -10,8 +10,8 @@
 #include <string>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-
 typedef OpenMesh::PolyMesh_ArrayKernelT< > Mesh;
+
 int main(int argc, char* argv[])
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/quads_to_stitch.off");

Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp

@@ -11,7 +11,7 @@
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef Kernel::Point_3                                     Point;
-typedef CGAL::Surface_mesh<Point>                           Surface_mesh;
+typedef CGAL::Surface_mesh<Point>                           Mesh;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 
@@ -20,7 +20,7 @@ int main(int argc, char* argv[])
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/P.off");
   const char* outfilename = (argc > 2) ? argv[2] : "P_tri.off";
 
-  Surface_mesh mesh;
+  Mesh mesh;
   if(!PMP::IO::read_polygon_mesh(filename, mesh))
   {
     std::cerr << "Error: Invalid input." << std::endl;
@@ -41,7 +41,7 @@ int main(int argc, char* argv[])
   PMP::triangulate_faces(mesh);
 
   // Confirm that all faces are triangles.
-  for(boost::graph_traits<Surface_mesh>::face_descriptor f : faces(mesh))
+  for(boost::graph_traits<Mesh>::face_descriptor f : faces(mesh))
   {
     if(!CGAL::is_triangle(halfedge(f, mesh), mesh))
       std::cerr << "Error: non-triangular face left in mesh." << std::endl;

Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp

@@ -6,15 +6,14 @@
 
 #include <iostream>
 #include <fstream>
-#include <map>
 #include <string>
 #include <unordered_map>
 #include <utility>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef Kernel::Point_3                    Point;
-typedef CGAL::Surface_mesh<Point>          Surface_mesh;
-typedef boost::graph_traits<Surface_mesh>::face_descriptor face_descriptor;
+typedef Kernel::Point_3                                     Point;
+typedef CGAL::Surface_mesh<Point>                           Mesh;
+typedef boost::graph_traits<Mesh>::face_descriptor          face_descriptor;
 
 class Insert_iterator
 {
@@ -41,7 +40,7 @@ public:
 };
 
 
-struct Visitor : public CGAL::Polygon_mesh_processing::Triangulate_faces::Default_visitor<Surface_mesh>
+struct Visitor : public CGAL::Polygon_mesh_processing::Triangulate_faces::Default_visitor<Mesh>
 {
    typedef std::unordered_map<face_descriptor,face_descriptor> Container;
 
@@ -78,7 +77,7 @@ int main(int argc, char* argv[])
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/P.off");
   std::ifstream input(filename);
 
-  Surface_mesh mesh;
+  Mesh mesh;
   if (!input || !(input >> mesh) || mesh.is_empty())
   {
     std::cerr << "Not a valid off file." << std::endl;
@@ -87,7 +86,7 @@ int main(int argc, char* argv[])
 
   std::unordered_map<face_descriptor,face_descriptor> t2q;
 
-  Surface_mesh copy;
+  Mesh copy;
 
   CGAL::copy_face_graph(mesh, copy, CGAL::parameters::face_to_face_output_iterator(Insert_iterator(t2q)));
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/volume_connected_components.cpp

@@ -12,7 +12,7 @@
 #include <string>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel   Kernel;
-typedef CGAL::Surface_mesh<Kernel::Point_3>                   Surface_mesh;
+typedef CGAL::Surface_mesh<Kernel::Point_3>                   Mesh;
 
 namespace PMP = CGAL::Polygon_mesh_processing;
 namespace params = CGAL::parameters;
@@ -21,7 +21,7 @@ int main(int argc, char** argv)
 {
   const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off");
 
-  Surface_mesh mesh;
+  Mesh mesh;
   if(!PMP::IO::read_polygon_mesh(filename, mesh))
   {
     std::cerr << "Invalid input." << std::endl;
@@ -29,8 +29,8 @@ int main(int argc, char** argv)
   }
 
   // property map to assign a volume id for each face
-  Surface_mesh::Property_map<Surface_mesh::Face_index, std::size_t> vol_id_map =
-    mesh.add_property_map<Surface_mesh::Face_index, std::size_t>().first;
+  Mesh::Property_map<Mesh::Face_index, std::size_t> vol_id_map =
+    mesh.add_property_map<Mesh::Face_index, std::size_t>().first;
 
   // fill the volume id map
   std::vector<PMP::Volume_error_code> err_codes;
@@ -39,7 +39,7 @@ int main(int argc, char** argv)
   std::cout << "Found " << nb_vol << " volumes\n";
 
   // write each volume in an OFF file
-  typedef CGAL::Face_filtered_graph<Surface_mesh> Filtered_graph;
+  typedef CGAL::Face_filtered_graph<Mesh> Filtered_graph;
   Filtered_graph vol_mesh(mesh, 0, vol_id_map);
   for(std::size_t id = 0; id < nb_vol; ++id)
   {
@@ -49,11 +49,12 @@ int main(int argc, char** argv)
     if(id > 0)
       vol_mesh.set_selected_faces(id, vol_id_map);
 
-    Surface_mesh out;
+    Mesh out;
     CGAL::copy_face_graph(vol_mesh, out);
     std::ostringstream oss;
     oss << "vol_" << id <<".off";
     std::ofstream os(oss.str().data());
     os << out;
   }
+  return 0;
 }

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Adaptive_sizing_field.h

@@ -38,6 +38,8 @@ namespace Polygon_mesh_processing
 * Edges too long with respect to the local target edge length are split in two, while
 * edges that are too short are collapsed.
 *
+* This class depends on the \eigen library.
+*
 * \cgalModels{PMPSizingField}
 *
 * \sa `isotropic_remeshing()`
@@ -213,13 +215,13 @@ private:
   }
 
 public:
-  FT at(const vertex_descriptor v) const
+  FT at(const vertex_descriptor v, const PolygonMesh& /* pmesh */) const
   {
     CGAL_assertion(get(m_vertex_sizing_map, v));
     return get(m_vertex_sizing_map, v);
   }
 
-  std::optional<FT> is_too_long(const vertex_descriptor va, const vertex_descriptor vb) const
+  std::optional<FT> is_too_long(const vertex_descriptor va, const vertex_descriptor vb, const PolygonMesh& /* pmesh */) const
   {
     const FT sqlen = sqlength(va, vb);
     FT sqtarg_len = CGAL::square(4./3. * (CGAL::min)(get(m_vertex_sizing_map, va),
@@ -251,7 +253,7 @@ public:
                     get(m_vpmap, source(h, pmesh)));
   }
 
-  void update(const vertex_descriptor v, const PolygonMesh& pmesh)
+  void register_split_vertex(const vertex_descriptor v, const PolygonMesh& pmesh)
   {
     // calculating it as the average of two vertices on other ends
     // of halfedges as updating is done during an edge split

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Uniform_sizing_field.h

@@ -102,12 +102,12 @@ private:
   }
 
 public:
-  FT at(const vertex_descriptor /* v */) const
+  FT at(const vertex_descriptor /* v */, const PolygonMesh& /* pmesh */) const
   {
     return m_size;
   }
 
-  std::optional<FT> is_too_long(const vertex_descriptor va, const vertex_descriptor vb) const
+  std::optional<FT> is_too_long(const vertex_descriptor va, const vertex_descriptor vb, const PolygonMesh& /* pmesh */) const
   {
     const FT sqlen = sqlength(va, vb);
     if (sqlen > m_sq_long)
@@ -133,7 +133,7 @@ public:
                     get(m_vpmap, source(h, pmesh)));
   }
 
-  void update(const vertex_descriptor /* v */, const PolygonMesh& /* pmesh */)
+  void register_split_vertex(const vertex_descriptor /* v */, const PolygonMesh& /* pmesh */)
   {}
 
 private:

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h

@@ -246,7 +246,7 @@ namespace internal {
             get(ecmap, e) ||
             get(fpm, face(h,pmesh))!=get(fpm, face(opposite(h,pmesh),pmesh)) )
       {
-        if (sizing.is_too_long(source(h, pmesh), target(h, pmesh)))
+        if (sizing.is_too_long(source(h, pmesh), target(h, pmesh), pmesh))
         {
           return false;
         }
@@ -400,7 +400,7 @@ namespace internal {
       for(edge_descriptor e : edge_range)
       {
         const halfedge_descriptor he = halfedge(e, mesh_);
-        std::optional<double> sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_));
+        std::optional<double> sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_), mesh_);
         if(sqlen != std::nullopt)
           long_edges.emplace(he, sqlen.value());
       }
@@ -433,16 +433,16 @@ namespace internal {
         std::cout << "   refinement point : " << refinement_point << std::endl;
 #endif
         //update sizing field with the new point
-        sizing.update(vnew, mesh_);
+        sizing.register_split_vertex(vnew, mesh_);
 
         //check sub-edges
         //if it was more than twice the "long" threshold, insert them
-        std::optional<double> sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_));
+        std::optional<double> sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_), mesh_);
         if(sqlen_new != std::nullopt)
           long_edges.emplace(hnew, sqlen_new.value());
 
         const halfedge_descriptor hnext = next(hnew, mesh_);
-        sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_));
+        sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_), mesh_);
         if (sqlen_new != std::nullopt)
           long_edges.emplace(hnext, sqlen_new.value());
 
@@ -500,7 +500,7 @@ namespace internal {
         if (!is_split_allowed(e))
           continue;
         const halfedge_descriptor he = halfedge(e, mesh_);
-        std::optional<double> sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_));
+        std::optional<double> sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_), mesh_);
         if(sqlen != std::nullopt)
           long_edges.emplace(halfedge(e, mesh_), sqlen.value());
       }
@@ -550,16 +550,16 @@ namespace internal {
         halfedge_added(hnew_opp, status(opposite(he, mesh_)));
 
         //update sizing field with the new point
-        sizing.update(vnew, mesh_);
+        sizing.register_split_vertex(vnew, mesh_);
 
         //check sub-edges
         //if it was more than twice the "long" threshold, insert them
-        std::optional<double> sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_));
+        std::optional<double> sqlen_new = sizing.is_too_long(source(hnew, mesh_), target(hnew, mesh_), mesh_);
         if(sqlen_new != std::nullopt)
           long_edges.emplace(hnew, sqlen_new.value());
 
         const halfedge_descriptor hnext = next(hnew, mesh_);
-        sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_));
+        sqlen_new = sizing.is_too_long(source(hnext, mesh_), target(hnext, mesh_), mesh_);
         if (sqlen_new != std::nullopt)
           long_edges.emplace(hnext, sqlen_new.value());
 
@@ -580,7 +580,7 @@ namespace internal {
 
           if (snew == PATCH)
           {
-            std::optional<double> sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_));
+            std::optional<double> sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_), mesh_);
             if(sql != std::nullopt)
               long_edges.emplace(hnew2, sql.value());
           }
@@ -603,7 +603,7 @@ namespace internal {
 
           if (snew == PATCH)
           {
-            std::optional<double> sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_));
+            std::optional<double> sql = sizing.is_too_long(source(hnew2, mesh_), target(hnew2, mesh_), mesh_);
             if (sql != std::nullopt)
               long_edges.emplace(hnew2, sql.value());
           }
@@ -747,7 +747,7 @@ namespace internal {
         for(halfedge_descriptor ha : halfedges_around_target(va, mesh_))
         {
           vertex_descriptor va_i = source(ha, mesh_);
-          std::optional<double> sqha = sizing.is_too_long(vb, va_i);
+          std::optional<double> sqha = sizing.is_too_long(vb, va_i, mesh_);
           if (sqha != std::nullopt)
           {
             collapse_ok = false;

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Sizing_field_base.h

@@ -60,13 +60,14 @@ public:
   typedef typename K::FT                                        FT;
 
 public:
-  virtual FT at(const vertex_descriptor v) const = 0;
+  virtual FT at(const vertex_descriptor v, const PolygonMesh&) const = 0;
   virtual std::optional<FT> is_too_long(const vertex_descriptor va,
-                                        const vertex_descriptor vb) const = 0;
+                                        const vertex_descriptor vb,
+                                        const PolygonMesh&) const = 0;
   virtual std::optional<FT> is_too_short(const halfedge_descriptor h,
                                          const PolygonMesh& pmesh) const = 0;
   virtual Point_3 split_placement(const halfedge_descriptor h, const PolygonMesh& pmesh) const = 0;
-  virtual void update(const vertex_descriptor v, const PolygonMesh& pmesh) = 0;
+  virtual void register_split_vertex(const vertex_descriptor v, const PolygonMesh& pmesh) = 0;
 
 };
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h

@@ -0,0 +1,174 @@
+// Copyright (c) 2021-2023 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Sebastien Loriot
+
+#ifndef CGAL_POLYGON_MESH_PROCESSING_REFINE_MESH_AT_ISOLEVEL_H
+#define CGAL_POLYGON_MESH_PROCESSING_REFINE_MESH_AT_ISOLEVEL_H
+
+#include <CGAL/license/Polygon_mesh_processing/miscellaneous.h>
+
+#include <CGAL/boost/graph/Euler_operations.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <unordered_set>
+#include <unordered_map>
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+
+/*!
+ * \ingroup PkgPolygonMeshProcessingRef
+ *
+ * refines `pm` by adding new vertices on edges having their incident vertices associated with
+ * values respectively larger and smaller than `isovalue` in `value_map`.
+ * The placement of new vertices on edges will be done by linear interpolation
+ * using the aforementioned values.
+ * New vertices will be associated `isovalue` in `value_map` when created.
+ * Additionally, new edges will be added by connecting new vertices created sharing
+ * a common incident face. Note that in case more than two new vertices are added
+ * on a face boundary, no edges will be created in that face.
+ *
+ * @tparam PolygonMesh a model of the concepts `EdgeListGraph` and `FaceListGraph`
+ * @tparam ValueMap a model of the concept `ReadWritePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+ *                  as key type and with its value type being the type of the coordinates of points associated with vertices
+ *                  in the vertex map provided to the `vertex_point_map()` named parameter.
+ * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" for `pm`
+ *
+ * @param pm the polygon mesh to be refined.
+ * @param value_map the property map containing a value at each vertex for a given function defined over the mesh.
+ * @param isovalue the value used to refine
+ * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{edge_is_constrained_map}
+ *     \cgalParamDescription{an output property map associating `true` to all edges connecting vertices on the isolevel,
+ *                           and `false` for all other edges.}
+ *     \cgalParamType{a class model of `WritablePropertyMap` with `boost::graph_traits<PolygonMesh>::%edge_descriptor`
+ *                    as key type and `bool` as value type}
+ *     \cgalParamDefault{No marks on edges will be put}
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{vertex_point_map}
+ *     \cgalParamDescription{a property map associating points to the vertices of `pm`}
+ *     \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+ *                    as key type and `%Point_3` as value type}
+ *     \cgalParamDefault{`boost::get(CGAL::vertex_point, pm)`}
+ *     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+ *                     must be available in `PolygonMesh`.}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ */
+template <class PolygonMesh, class ValueMap, class NamedParameters = parameters::Default_named_parameters>
+void refine_mesh_at_isolevel(PolygonMesh& pm,
+                             ValueMap value_map,
+                             typename boost::property_traits<ValueMap>::value_type isovalue,
+                             const NamedParameters& np = parameters::default_values())
+{
+  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<PolygonMesh>::edge_descriptor edge_descriptor;
+  typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
+  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
+  typedef typename boost::property_traits<ValueMap>::value_type FT;
+
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+  using parameters::is_default_parameter;
+
+  typedef Static_boolean_property_map<edge_descriptor, false>  Default_ECM;
+  typedef typename internal_np::Lookup_named_param_def<internal_np::edge_is_constrained_t,
+                                                       NamedParameters,
+                                                       Default_ECM>::type       ECM;
+  typedef typename GetVertexPointMap < PolygonMesh, NamedParameters>::type VPM;
+
+  VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                             get_property_map(vertex_point, pm));
+
+  ECM ecm = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), Default_ECM());
+
+  std::unordered_map<face_descriptor, std::vector<halfedge_descriptor> > faces_to_split;
+  std::vector<edge_descriptor> to_split;
+  std::unordered_set<vertex_descriptor> vertices_on_isoline;
+  for (edge_descriptor e : edges(pm))
+  {
+    vertex_descriptor src = source(e, pm), tgt = target(e, pm);
+
+    if (get(value_map, src)==isovalue)
+    {
+      vertices_on_isoline.insert(src);
+      if (get(value_map, tgt)==isovalue)
+      {
+        put(ecm, e, true); // special case for faces entirely on an isovalue
+        continue;
+      }
+      continue;
+    }
+    if (get(value_map, tgt)==isovalue)
+    {
+      vertices_on_isoline.insert(tgt);
+      continue;
+    }
+    if ( (get(value_map, tgt) < isovalue) != (get(value_map, src) < isovalue) )
+    {
+      to_split.push_back(e);
+    }
+  }
+
+  for (edge_descriptor e : to_split)
+  {
+    vertex_descriptor src = source(e, pm), tgt = target(e, pm);
+    FT ds = get(value_map, src);
+    FT dt = get(value_map, tgt);
+    FT alpha = (isovalue - dt) / (ds - dt);
+    halfedge_descriptor hnew = CGAL::Euler::split_edge(halfedge(e, pm), pm);
+    put(vpm, target(hnew, pm), barycenter(get(vpm,src), alpha, get(vpm, tgt), 1-alpha));
+    put(value_map,  target(hnew, pm) , isovalue);
+    face_descriptor f = face(hnew, pm);
+    if (f!=boost::graph_traits<PolygonMesh>::null_face())
+      faces_to_split[f].push_back(hnew);
+    hnew=pm.prev(opposite(hnew, pm));
+    f = face(hnew, pm);
+    if (f!=boost::graph_traits<PolygonMesh>::null_face())
+      faces_to_split[f].push_back(hnew);
+  }
+
+  for (vertex_descriptor vh : vertices_on_isoline)
+  {
+    for (halfedge_descriptor h : halfedges_around_target(vh, pm))
+    {
+      face_descriptor f = face(h, pm);
+      if (f!=boost::graph_traits<PolygonMesh>::null_face())
+        faces_to_split[f].push_back(h);
+    }
+  }
+
+  for (const auto& p : faces_to_split)
+  {
+    if(p.second.size()!=2) continue;
+
+    std::pair<edge_descriptor ,bool> res = edge(target(p.second[0],pm),
+                                                target(p.second[1],pm), pm);
+    if (res.second)
+    {
+      // no split as the edge already exists (the two vertices are on the isolevel)
+      put(ecm, res.first, true);
+      continue;
+    }
+
+    halfedge_descriptor hnew = CGAL::Euler::split_face(p.second[0], p.second[1], pm);
+    put(ecm, edge(hnew, pm), true);
+  }
+}
+
+} } // end of CGAL::Polygon_mesh_processing
+
+
+#endif

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh_planar_patches.h

@@ -299,7 +299,7 @@ template <typename Kernel,
           typename VertexPointMap>
 void
 mark_constrained_edges(
-  TriangleMesh& tm,
+  const TriangleMesh& tm,
   EdgeIsConstrainedMap edge_is_constrained,
   double coplanar_cos_threshold,
   const VertexPointMap& vpm)
@@ -319,7 +319,7 @@ template <typename Kernel,
           typename VertexCornerIdMap>
 std::size_t
 mark_corner_vertices(
-  TriangleMesh& tm,
+  const TriangleMesh& tm,
   EdgeIsConstrainedMap& edge_is_constrained,
   VertexCornerIdMap& vertex_corner_id,
   double coplanar_cos_threshold,
@@ -546,7 +546,7 @@ template <typename Kernel,
           typename FaceCCIdMap,
           typename VertexPointMap>
 std::pair<std::size_t, std::size_t>
-tag_corners_and_constrained_edges(TriangleMesh& tm,
+tag_corners_and_constrained_edges(const TriangleMesh& tm,
                                   double coplanar_cos_threshold,
                                   VertexCornerIdMap& vertex_corner_id,
                                   EdgeIsConstrainedMap& edge_is_constrained,

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h

@@ -288,9 +288,9 @@ void tangential_relaxation(const VertexRange& vertices,
 
             const double tri_area = gt_area(get(vpm, v), get(vpm, v1), get(vpm, v2));
             const double face_weight = tri_area
-                                       / (1. / 3. * (sizing.at(v)
-                                                   + sizing.at(v1)
-                                                   + sizing.at(v2)));
+                                       / (1. / 3. * (sizing.at(v, tm)
+                                                   + sizing.at(v1, tm)
+                                                   + sizing.at(v2, tm)));
             weight += face_weight;
 
             const Point_3 centroid = gt_centroid(get(vpm, v), get(vpm, v1), get(vpm, v2));

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -65,6 +65,7 @@ create_single_source_cgal_program("triangulate_hole_with_cdt_2_test.cpp")
 create_single_source_cgal_program("test_pmp_polyhedral_envelope.cpp")
 create_single_source_cgal_program("test_pmp_np_function.cpp")
 create_single_source_cgal_program("test_degenerate_pmp_clip_split_corefine.cpp")
+create_single_source_cgal_program("test_isolevel_refinement.cpp")
 #  create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp")
 
 find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater)

Polygon_mesh_processing/test/Polygon_mesh_processing/test_isolevel_refinement.cpp

@@ -0,0 +1,68 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h>
+#include <CGAL/Polygon_mesh_processing/connected_components.h>
+
+#include <fstream>
+#include <iostream>
+
+typedef CGAL::Simple_cartesian<double>                       Kernel;
+typedef Kernel::Point_3                                      Point_3;
+typedef CGAL::Surface_mesh<Point_3>                          Triangle_mesh;
+
+typedef boost::graph_traits<Triangle_mesh>::vertex_descriptor vertex_descriptor;
+typedef boost::graph_traits<Triangle_mesh>::edge_descriptor edge_descriptor;
+
+typedef Triangle_mesh::Property_map<vertex_descriptor,double> Vertex_distance_map;
+
+int main()
+{
+  const std::string filename = CGAL::data_file_path("meshes/elephant.off");
+
+  Triangle_mesh tm;
+  if(!CGAL::IO::read_polygon_mesh(filename, tm) ||
+     CGAL::is_empty(tm) || !CGAL::is_triangle_mesh(tm))
+  {
+    std::cerr << "Invalid input file." << std::endl;
+    return EXIT_FAILURE;
+  }
+  //property map for the distance values to the source set
+  Vertex_distance_map vertex_distance = tm.add_property_map<vertex_descriptor, double>("v:distance", 5).first;
+  std::vector<int> zero_vids={
+    /*a closed polyline*/ 2144,145,2690,1752,339,215,1395,338,77,2145,2052,2054,343,1936,22,1751,214,1499,142,358,2694,1750,301,65,59,2650,2060,205,2651,2061,2490,1939,898,13,298,
+    /*two adjacent triangles*/ 532, 185, 534, 2735,
+    /*another patch with missing crossed edges*/134,73,1883,2533,72,532,185,131,534
+  };
+
+  std::vector<int> minus_vids = {132, 364};
+
+  for (int i : zero_vids)
+    put(vertex_distance, Triangle_mesh::Vertex_index(i), 0);
+  for (int i : minus_vids)
+    put(vertex_distance, Triangle_mesh::Vertex_index(i), -5);
+
+  // property map to flag new cut edge added in the mesh
+  auto ecm = tm.add_property_map<edge_descriptor, bool>("e:is_constrained", 0).first;
+
+  CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(tm, vertex_distance, 0, CGAL::parameters::edge_is_constrained_map(ecm));
+
+  std::ofstream debug("edges.polylines.txt");
+  for (Triangle_mesh::Edge_index e : edges(tm))
+    if (get(ecm, e))
+      debug << "2 " << tm.point(source(e, tm)) << " " <<  tm.point(target(e, tm)) << "\n";
+  debug.close();
+
+  // split the mesh in connected components bounded by the isocurves
+  std::vector<Triangle_mesh> edges_split;
+  CGAL::Polygon_mesh_processing::split_connected_components(tm, edges_split, CGAL::parameters::edge_is_constrained_map(ecm));
+
+  // export each submesh in a file
+  for(std::size_t i=0; i<edges_split.size(); ++i)
+    std::ofstream("out_"+std::to_string(i)+".off") << std::setprecision(17) << edges_split[i];
+
+  assert(edges_split.size() == 5);
+
+
+  return 0;
+}

Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp

@@ -335,7 +335,7 @@ void test_triangulate_refine_and_fair_hole_compile() {
   std::vector<Vertex_handle> patch_vertices;
 
   // use all param
-  read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
+  read_poly_with_borders("elephant_quad_hole_no_DT3.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
   (poly, border_reps[0],
   CGAL::parameters::
@@ -345,7 +345,7 @@ void test_triangulate_refine_and_fair_hole_compile() {
     sparse_linear_solver(Default_solver()));
 
   // default solver
-  read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
+  read_poly_with_borders("elephant_quad_hole_no_DT3.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
     (poly, border_reps[0],
     CGAL::parameters::
@@ -354,7 +354,7 @@ void test_triangulate_refine_and_fair_hole_compile() {
       weight_calculator(CGAL::Weights::Uniform_weight<Polyhedron>()));
 
   // default solver and weight
-  read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
+  read_poly_with_borders("elephant_quad_hole_no_DT3.off", poly, border_reps);
   CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole
     (poly, border_reps[0],
         CGAL::parameters::
@@ -376,11 +376,11 @@ void generate_elephant_with_hole()
     {
       Halfedge_handle nh=opposite(halfedge(fd,poly), poly);
       CGAL::Euler::remove_face(halfedge(fd, poly), poly);
-      std::ofstream output("elephant_triangle_hole.off");
+      std::ofstream output("elephant_triangle_hole_no_DT3.off");
       output << poly;
       output.close();
       CGAL::Euler::remove_face(nh, poly);
-      output.open("elephant_quad_hole.off");
+      output.open("elephant_quad_hole_no_DT3.off");
       output << poly;
       return;
     }
@@ -396,8 +396,8 @@ typedef CGAL::Surface_mesh<typename Kernel::Point_3> Polyhedron;
   generate_elephant_with_hole<Polyhedron>();
 
   std::vector<std::string> input_files;
-  input_files.push_back("elephant_triangle_hole.off");
-  input_files.push_back("elephant_quad_hole.off");
+  input_files.push_back("elephant_triangle_hole_no_DT3.off");
+  input_files.push_back("elephant_quad_hole_no_DT3.off");
   input_files.push_back(CGAL::data_file_path("meshes/mech-holes-shark.off"));
   // std::cerr.precision(15);
   for(std::vector<std::string>::iterator it = input_files.begin(); it != input_files.end(); ++it) {

Polyhedron/demo/Polyhedron/CGAL_double_edit.cpp

@@ -19,7 +19,6 @@ public:
   }
   QValidator::State validate ( QString & input, int & pos ) const
   {
-    fixup(input);
     return QDoubleValidator::validate(input, pos);
   }
 };
@@ -59,7 +58,7 @@ public:
 
   void DoubleEdit::setRange(double rmin, double rmax)
   {
-    this->validator->setRange(rmin, rmax, this->validator->decimals());
+    this->validator->setRange(rmin, rmax, -1);
   }
 
   double DoubleEdit::getValue()

Polyhedron/demo/Polyhedron/Plugins/Alpha_wrap_3/Alpha_wrap_3_plugin.cpp

@@ -9,6 +9,10 @@
 #include "Scene_polylines_item.h"
 #include "Scene_points_with_normal_item.h"
 
+// Since we want to do visualization and interruption, it's better to use the sorted priority queue,
+// even if it is slower
+#define CGAL_AW3_USE_SORTED_PRIORITY_QUEUE
+
 #include <CGAL/alpha_wrap_3.h>
 
 #include <QAction>
@@ -35,13 +39,13 @@
 using TS_Oracle = CGAL::Alpha_wraps_3::internal::Triangle_soup_oracle<Kernel>;
 using SS_Oracle = CGAL::Alpha_wraps_3::internal::Segment_soup_oracle<Kernel, TS_Oracle>;
 using Oracle = CGAL::Alpha_wraps_3::internal::Point_set_oracle<Kernel, SS_Oracle>;
-using Wrapper = CGAL::Alpha_wraps_3::internal::Alpha_wrap_3<Oracle>;
+using Wrapper = CGAL::Alpha_wraps_3::internal::Alpha_wrapper_3<Oracle>;
 
 // Here is the pipeline for the interruption box:
 // - The main window is connected to a wrapping thread, which performs the wrapping.
 // - The wrapping has a visitor, AW3_interrupter_visitor, which has a shared_ptr to a Boolean
-// - When the user clicks the box, the Boolean is switched to *false*, and the visitor throws
-// - The wrapping thread catches the exception, and creates the wip mesh
+// - When the user clicks the box, the Boolean is switched to *false*
+// - The wrapping thread creates the wip mesh
 
 // Here is the pipeline for the iterative visualization:
 // - The main window is connected to a wrapping thread, which performs the wrapping.
@@ -113,10 +117,10 @@ public:
     if(!points || !faces || !fcolors || !vcolors)
       return;
 
-    // If the next top of the queue has vertices on the bbox, don't draw (as to avoid producing
+    // If the next top of the queue has vertices on the bbox, don't draw (try to avoid producing
     // spikes in the visualization)
 //    const auto& gate = wrapper.queue().top();
-//    if(wrapper.triangulation().number_of_vertices() > 500 && gate.is_artificial_facet())
+//    if(wrapper.triangulation().number_of_vertices() > 500 && gate.is_permissive_facet())
 //      return;
 
     // Skip some...
@@ -201,9 +205,6 @@ public:
   }
 };
 
-// Use a throw to get out of the AW3 refinement loop
-class Out_of_patience_exception : public std::exception { };
-
 template <typename BaseVisitor>
 struct AW3_interrupter_visitor
   : BaseVisitor
@@ -215,15 +216,10 @@ struct AW3_interrupter_visitor
     : BaseVisitor(base)
   { }
 
-  // Only overload this one because it gives a better state of the wrap (for other visitor calls,
-  // we often get tetrahedral spikes because there are artificial gates in the queue)
-  template <typename Wrapper, typename Point>
-  void before_Steiner_point_insertion(const Wrapper& wrapper, const Point& p)
+  template <typename Wrapper>
+  constexpr bool go_further(const Wrapper&)
   {
-    if(*should_stop)
-      throw Out_of_patience_exception();
-
-    return BaseVisitor::before_Steiner_point_insertion(wrapper, p);
+    return !(*should_stop);
   }
 };
 
@@ -273,25 +269,14 @@ public:
     QElapsedTimer elapsed_timer;
     elapsed_timer.start();
 
-    // try-catch because the stop visitor currently uses a throw
-    try
-    {
-      wrapper(alpha, offset, wrap,
-              CGAL::parameters::do_enforce_manifoldness(enforce_manifoldness)
-                               .visitor(visitor));
+    wrapper(alpha, offset, wrap,
+            CGAL::parameters::do_enforce_manifoldness(enforce_manifoldness)
+                             .visitor(visitor));
 
+    if(wrapper.queue().empty())
       Q_EMIT done(this);
-    }
-    catch(const Out_of_patience_exception&)
-    {
-      if(enforce_manifoldness)
-        wrapper.make_manifold();
-
-      // extract the wrap in its current state
-      wrapper.extract_surface(wrap, CGAL::get(CGAL::vertex_point, wrap), !enforce_manifoldness);
-
+    else
       Q_EMIT interrupted(this);
-    }
 
     std::cout << "Wrapping took " << elapsed_timer.elapsed() / 1000. << "s" << std::endl;
   }
@@ -562,8 +547,6 @@ public Q_SLOTS:
           triangles.emplace_back(get(vpm, target(h, *pMesh)),
                                  get(vpm, target(next(h, *pMesh), *pMesh)),
                                  get(vpm, source(h, *pMesh)));
-
-          m_wrap_bbox += triangles.back().bbox();
         }
 
         continue;
@@ -586,8 +569,6 @@ public Q_SLOTS:
           triangles.emplace_back(soup_item->points()[p[0]],
                                  soup_item->points()[p[1]],
                                  soup_item->points()[p[2]]);
-
-          m_wrap_bbox += triangles.back().bbox();
         }
 
         continue;
@@ -614,8 +595,6 @@ public Q_SLOTS:
           triangles.emplace_back(get(vpm, target(h, *pMesh)),
                                  get(vpm, target(next(h, *pMesh), *pMesh)),
                                  get(vpm, source(h, *pMesh)));
-
-          m_wrap_bbox += triangles.back().bbox();
         }
 
         segments.reserve(segments.size() + selection_item->selected_edges.size());
@@ -623,16 +602,12 @@ public Q_SLOTS:
         {
           segments.emplace_back(get(vpm, target(halfedge(e, *pMesh), *pMesh)),
                                 get(vpm, target(opposite(halfedge(e, *pMesh), *pMesh), *pMesh)));
-
-          m_wrap_bbox += segments.back().bbox();
         }
 
         points.reserve(points.size() + selection_item->selected_vertices.size());
         for(const auto& v : selection_item->selected_vertices)
         {
           points.push_back(get(vpm, v));
-
-          m_wrap_bbox += points.back().bbox();
         }
 
         continue;
@@ -671,6 +646,22 @@ public Q_SLOTS:
     std::cout << segments.size() << " edges" << std::endl;
     std::cout << points.size() << " points" << std::endl;
 
+    if(!triangles.empty())
+      m_wrap_bbox = triangles.front().bbox();
+    else if(!segments.empty())
+      m_wrap_bbox = segments.front().bbox();
+    else if(!points.empty())
+      m_wrap_bbox = points.front().bbox();
+
+    for(const Kernel::Triangle_3& tr : triangles)
+      m_wrap_bbox += tr.bbox();
+    for(const Kernel::Segment_3& sg : segments)
+      m_wrap_bbox += sg.bbox();
+    for(const Kernel::Point_3& pt : points)
+      m_wrap_bbox += pt.bbox();
+
+    std::cout << "Bbox:\n" << m_wrap_bbox << std::endl;
+
     // The relative value uses the bbox of the full scene and not that of selected items to wrap
     // This is intentional, both because it's tedious to make it otherwise, and because it seems
     // to be simpler to compare between "all wrapped" / "some wrapped"
@@ -703,8 +694,8 @@ public Q_SLOTS:
 
     const bool use_message_box = ui.enableMessageBox->isChecked();
 
-    std::cout << "Wrapping edges? " << std::boolalpha << wrap_segments << std::endl;
     std::cout << "Wrapping faces? " << std::boolalpha << wrap_triangles << std::endl;
+    std::cout << "Wrapping edges? " << std::boolalpha << wrap_segments << std::endl;
 
     if(!wrap_triangles)
     {
@@ -756,15 +747,11 @@ public Q_SLOTS:
 
     if(alpha <= 0. || offset <= 0.)
     {
-      print_message("Warning: alpha/offset must be strictly positive - nothing to wrap");
+      print_message("Warning: alpha/offset must be strictly positive");
       QApplication::restoreOverrideCursor();
       return;
     }
 
-    // Switch from 'wait' to 'busy'
-    QApplication::restoreOverrideCursor();
-    QApplication::setOverrideCursor(Qt::BusyCursor);
-
     for(int index : this->scene->selectionIndices())
     {
       Scene_surface_mesh_item* sm_item = qobject_cast<Scene_surface_mesh_item*>(this->scene->item(index));
@@ -824,6 +811,10 @@ public Q_SLOTS:
     // Create message box with stop button
     if(use_message_box)
     {
+      // Switch from 'wait' to 'busy'
+      QApplication::restoreOverrideCursor();
+      QApplication::setOverrideCursor(Qt::BusyCursor);
+
       m_message_box = new QMessageBox(QMessageBox::NoIcon,
                                      "Wrapping",
                                      "Wrapping in progress...",
@@ -841,6 +832,8 @@ public Q_SLOTS:
     }
 
     // Actual start
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+
     wrapper_thread->start();
 
     CGAL::Three::Three::getMutex()->lock();

Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin.cpp

@@ -624,6 +624,8 @@ void Mesh_3_plugin::mesh_3(const Mesh_type mesh_type,
   ui.approx->setRange(diag * 10e-7, // min
                       diag);        // max
   ui.approx->setValue(approx);
+  ui.approx->setToolTip(tr("Approximation error: in [%1; %2]")
+                       .arg(diag * 10e-7).arg(diag));
 
   ui.protect->setEnabled(features_protection_available);
   ui.protect->setChecked(features_protection_available);

Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_3_plugin_cgal_code.cpp

@@ -134,6 +134,8 @@ Meshing_thread* cgal_code_mesh_3(QList<const SMesh*> pMeshes,
   param.manifold = manifold;
   param.protect_features = protect_features || protect_borders;
   param.use_sizing_field_with_aabb_tree = polylines.empty() && protect_features;
+  param.image_3_ptr = nullptr;
+  param.weights_ptr = nullptr;
 
   typedef ::Mesh_function<Polyhedral_mesh_domain,
                           Mesh_fnt::Polyhedral_domain_tag> Mesh_function;
@@ -237,6 +239,8 @@ Meshing_thread* cgal_code_mesh_3(const QList<const SMesh*> pMeshes,
   param.manifold = manifold;
   param.protect_features = protect_features || protect_borders;
   param.use_sizing_field_with_aabb_tree = protect_features;
+  param.image_3_ptr = nullptr;
+  param.weights_ptr = nullptr;
 
   typedef ::Mesh_function<Polyhedral_complex_mesh_domain,
                           Mesh_fnt::Polyhedral_domain_tag> Mesh_function;
@@ -292,7 +296,8 @@ Meshing_thread* cgal_code_mesh_3(const Implicit_function_interface* pfunction,
   param.manifold = manifold;
   param.detect_connected_components = false; // to avoid random values
                                              // in the debug displays
-
+  param.image_3_ptr = nullptr;
+  param.weights_ptr = nullptr;
 
   typedef ::Mesh_function<Function_mesh_domain,
                           Mesh_fnt::Implicit_domain_tag> Mesh_function;

Polyhedron/demo/Polyhedron/Plugins/Mesh_3/Mesh_function.h

@@ -138,20 +138,40 @@ QStringList
 Mesh_parameters::
 log() const
 {
-  return QStringList()
-  << QString("edge max size: %1").arg(edge_sizing)
-  << QString("edge min size: %1").arg(edge_min_sizing)
-  << QString("facet min angle: %1").arg(facet_angle)
-  << QString("facet max size: %1").arg(facet_sizing)
-  << QString("facet min size: %1").arg(facet_min_sizing)
-  << QString("facet approx error: %1").arg(facet_approx)
-  << QString("tet shape (radius-edge): %1").arg(tet_shape)
-  << QString("tet max size: %1").arg(tet_sizing)
-  << QString("tet min size: %1").arg(tet_min_sizing)
-  << QString("detect connected components: %1")
-    .arg(detect_connected_components)
-  << QString("use weights: %1").arg(weights_ptr != nullptr)
-  << QString("protect features: %1").arg(protect_features);
+  QStringList res("Mesh criteria");
+
+  // doubles
+  if(edge_sizing > 0)
+    res << QString("edge max size: %1").arg(edge_sizing);
+  if(edge_min_sizing > 0)
+    res << QString("edge min size: %1").arg(edge_min_sizing);
+  if(facet_angle > 0)
+    res << QString("facet min angle: %1").arg(facet_angle);
+  if(facet_sizing > 0)
+    res << QString("facet max size: %1").arg(facet_sizing);
+  if(facet_min_sizing > 0)
+    res << QString("facet min size: %1").arg(facet_min_sizing);
+  if(facet_approx > 0)
+    res << QString("facet approx error: %1").arg(facet_approx);
+  if(tet_shape > 0)
+    res << QString("tet shape (radius-edge): %1").arg(tet_shape);
+  if(tet_sizing > 0)
+    res << QString("tet max size: %1").arg(tet_sizing);
+  if(tet_min_sizing > 0)
+    res << QString("tet min size: %1").arg(tet_min_sizing);
+
+  // booleans
+  res << QString("protect features: %1").arg(protect_features);
+  if(image_3_ptr != nullptr)
+  {
+    res << QString("detect connected components: %1")
+             .arg(detect_connected_components);
+    res << QString("use weights: %1").arg(weights_ptr != nullptr);
+  }
+  res << QString("use aabb tree: %1").arg(use_sizing_field_with_aabb_tree);
+  res << QString("manifold: %1").arg(manifold);
+
+  return res;
 }
 
 

STL_Extension/include/CGAL/STL_Extension/internal/mesh_parameters_interface.h

@@ -55,7 +55,7 @@ perturb(const CGAL_NP_CLASS& np = parameters::default_values())
 }
 
 template<typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_1, typename CGAL_NP_TEMPLATE_PARAMETERS_NO_DEFAULT_2, typename ... NP>
-Named_function_parameters<::CGAL::parameters::internal::Exude_options, ::CGAL::internal_np::exude_options_param_t, CGAL_NP_BASE>
+Named_function_parameters<::CGAL::parameters::internal::Perturb_options, ::CGAL::internal_np::perturb_options_param_t, CGAL_NP_BASE>
 perturb(const CGAL_NP_CLASS_1&  np1, const CGAL_NP_CLASS_2&  np2, const NP& ... nps)
 {
   return perturb(::CGAL::internal_np::combine_named_parameters(np1, np2, nps...));
@@ -91,7 +91,7 @@ exude(const CGAL_NP_CLASS& np = parameters::default_values())
   using ::CGAL::parameters::choose_parameter;
   using ::CGAL::parameters::get_parameter;
   double time_limit = choose_parameter(get_parameter(np,::CGAL::internal_np::maximum_running_time),::CGAL::parameters::internal::undef_parameter);
-  double sliver_bound = choose_parameter(get_parameter(np,::CGAL::internal_np::lower_sliver_bound),::CGAL::parameters::default_values_for_mesh_3::perturb_sliver_bound);
+  double sliver_bound = choose_parameter(get_parameter(np,::CGAL::internal_np::lower_sliver_bound),::CGAL::parameters::default_values_for_mesh_3::exude_sliver_bound);
 
   ::CGAL::parameters::internal::Exude_options options(true);
 
@@ -304,7 +304,7 @@ mesh_3_options(const CGAL_NP_CLASS& np = parameters::default_values())
   options.dump_after_refine_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_refine_prefix_param), "");
   options.dump_after_glob_opt_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_glob_opt_prefix_param), "");
   options.dump_after_perturb_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_perturb_prefix_param), "");
-  options.dump_after_exude_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_refine_surface_prefix_param), "");
+  options.dump_after_exude_prefix=choose_parameter(get_parameter(np, ::CGAL::internal_np::dump_after_exude_prefix_param), "");
   options.number_of_initial_points=choose_parameter(get_parameter(np, ::CGAL::internal_np::number_of_initial_points_param), -1);
   options.nonlinear_growth_of_balls = choose_parameter(get_parameter(np, ::CGAL::internal_np::nonlinear_growth_of_balls_param), false);
   options.maximal_number_of_vertices=choose_parameter(get_parameter(np, ::CGAL::internal_np::maximal_number_of_vertices_param), 0);

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -248,6 +248,8 @@ CGAL_add_named_parameter(smooth_constrained_edges_t, smooth_constrained_edges, s
 // List of named parameters used in Alpha_wrap_3
 CGAL_add_named_parameter(do_enforce_manifoldness_t, do_enforce_manifoldness, do_enforce_manifoldness)
 CGAL_add_named_parameter(seed_points_t, seed_points, seed_points)
+CGAL_add_named_parameter(refine_triangulation_t, refine_triangulation, refine_triangulation)
+CGAL_add_named_parameter(keep_inner_connected_components_t, keep_inner_connected_components, keep_inner_connected_components)
 
 // SMDS_3 parameters
 CGAL_add_named_parameter(surface_facets_t, surface_facets, surface_facets)
@@ -348,19 +350,14 @@ CGAL_add_named_parameter_with_compatibility(features_detector_param_t, features_
 CGAL_add_named_parameter_with_compatibility(input_features_param_t, input_features_param, input_features)
 
 CGAL_add_named_parameter_with_compatibility(edge_size_param_t, edge_size_param, edge_size)
-CGAL_add_named_parameter_with_compatibility_cref_only(edge_sizing_field_param_t, edge_sizing_field_param, edge_sizing_field)
 CGAL_add_named_parameter_with_compatibility(edge_min_size_param_t, edge_min_size_param, edge_min_size)
 CGAL_add_named_parameter_with_compatibility(facet_angle_param_t, facet_angle_param, facet_angle)
 CGAL_add_named_parameter_with_compatibility(facet_size_param_t, facet_size_param, facet_size)
-CGAL_add_named_parameter_with_compatibility_cref_only(facet_sizing_field_param_t, facet_sizing_field_param, facet_sizing_field)
 CGAL_add_named_parameter_with_compatibility_cref_only(facet_distance_param_t, facet_distance_param, facet_distance)
 CGAL_add_named_parameter_with_compatibility(facet_min_size_param_t, facet_min_size_param, facet_min_size)
 CGAL_add_named_parameter_with_compatibility(facet_topology_param_t, facet_topology_param, facet_topology)
-CGAL_add_named_parameter_with_compatibility(cell_radius_edge_param_t, cell_radius_edge_param, cell_radius_edge)
 CGAL_add_named_parameter_with_compatibility(cell_radius_edge_ratio_param_t, cell_radius_edge_ratio_param, cell_radius_edge_ratio)
 CGAL_add_named_parameter_with_compatibility(cell_size_param_t, cell_size_param, cell_size)
-CGAL_add_named_parameter_with_compatibility_cref_only(cell_sizing_field_param_t, cell_sizing_field_param, cell_sizing_field)
-CGAL_add_named_parameter_with_compatibility_cref_only(sizing_field_param_t, sizing_field_param, sizing_field)
 CGAL_add_named_parameter_with_compatibility(cell_min_size_param_t, cell_min_size_param, cell_min_size)
 
 CGAL_add_named_parameter_with_compatibility(function_param_t, function_param, function)

Scripts/developer_scripts/add_toc_to_github_wiki_page.py

@@ -1,125 +1,134 @@
-from sys import argv
-from sys import exit
 import codecs
 import re
 import argparse
+import sys
+from urllib.parse import quote
 
 parser = argparse.ArgumentParser()
-parser.add_argument("filename",
-                    help="the Mardown file to process")
-parser.add_argument("--codebase",
-                    help="for a Markdown file of Codebase instead of Github",
-                    action="store_true")
-parser.add_argument("--h1",
-                    help="support level one sections (h1)",
-                    action="store_true")
+parser.add_argument("filename", help="the Markdown file to process")
+parser.add_argument(
+    "--codebase",
+    help="for a Markdown file of Codebase instead of Github",
+    action="store_true",
+)
+parser.add_argument("--h1", help="support level one sections (h1)", action="store_true")
+parser.add_argument("--max-level", help="maximum level of sections", type=int, default = 5)
 args = parser.parse_args()
 
+
 # a probably incomplete version to generate an anchor from a section name
 def get_anchor(s):
-  s = s.replace("`","")
-  s = s.replace("(","")
-  s = s.replace(")","")
-  s = s.replace(".","")
-  s = s.replace("#","")
-  s = s.replace(":","")
-  s = s.replace(",","")
-  s = s.replace(";","")
-  if args.codebase:
-    s = s.replace("/","-")
-  else:
-    s = s.replace("/","")
-  s = s.replace("<","")
-  s = s.replace(">","")
-  s = s.replace("+","")
-  s = s.replace("=","")
-  s = s.replace("?","")
-  s = s.replace("@","")
-  s = s.lstrip(" ")
-  s = s.rstrip("\n")
-  s = s.rstrip(" ")
-  s = re.sub(r'\s+','-',s)
-  if not args.codebase:
-    s = s.lower()
-  if args.codebase:
-    s = s.replace("'","-and-39-")
-  return "#"+s
+    s = s.replace("`", "")
+    s = s.replace("(", "")
+    s = s.replace(")", "")
+    s = s.replace(".", "")
+    s = s.replace("#", "")
+    s = s.replace(":", "")
+    s = s.replace(",", "")
+    s = s.replace(";", "")
+    if args.codebase:
+        s = s.replace("/", "-")
+    else:
+        s = s.replace("/", "")
+    s = s.replace("<", "")
+    s = s.replace(">", "")
+    s = s.replace("+", "")
+    s = s.replace("=", "")
+    s = s.replace("?", "")
+    s = s.replace("@", "")
+    s = s.lstrip(" ")
+    s = s.rstrip("\n")
+    s = s.rstrip(" ")
+    s = re.sub(r"\s+", "-", s)
+    if not args.codebase:
+        s = s.lower()
+    if args.codebase:
+        s = s.replace("'", "-and-39-")
+    return "#" + quote(s)
+
 
 # indices the nesting level (first level allowed is ##)
 def get_level(s):
-  m = re.search('^(#+)\s', s)
-  if m:
-    return len(m.group(1))
-  else:
-    return 0
+    m = re.search(r"^(#+)\s", s)
+    if m:
+        return len(m.group(1))
+    else:
+        return 0
+
 
 def get_name(s):
-  m = re.search('^#+\s+(.*)\s*$', s)
-  if m:
-    return m.group(1)
-  else:
-    return "ERROR: Section name extraction"
-
-#generate the entry for one section
-def get_toc_entry(s):
-  name = get_name(s)
-  if args.h1:
-    level = get_level(s)-1
-  else:
-    level = get_level(s)-2
-  anchor = get_anchor(s)
-
-  if level<0:
-    return "ERROR: h1 sections are not allowed"
-
-  res="* ["+name+"]("+anchor+")"
-  for i in range(0,level):
-    res="  "+res
-  return res
-
-#now the main
-input = args.filename
-
-f = codecs.open(input, 'r', encoding='utf-8')
-
-if not f:
-  print("Cannot open "+input+"\n")
-  exit()
-
-#look for <!--TOC--> the begin of the file
-line=f.readline()
-if line.find("<!--TOC-->")==-1:
-  exit()
-
-#skip current TOC
-line=f.readline()
-while line and line.find("<!--TOC-->")==-1:
-  line=f.readline()
-
-if not line:
-  exit()
-
-buffer=""
-TOC="<!--TOC-->\n\n# Table of Contents\n"
-
-verbatim_mode=False # to ignore verbatim mode while looking for sections
-TOC_empty=True
-for line in f.readlines():
-  buffer+=line
-  if verbatim_mode:
-    if line[:3]=="```":
-      verbatim_mode=False
-  else:
-    if line[:3]=="```":
-      verbatim_mode=True
+    m = re.search(r"^#+\s+(.*)\s*$", s)
+    if m:
+        return m.group(1)
     else:
-      if line[0]=="#":
-        TOC+=(get_toc_entry(line)+"\n")
-        TOC_empty=False
-TOC+="\n<!--TOC-->\n"
+        return "ERROR: Section name extraction"
 
-if not TOC_empty:
-  f.close()
-  f = codecs.open(input, 'w', encoding='utf-8')
-  f.write(TOC)
-  f.write(buffer)
+
+# generate the entry for one section
+def get_toc_entry(s):
+    name = get_name(s)
+    if args.h1:
+        level = get_level(s) - 1
+    else:
+        level = get_level(s) - 2
+    anchor = get_anchor(s)
+
+    if level < 0:
+        return "ERROR: h1 sections are not allowed"
+
+    res = "* [" + name + "](" + anchor + ")"
+    for _ in range(0, level):
+        res = "  " + res
+    return res
+
+
+# now the main
+def main():
+    filename = args.filename
+
+    f = codecs.open(filename, "r", encoding="utf-8")
+
+    if not f:
+        print("Cannot open " + input + "\n")
+        sys.exit()
+
+    # look for <!--TOC--> the begin of the file
+    line = f.readline()
+    if line.find("<!--TOC-->") == -1:
+        sys.exit()
+
+    # skip current TOC
+    line = f.readline()
+    while line and line.find("<!--TOC-->") == -1:
+        line = f.readline()
+
+    if not line:
+        sys.exit()
+
+    buffer = ""
+    toc = "<!--TOC-->\n\n# Table of Contents\n"
+
+    verbatim_mode = False  # to ignore verbatim mode while looking for sections
+    toc_empty = True
+    for line in f.readlines():
+        buffer += line
+        if verbatim_mode:
+            if line[:3] == "```":
+                verbatim_mode = False
+        else:
+            if line[:3] == "```":
+                verbatim_mode = True
+            else:
+                if line[0] == "#" and get_level(line) <= args.max_level:
+                    toc += get_toc_entry(line) + "\n"
+                    toc_empty = False
+    toc += "\n<!--TOC-->\n"
+
+    if not toc_empty:
+        f.close()
+        f = codecs.open(filename, "w", encoding="utf-8")
+        f.write(toc)
+        f.write(buffer)
+
+if __name__ == "__main__":
+    main()

Shape_regularization/include/CGAL/Shape_regularization/regularize_segments.h

@@ -70,6 +70,9 @@ namespace Segments {
     that should be addressed. The function is based on the class `QP_regularization`.
     Please refer to that class and these concepts for more information.
 
+    This class requires a `QPSolver` model which defaults to the \ref thirdpartyOSQP "OSQP"
+    library, which must be available on the system.
+
     \tparam InputRange
     a model of `ConstRange` whose iterator type is `RandomAccessIterator`
 

Spatial_searching/doc/Spatial_searching/examples.txt

@@ -11,6 +11,7 @@
 \example Spatial_searching/searching_with_circular_query.cpp
 \example Spatial_searching/searching_surface_mesh_vertices.cpp
 \example Spatial_searching/searching_polyhedron_vertices.cpp
+\example Spatial_searching/searching_triangulation_vertices.cpp
 \example Spatial_searching/searching_with_point_with_info.cpp
 \example Spatial_searching/searching_with_point_with_info_inplace.cpp
 \example Spatial_searching/searching_with_point_with_info_pmap.cpp

Spatial_searching/examples/Spatial_searching/CMakeLists.txt

@@ -23,6 +23,7 @@ create_single_source_cgal_program("searching_with_point_with_info_inplace.cpp")
 create_single_source_cgal_program("searching_with_point_with_info_pmap.cpp")
 create_single_source_cgal_program("searching_surface_mesh_vertices.cpp")
 create_single_source_cgal_program("searching_polyhedron_vertices.cpp")
+create_single_source_cgal_program("searching_triangulation_vertices.cpp")
 create_single_source_cgal_program("searching_polyhedron_vertices_with_fuzzy_sphere.cpp")
 create_single_source_cgal_program("user_defined_point_and_distance.cpp")
 create_single_source_cgal_program("using_fair_splitting_rule.cpp")

Spatial_searching/examples/Spatial_searching/searching_triangulation_vertices.cpp

@@ -0,0 +1,65 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Search_traits_3.h>
+#include <CGAL/Search_traits_adapter.h>
+#include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+#include <CGAL/function_objects.h>
+
+#include <boost/property_map/function_property_map.hpp>
+
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+typedef Kernel::Point_3 Point_3;
+typedef CGAL::Delaunay_triangulation_3<Kernel> Delaunay;
+typedef Delaunay::Vertex_handle Vertex_handle;
+
+
+struct Project {
+  typedef Vertex_handle   argument_type;
+  typedef const Point_3&  Point;
+  typedef Point           result_type;
+  Point operator()( Vertex_handle  v) const { return v->point(); }
+};
+
+
+typedef boost::function_property_map<Project, Vertex_handle> Vertex_point_pmap;
+typedef CGAL::Search_traits_3<Kernel>                        Traits_base;
+typedef CGAL::Search_traits_adapter<Vertex_handle,Vertex_point_pmap,Traits_base> Traits;
+typedef CGAL::Orthogonal_k_neighbor_search<Traits>                               K_neighbor_search;
+typedef K_neighbor_search::Tree                                                  Tree;
+typedef Tree::Splitter                                                           Splitter;
+typedef K_neighbor_search::Distance                                              Distance;
+
+int main(int argc, char* argv[])
+{
+  std::ifstream in((argc>1)?argv[1]:CGAL::data_file_path("points_3/kitten.xyz"));
+  std::vector<Point_3> points;
+  Point_3 p;
+  while(in >> p){
+    points.push_back(p);
+  }
+
+  Delaunay dt(points.begin(), points.end());
+
+  const unsigned int K = 5;
+
+  Project project;
+  Vertex_point_pmap vppmap(project);
+
+  // Insert number_of_data_points in the tree
+  Tree tree(dt.finite_vertex_handles().begin(), dt.finite_vertex_handles().end(), Splitter(), Traits(vppmap));
+
+  // search K nearest neighbors
+  Point_3 query(0.0, 0.0, 0.0);
+  Distance tr_dist(vppmap);
+
+  K_neighbor_search search(tree, query, K,0,true,tr_dist);
+  std::cout <<"The "<< K << " nearest vertices to the query point at (0,0,0) are:" << std::endl;
+  for(K_neighbor_search::iterator it = search.begin(); it != search.end(); it++){
+    std::cout << "vertex " << &*(it->first) << " : " << vppmap[it->first] << " at distance "
+              << tr_dist.inverse_of_transformed_distance(it->second) << std::endl;
+  }
+
+  return 0;
+}

Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h

@@ -324,8 +324,8 @@ private:
 
     CGAL_precondition( (! (face(edge,mesh)==boost::graph_traits<Polyhedron>::null_face()))
                        && (! (face(opposite(edge,mesh),mesh)==boost::graph_traits<Polyhedron>::null_face())) );
-    const Point a = get(vertex_point_pmap,target(edge,mesh));
-    const Point b = get(vertex_point_pmap,target(prev(edge,mesh),mesh));
+    const Point a = get(vertex_point_pmap,source(edge,mesh));
+    const Point b = get(vertex_point_pmap,target(edge,mesh));
     const Point c = get(vertex_point_pmap,target(next(edge,mesh),mesh));
     const Point d = get(vertex_point_pmap,target(next(opposite(edge,mesh),mesh),mesh));
     // As far as I check: if, say, dihedral angle is 5, this returns 175,

TDS_3/include/CGAL/Triangulation_utils_3.h

@@ -80,6 +80,7 @@ struct Triangulation_utils_3
     CGAL_precondition( ( i >= 0 && i < 4 ) &&
                        ( j >= 0 && j < 4 ) &&
                        ( i != j ) );
+    CGAL_assume(i!=j);
     return tab_next_around_edge[i][j];
   }