Merge branch 'master' into Region_growing-revision-soesau

.github/dependabot.yml

@@ -0,0 +1,6 @@
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/"
+    schedule:
+      interval: "weekly"

.github/workflows/build_doc.yml

@@ -57,11 +57,11 @@ jobs:
           sudo apt-get update && sudo apt-get install -y graphviz ssh bibtex2html
           sudo pip install lxml
           sudo pip install 'pyquery==1.4.1' # it seems to be the last py2 compatible version
-          wget --no-verbose -O doxygen_exe https://cgal.geometryfactory.com/~mgimeno/doxygen/build_1_8_13/bin/doxygen
+          wget --no-verbose -O doxygen_exe https://cgal.geometryfactory.com/~cgaltest/doxygen_1_8_13_patched/doxygen
           sudo mv doxygen_exe /usr/bin/doxygen
           sudo chmod +x /usr/bin/doxygen
-          git config --global user.email "maxime.gimeno@geometryfactory.com"
-          git config --global user.name "Maxime Gimeno"
+          git config --global user.email "cgal@geometryfactory.com"
+          git config --global user.name "cgaltest"
 
       - name: configure all
         if: steps.get_round.outputs.result != 'stop'

.github/workflows/checks.yml

@@ -2,6 +2,9 @@ name: CMake Test Merge Branch
 
 on: [push, pull_request]
 
+permissions:
+  contents: read
+
 jobs:
   build:
 
@@ -16,12 +19,12 @@ jobs:
         sudo apt-get update && sudo apt-get install -y graphviz ssh bibtex2html
         sudo pip install lxml 
         sudo pip install 'pyquery==1.4.1' # it seems to be the last py2 compatible version
-        wget --no-verbose -O doxygen_exe https://cgal.geometryfactory.com/~mgimeno/doxygen/build_1_8_13/bin/doxygen
+        wget --no-verbose -O doxygen_exe https://cgal.geometryfactory.com/~cgaltest/doxygen_1_8_13_patched/doxygen
         sudo mv doxygen_exe /usr/bin/doxygen
         sudo chmod +x /usr/bin/doxygen
-        git config --global user.email "maxime.gimeno@geometryfactory.com"
-        git config --global user.name "Maxime Gimeno"
-        
+        git config --global user.email "cgal@geometryfactory.com"
+        git config --global user.name "cgaltest"
+
     - name: Run checks
       run: |
         zsh Scripts/developer_scripts/test_merge_of_branch HEAD

.github/workflows/cmake-all.yml

@@ -2,6 +2,9 @@ name: CMake Testsuite
 
 on: [push, pull_request]
 
+permissions:
+  contents: read
+
 jobs:
   cmake-testsuite:
 

.github/workflows/delete_doc.yml

@@ -3,9 +3,14 @@ name: Documentation Removal
 on:
   pull_request_target:
     types: [closed, removed]
+permissions:
+  contents: read
+
 jobs:
   build:
 
+    permissions:
+      contents: write  # for Git to git push
     runs-on: ubuntu-latest
 
     steps:
@@ -13,8 +18,8 @@ jobs:
     - name: delete directory
       run: |
          set -x
-         git config --global user.email "maxime.gimeno@geometryfactory.com"
-         git config --global user.name "Maxime Gimeno"
+         git config --global user.email "cgal@geometryfactory.com"
+         git config --global user.name "cgaltest"
          git clone https://maxGimeno:${{ secrets.PUSH_TO_CGAL_GITHUB_IO_TOKEN }}@github.com/CGAL/cgal.github.io.git --depth=5
          PR_NUMBER=$(python -c "import json; import os; y = json.load(open(os.environ['GITHUB_EVENT_PATH'])); print(y[\"number\"])")
          cd cgal.github.io/

.github/workflows/demo.yml

@@ -2,6 +2,9 @@ name: Test Polyhedron Demo
 
 on: [push, pull_request]
 
+permissions:
+  contents: read
+
 jobs:
   batch_1:
     runs-on: ubuntu-latest

.gitignore

@@ -2,6 +2,7 @@
 /*/*/*/build
 /*/*/*/VC*
 /*/*/*/GCC
+.vscode
 AABB_tree/demo/AABB_tree/AABB_demo
 AABB_tree/demo/AABB_tree/Makefile
 AABB_tree/examples/AABB_tree/*.kdev*

AABB_tree/benchmark/AABB_tree/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(AABB_traits_benchmark)
 
 find_package(CGAL REQUIRED OPTIONAL_COMPONENTS Core)

AABB_tree/demo/AABB_tree/CMakeLists.txt

@@ -1,6 +1,6 @@
 # This is the CMake script for compiling the AABB tree demo.
 
-cmake_minimum_required(VERSION 3.1...3.20)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(AABB_tree_Demo)
 
 # Find includes in corresponding build directories

AABB_tree/examples/AABB_tree/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(AABB_tree_Examples)
 
 find_package(CGAL REQUIRED)

AABB_tree/include/CGAL/AABB_primitive.h

@@ -53,9 +53,8 @@ public:
  * The two property maps which are template parameters of the class enable to get the datum and the reference point of
  * the primitive from the identifier. The last template parameter controls whether the primitive class holds a copy of the datum.
  *
- * \cgalModels `AABBPrimitive` if `ExternalPropertyMaps` is `CGAL::Tag_false`,
- *    and `AABBPrimitiveWithSharedData` if `ExternalPropertyMaps` is `CGAL::Tag_true`.
- *
+ * \cgalModels `AABBPrimitive` if `ExternalPropertyMaps` is `CGAL::Tag_false`.
+ * \cgalModels `AABBPrimitiveWithSharedData` if `ExternalPropertyMaps` is `CGAL::Tag_true`.
  *
  * \tparam ObjectPropertyMap is a model of `ReadablePropertyMap` with `Id` as
  *           `key_type`. It must be a model of `CopyConstructible`, `DefaultConstructible`, and `CopyAssignable`.
@@ -70,7 +69,6 @@ public:
  *           it is constructed on the fly to reduce the memory footprint.
  *           The default is `CGAL::Tag_false` (datum is not stored).
  *
- * \sa `AABBPrimitive`
  * \sa `AABB_segment_primitive<Iterator,CacheDatum>`
  * \sa `AABB_triangle_primitive<Iterator,CacheDatum>`
  * \sa `AABB_halfedge_graph_segment_primitive<HalfedgeGraph,OneHalfedgeGraphPerTree,CacheDatum>`

AABB_tree/include/CGAL/AABB_tree.h

@@ -562,11 +562,14 @@ public:
 
     /**
      * @brief Builds the tree by recursive expansion.
+     * @param node the root node of the subtree to generate
      * @param first the first primitive to insert
-     * @param last the last primitive to insert
+     * @param beyond the last primitive to insert
      * @param range the number of primitive of the range
+     * @param compute_bbox a functor
+     * @param split_primitives a functor
      *
-     * [first,last[ is the range of primitives to be added to the tree.
+     * [first,beyond[ is the range of primitives to be added to the tree.
      */
     template<typename ConstPrimitiveIterator, typename ComputeBbox, typename SplitPrimitives>
     void expand(Node& node,
@@ -574,8 +577,7 @@ public:
                 ConstPrimitiveIterator beyond,
                 const std::size_t range,
                 const ComputeBbox& compute_bbox,
-                const SplitPrimitives& split_primitives,
-                const AABBTraits&);
+                const SplitPrimitives& split_primitives);
 
   public:
     // returns a point which must be on one primitive
@@ -791,8 +793,7 @@ public:
                         ConstPrimitiveIterator beyond,
                         const std::size_t range,
                         const ComputeBbox& compute_bbox,
-                        const SplitPrimitives& split_primitives,
-                        const Tr& traits)
+                        const SplitPrimitives& split_primitives)
   {
     node.set_bbox(compute_bbox(first, beyond));
 
@@ -806,13 +807,13 @@ public:
       break;
     case 3:
       node.set_children(*first, new_node());
-      expand(node.right_child(), first+1, beyond, 2, compute_bbox, split_primitives, traits);
+      expand(node.right_child(), first+1, beyond, 2, compute_bbox, split_primitives);
       break;
     default:
       const std::size_t new_range = range/2;
       node.set_children(new_node(), new_node());
-      expand(node.left_child(), first, first + new_range, new_range, compute_bbox, split_primitives, traits);
-      expand(node.right_child(), first + new_range, beyond, range - new_range, compute_bbox, split_primitives, traits);
+      expand(node.left_child(), first, first + new_range, new_range, compute_bbox, split_primitives);
+      expand(node.right_child(), first + new_range, beyond, range - new_range, compute_bbox, split_primitives);
     }
   }
 
@@ -844,8 +845,7 @@ public:
              m_primitives.begin(), m_primitives.end(),
              m_primitives.size(),
              compute_bbox,
-             split_primitives,
-             m_traits);
+             split_primitives);
     }
 #ifdef CGAL_HAS_THREADS
     m_atomic_need_build.store(false, std::memory_order_release); // in case build() is triggered by a call to root_node()

AABB_tree/include/CGAL/AABB_tree/internal/triangle_datum_covering.h

@@ -0,0 +1,406 @@
+// Copyright (c) 2021 INRIA Sophia-Antipolis (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) : Mael Rouxel-Labbé
+
+#ifndef CGAL_AABB_TREE_INTERNAL_TRIANGLE_DATUM_COVERING_H
+#define CGAL_AABB_TREE_INTERNAL_TRIANGLE_DATUM_COVERING_H
+
+#include <CGAL/license/AABB_tree.h>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Cartesian_converter.h>
+
+#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
+#include <CGAL/AABB_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/array.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Container_helper.h>
+#include <CGAL/property_map.h>
+
+#include <queue>
+#include <unordered_set>
+#include <utility>
+#include <vector>
+
+namespace CGAL {
+namespace AABB_trees {
+namespace internal {
+
+// std::vector to property map
+// Not using Pointer_property_map because the underlying range is empty initially and can change
+template <typename T>
+struct Vector_property_map
+{
+  using Range = std::vector<T>;
+
+  using key_type = std::size_t;
+  using value_type = T;
+  using reference = value_type&;
+  using category = boost::read_write_property_map_tag;
+
+  Vector_property_map() : m_range_ptr(std::make_shared<Range>()) { }
+
+  inline friend void put(const Vector_property_map& map, const key_type& k, const value_type& v)
+  {
+    CGAL_precondition(map.m_range_ptr != nullptr);
+
+    if(k >= map.m_range_ptr->size())
+      map.m_range_ptr->resize(k+1);
+
+    map.m_range_ptr->operator[](k) = v;
+  }
+
+  inline friend reference get(const Vector_property_map& map, const key_type& k)
+  {
+    CGAL_precondition(map.m_range_ptr != nullptr);
+    return map.m_range_ptr->operator[](k);
+  }
+
+  Range& range() { return *m_range_ptr; }
+  const Range& range() const { return *m_range_ptr; }
+
+private:
+  std::shared_ptr<Range> m_range_ptr;
+};
+
+// Same as the standard traversal traits, but for multiple primitives per datum,
+// such that the final operation on the datum is only performed once.
+template <typename AABBTraits, typename BaseTraversalTraits>
+struct Covered_traversal_traits
+  : BaseTraversalTraits
+{
+  using Base = BaseTraversalTraits;
+  using Primitive = typename AABBTraits::Primitive;
+
+  std::unordered_set<std::size_t> visited_data;
+
+public:
+  template <typename ... Args>
+  Covered_traversal_traits(Args&&... args) : Base(std::forward<Args>(args)...) { }
+
+  template <typename Query>
+  void intersection(const Query& query, const Primitive& primitive)
+  {
+    // check a datum only once
+    auto is_insert_successful = visited_data.insert(primitive.id().second/*unique input face ID*/);
+    if(!is_insert_successful.second)
+      return;
+
+    return Base::intersection(query, primitive);
+  }
+};
+
+template <typename AABBTraits>
+struct Covered_tree_traversal_traits
+{
+  using Base_projection_traits = CGAL::internal::AABB_tree::Projection_traits<AABBTraits>;
+  using Projection_traits = Covered_traversal_traits<AABBTraits, Base_projection_traits>;
+
+  template <typename Query>
+  using Do_intersect_traits_base = CGAL::internal::AABB_tree::Do_intersect_traits<AABBTraits, Query>;
+  template <typename Query>
+  using Do_intersect_traits = Covered_traversal_traits<AABBTraits, Do_intersect_traits_base<Query> >;
+
+  template <typename Query>
+  using First_intersection_traits_base = CGAL::internal::AABB_tree::First_intersection_traits<AABBTraits, Query>;
+  template <typename Query>
+  using First_intersection_traits = Covered_traversal_traits<AABBTraits, First_intersection_traits_base<Query> >;
+
+  template <typename Query>
+  using First_primitive_traits_base = CGAL::internal::AABB_tree::First_primitive_traits<AABBTraits, Query>;
+  template <typename Query>
+  using First_primitive_traits = Covered_traversal_traits<AABBTraits, First_primitive_traits_base<Query> >;
+
+  template <typename Query, typename CountingIterator>
+  using Listing_primitive_traits_base = CGAL::internal::AABB_tree::Listing_primitive_traits<AABBTraits, Query, CountingIterator>;
+  template <typename Query, typename CountingIterator>
+  using Listing_primitive_traits = Covered_traversal_traits<AABBTraits, Listing_primitive_traits_base<Query, CountingIterator> >;
+
+  template <typename Query, typename CountingIterator>
+  using Listing_intersection_traits_base = CGAL::internal::AABB_tree::Listing_intersection_traits<AABBTraits, Query, CountingIterator>;
+  template <typename Query, typename CountingIterator>
+  using Listing_intersection_traits = Covered_traversal_traits<AABBTraits, Listing_intersection_traits_base<Query, CountingIterator> >;
+};
+
+// Dissociated from the class `AABB_covered_triangle_tree` for clarity
+template <typename Kernel, typename Point>
+struct AABB_covered_triangle_tree_traits
+{
+  using Triangle_3 = typename Kernel::Triangle_3;
+
+  // Below is a lot of trouble to cover a single datum with multiple primitives using smaller bboxes
+  using ID = std::pair<std::size_t /*primitive ID*/, std::size_t /*input face ID*/>;
+  using IDPM = CGAL::First_of_pair_property_map<ID>;
+
+  // Primitive ID --> box vector pos --> Bounding Box
+  using BPMB = internal::Vector_property_map<CGAL::Bbox_3>;
+  using BPM = CGAL::Property_map_binder<IDPM, BPMB>;
+
+  // Primitive ID --> point vector pos --> Reference Point
+  using RPPMB = internal::Vector_property_map<Point>;
+  using RPPM = CGAL::Property_map_binder<IDPM, RPPMB>;
+
+  // Primitive ID --> Datum pos vector pos --> Datum pos --> Datum
+  // The vector of data has size nf, but the vector of datum pos has size tree.size()
+  using DPPMB = internal::Vector_property_map<std::size_t>; // pos --> Datum pos
+  using DPPM = CGAL::Property_map_binder<IDPM, DPPMB>; // PID --> Datum pos
+  using DPMB = internal::Vector_property_map<Triangle_3>; // Datum pos --> Datum
+  using DPM = CGAL::Property_map_binder<DPPM, DPMB>; // PID --> Datum
+
+  using Primitive = CGAL::AABB_primitive<ID, DPM, RPPM,
+                                         CGAL::Tag_true /*external pmaps*/,
+                                         CGAL::Tag_false /*no caching*/>;
+
+  using AABB_geom_traits = Kernel;
+  using AABB_traits = CGAL::AABB_traits<AABB_geom_traits, Primitive, BPM>;
+  using AABB_tree = CGAL::AABB_tree<AABB_traits>;
+};
+
+template <typename Kernel,
+          typename Point = typename Kernel::Point_3>
+struct AABB_covered_triangle_tree
+  : public AABB_covered_triangle_tree_traits<Kernel, Point>::AABB_tree
+{
+  using FT = typename Kernel::FT;
+  using Triangle_3 = typename Kernel::Triangle_3;
+
+  using ACTTT = AABB_covered_triangle_tree_traits<Kernel, Point>;
+
+  using BPM = typename ACTTT::BPM;
+  using RPPM = typename ACTTT::RPPM;
+  using DPPMB = typename ACTTT::DPPMB;
+  using DPPM = typename ACTTT::DPPM;
+  using DPMB = typename ACTTT::DPMB;
+  using DPM = typename ACTTT::DPM;
+
+  using ID = typename ACTTT::ID;
+  using Primitive = typename ACTTT::Primitive;
+  using AABB_traits = typename ACTTT::AABB_traits;
+  using AABB_tree = typename ACTTT::AABB_tree;
+  using Base = AABB_tree;
+
+protected:
+  double m_sq_length;
+
+  DPPMB m_dppmb; // std::size_t (id) --> datum pos
+
+  BPM m_bpm; // std::size_t (id) --> bounding box
+  RPPM m_rppm; // std::size_t (id) --> reference point
+  DPMB m_dpmb; // std::size_t (datum pos) --> triangle datum
+
+  DPM m_dpm; // std::size_t (id) --> triangle (datum)
+
+  std::size_t fid = 0;
+
+public:
+  AABB_covered_triangle_tree(const double max_length,
+                             const AABB_traits& traits = AABB_traits())
+    : Base(traits),
+      m_sq_length(square(max_length)),
+      m_dppmb(), m_bpm(), m_rppm(), m_dpmb(),
+      m_dpm(DPPM(m_dppmb/*first binder's value_map*/)/*second binder's key map*/, m_dpmb)
+  {
+    initialize_tree_property_maps();
+  }
+
+private:
+  void initialize_tree_property_maps() const
+  {
+    // Can't be set in the default constructed traits that are passed to the base
+    // since m_bpm is a member of the derived class.
+    //
+    // 'const_cast' because CGAL::AABB_tree only gives a const& to its traits...
+    const_cast<AABB_traits&>(this->traits()).bbm = m_bpm;
+    const_cast<AABB_traits&>(this->traits()).set_shared_data(m_dpm, m_rppm);
+  }
+
+  template <typename P> // Kernel is Simple_Cartesian<Interval>
+  CGAL::Bbox_3 compute_bbox(const P& ap0, const P& ap1, const P& ap2)
+  {
+    double xmin = (CGAL::min)(ap0.x().inf(), (CGAL::min)(ap1.x().inf(), ap2.x().inf()));
+    double ymin = (CGAL::min)(ap0.y().inf(), (CGAL::min)(ap1.y().inf(), ap2.y().inf()));
+    double zmin = (CGAL::min)(ap0.z().inf(), (CGAL::min)(ap1.z().inf(), ap2.z().inf()));
+
+    double xmax = (CGAL::max)(ap0.x().sup(), (CGAL::max)(ap1.x().sup(), ap2.x().sup()));
+    double ymax = (CGAL::max)(ap0.y().sup(), (CGAL::max)(ap1.y().sup(), ap2.y().sup()));
+    double zmax = (CGAL::max)(ap0.z().sup(), (CGAL::max)(ap1.z().sup(), ap2.z().sup()));
+
+    return CGAL::Bbox_3(xmin, ymin, zmin, xmax, ymax, zmax);
+  }
+
+  template <typename P> // Kernel is Simple_Cartesian<Interval>
+  const P& compute_reference_point(const P&, const P&, const P& ap2)
+  {
+    return ap2; // ap2 is the midpoint when splitting
+  }
+
+public:
+  void reserve(std::size_t nf)
+  {
+    CGAL::internal::reserve(m_dpmb.range(), m_dpmb.range().size() + nf);
+
+    // Due to splitting, these might need more than 'nf'
+    CGAL::internal::reserve(m_dppmb.range(), m_dppmb.range().size() + nf);
+    CGAL::internal::reserve(m_rppm.value_map.range(), m_rppm.value_map.range().size() + nf);
+    CGAL::internal::reserve(m_bpm.value_map.range(), m_bpm.value_map.range().size() + nf);
+  }
+
+  void split_and_insert(const Triangle_3& tr)
+  {
+    // Convert to intervals to ensure that the bounding box fully covers the subdividing triangle
+    using IT = CGAL::Interval_nt<true>;
+    using NT = typename IT::value_type;
+    using AK = CGAL::Simple_cartesian<IT>;
+    using K2AK = CGAL::Cartesian_converter<Kernel, AK>;
+    using AFT = AK::FT;
+    using APoint_3 = AK::Point_3;
+
+    using APL = std::pair<APoint_3, NT>; // point and upper bound of the length of the opposite edge
+    using AT = std::array<APL, 3>;
+
+    const std::size_t data_size = m_dpmb.range().size();
+    put(m_dpmb, data_size, tr);
+
+    auto vertex = Kernel().construct_vertex_3_object();
+    const Point& p0 = vertex(tr, 0);
+    const Point& p1 = vertex(tr, 1);
+    const Point& p2 = vertex(tr, 2);
+
+    if(m_sq_length == FT(0)) // no splits
+    {
+      const std::size_t pid = this->size();
+      ID id = std::make_pair(pid, fid++);
+
+      put(m_dppmb, pid, data_size);
+      put(m_rppm, id, p1); // the ref point that `One_point_from_face_descriptor_map` would give
+      put(m_bpm, id, Kernel().construct_bbox_3_object()(tr));
+
+//      std::cout << "Primitive[" << id.first << " " << id.second << "]; "
+//                    << "Bbox: [" << get(m_bpm, id) << "] "
+//                    << "Point: (" << get(m_rppm, id) << ") "
+//                    << "Datum: [" << get(m_bpm, id) << "]" << std::endl;
+
+      Primitive p(id/*, m_dpm, m_rppm*/); // pmaps are external, shared data
+      this->insert(p);
+      return;
+    }
+
+    K2AK k2ak;
+    std::queue<AT> to_treat;
+
+    const APoint_3 ap0 = k2ak(p0);
+    const APoint_3 ap1 = k2ak(p1);
+    const APoint_3 ap2 = k2ak(p2);
+    const AFT sq_l0 = CGAL::squared_distance(ap1, ap2);
+    const AFT sq_l1 = CGAL::squared_distance(ap2, ap0);
+    const AFT sq_l2 = CGAL::squared_distance(ap0, ap1);
+
+    to_treat.push(CGAL::make_array(std::make_pair(ap0, sq_l0.sup()),
+                                   std::make_pair(ap1, sq_l1.sup()),
+                                   std::make_pair(ap2, sq_l2.sup())));
+
+    while(!to_treat.empty())
+    {
+      const AT t = std::move(to_treat.front());
+      to_treat.pop();
+
+      const APL& apl0 = t[0];
+      const APL& apl1 = t[1];
+      const APL& apl2 = t[2];
+
+      int i = (apl0.second >= apl1.second)
+                ? (apl0.second >= apl2.second) ? 0 : 2
+                : (apl1.second >= apl2.second) ? 1 : 2;
+
+      const NT max_sql = t[i].second;
+
+      // If the face is too big, do a split (two small bboxes rather than a big one)
+      if(max_sql > m_sq_length)
+      {
+        // Could be factorized, but this is simpler to read
+        if(i == 0)
+        {
+          // 0 1 2 into 0 1 A and 0 A 2
+          const APoint_3 amp = CGAL::midpoint(apl1.first, apl2.first);
+          const NT sq_half_length = apl0.second / NT(4);
+          const NT sq_diag_length = CGAL::squared_distance(amp, apl0.first).sup();
+
+          to_treat.push(CGAL::make_array(std::make_pair(apl0.first, sq_half_length),
+                                         std::make_pair(apl1.first, sq_diag_length),
+                                         std::make_pair(amp, apl2.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl2.first, sq_diag_length),
+                                         std::make_pair(apl0.first, sq_half_length),
+                                         std::make_pair(amp, apl1.second)));
+        }
+        else if(i == 1)
+        {
+          // 0 1 2 into 0 1 A and 1 2 A
+          const APoint_3 amp = CGAL::midpoint(apl2.first, apl0.first);
+          const NT sq_half_length = apl1.second / NT(4);
+          const NT sq_diag_length = CGAL::squared_distance(amp, apl1.first).sup();
+
+          to_treat.push(CGAL::make_array(std::make_pair(apl0.first, sq_diag_length),
+                                         std::make_pair(apl1.first, sq_half_length),
+                                         std::make_pair(amp, apl2.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl1.first, sq_half_length),
+                                         std::make_pair(apl2.first, sq_diag_length),
+                                         std::make_pair(amp, apl0.second)));
+        }
+        else // i == 2
+        {
+          // 0 1 2 into 0 A 2 and 2 A 1
+          const APoint_3 amp = CGAL::midpoint(apl0.first, apl1.first);
+          const NT sq_half_length = apl2.second / NT(4);
+          const NT sq_diag_length = CGAL::squared_distance(amp, apl2.first).sup();
+
+          to_treat.push(CGAL::make_array(std::make_pair(apl2.first, sq_half_length),
+                                         std::make_pair(apl0.first, sq_diag_length),
+                                         std::make_pair(amp, apl1.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl1.first, sq_diag_length),
+                                         std::make_pair(apl2.first, sq_half_length),
+                                         std::make_pair(amp, apl0.second)));
+        }
+      }
+      else // all edges have length below the threshold, create a primitive
+      {
+        const std::size_t pid = this->size();
+        ID id = std::make_pair(pid, fid);
+
+        put(m_dppmb, pid, data_size);
+
+        // this is basically to_double() of an APoint_3, but FT is not necessarily 'double'
+        const APoint_3& apt = compute_reference_point(apl0.first, apl1.first, apl2.first);
+        put(m_rppm, id, Point((FT(apt.x().sup()) + FT(apt.x().inf())) / FT(2),
+                              (FT(apt.y().sup()) + FT(apt.y().inf())) / FT(2),
+                              (FT(apt.z().sup()) + FT(apt.z().inf())) / FT(2)));
+
+        put(m_bpm, id, compute_bbox(apl0.first, apl1.first, apl2.first));
+
+//        std::cout << "Primitive[" << std::get<0>(id) << " " << std::get<1>(id) << " " << std::get<2>(id) << "]; "
+//                  << "Bbox: [" << get(m_bpm, id) << "] "
+//                  << "Point: (" << get(m_rppm, id) << ") "
+//                  << "Datum: [" << get(m_dpm, id) << "]" << std::endl;
+
+        Primitive p(id/*, m_dpm, m_rppm*/); // pmaps are external, shared data
+        this->insert(p);
+      }
+    }
+
+    ++fid;
+  }
+};
+
+} // namespace internal
+} // namespace AABB_trees
+} // namespace CGAL
+
+#endif // CGAL_AABB_TREE_INTERNAL_TRIANGLE_DATUM_COVERING_H

AABB_tree/test/AABB_tree/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(AABB_tree_Tests)
 
 find_package(CGAL REQUIRED)

AABB_tree/test/AABB_tree/aabb_any_all_benchmark.cpp

@@ -1,7 +1,7 @@
 #include <iostream>
 #include <fstream>
 #include <algorithm>
-#include <boost/tuple/tuple.hpp>
+#include <tuple>
 
 #include <CGAL/Cartesian.h>
 #include <CGAL/Simple_cartesian.h>
@@ -55,7 +55,7 @@ std::size_t intersect(ForwardIterator b, ForwardIterator e, const Tree& tree, lo
 }
 
 template<typename K>
-boost::tuple<std::size_t, std::size_t, std::size_t, long> test(const char* name) {
+std::tuple<std::size_t, std::size_t, std::size_t, long> test(const char* name) {
   typedef typename K::FT FT;
   typedef typename K::Ray_3 Ray;
   typedef typename K::Line_3 Line;
@@ -119,7 +119,7 @@ boost::tuple<std::size_t, std::size_t, std::size_t, long> test(const char* name)
 
   segments.erase(std::remove_if(segments.begin(), segments.end(), p), segments.end());
 
-  boost::tuple<std::size_t, std::size_t, std::size_t, long> tu;
+  std::tuple<std::size_t, std::size_t, std::size_t, long> tu;
 
     {
       CGAL::Timer t;
@@ -128,12 +128,12 @@ boost::tuple<std::size_t, std::size_t, std::size_t, long> test(const char* name)
       for(int i = 0; i < runs; ++i)
       {
         long counter = 0L;
-        tu = boost::make_tuple(intersect(lines.begin(), lines.end(), tree, counter),
+        tu = std::make_tuple(intersect(lines.begin(), lines.end(), tree, counter),
                                intersect(rays.begin(), rays.end(), tree, counter),
                                intersect(segments.begin(), segments.end(), tree, counter),
                                // cant use counter here
                                0);
-        boost::get<3>(tu) = counter;
+        std::get<3>(tu) = counter;
       }
       std::cout << t.time();
     }
@@ -146,41 +146,41 @@ int main()
   const char* filename = "data/finger.off";
 
   std::cout << "| Simple cartesian float kernel | ";
-  boost::tuple<std::size_t, std::size_t, std::size_t, long> t1 = test<CGAL::Simple_cartesian<float> >(filename);
+  std::tuple<std::size_t, std::size_t, std::size_t, long> t1 = test<CGAL::Simple_cartesian<float> >(filename);
   std::cout << " | " << std::endl;
 
   std::cout << "| Cartesian float kernel | ";
-  boost::tuple<std::size_t, std::size_t, std::size_t, long> t2 = test<CGAL::Cartesian<float> >(filename);
+  std::tuple<std::size_t, std::size_t, std::size_t, long> t2 = test<CGAL::Cartesian<float> >(filename);
   std::cout << " | " << std::endl;
 
   std::cout << "| Simple cartesian double kernel |";
-  boost::tuple<std::size_t, std::size_t, std::size_t, long> t3 = test<CGAL::Simple_cartesian<double> >(filename);
+  std::tuple<std::size_t, std::size_t, std::size_t, long> t3 = test<CGAL::Simple_cartesian<double> >(filename);
   std::cout << " | " << std::endl;
 
   std::cout << "| Cartesian double kernel |";
-  boost::tuple<std::size_t, std::size_t, std::size_t, long> t4 = test<CGAL::Cartesian<double> >(filename);
+  std::tuple<std::size_t, std::size_t, std::size_t, long> t4 = test<CGAL::Cartesian<double> >(filename);
   std::cout << " | " << std::endl;
 
   std::cout << "| Epic kernel |";
-  boost::tuple<std::size_t, std::size_t, std::size_t, long> t5 = test<CGAL::Exact_predicates_inexact_constructions_kernel>(filename);
+  std::tuple<std::size_t, std::size_t, std::size_t, long> t5 = test<CGAL::Exact_predicates_inexact_constructions_kernel>(filename);
   std::cout << " | " << std::endl;
 
   std::size_t a, b, c;
   long d;
 
-  boost::tie(a, b, c, d) = t5;
+  std::tie(a, b, c, d) = t5;
   std::cout << a << " " << b << " " << c << " " << d << std::endl;
 
-  boost::tie(a, b, c, d) = t4;
+  std::tie(a, b, c, d) = t4;
   std::cout << a << " " << b << " " << c << " " << d << std::endl;
 
-  boost::tie(a, b, c, d) = t3;
+  std::tie(a, b, c, d) = t3;
   std::cout << a << " " << b << " " << c << " " << d << std::endl;
 
-  boost::tie(a, b, c, d) = t2;
+  std::tie(a, b, c, d) = t2;
   std::cout << a << " " << b << " " << c << " " << d << std::endl;
 
-  boost::tie(a, b, c, d) = t1;
+  std::tie(a, b, c, d) = t1;
   std::cout << a << " " << b << " " << c << " " << d << std::endl;
   return 0;
 }

AABB_tree/test/AABB_tree/aabb_test_all_intersected_primitives_tree.cpp

@@ -1,7 +1,6 @@
 #include <fstream>
 #include <iterator>
-
-#include <CGAL/assertions.h>
+#include <cassert>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
@@ -59,11 +58,11 @@ int main()
   std::list<T_Tree::Primitive::Id> t_primitives;
   std::list<S_Tree::Primitive::Id> s_primitives;
   cube_tree.all_intersected_primitives(tet_tree,std::back_inserter(t_primitives));
-  CGAL_assertion(t_primitives.size() == 6);
+  assert(t_primitives.size() == 6);
   tet_tree.all_intersected_primitives(cube_tree,std::back_inserter(s_primitives));
-  CGAL_assertion(s_primitives.size() == 6);
-  CGAL_assertion(tet_tree.do_intersect(cube_tree));
-  CGAL_assertion(cube_tree.do_intersect(tet_tree));
+  assert(s_primitives.size() == 6);
+  assert(tet_tree.do_intersect(cube_tree));
+  assert(cube_tree.do_intersect(tet_tree));
 
   std::vector<T_Tree::Primitive::Id> all_primitives;
   cube_tree.all_intersected_primitives(tet_tree, std::back_inserter(all_primitives));
@@ -73,7 +72,7 @@ int main()
     if((int)prim.first == 5)
       found_f5 = true;
   }
-  CGAL_assertion(found_f5);
-  CGAL_USE(found_f5);
+  assert(found_f5);
+
   return 0;
 }

AABB_tree/test/AABB_tree/aabb_test_move_constructor.cpp

@@ -2,6 +2,7 @@
 #include <iostream>
 #include <list>
 #include <utility>
+#include <cassert>
 
 #include <CGAL/AABB_face_graph_triangle_primitive.h>
 #include <CGAL/AABB_halfedge_graph_segment_primitive.h>
@@ -149,12 +150,12 @@ public:
     std::list<S_Tree::Primitive::Id> s_primitives;
     cube_tree.all_intersected_primitives(tet_tree,
                                          std::back_inserter(t_primitives));
-    CGAL_assertion(t_primitives.size() == 6);
+    assert(t_primitives.size() == 6);
     tet_tree.all_intersected_primitives(cube_tree,
                                         std::back_inserter(s_primitives));
-    CGAL_assertion(s_primitives.size() == 6);
-    CGAL_assertion(tet_tree.do_intersect(cube_tree));
-    CGAL_assertion(cube_tree.do_intersect(tet_tree));
+    assert(s_primitives.size() == 6);
+    assert(tet_tree.do_intersect(cube_tree));
+    assert(cube_tree.do_intersect(tet_tree));
 
     std::vector<T_Tree::Primitive::Id> all_primitives;
     cube_tree.all_intersected_primitives(tet_tree,
@@ -164,8 +165,8 @@ public:
       if ((int)prim.first == 5)
         found_f5 = true;
     }
-    CGAL_assertion(found_f5);
-    CGAL_USE(found_f5);
+    assert(found_f5);
+
     return true;
   }
 };

AABB_tree/test/AABB_tree/aabb_test_multi_mesh.cpp

@@ -5,7 +5,6 @@
 #include <boost/functional/value_factory.hpp>
 #include <boost/array.hpp>
 
-#include <CGAL/assertions.h>
 #include <CGAL/algorithm.h>
 #include <CGAL/point_generators_3.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

AABB_tree/test/AABB_tree/aabb_test_ray_intersection.cpp

@@ -138,9 +138,8 @@ int main()
   for(std::vector<Ray>::iterator it = rays.begin(); it != rays.end(); ++it) {
     primitives2.push_back(tree.first_intersection(*it));
   }
-  CGAL_assertion_msg(primitives1.size() == primitives2.size(), "Different amount of primitives intersected.");
-  CGAL_assertion_msg(std::equal(primitives1.begin(), primitives1.end(), primitives2.begin()),
-                     "Primitives mismatch.");
+  assert(primitives1.size() == primitives2.size()); //  Different amount of primitives intersected
+  assert(std::equal(primitives1.begin(), primitives1.end(), primitives2.begin())); //  Primitives mismatch
   std::size_t c = primitives1.size() - std::count(primitives1.begin(), primitives1.end(), boost::none);
   std::cout << "Intersected " << c << " primitives with " << NB_RAYS << " rays" << std::endl;
   std::cout << "Primitive method had to sort " << accum/NB_RAYS

AABB_tree/test/AABB_tree/aabb_triangle_datum_covering.cpp

@@ -0,0 +1,262 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/AABB_tree/internal/triangle_datum_covering.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+
+#include <CGAL/point_generators_3.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
+#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
+#include <CGAL/Real_timer.h>
+
+#include <fstream>
+#include <iostream>
+
+constexpr int N = 1000;
+
+template <typename Mesh>
+void test_no_cover(const Mesh& mesh)
+{
+  std::cout << "Test WITHOUT cover" << std::endl;
+
+  using GT = typename CGAL::GetGeomTraits<Mesh>::type;
+  using FT = typename GT::FT;
+  using Point = typename GT::Point_3;
+  using Ray_3 = typename GT::Ray_3;
+  using Line_3 = typename GT::Line_3;
+
+  using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh, CGAL::Default, CGAL::Tag_true, CGAL::Tag_true>;
+  using Traits = CGAL::AABB_traits<GT, Primitive>;
+  using Tree = CGAL::AABB_tree<Traits>;
+
+  // Build
+  Tree tree(faces(mesh).first, faces(mesh).second, mesh);
+  std::cout << faces(mesh).size() << " input faces and " << tree.size() << " primitives" << std::endl;
+
+  // Usage
+  std::cout << "Traverse..." << std::endl;
+
+  CGAL::Random r; // = CGAL::Random(1337);
+  CGAL::Random_points_in_cube_3<Point, CGAL::Creator_uniform_3<FT, Point> > g(1., r);
+
+  std::vector<Point> points;
+  points.reserve(N);
+  std::copy_n(g, N, std::back_inserter(points));
+
+  CGAL::Real_timer timer;
+  timer.start();
+
+  for(int i=0; i<N; ++i)
+  {
+    bool does_intersect = tree.do_intersect(Line_3(CGAL::ORIGIN, points[i]));
+    assert(does_intersect);
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    auto inter_res = tree.any_intersection(Ray_3(CGAL::ORIGIN, points[i]));
+    assert(inter_res);
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    auto inter_res = tree.first_intersection(Ray_3(CGAL::ORIGIN, points[i]));
+    assert(inter_res);
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    Point res = tree.closest_point(points[i]);
+    CGAL_USE(res);
+  }
+
+  for(int i=0; i<N; ++i)
+    tree.all_intersected_primitives(Line_3(CGAL::ORIGIN, points[i]), CGAL::Emptyset_iterator());
+
+  timer.stop();
+  std::cout << "Elapsed: " << timer.time() << " s." << std::endl;
+}
+
+template <typename Mesh, typename NamedParameters>
+void test_cover(const Mesh& mesh,
+                const NamedParameters& np)
+{
+  std::cout << "Test WITH cover" << std::endl;
+
+  constexpr bool check_with_standard_tree = true;
+
+  using VPM = typename CGAL::GetVertexPointMap<Mesh, NamedParameters>::const_type;
+  using GT = typename CGAL::GetGeomTraits<Mesh, NamedParameters>::type;
+  using FT = typename GT::FT;
+  using Point = typename boost::property_traits<VPM>::value_type;
+  using Ray_3 = typename GT::Ray_3;
+  using Line_3 = typename GT::Line_3;
+  using Triangle_3 = typename GT::Triangle_3;
+
+  using AABB_tree = CGAL::AABB_trees::internal::AABB_covered_triangle_tree<GT, Point>;
+
+  using FG_Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh, CGAL::Default, CGAL::Tag_true, CGAL::Tag_true>;
+  using FG_Traits = CGAL::AABB_traits<GT, FG_Primitive>;
+  using FG_Tree = CGAL::AABB_tree<FG_Traits>;
+
+  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()));
+
+  // Build -----------------------------------------------------------------------------------------
+  using CGAL::parameters::get_parameter;
+  using CGAL::parameters::choose_parameter;
+
+  VPM vpm = choose_parameter(get_parameter(np, CGAL::internal_np::vertex_point),
+                             get_const_property_map(CGAL::vertex_point, mesh));
+  GT gt = choose_parameter<GT>(get_parameter(np, CGAL::internal_np::geom_traits));
+
+  AABB_tree tree(diag_length / 100.);
+  tree.reserve(num_faces(mesh));
+
+  for(auto f : faces(mesh))
+  {
+    if(CGAL::Polygon_mesh_processing::is_degenerate_triangle_face(f, mesh, np))
+      continue;
+
+    const Point& p0 = get(vpm, target(halfedge(f, mesh), mesh));
+    const Point& p1 = get(vpm, target(next(halfedge(f, mesh), mesh), mesh));
+    const Point& p2 = get(vpm, source(halfedge(f, mesh), mesh));
+
+    const Triangle_3 tr = gt.construct_triangle_3_object()(p0, p1, p2);
+    tree.split_and_insert(tr);
+  }
+
+  std::cout << faces(mesh).size() << " input faces and " << tree.size() << " primitives" << std::endl;
+
+  FG_Tree fg_tree(faces(mesh).first, faces(mesh).second, mesh);
+
+  // Usage -----------------------------------------------------------------------------------------
+  std::cout << "Traverse..." << std::endl;
+
+  CGAL::Random r; // = CGAL::Random(1337);
+  CGAL::Random_points_in_cube_3<Point, CGAL::Creator_uniform_3<FT, Point> > g(1., r);
+
+  std::vector<Point> points;
+  points.reserve(N);
+  std::copy_n(g, N, std::back_inserter(points));
+
+  using AABB_traits = typename AABB_tree::AABB_traits;
+  using AABB_traversal_traits = CGAL::AABB_trees::internal::Covered_tree_traversal_traits<AABB_traits>;
+
+  using Do_intersect_traits = typename AABB_traversal_traits::template Do_intersect_traits<Line_3>;
+  using First_intersection_traits = typename AABB_traversal_traits::template First_intersection_traits<Ray_3>;
+  using Projection_traits = typename AABB_traversal_traits::Projection_traits;
+
+  using size_type = typename AABB_tree::size_type;
+  using Primitive_id = typename AABB_tree::Primitive_id;
+  using Counting_iterator = CGAL::internal::AABB_tree::Counting_output_iterator<Primitive_id, size_type>;
+  using Listing_primitive_traits = typename AABB_traversal_traits::template Listing_primitive_traits<Line_3, Counting_iterator>;
+
+  CGAL::Real_timer timer;
+  timer.start();
+
+  for(int i=0; i<N; ++i)
+  {
+    Do_intersect_traits do_intersect_traversal_traits(tree.traits());
+    Line_3 query(CGAL::ORIGIN, points[i]);
+    tree.traversal(query, do_intersect_traversal_traits);
+    bool does_intersect = do_intersect_traversal_traits.is_intersection_found();
+    assert(does_intersect);
+
+    if(check_with_standard_tree)
+      assert(does_intersect == fg_tree.do_intersect(query));
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    Ray_3 query(CGAL::ORIGIN, points[i]);
+    First_intersection_traits first_intersection_traversal_traits(tree.traits());
+    tree.traversal(query, first_intersection_traversal_traits);
+    auto inter_res = first_intersection_traversal_traits.result();
+    assert(inter_res);
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    Ray_3 query(CGAL::ORIGIN, points[i]);
+    auto inter_res = tree.first_intersection(query);
+    assert(inter_res);
+
+    if(check_with_standard_tree)
+    {
+      auto fg_inter_res = fg_tree.first_intersection(query);
+      assert(inter_res->first == fg_inter_res->first);
+    }
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    auto hint = tree.best_hint(points[i]);
+    Projection_traits projection_traits(hint.first, hint.second, tree.traits());
+    Point res = tree.closest_point(points[i]);
+
+    if(check_with_standard_tree)
+    {
+      Point fg_res = fg_tree.closest_point(points[i]);
+      assert(res == fg_res);
+    }
+  }
+
+  for(int i=0; i<N; ++i)
+  {
+    Line_3 query(CGAL::ORIGIN, points[i]);
+
+    size_type counter = 0;
+    Counting_iterator out(&counter);
+    Listing_primitive_traits listing_traversal_traits(out, tree.traits());
+    tree.traversal(query, listing_traversal_traits);
+    assert(counter > 0);
+
+    if(check_with_standard_tree)
+    {
+      typename FG_Tree::size_type fg_counter = 0;
+      CGAL::internal::AABB_tree::Counting_output_iterator<typename FG_Tree::Primitive_id,
+                                                          typename FG_Tree::size_type> fg_out(&fg_counter);
+      fg_tree.all_intersected_primitives(query, fg_out);
+      assert(counter == fg_counter);
+    }
+  }
+
+  timer.stop();
+  std::cout << "Elapsed: " << timer.time() << " s." << std::endl;
+}
+
+template <typename K>
+void test(const char* filename)
+{
+  std::cout << " === Test Kernel " << typeid(K).name() << " === " << std::endl;
+
+  using Point_3 = typename K::Point_3;
+  using Mesh = CGAL::Surface_mesh<Point_3>;
+
+  Mesh mesh;
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh) ||
+     !is_triangle_mesh(mesh))
+  {
+    std::cerr << "Error: invalid input" << std::endl;
+    return;
+  }
+
+  test_no_cover(mesh);
+  test_cover(mesh, CGAL::parameters::default_values());
+}
+
+int main(int, char**)
+{
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  test<CGAL::Exact_predicates_inexact_constructions_kernel>("data/bunny00.off");
+  test<CGAL::Exact_predicates_exact_constructions_kernel>("data/bunny00.off");
+}
+

Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Advancing_front_surface_reconstruction_Examples)
 
 find_package(CGAL REQUIRED)

Advancing_front_surface_reconstruction/examples/Advancing_front_surface_reconstruction/reconstruction_surface_mesh.cpp

@@ -53,7 +53,6 @@ int main(int argc, char* argv[])
 {
   std::ifstream in((argc>1)?argv[1]:CGAL::data_file_path("points_3/half.xyz"));
   std::vector<Point_3> points;
-  std::vector<Facet> facets;
   Mesh m;
 
   std::copy(std::istream_iterator<Point_3>(in),

Advancing_front_surface_reconstruction/include/CGAL/Advancing_front_surface_reconstruction.h

@@ -394,7 +394,7 @@ namespace CGAL {
     bool deal_with_2d;
     Priority priority;
     int max_connected_component;
-    double K_init, K_step;
+    coord_type K_init, K_step;
     std::list<Vertex_handle> interior_edges;
     std::list< Incidence_request_elt > incidence_requests;
     typename std::list< Incidence_request_elt >::iterator sentinel;

Advancing_front_surface_reconstruction/test/Advancing_front_surface_reconstruction/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Advancing_front_surface_reconstruction_Tests)
 
 find_package(CGAL REQUIRED)

Algebraic_foundations/examples/Algebraic_foundations/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Algebraic_foundations_Examples)
 
 find_package(CGAL REQUIRED)

Algebraic_foundations/test/Algebraic_foundations/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Algebraic_foundations_Tests)
 
 find_package(CGAL REQUIRED COMPONENTS Core)

Algebraic_foundations/test/Algebraic_foundations/Chinese_remainder_traits.cpp

@@ -103,7 +103,7 @@ void test_CR_for(const NT& f){
     // try chinese remainder
     do{
         // chinese remainder failed if q > 2*max_coeff(f)
-         CGAL_postcondition_msg(pq < p*(2*max_coeff), " chinese remainder failed ");
+         assert(pq < p*(2*max_coeff));
 
         prime_index++;
         if(prime_index < 1000){

Algebraic_kernel_d/doc/Algebraic_kernel_d/CGAL/Algebraic_kernel_d_1.h

@@ -10,12 +10,10 @@ The template argument `Coeff` determines the coefficient type of the
 kernel, which is also the coefficient type of the supported polynomials.
 
 Currently, the following coefficient types are supported:
-
 - `Gmpz`, `Gmpq`, (requires configuration with external libraries GMP, MPFR and MPFI)
-
 - `CORE::BigInt`, `CORE::BigRat`, (requires configuration with external library GMP)
-
 - `leda_integer`, `leda_rational`. (requires configuration with external library LEDA)
+.
 
 \cgalAdvancedBegin
 The template argument type can also be set to `Sqrt_extension<NT,ROOT>`, where `NT` is

Algebraic_kernel_d/doc/Algebraic_kernel_d/CGAL/Algebraic_kernel_d_2.h

@@ -35,10 +35,10 @@ The template argument `Coeff` determines the coefficient type of the
 kernel, which is also the innermost coefficient type of the supported polynomials.
 
 Currently, the following coefficient types are supported:
-
 - `Gmpz`, `Gmpq`, (requires configuration with external libraries GMP, MPFR and MPFI)
 - `CORE::BigInt`, `CORE::BigRat`, (requires configuration with external library GMP)
 - `leda_integer`, `leda_rational`. (requires configuration with external library LEDA)
+.
 
 \cgalAdvancedBegin
 The template argument type can also be set to

Algebraic_kernel_d/examples/Algebraic_kernel_d/CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.1...3.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Algebraic_kernel_d_Examples)
 
 find_package(CGAL REQUIRED COMPONENTS Core)

Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/refine_zero_against.h

@@ -100,7 +100,7 @@ bool refine_zero_against(Field& low, Field& high, Polynomial p, Polynomial q) {
     if (CGAL::degree(q) == 0) return false;
 
     CGAL::Sign sign_p_low  = p.sign_at(low);
-    CGAL::Sign sign_p_high = p.sign_at(high);
+    CGAL_assertion_code(CGAL::Sign sign_p_high = p.sign_at(high));
     CGAL_precondition(sign_p_low  != CGAL::ZERO);
     CGAL_precondition(sign_p_high != CGAL::ZERO);
     CGAL_precondition(sign_p_high != sign_p_low);
@@ -148,9 +148,9 @@ bool refine_zero_against(Field& low, Field& high, Polynomial p, Polynomial q) {
             low = mid;
             sign_p_low = s;
         } else {
-            CGAL_postcondition(s == sign_p_high);
+            CGAL_assertion(s == sign_p_high);
             high = mid;
-            sign_p_high = s;
+            CGAL_assertion_code(sign_p_high = s);
         }
     }
 }
@@ -168,7 +168,7 @@ static bool strong_refine_zero_against(Field& low, Field& high,
     std::cout << "done, " << has_common_root << std::endl;
 
     CGAL::Sign sign_p_low = p.sign_at(low);
-    CGAL::Sign sign_p_high = p.sign_at(high);
+    CGAL_assertion_code(CGAL::Sign sign_p_high = p.sign_at(high));
 
     Field mid;
     CGAL::Sign s;
@@ -191,7 +191,7 @@ static bool strong_refine_zero_against(Field& low, Field& high,
         else {
             CGAL_assertion(s == sign_p_high);
             high = mid;
-            sign_p_high = s; //bogus?
+            CGAL_assertion_code(sign_p_high = s); //bogus?
         }
     }
 

Algebraic_kernel_d/test/Algebraic_kernel_d/CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.1...3.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Algebraic_kernel_d_Tests)
 
 # CGAL and its components

Algebraic_kernel_d/test/Algebraic_kernel_d/include/CGAL/_test_algebraic_kernel_1.h

@@ -409,9 +409,9 @@ void test_algebraic_kernel_1(const AlgebraicKernel_d_1& ak_1){
 #define CGAL_TEST_ALGEBRAIC_REAL_IO(_f)         \
     alg1=_f;                                    \
     ss<<CGAL::IO::oformat(alg1);                \
-    CGAL_assertion(ss.good());                  \
-    ss>>CGAL::IO::iformat(alg2);                    \
-    CGAL_assertion(!ss.fail());                 \
+    assert(ss.good());                          \
+    ss>>CGAL::IO::iformat(alg2);                \
+    assert(!ss.fail());                         \
     ss.clear();                                 \
     assert(alg1==alg2)
     // Note: after the reading ss>>CGAL::IO::iformat(alg2) the state of ss can

Algebraic_kernel_for_circles/test/Algebraic_kernel_for_circles/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Algebraic_kernel_for_circles_Tests)
 
 find_package(CGAL REQUIRED)

Algebraic_kernel_for_spheres/test/Algebraic_kernel_for_spheres/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Algebraic_kernel_for_spheres_Tests)
 
 find_package(CGAL REQUIRED)

Alpha_shapes_2/doc/Alpha_shapes_2/Alpha_shapes_2.txt

@@ -45,7 +45,7 @@ inside of the convex hull of \f$ S\f$. Hence, the alpha shape becomes
 the convex hull of \f$ S\f$ as \f$ \alpha \rightarrow \infty\f$.
 
 \cgal offers 2D and 3D alpha shapes. The GUDHI library offers a
-<a href="http://gudhi.gforge.inria.fr/doc/latest/group__alpha__complex.html"> dD Alpha complex</a>.
+<a href="https://gudhi.inria.fr/doc/latest/group__alpha__complex.html"> dD Alpha complex</a>.
 
 \section Alpha_shapes_2Definitions Definitions
 

Alpha_shapes_2/examples/Alpha_shapes_2/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Alpha_shapes_2_Examples)
 
 find_package(CGAL REQUIRED)

Alpha_shapes_2/examples/Alpha_shapes_2/ex_alpha_shapes_2.cpp

@@ -6,7 +6,6 @@
 #include <CGAL/Delaunay_triangulation_2.h>
 
 #include <CGAL/algorithm.h>
-#include <CGAL/assertions.h>
 
 #include <fstream>
 #include <iostream>

Alpha_shapes_2/test/Alpha_shapes_2/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Alpha_shapes_2_Tests)
 
 find_package(CGAL REQUIRED)

Alpha_shapes_3/demo/Alpha_shapes_3/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Alpha_shapes_3_Demo)
 
 # Find includes in corresponding build directories

Alpha_shapes_3/doc/Alpha_shapes_3/Alpha_shapes_3.txt

@@ -44,7 +44,7 @@ inside of the convex hull of \f$ S\f$. Hence, the alpha shape becomes
 the convex hull of \f$ S\f$ as \f$ \alpha \rightarrow \infty\f$.
 
 \cgal offers 2D and 3D alpha shapes. The GUDHI library offers a
-<a href="http://gudhi.gforge.inria.fr/doc/latest/group__alpha__complex.html"> dD Alpha complex</a>.
+<a href="https://gudhi.inria.fr/doc/latest/group__alpha__complex.html"> dD Alpha complex</a>.
 
 
 \section Alpha_shapes_3Definitions Definitions

Alpha_shapes_3/doc/Alpha_shapes_3/Doxyfile.in

@@ -2,5 +2,4 @@
 
 PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Alpha Shapes"
 
-EXAMPLE_PATH =  ${CGAL_Alpha_shapes_2_EXAMPLE_DIR} \
-                ${CGAL_Alpha_shapes_3_EXAMPLE_DIR}
+EXAMPLE_PATH +=  ${CGAL_Alpha_shapes_2_EXAMPLE_DIR}

Alpha_shapes_3/examples/Alpha_shapes_3/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Alpha_shapes_3_Examples)
 
 find_package(CGAL REQUIRED)

Alpha_shapes_3/examples/Alpha_shapes_3/visible_alpha_shape_facets_to_OFF.cpp

@@ -8,8 +8,8 @@
 #include <fstream>
 #include <vector>
 
-#include <boost/unordered_set.hpp>
-#include <boost/unordered_map.hpp>
+#include <unordered_set>
+#include <unordered_map>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel Gt;
 
@@ -48,7 +48,7 @@ int main()
 
 // collect alpha-shape facets accessible from the infinity
   // marks the cells that are in the same component as the infinite vertex by flooding
-  boost::unordered_set< Alpha_shape_3::Cell_handle > marked_cells;
+  std::unordered_set< Alpha_shape_3::Cell_handle > marked_cells;
   std::vector< Alpha_shape_3::Cell_handle > queue;
   queue.push_back( as.infinite_cell() );
 
@@ -86,7 +86,7 @@ int main()
 
 // dump into OFF format
   // assign an id per vertex
-  boost::unordered_map< Alpha_shape_3::Vertex_handle, std::size_t> vids;
+  std::unordered_map< Alpha_shape_3::Vertex_handle, std::size_t> vids;
   points.clear();
 
   for(Alpha_shape_3::Facet f : filtered_regular_facets)

Alpha_shapes_3/test/Alpha_shapes_3/CMakeLists.txt

@@ -1,7 +1,7 @@
 # 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.22)
+cmake_minimum_required(VERSION 3.1...3.23)
 project(Alpha_shapes_3_Tests)
 
 find_package(CGAL REQUIRED)

Alpha_wrap_3/doc/Alpha_wrap_3/Concepts/AlphaWrapOracle.h

@@ -0,0 +1,59 @@
+// Undocumented for now
+#ifndef DOXYGEN_RUNNING
+
+/*!
+\ingroup PkgAlphaWrap3Concepts
+\cgalConcept
+
+The concept `AlphaWrapOracle` defines the requirements for an Alpha Wrap <em>Oracle</em>, that is a class
+that answers a number of queries over the input of the algorithm.
+The oracle is the template parameter of the class `CGAL::Alpha_wraps_3_::Alpha_wrap_3`.
+
+\cgalHasModel `CGAL::Alpha_wraps_3_::Point_set_oracle`
+\cgalHasModel `CGAL::Alpha_wraps_3_::Segment_soup_oracle`
+\cgalHasModel `CGAL::Alpha_wraps_3_::Triangle_mesh_oracle`
+\cgalHasModel `CGAL::Alpha_wraps_3_::Triangle_soup_oracle`
+
+*/
+template <typename GeomTraits>
+class AlphaWrapOracle
+{
+public:
+  /// The geometric traits, must be a model of `AlphaWrapOracleTraits_3`
+  typedef unspecified_type Geom_traits;
+
+  /// Field type
+  typedef Geom_traits::FT FT;
+
+  /// Point type
+  typedef Geom_traits::Point_3 Point_3;
+
+  /// Sphere type
+  typedef Geom_traits::Ball_3 Ball_3;
+
+  /// Returns the geometric traits
+  Geom_traits geom_traits();
+
+  /// Returns an axis-aligned box enclosing the input data.
+  CGAL::Bbox_3 bbox();
+
+  /// Returns whether the ball `b` intersects the input data.
+  bool do_intersect(Ball_3 b);
+
+  /// Returns whether the tetrahedron `tet` intersects the input data.
+  template <typename K>
+  bool do_intersect(Tetrahedron_with_outside_info<K> tet);
+
+  /// Returns the intersection `o` of the segment `[p;q]` with the offset isolevel at distance `os`
+  /// from the input. In case of multiple intersections, the intersection closest to `p` is returned.
+  /// Returns `true` if there is an intersection, and `false` otherwise.
+  bool first_intersection(Point_3 p, Point_3 q, Point_3& o, double os);
+
+  /// Returns the point on the input data closest to `q`.
+  Point_3 closest_point(Point_3 q);
+
+  /// Returns the smallest squared distance between `q` and the input data.
+  FT squared_distance(Point_3 q);
+};
+
+#endif // DOXYGEN_RUNNING

Alpha_wrap_3/doc/Alpha_wrap_3/Concepts/AlphaWrapTraits_3.h

@@ -0,0 +1,188 @@
+// Undocumented for now
+#ifndef DOXYGEN_RUNNING
+
+/*!
+\ingroup PkgAlphaWrap3Concepts
+\cgalConcept
+
+@fixme Because of a few calls to PMP::, if you only look at the doc, it's all pointless because
+you require Kernel. Stitch_borders doesn't even have clear geometric traits requirements...
+
+The concept `AlphaWrapTraits_3` defines the requirements for the geometric traits class
+of an alpha wrap oracle.
+
+\cgalHasModel Any 3D %kernel is a model of this traits concept.
+*/
+
+class AlphaWrapTraits_3
+{
+public:
+  /// The field type
+  typedef unspecified_type FT;
+
+  /// The point type
+  typedef unspecified_type Point_3;
+
+  /// The vector type
+  typedef unspecified_type Vector_3;
+
+  /// The triangle type
+  typedef unspecified_type Triangle_3;
+
+  /// The tetrahedron type
+  typedef unspecified_type Tetrahedron_3;
+
+  /// The ball type
+  typedef unspecified_type Ball_3;
+
+  /*!
+  A predicate object that must provide the following function operator:
+
+  `CGAL::Comparison_result operator()(Point_3 p, Point_3 q, FT sqd)`,
+
+  which compares the squared distance between the two points `p` and `q` to the value `sqd`.
+  */
+  typedef unspecified_type Compare_squared_distance_3;
+
+  /*!
+  A construction object that must provide the following function operator:
+
+  `FT operator()(Vector_3 v)`,
+
+  which returns the squared length of the vector `v`.
+  */
+  typedef unspecified_type Compute_squared_length_3;
+
+  /*!
+  A construction object that must provide the following function operators:
+
+  `FT operator()(Point_3 p, Point_3 q, Point_3 r)`,
+
+  and
+
+  `FT operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s)`,
+
+  which return the squared radius of the smallest sphere enclosing the points.
+  */
+  typedef unspecified_type Compute_squared_radius_3;
+
+  /*!
+  A predicate object that must provide the following function operator:
+
+  `bool operator()(Point_3 p, Point_3 q, Point_3 r)`,
+
+  which returns `true`, iff `p`, `q`, and `r` are collinear.
+  */
+  typedef unspecified_type Collinear_3;
+
+  /*!
+  A construction object that must provide the following function operator:
+
+  `Ball_3 operator()(Point_3 p, FT sqr)`,
+
+  which returns the ball centered at `p` with squared radius `sqr`.
+  */
+  typedef unspecified_type Construct_ball_3;
+
+  /*!
+  A construction object that must provide the following function operators:
+
+  `CGAL::Bbox_3 operator()(Triangle_3 tr)`,
+
+  and
+
+  `CGAL::Bbox_3 operator()(Tetrahedron_3 tet)`,
+
+  which return an axis-aligned bounding box that encloses the object.
+  */
+  typedef unspecified_type Construct_bbox_3;
+
+  /*!
+  A construction object that must provide the following function operator:
+
+  `Vector_3 operator()(Vector_3 v, FT a)`,
+
+  which returns the vector `v` with length scaled by `a`.
+  */
+  typedef unspecified_type Construct_scaled_vector_3;
+
+  /*!
+  A construction object that must provide the following function operator:
+
+  `Point_3 operator()(Point_3 p, Vector_3 v)`,
+
+  which returns the point `p` translated by the vector `v`.
+  */
+  typedef unspecified_type Construct_translated_point_3;
+
+  /*!
+  A predicate object that must provide the following function operators:
+
+  `bool operator()(Tetrahedron_3 tet, Point_3 p)`,
+
+  and
+
+  `bool operator()(Sphere_3 s, Point_3 p)`,
+
+  which return `true` iff `p` is on the bounded side of the respective kernel objects.
+  */
+  typedef unspecified_type Has_on_bounded_side_3;
+
+  /*!
+  A predicate object that must provide the following function operator:
+
+  `bool operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s)`,
+
+  which returns `true` iff `s` is on the bounded side of the smallest sphere enclosing `p`, `q`, and `r`.
+  */
+  typedef unspecified_type Side_of_bounded_sphere_3;
+
+  // ===
+
+  /*!
+  returns the `Compare_squared_distance_3` predicate.
+  */
+  Compare_squared_distance_3 compare_squared_distance_3_object();
+
+  /*!
+  returns the `Compute_squared_length_3` predicate.
+  */
+  Compute_squared_length_3 compute_squared_length_3_object();
+
+  /*!
+  returns the `Compute_squared_radius_3` predicate.
+  */
+  Compute_squared_radius_3 compute_squared_radius_3_object();
+
+  /*!
+  returns the `Collinear_3` predicate.
+  */
+  Collinear_3 collinear_3_object();
+
+  /*!
+  returns the `Construct_ball_3` construction.
+  */
+  Construct_ball_3 construct_ball_3_object();
+
+  /*!
+  returns the `Construct_scaled_vector_3` construction.
+  */
+  Construct_scaled_vector_3 construct_scaled_vector_3_object();
+
+  /*!
+  returns the `Construct_translated_point_3` construction.
+  */
+  Construct_translated_point_3 construct_translated_point_3_object();
+
+  /*!
+  returns the `Has_on_bounded_side_3` predicate.
+  */
+  Has_on_bounded_side_3 has_on_bounded_side_3_object();
+
+  /*!
+  returns the `Side_of_bounded_sphere_3` predicate.
+  */
+  Side_of_bounded_sphere_3 side_of_bounded_sphere_3_object();
+};
+
+#endif // DOXYGEN_RUNNING

Alpha_wrap_3/doc/Alpha_wrap_3/Doxyfile.in

@@ -0,0 +1,20 @@
+@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
+
+PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - 3D Alpha Wrapping"
+
+# custom options for this package
+EXTRACT_ALL                =  false
+HIDE_UNDOC_MEMBERS         =  true
+HIDE_UNDOC_CLASSES         =  true
+
+HTML_EXTRA_FILES       = ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_bike_lod.jpg \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_triangle_soup.jpg \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_non_manifold_cases.jpg \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_pencil.png \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_steps.jpg \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_church_lod.jpg \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_sharp_feature.jpg \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_steiner.jpg \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_alpha_offset_bike.jpg \
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_double_sided.jpg \ 
+                         ${CGAL_PACKAGE_DOC_DIR}/fig/aw3_thingi10k_benchmark.jpg

Alpha_wrap_3/doc/Alpha_wrap_3/PackageDescription.txt

@@ -0,0 +1,35 @@
+/// \defgroup PkgAlphaWrap3Ref 3D Alpha Wrapping
+
+/// \defgroup AW3_free_functions_grp Free Functions
+/// Functions to create a wrap from point clouds, triangle soups, and triangle meshes.
+/// \ingroup PkgAlphaWrap3Ref
+
+/*!
+\addtogroup PkgAlphaWrap3Ref
+\cgalPkgDescriptionBegin{3D Alpha Wrapping,PkgAlphaWrap3}
+\cgalPkgPicture{alpha_wrap_3.jpg}
+\cgalPkgSummaryBegin
+\cgalPkgAuthors{Pierre Alliez, David Cohen-Steiner, Michael Hemmer, Cédric Portaneri, Mael Rouxel-Labbé}
+\cgalPkgDesc{This component takes a 3D triangle mesh, a triangle soup, or a point set as input, and generates
+             a valid triangulated surface mesh that strictly contains the input (watertight, intersection-free
+             and 2-manifold). The algorithm proceeds by shrink-wrapping
+             and refining a 3D Delaunay triangulation starting from a loose bounding box of the input.
+             Two user-defined parameters, alpha and offset, offer control over the maximum size
+             of cavities where the shrink-wrapping process can enter, and the tightness
+             of the final surface mesh to the input, respectively. Once combined, these parameters
+             provide a means to trade fidelity to the input for complexity of the output.}
+\cgalPkgManuals{Chapter_3D_Alpha_wrapping,PkgAlphaWrap3Ref}
+\cgalPkgSummaryEnd
+\cgalPkgShortInfoBegin
+\cgalPkgSince{5.5}
+\cgalPkgDependsOn{\ref PkgTriangulation3 and \ref PkgPolygonMeshProcessing}
+\cgalPkgBib{cgal:achpr-aw3}
+\cgalPkgLicense{\ref licensesGPL "GPL"}
+\cgalPkgShortInfoEnd
+\cgalPkgDescriptionEnd
+
+\cgalClassifedRefPages
+
+\cgalCRPSection{Functions}
+- \link AW3_free_functions_grp `CGAL::alpha_wrap_3()` \endlink
+*/

Alpha_wrap_3/doc/Alpha_wrap_3/alpha_wrap_3.txt

@@ -0,0 +1,368 @@
+namespace CGAL {
+
+/*!
+
+\mainpage User Manual
+\anchor Chapter_3D_Alpha_wrapping
+\cgalAutoToc
+
+\authors Pierre Alliez, David Cohen-Steiner, Michael Hemmer, Cédric Portaneri, Mael Rouxel-Labbé
+
+<center>
+<img src="aw3_bike_lod.jpg" style="max-width:70%;"/>
+</center>
+
+\section aw3_introduction Introduction
+
+Various tasks in geometric modeling and processing require 3D objects represented as valid surface meshes,
+where "valid" refers to meshes that are watertight, intersection-free, orientable, and 2-manifold.
+Such representations offer well-defined notions of interior/exterior and geodesic neighborhoods.
+
+3D data are usually acquired through measurements followed by reconstruction, designed by humans,
+or generated through imperfect automated processes.
+As a result, they can exhibit a wide variety of defects including gaps, missing data,
+self-intersections, degeneracies such as zero-volume structures, and non-manifold features.
+
+Given the large repertoire of possible defects, many methods and data structures have been proposed
+to repair specific defects, usually with the goal of guaranteeing specific properties in the repaired 3D model.
+Reliably repairing all types of defects is notoriously difficult and is often an ill-posed problem
+as many valid solutions exist for a given 3D model with defects.
+In addition, the input model can be overly complex with unnecessary geometric details,
+spurious topological structures, nonessential inner components, or excessively fine discretizations.
+For applications such as collision avoidance, path planning, or simulation,
+getting an approximation of the input can be more relevant than repairing it.
+Approximation herein refers to an approach capable of filtering out inner structures,
+fine details and cavities, as well as wrapping the input within a user-defined offset margin.
+
+Given an input 3D geometry, we address the problem of computing a conservative approximation,
+where conservative means that the output is guaranteed to strictly enclose the input.
+We seek unconditional robustness in the sense that the output mesh should be valid (oriented,
+2-manifold, and without self-intersections), even for raw input with many defects
+and degeneracies.
+The default input is a soup of 3D triangles, but the generic interface leaves the door open
+to other types of finite 3D primitives such as triangle soups and point sets.
+
+\cgalFigureAnchor{1}
+<center>
+<img src="aw3_triangle_soup.jpg" style="max-width:70%;"/>
+</center>
+\cgalFigureCaptionBegin{1}
+Shrink-wrapping output from a triangle soup, with many intersections and gaps.
+From left to right, input model, output wrap, and superposition.
+\cgalFigureCaptionEnd
+
+\cgalFigureAnchor{2}
+<center>
+<img src="aw3_non_manifold_cases.jpg" style="max-width:75%;"/>
+</center>
+\cgalFigureCaptionBegin{2}
+Dealing with non-manifold features and degeneracies.
+From left to right, a non-manifold vertex, self-intersecting faces and two adjacent triangles
+representing a zero-volume structure.
+The algorithm handles these cases by wrapping an offset of the input.
+\cgalFigureCaptionEnd
+
+\section aw3_definition Approach
+
+Many approaches have been devised to enclose a 3D model within a volume, featuring different balances
+between the runtime and quality (i.e., tightness) of the approximation.
+Within the simplest cases, an axis-aligned or oriented bounding box clearly satisfies some desired properties;
+however, the approximation error is uncontrollable and often very large.
+Computing the convex hull of the input also matches some of the desired properties
+and improves the quality of the result, albeit at the price of increasing the runtime.
+However, the approximation remains crude, especially in the case of several components.
+
+The convex hull is, in fact, a special case of alpha shapes (\ref Chapter_3D_Alpha_Shapes).
+Mathematically, the alpha shape is a subcomplex of the Delaunay triangulation, with simplicies
+being part of the complex depending on the size of their minimal (empty) Delaunay ball.
+Intuitively, constructing 3D alpha shapes can be thought of as carving 3D space with an empty ball
+of user-defined radius alpha.
+Alpha shapes yield provable, good piecewise-linear approximations of a shape \cgalCite{bb-srmua-97t},
+but are defined on point sets, whereas we wish to deal with more general input data, such as triangle soups.
+Even after sampling the triangle soup, alpha shapes do not guarantee to be conservative for any alpha.
+Finally, inner structures are also carved within the volumes, instead of being filtered out.
+
+Inspired by alpha shapes, we replace the above notion of carving by <em>shrink-wrapping</em>:
+we iteratively construct a subcomplex of a 3D Delaunay triangulation by starting from
+a simple 3D Delaunay triangulation enclosing the input, and then iteratively removing eligible tetrahedra
+that lie on the boundary of the complex.
+In addition, the underlying triangulation---and thus the complex incidentally---is refined
+as shrinking proceeds.
+Thus, instead of carving from the convex hull of the input data as in alpha shapes, we construct
+an entirely new mesh through a Delaunay refinement-like algorithm. The refinement algorithm inserts Steiner points
+on the boundary of an offset volume, defined as a level set of the unsigned distance field to the input.
+
+This process both prevents the creation of inner structures within the output and avoids superfluous computations.
+In addition, detaching our mesh construction from the geometry and discretization of the input has several advantages:
+(1) the underlying data is not restricted to a specific format (triangle soups, polygon soups, point clouds, etc.)
+as all that is required is answering three basic geometric queries: (a) the distance between a point
+and the input, (b) the projection of a query point onto the input, (c) an intersection test
+between a tetrahedron and the input, and (2) The user has more freedom to trade tightness
+to the input for final mesh complexity, as constructing a conservative approximation on a large offset
+of the input requires fewer mesh elements.
+
+\subsection aw3_algorithm Algorithm
+
+<b>Initialization</b>. The algorithm is initialized by inserting the eight corner vertices
+of a loose bounding box into a 3D Delaunay triangulation.
+In the 3D Delaunay triangulation of \cgal, all triangle facets are adjacent to two tetrahedron cells.
+Each facet of the boundary of the Delaunay triangulation, which coincides with one facet
+of the convex hull of the triangulation vertices, is adjacent to a so-called <em>infinite</em> tetrahedron cell,
+an abstract cell connected to the so-called <em>infinite vertex</em> to ensure the aforementioned double-facet adjacency.
+Initially, all infinite cells are tagged as outside, and all finite tetrahedron cells are tagged as inside.
+
+<b>Shrink-wrapping</b>. The shrink-wrapping algorithm proceeds by traversing the cells
+of the Delaunay triangulation from outside to inside, flood-filling from one cell to its adjacent cell,
+and tagging the adjacent cell as outside whenever possible (the term possible is specified later).
+Flood filling is implemented via a priority queue of Delaunay triangle facets representing
+the traversal between the two adjacent cells of the facet, from outside to inside.
+These triangle facets are referred to as <em>gates</em> in the following.
+
+Given an outside cell and its adjacent inside cell, the common facet (i.e., a gate) is said
+to be <em>alpha-traversable</em> if its circumradius is larger than the user-defined parameter alpha,
+where circumradius refers to the radius of the relating triangle's Delaunay ball.
+Intuitively, cavities smaller than alpha are not accessible as their gates are not alpha-traversable.
+
+Initialized by the alpha-traversable gates on the convex hull, the priority queue contains only
+alpha-traversable gates and is sorted by decreasing order
+of the circumradius of the gate.
+Traversal can be seen as a continuous process that advances along dual Voronoi edges of the gates,
+with a pencil of empty balls circumscribing the gate.
+
+\cgalFigureAnchor{3}
+<center>
+<img src="aw3_pencil.png" style="max-width:40%;"/>
+</center>
+\cgalFigureCaptionBegin{3}
+(Left) Pencil of empty circles (blue) circumscribing a Delaunay edge (green) in a 2D Delaunay triangulation (black).
+From the top triangle circumcenter <em>c1</em> to the bottom triangle circumcenter <em>c2</em>, the dual Voronoi edge denoted by <em>e</em> (doted red) is the trace of centers of the largest circles that are empty of Delaunay vertex.
+(Right) The graph corresponding to the left example. The x-axis corresponds to the position of empty circle centers located on the Voronoi edge <em>e</em>, from <em>c1</em> to <em>c2</em>. The y-axis is the radius value of the corresponding empty circles. In this case, the minimum radius of this pencil of empty circle is located at the midpoint of the green Delaunay edge.
+In our algorithm, a gate (green Delaunay edge) is said to be not alpha-traversable when the minimum radius of the pencil of empty circle is smaller than alpha.
+\cgalFigureCaptionEnd
+
+When traversing from an outside cell \f$ c_o \f$ to an inside cell \f$ c_i \f$ through an alpha-traversable facet \f$ f \f$,
+two criteria are tested to prevent the wrapping process from colliding with the input:
+
+(1) We check for an intersection between the dual Voronoi edge of \f$ f \f$, i.e. the segment between
+the circumcenters of the two incident cells, and the <em>offset surface</em>, defined as the level set
+of unsigned isosurface to the input.
+If one or several intersections exists, the first intersection point, along the dual Voronoi edge
+oriented from outside to inside is inserted into the triangulation as a Steiner point.
+
+(2) If the dual Voronoi edge does not intersect the offset surface but the neighboring cell \f$ c_i \f$
+intersects the input, we compute the projection of the circumcenter of \f$ c_i \f$
+onto the offset surface, and insert it into the triangulation as a Steiner point (which destroys \f$ c_i \f$).
+
+After each of the above Steiner point insertions, all new incident cells are tagged as inside,
+and the newly alpha-traversable gates are pushed into the priority queue.
+
+If none of the above two criteria are met, the neighboring cell \f$ c_i \f$ is traversed and tagged as outside.
+Alpha-Traversable facets of \f$ c_i \f$ that are separating inside from outside cells are pushed as new gates into the priority queue.
+
+Once the queue empties---a process that is guaranteed as facets (and their circumradii) become smaller
+due to the insertion of new Steiner points---the construction phase terminates.
+The output triangle surface mesh is extracted from the Delaunay triangulation as the set of facets
+separating inside from outside cells.
+
+The figure below depicts the steps of the algorithm in 2D.
+
+\cgalFigureAnchor{4}
+<center>
+<img src="aw3_steps.jpg" style="max-width:95%;"/>
+</center>
+\cgalFigureCaptionBegin{4}
+Steps of the shrink-wrapping algorithm in 2D.
+The algorithm is initialized by inserting the corners of the loose bounding box of the input (red)
+into a Delaunay triangulation, and all finite triangles are tagged inside (grey).
+The current gate (green edge) popped out from the queue is alpha-traversable. The triangle adjacent
+to the gate is tagged outside when it does not intersect the input, and new alpha-traversable gates
+are pushed to the queue. When the adjacent triangle intersects the input, a new Steiner point (large green disc)
+is computed and inserted into the triangulation, all neighboring triangles are tagged inside,
+new alpha-traversable gates are pushed to the queue, and traversal is resumed.
+Grey edges depict the Delaunay triangulation. Blue edges depict the Voronoi diagram.
+Pink circles depict the empty circle of radius alpha. The output edges (dark blue) separate inside from outside triangles.
+\cgalFigureCaptionEnd
+
+\subsection aw3_guarantees Guarantees
+
+The algorithm is proven to terminate and to produce a 2-manifold triangulated surface mesh
+that strictly encloses the input data.
+The key element to the proof is that we wrap from outside to inside and never allow a cell
+that intersects the input to be flagged inside.
+Furthermore, both criteria that lead to refinement of the triangulation insert Steiner points
+that are guaranteed to break the cells in need of refinement and reduce the neighbor facets circumradii.
+
+Because the main refinement criterion is the insertion of an intersection between a dual Voronoi edge
+and an offset of the input, or the projection of a Voronoi vertex onto the offset of the input,
+the algorithm has similarities to popular meshing algorithms based on Delaunay filtering
+and refinement (see \ref Chapter_3D_Mesh_Generation).
+
+\section aw3_interface Interface
+
+Our algorithm takes as input a set of triangles in 3D, provided either as a triangle soup or
+as a triangle surface mesh, and two user-defined scalar parameters: the <em>alpha</em> and the <em>offset</em> values.
+It proceeds by shrink-wrapping and refining a 3D Delaunay triangulation starting from a loose bounding box of the input.
+The parameter <em>alpha</em> refers to the size of cavities or holes that cannot be traversed during wrapping,
+and hence to the final level of detail, as alpha acts like a sizing field in a common Delaunay
+refinement algorithm (\ref Chapter_3D_Mesh_Generation).
+The parameter <em>offset</em> refers to the distance between the vertices of the refined triangulation
+and the input, so that a large offset translates into a loose enclosing of the input.
+This second parameter offers a means to control the trade-off between tightness and complexity.
+
+The main entry point of the component is the global function `CGAL::alpha_wrap_3()` that generates the alpha wrap;
+this function takes as input a polygon soup or a polygon mesh.
+There is no prerequisite on the input connectivity so that it can take an arbitrary triangle soup,
+with islands, self-intersections, or overlaps, as well as combinatorial or geometrical degeneracies.
+
+The underlying traits class must be a model of the `Kernel` concept. It should use
+a floating point number type as inexactness is inherent to the algorithm since there is no closed
+form description of new vertices on the offset surface.
+
+The output is a triangle surface mesh whose type is chosen by the user, under the constraint
+that it must be a model of the `MutableFaceGraph` concept.
+
+\section aw3_parameters Choosing Parameters
+
+The two parameters of the algorithm impact both the level of detail and complexity of the output mesh.
+
+\subsection aw3_alpha Alpha
+
+The main parameter, alpha, controls whether a Delaunay facet is traversable during shrink-wrapping.
+Alpha's main purpose is to control the size
+of the empty balls used during wrapping, and thus to determine which features will appear in the output:
+indeed, a facet is alpha-traversable if its circumradius is larger than alpha; hence, the algorithm
+can only shrink-wrap through straits or holes with diameters larger than alpha.
+A second, less direct consequence is that as long as a facet has a circumradius larger than alpha,
+the incident inside cell will be visited and possibly refined.
+Therefore, when the algorithm terminates, all facets have a circumradius smaller than alpha.
+This parameter thus also behaves like a sizing criterion on the triangle facets of the output.
+
+\cgalFigureAnchor{5}
+<center>
+<img src="aw3_church_lod.jpg" style="max-width:90%;"/>
+</center>
+\cgalFigureCaptionBegin{5}
+Impact of the alpha parameter on the output.
+(Left) The input triangle mesh, generated by surface reconstruction from a raw point cloud,
+has many non-manifold edges and vertices, superfluous geometric details and spurious topological structures.
+(Right) This component approximates the input conservatively and produces valid meshes with different
+complexity and fidelity to the input, depending on the alpha parameter.
+The smaller the alpha, the deeper the shrink-wrapping process will enter cavities.
+The alpha parameter is decreasing from left to right, to respectively 1/50, 1/100 and 1/300 of the longest diagonal of the input bounding box.
+A large alpha will produce an output less complex but less faithful to the input.
+\cgalFigureCaptionEnd
+
+\subsection aw3_offset Offset
+
+The second parameter, the offset distance, controls the distance from the input and thus the definition
+of the offset isosurface onto which the vertices of the output mesh are located.
+This parameter controls the tightness of the result, which has, in turn, a few consequences.
+Firstly, locating vertices away from the input enables the algorithm to generate
+a less complex mesh, especially in convex areas. A trivial example of this behavior would be a very dense mesh
+of a sphere, for which an as-tight-as-possible envelope would also be very dense.
+Secondly, the farther the isosurface is from the input, the more new points are inserted
+through the first criterion (i.e., through intersection with dual Voronoi edge, see Section \ref aw3_algorithm);
+thus, the quality of the output improves in terms of angles of the triangle elements.
+Finally, and depending on the value of the alpha parameter, a large offset can also offer defeaturing capabilities.
+However using a small offset parameter will tend to better preserve sharp features as projection
+Steiner points tend to project onto convex sharp features.
+
+\cgalFigureAnchor{6}
+<center>
+<img src="aw3_sharp_feature.jpg" style="max-width:90%;"/>
+</center>
+\cgalFigureCaptionBegin{6}
+Impact of the offset parameter on the output.
+(Left) Input mesh generated by meshing a NURBS CAD model in parameter space.
+(Right) The smaller the offset, the closest sample points are to the input.
+The offset parameter is decreasing from left to right, to respectively 1/50, 1/200 and 1/1000 of the longest diagonal of the input bounding box.
+The alpha parameter is equal to 1/50 of the longest diagonal of the input bounding box for all level of details.
+A larger offset will produce an output less complex with better triangle quality.
+However the sharp features (red edges) are well preserved when the offset parameter is small.
+\cgalFigureCaptionEnd
+
+\cgalFigureAnchor{7}
+<center>
+<img src="aw3_steiner.jpg" style="max-width:90%;"/>
+</center>
+\cgalFigureCaptionBegin{7}
+Steiner points.
+The projection Steiner points (green) are computed by projecting the triangle circumcenter onto the offset.
+The intersection Steiner points (blue) are computed as the first intersection point between the Voronoi edge and the offset.
+(Left) When the offset parameter is small, the algorithm produces more projection Steiner points,
+which tends to improve the preservation of convex sharp features.
+(Right) When the offset parameter is large, the algorithm produces more intersection Steiner points,
+which tends to generate triangles with better quality in terms of angles, in 3D.
+\cgalFigureCaptionEnd
+
+By default, we recommend to set the offset parameter to a small fraction of alpha, so that alpha
+becomes the main parameter that controls the final level of detail.
+
+The image below illustrates the impact of both parameters.
+
+\cgalFigureAnchor{8}
+<center>
+<img src="aw3_alpha_offset_bike.jpg" style="max-width:80%;"/>
+</center>
+\cgalFigureCaptionBegin{8}
+Different alpha and offset values on the bike model (533,000 triangles).
+The x-axis represents the offset value equal to 1/5000, 1/2000, 1/500, 1/200, 1/50, 1/20 and 1/5 of the longest diagonal of the input bounding box, from left to right.
+The y-axis represents the alpha value equal to 1/300, 1/100, 1/50, 1/20 and 1/5 of the longest diagonal of the input bounding box, from bottom to top.
+The numbers below each level of detail represents their number of triangles.
+Depending on the alpha value, an offset too small or too large will produce output mesh with higher complexity.
+For each alpha, the models with lower complexity can be used as a scale-space representations for collision detection, from near to far distances.
+\cgalFigureCaptionEnd
+
+\subsection aw3_two_side A Note on "Two-Sided" Wraps
+
+The offset parameter is crucial to our approach because it guarantees that the output is a closed,
+2-manifold surface mesh.
+Indeed, and even when the input is a zero-volume structure such as a single 3D triangle,
+the output wrap is a thin volume enclosing the said triangle \cgalFigureRef{2}.
+
+Users should keep in mind that the wrapping algorithm has no means of determining whether it is acting on the inside or the outside
+of the unsigned distance field, and will thus produce two-sided wraps in the case of holes in the input
+and values of alpha smaller than the size of the holes.
+
+\cgalFigureAnchor{9}
+<center>
+<img src="aw3_double_sided.jpg" style="max-width:80%;"/>
+</center>
+\cgalFigureCaptionBegin{9}
+Two-sided wrap.
+(Left) Wrapping a Bunny in 2D, with decreasing values for alpha.
+(Right) Wrapping a defect-laden Bunny in 3D. The rightmost column depicts a clipped visualization
+of the inside. When alpha is small enough with respect the diamater of the holes, the algorithm generates a two-sided wrap.
+\cgalFigureCaptionEnd
+
+\section aw3_performance Performance
+
+The charts below plots the computation times of the wrapping algorithm on the Thingi10k dataset, as well as the complexity of the output triangle mesh.
+
+\cgalFigureAnchor{10}
+<center>
+<img src="aw3_thingi10k_benchmark.jpg" style="max-width:80%;"/>
+</center>
+\cgalFigureCaptionBegin{9}
+Execution times and output complexity for different values of alpha on the Thingi10k data set.
+Alpha increases from 1/20 to 1/200 of the length of bounding box diagonal.
+The x axis represents the complexity of the output wrap mesh in number of triangle facets.
+The y axis represents the total computation time, in seconds.
+The color and diameter of the dots represent the number of faces in the input triangle soup,
+ranging from 10 (green) to 3154000 (blue).
+\cgalFigureCaptionEnd
+
+\section aw3_examples Examples
+
+Here is an example with an input triangle mesh, with alpha set to 1/20 of the bounding box longest diagonal edge length,
+and offset set to 1/30 of alpha (i.e., 1/600 of the bounding box diagonal edge length).
+
+\cgalExample{Alpha_wrap_3/triangle_mesh_wrap.cpp}
+
+Here is an example with a point cloud.
+
+\cgalExample{Alpha_wrap_3/point_set_wrap.cpp}
+
+*/
+
+}

Alpha_wrap_3/doc/Alpha_wrap_3/dependencies

@@ -0,0 +1,11 @@
+Manual
+Kernel_23
+STL_Extension
+Algebraic_foundations
+Circulator
+Stream_support
+AABB_tree
+Alpha_shapes_3
+BGL
+Mesh_3
+Triangulation_3

Alpha_wrap_3/doc/Alpha_wrap_3/examples.txt

@@ -0,0 +1,4 @@
+/*!
+\example Alpha_wrap_3/triangle_mesh_wrap.cpp
+\example Alpha_wrap_3/point_set_wrap.cpp
+*/

Alpha_wrap_3/examples/Alpha_wrap_3/CMakeLists.txt

@@ -0,0 +1,13 @@
+# 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.23)
+project(Alpha_wrap_3_Examples)
+
+find_package(CGAL REQUIRED)
+
+# create a target per cppfile
+create_single_source_cgal_program("triangle_mesh_wrap.cpp")
+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")

Alpha_wrap_3/examples/Alpha_wrap_3/mixed_inputs_wrap.cpp

@@ -0,0 +1,131 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/alpha_wrap_3.h>
+#include <CGAL/Real_timer.h>
+#include <CGAL/IO/read_points.h>
+#include <CGAL/IO/polygon_soup_io.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 Segment_3 = K::Segment_3;
+
+using Face = std::array<std::size_t, 3>;
+using Segments = std::vector<Segment_3>;
+using Points = std::vector<Point_3>;
+using Faces = std::vector<Face>;
+
+using Mesh = CGAL::Surface_mesh<Point_3>;
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+
+  // Read the inputs
+  const std::string ts_filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off"); // triangle soup
+  const std::string ss_filename = (argc > 2) ? argv[2] : CGAL::data_file_path("images/420.polylines.txt"); // segment soup
+  const std::string ps_filename = (argc > 3) ? argv[3] : CGAL::data_file_path("points_3/ball.ply"); // point set
+
+  std::cout << "Triangle soup input: " << ts_filename << std::endl;
+  std::cout << "Segment soup input: " << ss_filename << std::endl;
+  std::cout << "Point set input: " << ps_filename << std::endl;
+
+  // = read the soup
+  Points points;
+  Faces faces;
+  if(!CGAL::IO::read_polygon_soup(ts_filename, points, faces) || points.empty() || faces.empty())
+  {
+    std::cerr << "Invalid soup input: " << ts_filename << std::endl;
+    return EXIT_FAILURE;
+  }
+  std::cout << points.size() << " points (triangle soup)" << std::endl;
+
+  // = read the polylines
+  std::ifstream ifs(ss_filename);
+  Segments segments;
+  int len = 0;
+  while(ifs >> len)
+  {
+    std::vector<Point_3> polyline;
+    while(len--)
+    {
+      Point_3 point;
+      ifs >> point;
+      polyline.push_back(point);
+    }
+
+    if(polyline.size() >= 2)
+    {
+      for(std::size_t i=0; i<polyline.size() - 1; ++i)
+        segments.emplace_back(polyline[i], polyline[i+1]);
+    }
+  }
+  std::cout << segments.size() << " segments (segment soup)" << std::endl;
+
+  // = read the points
+  Points ps_points;
+  if(!CGAL::IO::read_points(ps_filename, std::back_inserter(ps_points)))
+  {
+    std::cerr << "Invalid point set input: " << ps_filename << std::endl;
+    return EXIT_FAILURE;
+  }
+  std::cout << ps_points.size() << " points (point set)" << std::endl;
+
+  const double relative_alpha = (argc > 4) ? std::stod(argv[4]) : 15.;
+  const double relative_offset = (argc > 5) ? std::stod(argv[5]) : 450.;
+
+  CGAL::Bbox_3 bbox = bbox_3(std::cbegin(points), std::cend(points));
+  CGAL::Bbox_3 ps_bbox = bbox_3(std::cbegin(ps_points), std::cend(ps_points));
+  bbox += ps_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;
+
+  CGAL::Real_timer t;
+  t.start();
+
+  using TS_Oracle = CGAL::Alpha_wraps_3::internal::Triangle_soup_oracle<K>;
+  using SS_Oracle = CGAL::Alpha_wraps_3::internal::Segment_soup_oracle<K, TS_Oracle>;
+  using Oracle = CGAL::Alpha_wraps_3::internal::Point_set_oracle<K, SS_Oracle>;
+
+  TS_Oracle ts_oracle(K{});
+  SS_Oracle ss_oracle(ts_oracle);
+  Oracle oracle(ss_oracle);
+
+  oracle.add_triangle_soup(points, faces, CGAL::parameters::default_values());
+  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);
+
+  Mesh output_mesh;
+  aw3(alpha, offset, output_mesh);
+
+  t.stop();
+  std::cout << "Took " << t.time() << std::endl;
+
+  std::string ts_name = std::string(ts_filename);
+  ts_name = ts_name.substr(ts_name.find_last_of("/") + 1, ts_name.length() - 1);
+  ts_name = ts_name.substr(0, ts_name.find_last_of("."));
+  std::string ss_name = std::string(ss_filename);
+  ss_name = ss_name.substr(ss_name.find_last_of("/") + 1, ss_name.length() - 1);
+  ss_name = ss_name.substr(0, ss_name.find_last_of("."));
+  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::cout << "Writing to " << output_name << std::endl;
+  CGAL::IO::write_polygon_mesh(output_name, output_mesh, CGAL::parameters::stream_precision(17));
+
+  return EXIT_SUCCESS;
+}

Alpha_wrap_3/examples/Alpha_wrap_3/point_set_wrap.cpp

@@ -0,0 +1,69 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/alpha_wrap_3.h>
+#include <CGAL/IO/read_points.h>
+#include <CGAL/Real_timer.h>
+
+#include <iostream>
+#include <string>
+
+namespace AW3 = CGAL::Alpha_wraps_3;
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point_3 = K::Point_3;
+
+using Point_container = std::vector<Point_3>;
+using Mesh = CGAL::Surface_mesh<Point_3>;
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+
+  // Read the input
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/oni.pwn");
+  std::cout << "Reading " << filename << "..." << std::endl;
+
+  Point_container points;
+  if(!CGAL::IO::read_points(filename, std::back_inserter(points)) || points.empty())
+  {
+    std::cerr << "Invalid input." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << points.size() << " points" << std::endl;
+
+  // Compute the alpha and offset values
+  const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 10.;
+  const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 300.;
+
+  CGAL::Bbox_3 bbox = CGAL::bbox_3(std::cbegin(points), std::cend(points));
+  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;
+  std::cout << "absolute alpha = " << alpha << " absolute offset = " << offset << std::endl;
+
+  // Construct the wrap
+  CGAL::Real_timer t;
+  t.start();
+
+  Mesh wrap;
+  CGAL::alpha_wrap_3(points, alpha, offset, wrap);
+
+  t.stop();
+  std::cout << "Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;
+  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";
+  std::cout << "Writing to " << output_name << std::endl;
+  CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
+
+  return EXIT_SUCCESS;
+}

Alpha_wrap_3/examples/Alpha_wrap_3/triangle_mesh_wrap.cpp

@@ -0,0 +1,71 @@
+#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 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 Mesh = CGAL::Surface_mesh<Point_3>;
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+
+  // Read the input
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.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." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  std::cout << "Input: " << num_vertices(mesh) << " vertices, " << num_faces(mesh) << " faces" << std::endl;
+
+  // Compute the alpha and offset values
+  const double relative_alpha = (argc > 2) ? std::stod(argv[2]) : 20.;
+  const double relative_offset = (argc > 3) ? std::stod(argv[3]) : 600.;
+
+  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()));
+
+  const double alpha = diag_length / relative_alpha;
+  const double offset = diag_length / relative_offset;
+
+  // Construct the wrap
+  CGAL::Real_timer t;
+  t.start();
+
+  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 << "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";
+  std::cout << "Writing to " << output_name << std::endl;
+  CGAL::IO::write_polygon_mesh(output_name, wrap, CGAL::parameters::stream_precision(17));
+
+  return EXIT_SUCCESS;
+}

Alpha_wrap_3/examples/Alpha_wrap_3/wrap_from_cavity.cpp

@@ -0,0 +1,85 @@
+#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 <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 Mesh = CGAL::Surface_mesh<Point_3>;
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+
+  // Read the input
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/armadillo.off");
+  std::cout << "Reading " << filename << "..." << std::endl;
+
+  Mesh input;
+  if(!PMP::IO::read_polygon_mesh(filename, input) ||
+     is_empty(input) || !is_triangle_mesh(input))
+  {
+    std::cerr << "Invalid input." << 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_offset = (argc > 3) ? std::stod(argv[3]) : 600.;
+
+  // Compute the alpha and offset values
+  CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input);
+  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;
+
+  // Construct the wrap
+  using Oracle = CGAL::Alpha_wraps_3::internal::Triangle_mesh_oracle<K>;
+  Oracle oracle;
+  oracle.add_triangle_mesh(input);
+
+  CGAL::Real_timer t;
+  t.start();
+
+  Mesh wrap;
+  CGAL::Alpha_wraps_3::internal::Alpha_wrap_3<Oracle> aw3(oracle);
+
+  // There is no limit on how many seeds can be used.
+  // However, the algorithm automatically determines whether a seed can be used
+  // to initialize the refinement based on a few conditions (distance to the offset, value of alpha, etc.)
+  // See internal function Alpha_wrap_3::initialize_from_cavities() for more information
+  std::vector<Point_3> seeds =
+  {
+    Point_3(0, 50, 0) // a point within the armadillo surface
+  };
+
+  aw3(alpha, offset, wrap, CGAL::parameters::seed_points(std::ref(seeds)));
+
+  t.stop();
+  std::cout << "Result: " << num_vertices(wrap) << " vertices, " << num_faces(wrap) << " faces" << std::endl;
+  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 + "_cavity_" + 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, wrap, CGAL::parameters::stream_precision(17));
+
+  return EXIT_SUCCESS;
+}

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

@@ -0,0 +1,298 @@
+// 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
+//                 Cedric Portaneri,
+//                 Mael Rouxel-Labbé
+//                 Andreas Fabri
+//                 Michael Hemmer
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_ALPHA_WRAP_AABB_TRAITS_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_ALPHA_WRAP_AABB_TRAITS_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Bbox_3.h>
+
+#include <array>
+#include <iostream>
+#include <bitset>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+template <typename K>
+struct Tetrahedron_with_outside_info
+{
+  using Kernel = K;
+  using Tetrahedron_3 = typename Kernel::Tetrahedron_3;
+  using Triangle_3 = typename Kernel::Triangle_3;
+
+  template <typename CellHandle>
+  Tetrahedron_with_outside_info(const CellHandle ch, const K& k)
+  {
+    typename K::Construct_bbox_3 bbox = k.construct_bbox_3_object();
+    typename K::Construct_tetrahedron_3 tetrahedron = k.construct_tetrahedron_3_object();
+    typename K::Construct_triangle_3 triangle = k.construct_triangle_3_object();
+
+    m_tet = tetrahedron(ch->vertex(0)->point(), ch->vertex(1)->point(),
+                        ch->vertex(2)->point(), ch->vertex(3)->point());
+    m_bbox = bbox(m_tet);
+
+    for(int i=0; i<4; ++i)
+    {
+      if(ch->neighbor(i)->info().is_outside)
+        m_b.set(i, true);
+
+      m_triangles[i] = triangle(ch->vertex((i+1)& 3)->point(),
+                                ch->vertex((i+2)& 3)->point(),
+                                ch->vertex((i+3)& 3)->point());
+      m_tbox[i] = bbox(m_triangles[i]);
+    }
+  }
+
+  Tetrahedron_with_outside_info& operator=(const Tetrahedron_with_outside_info& rhs) = default;
+
+  friend std::ostream& operator<<(std::ostream& os, const Tetrahedron_with_outside_info& t)
+  {
+    os << t.m_tet;
+    return os;
+  }
+
+  Tetrahedron_3 m_tet;
+  Bbox_3 m_bbox;
+  std::array<Bbox_3, 4> m_tbox;
+  std::array<Triangle_3, 4> m_triangles;
+  std::bitset<4> m_b;
+};
+
+template <typename K>
+class Ball_3
+  : private K::Sphere_3
+{
+  using FT = typename K::FT;
+  using Point_3 = typename K::Point_3;
+  using Sphere_3 = typename K::Sphere_3;
+
+public:
+  explicit Ball_3(const Sphere_3& s) : Sphere_3(s) { }
+
+public:
+  const Sphere_3& boundary() const { return static_cast<const Sphere_3&>(*this); }
+};
+
+template <typename GT>
+class Alpha_wrap_AABB_traits
+  : public GT
+{
+public:
+  using Ball_3 = internal::Ball_3<GT>;
+
+public:
+  Alpha_wrap_AABB_traits(const GT& gt = GT()) : GT(gt) { }
+
+public:
+  class Construct_ball_3
+  {
+    using FT = typename GT::FT;
+    using Point_3 = typename GT::Point_3;
+    using Ball_3 = internal::Ball_3<GT>;
+
+    const GT& m_base_traits;
+
+public:
+    Construct_ball_3(const GT& base_traits) : m_base_traits(base_traits) { }
+
+    Ball_3 operator()(const Point_3& p, const FT sqr)
+    {
+      return Ball_3(m_base_traits.construct_sphere_3_object()(p, sqr));
+    }
+  };
+
+  // Enrich the base's Do_intersect_3 with Tetrahedron_with_outside_info<K> and Ball_3<K> overloads
+  class Do_intersect_3
+    : public GT::Do_intersect_3
+  {
+    using Base = typename GT::Do_intersect_3;
+
+    using Point_3 = typename GT::Point_3;
+    using Segment_3 = typename GT::Segment_3;
+    using Triangle_3 = typename GT::Triangle_3;
+    using Ball_3 = internal::Ball_3<GT>;
+
+    const GT& m_base_traits;
+
+public:
+    Do_intersect_3(const GT& base_traits)
+      : Base(base_traits.do_intersect_3_object()),
+        m_base_traits(base_traits)
+    { }
+
+    using Base::operator();
+
+    // ======
+    // Ball_3
+
+    bool operator()(const Ball_3& ball,
+                    const Point_3& p) const
+    {
+      return !m_base_traits.has_on_unbounded_side_3_object()(ball.boundary(), p);
+    }
+
+    bool operator()(const Ball_3& ball,
+                    const Segment_3& s) const
+    {
+      typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
+      typename GT::Construct_vertex_3 vertex = m_base_traits.construct_vertex_3_object();
+      typename GT::Has_on_bounded_side_3 has_on_bounded_side = m_base_traits.has_on_bounded_side_3_object();
+
+      if(!do_overlap(bbox(ball.boundary()), bbox(s)))
+         return false;
+
+      if(Base::operator()(ball.boundary(), s))
+        return true;
+
+      return has_on_bounded_side(ball.boundary(), vertex(s, 0));
+    }
+
+    bool operator()(const Ball_3& ball,
+                    const Triangle_3& tr) const
+    {
+      typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
+      typename GT::Construct_vertex_3 vertex = m_base_traits.construct_vertex_3_object();
+      typename GT::Has_on_bounded_side_3 has_on_bounded_side = m_base_traits.has_on_bounded_side_3_object();
+
+      if(!do_overlap(bbox(ball.boundary()), bbox(tr)))
+         return false;
+
+      if(Base::operator()(ball.boundary(), tr))
+        return true;
+
+      return has_on_bounded_side(ball.boundary(), vertex(tr, 0));
+    }
+
+    bool operator()(const Ball_3& ball,
+                    const CGAL::Bbox_3& bb) const
+    {
+      typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
+      typename GT::Construct_point_3 point = m_base_traits.construct_point_3_object();
+      typename GT::Has_on_bounded_side_3 has_on_bounded_side = m_base_traits.has_on_bounded_side_3_object();
+
+      if(!do_overlap(bbox(ball.boundary()), bb))
+         return false;
+
+      if(Base::operator()(ball.boundary(), bb))
+        return true;
+
+      const Point_3 bbp = point(bb.xmin(), bb.ymin(), bb.zmin());
+
+      return has_on_bounded_side(ball.boundary(), bbp);
+    }
+
+    // ======
+    // Tetrahedron_with_outside_info
+
+    template <typename K>
+    bool operator()(const Tetrahedron_with_outside_info<K>& q,
+                    const Point_3& p) const
+    {
+      typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
+
+      const CGAL::Bbox_3 pbox = bbox(p);
+      if(!do_overlap(q.m_bbox, pbox))
+        return false;
+
+      for(int i=0; i<4; ++i)
+      {
+        if(!q.m_b.test(i) && do_overlap(q.m_tbox[i], pbox) && Base::operator()(q.m_triangles[i], p))
+          return true;
+      }
+
+      return m_base_traits.has_on_bounded_side_3_object()(q.m_tet, p);
+    }
+
+    template <typename K>
+    bool operator()(const Tetrahedron_with_outside_info<K>& q,
+                    const Segment_3& s) const
+    {
+      typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
+      typename GT::Construct_vertex_3 vertex = m_base_traits.construct_vertex_3_object();
+
+      const CGAL::Bbox_3 sbox = bbox(s);
+      if(!do_overlap(q.m_bbox, sbox))
+        return false;
+
+      for(int i=0; i<4; ++i)
+      {
+        if(!q.m_b.test(i) && do_overlap(q.m_tbox[i], sbox) && Base::operator()(q.m_triangles[i], s))
+          return true;
+      }
+
+      return m_base_traits.has_on_bounded_side_3_object()(q.m_tet, vertex(s, 0));
+    }
+
+    template <typename K>
+    bool operator()(const Tetrahedron_with_outside_info<K>& q,
+                    const Triangle_3& tr) const
+    {
+      typename GT::Construct_bbox_3 bbox = m_base_traits.construct_bbox_3_object();
+      typename GT::Construct_vertex_3 vertex = m_base_traits.construct_vertex_3_object();
+
+      const CGAL::Bbox_3 tbox = bbox(tr);
+      if(!do_overlap(q.m_bbox, tbox))
+        return false;
+
+      for(int i=0; i<4; ++i)
+      {
+        if(!q.m_b.test(i) && do_overlap(q.m_tbox[i], tbox) && Base::operator()(q.m_triangles[i], tr))
+          return true;
+      }
+
+      return m_base_traits.has_on_bounded_side_3_object()(q.m_tet, vertex(tr, 0));
+    }
+
+    template <typename K>
+    bool operator()(const Tetrahedron_with_outside_info<K>& q,
+                    const CGAL::Bbox_3& bb) const
+    {
+      if(!do_overlap(q.m_bbox, bb))
+        return false;
+
+      for(int i=0; i<4; ++i)
+      {
+        // this overload of do_intersect() must not filter based on q.m_b because
+        // it is called from the AABB_tree's traversal with a node's bounding box,
+        // and the fact that a facet is incident to an outside cell is irrelevant
+        // for the hierarchy of bounding boxes of the primitives.
+        if(do_overlap(q.m_tbox[i], bb) && Base::operator()(q.m_triangles[i], bb))
+          return true;
+      }
+
+      const Point_3 bbp = m_base_traits.construct_point_3_object()(bb.xmin(), bb.ymin(), bb.zmin());
+      return m_base_traits.has_on_bounded_side_3_object()(q.m_tet, bbp);
+    }
+  };
+
+  Construct_ball_3 construct_ball_3_object() const
+  {
+    return Construct_ball_3(static_cast<const GT&>(*this));
+  }
+
+  Do_intersect_3 do_intersect_3_object() const
+  {
+    return Do_intersect_3(static_cast<const GT&>(*this));
+  }
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_ALPHA_WRAP_AABB_TRAITS_H

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

@@ -0,0 +1,365 @@
+// 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)     : Mael Rouxel-Labbé
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_ORACLE_BASE_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_ORACLE_BASE_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Alpha_wrap_3/internal/offset_intersection.h>
+
+#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
+#include <CGAL/Default.h>
+
+#include <algorithm>
+#include <memory>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+template <typename AABBTraits>
+struct Default_traversal_traits
+{
+  using Projection_traits = CGAL::internal::AABB_tree::Projection_traits<AABBTraits>;
+
+  template <typename Query>
+  using Do_intersect_traits = CGAL::internal::AABB_tree::Do_intersect_traits<AABBTraits, Query>;
+
+  template <typename Query>
+  using First_intersection_traits = CGAL::internal::AABB_tree::First_intersection_traits<AABBTraits, Query>;
+};
+
+// Factorize the implementation of the functions calling the AABB tree
+template <typename AABBTree,
+          typename AABBTraversalTraits>
+struct AABB_tree_oracle_helper
+{
+  using Self = AABB_tree_oracle_helper<AABBTree, AABBTraversalTraits>;
+
+  using AABB_traits = typename AABBTree::AABB_traits;
+  using GT = typename AABB_traits::Geom_traits;
+
+  using FT = typename AABB_traits::FT;
+  using Point_3 = typename AABB_traits::Point_3;
+
+  template <typename Query>
+  static bool do_intersect(const Query& query,
+                           const AABBTree& tree)
+  {
+    CGAL_precondition(!tree.empty());
+
+    using Do_intersect_traits = typename AABBTraversalTraits::template Do_intersect_traits<Query>;
+
+    Do_intersect_traits traversal_traits(tree.traits());
+    tree.traversal(query, traversal_traits);
+    return traversal_traits.is_intersection_found();
+  }
+
+  static Point_3 closest_point(const Point_3& p,
+                               const AABBTree& tree)
+  {
+    CGAL_precondition(!tree.empty());
+
+    using Projection_traits = typename AABBTraversalTraits::Projection_traits;
+
+    const auto& hint = tree.best_hint(p);
+
+    Projection_traits projection_traits(hint.first, hint.second, tree.traits());
+    tree.traversal(p, projection_traits);
+    return projection_traits.closest_point();
+  }
+
+  static FT squared_distance(const Point_3& p,
+                             const AABBTree& tree)
+  {
+    CGAL_precondition(!tree.empty());
+
+    const Point_3 closest = Self::closest_point(p, tree);
+    return tree.traits().squared_distance_object()(p, closest);
+  }
+
+  static bool first_intersection(const Point_3& p, const Point_3& q, Point_3& o,
+                                 const FT offset_size,
+                                 const FT intersection_precision,
+                                 const AABBTree& tree)
+  {
+    CGAL_precondition(!tree.empty());
+
+    using AABB_distance_oracle = internal::AABB_distance_oracle<AABBTree, AABBTraversalTraits>;
+    using Offset_intersection = internal::Offset_intersection<GT, AABB_distance_oracle>;
+
+    AABB_distance_oracle dist_oracle(tree);
+    Offset_intersection offset_intersection(dist_oracle, offset_size, intersection_precision, 1 /*lip*/);
+    return offset_intersection.first_intersection(p, q, o);
+  }
+};
+
+template <typename GT,
+          typename AABBTree,
+          typename AABBTraversalTraits = CGAL::Default,
+          typename BaseOracle = int> // base oracle
+class AABB_tree_oracle
+  : public BaseOracle
+{
+protected:
+  using Geom_traits = GT;
+  using FT = typename GT::FT;
+  using Point_3 = typename GT::Point_3;
+
+  // When building oracle stacks, there are copies of (empty) trees, which isn't possible, thus pointer
+  using AABB_tree = AABBTree;
+  using AABB_tree_ptr = std::shared_ptr<AABB_tree>;
+  using AABB_traits = typename AABB_tree::AABB_traits;
+  using AABB_traversal_traits = typename Default::Get<AABBTraversalTraits,
+                                                      Default_traversal_traits<AABB_traits> >::type;
+
+  using AABB_helper = AABB_tree_oracle_helper<AABB_tree, AABB_traversal_traits>;
+
+public:
+  AABB_tree_oracle(const BaseOracle& base,
+                   const GT& gt)
+    : BaseOracle(base),
+      m_gt(gt),
+      m_tree_ptr(std::make_shared<AABB_tree>())
+  { }
+
+  AABB_tree_oracle(const AABB_tree_oracle&) = default;
+
+public:
+  const Geom_traits& geom_traits() const { return m_gt; }
+
+  AABB_tree& tree() { return *m_tree_ptr; }
+  const AABB_tree& tree() const { return *m_tree_ptr; }
+  BaseOracle& base() { return static_cast<BaseOracle&>(*this); }
+  const BaseOracle& base() const { return static_cast<const BaseOracle&>(*this); }
+
+  bool empty() const { return m_tree_ptr->empty(); }
+  bool do_call() const { return (!empty() || base().do_call()); }
+
+public:
+  typename AABB_tree::Bounding_box bbox() const
+  {
+    CGAL_precondition(do_call());
+
+    typename AABB_tree::Bounding_box bb;
+
+    if(base().do_call())
+      bb += base().bbox();
+
+    if(!empty())
+      bb += tree().bbox();
+
+    return bb;
+  }
+
+  template <typename T>
+  bool do_intersect(const T& t) const
+  {
+    CGAL_precondition(do_call());
+
+    if(base().do_call() && base().do_intersect(t))
+      return true;
+
+    if(!empty())
+      return AABB_helper::do_intersect(t, tree());
+
+    return false;
+  }
+
+  FT squared_distance(const Point_3& p) const
+  {
+    CGAL_precondition(do_call());
+
+    if(base().do_call())
+    {
+      if(!empty()) // both non empty
+      {
+        const FT base_sqd = base().squared_distance(p);
+        // @speed could do a smarter traversal, no need to search deeper than the current best
+        const FT this_sqd = AABB_helper::squared_distance(p, tree());
+        return (std::min)(base_sqd, this_sqd);
+      }
+      else // this level is empty
+      {
+        return base().squared_distance(p);
+      }
+    }
+    else // empty base
+    {
+      return AABB_helper::squared_distance(p, tree());
+    }
+  }
+
+  Point_3 closest_point(const Point_3& p) const
+  {
+    CGAL_precondition(do_call());
+
+    if(base().do_call())
+    {
+      if(!empty()) // both non empty
+      {
+        const Point_3 base_c = base().closest_point(p);
+        // @speed could do a smarter traversal, no need to search deeper than the current best
+        const Point_3 this_c = AABB_helper::closest_point(p, tree());
+        return (compare_distance_to_point(p, base_c, this_c) == CGAL::SMALLER) ? base_c : this_c;
+      }
+      else // this level is empty
+      {
+        return base().closest_point(p);
+      }
+    }
+    else // empty base
+    {
+      return AABB_helper::closest_point(p, tree());
+    }
+  }
+
+  bool first_intersection(const Point_3& p, const Point_3& q,
+                          Point_3& o,
+                          const FT offset_size,
+                          const FT intersection_precision) const
+  {
+    CGAL_precondition(do_call());
+
+    if(base().do_call())
+    {
+      if(!empty()) // both non empty
+      {
+        Point_3 base_o;
+        bool base_b = base().first_intersection(p, q, base_o, offset_size, intersection_precision);
+
+        if(base_b) // intersection found in base
+        {
+          // @speed could do a smarter traversal, no need to search deeper than the current best
+          Point_3 this_o;
+          bool this_b = AABB_helper::first_intersection(p, q, this_o, offset_size, intersection_precision, tree());
+          if(this_b)
+            o = (compare_distance_to_point(p, base_o, this_o) == SMALLER) ? base_o : this_o;
+          else
+            o = base_o;
+
+          return true;
+        }
+        else // no intersection found in non-empty base
+        {
+          return AABB_helper::first_intersection(p, q, o, offset_size, intersection_precision, tree());
+        }
+      }
+      else // this level is empty
+      {
+        return base().first_intersection(p, q, o, offset_size, intersection_precision);
+      }
+    }
+    else // empty base
+    {
+      return AABB_helper::first_intersection(p, q, o, offset_size, intersection_precision, tree());
+    }
+  }
+
+  bool first_intersection(const Point_3& p, const Point_3& q,
+                          Point_3& o,
+                          const FT offset_size) const
+  {
+    return first_intersection(p, q, o, offset_size, 1e-2 * offset_size);
+  }
+
+protected:
+  Geom_traits m_gt;
+  AABB_tree_ptr m_tree_ptr;
+};
+
+// partial specialization, when there is no further oracle beneath in the stack.
+//
+// `int` is used to denote the absence of base rather than `void`,
+// as to use the same constructor for both versions (thus requires a default construction)
+template <typename GT,
+          typename AABBTree,
+          typename AABBTraversalTraits>
+class AABB_tree_oracle<GT, AABBTree, AABBTraversalTraits, int>
+{
+protected:
+  using Geom_traits = GT;
+  using FT = typename GT::FT;
+  using Point_3 = typename GT::Point_3;
+
+  using AABB_tree = AABBTree;
+  using AABB_tree_ptr = std::shared_ptr<AABB_tree>;
+  using AABB_traits = typename AABB_tree::AABB_traits;
+  using AABB_traversal_traits = typename Default::Get<AABBTraversalTraits,
+                                                      Default_traversal_traits<AABB_traits> >::type;
+
+  using AABB_helper = AABB_tree_oracle_helper<AABB_tree, AABB_traversal_traits>;
+
+public:
+  AABB_tree_oracle(const int, // to have a common constructor API between both versions
+                   const GT& gt)
+    : m_gt(gt), m_tree_ptr(std::make_shared<AABB_tree>())
+  { }
+
+public:
+  const Geom_traits& geom_traits() const { return m_gt; }
+  AABB_tree& tree() { return *m_tree_ptr; }
+  const AABB_tree& tree() const { return *m_tree_ptr; }
+
+  bool empty() const { return m_tree_ptr->empty(); }
+  bool do_call() const { return !empty(); }
+
+public:
+  typename AABB_tree::Bounding_box bbox() const
+  {
+    CGAL_precondition(!empty());
+
+    return tree().bbox();
+  }
+
+  template <typename T>
+  bool do_intersect(const T& t) const
+  {
+    CGAL_precondition(!empty());
+    return AABB_helper::do_intersect(t, tree());
+  }
+
+  FT squared_distance(const Point_3& p) const
+  {
+    CGAL_precondition(!empty());
+    return AABB_helper::squared_distance(p, tree());
+  }
+
+  Point_3 closest_point(const Point_3& p) const
+  {
+    CGAL_precondition(!empty());
+    return AABB_helper::closest_point(p, tree());
+  }
+
+  bool first_intersection(const Point_3& p, const Point_3& q, Point_3& o,
+                          const FT offset_size, const FT intersection_precision) const
+  {
+    CGAL_precondition(!empty());
+    return AABB_helper::first_intersection(p, q, o, offset_size, intersection_precision, tree());
+  }
+
+  bool first_intersection(const Point_3& p, const Point_3& q, Point_3& o, const FT offset_size) const
+  {
+    CGAL_precondition(!empty());
+    return AABB_helper::first_intersection(p, q, o, offset_size, 1e-2 * offset_size, tree());
+  }
+
+private:
+  Geom_traits m_gt;
+  AABB_tree_ptr m_tree_ptr;
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_ORACLE_BASE_H

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

@@ -0,0 +1,126 @@
+// 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)     : Mael Rouxel-Labbé
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_POINT_SET_ORACLE_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_POINT_SET_ORACLE_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
+#include <CGAL/Alpha_wrap_3/internal/Oracle_base.h>
+
+#include <CGAL/AABB_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/property_map.h>
+
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <functional>
+#include <vector>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+// Just some typedefs for readability
+template <typename GT_>
+struct PS_oracle_traits
+{
+  using Geom_traits = Alpha_wrap_AABB_traits<GT_>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
+
+  using Points = std::vector<typename GT_::Point_3>;
+  using PR_iterator = typename Points::const_iterator;
+
+  using Primitive = AABB_primitive<PR_iterator,
+                                   Input_iterator_property_map<PR_iterator> /*DPM*/,
+                                   Input_iterator_property_map<PR_iterator> /*RPM*/,
+                                   CGAL::Tag_false, // not external
+                                   CGAL::Tag_false>; // no caching
+
+  using AABB_traits = CGAL::AABB_traits<Geom_traits, Primitive>;
+  using AABB_tree = CGAL::AABB_tree<AABB_traits>;
+};
+
+template <typename GT_,
+          typename BaseOracle = int>
+class Point_set_oracle
+  : public AABB_tree_oracle<typename PS_oracle_traits<GT_>::Geom_traits,
+                            typename PS_oracle_traits<GT_>::AABB_tree,
+                            CGAL::Default, // Default_traversal_traits<AABB_traits>
+                            BaseOracle>
+{
+  using PSOT = PS_oracle_traits<GT_>;
+  using Base_GT = GT_;
+
+public:
+  using Geom_traits = typename PSOT::Geom_traits;
+
+private:
+  using Points = typename PSOT::Points;
+  using AABB_tree = typename PSOT::AABB_tree;
+  using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, CGAL::Default, BaseOracle>;
+
+private:
+  Points m_points;
+
+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)
+  { }
+
+  Point_set_oracle(const Base_GT& gt,
+                   const BaseOracle& base_oracle = BaseOracle())
+    : Oracle_base(base_oracle, gt)
+  { }
+
+public:
+  // adds a range of points to the oracle
+  template <typename PointRange,
+            typename CGAL_NP_TEMPLATE_PARAMETERS>
+  void add_point_set(const PointRange& points,
+                     const CGAL_NP_CLASS& /*np*/ = CGAL::parameters::default_values())
+  {
+    if(points.empty())
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cout << "Warning: Input is empty " << 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));
+
+#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));
+
+    CGAL_postcondition(this->tree().size() == m_points.size());
+  }
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_POINT_SET_ORACLE_H

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

@@ -0,0 +1,124 @@
+// 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)     : Mael Rouxel-Labbé
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_SEGMENT_SOUP_ORACLE_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_SEGMENT_SOUP_ORACLE_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
+#include <CGAL/Alpha_wrap_3/internal/Oracle_base.h>
+
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_segment_primitive.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/property_map.h>
+
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <functional>
+#include <vector>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+// Just some typedefs for readability
+template <typename GT_>
+struct SS_oracle_traits
+{
+  using Geom_traits = Alpha_wrap_AABB_traits<GT_>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
+
+  using Segments = std::vector<typename GT_::Segment_3>;
+  using SR_iterator = typename Segments::const_iterator;
+
+  using Primitive = AABB_primitive<SR_iterator,
+                                   Input_iterator_property_map<SR_iterator>, // DPM
+                                   CGAL::internal::Source_of_segment_3_iterator_property_map<Geom_traits, SR_iterator>, // RPM
+                                   CGAL::Tag_false, // not external
+                                   CGAL::Tag_false>; // no caching
+
+  using AABB_traits = CGAL::AABB_traits<Geom_traits, Primitive>;
+  using AABB_tree = CGAL::AABB_tree<AABB_traits>;
+};
+
+template <typename GT_,
+          typename BaseOracle = int>
+class Segment_soup_oracle
+  : public AABB_tree_oracle<typename SS_oracle_traits<GT_>::Geom_traits,
+                            typename SS_oracle_traits<GT_>::AABB_tree,
+                            CGAL::Default, // Default_traversal_traits<AABB_traits>
+                            BaseOracle>
+{
+  using SSOT = SS_oracle_traits<GT_>;
+  using Base_GT = GT_;
+
+public:
+  using Geom_traits = typename SSOT::Geom_traits;
+
+private:
+  using Segments = typename SSOT::Segments;
+  using AABB_tree = typename SSOT::AABB_tree;
+  using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, CGAL::Default, BaseOracle>;
+
+private:
+  Segments m_segments;
+
+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)
+  { }
+
+  Segment_soup_oracle(const Base_GT& gt,
+                      const BaseOracle& base_oracle = BaseOracle())
+    : Oracle_base(base_oracle, gt)
+  { }
+
+public:
+  template <typename SegmentRange,
+            typename CGAL_NP_TEMPLATE_PARAMETERS>
+  void add_segment_soup(const SegmentRange& segments,
+                        const CGAL_NP_CLASS& /*np*/ = CGAL::parameters::default_values())
+  {
+    if(segments.empty())
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cout << "Warning: Input is empty " << 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));
+
+#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));
+
+    CGAL_postcondition(this->tree().size() == m_segments.size());
+  }
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_SEGMENT_SOUP_ORACLE_H

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

@@ -0,0 +1,177 @@
+// 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)     : Mael Rouxel-Labbé
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_MESH_ORACLE_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_MESH_ORACLE_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
+#include <CGAL/Alpha_wrap_3/internal/Oracle_base.h>
+#include <CGAL/Alpha_wrap_3/internal/splitting_helper.h>
+
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
+
+#include <algorithm>
+#include <iostream>
+#include <functional>
+#include <memory>
+#include <type_traits>
+#include <vector>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+// Just some typedefs for readability in the main oracle class
+template <typename GT_>
+struct TM_oracle_traits
+{
+  using Geom_traits = Alpha_wrap_AABB_traits<GT_>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
+
+  using Point_3 = typename Geom_traits::Point_3;
+  using AABB_traits = typename AABB_tree_splitter_traits<Point_3, Geom_traits>::AABB_traits;
+  using AABB_tree = typename AABB_tree_splitter_traits<Point_3, Geom_traits>::AABB_tree;
+};
+
+// @speed could do a partial specialization 'subdivide = false' with simpler code for speed?
+template <typename GT_,
+          typename BaseOracle = int,
+          bool subdivide = true>
+class Triangle_mesh_oracle
+  : // this is the base that handles calls to the AABB tree
+    public AABB_tree_oracle<typename TM_oracle_traits<GT_>::Geom_traits,
+                            typename TM_oracle_traits<GT_>::AABB_tree,
+                            typename std::conditional<
+                                      /*condition*/subdivide,
+                                      /*true*/Splitter_traversal_traits<typename TM_oracle_traits<GT_>::AABB_traits>,
+                                      /*false*/Default_traversal_traits<typename TM_oracle_traits<GT_>::AABB_traits> >::type,
+                            BaseOracle>,
+    // this is the base that handles splitting input faces and inserting them into the AABB tree
+    public AABB_tree_oracle_splitter<subdivide,
+                                     typename TM_oracle_traits<GT_>::Point_3,
+                                     typename TM_oracle_traits<GT_>::Geom_traits>
+{
+  using TMOT = TM_oracle_traits<GT_>;
+  using Base_GT = GT_;
+
+public:
+  using Geom_traits = typename TMOT::Geom_traits;
+
+private:
+  using Point_3 = typename Geom_traits::Point_3;
+  using Triangle_3 = typename Geom_traits::Triangle_3;
+
+  using AABB_traits = typename TMOT::AABB_traits;
+  using AABB_tree = typename TMOT::AABB_tree;
+  using AABB_traversal_traits = typename std::conditional<
+                                  /*condition*/subdivide,
+                                  /*true*/Splitter_traversal_traits<AABB_traits>,
+                                  /*false*/Default_traversal_traits<AABB_traits> >::type;
+
+  using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, AABB_traversal_traits, BaseOracle>;
+  using Splitter_base = AABB_tree_oracle_splitter<subdivide, Point_3, Geom_traits>;
+
+public:
+  // Constructors
+  //
+  // When using this constructor (and thus doing actual splitting), note that the oracle
+  // will be adapted to this particular 'alpha', and so when calling again AW3(other_alpha)
+  // the oracle might not have performed a split that is adapted to this other alpha value.
+  Triangle_mesh_oracle(const double alpha,
+                       const BaseOracle& base_oracle = BaseOracle(),
+                       const Base_GT& gt = Base_GT())
+    : Oracle_base(base_oracle, gt), Splitter_base(alpha)
+  {
+    Splitter_base::initialize_tree_property_maps(this->tree());
+  }
+
+  Triangle_mesh_oracle(const double alpha,
+                       const Base_GT& gt,
+                       const BaseOracle& base_oracle = BaseOracle())
+    : Triangle_mesh_oracle(alpha, base_oracle, gt)
+  { }
+
+ Triangle_mesh_oracle(const BaseOracle& base_oracle,
+                      const Base_GT& gt = Base_GT())
+   : Triangle_mesh_oracle(0. /*alpha*/, base_oracle, gt)
+ { }
+
+ Triangle_mesh_oracle(const Base_GT& gt,
+                      const BaseOracle& base_oracle = BaseOracle())
+   : Triangle_mesh_oracle(0. /*alpha*/, base_oracle, gt)
+ { }
+
+ Triangle_mesh_oracle()
+   : Triangle_mesh_oracle(0. /*alpha*/, BaseOracle(), Base_GT())
+ { }
+
+public:
+  template <typename TriangleMesh,
+            typename CGAL_NP_TEMPLATE_PARAMETERS>
+  void add_triangle_mesh(const TriangleMesh& tmesh,
+                         const CGAL_NP_CLASS& np = CGAL::parameters::default_values())
+  {
+    using parameters::get_parameter;
+    using parameters::choose_parameter;
+
+    using face_descriptor = typename boost::graph_traits<TriangleMesh>::face_descriptor;
+
+    using VPM = typename GetVertexPointMap<TriangleMesh>::const_type;
+    using Point_ref = typename boost::property_traits<VPM>::reference;
+
+    CGAL_precondition(CGAL::is_triangle_mesh(tmesh));
+
+    if(is_empty(tmesh))
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cout << "Warning: Input is empty " << std::endl;
+#endif
+      return;
+    }
+
+#ifdef CGAL_AW3_DEBUG
+    std::cout << "Insert into AABB tree (faces)..." << std::endl;
+#endif
+
+    VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                               get_const_property_map(vertex_point, tmesh));
+    CGAL_static_assertion((std::is_same<typename boost::property_traits<VPM>::value_type, Point_3>::value));
+
+    Splitter_base::reserve(num_faces(tmesh));
+
+    for(face_descriptor f : faces(tmesh))
+    {
+      if(Polygon_mesh_processing::is_degenerate_triangle_face(f, tmesh, np))
+        continue;
+
+      const Point_ref p0 = get(vpm, source(halfedge(f, tmesh), tmesh));
+      const Point_ref p1 = get(vpm, target(halfedge(f, tmesh), tmesh));
+      const Point_ref p2 = get(vpm, target(next(halfedge(f, tmesh), tmesh), tmesh));
+
+      const Triangle_3 tr = this->geom_traits().construct_triangle_3_object()(p0, p1, p2);
+
+      Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits());
+    }
+
+#ifdef CGAL_AW3_DEBUG
+    std::cout << "Tree: " << this->tree().size() << " primitives (" << num_faces(tmesh) << " faces in input)" << std::endl;
+#endif
+  }
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_MESH_ORACLE_H

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

@@ -0,0 +1,198 @@
+// 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)     : Mael Rouxel-Labbé
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_SOUP_ORACLE_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_SOUP_ORACLE_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
+#include <CGAL/Alpha_wrap_3/internal/Oracle_base.h>
+#include <CGAL/Alpha_wrap_3/internal/splitting_helper.h>
+
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Named_function_parameters.h>
+
+#include <boost/range/value_type.hpp>
+
+#include <algorithm>
+#include <iostream>
+#include <functional>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+// Just some typedefs for readability
+template <typename GT_>
+struct TS_oracle_traits
+{
+  using Geom_traits = Alpha_wrap_AABB_traits<GT_>; // Wrap the kernel to add Ball_3 + custom Do_intersect_3
+  using Point_3 = typename Geom_traits::Point_3;
+  using AABB_traits = typename AABB_tree_splitter_traits<Point_3, Geom_traits>::AABB_traits;
+  using AABB_tree = typename AABB_tree_splitter_traits<Point_3, Geom_traits>::AABB_tree;
+};
+
+template <typename GT_,
+          typename BaseOracle = int,
+          bool subdivide = true>
+class Triangle_soup_oracle
+  : // this is the base that handles calls to the AABB tree
+    public AABB_tree_oracle<typename TS_oracle_traits<GT_>::Geom_traits,
+                            typename TS_oracle_traits<GT_>::AABB_tree,
+                            typename std::conditional<
+                              /*condition*/subdivide,
+                              /*true*/Splitter_traversal_traits<typename TS_oracle_traits<GT_>::AABB_traits>,
+                              /*false*/Default_traversal_traits<typename TS_oracle_traits<GT_>::AABB_traits> >::type,
+                            BaseOracle>,
+    // this is the base that handles splitting input faces and inserting them into the AABB tree
+    public AABB_tree_oracle_splitter<subdivide,
+                                     typename TS_oracle_traits<GT_>::Point_3,
+                                     typename TS_oracle_traits<GT_>::Geom_traits>
+{
+  using TSOT = TS_oracle_traits<GT_>;
+  using Base_GT = GT_;
+
+public:
+  using Geom_traits = typename TSOT::Geom_traits;
+
+private:
+  using Point_3 = typename Geom_traits::Point_3;
+  using Triangle_3 = typename Geom_traits::Triangle_3;
+
+  using AABB_traits = typename TSOT::AABB_traits;
+  using AABB_tree = typename TSOT::AABB_tree;
+  using AABB_traversal_traits = typename std::conditional<
+                                  /*condition*/subdivide,
+                                  /*true*/Splitter_traversal_traits<AABB_traits>,
+                                  /*false*/Default_traversal_traits<AABB_traits> >::type;
+
+  using Oracle_base = AABB_tree_oracle<Geom_traits, AABB_tree, AABB_traversal_traits, BaseOracle>;
+  using Splitter_base = AABB_tree_oracle_splitter<subdivide, Point_3, Geom_traits>;
+
+public:
+  // Constructors
+  //
+  // When using this constructor (and thus doing actual splitting), note that the oracle
+  // will be adapted to this particular 'alpha', and so when calling again AW3(other_alpha)
+  // the oracle might not have performed a split that is adapted to this other alpha value.
+  Triangle_soup_oracle(const double alpha,
+                       const BaseOracle& base_oracle = BaseOracle(),
+                       const Base_GT& gt = Base_GT())
+    : Oracle_base(base_oracle, gt), Splitter_base(alpha)
+  {
+    Splitter_base::initialize_tree_property_maps(this->tree());
+  }
+
+  Triangle_soup_oracle(const double alpha,
+                       const Base_GT& gt,
+                       const BaseOracle& base_oracle = BaseOracle())
+    : Triangle_soup_oracle(alpha, base_oracle, gt)
+  { }
+
+ Triangle_soup_oracle(const BaseOracle& base_oracle,
+                      const Base_GT& gt = Base_GT())
+   : Triangle_soup_oracle(0. /*alpha*/, base_oracle, gt)
+ { }
+
+ Triangle_soup_oracle(const Base_GT& gt,
+                      const BaseOracle& base_oracle = BaseOracle())
+   : Triangle_soup_oracle(0. /*alpha*/, base_oracle, gt)
+ { }
+
+ Triangle_soup_oracle()
+   : Triangle_soup_oracle(0. /*alpha*/, BaseOracle(), Base_GT())
+ { }
+
+public:
+  template <typename PointRange, typename FaceRange,
+            typename CGAL_NP_TEMPLATE_PARAMETERS>
+  void add_triangle_soup(const PointRange& points,
+                         const FaceRange& faces,
+                         const CGAL_NP_CLASS& np = CGAL::parameters::default_values())
+  {
+    using parameters::choose_parameter;
+    using parameters::get_parameter;
+
+    using PPM = typename GetPointMap<PointRange, CGAL_NP_CLASS>::const_type;
+    using Point_ref = typename boost::property_traits<PPM>::reference;
+
+    using Face = typename boost::range_value<FaceRange>::type;
+
+    if(points.empty() || faces.empty())
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cout << "Warning: Input is empty " << std::endl;
+#endif
+      return;
+    }
+
+#ifdef CGAL_AW3_DEBUG
+    std::cout << "Insert into AABB Tree (triangles)..." << std::endl;
+#endif
+
+    PPM pm = choose_parameter<PPM>(get_parameter(np, internal_np::point_map));
+    CGAL_static_assertion((std::is_same<typename boost::property_traits<PPM>::value_type, Point_3>::value));
+
+    Splitter_base::reserve(faces.size());
+
+    typename Geom_traits::Construct_triangle_3 triangle = this->geom_traits().construct_triangle_3_object();
+    typename Geom_traits::Is_degenerate_3 is_degenerate = this->geom_traits().is_degenerate_3_object();
+
+    for(const Face& f : faces)
+    {
+      CGAL_precondition(std::distance(std::cbegin(f), std::cend(f)) == 3);
+
+      auto vi = std::cbegin(f);
+      CGAL_assertion(*vi < points.size());
+      Point_ref p0 = get(pm, points[*vi++]);
+      CGAL_assertion(*vi < points.size());
+      Point_ref p1 = get(pm, points[*vi++]);
+      CGAL_assertion(*vi < points.size());
+      Point_ref p2 = get(pm, points[*vi]);
+
+      const Triangle_3 tr = triangle(p0, p1, p2);
+      if(is_degenerate(tr))
+        continue;
+
+      Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits());
+    }
+
+#ifdef CGAL_AW3_DEBUG
+    std::cout << "Tree: " << this->tree().size() << " primitives (" << faces.size() << " faces in input)" << std::endl;
+#endif
+  }
+
+  template <typename TriangleRange,
+            typename CGAL_NP_TEMPLATE_PARAMETERS>
+  void add_triangle_soup(const TriangleRange& triangles,
+                         const CGAL_NP_CLASS& /*np*/ = CGAL::parameters::default_values())
+  {
+    typename Geom_traits::Is_degenerate_3 is_degenerate = this->geom_traits().is_degenerate_3_object();
+
+    for(const Triangle_3& tr : triangles)
+    {
+      if(is_degenerate(tr))
+        continue;
+
+      Splitter_base::split_and_insert_datum(tr, this->tree(), this->geom_traits());
+    }
+  }
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_TRIANGLE_SOUP_ORACLE_H

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

@@ -0,0 +1,101 @@
+// 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
+//                 Cedric Portaneri,
+//                 Mael Rouxel-Labbé
+//                 Andreas Fabri
+//                 Michael Hemmer
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_GATE_PRIORITY_QUEUE_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_GATE_PRIORITY_QUEUE_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <boost/property_map/property_map.hpp>
+
+#include <cassert>
+#include <iostream>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+// Represents an alpha-traversable facet in the mutable priority queue
+template <typename DT3>
+class Gate
+{
+  using Facet = typename DT3::Facet;
+  using FT = typename DT3::Geom_traits::FT;
+
+private:
+  Facet m_facet;
+  FT m_priority; // circumsphere sq_radius
+  bool m_is_artificial_facet;
+
+public:
+  // Constructors
+  Gate(const Facet& facet,
+       const FT& priority,
+       const bool is_artificial_facet)
+    :
+      m_facet(facet),
+      m_priority(priority),
+      m_is_artificial_facet(is_artificial_facet)
+  {
+    CGAL_assertion(priority >= 0);
+  }
+
+  // This overload is only used for contains() and erase(), priority and bbox flag are dummy value
+  Gate(const Facet& facet)
+    : Gate(facet, 0, false)
+  { }
+
+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; }
+};
+
+struct Less_gate
+{
+  template <typename DT3>
+  bool operator()(const Gate<DT3>& a, const Gate<DT3>& 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;
+    return a.priority() > b.priority();
+  }
+};
+
+template <typename DT3>
+struct Gate_ID_PM
+{
+  using key_type = Gate<DT3>;
+  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;
+
+    const Facet& f = k.facet();
+    return (4 * f.first->time_stamp() + f.second);
+  }
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_GATE_PRIORITY_QUEUE_H

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

@@ -0,0 +1,161 @@
+// 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
+//                 Cedric Portaneri,
+//                 Mael Rouxel-Labbé
+//                 Andreas Fabri
+//                 Michael Hemmer
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_GEOMETRY_UTILS_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_GEOMETRY_UTILS_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Simple_cartesian.h>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+template <typename K>
+struct Orientation_of_circumcenter
+{
+  typedef typename K::Point_3 Point_3;
+
+  typedef Orientation result_type;
+
+  Orientation operator()(const Point_3& p, const Point_3& q, const Point_3& r,
+                         const Point_3& ccp, const Point_3& ccq, const Point_3& ccr, const Point_3& ccs) const
+  {
+    Point_3 cc = circumcenter(ccp, ccq, ccr, ccs);
+    return orientation(p, q, r, cc);
+  }
+};
+
+template <typename Dt>
+bool
+less_squared_radius_of_min_empty_sphere(typename Dt::Geom_traits::FT sq_alpha,
+                                        const typename Dt::Facet& fh,
+                                        const Dt& dt)
+{
+  using Cell_handle = typename Dt::Cell_handle;
+  using Point = typename Dt::Point;
+
+  using CK = typename Dt::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;
+
+  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));
+
+  // 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))
+  {
+    Orientation ori = orientation_of_circumcenter(p1, p2, p3,
+                                                  dt.point(c, 0), dt.point(c, 1),
+                                                  dt.point(c, 2), dt.point(c, 3));
+
+    if(ori == POSITIVE)
+    {
+      Comparison_result cr = compare_squared_radius(p1, p2, p3, sq_alpha);
+      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),
+                                                    sq_alpha);
+      return cr == LARGER;
+    }
+  }
+
+  if(dt.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));
+
+    if(ori == NEGATIVE)
+    {
+      Comparison_result cr = compare_squared_radius(p1, p2, p3, sq_alpha);
+      return cr == LARGER;
+    }
+    else
+    {
+      Comparison_result cr = compare_squared_radius(dt.point(n, 0), dt.point(n, 1),
+                                                    dt.point(n, 2), dt.point(n, 3),
+                                                    sq_alpha);
+      return cr == LARGER;
+    }
+  }
+
+  // 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)) !=
+     orientation_of_circumcenter(p1, p2, p3,
+                                 dt.point(n, 0), dt.point(n, 1), dt.point(n, 2), dt.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)
+              << " 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),
+                                                  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))
+              << " 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),
+                                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))
+              << " sq alpha " << sq_alpha << std::endl;
+#endif
+
+    return cr == LARGER;
+  }
+}
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_GEOMETRY_UTILS_H

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

@@ -0,0 +1,220 @@
+// 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
+//                 Cedric Portaneri,
+//                 Mael Rouxel-Labbé
+//                 Andreas Fabri
+//                 Michael Hemmer
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_OFFSET_INTERSECTION_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_OFFSET_INTERSECTION_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/number_utils.h>
+
+#include <boost/algorithm/clamp.hpp>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+template <typename AABBTree,
+          typename AABBTraversalTraits>
+struct AABB_tree_oracle_helper;
+
+template <typename AABBTree, typename AABBTraversalTraits>
+struct AABB_distance_oracle
+{
+  using FT = typename AABBTree::FT;
+  using Point_3 = typename AABBTree::Point;
+
+  using AABB_helper = AABB_tree_oracle_helper<AABBTree, AABBTraversalTraits>;
+
+  AABB_distance_oracle(const AABBTree& tree) : tree(tree) { }
+
+  FT operator()(const Point_3& p) const
+  {
+    return approximate_sqrt(AABB_helper::squared_distance(p, tree));
+  }
+
+private:
+  const AABBTree& tree;
+};
+
+// @todo even with EPECK, the precision cannot be 0 (otherwise it will not converge),
+// thus exactness is pointless. Might as well use a cheap kernel (e.g. SC<double>), as long
+// as there exists a mechanism to catch when the cheap kernel fails to converge (iterations?
+// see also Tr_3::locate() or Mesh_3::Robust_intersection_traits_3.h)
+template <class Kernel, class DistanceOracle>
+class Offset_intersection
+{
+  using FT = typename Kernel::FT;
+  using Point_2 = typename Kernel::Point_2;
+  using Point_3 = typename Kernel::Point_3;
+  using Vector_3 = typename Kernel::Vector_3;
+
+public:
+  Offset_intersection(const DistanceOracle& oracle,
+                      const FT& off,
+                      const FT& prec,
+                      const FT& lip)
+    : dist_oracle(oracle), offset(off), precision(prec), lipschitz(lip)
+  { }
+
+  bool first_intersection(const Point_3& s,
+                          const Point_3& t,
+                          Point_3& output_pt)
+  {
+    source = s;
+    target = t;
+    seg_length = approximate_sqrt(squared_distance(s, t));
+    seg_unit_v = (t - s) / seg_length;
+    const Point_2 p0 { 0, dist_oracle(source) };
+    const Point_2 p1 { seg_length, dist_oracle(target) };
+
+    return recursive_dichotomic_search(p0, p1, output_pt);
+  }
+
+private:
+  Point_3 source;
+  Point_3 target;
+  FT seg_length;
+  Vector_3 seg_unit_v;
+  DistanceOracle dist_oracle;
+  FT offset;
+  FT precision;
+  FT lipschitz;
+
+  template <class Point>
+  bool recursive_dichotomic_search(const Point_2& s, const Point_2& t,
+                                   Point& output_pt)
+  {
+    if(CGAL::abs(s.x() - t.x()) < precision)
+    {
+      if(CGAL::abs(s.y() - offset) < precision)
+      {
+        const FT x_clamp = boost::algorithm::clamp<FT>(s.x(), FT{0}, seg_length);
+        output_pt = source + (seg_unit_v * x_clamp);
+        return true;
+      }
+
+      return false;
+    }
+
+    const bool sign_s = (s.y() > offset);
+    const bool sign_t = (t.y() > offset);
+    const FT gs_a = (sign_s) ? -lipschitz : lipschitz;
+    const FT gs_b = s.y() - (gs_a * s.x());
+    const FT gt_a = (sign_t) ? lipschitz : -lipschitz;
+
+    const FT gt_b = t.y() - (gt_a * t.x());
+    FT ms = (offset - gs_b) / gs_a;
+    FT mt = (offset - gt_b) / gt_a;
+
+    // early exit if there is no intersection
+    if(sign_s == sign_t)
+    {
+      FT ui = (gt_b - gs_b) / (gs_a - gt_a);
+      const FT gs_ui = (gs_a * ui) + gs_b;
+      if((sign_s && (gs_ui > offset)) || (!sign_s && (gs_ui < offset)))
+      {
+        if(CGAL::abs(s.y() - offset) < precision)
+        {
+          const FT x_clamp = boost::algorithm::clamp<FT>(s.x(), FT{0}, seg_length);
+          output_pt = source + (seg_unit_v * x_clamp);
+          return true;
+        }
+        else if(CGAL::abs(t.y() - offset) < precision)
+        {
+          const FT x_clamp = boost::algorithm::clamp<FT>(t.x(), FT{0}, seg_length);
+          output_pt = source + (seg_unit_v * x_clamp);
+          return true;
+        }
+
+        return false;
+      }
+      else
+      {
+        ms = boost::algorithm::clamp<FT>(ms, FT{0}, seg_length);
+        ui = boost::algorithm::clamp<FT>(ui, FT{0}, seg_length);
+        mt = boost::algorithm::clamp<FT>(mt, FT{0}, seg_length);
+        const Point_2 ms_pt { ms, dist_oracle(source + (seg_unit_v * ms)) };
+        const Point_2 ui_pt { ui, dist_oracle(source + (seg_unit_v * ui)) };
+        const Point_2 mt_pt { mt, dist_oracle(source + (seg_unit_v * mt)) };
+
+        if(CGAL::abs(ms_pt.y() - offset) < precision)
+        {
+          const FT x_clamp = boost::algorithm::clamp<FT>(ms_pt.x(), FT{0}, seg_length);
+          output_pt = source + (seg_unit_v * x_clamp);
+          return true;
+        }
+        else if(CGAL::abs(ui_pt.y() - offset) < precision)
+        {
+          const FT x_clamp = boost::algorithm::clamp<FT>(ui_pt.x(), FT{0}, seg_length);
+          output_pt = source + (seg_unit_v * x_clamp);
+          return true;
+        }
+        else if(CGAL::abs(mt_pt.y() - offset) < precision)
+        {
+          const FT x_clamp = boost::algorithm::clamp<FT>(mt_pt.x(), FT{0}, seg_length);
+          output_pt = source + (seg_unit_v * x_clamp);
+          return true;
+        }
+
+        return (recursive_dichotomic_search(ms_pt, ui_pt, output_pt) ||
+                recursive_dichotomic_search(ui_pt, mt_pt, output_pt));
+      }
+    }
+    else // there is an intersection
+    {
+      if(CGAL::abs(mt - ms) <= precision) // linear approximation
+      {
+        const FT fsft_a = (t.y() - s.y()) / (t.x() - s.x());
+        const FT fsft_b = s.y() - fsft_a * s.x();
+        FT m_fsft;
+        if(fsft_a == FT{0})
+        {
+          if(CGAL::abs(s.y() - offset) < precision)
+            m_fsft = s.x();
+          else
+            return false;
+        }
+        else
+        {
+          m_fsft = (offset - fsft_b) / fsft_a;
+        }
+        m_fsft = boost::algorithm::clamp<FT>(m_fsft, FT{0}, seg_length);
+        output_pt = source + (seg_unit_v * m_fsft);
+        return true;
+      }
+      else
+      {
+        FT m = (ms + mt) / FT{2};
+        ms = boost::algorithm::clamp<FT>(ms, FT{0}, seg_length);
+        m = boost::algorithm::clamp<FT>(m, FT{0}, seg_length);
+        mt = boost::algorithm::clamp<FT>(mt, FT{0}, seg_length);
+
+        const Point_2 ms_pt { ms, dist_oracle(source + (seg_unit_v * ms)) };
+        const Point_2 m_pt { m, dist_oracle(source + (seg_unit_v * m)) };
+        const Point_2 mt_pt { mt, dist_oracle(source + (seg_unit_v * mt)) };
+
+        return (recursive_dichotomic_search(ms_pt, m_pt, output_pt) ||
+                recursive_dichotomic_search(m_pt, mt_pt, output_pt));
+      }
+    }
+  }
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_OFFSET_INTERSECTION_H

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

@@ -0,0 +1,25 @@
+// 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)     : Mael Rouxel-Labbé
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_ORACLES_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_ORACLES_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
+#include <CGAL/Alpha_wrap_3/internal/offset_intersection.h>
+
+#include <CGAL/Alpha_wrap_3/internal/Triangle_mesh_oracle.h>
+#include <CGAL/Alpha_wrap_3/internal/Triangle_soup_oracle.h>
+#include <CGAL/Alpha_wrap_3/internal/Segment_soup_oracle.h>
+#include <CGAL/Alpha_wrap_3/internal/Point_set_oracle.h>
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_ORACLES_H

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

@@ -0,0 +1,424 @@
+// 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)     : Mael Rouxel-Labbé
+//
+#ifndef CGAL_ALPHA_WRAP_3_INTERNAL_SPLITTING_HELPER_H
+#define CGAL_ALPHA_WRAP_3_INTERNAL_SPLITTING_HELPER_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_AABB_traits.h>
+
+#include <CGAL/AABB_tree/internal/AABB_traversal_traits.h>
+#include <CGAL/AABB_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/array.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Container_helper.h>
+#include <CGAL/property_map.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Cartesian_converter.h>
+
+#include <queue>
+#include <unordered_set>
+#include <utility>
+#include <vector>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+// std::vector to property map
+// Not using Pointer_property_map because the underlying range is empty initially and can change
+template <typename T>
+struct Vector_property_map
+{
+  using Range = std::vector<T>;
+
+  using key_type = std::size_t;
+  using value_type = T;
+  using reference = value_type&;
+  using category = boost::read_write_property_map_tag;
+
+  Vector_property_map() : m_range_ptr(std::make_shared<Range>()) { }
+
+  inline friend void put(const Vector_property_map& map, const key_type& k, const value_type& v)
+  {
+    CGAL_precondition(map.m_range_ptr != nullptr);
+
+    if(k >= map.m_range_ptr->size())
+      map.m_range_ptr->resize(k+1);
+
+    map.m_range_ptr->operator[](k) = v;
+  }
+
+  inline friend reference get(const Vector_property_map& map, const key_type& k)
+  {
+    CGAL_precondition(map.m_range_ptr != nullptr);
+    return map.m_range_ptr->operator[](k);
+  }
+
+  Range& range() { return *m_range_ptr; }
+  const Range& range() const { return *m_range_ptr; }
+
+private:
+  std::shared_ptr<Range> m_range_ptr;
+};
+
+// Same as the standard traversal traits, but for multiple primitives per datum: the final operation
+// done on the datum is only performed once.
+template <typename AABBTraits>
+struct Splitter_traversal_traits
+{
+  // -----------------------------------------------------------------------------------------------
+  class Projection_traits
+    : public CGAL::internal::AABB_tree::Projection_traits<AABBTraits>
+  {
+    using Base = CGAL::internal::AABB_tree::Projection_traits<AABBTraits>;
+    using Point = typename AABBTraits::Point_3;
+    using Primitive = typename AABBTraits::Primitive;
+
+    std::unordered_set<std::size_t> visited_data;
+
+public:
+    Projection_traits(const Point& hint,
+                      const typename Primitive::Id& hint_primitive,
+                      const AABBTraits& traits)
+      : Base(hint, hint_primitive, traits)
+    { }
+
+    void intersection(const Point& query, const Primitive& primitive)
+    {
+      // check a datum only once
+      auto is_insert_successful = visited_data.insert(primitive.id().second/*unique input face ID*/);
+      if(!is_insert_successful.second)
+        return;
+
+      return Base::intersection(query, primitive);
+    }
+  };
+
+  // -----------------------------------------------------------------------------------------------
+  template <typename Query>
+  class Do_intersect_traits
+    : public CGAL::internal::AABB_tree::Do_intersect_traits<AABBTraits, Query>
+  {
+    using Base = CGAL::internal::AABB_tree::Do_intersect_traits<AABBTraits, Query>;
+    using Primitive = typename AABBTraits::Primitive;
+
+    std::unordered_set<std::size_t> visited_data;
+
+public:
+    Do_intersect_traits(const AABBTraits& traits) : Base(traits) { }
+
+    void intersection(const Query& query, const Primitive& primitive)
+    {
+      // check a datum only once
+      auto is_insert_successful = visited_data.insert(primitive.id().second/*unique input face ID*/);
+      if(!is_insert_successful.second)
+        return;
+
+      return Base::intersection(query, primitive);
+    }
+  };
+
+  // -----------------------------------------------------------------------------------------------
+  template <typename Query>
+  class First_intersection_traits
+    : public CGAL::internal::AABB_tree::First_intersection_traits<AABBTraits, Query>
+  {
+    using Base = CGAL::internal::AABB_tree::First_intersection_traits<AABBTraits, Query>;
+    using Primitive = typename AABBTraits::Primitive;
+
+    std::unordered_set<std::size_t> visited_data;
+
+public:
+    First_intersection_traits(const AABBTraits& traits) : Base(traits) { }
+
+    void intersection(const Query& query, const Primitive& primitive)
+    {
+      // check a datum only once
+      auto is_insert_successful = visited_data.insert(primitive.id().second/*unique input face ID*/);
+      if(!is_insert_successful.second)
+        return;
+
+      return Base::intersection(query, primitive);
+    }
+  };
+};
+
+// Dissociated from the class `AABB_tree_oracle_splitter` for clarity (the AABB_tree is a template
+// parameter of the other base of the oracle too)
+template <typename Point, typename GT>
+struct AABB_tree_splitter_traits
+{
+  using Triangle_3 = typename GT::Triangle_3;
+
+  // Below is a lot of trouble to cover a single datum with multiple primitives using smaller bboxes
+
+  // The input face ID serves when traversing the tree, to avoid doing the same intersection()
+  // on the same datum seen from different primitives.
+  //
+  // Technically, FPM could type-erase the mesh and the VPM, as it currently forces all independant
+  // inputs to have the same types. This is not such much of an issue for the mesh type,
+  // but it can be annoying for the VPM type.
+  using ID = std::pair<std::size_t /*primitive ID*/, std::size_t /*input face ID*/>;
+  using IDPM = CGAL::First_of_pair_property_map<ID>;
+
+  // Primitive ID --> box vector pos --> Bounding Box
+  using BPMB = internal::Vector_property_map<CGAL::Bbox_3>;
+  using BPM = CGAL::Property_map_binder<IDPM, BPMB>;
+
+  // Primitive ID --> point vector pos --> Reference Point
+  using RPPMB = internal::Vector_property_map<Point>;
+  using RPPM = CGAL::Property_map_binder<IDPM, RPPMB>;
+
+  // Primitive ID --> Datum pos vector pos --> Datum pos --> Datum
+  // The vector of data has size nf, but the vector of datum pos has size tree.size()
+  using DPPMB = internal::Vector_property_map<std::size_t>; // pos --> Datum pos
+  using DPPM = CGAL::Property_map_binder<IDPM, DPPMB>; // PID --> Datum pos
+  using DPMB = internal::Vector_property_map<Triangle_3>; // Datum pos --> Datum
+  using DPM = CGAL::Property_map_binder<DPPM, DPMB>; // PID --> Datum
+
+  using Primitive = CGAL::AABB_primitive<ID, DPM, RPPM,
+                                         CGAL::Tag_true /*external pmaps*/,
+                                         CGAL::Tag_false /*no caching*/>;
+
+  using AABB_traits = CGAL::AABB_traits<GT, Primitive, BPM>;
+  using AABB_tree = CGAL::AABB_tree<AABB_traits>;
+};
+
+template <bool subdivide, typename Point, typename GT>
+struct AABB_tree_oracle_splitter
+{
+  using FT = typename GT::FT;
+  using Triangle_3 = typename GT::Triangle_3;
+
+  using ATST = AABB_tree_splitter_traits<Point, GT>;
+
+  using BPM = typename ATST::BPM;
+  using RPPM = typename ATST::RPPM;
+  using DPPMB = typename ATST::DPPMB;
+  using DPPM = typename ATST::DPPM;
+  using DPMB = typename ATST::DPMB;
+  using DPM = typename ATST::DPM;
+
+  using ID = typename ATST::ID;
+  using Primitive = typename ATST::Primitive;
+  using AABB_traits = typename ATST::AABB_traits;
+  using AABB_tree = typename ATST::AABB_tree;
+
+protected:
+  double m_sq_alpha;
+
+  // one per face
+  DPPMB m_dppmb; // std::size_t (id) --> datum pos
+
+  // possibly more than one per face
+  BPM m_bpm; // std::size_t (id) --> bounding box
+  RPPM m_rppm; // std::size_t (id) --> reference point
+  DPMB m_dpmb; // std::size_t (datum pos) --> triangle datum
+
+  DPM m_dpm; // std::size_t (id) --> triangle (datum)
+
+  std::size_t fid = 0;
+
+public:
+  AABB_tree_oracle_splitter(const double alpha)
+    :
+      m_sq_alpha(square(alpha)),
+      m_dppmb(), m_bpm(), m_rppm(), m_dpmb(),
+      m_dpm(DPPM(m_dppmb/*first binder's value_map*/)/*second binder's key map*/, m_dpmb)
+  { }
+
+protected:
+  void initialize_tree_property_maps(const AABB_tree& tree) const
+  {
+    // Can't be set in the default constructed traits that are passed to the base
+    // since m_bpm is a member of the derived class.
+    //
+    // 'const_cast' because CGAL::AABB_tree only gives a const& to its traits.
+    const_cast<AABB_traits&>(tree.traits()).bbm = m_bpm;
+
+    const_cast<AABB_traits&>(tree.traits()).set_shared_data(m_dpm, m_rppm);
+  }
+
+  void reserve(std::size_t nf)
+  {
+    CGAL::internal::reserve(m_dpmb.range(), m_dpmb.range().size() + nf);
+
+    // Due to splitting, these might need more than 'nf'
+    CGAL::internal::reserve(m_dppmb.range(), m_dppmb.range().size() + nf);
+    CGAL::internal::reserve(m_rppm.value_map.range(), m_rppm.value_map.range().size() + nf);
+    CGAL::internal::reserve(m_bpm.value_map.range(), m_bpm.value_map.range().size() + nf);
+  }
+
+  template <typename P> // Kernel is Simple_Cartesian<Interval>
+  CGAL::Bbox_3 compute_bbox(const P& ap0, const P& ap1, const P& ap2)
+  {
+    double xmin = (CGAL::min)(ap0.x().inf(), (CGAL::min)(ap1.x().inf(), ap2.x().inf()));
+    double ymin = (CGAL::min)(ap0.y().inf(), (CGAL::min)(ap1.y().inf(), ap2.y().inf()));
+    double zmin = (CGAL::min)(ap0.z().inf(), (CGAL::min)(ap1.z().inf(), ap2.z().inf()));
+
+    double xmax = (CGAL::max)(ap0.x().sup(), (CGAL::max)(ap1.x().sup(), ap2.x().sup()));
+    double ymax = (CGAL::max)(ap0.y().sup(), (CGAL::max)(ap1.y().sup(), ap2.y().sup()));
+    double zmax = (CGAL::max)(ap0.z().sup(), (CGAL::max)(ap1.z().sup(), ap2.z().sup()));
+
+    return CGAL::Bbox_3(xmin, ymin, zmin, xmax, ymax, zmax);
+  }
+
+  template <typename P> // Kernel is Simple_Cartesian<Interval>
+  const P& compute_reference_point(const P&, const P&, const P& ap2)
+  {
+    return ap2; // ap2 is the midpoint when splitting
+  }
+
+  void split_and_insert_datum(const Triangle_3& tr,
+                              AABB_tree& tree,
+                              const GT& gt)
+  {
+    // Convert to intervals to ensure that the bounding box fully covers the subdividing triangle
+    using AK = CGAL::Simple_cartesian<CGAL::Interval_nt<true> >;
+    using K2AK = CGAL::Cartesian_converter<GT, AK>;
+    using AK2K = CGAL::Cartesian_converter<AK, GT>;
+    using AFT = AK::FT;
+    using APoint_3 = AK::Point_3;
+
+    using APL = std::pair<APoint_3, AFT>; // point and length of the opposite edge
+    using AT = std::array<APL, 3>;
+
+    const std::size_t data_size = m_dpmb.range().size();
+    put(m_dpmb, data_size, tr);
+
+    auto vertex = gt.construct_vertex_3_object();
+
+    const Point& p0 = vertex(tr, 0);
+    const Point& p1 = vertex(tr, 1);
+    const Point& p2 = vertex(tr, 2);
+
+    if(!subdivide || m_sq_alpha == 0.) // no splits
+    {
+      const std::size_t pid = tree.size();
+      ID id = std::make_pair(pid, fid++);
+
+      put(m_dppmb, pid, data_size);
+      put(m_rppm, id, p1); // the ref point that `One_point_from_face_descriptor_map` would give
+      put(m_bpm, id, gt.construct_bbox_3_object()(tr));
+
+//      std::cout << "Primitive[" << id.first << " " << id.second << "]; "
+//                    << "Bbox: [" << get(m_bpm, id) << "] "
+//                    << "Point: (" << get(m_rppm, id) << ") "
+//                    << "Datum: [" << get(m_bpm, id) << "]" << std::endl;
+
+      Primitive p(id/*, m_dpm, m_rppm*/); // pmaps are external, shared data
+      tree.insert(p);
+      return;
+    }
+
+    K2AK k2ak;
+    AK2K ak2k;
+    std::queue<AT> to_treat;
+
+    const APoint_3 ap0 = k2ak(p0);
+    const APoint_3 ap1 = k2ak(p1);
+    const APoint_3 ap2 = k2ak(p2);
+    const AFT sq_l0 = CGAL::squared_distance(ap1, ap2);
+    const AFT sq_l1 = CGAL::squared_distance(ap2, ap0);
+    const AFT sq_l2 = CGAL::squared_distance(ap0, ap1);
+
+    to_treat.push(CGAL::make_array(std::make_pair(ap0, sq_l0),
+                                   std::make_pair(ap1, sq_l1),
+                                   std::make_pair(ap2, sq_l2)));
+
+    while(!to_treat.empty())
+    {
+      const AT t = std::move(to_treat.front());
+      to_treat.pop();
+
+      const APL& apl0 = t[0];
+      const APL& apl1 = t[1];
+      const APL& apl2 = t[2];
+
+      int i = (apl0.second.sup() >= apl1.second.sup())
+                ? (apl0.second.sup() >= apl2.second.sup()) ? 0 : 2
+                : (apl1.second.sup() >= apl2.second.sup()) ? 1 : 2;
+
+      const FT max_sql = t[i].second.sup();
+
+      // The '3 * alpha' is some empirically-determined value
+      const FT sq_bound = 9. * m_sq_alpha;
+
+      // If the face is too big, do a fake split (two small bboxes rather than a big one)
+      if(max_sql > sq_bound)
+      {
+        // Could be factorized, but this is simpler to read
+        if(i == 0)
+        {
+          // 0 1 2 into 0 1 A and 0 A 2
+          const APoint_3 amp = CGAL::midpoint(apl1.first, apl2.first);
+          to_treat.push(CGAL::make_array(std::make_pair(apl0.first, CGAL::squared_distance(apl1.first, amp)),
+                                         std::make_pair(apl1.first, CGAL::squared_distance(amp, apl0.first)),
+                                         std::make_pair(amp, apl2.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl2.first, CGAL::squared_distance(apl0.first, amp)),
+                                         std::make_pair(apl0.first, CGAL::squared_distance(amp, apl2.first)),
+                                         std::make_pair(amp, apl1.second)));
+        }
+        else if(i == 1)
+        {
+          // 0 1 2 into 0 1 A and 1 2 A
+          const APoint_3 amp = CGAL::midpoint(apl2.first, apl0.first);
+          to_treat.push(CGAL::make_array(std::make_pair(apl0.first, CGAL::squared_distance(apl1.first, amp)),
+                                         std::make_pair(apl1.first, CGAL::squared_distance(amp, apl0.first)),
+                                         std::make_pair(amp, apl2.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl1.first, CGAL::squared_distance(apl2.first, amp)),
+                                         std::make_pair(apl2.first, CGAL::squared_distance(amp, apl1.first)),
+                                         std::make_pair(amp, apl0.second)));
+        }
+        else // i == 2
+        {
+          // 0 1 2 into 0 A 2 and 2 A 1
+          const APoint_3 amp = CGAL::midpoint(apl0.first, apl1.first);
+          to_treat.push(CGAL::make_array(std::make_pair(apl2.first, CGAL::squared_distance(apl0.first, amp)),
+                                         std::make_pair(apl0.first, CGAL::squared_distance(amp, apl2.first)),
+                                         std::make_pair(amp, apl1.second)));
+          to_treat.push(CGAL::make_array(std::make_pair(apl1.first, CGAL::squared_distance(apl2.first, amp)),
+                                         std::make_pair(apl2.first, CGAL::squared_distance(amp, apl1.first)),
+                                         std::make_pair(amp, apl0.second)));
+        }
+      }
+      else // all edges have length below the threshold, create a primitive
+      {
+        const std::size_t pid = tree.size();
+        ID id = std::make_pair(pid, fid);
+
+        put(m_dppmb, pid, data_size);
+        put(m_rppm, id, ak2k(compute_reference_point(apl0.first, apl1.first, apl2.first)));
+        put(m_bpm, id, compute_bbox(apl0.first, apl1.first, apl2.first));
+
+//        std::cout << "Primitive[" << id.first << " " << id.second << "]; "
+//                  << "Bbox: [" << get(m_bpm, id) << "] "
+//                  << "Point: (" << get(m_rppm, id) << ") "
+//                  << "Datum: [" << get(m_dpm, id) << "]" << std::endl;
+
+        Primitive p(id/*, m_dpm, m_rppm*/); // pmaps are external, shared data
+        tree.insert(p);
+      }
+    }
+
+    ++fid;
+  }
+};
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_INTERNAL_SPLITTING_HELPER_H

Alpha_wrap_3/include/CGAL/alpha_wrap_3.h

@@ -0,0 +1,417 @@
+// 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
+//                 Cedric Portaneri,
+//                 Mael Rouxel-Labbé
+//                 Andreas Fabri
+//                 Michael Hemmer
+//
+#ifndef CGAL_ALPHA_WRAP_3_H
+#define CGAL_ALPHA_WRAP_3_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h>
+
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/property_map.h>
+
+#include <boost/range/has_range_iterator.hpp>
+
+#include <type_traits>
+
+namespace CGAL {
+
+// /////////////////////////////////////////////////////////////////////////////////////////////////
+// WITH A TRIANGLE SOUP ---------------------------------------------------------------------------
+// /////////////////////////////////////////////////////////////////////////////////////////////////
+
+/*!
+* \ingroup AW3_free_functions_grp
+*
+* \brief computes a watertight, 2-manifold, and intersection-free triangulated surface mesh
+* that strictly contains an input triangle soup.
+*
+* The parameters `alpha` and `offset` respectively control which features will appear in the output,
+* and the distance from the input. See Section \ref aw3_parameters for a detailed breakdown of their influence.
+*
+* \tparam PointRange a model of `Range` whose value type is the point type
+* \tparam FaceRange a model of `RandomAccessContainer` whose value type is a model of `RandomAccessContainer` whose value type is an integral type
+* \tparam OutputMesh model of `MutableFaceGraph`.
+* \tparam InputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+* \tparam OutputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* \param points the input points
+* \param faces the input faces, with each element of the range being a range of indices corresponding to points in `points`
+* \param alpha the value of the parameter `alpha`
+* \param offset the value of the parameter `offset`
+* \param alpha_wrap the output surface mesh
+* \param in_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{a property map associating points to the elements of the point set `points`}
+*     \cgalParamType{a model of `ReadablePropertyMap` whose key type is the value type
+*                    of the iterator of `PointRange` and whose value type is `geom_traits::Point_3`}
+*     \cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{an instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
+*     \cgalParamExtra{<ul><li>The geometric traits class must be compatible with the point type.</li>
+*                         <li>The geometric traits should use a floating point number type (see \ref aw3_interface).</li></ul>}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* \param out_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{vertex_point_map}
+*     \cgalParamDescription{a property map associating points to the vertices of `alpha_wrap`}
+*     \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<OutputMesh>::%vertex_descriptor`
+*                    as key type and `%Point_3` as value type}
+*     \cgalParamDefault{`boost::get(CGAL::vertex_point, alpha_wrap)`}
+*     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+*                     must be available in `OutputMesh`.}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* \pre The elements of `faces` are triangles.
+* \pre `alpha` and `offset` are strictly positive values.
+*/
+template <typename PointRange, typename FaceRange, typename OutputMesh,
+          typename InputNamedParameters, typename OutputNamedParameters>
+void alpha_wrap_3(const PointRange& points,
+                  const FaceRange& faces,
+                  const double alpha,
+                  const double offset,
+                  OutputMesh& alpha_wrap,
+                  const InputNamedParameters& in_np,
+                  const OutputNamedParameters& out_np)
+{
+  using parameters::get_parameter;
+  using parameters::choose_parameter;
+
+  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>;
+
+  Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
+
+  Oracle oracle(alpha, gt);
+  oracle.add_triangle_soup(points, faces, in_np);
+  AW3 alpha_wrap_builder(oracle);
+  alpha_wrap_builder(alpha, offset, alpha_wrap, out_np);
+}
+
+// Convenience overloads
+template <typename PointRange, typename FaceRange, typename OutputMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+void alpha_wrap_3(const PointRange& points,
+                  const FaceRange& faces,
+                  const double alpha,
+                  const double offset,
+                  OutputMesh& alpha_wrap,
+                  const CGAL_NP_CLASS& in_np)
+{
+  return alpha_wrap_3(points, faces, alpha, offset, alpha_wrap, in_np, CGAL::parameters::default_values());
+}
+
+template <typename PointRange, typename FaceRange, typename OutputMesh>
+void alpha_wrap_3(const PointRange& points,
+                  const FaceRange& faces,
+                  const double alpha,
+                  const double offset,
+                  OutputMesh& alpha_wrap)
+{
+  return alpha_wrap_3(points, faces, alpha, offset, alpha_wrap, CGAL::parameters::default_values());
+}
+
+// without offset
+template <typename PointRange, typename FaceRange, typename OutputMesh,
+          typename T_I, typename Tag_I, typename Base_I,
+          typename T_O, typename Tag_O, typename Base_O>
+void alpha_wrap_3(const PointRange& points,
+                  const FaceRange& faces,
+                  const double alpha,
+                  OutputMesh& alpha_wrap,
+                  const CGAL::Named_function_parameters<T_I, Tag_I, Base_I>& in_np,
+                  const CGAL::Named_function_parameters<T_O, Tag_O, Base_O>& out_np,
+                  typename std::enable_if<boost::has_range_const_iterator<FaceRange>::value>::type* = nullptr)
+{
+  return alpha_wrap_3(points, faces, alpha, alpha / 30., alpha_wrap, in_np, out_np);
+}
+
+template <typename PointRange, typename FaceRange, typename OutputMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+void alpha_wrap_3(const PointRange& points,
+                  const FaceRange& faces,
+                  const double alpha,
+                  OutputMesh& alpha_wrap,
+                  const CGAL_NP_CLASS& in_np,
+                  typename std::enable_if<boost::has_range_const_iterator<FaceRange>::value>::type* = nullptr)
+{
+  return alpha_wrap_3(points, faces, alpha, alpha / 30., alpha_wrap, in_np,
+                      CGAL::parameters::default_values());
+}
+
+template <typename PointRange, typename FaceRange, typename OutputMesh>
+void alpha_wrap_3(const PointRange& points,
+                  const FaceRange& faces,
+                  const double alpha,
+                  OutputMesh& alpha_wrap,
+                  typename std::enable_if<boost::has_range_const_iterator<FaceRange>::value>::type* = nullptr)
+{
+  return alpha_wrap_3(points, faces, alpha, alpha / 30., alpha_wrap,
+                      CGAL::parameters::default_values(), CGAL::parameters::default_values());
+}
+
+// /////////////////////////////////////////////////////////////////////////////////////////////////
+// WITH A TRIANGLE MESH ----------------------------------------------------------------------------
+// /////////////////////////////////////////////////////////////////////////////////////////////////
+
+/*!
+* \ingroup AW3_free_functions_grp
+*
+* \brief computes a watertight, 2-manifold, and intersection-free triangulated surface mesh
+* that strictly contains an input triangle mesh.
+*
+* The parameters `alpha` and `offset` respectively control which features will appear in the output,
+* and the distance from the input. See Section \ref aw3_parameters for a detailed breakdown of their influence.
+*
+* \tparam TriangleMesh model of `FaceListGraph`.
+* \tparam OutputMesh model of `MutableFaceGraph`.
+* \tparam InputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+* \tparam OutputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* \param tmesh a triangle mesh
+* \param alpha the value of the parameter `alpha`
+* \param offset the value of the parameter `offset`
+* \param alpha_wrap the output surface mesh
+* \param in_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{vertex_point_map}
+*     \cgalParamDescription{a property map associating points to the vertices of `tmesh`}
+*     \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<TriangleMesh>::%vertex_descriptor`
+*                    as key type and `%Point_3` as value type}
+*     \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`}
+*     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+*                     must be available in `TriangleMesh`.}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{an instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
+*     \cgalParamExtra{<ul><li>The geometric traits class must be compatible with the point type.</li>
+*                         <li>The geometric traits should use a floating point number type (see \ref aw3_interface).</li></ul>}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* \param out_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{vertex_point_map}
+*     \cgalParamDescription{a property map associating points to the vertices of `alpha_wrap`}
+*     \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<OutputMesh>::%vertex_descriptor`
+*                    as key type and `%Point_3` as value type}
+*     \cgalParamDefault{`boost::get(CGAL::vertex_point, alpha_wrap)`}
+*     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+*                     must be available in `OutputMesh`.}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* \pre `tmesh` is a triangle mesh.
+* \pre `alpha` and `offset` are strictly positive values.
+*/
+template <typename TriangleMesh, typename OutputMesh,
+          typename InputNamedParameters, typename OutputNamedParameters>
+void alpha_wrap_3(const TriangleMesh& tmesh,
+                  const double alpha,
+                  const double offset,
+                  OutputMesh& alpha_wrap,
+                  const InputNamedParameters& in_np,
+                  const OutputNamedParameters& out_np
+#ifndef DOXYGEN_RUNNING
+                  , typename std::enable_if<! boost::has_range_const_iterator<TriangleMesh>::value>::type* = nullptr
+#endif
+                  )
+{
+  using parameters::get_parameter;
+  using parameters::choose_parameter;
+
+  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>;
+
+  Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
+
+  Oracle oracle(alpha, gt);
+  oracle.add_triangle_mesh(tmesh, in_np);
+  AW3 alpha_wrap_builder(oracle);
+  alpha_wrap_builder(alpha, offset, alpha_wrap, out_np);
+}
+
+// The convenience overloads are the same for triangle mesh & point set
+
+// /////////////////////////////////////////////////////////////////////////////////////////////////
+//  WITH A POINT SET -------------------------------------------------------------------------------
+// /////////////////////////////////////////////////////////////////////////////////////////////////
+
+/*!
+* \ingroup AW3_free_functions_grp
+*
+* \brief computes a watertight, 2-manifold, and intersection-free triangulated surface mesh
+* that strictly contains an input point set.
+*
+* The parameters `alpha` and `offset` respectively control which features will appear in the output,
+* and the distance from the input. See Section \ref aw3_parameters for a detailed breakdown of their influence.
+*
+* \tparam PointRange model of `Range` whose value type is a point type.
+* \tparam OutputMesh model of `MutableFaceGraph`.
+* \tparam InputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+* \tparam OutputNamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
+*
+* \param points the input points
+* \param alpha the value of the parameter `alpha`
+* \param offset the value of the parameter `offset`
+* \param alpha_wrap the output surface mesh
+* \param in_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{point_map}
+*     \cgalParamDescription{a property map associating points to the elements of the point range}
+*     \cgalParamType{a model of `ReadablePropertyMap` with value type `geom_traits::Point_3`}
+*     \cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
+*   \cgalParamNEnd
+*
+*   \cgalParamNBegin{geom_traits}
+*     \cgalParamDescription{an instance of a geometric traits class}
+*     \cgalParamType{a class model of `Kernel`}
+*     \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
+*     \cgalParamExtra{<ul><li>The geometric traits class must be compatible with the point type.</li>
+*                         <li>The geometric traits should use a floating point number type (see \ref aw3_interface).</li></ul>}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* \param out_np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{vertex_point_map}
+*     \cgalParamDescription{a property map associating points to the vertices of `alpha_wrap`}
+*     \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<OutputMesh>::%vertex_descriptor`
+*                    as key type and `%Point_3` as value type}
+*     \cgalParamDefault{`boost::get(CGAL::vertex_point, alpha_wrap)`}
+*     \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+*                     must be available in `OutputMesh`.}
+*   \cgalParamNEnd
+* \cgalNamedParamsEnd
+*
+* \pre `alpha` and `offset` are strictly positive values.
+*/
+template <typename PointRange, typename OutputMesh,
+#ifdef DOXYGEN_RUNNING
+          typename InputNamedParameters, typename OutputNamedParameters>
+#else
+          typename T_I, typename Tag_I, typename Base_I,
+          typename T_O, typename Tag_O, typename Base_O>
+#endif
+void alpha_wrap_3(const PointRange& points,
+                  const double alpha,
+                  const double offset,
+                  OutputMesh& alpha_wrap,
+#ifdef DOXYGEN_RUNNING
+                  const InputNamedParameters& in_np,
+                  const OutputNamedParameters& out_np
+#else
+                  const CGAL::Named_function_parameters<T_I, Tag_I, Base_I>& in_np,
+                  const CGAL::Named_function_parameters<T_O, Tag_O, Base_O>& out_np,
+                  typename std::enable_if<boost::has_range_const_iterator<PointRange>::value>::type* = nullptr
+#endif
+                  )
+{
+  using parameters::get_parameter;
+  using parameters::choose_parameter;
+
+  using InputNamedParameters = CGAL::Named_function_parameters<T_I, Tag_I, Base_I>;
+
+  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>;
+
+  Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(in_np, internal_np::geom_traits));
+
+  Oracle oracle(gt);
+  oracle.add_point_set(points, in_np);
+  AW3 alpha_wrap_builder(oracle);
+  alpha_wrap_builder(alpha, offset, alpha_wrap, out_np);
+}
+
+// Convenience overloads, common to both mesh and point set
+template <typename Input, typename OutputMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+void alpha_wrap_3(const Input& input,
+                  const double alpha,
+                  const double offset,
+                  OutputMesh& alpha_wrap,
+                  const CGAL_NP_CLASS& in_np)
+{
+  return alpha_wrap_3(input, alpha, offset, alpha_wrap, in_np, CGAL::parameters::default_values());
+}
+
+template <typename Input, typename OutputMesh>
+void alpha_wrap_3(const Input& input,
+                  const double alpha,
+                  const double offset,
+                  OutputMesh& alpha_wrap)
+{
+  return alpha_wrap_3(input, alpha, offset, alpha_wrap, CGAL::parameters::default_values());
+}
+
+// without offset
+template <typename Input, typename OutputMesh,
+          typename T_I, typename Tag_I, typename Base_I,
+          typename T_O, typename Tag_O, typename Base_O>
+void alpha_wrap_3(const Input& input,
+                  const double alpha,
+                  OutputMesh& alpha_wrap,
+                  const CGAL::Named_function_parameters<T_I, Tag_I, Base_I>& in_np,
+                  const CGAL::Named_function_parameters<T_O, Tag_O, Base_O>& out_np)
+{
+  return alpha_wrap_3(input, alpha, alpha / 30., alpha_wrap, in_np, out_np);
+}
+
+template <typename Input, typename OutputMesh, typename CGAL_NP_TEMPLATE_PARAMETERS>
+void alpha_wrap_3(const Input& input,
+                  const double alpha,
+                  OutputMesh& alpha_wrap,
+                  const CGAL_NP_CLASS& in_np)
+{
+  return alpha_wrap_3(input, alpha, alpha / 30., alpha_wrap, in_np, CGAL::parameters::default_values());
+}
+
+template <typename Input, typename OutputMesh>
+void alpha_wrap_3(const Input& input,
+                  const double alpha,
+                  OutputMesh& alpha_wrap)
+{
+  return alpha_wrap_3(input, alpha, alpha / 30., alpha_wrap,
+                      CGAL::parameters::default_values(), CGAL::parameters::default_values());
+}
+
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_H

Alpha_wrap_3/package_info/Alpha_wrap_3/copyright

@@ -0,0 +1 @@
+Google LLC (USA).

Alpha_wrap_3/package_info/Alpha_wrap_3/dependencies

@@ -0,0 +1,33 @@
+AABB_tree
+Algebraic_foundations
+Alpha_wrap_3
+Arithmetic_kernel
+BGL
+Box_intersection_d
+Cartesian_kernel
+Circulator
+Distance_2
+Distance_3
+Filtered_kernel
+Generator
+Hash_map
+Homogeneous_kernel
+Installation
+Intersections_2
+Intersections_3
+Interval_support
+Kernel_23
+Kernel_d
+Modular_arithmetic
+Number_types
+Polygon_mesh_processing
+Profiling_tools
+Property_map
+Random_numbers
+STL_Extension
+Spatial_searching
+Spatial_sorting
+Stream_support
+TDS_3
+Triangulation_3
+Union_find

Alpha_wrap_3/package_info/Alpha_wrap_3/description.txt

@@ -0,0 +1 @@
+The 3D Alpha Wrapping package is a generic meshing algorithm to create valid and input-enclosing meshes.

Alpha_wrap_3/package_info/Alpha_wrap_3/license.txt

@@ -0,0 +1 @@
+GPL (v3 or later)

Alpha_wrap_3/package_info/Alpha_wrap_3/long_description.txt

@@ -0,0 +1,7 @@
+This component takes a 3D triangle mesh, a triangle soup, or a point set as input, and generates
+a valid triangulated surface mesh that strictly contains the input (watertight, intersection-free
+and 2-manifold). The algorithm proceeds by shrink-wrapping and refining a 3D Delaunay triangulation
+loosely bounding the input. Two user-defined parameters, alpha and offset, offer control over the maximum size
+of cavities where the shrink-wrapping process can enter, and the tightness of the final surface mesh
+to the input, respectively. Once combined, these parameters provide a means to trade fidelity
+to the input for complexity of the output.

Alpha_wrap_3/package_info/Alpha_wrap_3/maintainer

@@ -0,0 +1 @@
+GeometryFactory

Alpha_wrap_3/test/Alpha_wrap_3/CMakeLists.txt

@@ -0,0 +1,14 @@
+# 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.23)
+project(Alpha_wrap_3_Tests)
+
+find_package(CGAL REQUIRED)
+
+# create a target per cppfile
+create_single_source_cgal_program("test_alpha_wrap_3_mesh.cpp")
+create_single_source_cgal_program("test_AW3_cavity_initializations.cpp")
+create_single_source_cgal_program("test_AW3_manifoldness.cpp")
+create_single_source_cgal_program("test_AW3_multiple_calls.cpp")
+create_single_source_cgal_program("test_AW3_compilation.cpp")

Alpha_wrap_3/test/Alpha_wrap_3/alpha_wrap_validation.h

@@ -0,0 +1,378 @@
+// 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)     : Cedric Portaneri,
+//                 Mael Rouxel-Labbé
+//
+#ifndef CGAL_ALPHA_WRAP_3_TEST_ALPHA_WRAP_VALIDATION_H
+#define CGAL_ALPHA_WRAP_3_TEST_ALPHA_WRAP_VALIDATION_H
+
+#include <CGAL/license/Alpha_wrap_3.h>
+
+#include <CGAL/boost/graph/helpers.h>
+#include <CGAL/Kernel_traits.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+#include <CGAL/Polygon_mesh_processing/intersection.h>
+#include <CGAL/Polygon_mesh_processing/manifoldness.h>
+#include <CGAL/Polygon_mesh_processing/measure.h>
+#include <CGAL/Polygon_mesh_processing/orientation.h>
+#include <CGAL/Polygon_mesh_processing/repair.h>
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Side_of_triangle_mesh.h>
+
+namespace CGAL {
+namespace Alpha_wraps_3 {
+namespace internal {
+
+// @todo this does not detect non-manifold edges with manifold vertices
+// (but PMP::self_intersections will catch it)
+template <typename TriangleMesh>
+bool is_combinatorially_non_manifold(const TriangleMesh& mesh)
+{
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  using halfedge_descriptor = typename boost::graph_traits<TriangleMesh>::halfedge_descriptor;
+
+  std::vector<halfedge_descriptor> non_manifold_cones;
+  PMP::non_manifold_vertices(mesh, std::back_inserter(non_manifold_cones));
+
+  if(!non_manifold_cones.empty())
+    return true;
+
+  return false;
+}
+
+template <typename TriangleMesh,
+          typename NamedParameters = parameters::Default_named_parameters>
+bool has_degenerated_faces(const TriangleMesh& mesh,
+                           const NamedParameters& np = CGAL::parameters::default_values())
+{
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  for(auto f : faces(mesh))
+    if(PMP::is_degenerate_triangle_face(f, mesh, np))
+      return true;
+
+  return false;
+}
+
+// Edge length is bounded by twice the circumradius
+template <typename TriangleMesh,
+          typename NamedParameters = parameters::Default_named_parameters>
+bool check_edge_length(const TriangleMesh& output_mesh,
+                       const double alpha,
+                       const NamedParameters& np = CGAL::parameters::default_values())
+{
+  const auto sq_alpha_bound = 4 * square(alpha);
+  for(auto e : edges(output_mesh))
+  {
+    const auto sqd = Polygon_mesh_processing::squared_edge_length(e, output_mesh, np);
+    if(sqd > sq_alpha_bound) // alpha is the circumradius
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cerr << "Error: " << sqd << " greater than " << sq_alpha_bound << std::endl;
+      std::cerr << get(CGAL::vertex_point, output_mesh, source(e, output_mesh)) << std::endl;
+      std::cerr << get(CGAL::vertex_point, output_mesh, target(e, output_mesh)) << std::endl;
+#endif
+      return false;
+    }
+  }
+
+  return true;
+}
+
+template <typename ConcurrencyTag = CGAL::Sequential_tag,
+          typename TriangleMesh, typename FT,
+          typename InputNamedParameters = parameters::Default_named_parameters,
+          typename OutputNamedParameters = parameters::Default_named_parameters>
+bool has_expected_Hausdorff_distance(const TriangleMesh& wrap,
+                                     const TriangleMesh& input,
+                                     const FT alpha, const FT offset,
+                                     const InputNamedParameters& in_np = parameters::default_values(),
+                                     const OutputNamedParameters& out_np = parameters::default_values())
+{
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  using face_descriptor = typename boost::graph_traits<TriangleMesh>::face_descriptor;
+
+  std::vector<std::pair<face_descriptor, face_descriptor> > fpairs;
+
+  const FT bound = 0.01 * (std::min)(alpha, offset);
+  const FT d = PMP::bounded_error_Hausdorff_distance<ConcurrencyTag>(
+                 wrap, input, bound, in_np.output_iterator(std::back_inserter(fpairs)), out_np);
+
+#ifdef CGAL_AW3_DEBUG
+  std::cout << "Alpha: " << alpha << " Offset: " << offset << " Bound: " << bound << " Hausdorff_distance: " << d << std::endl;
+  std::cout << "Maximum distance on faces " << fpairs.back().first << " " << fpairs.back().second << std::endl;
+#endif
+
+  return (d < (alpha + offset + bound));
+}
+
+template <typename TriangleMesh, typename NamedParameters = parameters::Default_named_parameters>
+bool is_valid_wrap(const TriangleMesh& wrap,
+                   const bool check_manifoldness = true,
+                   const NamedParameters& np = parameters::default_values())
+{
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  if(is_empty(wrap))
+  {
+#ifdef CGAL_AW3_DEBUG
+    std::cerr << "Error: empty wrap" << std::endl;
+#endif
+    return false;
+  }
+
+  if(!is_valid_polygon_mesh(wrap))
+  {
+#ifdef CGAL_AW3_DEBUG
+    std::cerr << "Error: Invalid wrap mesh" << std::endl;
+#endif
+    return false;
+  }
+
+  if(!is_triangle_mesh(wrap))
+  {
+#ifdef CGAL_AW3_DEBUG
+    std::cerr << "Error: Wrap is not triangulated" << std::endl;
+#endif
+    return false;
+  }
+
+  if(!is_closed(wrap))
+  {
+    if(check_manifoldness)
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cerr << "Error: Wrap is not closed" << std::endl;
+#endif
+      return false;
+    }
+    else
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cerr << "W: Wrap is not closed" << std::endl;
+#endif
+    }
+  }
+  else if(!PMP::does_bound_a_volume(wrap, np))
+  {
+#ifdef CGAL_AW3_DEBUG
+    std::cerr << "Error: Wrap does not bound a volume" << std::endl;
+#endif
+    return false;
+  }
+
+  if(has_degenerated_faces(wrap, np))
+  {
+#ifdef CGAL_AW3_DEBUG
+    std::cerr << "Error: Wrap has degenerate faces" << std::endl;
+#endif
+    return false;
+  }
+
+  if(is_combinatorially_non_manifold(wrap))
+  {
+#ifdef CGAL_AW3_DEBUG
+    std::cerr << "Error: Wrap is combinatorially non-manifold" << std::endl;
+#endif
+    return false;
+  }
+
+  if(PMP::does_self_intersect(wrap, np))
+  {
+    if(check_manifoldness)
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cerr << "Error: Wrap self-intersects" << std::endl;
+#endif
+      return false;
+    }
+#ifdef CGAL_AW3_DEBUG
+    std::cerr << "W: Wrap self-intersects" << std::endl;
+#endif
+  }
+
+  return true;
+}
+
+template <typename InputTriangleMesh, typename OutputTriangleMesh,
+          typename InputNamedParameters = parameters::Default_named_parameters,
+          typename OutputNamedParameters = parameters::Default_named_parameters>
+bool is_outer_wrap_of_triangle_mesh(const OutputTriangleMesh& wrap,
+                                    const InputTriangleMesh& input,
+                                    const OutputNamedParameters& out_np = parameters::default_values(),
+                                    const InputNamedParameters& in_np = parameters::default_values())
+{
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  using parameters::get_parameter;
+  using parameters::choose_parameter;
+
+  using IVPM = typename GetVertexPointMap<InputTriangleMesh, InputNamedParameters>::const_type;
+  using OVPM = typename GetVertexPointMap<OutputTriangleMesh, OutputNamedParameters>::const_type;
+  using K = typename GetGeomTraits<InputTriangleMesh, InputNamedParameters>::type;
+
+//  CGAL::Rigid_triangle_mesh_collision_detection<TriangleMesh> collision_detection;
+//  collision_detection.add_mesh(input);
+//  collision_detection.add_mesh(wrap);
+
+//  auto res = collision_detection.get_all_intersections(1);
+//  if(res.size() != 0)
+//  {
+//    std::cerr << "Error: The wrap intersects the input mesh" << std::endl;
+//    return EXIT_FAILURE;
+//  }
+
+  if(PMP::do_intersect(input, wrap, in_np, out_np))
+  {
+#ifdef CGAL_AW3_DEBUG
+    std::cerr << "Error: The wrap intersects the input mesh" << std::endl;
+#endif
+    return false;
+  }
+
+  IVPM in_vpm = choose_parameter(get_parameter(in_np, internal_np::vertex_point),
+                                 get_const_property_map(vertex_point, input));
+  OVPM out_vpm = choose_parameter(get_parameter(out_np, internal_np::vertex_point),
+                                  get_const_property_map(vertex_point, wrap));
+
+  CGAL::Side_of_triangle_mesh<OutputTriangleMesh, K, OVPM> side_of_wrap(wrap, out_vpm);
+
+  // @speed a single vertex per CC would be sufficient
+  for(auto v : vertices(input))
+  {
+    if(side_of_wrap(get(in_vpm, v)) != CGAL::ON_BOUNDED_SIDE)
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cerr << "Error: Part(s) of the input mesh are outside the wrap: " << get(in_vpm, v) << std::endl;
+#endif
+      return false;
+    }
+  }
+
+  return true;
+}
+
+template <typename OutputTriangleMesh, typename InputTriangleMesh,
+          typename OutputNamedParameters = parameters::Default_named_parameters,
+          typename InputNamedParameters = parameters::Default_named_parameters>
+bool is_valid_wrap_of_triangle_mesh(const OutputTriangleMesh& wrap,
+                                    const InputTriangleMesh& input,
+                                    const OutputNamedParameters& out_np = parameters::default_values(),
+                                    const InputNamedParameters& in_np = parameters::default_values())
+{
+  if(!is_valid_wrap(wrap, out_np))
+    return false;
+
+  if(!is_outer_wrap_of_triangle_mesh(wrap, input, out_np, in_np))
+    return false;
+
+  return true;
+}
+
+template <typename TriangleMesh, typename PointRange, typename FaceRange,
+          typename OutputNamedParameters = parameters::Default_named_parameters,
+          typename InputNamedParameters = parameters::Default_named_parameters>
+bool is_outer_wrap_of_triangle_soup(const TriangleMesh& wrap,
+                                    PointRange points, // intentional copies
+                                    FaceRange faces,
+                                    const OutputNamedParameters& out_np = parameters::default_values(),
+                                    const InputNamedParameters& in_np = parameters::default_values())
+{
+  namespace PMP = CGAL::Polygon_mesh_processing;
+
+  // Make a mesh out of the soup
+  PMP::repair_polygon_soup(points, faces);
+  PMP::orient_polygon_soup(points, faces);
+  CGAL_assertion(PMP::is_polygon_soup_a_polygon_mesh(faces));
+
+  TriangleMesh mesh;
+  PMP::polygon_soup_to_polygon_mesh(points, faces, mesh, in_np);
+
+  return is_outer_wrap_of_triangle_mesh(wrap, mesh, out_np);
+}
+
+template <typename TriangleMesh, typename PointRange, typename FaceRange,
+          typename OutputNamedParameters = parameters::Default_named_parameters,
+          typename InputNamedParameters = parameters::Default_named_parameters>
+bool is_valid_wrap_of_triangle_soup(const TriangleMesh& wrap,
+                                    const PointRange& points,
+                                    const FaceRange& faces,
+                                    const OutputNamedParameters& out_np = parameters::default_values(),
+                                    const InputNamedParameters& in_np = parameters::default_values())
+{
+  if(!is_valid_wrap(wrap, out_np))
+    return false;
+
+  if(!is_outer_wrap_of_triangle_soup(wrap, points, faces, out_np, in_np))
+    return false;
+
+  return true;
+}
+
+template <typename TriangleMesh, typename PointRange,
+          typename OutputNamedParameters = parameters::Default_named_parameters,
+          typename InputNamedParameters = parameters::Default_named_parameters>
+bool is_outer_wrap_of_point_set(const TriangleMesh& wrap,
+                                const PointRange& points,
+                                const OutputNamedParameters& out_np = parameters::default_values(),
+                                const InputNamedParameters& in_np = parameters::default_values())
+{
+  using parameters::get_parameter;
+  using parameters::choose_parameter;
+
+  using OVPM = typename GetVertexPointMap<TriangleMesh, OutputNamedParameters>::const_type;
+  using IPM = typename GetPointMap<PointRange, InputNamedParameters>::const_type;
+  using K = typename Kernel_traits<typename boost::property_traits<IPM>::value_type>::Kernel;
+
+  OVPM out_vpm = choose_parameter(get_parameter(out_np, internal_np::vertex_point),
+                                  get_const_property_map(vertex_point, wrap));
+  IPM in_pm = choose_parameter<IPM>(get_parameter(in_np, internal_np::point_map));
+
+  CGAL::Side_of_triangle_mesh<TriangleMesh, K, OVPM> side_of_wrap(wrap, out_vpm);
+
+  // @speed a single vertex per CC would be sufficient
+  for(const auto& p : points)
+  {
+    if(side_of_wrap(get(in_pm, p)) != CGAL::ON_BOUNDED_SIDE)
+    {
+#ifdef CGAL_AW3_DEBUG
+      std::cerr << "Part(s) of the input mesh are outside the wrap: " << get(in_pm, p) << std::endl;
+#endif
+      return false;
+    }
+  }
+
+  return true;
+}
+
+template <typename TriangleMesh, typename PointRange,
+          typename OutputNamedParameters = parameters::Default_named_parameters,
+          typename InputNamedParameters = parameters::Default_named_parameters>
+bool is_valid_wrap_of_point_set(const TriangleMesh& wrap,
+                                const PointRange& points,
+                                const OutputNamedParameters& out_np = parameters::default_values(),
+                                const InputNamedParameters& in_np = parameters::default_values())
+{
+  if(!is_valid_wrap(wrap, out_np))
+    return false;
+
+  if(!is_outer_wrap_of_point_set(wrap, points, out_np, in_np))
+    return false;
+
+  return true;
+}
+
+} // namespace internal
+} // namespace Alpha_wraps_3
+} // namespace CGAL
+
+#endif // CGAL_ALPHA_WRAP_3_TEST_ALPHA_WRAP_VALIDATION_H

Alpha_wrap_3/test/Alpha_wrap_3/data/combinatorial_manifold.off

@@ -0,0 +1,84 @@
+OFF
+28 54 0
+0.553363 0.305 0 
+0 0 0 
+0 0.61 0 
+0 0.305 0.5901607 
+0.553363 0.5409906 0.7191805 
+0.553363 0.07099058 0.7191805 
+-0.1135536 0.007559509 0.07658271 
+-0.1997082 0.0180451 0.1828085 
+-0.250964 0.03054399 0.3094303 
+-0.2628592 0.04396815 0.4454258 
+-0.2343584 0.05714902 0.5789564 
+-0.1679425 0.06893919 0.6983985 
+-0.06939305 0.07831232 0.7933543 
+0.4494582 0.07978794 0.8083033 
+0.05271121 0.08445249 0.8555582 
+0.3239422 0.08522387 0.8633728 
+0.1877411 0.08682518 0.8795952 
+-0.09665208 0.6097971 0.08148474 
+-0.1673462 0.6069757 0.1862502 
+-0.206674 0.6017519 0.3062813 
+-0.2116265 0.594525 0.432395 
+-0.181825 0.5858482 0.554943 
+-0.1195494 0.5763851 0.6645497 
+-0.02956408 0.5668598 0.7528296 
+0.4498522 0.541643 0.7917513 
+0.08124659 0.5580011 0.8130288 
+0.3304889 0.5448912 0.8332637 
+0.204405 0.5504866 0.8405418 
+3 0 1 2
+3 3 0 2
+3 3 1 0
+3 4 0 3
+3 5 3 0
+3 4 5 0
+3 5 6 1
+3 6 5 7
+3 7 5 8
+3 8 5 9
+3 9 5 10
+3 10 5 11
+3 11 5 12
+3 12 5 13
+3 12 13 14
+3 14 13 15
+3 14 15 16
+3 17 4 2
+3 4 17 18
+3 4 18 19
+3 4 19 20
+3 4 20 21
+3 4 21 22
+3 4 22 23
+3 4 23 24
+3 24 23 25
+3 24 25 26
+3 26 25 27
+3 5 1 3
+3 4 3 2
+3 24 26 13
+3 4 24 5
+3 5 24 13
+3 13 26 15
+3 26 27 15
+3 15 27 16
+3 27 25 16
+3 16 25 14
+3 25 23 14
+3 14 23 12
+3 23 22 12
+3 12 22 11
+3 22 10 11
+3 21 10 22
+3 21 9 10
+3 20 9 21
+3 19 9 20
+3 19 8 9
+3 18 8 19
+3 18 7 8
+3 17 7 18
+3 17 6 7
+3 1 6 2
+3 2 6 17

Alpha_wrap_3/test/Alpha_wrap_3/data/combinatorial_manifold_multiple_components.off

@@ -0,0 +1,16 @@
+OFF
+6 8 0
+0 0.61 0 
+0 0 0 
+0 0.305 0.5901607 
+0.553363 0.305 0 
+0.553363 0.5409906 0.7191805 
+0.553363 0.07099058 0.7191805 
+3 0 1 2
+3 3 0 1
+3 2 3 0
+3 2 1 3
+3 4 2 5
+3 4 3 2
+3 5 2 3
+3 4 5 3

Alpha_wrap_3/test/Alpha_wrap_3/data/non_manifold.off

@@ -0,0 +1,21 @@
+OFF
+7 8 0
+
+0 0 0
+1 1 0
+1 0 1
+0 1 1
+
+1 1 10
+1 0 11
+0 1 11
+
+3 0 1 2
+3 3 1 0 
+3 3 2 1
+3 3 0 2
+
+3 0 4 5
+3 6 4 0 
+3 6 5 4
+3 6 0 5

Alpha_wrap_3/test/Alpha_wrap_3/data/open_tetrahedron.off

@@ -0,0 +1,9 @@
+OFF
+4 3 0
+0 0 0 
+1 1 0 
+1 0 1 
+0 1 1 
+3 0 1 2
+3 3 1 0
+3 3 2 1

Alpha_wrap_3/test/Alpha_wrap_3/data/plane.off

@@ -0,0 +1,8 @@
+OFF
+4 2 0
+0 0 0
+1 1 0
+1 0 0
+0 1 0
+3 0 1 2
+3 3 1 0