Merge pull request #7672 from JacksonCampolattaro/orthtree-generalization

Orthtree generalization

Installation/CHANGES.md

@@ -44,6 +44,12 @@ Release date: October 2023
 -   Removed the class templates `Gray_image_mesh_domain_3`, `Implicit_mesh_domain_3`, and `Labeled_image_mesh_domain_3`
     which are deprecated since CGAL-4.13.
 
+### [Quadtrees, Octrees, and Orthtrees](https://doc.cgal.org/6.0/Manual/packages.html#PkgOrthtree)
+- **Breaking change**:
+  - Node splitting behavior and per-node data are now customizable via the Traits class.
+  - Nodes are now stored as a property map, with properties of each node accessed by index.
+  - Nearest neighbors functions only work for Orthtrees which provide the necessary functionality.
+
 ### [Polygon Mesh Processing](https://doc.cgal.org/6.0/Manual/packages.html#PkgPolygonMeshProcessing)
 
 -   Added the function `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures()` which can be used to compute

Orthtree/benchmark/Orthtree/construction.cpp

@@ -12,16 +12,14 @@
 #include <iostream>
 #include <chrono>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef Point_set::Point_map Point_map;
-
-typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
-
-typedef CGAL::Search_traits_3<Kernel> Kd_tree_traits;
-typedef CGAL::Orthogonal_k_neighbor_search<Kd_tree_traits> Kd_tree_search;
-typedef Kd_tree_search::Tree Kdtree;
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Point_map = Point_set::Point_map;
+using Octree = CGAL::Octree<Kernel, Point_set, Point_map>;
+using Kd_tree_traits = CGAL::Search_traits_3<Kernel>;
+using Kd_tree_search = CGAL::Orthogonal_k_neighbor_search<Kd_tree_traits>;
+using Kdtree = Kd_tree_search::Tree;
 
 using std::chrono::high_resolution_clock;
 using std::chrono::duration_cast;

Orthtree/benchmark/Orthtree/nearest_neighbor.cpp

@@ -12,21 +12,19 @@
 #include <iostream>
 #include <chrono>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef Point_set::Point_map Point_map;
-
-typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
-
-typedef CGAL::Search_traits_3<Kernel> Kd_tree_traits;
-typedef CGAL::Orthogonal_k_neighbor_search<Kd_tree_traits> Kd_tree_search;
-typedef Kd_tree_search::Tree Kdtree;
-
 using std::chrono::high_resolution_clock;
 using std::chrono::duration_cast;
 using std::chrono::microseconds;
 
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Point_map = Point_set::Point_map;
+using Octree = CGAL::Octree<Kernel, Point_set, Point_map>;
+using Kd_tree_traits = CGAL::Search_traits_3<Kernel>;
+using Kd_tree_search = CGAL::Orthogonal_k_neighbor_search<Kd_tree_traits>;
+using Kdtree = Kd_tree_search::Tree;
+
 int main(int argc, char **argv) {
 
   int num_runs = 100;
@@ -103,8 +101,8 @@ int main(int argc, char **argv) {
       // Time how long it takes to find neighbors using the octree
       auto octreeTime = bench<microseconds>(
               [&] {
-                std::vector<Point> nearest_neighbors;
-                octree.nearest_neighbors(search_point, k, std::back_inserter(nearest_neighbors));
+                std::vector<Point_set::Index> nearest_neighbors;
+                octree.nearest_k_neighbors(search_point, k, std::back_inserter(nearest_neighbors));
               }
       );
 

Orthtree/benchmark/Orthtree/util.h

@@ -10,8 +10,8 @@
 template<class Kernel>
 CGAL::Point_set_3<typename Kernel::Point_3> generate(size_t num_points = 1) {
 
-  typedef typename Kernel::Point_3 Point;
-  typedef CGAL::Point_set_3<Point> Point_set;
+  using Point = typename Kernel::Point_3;
+  using Point_set = CGAL::Point_set_3<Point>;
 
   // Create an empty point set
   Point_set points;

Orthtree/doc/Orthtree/Concepts/CollectionPartitioningOrthtreeTraits.h

@@ -0,0 +1,101 @@
+/*!
+  \ingroup PkgOrthtreeConcepts
+  \cgalConcept
+
+  Refinement of the `OrthtreeTraitsWithData` concept, adding requirements for the
+  traits class of a `CGAL::Orthtree` in order to supports nearest-neighbor searching.
+
+  Nearest neighbor searches expect a tree where `Node_data` is a model of `ForwardRange`.
+  The leaf nodes of the tree represent an exclusive partition of the elements contained in the tree.
+  This means that no element should be contained by more than one node.
+
+  \cgalRefines{OrthtreeTraitsWithData}
+
+  \cgalHasModelsBegin
+  \cgalHasModels{CGAL::Orthtree_traits_point<GeomTraits, PointRange, PointMap, DimensionTag>}
+  \cgalHasModelsEnd
+*/
+class CollectionPartitioningOrthtreeTraits {
+public:
+
+
+  /// \name Types
+  /// @{
+
+  /*!
+   * Sphere Type used for the shrinking-sphere approach for neighbor queries; needs to be copy assignable.
+   */
+  using Sphere_d = unspecified_type;
+
+  /*!
+   * \brief The data type contained by each node; must be a model of `ForwardRange` in addition to default constructible, copy constructible and copy assignable.
+   */
+  using Node_data = unspecified_type;
+
+  /*!
+   * \brief An element of the `Node_data` list-like type.
+   *
+   * Must be constructible from the value type of a `Node_data::iterator`.
+   * Typically the same as that type, but can also be an `std::reference_wrapper<>` if the type is not copyable.
+   */
+  using Node_data_element = unspecified_type;
+
+  /*!
+   * \brief Functor with an operator that calculates the squared distance of a `Node_data_element` from a point.
+   *
+   * Provides the operator:
+   * `FT operator()(const Node_data_element&, const Point_d&)`
+   */
+  using Squared_distance_of_element = unspecified_type;
+
+  /*!
+   * \brief Functor with an operator that constructs a `Sphere_d` from a provided center and squared radius.
+   *
+   * Provides the operator:
+   * `Sphere_d operator()(const Point_d&, const FT&)`
+  */
+  using Construct_sphere_d = unspecified_type;
+
+  /*!
+   * \brief Functor with an operator that provides the center of a `Sphere_d`.
+   *
+   * Provides the operator:
+   * `Point_d operator()(const Sphere_d&)`
+  */
+  using Construct_center_d = unspecified_type;
+
+  /*!
+   * \brief Functor with an operator that provides the squared radius of a `Sphere_d`.
+   *
+   * Provides the operator:
+   * `FT operator()(const Sphere_d&)`
+  */
+  using Compute_squared_radius_d = unspecified_type;
+
+  /// @}
+
+  /// \name Operations
+  /// @{
+
+  /*!
+   * constructs an object of type `ConstructSphere_d`.
+   */
+  Construct_sphere_d construct_sphere_d_object() const;
+
+  /*!
+   * constructs an object of type `ConstructCenter_d`.
+   */
+  Construct_center_d construct_center_d_object() const;
+
+  /*!
+   * constructs an object of type `ComputeSquaredRadius_d`.
+   */
+  Compute_squared_radius_d compute_squared_radius_d_object() const;
+
+  /*!
+   * constructs an object of type `Squared_distance_of_element`.
+   */
+  Squared_distance_of_element squared_distance_of_element_object() const;
+
+  /// @}
+};

Orthtree/doc/Orthtree/Concepts/OrthtreeTraits.h

@@ -6,9 +6,9 @@
   template parameter of the `CGAL::Orthtree` class.
 
   \cgalHasModelsBegin
-  \cgalHasModels{CGAL::Orthtree_traits_2<GeomTraits>}
-  \cgalHasModels{CGAL::Orthtree_traits_3<GeomTraits>}
-  \cgalHasModels{CGAL::Orthtree_traits_d<GeomTraits,Dimension>}
+  \cgalHasModels{CGAL::Orthtree_traits_point<GeomTraits, PointRange, PointMap, dimension>}
+  \cgalHasModels{CGAL::Orthtree_traits_face_graph<PolygonMesh, VPM>}
+  \cgalHasModels{CGAL::Orthtree_traits_base<K, dimension>}
   \cgalHasModelsEnd
 */
 class OrthtreeTraits
@@ -17,31 +17,47 @@ public:
 
   /// \name Types
   /// @{
-
-  typedef unspecified_type Dimension; ///< Dimension type (see `CGAL::Dimension_tag`).
-  typedef unspecified_type Bbox_d; ///< Bounding box type.
-  typedef unspecified_type FT; ///< The number type of the %Cartesian coordinates of types `Point_d`
-  typedef unspecified_type Point_d; ///< Point type.
-  typedef unspecified_type Sphere_d; ///< The sphere type for neighbor queries.
+  using Node_index = unspecified_type; ///< An integer type for nodes
+  constexpr int dimension; ///< Dimension.
+  using FT = unspecified_type; ///< The number type of the %Cartesian coordinates of types `Point_d`
+  using Point_d = unspecified_type; ///< Point type.
+  using Bbox_d = unspecified_type; ///< Bounding box type. Must be constructible from a pair of `Point_d` objects.
 
   /*!
     A random access iterator type to enumerate the
-    %Cartesian coordinates of a point.
+    %Cartesian coordinates of a point of type `Point_d`.
   */
-  typedef unspecified_type Cartesian_const_iterator_d;
-  typedef std::array<FT, Dimension::value> Array; ///< Array used for easy point constructions.
-
-  typedef unspecified_type Adjacency; ///< Specify the adjacency directions
+  using Cartesian_const_iterator_d = unspecified_type;
 
   /*!
-    Functor with an operator to construct a `Point_d` from an `Array` object.
-  */
-  typedef unspecified_type Construct_point_d_from_array;
+   * \brief Integral number type which can take on values indicating adjacency directions.
+   *
+   * Must be able to take on values ranging from 0 to the number of faces of the (hyper)rectangle type, equivalent to 2 * D.
+   */
+  using Adjacency = unspecified_type; ///< Specify the adjacency directions
 
   /*!
-    Functor with an operator to construct a `Bbox_d` from two `Array` objects (coordinates of minimum and maximum points).
-  */
-  typedef unspecified_type Construct_bbox_d;
+   * \brief Functor with an operator to create the bounding box of the root node.
+   *
+   * Provides the operator:
+   * `Bbox_d operator()()`
+   *
+   * The bounding box must enclose all elements contained by the tree.
+   * It may be tight-fitting. The orthtree constructor produces a bounding box surrounding this region.
+   * For traits which assign no data to each node, this can be defined to return a fixed region.
+   */
+  using Construct_root_node_bbox = unspecified_type;
+
+  /*!
+   * \brief Functor with an operator to construct a `Point_d` from an initializer list of type `FT`.
+   *
+   * Provides the operator:
+   * `Point_d operator()(arg1, arg2,...)`
+   *
+   * For trees which use a different kernel for the bounding box type,
+   * the return type of this functor must match the kernel used by the bounding box type and not that of the contents.
+   */
+  using Construct_point_d = unspecified_type;
 
   /// @}
 
@@ -49,14 +65,14 @@ public:
   /// @{
 
   /*!
-    Function used to construct an object of type `Construct_point_d_from_array`.
-  */
-  Construct_point_d_from_array construct_point_d_from_array_object() const;
+   * constructs an object of type `Construct_root_node_bbox`.
+   */
+  Construct_root_node_bbox construct_root_node_bbox_object() const;
 
   /*!
-    Function used to construct an object of type `Construct_bbox_d`.
-  */
-  Construct_bbox_d construct_bbox_d_object() const;
+   * constructs an object of type `Construct_point_d`.
+   */
+  Construct_point_d construct_point_d_object() const;
 
   /// @}
 };

Orthtree/doc/Orthtree/Concepts/OrthtreeTraitsWithData.h

@@ -0,0 +1,76 @@
+/*!
+  \ingroup PkgOrthtreeConcepts
+  \cgalConcept
+
+  The concept `OrthtreeTraitsWithData` defines the requirements for the
+  template parameter of the `CGAL::Orthtree` class for a node type that stores data.
+
+  \cgalRefines{OrthtreeTraits}
+
+  \cgalHasModelsBegin
+  \cgalHasModels{CGAL::Orthtree_traits_point<GeomTraits, PointRange, PointMap, dimension>}
+  \cgalHasModels{CGAL::Orthtree_traits_face_graph<PolygonMesh, VPM>}
+  \cgalHasModels{CGAL::Orthtree_traits_base<K, dimension>}
+  \cgalHasModelsEnd
+*/
+class OrthtreeTraitsWithData
+{
+public:
+
+  /// \name Types
+  /// @{
+  /*!
+   * \brief The data type contained by each node. Must be default constructible, copy constructible and copy assignable.
+   */
+  using Node_data = unspecified_type;
+
+  /*!
+   * \brief Functor which initializes elements contained by the root node.
+   *
+   * Each node of a tree has an associated `Node_data` value.
+   * This functor initializes the `Node_data` of the root node.
+   * It takes no arguments, and returns an instance of `Node_data`.
+   *
+   * Provides the operator:
+   * `Node_data operator()()`
+   *
+   * Typically, the `Node_data` of the root node contains all the elements in the tree.
+   * For a tree in which each node contains a span (such as `std::span()`) this function would return the span containing all items.
+   *
+   */
+  using Construct_root_node_contents = unspecified_type;
+
+  /*!
+   * \brief Functor which fills the contents of the nodes children.
+   *
+   * Provides the operator:
+   * `void operator()(Node_index, Orthtree<Traits>&, const Point_d&)`
+   *
+   * The functor is called during refinement of the `Orthtree` on a node after it has been split. The purpose of the functor is to
+   * distribute the `Node_data`, accessible via `tree.data()`, to the data of the nodes children, accessible via `tree.children()`.
+   * The first parameter is the `Node_index` of the node. The second parameter provides the instance of the `Orthtree`
+   * and the last parameter is the barycenter of the node which will be used as shared corner amongst the children of the node.
+   *
+   * For a tree in which each node contains a span, this may mean rearranging the contents of the original node
+   * and producing spans containing a subset of its contents for each of its children.
+   * For compatibility with locate, the center of the node is considered to be part of the upper half.
+   */
+  using Distribute_node_contents = unspecified_type;
+
+  /// @}
+
+  /// \name Operations
+  /// @{
+
+  /*!
+   * constructs an object of type `Construct_root_node_contents`.
+   */
+  Construct_root_node_contents construct_root_node_contents_object() const;
+
+  /*!
+   * constructs an object of type `Distribute_node_contents`.
+   */
+  Distribute_node_contents distribute_node_contents_object() const;
+
+  /// @}
+};

Orthtree/doc/Orthtree/Concepts/OrthtreeTraversal.h

@@ -1,13 +1,8 @@
-
 /*!
   \ingroup PkgOrthtreeConcepts
   \cgalConcept
 
-  \brief a traversal provides the functions needed to traverse the
-  nodes of an orthtree.
-
-  A traversal is used to iterate on a tree with a user-selected order
-  (e.g. preorder, postorder).
+  \brief Requirements for defining a traversal strategy of an orthtree.
 
   \cgalHasModelsBegin
   \cgalHasModels{CGAL::Orthtrees::Preorder_traversal}
@@ -20,15 +15,15 @@ class OrthtreeTraversal {
 
 public:
 
-  using Node = unspecified_type; ///< A specialization of [CGAL::Orthtree<Traits,PointRange,PointMap>::Node](@ref CGAL::Orthtree::Node)
+  using Node_index = unspecified_type; ///< Index type of the orthtree to be traversed
 
   /*!
-    \brief returns the first node to iterate to, given the root of the Orthtree.
+    \brief returns the first node of the traversal.
    */
-  Node first (Node root) const;
+  Node_index first_index() const;
 
   /*!
-    \brief returns the next node to iterate to, given the previous node.
+    \brief returns the next node to be traversed given the previous node `n`.
    */
-  Node next(Node n) const;
+  Node_index next(Node_index n) const;
 };

Orthtree/doc/Orthtree/Doxyfile.in

@@ -5,3 +5,7 @@ PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Quadtrees, Octrees, and Orthtrees"
 EXTRACT_ALL                =  false
 HIDE_UNDOC_MEMBERS         =  true
 HIDE_UNDOC_CLASSES         =  true
+WARN_IF_UNDOCUMENTED       =  false
+
+HTML_EXTRA_FILES           =  ${CGAL_PACKAGE_DOC_DIR}/fig/orthtree.png \
+                              ${CGAL_PACKAGE_DOC_DIR}/fig/orthtree_mesh.png

Orthtree/doc/Orthtree/Orthtree.txt

@@ -10,9 +10,9 @@ namespace CGAL {
 \section Section_Orthtree_Introduction Introduction
 
 Quadtrees are tree data structures in which each node encloses a
-square section of space, and each internal node has exactly 4
+rectangular section of space, and each internal node has exactly 4
 children. Octrees are a similar data structure in 3D in which each
-node encloses a cubic section of space, and each internal node has
+node encloses a rectangular cuboid section of space, and each internal node has
 exactly 8 children.
 
 We call the generalization of such data structure "orthtrees", as
@@ -22,20 +22,32 @@ structures in dimensions 4 and higher.
 
 This package provides a general data structure `Orthtree` along with
 aliases for `Quadtree` and `Octree`. These trees can be constructed
-with custom point ranges and split predicates, and iterated on with
-various traversal methods.
-
-\cgalFigureBegin{Orthtree_fig, orthtree.png}
-Building an %orthtree in 3D (%octree) from a point cloud.
-\cgalFigureEnd
+with custom contents and split predicates, and iterated on with
+various traversal methods. Orthants can be orthotopes and not only hypercubes.
 
+\cgalFigureAnchor{Orthtree_fig  }
+<center>
+<table border=0>
+<tr>
+  <td><img border="0" src="orthtree.png"></td>
+</tr>
+<tr>
+  <td><img border="0" src="orthtree_mesh.png"></td>
+</tr>
+</table>
+</center>
+\cgalFigureCaptionBegin{Orthtree_fig}
+Top: an %orthtree in 3D (%octree) from a point cloud (top);
+Bottom: an %orthtree in 3D (%octree) from the triangle faces of a mesh.
+\cgalFigureCaptionEnd
 
 \section Section_Orthtree_Building Building
 
-An orthtree is created using a set of points. The points are not
-copied: the provided point range is used directly and is rearranged by
+A common purpose for an orthtree is to subdivide a collection of points,
+and the Orthtree package provides a traits class for this purpose.
+The points are not copied: the provided point range is used directly and rearranged by
 the orthtree. Altering the point range after creating the orthtree
-might leave it in an invalid state. The constructor returns a tree
+may leave the tree in an invalid state. The constructor returns a tree
 with a single (root) node that contains all the points.
 
 The method [refine()](@ref CGAL::Orthtree::refine) must be called to
@@ -43,9 +55,11 @@ subdivide space further. This method uses a split predicate which
 takes a node as input and returns `true` if this node should be
 split, `false` otherwise: this enables users to choose on what
 criterion should the orthtree be refined. Predefined predicates are
-provided such as [Maximum_depth](@ref CGAL::Orthtrees::Maximum_depth) or [Maximum_number_of_inliers](@ref CGAL::Orthtrees::Maximum_number_of_inliers).
+provided for the point-set orthtree,
+including [Maximum_depth](@ref CGAL::Orthtrees::Maximum_depth)
+and [Maximum_contained_elements](@ref CGAL::Orthtrees::Maximum_contained_elements).
 
-The simplest way to create an orthtree is using a vector of points.
+The simplest way to create a point-set orthtree is from a vector of points.
 The constructor generally expects a separate point range and map,
 but the point map defaults to `Identity_property_map` if none is provided.
 
@@ -55,29 +69,27 @@ the existing ones do not match users' needs.
 
 \subsection Section_Orthtree_Quadtree Building a Quadtree
 
-The `Orthtree` class may be templated with `Orthtree_traits_2` and thus
-behave as a %quadtree. For convenience, the alias `Quadtree` is provided.
+The `Orthtree` class may be templated with `Orthtree_traits_point<>`
+while specifying a 2d ambient space and thus behave as a %quadtree.
+For convenience, the alias `CGAL::Quadtree` is provided.
 
-The following example shows how to create a %quadtree object from a
-vector of `Point_2` objects and refine it, which means constructing
-the tree's space subdivision itself, using a maximum depth of 10 and a
-maximum number of inliers per node (bucket size) of 5. The refinement
-is stopped as soon as one of the conditions is violated: if a node has
-more inliers than `bucket_size` but is already at `max_depth`, it is
-not split. Similarly, a node that is at a depth smaller than
-`max_depth` but already has fewer inliers than `bucket_size` is not
-split.
+The following example shows the construction of %quadtree from a vector of `Point_2` objects.
+`quadtree.refine(10, 5)` uses the default split predicate, which
+enforces a max-depth and a maximum of elements contained per node ("bucket size").
+Nodes are split if their depth is less than 10, and they contain more than 5 points.
 
 
 \cgalExample{Orthtree/quadtree_build_from_point_vector.cpp}
 
 \subsection Section_Orthtree_Point_Vector Building an Octree
 
-The `Orthtree` class may be templated with `Orthtree_traits_3` and thus
-behave as an %octree. For convenience, the alias `Octree` is provided.
+`Orthtree_traits_point<>` can also be templated with dimension 3 and thus
+behave as an %octree. For convenience, the alias `CGAL::Octree` is provided.
+
+The following example shows how to create an %octree from a vector of `Point_3` objects.
+As with the %quadtree example, we use the default split predicate.
+In this case the maximum depth is 10, and the bucket size is 1.
 
-The following example shows how to create an %octree from a vector of
-`Point_3` objects:
 
 \cgalExample{Orthtree/octree_build_from_point_vector.cpp}
 
@@ -87,12 +99,12 @@ Some data structures such as `Point_set_3` require a non-default point
 map type and object. This example illustrates how to create an octree from a `Point_set_3` loaded from a file.
 It also shows a more explicit way of setting the split predicate when refining the tree.
 
-An octree is constructed from the point set and its map.
-The tree is refined with a maximum depth (deepest node allowed) of 10,
-and a bucket size (maximum number of points contained by a single node) of 20.
-The tree is then written to the standard output.
+An %octree is constructed from the point set and its corresponding map,
+and then written to the standard output.
 
 The split predicate is manually constructed and passed to the refine method.
+In this case, we use a maximum number of contained elements with no corresponding maximum depth.
+This means that nodes will continue to be split until none contain more than 10 points.
 
 \cgalExample{Orthtree/octree_build_from_point_set.cpp}
 
@@ -108,38 +120,47 @@ at depth 7 can hold 14.
 
 \subsection Section_Orthtree_Orthtree_Point_Vector Building an Orthtree
 
-The following example shows how to build an generalized orthtree in dimension 4.
-A `std::vector<Point_d>` is manually filled with points.
-The vector is used as the point set,
-an `Identity_property_map` is automatically set as the orthtree's map type, so a map does not need to be provided.
+An orthtree can also be used with an arbitrary number of dimensions.
+The `Orthtree_traits_point` template can infer the arbitrary dimension count from the d-dimensional kernel.
+
+The following example shows how to build a generalized orthtree in dimension 4.
+As `std::vector<Point_d>` is manually filled with 4-dimensional points.
+The vector is used as the point set, and an `Identity_property_map` is automatically
+set as the orthtree's map type, so a map does not need to be provided.
 
 \cgalExample{Orthtree/orthtree_build.cpp}
 
+\section Section_Orthtree_Properties Properties
+The Orthtree uses a mechanism to attach properties to nodes at run-time which follows \ref sectionSurfaceMesh_properties "Surface mesh properties". Each property is identified by a string and its value type and stored in a consecutive block of memory.
+
 \section Section_Orthtree_Traversal Traversal
 
-\note For simplicity, the rest of the user manual will only use
-octrees, but all the presented features also apply to quadtrees and
-higher dimension orthtrees.
-
 %Traversal is the act of navigating among the nodes of the tree.
-The `Orthtree` and [Node](@ref CGAL::Orthtree::Node) classes provide a
+The `Orthtree` class provides a
 number of different solutions for traversing the tree.
 
 \subsection Section_Orthtree_Manual_Traveral Manual Traversal
 
 Because our orthtree is a form of connected acyclic undirected graph, it is possible to navigate between any two nodes.
-What that means in practice, is that given a node on the tree, it is possible to
+What that means in practice is that given a node on the tree, it is possible to
 access any other node using the right set of operations.
-The `Node` class provides functions that enable the user to access each of its children, as well as its parent (if it exists).
+The `Node_index` type provides a handle on a node, and the `orthtree` class provides methods
+that enable the user to retrieve the indices of each of its children as well as its parent (if they exist).
 
-The following example demonstrates ways of accessing different nodes of a tree, given a reference to one.
+Given any node index `nid`, the `n`th child of that node can be found with `orthtree.child(nid, n)`.
+For an octree, values of `n` from 0-7 provide access to the different children.
+For non-root nodes, it is possible to access parent nodes using the `orthtree.parent()` accessor.
 
-From the root node, children can be accessed using the subscript operator `CGAL::Orthtree::Node::operator[]()`.
-For an octree, values from 0-7 provide access to the different children.
+To access grandchildren, it isn't necessary to use nested `orthtree.child()` calls.
+Instead, the `orthtree.descendant(node, a, b, ...)` convenience method is provided.
+This retrieves the `b`th child of the `a`th child, to any depth.
 
-For non-root nodes, it is possible to access parent nodes using the [parent()](@ref CGAL::Orthtree::Node::parent) accessor.
+In most cases, we want to find the descendants of the root node.
+For this case, there is another convenience method `orthtree.node(a, b, c, ...)`.
+This retrieves the node specified by the descent `a`, `b`, `c`.
+It is equivalent to `orthtree.descendant(orthtree.root(), a, b, c, ...)`.
 
-These accessors and operators can be chained to access any node in the tree in a single line of code, as shown in the following example:
+The following example demonstrates the use of several of these accessors.
 
 \cgalExample{Orthtree/octree_traversal_manual.cpp}
 
@@ -147,9 +168,10 @@ These accessors and operators can be chained to access any node in the tree in a
 
 It is often useful to be able to iterate over the nodes of the tree in a particular order.
 For example, the stream operator `<<` uses a traversal to print out each node.
-A few traversals are provided, among them [Preorder_traversal](@ref CGAL::Orthtrees::Preorder_traversal) and [Postorder_traversal](@ref CGAL::Orthtrees::Postorder_traversal).
-To traverse a tree in preorder is to visit each parent immediately followed by its children,
-whereas in postorder, traversal the children are visited first.
+A few traversals are provided, among them [Preorder_traversal](@ref CGAL::Orthtrees::Preorder_traversal)
+and [Postorder_traversal](@ref CGAL::Orthtrees::Postorder_traversal).
+Traversing a tree in preorder means to visit each parent immediately followed by its children,
+whereas in postorder traversal the children are visited first.
 
 The following example illustrates how to use the provided traversals.
 
@@ -161,13 +183,15 @@ In this case, we print out the nodes of the tree without indentation instead.
 
 \subsection Section_Orthtree_Custom_Traversal Custom Traversal
 
-Users can define their own traversal methods by creating models of the
-`OrthtreeTraversal` concept. The following example shows how to define a
-custom traversal that only traverses the first branch of the octree:
+Users can define their own traversal methods by creating models of the `OrthtreeTraversal` concept.
+The custom traversal must provide a method which returns the starting point of the traversal (`first_index()`)
+and another method which returns the next node in the sequence (`next_index()`).
+`next_index()` returns an empty optional when the end of the traversal is reached.
+The following example shows how to define a custom traversal that only traverses the first branch an octree:
 
 \cgalExample{Orthtree/octree_traversal_custom.cpp}
 
-\subsection Comparison of Traversals
+\subsection Section_Orthtree_Cmp_Trav Comparison of Traversals
 
 Figure \cgalFigureRef{Orthtree_traversal_fig} shows in which order
 nodes are visited depending on the traversal method used.
@@ -184,12 +208,11 @@ Once an orthtree is built, its structure can be used to accelerate different tas
 
 \subsection Section_Orthtree_Nearest_Neighbor Finding the Nearest Neighbor of a Point
 
-The naive way of finding the nearest neighbor of a point requires finding the distance to every other point.
+The naive way of finding the nearest neighbor of a point requires finding the distance to all elements.
 An orthtree can be used to perform the same task in significantly less time.
-For large numbers of points, this can be a large enough difference to outweigh the time spent building the tree.
+For large numbers of elements, this can be a large enough difference to outweigh the time spent building the tree.
 
-Note that a kd-tree is expected to outperform the orthtree for this task,
-it should be preferred unless features specific to the orthtree are needed.
+Note that a kd-tree is expected to outperform the orthtree for this task on points, it should be preferred unless features specific to the orthtree are needed.
 
 The following example illustrates how to use an octree to accelerate the search for points close to a location.
 
@@ -200,9 +223,14 @@ Results are put in a vector, and then printed.
 
 \cgalExample{Orthtree/octree_find_nearest_neighbor.cpp}
 
+Not all octrees are compatible with nearest neighbor functionality,
+as the idea of a nearest neighbor may not make sense for some tree contents.
+For the nearest neighbor methods to work, the traits class must implement the
+`CollectionPartitioningOrthtreeTraits` concept.
+
 \subsection Section_Orthtree_Grade Grading
 
-An orthtree is graded if the difference of depth between two adjacent
+An orthtree is considered "graded" if the difference of depth between two adjacent
 leaves is at most 1 for every pair of leaves.
 
 \cgalFigureBegin{Orthtree_quadree_graded_fig, quadtree_graded.png}
@@ -256,15 +284,90 @@ purposes.
 
 For nontrivial point counts, the naive approach's calculation time dwarfs that of either the %orthtree or kd-tree.
 
+\section Section_Orthtree_Migration Migrating Code Written Before Release 6.0
+
+The orthtree package changed to allow for custom data stored per node in the orthtree. To migrate existing code written before \cgal 6.0 `Orthtree_traits_point` can be used for a point-based orthtrees. The aliases `CGAL::Quadtree` and `CGAL::Octree` have been extended by a boolean template parameter to allow for non-cubic cells, which is the default. The data is passed via the traits class and no longer directly to the orthtree.
+
+Former code to declare and define an Octree with cubic cells was as follows:
+\code{.cpp}
+typedef CGAL::Point_set_3<Kernel::Point_3> Point_set;
+typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
+...
+Octree octree(points, points.point_map());
+\endcode
+
+\cgal 6.0 code with identical behavior is now:
+\code{.cpp}
+typedef CGAL::Point_set_3<Kernel::Point_3> Point_set;
+typedef CGAL::Octree<Kernel, Point_set, Point_map, false> Octree;
+...
+Octree octree(points, points.point_map());
+\endcode
+
+The node class does not exist anymore and has been replaced by the lightweight type `Node_index`. All information formerly contained in the node class is now accessible via the `Orthtree` interface.
+
+Former node access was as follows:
+\code{.cpp}
+Orthtree::Node root = orthtree.root();
+Orthtree::Node child = root[0];
+bool is_leaf = child.is_leaf();
+
+for (Orthtree::Node::const_iterator it = child.begin(); it != child.end(); it++) {
+  const Orthtree::Point &p = *it;
+  ...
+}
+\endcode
+
+\cgal 6.0 node access is now:
+\code{.cpp}
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+Point_set points;
+...
+Orthtree::Node_index root = orthtree.root();
+Orthtree::Node_index child = orthtree.child(root, 0);
+bool is_leaf = orthtree.is_leaf(child);
+
+for (Octree::Traits::Node_data_element& e : octree.data(child)) {
+  // Using a pointmap is necessary when using a Point_set_3.
+  Point& p = points.point(e);
+
+  // If the orthtree is used with a std::vector<Point>, Node_data_element is identical to Point.
+  // Point& p = e;
+  ...
+}
+\endcode
+
+The nearest neighbor search behaves as before, however, the output iterator will be required to store `CollectionPartitioningOrthtreeTraits::Node_data_element`. This may be the actual point or, e.g., in case a `Point_set_3` is used, an index which requires the use of a point map to retrieve the point.
+
+The provided traversals, i.e., `CGAL::Orthtrees::Leaves_traversal`, `CGAL::Orthtrees::Preorder_traversal`, `CGAL::Orthtrees::Postorder_traversal`, require the orthtree as template parameter now.
+
+Former traversal use was as follows:
+\code{.cpp}
+for (Orthtree::Node node : orthtree.traverse<CGAL::Orthtrees::Leaves_traversal>())
+  ...
+\endcode
+
+\cgal 6.0 traversal use is now:
+\code{.cpp}
+for (Orthtree::Node_index node : orthtree.traverse<CGAL::Orthtrees::Leaves_traversal<Orthtree>>())
+  ...
+\endcode
+
 \section Section_Orthtree_History History
 
 A prototype code was implemented by Pierre Alliez and improved by Tong
 Zhao and Cédric Portaneri. From this prototype code, the package was
-developed by Jackson Campolatarro as part of the Google Summer of Code
-2020. Simon Giraudot, supervisor of the GSoC internship, completed and
+developed by Jackson Campolattaro as part of the Google Summer of Code 2020.
+Simon Giraudot, supervisor of the GSoC internship, completed and
 finalized the package for integration in CGAL 5.3. Pierre Alliez
 provided kind help and advice all the way through.
+Starting with CGAL 6.0 an API improvement of the Orthtree package was released.
+Most notably, the orthtree nodes do not need to store anything. Data to be stored
+by the node are managed through a mechanism of dynamic properties.
+This improvement was done by Jackson Campolattaro thanks to a funding provided by
+INRIA, together with GeometryFactory.
 
 */
 
-}
+}

Orthtree/doc/Orthtree/PackageDescription.txt

@@ -1,20 +1,20 @@
-/// \defgroup PkgOrthtreeRef Quadtree\, Octree and Orthtree Reference
-
-/// \defgroup PkgOrthtreeClasses Classes
-/// \ingroup PkgOrthtreeRef
-
-/// \defgroup PkgOrthtreeTraits Traits
-/// \ingroup PkgOrthtreeClasses
-
-/// \defgroup PkgOrthtreeSplitPredicates Split Predicates
-/// \ingroup PkgOrthtreeClasses
-
-/// \defgroup PkgOrthtreeTraversal Traversal
-/// \ingroup PkgOrthtreeClasses
+/// \defgroup PkgOrthtreeRef Quadtrees, Octrees and Orthtrees Reference
+Quadtree, Octree and Orthtree Reference
 
 /// \defgroup PkgOrthtreeConcepts Concepts
 /// \ingroup PkgOrthtreeRef
 
+/// \defgroup PkgOrthtreeTraits Traits
+/// \ingroup PkgOrthtreeRef
+
+/// \defgroup PkgOrthtreeSplitPredicates Split Predicates
+/// \ingroup PkgOrthtreeRef
+
+/// \defgroup PkgOrthtreeTraversal Traversal
+/// \ingroup PkgOrthtreeRef
+
+
+
 /*!
 \addtogroup PkgOrthtreeRef
 
@@ -39,22 +39,25 @@
 
 \cgalCRPSection{Concepts}
 - `OrthtreeTraits`
+- `OrthtreeTraitsWithData`
 - `OrthtreeTraversal`
+- `CollectionPartitioningOrthtreeTraits`
 
 \cgalCRPSection{Classes}
 - `CGAL::Quadtree<GeomTraits, PointRange, PointMap>`
 - `CGAL::Octree<GeomTraits, PointRange, PointMap>`
-- `CGAL::Orthtree<Traits, PointRange, PointMap>`
+- `CGAL::Orthtree<GeomTraits, PointRange, PointMap>`
 
 \cgalCRPSection{Traits}
-- `CGAL::Orthtree_traits_2<GeomTraits>`
-- `CGAL::Orthtree_traits_3<GeomTraits>`
-- `CGAL::Orthtree_traits_d<GeomTraits, Dimension>`
+- `CGAL::Orthtree_traits<GeomTraits, dimension>`
+- `CGAL::Orthtree_traits_point<GeomTraits, PointRange, PointMap, dimension>`
+- `CGAL::Orthtree_traits_base<GeomTraits>`
+- `CGAL::Orthtree_traits_face_graph<TriangleMesh, VertexPointMap>`
 
 \cgalCRPSection{Split Predicates}
-- `CGAL::Orthtrees::Maximum_number_of_inliers`
+- `CGAL::Orthtrees::Maximum_contained_elements`
 - `CGAL::Orthtrees::Maximum_depth`
-- `CGAL::Orthtrees::Maximum_depth_and_maximum_number_of_inliers`
+- `CGAL::Orthtrees::Maximum_depth_and_maximum_contained_elements`
 
 \cgalCRPSection{%Traversal}
 - `CGAL::Orthtrees::Preorder_traversal`

Orthtree/doc/Orthtree/dependencies

@@ -8,3 +8,4 @@ Property_map
 STL_Extension
 Spatial_searching
 Stream_support
+Surface_mesh

Orthtree/doc/Orthtree/examples.txt

@@ -3,6 +3,7 @@
 \example Orthtree/octree_build_from_point_vector.cpp
 \example Orthtree/octree_build_with_custom_split.cpp
 \example Orthtree/octree_find_nearest_neighbor.cpp
+\example Orthtree/octree_surface_mesh.cpp
 \example Orthtree/octree_traversal_manual.cpp
 \example Orthtree/octree_traversal_preorder.cpp
 \example Orthtree/octree_traversal_custom.cpp

Orthtree/examples/Orthtree/CMakeLists.txt

@@ -15,6 +15,8 @@ create_single_source_cgal_program("octree_traversal_manual.cpp")
 create_single_source_cgal_program("octree_traversal_preorder.cpp")
 create_single_source_cgal_program("octree_grade.cpp")
 create_single_source_cgal_program("quadtree_build_from_point_vector.cpp")
+create_single_source_cgal_program("octree_surface_mesh.cpp")
+create_single_source_cgal_program("quadtree_build_manually.cpp")
 
 find_package(Eigen3 3.1.91 QUIET) #(requires 3.1.91 or greater)
 include(CGAL_Eigen_support)

Orthtree/examples/Orthtree/octree_build_from_point_set.cpp

@@ -7,12 +7,11 @@
 #include <CGAL/Point_set_3/IO.h>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef Point_set::Point_map Point_map;
-
-typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Point_map = Point_set::Point_map;
+using Octree = CGAL::Octree<Kernel, Point_set, Point_map>;
 
 int main(int argc, char **argv) {
 
@@ -33,10 +32,10 @@ int main(int argc, char **argv) {
   Octree octree(points, points.point_map());
 
   // Build the octree with a small bucket size, using a more verbose method
-  octree.refine(CGAL::Orthtrees::Maximum_depth_and_maximum_number_of_inliers(5, 10));
+  octree.refine(CGAL::Orthtrees::Maximum_contained_elements(10));
 
   // Print out the tree
-  std::cout << octree;
+  std::cout << octree << std::endl;
 
   return EXIT_SUCCESS;
 }

Orthtree/examples/Orthtree/octree_build_from_point_vector.cpp

@@ -4,11 +4,10 @@
 #include <CGAL/Octree.h>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef std::list<Point> Point_vector;
-
-typedef CGAL::Octree<Kernel, Point_vector> Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_vector = std::vector<Point>;
+using Octree = CGAL::Octree<Kernel, Point_vector>;
 
 int main() {
 
@@ -27,7 +26,7 @@ int main() {
   Octree octree(points);
 
   // Build the octree
-  octree.refine(10, 20);
+  octree.refine(10, 1);
 
   // Print out the tree
   std::cout << octree;

Orthtree/examples/Orthtree/octree_build_with_custom_split.cpp

@@ -7,25 +7,26 @@
 #include <CGAL/Point_set_3/IO.h>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef Kernel::FT FT;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef Point_set::Point_map Point_map;
-
-typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using FT = Kernel::FT;
+using Point_set = CGAL::Point_set_3<Point>;
+using Point_map = Point_set::Point_map;
+using Octree = CGAL::Octree<Kernel, Point_set, Point_map>;
 
 // Split Predicate
-// The predicate is a functor which returns a boolean value, whether a node needs to be split or not
+// The predicate is a functor which returns a Boolean value, whether a node needs to be split or not
 struct Split_by_ratio {
 
   std::size_t ratio;
 
-  Split_by_ratio(std::size_t ratio) : ratio(ratio) {}
+  explicit Split_by_ratio(std::size_t ratio) : ratio(ratio) {}
 
-  template<class Node>
-  bool operator()(const Node &n) const {
-    return n.size() > (ratio * n.depth());
+  template<typename Node_index, typename Tree>
+  bool operator()(Node_index i, const Tree &tree) const {
+    std::size_t num_points = tree.data(i).size();
+    std::size_t depth = tree.depth(i);
+    return num_points > (ratio * depth);
   }
 };
 

Orthtree/examples/Orthtree/octree_find_nearest_neighbor.cpp

@@ -9,14 +9,13 @@
 #include <boost/iterator/function_output_iterator.hpp>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef Point_set::Point_map Point_map;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Point_map = Point_set::Point_map;
+using Octree = CGAL::Octree<Kernel, Point_set, Point_map>;
 
-typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
-
-int main(int argc, char **argv) {
+int main(int argc, char** argv) {
 
   // Point set will be used to hold our points
   Point_set points;
@@ -39,21 +38,42 @@ int main(int argc, char **argv) {
 
   // Find the nearest points to a few locations
   std::vector<Point> points_to_find = {
-          {0, 0, 0},
-          {1, 1, 1},
-          {-1, -1, -1},
-          {-0.46026, -0.25353, 0.32051},
-          {-0.460261, -0.253533, 0.320513}
+    {0,         0,         0},
+    {1,         1,         1},
+    {-1,        -1,        -1},
+    {-0.46026,  -0.25353,  0.32051},
+    {-0.460261, -0.253533, 0.320513}
   };
+
   for (const Point& p : points_to_find)
-    octree.nearest_neighbors
-      (p, 1, // k=1 to find the single closest point
-       boost::make_function_output_iterator
-       ([&](const Point& nearest)
+    octree.nearest_k_neighbors
+    (p, 1, // k=1 to find the single closest point
+      boost::make_function_output_iterator
+      ([&](const Point_set::Index& nearest)
         {
           std::cout << "the nearest point to (" << p <<
-            ") is (" << nearest << ")" << std::endl;
+            ") is (" << points.point(nearest) << ")" << std::endl;
         }));
 
+  typename Octree::Sphere s(points_to_find[0], 0.0375);
+  std::cout << std::endl << "Closest points within the sphere around " << s.center() << " with squared radius of " << s.squared_radius() << ":" << std::endl;
+  octree.neighbors_within_radius
+  (s,
+    boost::make_function_output_iterator
+    ([&](const Point_set::Index& nearest)
+      {
+        std::cout << points.point(nearest) << "    dist: " << (Point(0, 0, 0) - points.point(nearest)).squared_length() << std::endl;
+      }));
+
+  std::size_t k = 3;
+  std::cout << std::endl << "The up to " << k << " closest points to(" << s.center() << ") within a squared radius of " << s.squared_radius() << " are:" << std::endl;
+  octree.nearest_k_neighbors_within_radius
+  (s, k,
+    boost::make_function_output_iterator
+    ([&](const Point_set::Index& nearest)
+      {
+        std::cout << "(" << points.point(nearest) << "    dist: " << (Point(0, 0, 0) - points.point(nearest)).squared_length() << std::endl;
+      }));
+
   return EXIT_SUCCESS;
 }

Orthtree/examples/Orthtree/octree_grade.cpp

@@ -4,10 +4,10 @@
 #include <CGAL/Octree.h>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef std::vector<Point> Point_vector;
-typedef CGAL::Octree<Kernel, Point_vector> Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_vector = std::vector<Point>;
+using Octree = CGAL::Octree<Kernel, Point_vector>;
 
 int main() {
 

Orthtree/examples/Orthtree/octree_surface_mesh.cpp

@@ -0,0 +1,68 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Octree.h>
+#include <CGAL/Orthtree_traits_face_graph.h>
+
+#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Mesh = CGAL::Surface_mesh<K::Point_3>;
+
+using OTraits = CGAL::Orthtree_traits_face_graph<Mesh, Mesh::Property_map<Mesh::Vertex_index, K::Point_3>>;
+using Octree = CGAL::Orthtree<OTraits>;
+
+void dump_as_polylines(const Octree& ot)
+{
+    std::ofstream out("octree.polylines.txt");
+    for (Octree::Node_index node : ot.traverse(CGAL::Orthtrees::Leaves_traversal<Octree>(ot)))
+    {
+      if (!ot.is_leaf(node))
+        continue;
+      auto bb = ot.bbox(node);
+      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin()
+          << "  " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin() << "\n";
+      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin()
+          << "  " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin() << "\n";
+      out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin()
+          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin() << "\n";
+      out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin()
+          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin() << "\n";
+//
+      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin()
+          << "  " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax() << "\n";
+      out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin()
+          << "  " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax() << "\n";
+      out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin()
+          << "  " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax() << "\n";
+      out << "2 " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin()
+          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n";
+//
+      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax()
+          << "  " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax() << "\n";
+      out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax()
+          << "  " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax() << "\n";
+      out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax()
+          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n";
+      out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax()
+          << "  " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n";
+    }
+}
+
+int main(int argc, char** argv)
+{
+  const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off");
+
+  Mesh mesh;
+  if(!CGAL::IO::read_polygon_mesh(filename, mesh))
+  {
+    std::cerr << "Error: cannot read file" << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  Octree tree(mesh, mesh.points());
+  OTraits::Split_predicate_node_min_extent sp(0.01);
+  tree.refine(sp);
+
+  dump_as_polylines(tree);
+}

Orthtree/examples/Orthtree/octree_traversal_custom.cpp

@@ -3,34 +3,42 @@
 
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Octree.h>
+#include <CGAL/Orthtree/IO.h>
 #include <CGAL/Point_set_3.h>
 #include <CGAL/Point_set_3/IO.h>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef Point_set::Point_map Point_map;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Point_map = Point_set::Point_map;
+using Octree = CGAL::Octree<Kernel, Point_set, Point_map>;
 
-typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
-typedef Octree::Node Node;
+template <typename Tree>
+struct First_branch_traversal {
 
-struct First_branch_traversal
-{
-  Node first (Node root) const
-  {
-    return root;
+  using Node_index = typename Tree::Node_index;
+
+  const Tree& m_orthtree;
+
+  explicit First_branch_traversal(const Tree& orthtree) : m_orthtree(orthtree) {}
+
+  Node_index first_index() const {
+    return m_orthtree.root();
   }
 
-  Node next (Node n) const
-  {
-    if (n.is_leaf())
-      return Node();
-    return n[0];
+  std::optional<Node_index> next_index(Node_index n) const {
+
+    // Stop when we reach the base of the tree
+    if (m_orthtree.is_leaf(n))
+      return {};
+
+    // Always descend the first child
+    return m_orthtree.child(n, 0);
   }
 };
 
-int main(int argc, char **argv) {
+int main(int argc, char** argv) {
 
   // Point set will be used to hold our points
   Point_set points;
@@ -52,8 +60,9 @@ int main(int argc, char **argv) {
   octree.refine();
 
   // Print out the first branch using custom traversal
-  for (auto &node : octree.traverse<First_branch_traversal>())
-    std::cout << node << std::endl;
+  for (auto node: octree.traverse<First_branch_traversal<Octree>>()) {
+    std::cout << octree.to_string(node) << std::endl;
+  }
 
   return EXIT_SUCCESS;
 }

Orthtree/examples/Orthtree/octree_traversal_manual.cpp

@@ -7,14 +7,13 @@
 #include <CGAL/Point_set_3/IO.h>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef Point_set::Point_map Point_map;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Point_map = Point_set::Point_map;
+using Octree = CGAL::Octree<Kernel, Point_set, Point_map>;
 
-typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
-
-int main(int argc, char **argv) {
+int main(int argc, char** argv) {
 
   // Point set will be used to hold our points
   Point_set points;
@@ -39,29 +38,30 @@ int main(int argc, char **argv) {
   std::cout << "Navigation relative to the root node" << std::endl;
   std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
   std::cout << "the root node: " << std::endl;
-  std::cout << octree.root() << std::endl;
+  std::cout << octree.to_string(octree.root()) << std::endl;
   std::cout << "the first child of the root node: " << std::endl;
-  std::cout << octree.root()[0] << std::endl;
+  std::cout << octree.to_string(octree.child(octree.root(), 0)) << std::endl;
   std::cout << "the fifth child: " << std::endl;
-  std::cout << octree.root()[4] << std::endl;
-  std::cout << "the fifth child, accessed without the root keyword: " << std::endl;
-  std::cout << octree[4] << std::endl;
+  std::cout << octree.to_string(octree.child(octree.root(), 4)) << std::endl;
+  std::cout << "the fifth child, accessed without going through root: " << std::endl;
+  std::cout << octree.to_string(octree.node(4)) << std::endl;
   std::cout << "the second child of the fourth child: " << std::endl;
-  std::cout << octree.root()[4][1] << std::endl;
-  std::cout << "the second child of the fourth child, accessed without the root keyword: " << std::endl;
-  std::cout << octree[4][1] << std::endl;
+  std::cout << octree.to_string(octree.child(octree.child(octree.root(), 4), 1)) << std::endl;
+  std::cout << "the second child of the fourth child, accessed without going through root: " << std::endl;
+  std::cout << octree.to_string(octree.node(4, 1)) << std::endl;
   std::cout << std::endl;
 
   // Retrieve one of the deeper children
-  Octree::Node cur = octree[3][2];
+  Octree::Node_index cur = octree.node(4, 3);
   std::cout << "Navigation relative to a child node" << std::endl;
   std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
   std::cout << "the third child of the fourth child: " << std::endl;
-  std::cout << cur << std::endl;
+  std::cout << octree.to_string(cur) << std::endl;
   std::cout << "the third child: " << std::endl;
-  std::cout << cur.parent() << std::endl;
+  std::cout << octree.to_string(octree.parent(cur)) << std::endl;
   std::cout << "the next sibling of the third child of the fourth child: " << std::endl;
-  std::cout << cur.parent()[cur.local_coordinates().to_ulong() + 1] << std::endl;
+  std::cout << octree.to_string(octree.child(octree.parent(cur), octree.local_coordinates(cur).to_ulong() + 1))
+            << std::endl;
 
   return EXIT_SUCCESS;
 }

Orthtree/examples/Orthtree/octree_traversal_preorder.cpp

@@ -7,13 +7,12 @@
 #include <CGAL/Point_set_3/IO.h>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef Point_set::Point_map Point_map;
-
-typedef CGAL::Octree<Kernel, Point_set, Point_map> Octree;
-typedef CGAL::Orthtrees::Preorder_traversal Preorder_traversal;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Point_map = Point_set::Point_map;
+using Octree = CGAL::Octree<Kernel, Point_set, Point_map>;
+using Preorder_traversal = CGAL::Orthtrees::Preorder_traversal<Octree>;
 
 int main(int argc, char **argv) {
 
@@ -37,9 +36,8 @@ int main(int argc, char **argv) {
   octree.refine();
 
   // Print out the octree using preorder traversal
-  for (Octree::Node node : octree.traverse<Preorder_traversal>()) {
-
-    std::cout << node << std::endl;
+  for (auto node : octree.traverse<Preorder_traversal>()) {
+    std::cout << octree.to_string(node) << std::endl;
   }
 
   return EXIT_SUCCESS;

Orthtree/examples/Orthtree/orthtree_build.cpp

@@ -1,27 +1,27 @@
 #include <iostream>
 
 #include <CGAL/Epick_d.h>
+#include <CGAL/Cartesian_d.h>
 #include <CGAL/Orthtree.h>
-#include <CGAL/Orthtree_traits_d.h>
+#include <CGAL/Orthtree_traits_point.h>
 #include <CGAL/Random.h>
 
 // Type Declarations
-typedef CGAL::Dimension_tag<4> Dimension;
-typedef CGAL::Epick_d<Dimension> Kernel;
-typedef Kernel::Point_d Point_d;
-typedef std::vector<Point_d> Point_vector;
-
-typedef CGAL::Orthtree_traits_d<Kernel, Dimension> Traits;
-typedef CGAL::Orthtree<Traits, Point_vector> Orthtree;
+const int dimension = 4;
+using Kernel = CGAL::Epick_d<CGAL::Dimension_tag<dimension> >;
+using Point_d = Kernel::Point_d;
+using Point_vector = std::vector<Point_d>;
+using Traits = CGAL::Orthtree_traits_point<Kernel, Point_vector>;
+using Orthtree = CGAL::Orthtree<Traits>;
 
 int main()
 {
   CGAL::Random r;
 
   Point_vector points_dd;
-  for (std::size_t i = 0; i < 5; ++ i)
+  for (std::size_t i = 0; i < 20; ++ i)
   {
-    std::array<double, Dimension::value> init;
+    std::array<double, dimension> init{};
     for (double& v : init)
       v = r.get_double(-1., 1.);
     points_dd.emplace_back (init.begin(), init.end());
@@ -30,5 +30,7 @@ int main()
   Orthtree orthtree(points_dd);
   orthtree.refine(10, 5);
 
+  std::cout << orthtree << std::endl;
+
   return EXIT_SUCCESS;
 }

Orthtree/examples/Orthtree/quadtree_build_from_point_vector.cpp

@@ -5,11 +5,10 @@
 #include <CGAL/Random.h>
 
 // Type Declarations
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_2 Point_2;
-typedef std::vector<Point_2> Point_vector;
-
-typedef CGAL::Quadtree<Kernel, Point_vector> Quadtree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point_2 = Kernel::Point_2;
+using Point_vector = std::vector<Point_2>;
+using Quadtree = CGAL::Quadtree<Kernel, Point_vector>;
 
 int main()
 {

Orthtree/examples/Orthtree/quadtree_build_manually.cpp

@@ -0,0 +1,59 @@
+#include <iostream>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Dimension.h>
+
+#include <CGAL/Orthtree.h>
+#include <CGAL/Orthtree_traits_base.h>
+
+using Kernel = CGAL::Simple_cartesian<double>;
+
+namespace CGAL {
+
+struct empty_type {
+};
+
+template <typename K, int dimension>
+struct Orthtree_traits_empty : public Orthtree_traits_base<K, dimension> {
+
+  using Self = Orthtree_traits_empty<K, dimension>;
+  using Tree = Orthtree<Self>;
+
+  using Node_data = std::array<empty_type, 0>;
+
+  Orthtree_traits_empty(typename Self::Bbox_d bbox) : m_bbox(bbox) {};
+
+  auto construct_root_node_bbox_object() const {
+    return [&]() -> typename Self::Bbox_d { return m_bbox; };
+  }
+
+  auto construct_root_node_contents_object() const { return [&]() -> Node_data { return {}; }; }
+
+  auto distribute_node_contents_object() {
+    return [&](typename Tree::Node_index /* n */, Tree& /* tree */, const typename Self::Point_d& /* center */) -> void {};
+  }
+
+private:
+
+  typename Self::Bbox_d m_bbox;
+
+};
+}
+
+using EmptyQuadtree = CGAL::Orthtree<CGAL::Orthtree_traits_empty<Kernel, 2 >>;
+
+int main() {
+
+  // Build an empty quadtree which covers the domain (-1, -1) to (1, 1)
+  EmptyQuadtree quadtree{{{-1, -1, 1, 1}}};
+
+  // Split several nodes of the tree
+  quadtree.split(quadtree.root());
+  quadtree.split(quadtree.node(0));
+  quadtree.split(quadtree.node(3));
+  quadtree.split(quadtree.node(3, 0));
+
+  std::cout << quadtree.bbox(quadtree.root()) << std::endl;
+  std::cout << quadtree << std::endl;
+  return EXIT_SUCCESS;
+}

Orthtree/include/CGAL/Octree.h

@@ -15,34 +15,27 @@
 #include <CGAL/license/Orthtree.h>
 
 #include <CGAL/Orthtree.h>
-#include <CGAL/Orthtree_traits_3.h>
+#include <CGAL/Orthtree_traits_point.h>
 
 namespace CGAL {
 
 /*!
-  \ingroup PkgOrthtreeClasses
+  \ingroup PkgOrthtreeRef
 
-  \brief Alias that specializes the `Orthtree` class to a 3D octree.
+  \brief Alias that specializes the `Orthtree` class to a 3D octree storing 3D points.
 
-  These two types are exactly equivalent:
-  - `Octree<GeomTraits, PointRange, PointMap>`
-  - `Orthtree<Orthtree_traits_3<GeomTraits>, PointRange, PointMap>`.
-
-  \warning This is a not a real class but an alias, please refer to
-  the documentation of `Orthtree`.
-
-  \tparam GeomTraits must be a model of `Kernel`
-  \tparam PointRange_ must be a model of range whose value type is the key type of `PointMap`
-  \tparam PointMap must be a model of `ReadablePropertyMap` whose value type is `GeomTraits::Point_3`
+  \tparam GeomTraits a model of `Kernel`
+  \tparam PointRange a model of `Range` whose value type is the key type of `PointMap`
+  \tparam PointMap a model of `ReadablePropertyMap` whose value type is `GeomTraits::Point_3`
+  \tparam cubic_nodes Boolean to enforce cubic nodes
  */
-template <typename GeomTraits, typename PointRange,
-          typename PointMap = Identity_property_map
-         <typename std::iterator_traits<typename PointRange::iterator>::value_type> >
-#ifdef DOXYGEN_RUNNING
-class Octree;
-#else
-using Octree = Orthtree<Orthtree_traits_3<GeomTraits>, PointRange, PointMap>;
-#endif
+template <
+  typename GeomTraits,
+  typename PointRange,
+  typename PointMap = Identity_property_map<typename std::iterator_traits<typename PointRange::iterator>::value_type>,
+  bool cubic_nodes = false
+>
+using Octree = Orthtree<Orthtree_traits_point<GeomTraits, PointRange, PointMap, cubic_nodes, 3>>;
 
 } // namespace CGAL
 

Orthtree/include/CGAL/Orthtree/Cartesian_ranges.h

@@ -29,8 +29,8 @@ namespace internal
 template <typename Traits>
 struct Cartesian_ranges
 {
-  typedef typename Traits::Point_d Point;
-  typedef typename Traits::Cartesian_const_iterator_d Cartesian_const_iterator;
+  using Point = typename Traits::Point_d;
+  using Cartesian_const_iterator = typename Traits::Cartesian_const_iterator_d;
 
   using Range_single = CGAL::Iterator_range<Cartesian_const_iterator>;
 

Orthtree/include/CGAL/Orthtree/IO.h

@@ -22,46 +22,39 @@ namespace CGAL
 namespace internal
 {
 
-template<typename Node>
-std::ostream& print_orthtree_node(std::ostream& os, const Node& node)
+template<typename Node_index, typename Tree>
+std::ostream& print_orthtree_node(std::ostream& os, const Node_index& node, const Tree &tree)
 {
-  // Show the depth of the node
-//  for (int i = 0; i < node.depth(); ++i)
-//    os << ". ";
 
   // Wrap information in brackets
   os << "{ ";
 
   // Index identifies which child this is
   os << "(";
-  for (std::size_t i = 0; i < node.local_coordinates().size(); ++ i)
-    os << node.local_coordinates()[i];
+  for (std::size_t i = 0; i < tree.local_coordinates(node).size(); ++ i)
+    os << tree.local_coordinates(node)[i];
   os << ") ";
 
   // Location
   os << "( ";
-  for (const auto& b : node.global_coordinates())
+  for (const auto& b : tree.global_coordinates(node))
     os << b << " ";
   os << ") ";
 
   // Depth
   os << "("
-     << +node.depth() // The + forces printing as an int instead of a char
+     << +tree.depth(node) // The + forces printing as an int instead of a char
      << ") ";
 
   os << "("
-     << node.size()
+     << tree.data(node).size()
      << ") ";
 
-//  // If a node has points, indicate how many
-//  if (!node.is_empty())
-//    os << "[" << node.num_points() << " points] ";
-
   // If a node is a leaf, mark it
-  os << (node.is_leaf() ? "[leaf] " : "");
+  os << (tree.is_leaf(node) ? "[leaf] " : "");
 
   // If a node is root, mark it
-  os << (node.is_root() ? "[root] " : "");
+  os << (tree.is_root(node) ? "[root] " : "");
 
   // Wrap information in brackets
   os << "}";

Orthtree/include/CGAL/Orthtree/Node.h

@@ -1,619 +0,0 @@
-// Copyright (c) 2007-2020  INRIA (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)     : Jackson Campolattaro, Cédric Portaneri, Tong Zhao
-
-#ifndef CGAL_ORTHTREE_NODE_H
-#define CGAL_ORTHTREE_NODE_H
-
-#include <CGAL/license/Orthtree.h>
-
-#include <CGAL/Orthtree.h>
-
-#include <boost/range/iterator_range.hpp>
-
-#include <array>
-#include <memory>
-#include <bitset>
-#include <cassert>
-#include <iostream>
-
-namespace CGAL {
-
-/// \cond SKIP_IN_MANUAL
-namespace Orthtrees
-{
-
-// Non-documented, or testing purpose only
-struct Node_access
-{
-  template <typename Node, typename LC>
-  static Node create_node (Node parent, LC local_coordinates)
-  {
-    return Node(parent, local_coordinates);
-  }
-
-  template <typename Node>
-  static typename Node::Point_range& points(Node node) { return node.points(); }
-
-  template <typename Node>
-  static void split(Node node) { return node.split(); }
-
-  template <typename Node>
-  static void free(Node node)
-  {
-    typedef Dimension_tag<(2 << (Node::Dimension::value - 1))> Degree;
-    std::queue<Node> nodes;
-    nodes.push(node);
-    while (!nodes.empty())
-    {
-      Node node = nodes.front();
-      nodes.pop();
-      if (!node.is_leaf()){
-        for (std::size_t i = 0; i < Degree::value; ++ i){
-          nodes.push (node[i]);
-        }
-      }
-      node.free();
-    }
-  }
-
-};
-
-} // namespace Orthtrees
-/// \endcond
-
-/*!
-
-  \brief represents a single node of the tree. Alternatively referred
-  to as a cell, orthant, or sub-tree.
-
-  A `Node` is a lightweight object and thus generally passed by
-  copy. It is also a model of `ConstRange` with value type `Traits::Point_d`.
-
-  \cgalModels{ConstRange}
- */
-template<typename Traits, typename PointRange, typename PointMap>
-class Orthtree<Traits, PointRange, PointMap>::Node
-{
-
-public:
-
-  /// \name Types
-  /// @{
-
-  typedef Orthtree<Traits, PointRange, PointMap> Enclosing; ///< Orthtree type (enclosing class).
-  typedef typename Enclosing::Dimension Dimension; ///< Dimension type.
-  typedef typename Enclosing::Degree Degree; ///< Degree type.
-
-  /*!
-    \brief Self typedef for convenience.
-   */
-  typedef typename Orthtree<Traits, PointRange, PointMap>::Node Self;
-
-
-  /// \cond SKIP_IN_MANUAL
-  /*!
-   * \brief Array for containing the child nodes of this node.
-   */
-  typedef std::array<Self, Degree::value> Children;
-  /// \endcond
-
-  /*!
-    \brief Set of bits representing this node's relationship to its parent.
-
-    Equivalent to an array of Booleans, where index[0] is whether `x`
-    is greater, index[1] is whether `y` is greater, index[2] is whether
-    `z` is greater, and so on for higher dimensions if needed.
-    Used to represent a node's relationship to the center of its parent.
-   */
-  typedef std::bitset<Dimension::value> Local_coordinates;
-
-  /*!
-    \brief Coordinates representing this node's relationship
-    with the rest of the tree.
-
-    Each value `(x, y, z, ...)` of global coordinates is calculated by doubling
-    the parent's global coordinates and adding the local coordinates.
-   */
-  typedef std::array<std::uint32_t, Dimension::value> Global_coordinates;
-
-
-  typedef typename PointRange::const_iterator const_iterator; ///< constant iterator type.
-
-  /*!
-    \brief Adjacency directions.
-   */
-  typedef typename Traits::Adjacency Adjacency;
-
-  /// @}
-
-private:
-
-  /// \cond SKIP_IN_MANUAL
-  /*!
-    \brief point range type.
-   */
-  typedef typename PointRange::iterator iterator; ///< constant iterator type.
-  typedef boost::iterator_range<iterator> Point_range;
-  /// \endcond
-
-  // make Node trivially copiabled
-  struct Data
-  {
-    Point_range points;
-    Self parent;
-    std::uint8_t depth;
-    Global_coordinates global_coordinates;
-    std::unique_ptr<Children> children;
-
-    Data (Self parent)
-      : parent (parent), depth (0) { }
-  };
-
-  Data* m_data;
-
-
-  /// \cond SKIP_IN_MANUAL
-
-  // Only the Orthtree class has access to the non-default
-  // constructor, mutators, etc.
-  friend Enclosing;
-
-  // Hidden class to access methods for testing purposes
-  friend Orthtrees::Node_access;
-
-  /*!
-   * \brief Access to the content held by this node
-   * \return a reference to the collection of point indices
-   */
-  Point_range &points() { return m_data->points; }
-  const Point_range &points() const { return m_data->points; }
-
-  /// \name Construction
-  /// @{
-
-  /*!
-    \brief creates a new node, optionally as the child of a parent
-
-    If no parent is provided, the node created is assumed to be the
-    root of a tree.  This means that `parent.is_null()` returns
-    `true`, and the depth is zero.  If a parent is provided, the node
-    becomes the child of that parent.  In that case, an index should
-    be passed, telling this node its relationship to its parent.
-    Depth and global coordinates are automatically determined in the
-    constructor, and should generally be considered immutable after
-    construction.
-
-    \param parent the node containing this one
-    \param index this node's relationship to its parent
-  */
-  explicit Node(Self parent, Local_coordinates local_coordinates)
-    : m_data (new Data(parent)) {
-
-    if (!parent.is_null()) {
-
-      m_data->depth = parent.m_data->depth + 1;
-
-      for (int i = 0; i < Dimension::value; i++)
-        m_data->global_coordinates[i] = (2 * parent.m_data->global_coordinates[i]) + local_coordinates[i];
-
-    }
-    else
-      for (int i = 0; i < Dimension::value; i++)
-        m_data->global_coordinates[i] = 0;
-  }
-
-  void free() { delete m_data; }
-
-  Node deep_copy(Self parent = Node()) const
-  {
-    if (is_null())
-      return Node();
-
-    Node out;
-    out.m_data = new Data(parent);
-
-    out.m_data->points = m_data->points;
-    out.m_data->depth = m_data->depth;
-    out.m_data->global_coordinates = m_data->global_coordinates;
-    std::unique_ptr<Children> children;
-    if (!is_leaf())
-    {
-      out.m_data->children = std::make_unique<Children>();
-      for (int index = 0; index < Degree::value; index++)
-        (*out.m_data->children)[index] = (*this)[index].deep_copy(out);
-    }
-    return out;
-  }
-
-  /// @}
-
-  /// \name Mutators
-  /// @{
-
-  /*!
-    \brief splits a node into subnodes.
-
-    Only leaf nodes should be split.
-    When a node is split it is no longer a leaf node.
-    A number of `Degree::value` children are constructed automatically, and their values are set.
-    Contents of this node are _not_ propagated automatically.
-    It is the responsibility of the caller to redistribute the points contained by a node after splitting
-   */
-  void split() {
-
-    CGAL_precondition (is_leaf());
-
-    m_data->children = std::make_unique<Children>();
-    for (int index = 0; index < Degree::value; index++) {
-
-      (*m_data->children)[index] = std::move(Self(*this, {Local_coordinates(index)}));
-
-    }
-  }
-
-  /*!
-   * \brief eliminates this node's children, making it a leaf node.
-   *
-   * When a node is un-split, its children are automatically deleted.
-   * After un-splitting a node it will be considered a leaf node.
-   */
-  void unsplit() {
-
-    m_data->children.reset();
-  }
-
-  /// @}
-
-  /// \endcond
-
-
-public:
-
-  /// \cond SKIP_IN_MANUAL
-  // Default creates null node
-  Node() : m_data(nullptr) { }
-
-  // Comparison operator
-  bool operator< (const Node& other) const { return m_data < other.m_data; }
-  /// \endcond
-
-  /// \name Type & Location
-  /// @{
-
-  /*!
-    \brief returns `true` if the node is null, `false` otherwise.
-  */
-  bool is_null() const { return (m_data == nullptr); }
-
-  /*!
-    \brief returns `true` if the node has no parent, `false` otherwise.
-    \pre `!is_null()`
-  */
-  bool is_root() const
-  {
-    CGAL_precondition(!is_null());
-    return m_data->parent.is_null();
-  }
-
-  /*!
-    \brief returns `true` if the node has no children, `false` otherwise.
-    \pre `!is_null()`
-  */
-  bool is_leaf() const
-  {
-    CGAL_precondition(!is_null());
-    return (!m_data->children);
-  }
-
-  /*!
-    \brief returns this node's depth.
-    \pre `!is_null()`
-  */
-  std::uint8_t depth() const
-  {
-    CGAL_precondition (!is_null());
-    return m_data->depth;
-  }
-
-  /*!
-    \brief returns this node's local coordinates (in relation to its parent).
-    \pre `!is_null()`
-  */
-  Local_coordinates local_coordinates() const {
-
-    CGAL_precondition (!is_null());
-    // TODO: There must be a better way of doing this!
-
-    Local_coordinates result;
-
-    for (std::size_t i = 0; i < Dimension::value; ++ i)
-      result[i] = global_coordinates()[i] & 1;
-
-    return result;
-  }
-
-  /*!
-    \brief returns this node's global coordinates.
-    \pre `!is_null()`
-  */
-  Global_coordinates global_coordinates() const
-  {
-    CGAL_precondition (!is_null());
-    return m_data->global_coordinates;
-  }
-
-
-  /// \name Adjacency
-  /// @{
-
-  /*!
-    \brief returns this node's parent.
-    \pre `!is_null()`
-  */
-  Self parent() const
-  {
-    CGAL_precondition (!is_null());
-    return m_data->parent;
-  }
-
-  /*!
-    \brief returns the nth child of this node.
-
-    \pre `!is_null()`
-    \pre `!is_leaf()`
-    \pre `0 <= index && index < Degree::value`
-
-    The orthtree subdivides the space in 2 on each dimension
-    available, so a child node can be accessed by selecting a Boolean
-    value for each dimension. The `index` parameter is thus
-    interpreted as a bitmap, where each bit matches a dimension
-    (starting by the least significant bit for coordinate X).
-
-    For example, in the case of an octree (dimension 3):
-
-    - index 0 (000 in binary) is the children on the "minimum corner" (xmin, ymin, zmin)
-    - index 1 (001 in binary) is on (xmax, ymin, zmin)
-    - index 2 (010 in binary) is on (xmin, ymax, zmin)
-    - index 6 (101 in binary) is on (xmax, ymin, zmax)
-
-    Diagram of a quadtree:
-
-    ```
-    Children of the current node:
-
-    Y axis
-    A
-    |  +-------------------+-------------------+
-    |  | Coord:  Ymax Xmin | Coord:  Ymax Xmax |
-    |  | Bitmap:    1    0 | Bitmap:    1    1 |
-    |  |                   |                   |
-    |  | -> index = 2      | -> index = 3      |
-    |  |                   |                   |
-    |  |                   |                   |
-    |  |                   |                   |
-    |  |                   |                   |
-    |  +-------------------+-------------------+
-    |  | Coord:  Ymin Xmin | Coord:  Ymin Xmax |
-    |  | Bitmap:    0    0 | Bitmap:    0    1 |
-    |  |                   |                   |
-    |  | -> index = 0      | -> index = 1      |
-    |  |                   |                   |
-    |  |                   |                   |
-    |  |                   |                   |
-    |  |                   |                   |
-    |  +-------------------+-------------------+
-    |
-    +--------------------------------------------> X axis
-    0
-    ```
-
-    The operator can be chained. For example, `n[5][2][3]` returns the
-    third child of the second child of the fifth child of a node `n`.
-  */
-  Self operator[](std::size_t index) const {
-
-    CGAL_precondition (!is_null());
-    CGAL_precondition (!is_leaf());
-    CGAL_precondition (index < Degree::value);
-
-    return (*m_data->children)[index];
-  }
-
-  /*!
-    \brief finds the directly adjacent node in a specific direction
-
-    \pre `!is_null()`
-    \pre `direction.to_ulong < 2 * Dimension::value`
-
-    Adjacent nodes are found according to several properties:
-    - adjacent nodes may be larger than the seek node, but never smaller
-    - a node has at most `2 * Dimension::value` different adjacent nodes (in 3D: left, right, up, down, front, back)
-    - adjacent nodes are not required to be leaf nodes
-
-    Here's a diagram demonstrating the concept for a Quadtree:
-
-    ```
-    +---------------+---------------+
-    |               |               |
-    |               |               |
-    |               |               |
-    |       A       |               |
-    |               |               |
-    |               |               |
-    |               |               |
-    +-------+-------+---+---+-------+
-    |       |       |   |   |       |
-    |   A   |  (S)  +---A---+       |
-    |       |       |   |   |       |
-    +---+---+-------+---+---+-------+
-    |   |   |       |       |       |
-    +---+---+   A   |       |       |
-    |   |   |       |       |       |
-    +---+---+-------+-------+-------+
-    ```
-
-    - (S) : Seek node
-    - A  : Adjacent node
-
-    Note how the top adjacent node is larger than the seek node.  The
-    right adjacent node is the same size, even though it contains
-    further subdivisions.
-
-    This implementation returns the adjacent node if it's found.  If
-    there is no adjacent node in that direction, it returns a null
-    node.
-
-    \param direction which way to find the adjacent node relative to
-    this one. Each successive bit selects the direction for the
-    corresponding dimension: for an Octree in 3D, 010 means: negative
-    direction in X, position direction in Y, negative direction in Z.
-
-    \return the adjacent node if it exists, a null node otherwise.
-  */
-  Self adjacent_node (Local_coordinates direction) const
-  {
-    CGAL_precondition(!is_null());
-
-    // Direction:   LEFT  RIGHT  DOWN    UP  BACK FRONT
-    // direction:    000    001   010   011   100   101
-
-    // Nodes only have up to 2*dim different adjacent nodes (since cubes have 6 sides)
-    CGAL_precondition(direction.to_ulong() < Dimension::value * 2);
-
-    // The root node has no adjacent nodes!
-    if (is_root())
-      return Self();
-
-    // The least significant bit indicates the sign (which side of the node)
-    bool sign = direction[0];
-
-    // The first two bits indicate the dimension/axis (x, y, z)
-    uint8_t dimension = uint8_t((direction >> 1).to_ulong());
-
-    // Create an offset so that the bit-significance lines up with the dimension (e.g. 1, 2, 4 --> 001, 010, 100)
-    int8_t offset = (uint8_t) 1 << dimension;
-
-    // Finally, apply the sign to the offset
-    offset = (sign ? offset : -offset);
-
-    // Check if this child has the opposite sign along the direction's axis
-    if (local_coordinates()[dimension] != sign) {
-
-      // This means the adjacent node is a direct sibling, the offset can be applied easily!
-      return parent()[local_coordinates().to_ulong() + offset];
-    }
-
-    // Find the parent's neighbor in that direction if it exists
-    Self adjacent_node_of_parent = parent().adjacent_node(direction);
-
-    // If the parent has no neighbor, then this node doesn't have one
-    if (adjacent_node_of_parent.is_null())
-      return Node();
-
-    // If the parent's adjacent node has no children, then it's this node's adjacent node
-    if (adjacent_node_of_parent.is_leaf())
-      return adjacent_node_of_parent;
-
-    // Return the nearest node of the parent by subtracting the offset instead of adding
-    return adjacent_node_of_parent[local_coordinates().to_ulong() - offset];
-
-  }
-
-  /*!
-    \brief equivalent to `adjacent_node()`, with an adjacency direction
-    rather than a bitset.
-   */
-  Self adjacent_node(Adjacency adjacency) const {
-    return adjacent_node(std::bitset<Dimension::value>(static_cast<int>(adjacency)));
-  }
-
-  /// @}
-
-  /// \name Point Range
-  /// @{
-
-  /*!
-    \brief checks whether the node is empty of points or not.
-   */
-  bool empty() const {
-    return m_data->points.empty();
-  }
-
-  /*!
-    \brief returns the number of points of this node.
-   */
-  std::size_t size() const {
-    return std::size_t(std::distance(m_data->points.begin(), m_data->points.end()));
-  }
-
-  /*!
-    \brief returns the iterator at the start of the collection of
-    points held by this node.
-   */
-  const_iterator begin() const { return m_data->points.begin(); }
-
-  /*!
-    \brief returns the iterator at the end of the collection of
-    points held by this node.
-   */
-  const_iterator end() const { return m_data->points.end(); }
-
-  /// @}
-
-  /// \cond SKIP_IN_MANUAL
-  /// \name Operators
-  /// @{
-
-  /*!
-   * \brief compares the topology of this node to another node.
-   *
-   * \todo
-   *
-   * \param rhs node to compare with
-   * \return whether the nodes have different topology.
-   */
-  bool operator==(const Self &rhs) const {
-    return m_data == rhs.m_data;
-  }
-
-  static bool is_topology_equal (const Self& a, const Self& b)
-  {
-    CGAL_assertion (!a.is_null() && !b.is_null());
-
-    // If one node is a leaf, and the other isn't, they're not the same
-    if (a.is_leaf() != b.is_leaf())
-      return false;
-
-    // If both nodes are non-leaf nodes
-    if (!a.is_leaf()) {
-
-      // Check all the children
-      for (int i = 0; i < Degree::value; ++i) {
-        // If any child cell is different, they're not the same
-        if (!is_topology_equal(a[i], b[i]))
-          return false;
-      }
-    }
-
-    // If both nodes are leaf nodes, they must be in the same location
-    return (a.global_coordinates() == b.global_coordinates());
-  }
-
-  friend std::ostream& operator<< (std::ostream& os, const Self& node)
-  {
-    return internal::print_orthtree_node(os, node);
-  }
-  /// \endcond
-};
-
-}
-
-#endif //CGAL_ORTHTREE_NODE_H

Orthtree/include/CGAL/Orthtree/Split_predicates.h

@@ -18,33 +18,40 @@
 
 namespace CGAL {
 
+template <typename GeomTraits>
+class Orthtree;
+
 namespace Orthtrees {
 
 /*!
   \ingroup PkgOrthtreeSplitPredicates
-  \brief A class used to choose when a node should be split depending on the number of inliers.
+  \brief A class used to choose when a node should be split depending on the number of contained elements.
 
   This is a bucket size predicate that considers a node should be
   split if it contains more than a certain number of items.
+
+  \warning This split predicate is only appropriate for trees with traits classes which are models of `OrthtreeTraitsWithData`
+  and where `Node_data` is a model of `Range`. `RandomAccessRange` is suggested for performance.
  */
-class Maximum_number_of_inliers {
+class Maximum_contained_elements {
 
   std::size_t m_bucket_size;
 
 public:
   /*!
-    \brief creates a predicate based on the number of inliers (bucket size).
+    \brief creates a predicate based on the number of contained elements.
    */
-  Maximum_number_of_inliers(std::size_t bucket_size) :
+  Maximum_contained_elements(std::size_t bucket_size) :
           m_bucket_size(bucket_size) {}
 
   /*!
-    \brief returns `true` if `n` should be split, `false` otherwise.
+    \brief returns `true` if the node with index `i` should be split, `false` otherwise.
    */
-  template<typename Node>
-  bool operator()(const Node &n) const {
-    return (n.size() > m_bucket_size);
+  template<typename GeomTraits>
+  bool operator()(typename Orthtree<GeomTraits>::Node_index i, const Orthtree<GeomTraits> &tree) const {
+    return (tree.data(i).size() > m_bucket_size);
   }
+
 };
 
 /*!
@@ -65,30 +72,33 @@ public:
   Maximum_depth(std::size_t max_depth) : m_max_depth(max_depth) {}
 
   /*!
-    \brief returns `true` if `n` should be split, `false` otherwise.
+    \brief returns `true` if the node with index `i` should be split, `false` otherwise.
    */
-  template<typename Node>
-  bool operator()(const Node &n) const {
-    return n.depth() < m_max_depth;
+  template<typename GeomTraits>
+  bool operator()(typename Orthtree<GeomTraits>::Node_index i, const Orthtree<GeomTraits> &tree) const {
+    return (tree.depth(i) < m_max_depth);
   }
+
 };
 
 /*!
   \ingroup PkgOrthtreeSplitPredicates
-  \brief A class used to choose when a node should be split depending on the depth and the number of inliers.
+  \brief A class used to choose when a node should be split depending on the depth and the number of contained elements.
 
   This predicate makes a note split if it contains more than a
   certain number of items and if its depth is lower than a certain
   limit.
 
   The refinement is stopped as soon as one of the conditions is
-  violated: if a node has more inliers than `bucket_size` but is
+  violated: if a node has more elements than `bucket_size` but is
   already at `max_depth`, it is not split. Similarly, a node that is
-  at a depth smaller than `max_depth` but already has fewer inliers
+  at a depth smaller than `max_depth` but already has fewer elements
   than `bucket_size`, it is not split.
 
+  \warning This split predicate is only appropriate for trees with traits classes which are models of `OrthtreeTraitsWithData`
+  and where `Node_data` is a model of `Range`. `RandomAccessRange` is suggested for performance.
  */
-class Maximum_depth_and_maximum_number_of_inliers {
+class Maximum_depth_and_maximum_contained_elements {
 
   std::size_t m_max_depth, m_bucket_size;
 
@@ -96,18 +106,19 @@ public:
 
   /*!  \brief creates a predicate using maximum depth or bucket size.
    */
-  Maximum_depth_and_maximum_number_of_inliers(std::size_t max_depth, std::size_t bucket_size) :
+  Maximum_depth_and_maximum_contained_elements(std::size_t max_depth, std::size_t bucket_size) :
           m_max_depth(max_depth), m_bucket_size(bucket_size) {}
 
   /*!
-    \brief returns `true` if `n` should be split, `false` otherwise.
+    \brief returns `true` if the node with index `i` should be split, `false` otherwise.
    */
-  template<typename Node>
-  bool operator()(const Node &n) const {
-    std::size_t num_points = n.size();
-    std::size_t depth = n.depth();
+  template<typename GeomTraits>
+  bool operator()(typename Orthtree<GeomTraits>::Node_index i, const Orthtree<GeomTraits> &tree) const {
+    std::size_t num_points = tree.data(i).size();
+    std::size_t depth = tree.depth(i);
     return (num_points > m_bucket_size && depth < m_max_depth);
   }
+
 };
 
 } // Orthtrees

Orthtree/include/CGAL/Orthtree/Traversal_iterator.h

@@ -14,7 +14,10 @@
 
 #include <CGAL/license/Orthtree.h>
 
+#include <optional>
+
 #include <boost/function.hpp>
+#include <boost/optional.hpp>
 #include <boost/iterator/iterator_facade.hpp>
 
 /// \cond SKIP_IN_MANUAL
@@ -22,18 +25,21 @@
 namespace CGAL {
 
 /*!
- * \ingroup PkgOrthtreeClasses
+ * \ingroup PkgOrthtreeTraversal
  *
- * \brief
+ * \brief Wraps a traversal definition to produce an iterator which traverses the tree when incremented.
  *
  * \todo
  *
- * \tparam Value
+ * \tparam Tree The orthtree type to iterate over
  */
-template<class Value>
-class Traversal_iterator :
-        public boost::iterator_facade<Traversal_iterator<Value>, Value, boost::forward_traversal_tag> {
-
+template <class Tree>
+class Index_traversal_iterator : public boost::iterator_facade<
+  Index_traversal_iterator<Tree>,
+  const typename Tree::Node_index,
+  boost::forward_traversal_tag,
+  const typename Tree::Node_index
+> {
 public:
 
   /// \name Types
@@ -44,7 +50,9 @@ public:
    *
    * \todo
    */
-  typedef std::function<Value(Value)> Traversal_function;
+  using Traversal_function = std::function<std::optional<std::size_t>(const Tree&, std::size_t)>;
+
+  using Node_index = typename Tree::Node_index;
 
   /// @}
 
@@ -54,44 +62,51 @@ public:
   /// @{
 
   /*!
-   * \brief
+   * \brief Default constructor, creates an end sentinel
    *
    * \todo
    */
-  Traversal_iterator() : m_value(), m_next() {}
+  Index_traversal_iterator() : m_next() {}
 
   /*!
    * \brief
    *
    * \todo
    *
+   * \param tree
    * \param first
    * \param next
    */
-  Traversal_iterator(Value first, const Traversal_function &next) : m_value(first), m_next(next) {}
+  Index_traversal_iterator(const Tree& tree, Node_index first, const Traversal_function& next) :
+    m_tree(&tree), m_index(first), m_next(next) {}
 
   /// @}
 
 private:
+
   friend class boost::iterator_core_access;
 
-  bool equal(Traversal_iterator<Value> const &other) const {
-    return m_value == other.m_value;
+  bool equal(Index_traversal_iterator<Tree> const& other) const {
+    return m_index == other.m_index;
   }
 
   void increment() {
-    m_value = m_next(m_value);
+    // invoking increment on the sentinel is undefined behavior
+    m_index = m_next(*m_tree, *m_index);
   }
 
-  Value &dereference() const {
-    return const_cast<Value&>(m_value);
+  Node_index dereference() const {
+    return *m_index;
   }
 
 private:
 
-  Value m_value;
+  const Tree* m_tree = nullptr;
+  std::optional<std::size_t> m_index;
   Traversal_function m_next;
+
 };
+
 }
 
 /// \endcond

Orthtree/include/CGAL/Orthtree/Traversals.h

@@ -22,159 +22,80 @@
 
 namespace CGAL {
 
-/// \cond SKIP_IN_MANUAL
-// Forward declaration
-template<typename T, typename PR, typename PM>
-class Orthtree;
-/// \endcond
-
 namespace Orthtrees {
 
-/// \cond SKIP_IN_MANUAL
-
-template <typename Node>
-Node next_sibling(Node n) {
-
-  // Passing null returns the first node
-  if (n.is_null())
-    return Node();
-
-  // If this node has no parent, it has no siblings
-  if (n.parent().is_null())
-    return Node();
-
-  // Find out which child this is
-  std::size_t index = n.local_coordinates().to_ulong();
-
-  constexpr static int degree = Node::Degree::value;
-  // Return null if this is the last child
-  if (int(index) == degree - 1)
-    return Node();
-
-  // Otherwise, return the next child
-  return n.parent()[index + 1];
-}
-
-template <typename Node>
-Node next_sibling_up(Node n) {
-
-  if (n.is_null())
-    return Node();
-
-  Node up = n.parent();
-
-  while (!up.is_null()) {
-
-    if (!next_sibling(up).is_null())
-      return next_sibling(up);
-
-    up = up.parent();
-  }
-
-  return Node();
-}
-
-template <typename Node>
-Node deepest_first_child(Node n) {
-
-  if (n.is_null())
-    return Node();
-
-  // Find the deepest child on the left
-  Node first = n;
-  while (!first.is_leaf())
-    first = first[0];
-  return first;
-}
-
-template <typename Node>
-Node first_child_at_depth(Node n, std::size_t depth) {
-
-  if (n.is_null())
-    return Node();
-
-  std::stack<Node> todo;
-  todo.push(n);
-
-  if (n.depth() == depth)
-    return n;
-
-  while (!todo.empty())
-  {
-    Node node = todo.top();
-    todo.pop();
-
-    if (node.depth() == depth)
-      return node;
-
-    if (!node.is_leaf())
-      for (int i = 0; i < Node::Degree::value; ++ i)
-        todo.push(node[std::size_t(Node::Degree::value - 1 - i)]);
-  }
-
-  return Node();
-}
-
-/// \endcond
-
 /*!
   \ingroup PkgOrthtreeTraversal
   \brief A class used for performing a preorder traversal.
 
+  \tparam Tree an instance of `Orthtree`
+
   A preorder traversal starts from the root towards the leaves.
 
   \cgalModels{OrthtreeTraversal}
  */
+template <typename Tree>
 struct Preorder_traversal {
+private:
 
-  template <typename Node>
-  Node first(Node root) const {
-    return root;
+  const Tree& m_orthtree;
+
+public:
+
+  using Node_index = typename Tree::Node_index;
+
+  Preorder_traversal(const Tree& orthtree) : m_orthtree(orthtree) {}
+
+  Node_index first_index() const {
+    return m_orthtree.root();
   }
 
-  template <typename Node>
-  Node next(Node n) const {
+  std::optional<Node_index> next_index(Node_index n) const {
 
-    if (n.is_leaf()) {
+    if (m_orthtree.is_leaf(n)) {
 
-      Node next = next_sibling(n);
+      auto next = m_orthtree.next_sibling(n);
 
-      if (next.is_null())
-        return next_sibling_up(n);
+      if (!next)
+        next = m_orthtree.next_sibling_up(n);
 
       return next;
 
+    } else {
+      return m_orthtree.child(n, 0);
     }
-    else // Return the first child of this node
-      return n[0];
   }
+
 };
 
 /*!
   \ingroup PkgOrthtreeTraversal
   \brief A class used for performing a postorder traversal.
 
+  \tparam Tree an instance of `Orthtree`
+
   A postorder traversal starts from the leaves towards the root.
 
   \cgalModels{OrthtreeTraversal}
  */
+template <typename Tree>
 struct Postorder_traversal {
+private:
 
-  template <typename Node>
-  Node first(Node root) const {
+  const Tree& m_orthtree;
 
-    return deepest_first_child(root);
+public:
+
+  using Node_index = typename Tree::Node_index;
+
+  Postorder_traversal(const Tree& orthtree) : m_orthtree(orthtree) {}
+
+  Node_index first_index() const {
+    return m_orthtree.deepest_first_child(m_orthtree.root());
   }
 
-  template <typename Node>
-  Node next(Node n) const {
-
-    Node next = deepest_first_child(next_sibling(n));
-
-    if (next.is_null())
-      next = n.parent();
-
-    return next;
+  std::optional<Node_index> next_index(Node_index n) const {
+    return m_orthtree.index(next(&m_orthtree[n]));
   }
 };
 
@@ -182,27 +103,37 @@ struct Postorder_traversal {
   \ingroup PkgOrthtreeTraversal
   \brief A class used for performing a traversal on leaves only.
 
-  All non-leave nodes are ignored.
+  \tparam Tree an instance of `Orthtree`
+
+  All non-leaf nodes are ignored.
 
   \cgalModels{OrthtreeTraversal}
  */
+template <typename Tree>
 struct Leaves_traversal {
+private:
 
-  template <typename Node>
-  Node first(Node root) const {
+  const Tree& m_orthtree;
 
-    return deepest_first_child(root);
+public:
+
+  using Node_index = typename Tree::Node_index;
+
+  Leaves_traversal(const Tree& orthtree) : m_orthtree(orthtree) {}
+
+  Node_index first_index() const {
+    return m_orthtree.deepest_first_child(m_orthtree.root());
   }
 
-  template <typename Node>
-  Node next(Node n) const {
+  std::optional<Node_index> next_index(Node_index n) const {
 
-    Node next = deepest_first_child(next_sibling(n));
+    if (m_orthtree.next_sibling(n))
+      return m_orthtree.deepest_first_child(*m_orthtree.next_sibling(n));
 
-    if (next.is_null())
-      next = deepest_first_child(next_sibling_up(n));
+    if (m_orthtree.next_sibling_up(n))
+      return m_orthtree.deepest_first_child(*m_orthtree.next_sibling_up(n));
 
-    return next;
+    return {};
   }
 };
 
@@ -210,52 +141,58 @@ struct Leaves_traversal {
   \ingroup PkgOrthtreeTraversal
   \brief A class used for performing a traversal of a specific depth level.
 
-  All trees at another depth are ignored. If the selected depth is
+  \tparam Tree an instance of `Orthtree`
+
+  All tree nodes at another depth are ignored. If the selected depth is
+  All tree nodes at another depth are ignored. If the selected depth is
   higher than the maximum depth of the orthtree, no node will be traversed.
 
   \cgalModels{OrthtreeTraversal}
  */
+template <typename Tree>
 struct Level_traversal {
-
 private:
 
-  const std::size_t depth;
+  const Tree& m_orthtree;
+  const std::size_t m_depth;
 
 public:
 
+  using Node_index = typename Tree::Node_index;
+
   /*!
     constructs a `depth`-level traversal.
   */
-  Level_traversal (std::size_t depth) : depth(depth) { }
+  Level_traversal(const Tree& orthtree, std::size_t depth) : m_orthtree(orthtree), m_depth(depth) {}
 
-  template <typename Node>
-  Node first(Node root) const {
-    return first_child_at_depth(root, depth);
+  Node_index first_index() const {
+    // assumes the tree has at least one child at m_depth
+    return *m_orthtree.first_child_at_depth(m_orthtree.root(), m_depth);
   }
 
-  template <typename Node>
-  Node next(Node n) const {
-    std::cerr << depth << " ";
-    Node next = next_sibling(n);
+  std::optional<Node_index> next_index(Node_index n) const {
 
-    if (next.is_null())
-    {
-      Node up = n;
-      do
-      {
-        up = next_sibling_up(up);
-        if (up.is_null())
-          return Node();
-        next = first_child_at_depth(up, depth);
-      }
-      while (next.is_null());
+    auto next = m_orthtree.next_sibling(n);
+
+    if (!next) {
+
+      auto up = n;
+      do {
+
+        if (!m_orthtree.next_sibling_up(up))
+          return {};
+
+        up = *m_orthtree.next_sibling_up(up);
+        next = m_orthtree.first_child_at_depth(up, m_depth);
+
+      } while (!next);
     }
 
     return next;
   }
 };
 
-} // Orthtree
+} // Orthtrees
 } // CGAL
 
 #endif //CGAL_ORTHTREE_TRAVERSALS_H

Orthtree/include/CGAL/Orthtree_traits.h

@@ -0,0 +1,61 @@
+// Copyright (c) 2023  INRIA (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)     : Sven Oesau
+
+#ifndef ORTHTREE_TESTS_ORTHTREE_TRAITS_H
+#define ORTHTREE_TESTS_ORTHTREE_TRAITS_H
+
+#include <CGAL/license/Orthtree.h>
+
+#include <CGAL/Dimension.h>
+#include <CGAL/Bbox_2.h>
+#include <CGAL/Point_set_2.h>
+#include <CGAL/Orthtree/Cartesian_ranges.h>
+
+#include <CGAL/Orthtree_traits_base.h>
+
+namespace CGAL {
+
+/*!
+  \ingroup PkgOrthtreeTraits
+
+  Traits class for defining an orthtree using the class `CGAL::Orthtree` without storing data in the nodes.
+
+  \tparam GeomTraits model of `Kernel`.
+  \tparam dimension the dimension of the ambient Euclidean space.
+
+  \cgalModels{OrthtreeTraits}
+  \sa `CGAL::Octree`
+  \sa `CGAL::Quadtree`
+  \sa `CGAL::Orthtree_traits_base<GeomTraits, DimensionTag>`
+*/
+template <typename GeomTraits, int dimension>
+struct Orthtree_traits : public Orthtree_traits_base<GeomTraits, dimension> {
+public:
+  using Base = Orthtree_traits_base<GeomTraits, dimension>;
+  using Self = Orthtree_traits<GeomTraits, dimension>;
+  using Tree = Orthtree<Self>;
+
+  using Node_index = typename Base::Node_index;
+
+  Orthtree_traits() {}
+
+  auto construct_root_node_bbox_object() const {
+    return [&]() -> typename Self::Bbox_d {
+      return {std::apply(Self::construct_point_d_object(), std::array<typename Self::FT, Self::dimension>{-1.0, -1.0, -1.0}),
+              std::apply(Self::construct_point_d_object(), std::array<typename Self::FT, Self::dimension>{1.0, 1.0, 1.0})};
+    };
+  }
+};
+
+}
+
+
+#endif //ORTHTREE_TESTS_ORTHTREE_TRAITS_H

Orthtree/include/CGAL/Orthtree_traits_2.h

@@ -1,143 +0,0 @@
-// Copyright (c) 2020  GeometryFactory (France).
-// All rights reserved.
-//
-// This file is part of CGAL (www.cgal.org).
-//
-// $URL$
-// $Id$
-// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
-//
-// Author(s)     : Simon Giraudot
-
-#ifndef CGAL_ORTHTREE_TRAITS_2_H
-#define CGAL_ORTHTREE_TRAITS_2_H
-
-#include <CGAL/license/Orthtree.h>
-
-#include <CGAL/Dimension.h>
-#include <CGAL/Bbox_2.h>
-
-namespace CGAL
-{
-
-/*!
-  \ingroup PkgOrthtreeTraits
-
-  The class `Orthtree_traits_2` can be used as a template parameter of
-  the `Orthtree` class.
-
-  \tparam GeomTraits model of `Kernel`.
-
-  \cgalModels{OrthtreeTraits}
-  \sa `CGAL::Quadtree`
-  \sa `CGAL::Orthtree_traits_3`
-  \sa `CGAL::Orthtree_traits_d`
-*/
-template <typename GeomTraits>
-struct Orthtree_traits_2
-{
-public:
-
-  /// \name Types
-  /// @{
-
-  typedef Dimension_tag<2> Dimension; ///< Dimension type.
-  typedef Bbox_2 Bbox_d; ///< Bounding box type.
-  typedef typename GeomTraits::FT FT; ///< Number type.
-  typedef typename GeomTraits::Point_2 Point_d; ///< Point type.
-  typedef typename GeomTraits::Circle_2 Sphere_d; ///< Sphere type.
-  typedef typename GeomTraits::Cartesian_const_iterator_2 Cartesian_const_iterator_d; ///< An iterator over the %Cartesian coordinates.
-  typedef std::array<FT, Dimension::value> Array; ///< Array type.
-
-  /*!
-   * \brief Two directions along each axis in Cartesian space, relative to a node.
-   *
-   * Directions are mapped to numbers as 2-bit integers.
-   *
-   * The first bit indicates the axis (0 = x, 1 = y),
-   * the second bit indicates the direction along that axis (0 = -, 1 = +).
-   *
-   * The following diagram may be a useful reference:
-   *
-   *            3 *
-   *              |
-   *              |                    y+
-   *              |                    *
-   *     0 *------+------* 1           |
-   *              |                    |
-   *              |                    +-----* x+
-   *              |
-   *              * 2
-   *
-   * This lookup table may also be helpful:
-   *
-   * | Direction | bitset | number | Enum  |
-   * | --------- | ------ | ------ | ----- |
-   * | `-x`      |  00    | 0      | LEFT  |
-   * | `+x`      |  01    | 1      | RIGHT |
-   * | `-y`      |  10    | 2      | DOWN  |
-   * | `+y`      |  11    | 3      | UP    |
-   */
-  enum Adjacency
-  {
-    LEFT,
-    RIGHT,
-    DOWN,
-    UP
-  };
-
-#ifdef DOXYGEN_RUNNING
-  /*!
-    Functor with an operator to construct a `Point_d` from an `Array` object.
-  */
-  typedef unspecified_type Construct_point_d_from_array;
-#else
-  struct Construct_point_d_from_array
-  {
-    Point_d operator() (const Array& array) const
-    {
-      return Point_d (array[0], array[1]);
-    }
-  };
-#endif
-
-
-#ifdef DOXYGEN_RUNNING
-  /*!
-    Functor with an operator to construct a `Bbox_d` from two `Array` objects (coordinates of minimum and maximum points).
-  */
-  typedef unspecified_type Construct_bbox_d;
-#else
-  struct Construct_bbox_d
-  {
-    Bbox_d operator() (const Array& min,
-                       const Array& max) const
-    {
-      return Bbox_d (min[0], min[1], max[0], max[1]);
-    }
-  };
-#endif
-
-  /// @}
-
-  /// \name Operations
-  /// @{
-
-  /*!
-    Function used to construct an object of type `Construct_point_d_from_array`.
-  */
-  Construct_point_d_from_array construct_point_d_from_array_object() const
-  { return Construct_point_d_from_array(); }
-
-  /*!
-    Function used to construct an object of type `Construct_bbox_d`.
-  */
-  Construct_bbox_d construct_bbox_d_object() const
-  { return Construct_bbox_d(); }
-
-  /// @}
-};
-
-}
-
-#endif // CGAL_ORTHTREE_TRAITS_2_H

Orthtree/include/CGAL/Orthtree_traits_3.h

@@ -1,160 +0,0 @@
-// Copyright (c) 2020  GeometryFactory (France).
-// All rights reserved.
-//
-// This file is part of CGAL (www.cgal.org).
-//
-// $URL$
-// $Id$
-// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
-//
-// Author(s)     : Simon Giraudot
-
-#ifndef CGAL_ORTHTREE_TRAITS_3_H
-#define CGAL_ORTHTREE_TRAITS_3_H
-
-#include <CGAL/license/Orthtree.h>
-
-#include <CGAL/Dimension.h>
-#include <CGAL/Bbox_3.h>
-
-namespace CGAL
-{
-
-/*!
-  \ingroup PkgOrthtreeTraits
-
-  The class `Orthtree_traits_3` can be used as a template parameter of
-  the `Orthtree` class.
-
-  \tparam GeomTraits model of `Kernel`.
-
-  \cgalModels{OrthtreeTraits}
-  \sa `CGAL::Octree`
-  \sa `CGAL::Orthtree_traits_2`
-  \sa `CGAL::Orthtree_traits_d`
-*/
-template <typename GeomTraits>
-struct Orthtree_traits_3
-{
-public:
-
-  /// \name Types
-  /// @{
-
-  typedef Dimension_tag<3> Dimension; ///< Dimension type.
-  typedef Bbox_3 Bbox_d; ///< Bounding box type.
-  typedef typename GeomTraits::FT FT; ///< Number type.
-  typedef typename GeomTraits::Point_3 Point_d; ///< Point type.
-  typedef typename GeomTraits::Sphere_3 Sphere_d; ///< Sphere type.
-  typedef typename GeomTraits::Cartesian_const_iterator_3 Cartesian_const_iterator_d; ///< An iterator over the %Cartesian coordinates.
-  typedef std::array<FT, Dimension::value> Array; ///< Array type.
-
-  /*!
-   * \brief Two directions along each axis in Cartesian space, relative to a node.
-   *
-   * Directions are mapped to numbers as 3-bit integers,
-   * though the numbers 6 and 7 are not used because there are only 6 different directions.
-   *
-   * The first two bits indicate the axis (00 = x, 01 = y, 10 = z),
-   * the third bit indicates the direction along that axis (0 = -, 1 = +).
-   *
-   * The following diagram may be a useful reference:
-   *
-   *            3 *
-   *              |  * 5
-   *              | /                  y+
-   *              |/                   *  z+
-   *     0 *------+------* 1           | *
-   *             /|                    |/
-   *            / |                    +-----* x+
-   *         4 *  |
-   *              * 2
-   *
-   * This lookup table may also be helpful:
-   *
-   * | Direction | bitset | number | Enum  |
-   * | --------- | ------ | ------ | ----- |
-   * | `-x`      | 000    | 0      | LEFT  |
-   * | `+x`      | 001    | 1      | RIGHT |
-   * | `-y`      | 010    | 2      | DOWN  |
-   * | `+y`      | 011    | 3      | UP    |
-   * | `-z`      | 100    | 4      | BACK  |
-   * | `+z`      | 101    | 5      | FRONT |
-   */
-  enum Adjacency
-  {
-    LEFT,
-    RIGHT,
-    DOWN,
-    UP,
-    BACK,
-    FRONT
-  };
-
-  /// \cond SKIP_IN_MANUAL
-  enum Child {
-    LEFT_BOTTOM_BACK,
-    RIGHT_BOTTOM_BACK,
-    LEFT_TOP_BACK,
-    RIGHT_TOP_BACK,
-    LEFT_BOTTOM_FRONT,
-    RIGHT_BOTTOM_FRONT,
-    LEFT_TOP_FRONT,
-    RIGHT_TOP_FRONT
-  };
-  /// \endcond
-
-#ifdef DOXYGEN_RUNNING
-  /*!
-    Functor with an operator to construct a `Point_d` from an `Array` object.
-  */
-  typedef unspecified_type Construct_point_d_from_array;
-#else
-  struct Construct_point_d_from_array
-  {
-    Point_d operator() (const Array& array) const
-    {
-      return Point_d (array[0], array[1], array[2]);
-    }
-  };
-#endif
-
-#ifdef DOXYGEN_RUNNING
-  /*!
-    Functor with an operator to construct a `Bbox_d` from two `Array` objects (coordinates of minimum and maximum points).
-  */
-  typedef unspecified_type Construct_bbox_d;
-#else
-  struct Construct_bbox_d
-  {
-    Bbox_d operator() (const Array& min,
-                       const Array& max) const
-    {
-      return Bbox_d (min[0], min[1], min[2], max[0], max[1], max[2]);
-    }
-  };
-#endif
-
-  /// @}
-
-  /// \name Operations
-  /// @{
-
-  /*!
-    Function used to construct an object of type `Construct_point_d_from_array`.
-  */
-  Construct_point_d_from_array construct_point_d_from_array_object() const
-  { return Construct_point_d_from_array(); }
-
-  /*!
-    Function used to construct an object of type `Construct_bbox_d`.
-  */
-  Construct_bbox_d construct_bbox_d_object() const
-  { return Construct_bbox_d(); }
-
-  /// @}
-};
-
-}
-
-#endif // CGAL_ORTHTREE_TRAITS_3_H

Orthtree/include/CGAL/Orthtree_traits_base.h

@@ -0,0 +1,167 @@
+// Copyright (c) 2023  INRIA (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)     : Jackson Campolattaro
+
+#ifndef ORTHTREE_ORTHTREE_TRAITS_BASE_H
+#define ORTHTREE_ORTHTREE_TRAITS_BASE_H
+
+#include <CGAL/license/Orthtree.h>
+
+#include <CGAL/Bbox_2.h>
+#include <CGAL/Point_set_2.h>
+#include <CGAL/Orthtree/Cartesian_ranges.h>
+#include <CGAL/Orthtree.h>
+
+namespace CGAL {
+
+/*!
+  \ingroup PkgOrthtreeTraits
+
+  The class `Orthtree_traits_base` is a base class providing common choices for types and functors.
+  The base class is extended by `CGAL::Orthtree_traits_point<GeomTraits, PointRange, PointMap, dimension>` and by `CGAL::Orthtree_traits_face_graph<PolygonMesh, VertexPointMap>`.
+
+  \tparam GeomTraits a model of `Kernel`.
+  \tparam dim dimension of the ambient Euclidean space.
+
+  \sa `CGAL::Orthtree_traits_point<GeomTraits, PointRange, PointMap, dim>`
+  \sa `CGAL::Orthtree_traits_face_graph<PolygonMesh, VertexPointMap>`
+*/
+
+template <typename GeomTraits, int dim>
+struct Orthtree_traits_base {
+  /// \name Types
+  /// @{
+  using Node_index = std::size_t;
+  using Kernel = GeomTraits;
+  static constexpr int dimension = dim;
+  using FT = typename GeomTraits::FT;
+  using Point_d = typename GeomTraits::Point_d;
+  using Bbox_d = typename GeomTraits::Iso_box_d;
+  using Sphere_d = typename GeomTraits::Sphere_d;
+  using Cartesian_const_iterator_d = typename GeomTraits::Cartesian_const_iterator_d;
+  /*!
+   * Adjacency type.
+   *
+   * \note This type is used to identify adjacency directions with
+   * easily understandable keywords (left, right, up, down, ...) and is thus
+   * mainly useful in 2D and 3D. In
+   * higher dimensions, such keywords do not exist and this type is
+   * simply an integer. Conversions from this integer to bitsets still
+   * work but do not provide any user-friendly API for adjacency selection.
+   *
+   * Two directions along each axis in %Cartesian space, relative to a node.
+   *
+   * Directions are mapped to numbers as 3-bit integers in the 3D case or as 2-bit integers in the 2D case.
+   * In the 3d case the numbers 6 and 7 are not used because there are only 6 different directions.
+   *
+   * The first two bits indicate the axis (00 = x, 01 = y, 10 = z),
+   * the third bit indicates the direction along that axis (0 = -, 1 = +).
+   *
+   * The following diagram and table showing the 3D case may be a useful reference (2D case is identical with one dimension less):
+   *
+   *            3 *
+   *              |  * 4
+   *              | /                  y+
+   *              |/                   *
+   *     0 *------+------* 1           |
+   *             /|                    |
+   *            / |                    +-----* x+
+   *         5 *  |                   /
+   *              * 2                /
+   *                                * z+
+   *
+   * This lookup table may also be helpful:
+   *
+   * | Direction | bitset | number | Enum  |
+   * | --------- | ------ | ------ | ----- |
+   * | `-x`      | 000    | 0      | LEFT  |
+   * | `+x`      | 001    | 1      | RIGHT |
+   * | `-y`      | 010    | 2      | DOWN  |
+   * | `+y`      | 011    | 3      | UP    |
+   * | `-z`      | 100    | 4      | BACK  |
+   * | `+z`      | 101    | 5      | FRONT |
+   */
+  using Adjacency = int;
+  /// @}
+
+  auto construct_point_d_object() const {
+    return [](auto... Args) -> Point_d {
+      std::initializer_list<FT> args_list{Args...};
+      return Point_d{static_cast<int>(args_list.size()), args_list.begin(), args_list.end()};
+    };
+  }
+};
+
+template <typename GeomTraits>
+struct Orthtree_traits_base<GeomTraits, 2> {
+  using Node_index = std::size_t;
+  using Kernel = GeomTraits;
+  static constexpr int dimension = 2;
+  using FT = typename GeomTraits::FT;
+  using Point_d = typename GeomTraits::Point_2;
+  using Bbox_d = typename GeomTraits::Iso_rectangle_2;
+  using Sphere_d = typename GeomTraits::Circle_2;
+  using Cartesian_const_iterator_d = typename GeomTraits::Cartesian_const_iterator_2;
+
+  enum Adjacency {
+    LEFT,
+    RIGHT,
+    DOWN,
+    UP
+  };
+
+  auto construct_point_d_object() const {
+    return [](const FT& x, const FT& y) -> Point_d {
+      return {x, y};
+    };
+  }
+};
+
+template <typename GeomTraits>
+struct Orthtree_traits_base<GeomTraits, 3> {
+  using Node_index = std::size_t;
+  using Kernel = GeomTraits;
+  static constexpr int dimension = 3;
+  using FT = typename GeomTraits::FT;
+  using Point_d = typename GeomTraits::Point_3;
+  using Bbox_d = typename GeomTraits::Iso_cuboid_3;
+  using Sphere_d = typename GeomTraits::Sphere_3;
+  using Cartesian_const_iterator_d = typename GeomTraits::Cartesian_const_iterator_3;
+
+  enum Adjacency {
+    LEFT,
+    RIGHT,
+    DOWN,
+    UP,
+    BACK,
+    FRONT
+  };
+
+  enum Child {
+    LEFT_BOTTOM_BACK,
+    RIGHT_BOTTOM_BACK,
+    LEFT_TOP_BACK,
+    RIGHT_TOP_BACK,
+    LEFT_BOTTOM_FRONT,
+    RIGHT_BOTTOM_FRONT,
+    LEFT_TOP_FRONT,
+    RIGHT_TOP_FRONT
+  };
+
+  auto construct_point_d_object() const {
+    return [](const FT& x, const FT& y, const FT& z) -> Point_d {
+      return {x, y, z};
+    };
+  }
+};
+
+}
+
+#endif //ORTHTREE_ORTHTREE_TRAITS_BASE_H

Orthtree/include/CGAL/Orthtree_traits_d.h

@@ -1,137 +0,0 @@
-// Copyright (c) 2020  GeometryFactory (France).
-// All rights reserved.
-//
-// This file is part of CGAL (www.cgal.org).
-//
-// $URL$
-// $Id$
-// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
-//
-// Author(s)     : Simon Giraudot
-
-#ifndef CGAL_ORTHTREE_TRAITS_D_H
-#define CGAL_ORTHTREE_TRAITS_D_H
-
-#include <CGAL/license/Orthtree.h>
-
-#include <CGAL/Dimension.h>
-
-namespace CGAL
-{
-
-/*!
-  \ingroup PkgOrthtreeTraits
-
-  The class `Orthtree_traits_d` can be used as a template parameter of
-  the `Orthtree` class.
-
-  \tparam GeomTraits model of `Kernel`.
-  \tparam DimensionTag specialization of `CGAL::Dimension_tag`.
-
-  \cgalModels{OrthtreeTraits}
-  \sa `CGAL::Orthtree`
-  \sa `CGAL::Orthtree_traits_2`
-  \sa `CGAL::Orthtree_traits_3`
-*/
-
-template <typename GeomTraits, typename DimensionTag>
-struct Orthtree_traits_d
-{
-public:
-
-  /// \name Types
-  /// @{
-
-  typedef DimensionTag Dimension; ///< Dimension type.
-  typedef typename GeomTraits::FT FT; ///< Number type.
-  typedef typename GeomTraits::Point_d Point_d; ///< Point type.
-  typedef typename GeomTraits::Sphere_d Sphere_d; ///< Sphere type.
-  typedef typename GeomTraits::Cartesian_const_iterator_d Cartesian_const_iterator_d; ///< An iterator over the %Cartesian coordinates.
-  typedef std::array<FT, Dimension::value> Array; ///< Array type.
-
-#ifdef DOXYGEN_RUNNING
-  typedef unspecified_type Bbox_d; ///< Bounding box type.
-#else
-  class Bbox_d
-  {
-    Point_d m_min, m_max;
-  public:
-
-    Bbox_d (const Point_d& pmin, const Point_d& pmax)
-      : m_min (pmin), m_max (pmax)
-    { }
-
-    const Point_d& min BOOST_PREVENT_MACRO_SUBSTITUTION () { return m_min; }
-    const Point_d& max BOOST_PREVENT_MACRO_SUBSTITUTION () { return m_max; }
-  };
-#endif
-
-  /*!
-    Adjacency type.
-
-    \note This type is used to identify adjacency directions with
-    easily understandable keywords (left, right, up, etc.) and is thus
-    mainly useful for `Orthtree_traits_2` and `Orthtree_traits_3`. In
-    higher dimensions, such keywords do not exist and this type is
-    simply an integer. Conversions from this integer to bitsets still
-    works but do not provide any easier API for adjacency selection.
-  */
-  typedef int Adjacency;
-
-
-#ifdef DOXYGEN_RUNNING
-  /*!
-    Functor with an operator to construct a `Point_d` from an `Array` object.
-  */
-  typedef unspecified_type Construct_point_d_from_array;
-#else
-  struct Construct_point_d_from_array
-  {
-    Point_d operator() (const Array& array) const
-    {
-      return Point_d (array.begin(), array.end());
-    }
-  };
-#endif
-
-
-#ifdef DOXYGEN_RUNNING
-  /*!
-    Functor with an operator to construct a `Bbox_d` from two `Array` objects (coordinates of minimum and maximum points).
-  */
-  typedef unspecified_type Construct_bbox_d;
-#else
-  struct Construct_bbox_d
-  {
-    Bbox_d operator() (const Array& min,
-                       const Array& max) const
-    {
-      return Bbox_d (Point_d (min.begin(), min.end()),
-                     Point_d (max.begin(), max.end()));
-    }
-  };
-#endif
-
-  /// @}
-
-  /// \name Operations
-  /// @{
-
-  /*!
-    Function used to construct an object of type `Construct_point_d_from_array`.
-  */
-  Construct_point_d_from_array construct_point_d_from_array_object() const
-  { return Construct_point_d_from_array(); }
-
-  /*!
-    Function used to construct an object of type `Construct_bbox_d`.
-  */
-  Construct_bbox_d construct_bbox_d_object() const
-  { return Construct_bbox_d(); }
-
-  /// @}
-};
-
-}
-
-#endif // CGAL_ORTHTREE_TRAITS_D_H

Orthtree/include/CGAL/Orthtree_traits_face_graph.h

@@ -0,0 +1,168 @@
+// Copyright (c) 2023  INRIA (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Sebastien LORIOT
+
+
+#ifndef CGAL_ORTHREE_TRAITS_FACE_GRAPH_H
+#define CGAL_ORTHREE_TRAITS_FACE_GRAPH_H
+
+#include <CGAL/license/Orthtree.h>
+
+#include <CGAL/Orthtree_traits_base.h>
+
+#include <boost/graph/graph_traits.hpp>
+#include <CGAL/Dimension.h>
+#include <CGAL/Bbox_3.h>
+
+namespace CGAL {
+
+/*!
+\ingroup PkgOrthtreeTraits
+
+Traits class for the `Orthtree` class to be used to construct a 3D octree around
+a triangulated surface mesh. Each node of the octree will store all the faces of the
+mesh intersected by its bounding box. The subdivision of the octree is controlled
+by the nested class `Orthtree_traits_face_graph::Split_predicate_node_min_extent`
+to which the minimal extent of a node should be provided.
+
+\tparam TriangleMesh a model of `FaceListGraph` with all faces being triangles
+\tparam VertexPointMap a property map associating points to the vertices of `TriangleMesh`
+
+\cgalModels{OrthtreeTraitsWithData}
+*/
+template <class TriangleMesh, class VertexPointMap>
+struct Orthtree_traits_face_graph : public Orthtree_traits_base<
+  typename Kernel_traits<typename boost::property_traits<VertexPointMap>::value_type>::type, 3 > {
+
+  Orthtree_traits_face_graph(const TriangleMesh& pm, VertexPointMap vpm)
+    : m_pm(pm), m_vpm(vpm) {}
+
+  /// \name Types
+  /// @{
+
+  using Base = Orthtree_traits_base<
+  typename Kernel_traits<typename boost::property_traits<VertexPointMap>::value_type>::type, 3 >;
+  using Self = Orthtree_traits_face_graph<TriangleMesh, VertexPointMap>;
+  using Tree = Orthtree<Self>;
+
+  using Point_d = typename Self::Point_d;
+  using Bbox_d = typename Self::Bbox_d;
+  using FT = typename Self::FT;
+  using Cartesian_const_iterator_d = typename Self::Cartesian_const_iterator_d;
+
+  using Node_index = typename Base::Node_index;
+  using Node_data = std::vector<typename boost::graph_traits<TriangleMesh>::face_descriptor>;
+
+  using Geom_traits = typename Kernel_traits<Point_d>::type;
+
+  using Construct_root_node_bbox = std::function<Bbox_d()>;
+  using Construct_root_node_contents = std::function<Node_data()>;
+  using Distribute_node_contents = std::function<void(Node_index, Tree&, const Point_d&)>;
+
+  /// @}
+
+  /// \name Operations
+  /// @{
+
+  auto construct_root_node_bbox_object() const {
+    return [&]() -> Bbox_d {
+
+      std::array<FT, Base::dimension> min = {0.0, 0}, max = {0.0, 0};
+      if (faces(m_pm).begin() != faces(m_pm).end()) {
+        bool first = true;
+        for (auto v: vertices(m_pm)) {
+          const Point_d& p_v = get(m_vpm, v);
+          for (int i = 0; i < 3; ++i) {
+            if (first || p_v[i] < min[i]) min[i] = p_v[i];
+            if (first || p_v[i] > max[i]) max[i] = p_v[i];
+          }
+          first=false;
+        }
+      }
+
+      return {std::apply(Self::construct_point_d_object(), min),
+              std::apply(Self::construct_point_d_object(), max)};
+    };
+  }
+
+  auto construct_root_node_contents_object() const {
+    return [&]() -> Node_data {
+      return {faces(m_pm).begin(), faces(m_pm).end()};
+    };
+  }
+
+  auto distribute_node_contents_object() const {
+    return [&](Node_index n, Tree& tree, const Point_d& /* center */) -> void {
+      Node_data& ndata = tree.data(n);
+      for (int i = 0; i < 8; ++i) {
+        Node_index child = tree.child(n, i);
+        Node_data& child_data = tree.data(child);
+        Bbox_d bbox = tree.bbox(child);
+        for (auto f : ndata) {
+          typename boost::graph_traits<TriangleMesh>::halfedge_descriptor
+            h = halfedge(f, m_pm);
+          typename Geom_traits::Triangle_3 t(get(m_vpm, source(h, m_pm)),
+            get(m_vpm, target(h, m_pm)),
+            get(m_vpm, target(next(h, m_pm), m_pm)));
+          if (do_intersect(t, bbox))
+            child_data.push_back(f);
+        }
+      }
+      };
+  }
+
+  /// @}
+
+  /// Recommended split predicate to pass to `Orthtree::refine()` function so
+  /// that the octree is refined until a node is either empty or has an extent
+  /// that would be smaller after split than the corresponding value provided to the constructor.
+  class Split_predicate_node_min_extent {
+
+    std::array<FT, 3> m_min_extent;
+
+  public:
+
+    /// constructor with `me` being the minimal value a node extent could be
+    /// (same value for all dimension).
+    Split_predicate_node_min_extent(const FT& me)
+      : m_min_extent({me, me, me}) {}
+
+    /// constructor with `me` being the minimal value a node extent could be
+    /// (one value per dimension).
+    Split_predicate_node_min_extent(const std::array<FT, 3>& me)
+      : m_min_extent(me) {}
+
+    /*!
+      \brief returns `true` if `ni` should be split, `false` otherwise.
+     */
+    template <typename NodeIndex, typename Tree>
+    bool operator()(NodeIndex ni, const Tree& tree) const {
+      if (tree.data(ni).empty()) return false;
+
+      Bbox_d bb = tree.bbox(ni);
+
+      for (int i = 0; i < 3; ++i)
+        if (((bb.max)()[i] - (bb.min)()[i]) < 2 * m_min_extent[i])
+          return false;
+      return true;
+    }
+  };
+
+
+private:
+
+  const TriangleMesh& m_pm;
+  VertexPointMap m_vpm;
+};
+
+} // end of CGAL namespace
+
+
+#endif // CGAL_ORTHREE_TRAITS_FACE_GRAPH_H

Orthtree/include/CGAL/Orthtree_traits_point.h

@@ -0,0 +1,201 @@
+// Copyright (c) 2023  INRIA (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)     : Jackson Campolattaro
+
+#ifndef ORTHTREE_TESTS_ORTHTREE_TRAITS_POINT_H
+#define ORTHTREE_TESTS_ORTHTREE_TRAITS_POINT_H
+
+#include <CGAL/license/Orthtree.h>
+
+#include <CGAL/Dimension.h>
+#include <CGAL/Bbox_2.h>
+#include <CGAL/Point_set_2.h>
+#include <CGAL/Orthtree/Cartesian_ranges.h>
+
+#include <CGAL/Orthtree_traits_base.h>
+
+namespace CGAL {
+
+template <typename Tree, typename PointMap>
+void reassign_points(
+  Tree& tree, PointMap& point_map,
+  typename Tree::Node_index n, const typename Tree::Point& center, typename Tree::Node_data points,
+  std::bitset<Tree::dimension> coord = {}, std::size_t dimension = 0
+) {
+
+  // Root case: reached the last dimension
+  if (dimension == Tree::dimension) {
+    tree.data(tree.child(n, coord.to_ulong())) = points;
+    return;
+  }
+
+  // Split the point collection around the center point on this dimension
+  auto split_point = std::partition(
+    points.begin(), points.end(),
+    [&](const auto& p) -> bool {
+      return (get(point_map, p)[int(dimension)] < center[int(dimension)]);
+    }
+  );
+
+  // Further subdivide the first side of the split
+  std::bitset<Tree::dimension> coord_left = coord;
+  coord_left[dimension] = false;
+  reassign_points(tree, point_map, n, center, {points.begin(), split_point}, coord_left, dimension + 1);
+
+  // Further subdivide the second side of the split
+  std::bitset<Tree::dimension> coord_right = coord;
+  coord_right[dimension] = true;
+  reassign_points(tree, point_map, n, center, {split_point, points.end()}, coord_right, dimension + 1);
+}
+
+/*!
+  \ingroup PkgOrthtreeTraits
+
+  Traits class for defining an orthtree of points using the class `CGAL::Orthtree`.
+
+  \tparam GeomTraits model of `Kernel`.
+  \tparam PointRange must be a model of `Range` whose value type is the key type of `PointMap` and whose iterator type is model of `RandomAccessIterator`
+  \tparam PointMap must be a model of `ReadablePropertyMap` whose value type is a point type from `GeomTraits` matching the current dimension
+  \tparam dimension the dimension of the ambient Euclidean space.
+
+  \warning The input point set is not copied. It is used directly
+  and is rearranged by the `Orthtree`. Altering the point range
+  after creating the orthtree will leave it in an invalid state.
+
+  \cgalModels{CollectionPartitioningOrthtreeTraits}
+  \sa `CGAL::Octree`
+  \sa `CGAL::Quadtree`
+  \sa `CGAL::Orthtree_traits_base<GeomTraits, dimension>`
+*/
+template <
+  typename GeomTraits,
+  typename PointRange,
+  typename PointMap = Identity_property_map<typename std::iterator_traits<typename PointRange::iterator>::value_type>,
+  bool hypercubic_nodes = false,
+  int dimension = Ambient_dimension<
+    typename std::iterator_traits<typename PointRange::iterator>::value_type,
+    GeomTraits
+  >::value
+>
+struct Orthtree_traits_point : public Orthtree_traits_base<GeomTraits, dimension> {
+public:
+  /// \name Types
+  /// @{
+  using Node_data = boost::iterator_range<typename PointRange::iterator>;
+  /// @}
+
+  using Base = Orthtree_traits_base<GeomTraits, dimension>;
+  using Self = Orthtree_traits_point<GeomTraits, PointRange, PointMap, hypercubic_nodes, dimension>;
+  using Tree = Orthtree<Self>;
+
+  using Node_index = typename Base::Node_index;
+  using Node_data_element = typename std::iterator_traits<typename PointRange::iterator>::value_type;
+
+  Orthtree_traits_point(
+    PointRange& points,
+    PointMap point_map = PointMap()
+  ) : m_points(points), m_point_map(point_map) {}
+
+  auto construct_root_node_bbox_object() const {
+    return [&]() -> typename Self::Bbox_d {
+
+      std::array<typename Self::FT, Self::dimension> bbox_min, bbox_max;
+      Orthtrees::internal::Cartesian_ranges<Self> cartesian_range;
+
+      // init bbox with first values found
+      {
+        const typename Self::Point_d& point = get(m_point_map, *(m_points.begin()));
+        std::size_t i = 0;
+        for (const typename Self::FT& x: cartesian_range(point)) {
+          bbox_min[i] = x;
+          bbox_max[i] = x;
+          ++i;
+        }
+      }
+      // Expand bbox to contain all points
+      for (const auto& p: m_points) {
+        const typename Self::Point_d& point = get(m_point_map, p);
+        std::size_t i = 0;
+        for (const typename Self::FT& x: cartesian_range(point)) {
+          bbox_min[i] = (std::min)(x, bbox_min[i]);
+          bbox_max[i] = (std::max)(x, bbox_max[i]);
+          ++i;
+        }
+      }
+
+#if !defined(_MSC_VER) || _MSC_VER > 1920
+      if constexpr (hypercubic_nodes) {
+#else
+      if (hypercubic_nodes) {
+#endif
+        std::array<typename Self::FT, Self::dimension> center;
+        typename Self::FT max_side = 0;
+        for (int i = 0; i < Self::dimension; i++) {
+          typename Self::FT side = bbox_max[i] - bbox_min[i];
+          max_side = (std::max<typename Self::FT>)(max_side, side);
+          center[i] = (bbox_min[i] + bbox_max[i]) * 0.5f;
+        }
+        max_side *= 0.5f;
+        for (int i = 0; i < Self::dimension; i++) {
+          bbox_min[i] = center[i] - max_side;
+          bbox_max[i] = center[i] + max_side;
+        }
+      }
+
+      return {std::apply(Self::construct_point_d_object(), bbox_min),
+              std::apply(Self::construct_point_d_object(), bbox_max)};
+    };
+  }
+
+  auto construct_root_node_contents_object() const {
+    return [&]() -> typename Self::Node_data {
+      return {m_points.begin(), m_points.end()};
+    };
+  }
+
+  auto distribute_node_contents_object() const {
+    return [&](Node_index n, Tree& tree, const typename Self::Point_d& center) {
+      CGAL_precondition(!tree.is_leaf(n));
+      reassign_points(tree, m_point_map, n, center, tree.data(n));
+    };
+  }
+
+  auto construct_sphere_d_object() const {
+    return [](const typename Self::Point_d& center, const typename Self::FT& squared_radius) -> typename Self::Sphere_d {
+      return typename Self::Sphere_d(center, squared_radius);
+      };
+  }
+
+  auto construct_center_d_object() const {
+    return [](const typename Self::Sphere_d& sphere) -> typename Self::Point_d {
+      return sphere.center();
+      };
+  }
+
+  auto compute_squared_radius_d_object() const {
+    return [](const typename Self::Sphere_d& sphere) -> typename Self::FT {
+      return sphere.squared_radius();
+      };
+  }
+
+  auto squared_distance_of_element_object() const {
+    return [&](const Node_data_element& index, const typename Self::Point_d& point) -> typename Self::FT {
+      return CGAL::squared_distance(get(m_point_map, index), point);
+      };
+  }
+
+  PointRange& m_points;
+  PointMap m_point_map;
+};
+
+}
+
+
+#endif //ORTHTREE_TESTS_ORTHTREE_TRAITS_POINT_H

Orthtree/include/CGAL/Quadtree.h

@@ -15,34 +15,28 @@
 #include <CGAL/license/Orthtree.h>
 
 #include <CGAL/Orthtree.h>
-#include <CGAL/Orthtree_traits_2.h>
+#include <CGAL/Orthtree_traits_point.h>
 
 namespace CGAL {
 
 /*!
-  \ingroup PkgOrthtreeClasses
+  \ingroup PkgOrthtreeRef
 
-  \brief Alias that specializes the `Orthtree` class to a 2D quadtree.
-
-  These two types are exactly equivalent:
-  - `Quadtree<GeomTraits, PointRange, PointMap>`
-  - `Orthtree<Orthtree_traits_2<GeomTraits>, PointRange, PointMap>`.
-
-  \warning This is a not a real class but an alias, please refer to
-  the documentation of `Orthtree`.
+  \brief Alias that specializes the `Orthtree` class to a 2D quadtree storing 2D points.
 
   \tparam GeomTraits must be a model of `Kernel`
-  \tparam PointRange_ must be a model of range whose value type is the key type of `PointMap`
+  \tparam PointRange must be a model of `Range` whose value type is the key type of `PointMap`
   \tparam PointMap must be a model of `ReadablePropertyMap` whose value type is `GeomTraits::Point_2`
+  \tparam square_nodes Boolean to enforce square nodes
 */
 template <typename GeomTraits, typename PointRange,
           typename PointMap = Identity_property_map
-         <typename std::iterator_traits<typename PointRange::iterator>::value_type> >
-#ifdef DOXYGEN_RUNNING
-class Quadtree;
-#else
-using Quadtree = Orthtree<Orthtree_traits_2<GeomTraits>, PointRange, PointMap>;
-#endif
+         <typename std::iterator_traits<typename PointRange::iterator>::value_type>,
+          bool squared_nodes = false
+>
+
+using Quadtree = Orthtree<Orthtree_traits_point<GeomTraits, PointRange, PointMap, squared_nodes, 2>>;
+
 
 } // namespace CGAL
 

Orthtree/package_info/Orthtree/dependencies

@@ -1,6 +1,10 @@
 Algebraic_foundations
+Cartesian_kernel
+Circulator
 Distance_2
 Distance_3
+Filtered_kernel
+Hash_map
 Installation
 Intersections_2
 Intersections_3
@@ -9,7 +13,11 @@ Kernel_23
 Modular_arithmetic
 Number_types
 Orthtree
+Point_set_2
 Profiling_tools
 Property_map
 STL_Extension
+Spatial_sorting
 Stream_support
+TDS_2
+Triangulation_2

Orthtree/test/Orthtree/CMakeLists.txt

@@ -16,6 +16,7 @@ create_single_source_cgal_program("test_octree_traverse.cpp")
 create_single_source_cgal_program("test_octree_intersecting.cpp")
 create_single_source_cgal_program("test_octree_copy_move_constructors.cpp")
 create_single_source_cgal_program("test_octree_kernels.cpp")
+create_single_source_cgal_program("test_octree_custom_properties.cpp")
 
 create_single_source_cgal_program("test_node_index.cpp")
 create_single_source_cgal_program("test_node_adjacent.cpp")

Orthtree/test/Orthtree/test_node.cpp

@@ -1,55 +0,0 @@
-
-#include <CGAL/Orthtree/Node.h>
-#include <CGAL/Octree/IO.h>
-#include <iostream>
-#include <cassert>
-
-typedef CGAL::Orthtree::Node<std::vector<int>::iterator> Node;
-
-int main(void) {
-
-  // Build a new node
-  Node n = Node();
-
-  // Check that its values are correct
-  assert(n.is_root());
-  assert(n.is_leaf());
-  assert(!n.parent());
-  assert(n.depth() == 0);
-  assert(n.location()[0] == 0 && n.location()[1] == 0 && n.location()[2] == 0);
-
-  // Split the node
-  n.split();
-
-  // Check that it's children's values are also correct
-  for (std::size_t i = 0; i < 8; ++i) {
-
-    assert(!n[i].is_root());
-    assert(n[i].is_leaf());
-    assert(*n[i].parent() == n);
-    assert(n[i].depth() == 1);
-  }
-
-  // Check that the parent has updated
-  assert(n.is_root());
-  assert(!n.is_leaf());
-
-  // Split one of the children
-  n[1].split();
-
-  // Check each of that child's children
-  for (std::size_t i = 0; i < 8; ++i) {
-
-    assert(!n[1][i].is_root());
-    assert(n[1][i].is_leaf());
-    assert(*n[1][i].parent() == n[1]);
-    assert(*n[1][i].parent()->parent() == n);
-    assert(n[1][i].depth() == 2);
-  }
-
-  // Check that the child's values have updated
-  assert(!n[1].is_root());
-  assert(!n[1].is_leaf());
-
-  return 0;
-}

Orthtree/test/Orthtree/test_node_adjacent.cpp

@@ -3,16 +3,14 @@
 #include <CGAL/Octree.h>
 #include <CGAL/Orthtree/Traversals.h>
 
-#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Point_set_3.h>
+#include <CGAL/Simple_cartesian.h>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>
-Octree;
-typedef Octree::Node Node;
-typedef Octree::Traits Traits;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
+using Traits = Octree::Traits;
 
 int main(void) {
 
@@ -42,33 +40,43 @@ int main(void) {
   std::cout << octree << std::endl;
 
   // Root node should have no siblings
-  assert(octree.root().adjacent_node(0).is_null());
-  assert(octree.root().adjacent_node(1).is_null());
-  assert(octree.root().adjacent_node(2).is_null());
-  assert(octree.root().adjacent_node(3).is_null());
-  assert(octree.root().adjacent_node(4).is_null());
-  assert(octree.root().adjacent_node(5).is_null());
+  assert(!octree.adjacent_node(octree.root(), 0));
+  assert(!octree.adjacent_node(octree.root(), 1));
+  assert(!octree.adjacent_node(octree.root(), 2));
+  assert(!octree.adjacent_node(octree.root(), 3));
+  assert(!octree.adjacent_node(octree.root(), 4));
+  assert(!octree.adjacent_node(octree.root(), 5));
 
   // Left Top Front node should have siblings to the Right, Down, and Back
-  auto left_top_back = octree.root()[Traits::LEFT_TOP_BACK];
+  auto left_top_back = octree.node(Traits::LEFT_TOP_BACK);
 
-  assert(octree.root()[Traits::RIGHT_TOP_BACK] == left_top_back.adjacent_node(Traits::RIGHT));
-  assert(octree.root()[Traits::LEFT_BOTTOM_BACK] == left_top_back.adjacent_node(Traits::DOWN));
-  assert(octree.root()[Traits::LEFT_TOP_FRONT] == left_top_back.adjacent_node(Traits::FRONT));
-  assert(left_top_back.adjacent_node(Traits::LEFT).is_null());
-  assert(left_top_back.adjacent_node(Traits::UP).is_null());
-  assert(left_top_back.adjacent_node(Traits::BACK).is_null());
+  assert(octree.node(Traits::RIGHT_TOP_BACK) ==
+         *octree.adjacent_node(left_top_back, Traits::RIGHT));
+  assert(octree.node(Traits::LEFT_BOTTOM_BACK) ==
+         *octree.adjacent_node(left_top_back, Traits::DOWN));
+  assert(octree.node(Traits::LEFT_TOP_FRONT) ==
+         *octree.adjacent_node(left_top_back, Traits::FRONT));
+  assert(!octree.adjacent_node(left_top_back, Traits::LEFT));
+  assert(!octree.adjacent_node(left_top_back, Traits::UP));
+  assert(!octree.adjacent_node(left_top_back, Traits::BACK));
 
-  std::cout << octree.root()[Traits::LEFT_BOTTOM_BACK] << std::endl;
+  auto right_top_back_of_left_bottom_back = octree.node(Traits::LEFT_BOTTOM_BACK, Traits::RIGHT_TOP_BACK);
 
-  auto right_top_back_of_left_bottom_back = octree.root()[Traits::LEFT_BOTTOM_BACK][Traits::RIGHT_TOP_BACK];
-  assert(octree.root()[Traits::LEFT_BOTTOM_BACK][Traits::LEFT_TOP_BACK] == right_top_back_of_left_bottom_back.adjacent_node(Traits::LEFT));
-  assert(octree.root()[Traits::RIGHT_BOTTOM_BACK] == right_top_back_of_left_bottom_back.adjacent_node(Traits::RIGHT));
-  assert(!right_top_back_of_left_bottom_back.adjacent_node(Traits::RIGHT).is_null());
-  assert(!right_top_back_of_left_bottom_back.adjacent_node(Traits::UP).is_null());
-  assert(!right_top_back_of_left_bottom_back.adjacent_node(Traits::DOWN).is_null());
-  assert(right_top_back_of_left_bottom_back.adjacent_node(Traits::BACK).is_null());
-  assert(!right_top_back_of_left_bottom_back.adjacent_node(Traits::FRONT).is_null());
+  assert(
+    octree.node(Traits::LEFT_BOTTOM_BACK, Traits::LEFT_TOP_BACK) ==
+    octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::LEFT)
+  );
+  assert(
+    octree.node(Traits::RIGHT_BOTTOM_BACK) ==
+    octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::RIGHT)
+  );
+  assert(octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::RIGHT).has_value());
+  assert(octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::UP).has_value());
+  assert(octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::DOWN).has_value());
+  assert(octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::FRONT).has_value());
+
+  // A node at the back of the tree should have no neighbor to its back
+  assert(!octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::BACK));
 
   return 0;
 }

Orthtree/test/Orthtree/test_node_index.cpp

@@ -3,14 +3,13 @@
 #include <CGAL/Octree.h>
 #include <CGAL/Orthtree/Traversals.h>
 
-#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Point_set_3.h>
+#include <CGAL/Simple_cartesian.h>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>
-Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
 
 int main(void) {
 
@@ -37,12 +36,11 @@ int main(void) {
   Octree octree(points, points.point_map());
   octree.refine(10, 1);
 
-  std::cout << "root: " << octree.root().local_coordinates() << std::endl;
-  std::cout << "first child: " << octree.root()[0].local_coordinates() << std::endl;
-  std::cout << "fifth child: " << octree.root()[4].local_coordinates() << std::endl;
-  std::cout << "fifth child of first child: " << octree.root()[0][4].local_coordinates() << std::endl;
-
-  // TODO
+  std::cout << "root: " << octree.local_coordinates(octree.root()) << std::endl;
+  std::cout << "first child: " << octree.local_coordinates(octree.child(octree.root(), 0)) << std::endl;
+  std::cout << "fifth child: " << octree.local_coordinates(octree.child(octree.root(), 4)) << std::endl;
+  std::cout << "fifth child of first child: "
+            << octree.local_coordinates(octree.child(octree.child(octree.root(), 0), 4)) << std::endl;
 
   return 0;
 }

Orthtree/test/Orthtree/test_octree_bbox.cpp

@@ -1,17 +1,14 @@
-
-#include <iostream>
 #include <CGAL/Octree.h>
-#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Point_set_3.h>
-
+#include <CGAL/Simple_cartesian.h>
+#include <iostream>
 #include <cassert>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef Kernel::FT FT;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>
-Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using FT = Kernel::FT;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
 
 void test_1_node() {
 
@@ -23,8 +20,10 @@ void test_1_node() {
   Octree octree(points, points.point_map());
   octree.refine(10, 1);
 
+  Octree::Bbox expected_bbox{-1, -1, -1, -1, -1, -1};
+
   // Compare the top (only) node
-  assert(octree.bbox(octree.root()) == CGAL::Bbox_3(-1, -1, -1, -1, -1, -1));
+  assert(octree.bbox(octree.root()) == Octree::Bbox(-1, -1, -1, -1, -1, -1));
 }
 
 void test_9_nodes() {
@@ -35,21 +34,21 @@ void test_9_nodes() {
   points.insert({1, 1, 1});
 
   // Create the octree
-  Octree octree(points, points.point_map(), 1.1);
+  Octree octree(points, points.point_map());
   octree.refine(10, 1);
 
   // Compare the top node
-  assert(octree.bbox(octree.root()) == CGAL::Bbox_3(-1.1, -1.1, -1.1, 1.1, 1.1, 1.1));
+  assert(octree.bbox(octree.root()) == Octree::Bbox(-1, -1, -1, 1, 1, 1));
 
   // Compare the child nodes
-  assert(octree.bbox(octree.root()[0]) == CGAL::Bbox_3(-1.1, -1.1, -1.1, 0, 0, 0));
-  assert(octree.bbox(octree.root()[1]) == CGAL::Bbox_3(0, -1.1, -1.1, 1.1, 0, 0));
-  assert(octree.bbox(octree.root()[2]) == CGAL::Bbox_3(-1.1, 0, -1.1, 0, 1.1, 0));
-  assert(octree.bbox(octree.root()[3]) == CGAL::Bbox_3(0, 0, -1.1, 1.1, 1.1, 0));
-  assert(octree.bbox(octree.root()[4]) == CGAL::Bbox_3(-1.1, -1.1, 0, 0, 0, 1.1));
-  assert(octree.bbox(octree.root()[5]) == CGAL::Bbox_3(0, -1.1, 0, 1.1, 0, 1.1));
-  assert(octree.bbox(octree.root()[6]) == CGAL::Bbox_3(-1.1, 0, 0, 0, 1.1, 1.1));
-  assert(octree.bbox(octree.root()[7]) == CGAL::Bbox_3(0, 0, 0, 1.1, 1.1, 1.1));
+  assert(octree.bbox(octree.node(0)) == Octree::Bbox(-1, -1, -1, 0, 0, 0));
+  assert(octree.bbox(octree.node(1)) == Octree::Bbox(0, -1, -1, 1, 0, 0));
+  assert(octree.bbox(octree.node(2)) == Octree::Bbox(-1, 0, -1, 0, 1, 0));
+  assert(octree.bbox(octree.node(3)) == Octree::Bbox(0, 0, -1, 1, 1, 0));
+  assert(octree.bbox(octree.node(4)) == Octree::Bbox(-1, -1, 0, 0, 0, 1));
+  assert(octree.bbox(octree.node(5)) == Octree::Bbox(0, -1, 0, 1, 0, 1));
+  assert(octree.bbox(octree.node(6)) == Octree::Bbox(-1, 0, 0, 0, 1, 1));
+  assert(octree.bbox(octree.node(7)) == Octree::Bbox(0, 0, 0, 1, 1, 1));
 }
 
 void test_25_nodes() {
@@ -62,41 +61,70 @@ void test_25_nodes() {
   points.insert({1, 0.5, 1});
 
   // Create the octree
-  Octree octree(points, points.point_map(), 1.5);
+  Octree octree(points, points.point_map());
   octree.refine(10, 1);
 
   // Compare the top node
-  assert(octree.bbox(octree.root()) == CGAL::Bbox_3(-1.5, -1.5, -1.5, 1.5, 1.5, 1.5));
+  assert(octree.bbox(octree.root()) == Octree::Bbox(-1, -1, -1, 1, 1, 1));
 
   // Compare the child nodes
-  assert(octree.bbox(octree.root()[0]) == CGAL::Bbox_3(-1.5, -1.5, -1.5, 0, 0, 0));
-  assert(octree.bbox(octree.root()[1]) == CGAL::Bbox_3(0, -1.5, -1.5, 1.5, 0, 0));
-  assert(octree.bbox(octree.root()[2]) == CGAL::Bbox_3(-1.5, 0, -1.5, 0, 1.5, 0));
-  assert(octree.bbox(octree.root()[3]) == CGAL::Bbox_3(0, 0, -1.5, 1.5, 1.5, 0));
-  assert(octree.bbox(octree.root()[4]) == CGAL::Bbox_3(-1.5, -1.5, 0, 0, 0, 1.5));
-  assert(octree.bbox(octree.root()[5]) == CGAL::Bbox_3(0, -1.5, 0, 1.5, 0, 1.5));
-  assert(octree.bbox(octree.root()[6]) == CGAL::Bbox_3(-1.5, 0, 0, 0, 1.5, 1.5));
-  assert(octree.bbox(octree.root()[7]) == CGAL::Bbox_3(0, 0, 0, 1.5, 1.5, 1.5));
+  assert(octree.bbox(octree.node(0)) == Octree::Bbox(-1, -1, -1, 0, 0, 0));
+  assert(octree.bbox(octree.node(1)) == Octree::Bbox(0, -1, -1, 1, 0, 0));
+  assert(octree.bbox(octree.node(2)) == Octree::Bbox(-1, 0, -1, 0, 1, 0));
+  assert(octree.bbox(octree.node(3)) == Octree::Bbox(0, 0, -1, 1, 1, 0));
+  assert(octree.bbox(octree.node(4)) == Octree::Bbox(-1, -1, 0, 0, 0, 1));
+  assert(octree.bbox(octree.node(5)) == Octree::Bbox(0, -1, 0, 1, 0, 1));
+  assert(octree.bbox(octree.node(6)) == Octree::Bbox(-1, 0, 0, 0, 1, 1));
+  assert(octree.bbox(octree.node(7)) == Octree::Bbox(0, 0, 0, 1, 1, 1));
 
   // Compare children of the first child
-  assert(octree.bbox(octree.root()[0][0]) == CGAL::Bbox_3(-1.5, -1.5, -1.5, -0.75, -0.75, -0.75));
-  assert(octree.bbox(octree.root()[0][1]) == CGAL::Bbox_3(-0.75, -1.5, -1.5, 0, -0.75, -0.75));
-  assert(octree.bbox(octree.root()[0][2]) == CGAL::Bbox_3(-1.5, -0.75, -1.5, -0.75, 0, -0.75));
-  assert(octree.bbox(octree.root()[0][3]) == CGAL::Bbox_3(-0.75, -0.75, -1.5, 0, 0, -0.75));
-  assert(octree.bbox(octree.root()[0][4]) == CGAL::Bbox_3(-1.5, -1.5, -0.75, -0.75, -0.75, 0));
-  assert(octree.bbox(octree.root()[0][5]) == CGAL::Bbox_3(-0.75, -1.5, -0.75, 0, -0.75, 0));
-  assert(octree.bbox(octree.root()[0][6]) == CGAL::Bbox_3(-1.5, -0.75, -0.75, -0.75, 0, 0));
-  assert(octree.bbox(octree.root()[0][7]) == CGAL::Bbox_3(-0.75, -0.75, -0.75, 0, 0, 0));
+  assert(octree.bbox(octree.node(0, 0)) ==
+         Octree::Bbox(-1, -1, -1, -0.5, -0.5, -0.5));
+  assert(octree.bbox(octree.node(0, 1)) ==
+         Octree::Bbox(-0.5, -1, -1, 0, -0.5, -0.5));
+  assert(octree.bbox(octree.node(0, 2)) ==
+         Octree::Bbox(-1, -0.5, -1, -0.5, 0, -0.5));
+  assert(octree.bbox(octree.node(0, 3)) ==
+         Octree::Bbox(-0.5, -0.5, -1, 0, 0, -0.5));
+  assert(octree.bbox(octree.node(0, 4)) ==
+         Octree::Bbox(-1, -1, -0.5, -0.5, -0.5, 0));
+  assert(octree.bbox(octree.node(0, 5)) ==
+         Octree::Bbox(-0.5, -1, -0.5, 0, -0.5, 0));
+  assert(octree.bbox(octree.node(0, 6)) ==
+         Octree::Bbox(-1, -0.5, -0.5, -0.5, 0, 0));
+  assert(octree.bbox(octree.node(0, 7)) ==
+         Octree::Bbox(-0.5, -0.5, -0.5, 0, 0, 0));
 
   // Compare children of the last child
-  assert(octree.bbox(octree.root()[7][0]) == CGAL::Bbox_3(0, 0, 0, 0.75, 0.75, 0.75));
-  assert(octree.bbox(octree.root()[7][1]) == CGAL::Bbox_3(0.75, 0, 0, 1.5, 0.75, 0.75));
-  assert(octree.bbox(octree.root()[7][2]) == CGAL::Bbox_3(0, 0.75, 0, 0.75, 1.5, 0.75));
-  assert(octree.bbox(octree.root()[7][3]) == CGAL::Bbox_3(0.75, 0.75, 0, 1.5, 1.5, 0.75));
-  assert(octree.bbox(octree.root()[7][4]) == CGAL::Bbox_3(0, 0, 0.75, 0.75, 0.75, 1.5));
-  assert(octree.bbox(octree.root()[7][5]) == CGAL::Bbox_3(0.75, 0, 0.75, 1.5, 0.75, 1.5));
-  assert(octree.bbox(octree.root()[7][6]) == CGAL::Bbox_3(0, 0.75, 0.75, 0.75, 1.5, 1.5));
-  assert(octree.bbox(octree.root()[7][7]) == CGAL::Bbox_3(0.75, 0.75, 0.75, 1.5, 1.5, 1.5));
+  assert(octree.bbox(octree.node(7, 0)) ==
+         Octree::Bbox(0, 0, 0, 0.5, 0.5, 0.5));
+  assert(octree.bbox(octree.node(7, 1)) ==
+         Octree::Bbox(0.5, 0, 0, 1, 0.5, 0.5));
+  assert(octree.bbox(octree.node(7, 2)) ==
+         Octree::Bbox(0, 0.5, 0, 0.5, 1, 0.5));
+  assert(octree.bbox(octree.node(7, 3)) ==
+         Octree::Bbox(0.5, 0.5, 0, 1, 1, 0.5));
+  assert(octree.bbox(octree.node(7, 4)) ==
+         Octree::Bbox(0, 0, 0.5, 0.5, 0.5, 1));
+  assert(octree.bbox(octree.node(7, 5)) ==
+         Octree::Bbox(0.5, 0, 0.5, 1, 0.5, 1));
+  assert(octree.bbox(octree.node(7, 6)) ==
+         Octree::Bbox(0, 0.5, 0.5, 0.5, 1, 1));
+  assert(octree.bbox(octree.node(7, 7)) ==
+         Octree::Bbox(0.5, 0.5, 0.5, 1, 1, 1));
+
+  // All child nodes should share a vertex
+  auto center_of_last_child = octree.bbox(octree.node(7, 7)).vertex(0);
+  assert(octree.bbox(octree.node(7, 0)).vertex(7) == center_of_last_child);
+  assert(octree.bbox(octree.node(7, 1)).vertex(4) == center_of_last_child);
+  assert(octree.bbox(octree.node(7, 2)).vertex(6) == center_of_last_child);
+  assert(octree.bbox(octree.node(7, 3)).vertex(5) == center_of_last_child);
+  assert(octree.bbox(octree.node(7, 4)).vertex(2) == center_of_last_child);
+  assert(octree.bbox(octree.node(7, 5)).vertex(3) == center_of_last_child);
+  assert(octree.bbox(octree.node(7, 6)).vertex(1) == center_of_last_child);
+
+  // Nodes of different sizes should also share vertices
+  assert(octree.bbox(octree.node(7, 0)).vertex(0) == octree.bbox(octree.node(0, 7)).vertex(7));
 }
 
 int main(void) {

Orthtree/test/Orthtree/test_octree_copy_move_constructors.cpp

@@ -1,21 +1,54 @@
 
 #define CGAL_TRACE_STREAM std::cerr
 
-#include <iostream>
 #include <CGAL/Octree.h>
-#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Orthtree_traits.h>
+#include <CGAL/Orthtree/Split_predicates.h>
+
 #include <CGAL/Point_set_3.h>
-
-#include <cassert>
 #include <CGAL/point_generators_3.h>
+#include <CGAL/Simple_cartesian.h>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef Kernel::FT FT;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map> Octree;
+#include <iostream>
+#include <cassert>
 
-int main(void)
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using FT = Kernel::FT;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
+using Octree_without_data = CGAL::Orthtree<CGAL::Orthtree_traits<Kernel, 3>>;
+
+template<typename Tree>
+int test(Tree &tree)
+{
+  assert (tree.is_leaf(tree.root())); // tree is not refined yet
+
+  Tree copy1 (tree);
+  assert (copy1.is_leaf(copy1.root())); // copy1 is thus not refined either
+  assert (tree == copy1); // tree should be equal to copy1
+
+  tree.refine(CGAL::Orthtrees::Maximum_depth(5));
+  assert (!tree.is_leaf(tree.root())); // tree is now refined
+  assert (copy1.is_leaf(copy1.root())); // copy1 should be unaffected and still unrefined
+  assert (tree != copy1); // tree should be different from copy1
+
+  Tree copy2 (tree);
+  assert (!copy2.is_leaf(copy2.root())); // copy2 should be refined
+  assert (tree == copy2); // tree should be equal to copy2
+
+  Tree move (std::move(tree));
+  assert (!move.is_leaf(move.root())); // move should be refined
+
+  assert (tree.is_leaf(tree.root())); // tree should be back to init state (unrefined)
+  assert (copy1.is_leaf(copy1.root())); // copy1 still unaffected and still unrefined
+  assert (!copy2.is_leaf(copy2.root())); // copy2 unaffected by move and still refined
+  assert (move == copy2); // move should be equal to copy2
+
+  return EXIT_SUCCESS;
+}
+
+int main()
 {
   std::size_t nb_pts = 100;
   Point_set points;
@@ -24,28 +57,9 @@ int main(void)
   for (std::size_t i = 0; i < nb_pts; ++i)
     points.insert(*(generator++));
 
-  Octree base (points, points.point_map());
-  assert (base.root().is_leaf()); // base is not refined yet
+  Octree base(points, points.point_map());
+  test(base);
 
-  Octree copy1 (base);
-  assert (copy1.root().is_leaf()); // copy1 is thus not refined either
-  assert (base == copy1); // base should be equal to copy1
-
-  base.refine();
-  assert (!base.root().is_leaf()); // base is now refined
-  assert (copy1.root().is_leaf()); // copy1 should be unaffected and still unrefined
-  assert (base != copy1); // base should be different from copy1
-
-  Octree copy2 (base);
-  assert (!copy2.root().is_leaf()); // copy2 should be refined
-  assert (base == copy2); // base should be equal to copy2
-
-  Octree move (std::move(base));
-  assert (!move.root().is_leaf()); // move should be refined
-  assert (base.root().is_leaf()); // base should be back to init state (unrefined)
-  assert (copy1.root().is_leaf()); // copy1 still unaffected and still unrefined
-  assert (!copy2.root().is_leaf()); // copy2 unaffected by move and still refined
-  assert (move == copy2); // move should be equal to copy2
-
-  return EXIT_SUCCESS;
+  Octree_without_data base2({});
+  test(base2);
 }

Orthtree/test/Orthtree/test_octree_custom_properties.cpp

@@ -0,0 +1,72 @@
+#define CGAL_TRACE_STREAM std::cerr
+
+#include <CGAL/Octree.h>
+
+#include <CGAL/Point_set_3.h>
+#include <CGAL/point_generators_3.h>
+#include <CGAL/Simple_cartesian.h>
+
+#include <iostream>
+#include <cassert>
+
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using FT = Kernel::FT;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
+
+int main(void) {
+
+  std::size_t nb_pts = 100;
+  Point_set points;
+  CGAL::Random_points_in_cube_3<Point> generator;
+  points.reserve(nb_pts);
+  for (std::size_t i = 0; i < nb_pts; ++i)
+    points.insert(*(generator++));
+
+  Octree tree(points, points.point_map());
+
+  // Testing built in node properties
+  typename Octree::Property_map<typename Octree::Node_data> data_prop = *tree.property<typename Octree::Node_data>("contents");
+  CGAL_USE(data_prop);
+
+  // list of properties
+  assert(tree.properties().size() == 5);
+  auto prop1 = tree.add_property("test", int(5));
+  assert(prop1.second);
+  // One property added
+  assert(tree.properties().size() == 6);
+  // Default value should be respected
+  assert(prop1.first[tree.root()] == 5);
+  // Changes to individual nodes should be respected
+  prop1.first[tree.root()] = 0;
+  assert(prop1.first[tree.root()] == 0);
+
+  auto prop2 = tree.add_property("test", int(0));
+  assert(!prop2.second);
+  assert(tree.properties().size() == 6);
+
+  auto prop3 = tree.add_property("test2", std::string());
+  assert(prop3.second);
+  assert(tree.properties().size() == 7);
+
+  auto prop4 = tree.property<int>("test");
+  assert(prop4.has_value());
+
+  auto prop5 = tree.property<std::string>("test");
+  assert(!prop5.has_value());
+
+  // Expanding the tree; new nodes should be assigned the default value
+  tree.refine(10, 1);
+  for (auto n : tree.traverse<CGAL::Orthtrees::Preorder_traversal<Octree>>()) {
+    // Everything but the root will have the default value
+    if (!tree.is_root(n)) assert(prop1.first[n] == 5);
+  }
+  // The root should have preserved its custom value
+  assert(prop1.first[tree.root()] == 0);
+
+  tree.remove_property(prop1.first);
+  assert(tree.properties().size() == 6);
+
+  return EXIT_SUCCESS;
+}

Orthtree/test/Orthtree/test_octree_equality.cpp

@@ -1,16 +1,15 @@
 
 #define CGAL_TRACE_STREAM std::cerr
 
-#include <iostream>
 #include <CGAL/Octree.h>
-#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Point_set_3.h>
+#include <CGAL/Simple_cartesian.h>
+#include <iostream>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>
-Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
 
 void test_identical_trees() {
 

Orthtree/test/Orthtree/test_octree_grade.cpp

@@ -1,36 +1,36 @@
 
 #define CGAL_TRACE_STREAM std::cerr
 
-#include <iostream>
 #include <CGAL/Octree.h>
 #include <CGAL/Orthtree/Traversals.h>
-#include <CGAL/Simple_cartesian.h>
+
 #include <CGAL/Point_set_3.h>
-
-#include <cassert>
 #include <CGAL/point_generators_3.h>
+#include <CGAL/Simple_cartesian.h>
+#include <iostream>
+#include <cassert>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef Kernel::FT FT;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map> Octree;
-typedef CGAL::Orthtrees::Leaves_traversal Leaves_traversal;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using FT = Kernel::FT;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
+using Leaves_traversal = CGAL::Orthtrees::Leaves_traversal<Octree>;
 
-std::size_t count_jumps(Octree &octree) {
+std::size_t count_jumps(Octree& octree) {
 
   std::size_t jumps = 0;
 
-  for (auto &node : octree.traverse(Leaves_traversal())) {
+  for (auto node: octree.traverse<Leaves_traversal>()) {
 
     for (int direction = 0; direction < 6; ++direction) {
 
-      auto adjacent_node = node.adjacent_node(direction);
+      auto adjacent_node = octree.adjacent_node(node, direction);
 
-      if (adjacent_node.is_null())
+      if (!adjacent_node)
         continue;
 
-      if ((node.depth() - adjacent_node.depth()) > 1)
+      if ((octree.depth(node) - octree.depth(*adjacent_node)) > 1)
         jumps++;
     }
   }

Orthtree/test/Orthtree/test_octree_intersecting.cpp

@@ -1,17 +1,17 @@
 #define CGAL_TRACE_STREAM std::cerr
 
-#include <iostream>
 #include <CGAL/Octree.h>
 #include <CGAL/Orthtree/Traversals.h>
-#include <CGAL/Simple_cartesian.h>
 
-#include <cassert>
 #include <CGAL/point_generators_3.h>
+#include <CGAL/Simple_cartesian.h>
+#include <iostream>
+#include <cassert>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef std::vector<Point> Point_vector;
-typedef CGAL::Octree<Kernel, Point_vector> Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_vector = std::vector<Point>;
+using Octree = CGAL::Octree<Kernel, Point_vector>;
 
 int main(void) {
 
@@ -34,7 +34,7 @@ int main(void) {
   points.emplace_back(-0.9, -1, -1);
 
   // Create an octree from the vector
-  Octree octree(points);
+  Octree octree(Octree::Traits{points});
 
   // Build the octree
   octree.refine(10, 2);
@@ -46,7 +46,7 @@ int main(void) {
     auto query = Point{1, 1, 1};
 
     // Get a list of nodes intersected
-    std::vector<Octree::Node> nodes{};
+    std::vector<Octree::Node_index> nodes{};
     octree.intersected_nodes(query, std::back_inserter(nodes));
 
     // A point should only intersect one node
@@ -63,15 +63,15 @@ int main(void) {
     auto query = Kernel::Sphere_3(Point{1, 0.5, 1}, 1.0);
 
     // Get a list of nodes intersected
-    std::vector<Octree::Node> nodes{};
+    std::vector<Octree::Node_index> nodes{};
     octree.intersected_nodes(query, std::back_inserter(nodes));
 
     // Check the results
     assert(4 == nodes.size());
-    assert(octree[Octree::Traits::RIGHT_TOP_BACK] == nodes[0]);
-    assert(octree[Octree::Traits::RIGHT_BOTTOM_FRONT] == nodes[1]);
-    assert(octree[Octree::Traits::LEFT_TOP_FRONT] == nodes[2]);
-    assert(octree[Octree::Traits::RIGHT_TOP_FRONT] == nodes[3]);
+    assert(octree.node(Octree::Traits::RIGHT_TOP_BACK) == nodes[0]);
+    assert(octree.node(Octree::Traits::RIGHT_BOTTOM_FRONT) == nodes[1]);
+    assert(octree.node(Octree::Traits::LEFT_TOP_FRONT) == nodes[2]);
+    assert(octree.node(Octree::Traits::RIGHT_TOP_FRONT) == nodes[3]);
   }
 
   // Intersection with a ray
@@ -81,19 +81,22 @@ int main(void) {
     auto query = Kernel::Ray_3(Point{1, 1, 1}, Point{0, 0, 0});
 
     // Get a list of nodes intersected
-    std::vector<Octree::Node> nodes{};
+    std::vector<Octree::Node_index> nodes{};
     octree.intersected_nodes(query, std::back_inserter(nodes));
 
     // Check the results
     assert(8 == nodes.size());
-    assert(octree[Octree::Traits::LEFT_BOTTOM_BACK] == nodes[0]);
-    assert(octree[Octree::Traits::RIGHT_BOTTOM_BACK][Octree::Traits::LEFT_TOP_FRONT] == nodes[1]);
-    assert(octree[Octree::Traits::LEFT_TOP_BACK] == nodes[2]);
-    assert(octree[Octree::Traits::RIGHT_TOP_BACK] == nodes[3]);
-    assert(octree[Octree::Traits::LEFT_BOTTOM_FRONT] == nodes[4]);
-    assert(octree[Octree::Traits::RIGHT_BOTTOM_FRONT] == nodes[5]);
-    assert(octree[Octree::Traits::LEFT_TOP_FRONT] == nodes[6]);
-    assert(octree[Octree::Traits::RIGHT_TOP_FRONT] == nodes[7]);
+    assert(octree.node(Octree::Traits::LEFT_BOTTOM_BACK) == nodes[0]);
+    assert(
+      octree.node(Octree::Traits::RIGHT_BOTTOM_BACK, Octree::Traits::LEFT_TOP_FRONT)
+      == nodes[1]
+    );
+    assert(octree.node(Octree::Traits::LEFT_TOP_BACK) == nodes[2]);
+    assert(octree.node(Octree::Traits::RIGHT_TOP_BACK) == nodes[3]);
+    assert(octree.node(Octree::Traits::LEFT_BOTTOM_FRONT) == nodes[4]);
+    assert(octree.node(Octree::Traits::RIGHT_BOTTOM_FRONT) == nodes[5]);
+    assert(octree.node(Octree::Traits::LEFT_TOP_FRONT) == nodes[6]);
+    assert(octree.node(Octree::Traits::RIGHT_TOP_FRONT) == nodes[7]);
   }
 
   return EXIT_SUCCESS;

Orthtree/test/Orthtree/test_octree_kernels.cpp

@@ -1,16 +1,16 @@
 #include <CGAL/Octree.h>
+#include <CGAL/Point_set_3.h>
+#include <CGAL/point_generators_3.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Point_set_3.h>
-#include <CGAL/point_generators_3.h>
 
 template <typename Kernel>
 void test()
 {
-  typedef typename Kernel::Point_3 Point;
-  typedef CGAL::Point_set_3<Point> Point_set;
-  typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map> Octree;
+  using Point = typename Kernel::Point_3;
+  using Point_set = CGAL::Point_set_3<Point>;
+  using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
 
   Point_set points;
   CGAL::Random_points_in_cube_3<Point> generator;

Orthtree/test/Orthtree/test_octree_locate.cpp

@@ -1,19 +1,18 @@
 
 #define CGAL_TRACE_STREAM std::cerr
 
-#include <iostream>
 #include <CGAL/Octree.h>
-#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Point_set_3.h>
 
+#include <CGAL/Simple_cartesian.h>
+#include <iostream>
 #include <cassert>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef Kernel::FT FT;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>
-Octree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using FT = Kernel::FT;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
 
 void test_1_point() {
 
@@ -52,24 +51,24 @@ void test_8_points() {
   octree.refine(10, 1);
 
   // Existing points should end up in the same place
-  assert(octree.root()[0] == octree.locate({-1, -1, -1}));
-  assert(octree.root()[1] == octree.locate({1, -1, -1}));
-  assert(octree.root()[2] == octree.locate({-1, 1, -1}));
-  assert(octree.root()[3] == octree.locate({1, 1, -1}));
-  assert(octree.root()[4] == octree.locate({-1, -1, 1}));
-  assert(octree.root()[5] == octree.locate({1, -1, 1}));
-  assert(octree.root()[6] == octree.locate({-1, 1, 1}));
-  assert(octree.root()[7] == octree.locate({1, 1, 1}));
+  assert(octree.node(0) == octree.locate({-1, -1, -1}));
+  assert(octree.node(1) == octree.locate({1, -1, -1}));
+  assert(octree.node(2) == octree.locate({-1, 1, -1}));
+  assert(octree.node(3) == octree.locate({1, 1, -1}));
+  assert(octree.node(4) == octree.locate({-1, -1, 1}));
+  assert(octree.node(5) == octree.locate({1, -1, 1}));
+  assert(octree.node(6) == octree.locate({-1, 1, 1}));
+  assert(octree.node(7) == octree.locate({1, 1, 1}));
 
   // Points adjacent to the existing points should also end up in the same place
-  assert(octree.root()[0] == octree.locate({-1.1, -1.1, -1.1}));
-  assert(octree.root()[1] == octree.locate({1.1, -1.1, -1.1}));
-  assert(octree.root()[2] == octree.locate({-1.1, 1.1, -1.1}));
-  assert(octree.root()[3] == octree.locate({1.1, 1.1, -1.1}));
-  assert(octree.root()[4] == octree.locate({-1.1, -1.1, 1.1}));
-  assert(octree.root()[5] == octree.locate({1.1, -1.1, 1.1}));
-  assert(octree.root()[6] == octree.locate({-1.1, 1.1, 1.1}));
-  assert(octree.root()[7] == octree.locate({1.1, 1.1, 1.1}));
+  assert(octree.node(0) == octree.locate({-0.99, -0.99, -0.99}));
+  assert(octree.node(1) == octree.locate({0.99, -0.99, -0.99}));
+  assert(octree.node(2) == octree.locate({-0.99, 0.99, -0.99}));
+  assert(octree.node(3) == octree.locate({0.99, 0.99, -0.99}));
+  assert(octree.node(4) == octree.locate({-0.99, -0.99, 0.99}));
+  assert(octree.node(5) == octree.locate({0.99, -0.99, 0.99}));
+  assert(octree.node(6) == octree.locate({-0.99, 0.99, 0.99}));
+  assert(octree.node(7) == octree.locate({0.99, 0.99, 0.99}));
 
 }
 
@@ -93,24 +92,24 @@ void test_10_points() {
   octree.refine(10, 1);
 
   // Existing points should end up in the same place
-  assert(octree.root()[0] == octree.locate({-1, -1, -1}));
-  assert(octree.root()[1] == octree.locate({1, -1, -1}));
-  assert(octree.root()[2] == octree.locate({-1, 1, -1}));
-  assert(octree.root()[3][3][3] == octree.locate({1, 1, -1}));
-  assert(octree.root()[4][4][4] == octree.locate({-1, -1, 1}));
-  assert(octree.root()[5] == octree.locate({1, -1, 1}));
-  assert(octree.root()[6] == octree.locate({-1, 1, 1}));
-  assert(octree.root()[7] == octree.locate({1, 1, 1}));
+  assert(octree.node(0) == octree.locate({-1, -1, -1}));
+  assert(octree.node(1) == octree.locate({1, -1, -1}));
+  assert(octree.node(2) == octree.locate({-1, 1, -1}));
+  assert(octree.node(3, 3, 3, 3, 3) == octree.locate({1, 1, -1}));
+  assert(octree.node(4, 4, 4) == octree.locate({-1, -1, 1}));
+  assert(octree.node(5) == octree.locate({1, -1, 1}));
+  assert(octree.node(6) == octree.locate({-1, 1, 1}));
+  assert(octree.node(7) == octree.locate({1, 1, 1}));
 
   // Points adjacent to the existing points might end up in different places
-  assert(octree.root()[0] == octree.locate({-1.1, -1.1, -1.1}));
-  assert(octree.root()[1] == octree.locate({1.1, -1.1, -1.1}));
-  assert(octree.root()[2] == octree.locate({-1.1, 1.1, -1.1}));
-  assert(octree.root()[3][3][3] == octree.locate({1.1, 1.1, -1.1}));
-  assert(octree.root()[4][4][4] == octree.locate({-1.1, -1.1, 1.1}));
-  assert(octree.root()[5] == octree.locate({1.1, -1.1, 1.1}));
-  assert(octree.root()[6] == octree.locate({-1.1, 1.1, 1.1}));
-  assert(octree.root()[7] == octree.locate({1.1, 1.1, 1.1}));
+  assert(octree.node(0) == octree.locate({-0.99, -0.99, -0.99}));
+  assert(octree.node(1) == octree.locate({0.99, -0.99, -0.99}));
+  assert(octree.node(2) == octree.locate({-0.99, 0.99, -0.99}));
+  assert(octree.node(3, 3, 3, 3, 3) == octree.locate({0.99, 0.99, -0.99}));
+  assert(octree.node(4, 4, 4) == octree.locate({-0.99, -0.99, 0.99}));
+  assert(octree.node(5) == octree.locate({0.99, -0.99, 0.99}));
+  assert(octree.node(6) == octree.locate({-0.99, 0.99, 0.99}));
+  assert(octree.node(7) == octree.locate({0.99, 0.99, 0.99}));
 
 }
 

Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp

@@ -1,30 +1,28 @@
 
 #define CGAL_TRACE_STREAM std::cerr
 
-#include <iostream>
 #include <CGAL/Octree.h>
-#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Point_set_3.h>
 #include <CGAL/point_generators_3.h>
 #include <CGAL/squared_distance_3.h>
 #include <CGAL/Orthogonal_k_neighbor_search.h>
 #include <CGAL/Search_traits_3.h>
 
+#include <CGAL/Simple_cartesian.h>
+#include <iostream>
 #include <chrono>
 #include <cassert>
 
 using namespace std::chrono;
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef Kernel::FT FT;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>
-Octree;
-
-typedef CGAL::Search_traits_3<Kernel> Kd_tree_traits;
-typedef CGAL::Orthogonal_k_neighbor_search<Kd_tree_traits> Kd_tree_search;
-typedef Kd_tree_search::Tree Kd_tree;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using FT = Kernel::FT;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
+using Kd_tree_traits = CGAL::Search_traits_3<Kernel>;
+using Kd_tree_search = CGAL::Orthogonal_k_neighbor_search<Kd_tree_traits>;
+using Kd_tree = Kd_tree_search::Tree;
 
 
 void naive_vs_octree(std::size_t dataset_size) {
@@ -47,7 +45,7 @@ void naive_vs_octree(std::size_t dataset_size) {
   {
 
     FT distance_nearest = (std::numeric_limits<FT>::max)();
-    for (auto &p : points.points()) {
+    for (auto& p: points.points()) {
 
       FT distance_current = CGAL::squared_distance(p, random_point);
       if (distance_current < distance_nearest) {
@@ -72,10 +70,9 @@ void naive_vs_octree(std::size_t dataset_size) {
   octree.refine(10, 20);
   auto octree_start_time = high_resolution_clock::now();
   {
-    // TODO: Write a nearest-neighbor implementation and use it here
-    std::vector<Point> k_neighbors;
-    octree.nearest_neighbors(random_point, 1, std::back_inserter(k_neighbors));
-    octree_nearest = *k_neighbors.begin();
+    std::vector<Point_set::Index> k_neighbors;
+    octree.nearest_k_neighbors(random_point, 1, std::back_inserter(k_neighbors));
+    octree_nearest = get(points.point_map(), *k_neighbors.begin());
   }
   duration<float> octree_elapsed_time = high_resolution_clock::now() - octree_start_time;
 
@@ -109,9 +106,9 @@ void kdtree_vs_octree(std::size_t dataset_size, std::size_t K) {
   Kd_tree kd_tree(points.points().begin(), points.points().end());
   kd_tree.build();
   auto kd_tree_start_time = high_resolution_clock::now();
-  Kd_tree_search search(kd_tree, random_point, (unsigned int)(K));
+  Kd_tree_search search(kd_tree, random_point, (unsigned int) (K));
   duration<float> kd_tree_elapsed_time = high_resolution_clock::now() - kd_tree_start_time;
-  for (auto p : search)
+  for (auto p: search)
     kd_tree_nearest_neighbors.push_back(p.first);
 
   std::cout << "Kd_tree --> "
@@ -120,11 +117,11 @@ void kdtree_vs_octree(std::size_t dataset_size, std::size_t K) {
             << std::endl;
 
   // Do the same using the octree
-  std::vector<Point> octree_nearest_neighbors;
+  std::vector<Point_set::Index> octree_nearest_neighbors;
   Octree octree(points, points.point_map());
   octree.refine(10, 20);
   auto octree_start_time = high_resolution_clock::now();
-  octree.nearest_neighbors(random_point, K, std::back_inserter(octree_nearest_neighbors));
+  octree.nearest_k_neighbors(random_point, K, std::back_inserter(octree_nearest_neighbors));
   duration<float> octree_elapsed_time = high_resolution_clock::now() - octree_start_time;
 
   std::cout << "Octree --> "
@@ -137,12 +134,13 @@ void kdtree_vs_octree(std::size_t dataset_size, std::size_t K) {
 
   // Check that they produce the same answer
   for (std::size_t j = 0; j < K; ++j)
-    assert(octree_nearest_neighbors[j] == kd_tree_nearest_neighbors[j]);
+    assert(get(points.point_map(), octree_nearest_neighbors[j]) == kd_tree_nearest_neighbors[j]);
 
 }
 
 int main(void) {
 
+  naive_vs_octree(21);
   naive_vs_octree(500);
   naive_vs_octree(1000);
   naive_vs_octree(10000);

Orthtree/test/Orthtree/test_octree_refine.cpp

@@ -1,18 +1,34 @@
 
 #define CGAL_TRACE_STREAM std::cerr
 
-#include <iostream>
 #include <CGAL/Octree.h>
-#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Point_set_3.h>
 
+#include <CGAL/Simple_cartesian.h>
+#include <iostream>
 #include <cassert>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map> Octree;
-typedef Octree::Node Node;
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
+
+class Split_nth_child_of_root {
+  std::size_t m_n;
+public:
+
+  Split_nth_child_of_root(std::size_t n) : m_n(n) {}
+
+  template <typename Node>
+  bool operator()(const Node& node) const {
+    return (node.depth() == 1 && node.local_coordinates().to_ulong() == m_n);
+  }
+
+  template<typename Node_index, typename Tree>
+  bool operator()(Node_index i, const Tree &tree) const {
+    return (tree.depth(i) == 1 && tree.local_coordinates(i).to_ulong() == m_n);
+  }
+};
 
 void test_1_point() {
 
@@ -24,13 +40,9 @@ void test_1_point() {
   Octree octree(points, points.point_map());
   octree.refine(10, 1);
 
-  // Check that the topology matches
-  Node single_node = CGAL::Orthtrees::Node_access::create_node(Node(), 0);
-  CGAL::Orthtrees::Node_access::points(single_node)
-    = CGAL::Orthtrees::Node_access::points(octree.root());
-  assert(Node::is_topology_equal(single_node, octree.root()));
+  // Check that the root node was never split
+  assert(octree.is_leaf(octree.root()));
   assert(0 == octree.depth());
-  CGAL::Orthtrees::Node_access::free(single_node);
 }
 
 void test_2_points() {
@@ -45,12 +57,10 @@ void test_2_points() {
   octree.refine(10, 1);
 
   // The octree should have been split once
-  Node other = CGAL::Orthtrees::Node_access::create_node(Node(), 0);
-  CGAL::Orthtrees::Node_access::split(other);
-  assert(Node::is_topology_equal(other, octree.root()));
+  Octree other(points, points.point_map());
+  other.split(other.root());
+  assert(Octree::is_topology_equal(other, octree));
   assert(1 == octree.depth());
-  CGAL::Orthtrees::Node_access::free(other);
-
 }
 
 void test_4_points() {
@@ -65,14 +75,18 @@ void test_4_points() {
   Octree octree(points, points.point_map());
   octree.refine(10, 1);
 
-  // The octree should have been split once on the first level, and twice on the second
-  Node other = CGAL::Orthtrees::Node_access::create_node(Node(), 0);
-  CGAL::Orthtrees::Node_access::split(other);
-  CGAL::Orthtrees::Node_access::split(other[3]);
-  CGAL::Orthtrees::Node_access::split(other[7]);
-  assert(Node::is_topology_equal(other, octree.root()));
+  Octree other(points, points.point_map());
+  other.split(other.root());
+  other.split(other.node(3));
+  other.split(other.node(7));
+  assert(Octree::is_topology_equal(other, octree));
   assert(2 == octree.depth());
-  CGAL::Orthtrees::Node_access::free(other);
+
+  // Applying another splitting criterion shouldn't reset the tree.
+  octree.refine(Split_nth_child_of_root(2));
+  other.split(other.node(2));
+  assert(Octree::is_topology_equal(other, octree));
+
 }
 
 int main(void) {

Orthtree/test/Orthtree/test_octree_traverse.cpp

@@ -2,17 +2,18 @@
 
 #include <CGAL/Octree.h>
 #include <CGAL/Orthtree/Traversals.h>
-
-#include <CGAL/Simple_cartesian.h>
 #include <CGAL/Point_set_3.h>
 
+#include <CGAL/Simple_cartesian.h>
 #include <cassert>
 
-typedef CGAL::Simple_cartesian<double> Kernel;
-typedef Kernel::Point_3 Point;
-typedef CGAL::Point_set_3<Point> Point_set;
-typedef CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map> Octree;
-typedef CGAL::Orthtrees::Preorder_traversal Preorder_traversal;
+
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point = Kernel::Point_3;
+using Point_set = CGAL::Point_set_3<Point>;
+using Octree = CGAL::Octree<Kernel, Point_set, typename Point_set::Point_map>;
+using Preorder_traversal = CGAL::Orthtrees::Preorder_traversal<Octree>;
+using Level_traversal = CGAL::Orthtrees::Level_traversal<Octree>;
 
 bool test_preorder_1_node() {
 
@@ -53,7 +54,31 @@ bool test_preorder_9_nodes() {
   assert(*iter == octree.root());
   for (int i = 0; i < 8; ++i) {
     iter++;
-    assert((*iter == octree.root()[i]));
+    assert(*iter == octree.node(i));
+  }
+
+  return true;
+}
+
+bool test_level_9_nodes() {
+
+  // Define the dataset
+  Point_set points;
+  points.insert({-1, -1, -1});
+  points.insert({1, -1, -1});
+
+  // Create the octree
+  Octree octree(points, points.point_map());
+  octree.refine(10, 1);
+
+  // Create the range
+  auto nodes = octree.traverse<Level_traversal>(static_cast<std::size_t>(1));
+
+  // Check each item in the range
+  auto iter = nodes.begin();
+  for (int i = 0; i < 8; ++i) {
+    assert(*iter == octree.node(i));
+    iter++;
   }
 
   return true;
@@ -71,6 +96,7 @@ bool test_preorder_25_nodes() {
   // Create the octree
   Octree octree(points, points.point_map());
   octree.refine(10, 1);
+  std::cout << octree << std::endl;
 
   // Create the range
   auto nodes = octree.traverse<Preorder_traversal>();
@@ -79,28 +105,28 @@ bool test_preorder_25_nodes() {
   auto iter = nodes.begin();
   assert(*iter == octree.root());
   iter++;
-  assert((*iter == octree.root()[0]));
+  assert(*iter == octree.node(0));
   iter++;
-  assert((*iter == octree.root()[1]));
+  assert(*iter == octree.node(1));
   iter++;
-  assert((*iter == octree.root()[2]));
+  assert((*iter == octree.node(2)));
   iter++;
-  assert((*iter == octree.root()[3]));
+  assert(*iter == octree.node(3));
   for (int i = 0; i < 8; ++i) {
     iter++;
-    assert((*iter == octree.root()[3][i]));
+    assert(*iter == octree.node(3, i));
   }
   iter++;
-  assert((*iter == octree.root()[4]));
+  assert((*iter == octree.node(4)));
   iter++;
-  assert((*iter == octree.root()[5]));
+  assert((*iter == octree.node(5)));
   iter++;
-  assert((*iter == octree.root()[6]));
+  assert((*iter == octree.node(6)));
   iter++;
-  assert((*iter == octree.root()[7]));
+  assert((*iter == octree.node(7)));
   for (int i = 0; i < 8; ++i) {
     iter++;
-    assert((*iter == octree.root()[7][i]));
+    assert(*iter == octree.node(7, i));
   }
 
   return true;
@@ -110,6 +136,7 @@ int main(void) {
 
   test_preorder_1_node();
   test_preorder_9_nodes();
+  test_level_9_nodes();
   test_preorder_25_nodes();
 
   return 0;

Property_map/include/CGAL/Property_container.h

@@ -0,0 +1,630 @@
+// Copyright (c) 2023 INRIA
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Jackson Campolattaro
+
+#ifndef CGAL_PROPERTY_CONTAINTER_H
+#define CGAL_PROPERTY_CONTAINTER_H
+
+#include <CGAL/assertions.h>
+
+#include <map>
+#include <optional>
+#include <iterator>
+
+#include <boost/property_map/property_map.hpp>
+
+#ifndef DOXYGEN_RUNNING
+
+namespace CGAL::Properties::Experimental {
+
+template <typename Index>
+class Property_array_base {
+public:
+
+  Property_array_base() = default;
+
+  Property_array_base(const Property_array_base<Index>& rhs) = delete;
+
+  virtual ~Property_array_base() = default;
+
+  // Declare virtual functions here, for things which need to be done within the Property container
+  // todo: these should mostly be private, and made available using friend
+
+  virtual std::shared_ptr<Property_array_base<Index>> empty_clone(const std::vector<bool>& active_indices) = 0;
+
+  virtual std::shared_ptr<Property_array_base<Index>> clone(const std::vector<bool>& active_indices) = 0;
+
+  virtual void copy(const Property_array_base<Index>& other) = 0;
+
+  // desactived as MSVC 2017 as an issue with that but it is not currently used.
+#if 0
+  virtual void move(Property_array_base<Index>&& other) = 0;
+#endif
+
+  virtual void append(const Property_array_base<Index>& other) = 0;
+
+  virtual void reserve(std::size_t n) = 0;
+
+  virtual void shrink_to_fit() = 0;
+
+  virtual void swap(Index a, Index b) = 0;
+
+  virtual void reset(Index i) = 0;
+
+  virtual const std::type_info& type() const = 0;
+
+  virtual void transfer_from(const Property_array_base<Index>& other_base, Index other_index, Index this_index) = 0;
+
+};
+
+/*!
+ * \brief Indexed storage for arbitrary types
+ *
+ * todo: make this effectively private, prioritize the use of Property_array_handle
+ *
+ * @tparam T
+ */
+template <typename Index, typename T>
+class Property_array : public Property_array_base<Index> {
+
+  std::vector<T> m_data;
+  const std::vector<bool>& m_active_indices;
+  T m_default_value;
+
+public:
+
+  using value_type = T;
+  using reference = typename std::vector<T>::reference;
+  using const_reference = typename std::vector<T>::const_reference;
+  using iterator = typename std::vector<T>::iterator;
+  using const_iterator = typename std::vector<T>::const_iterator;
+
+  Property_array(const std::vector<bool>& active_indices, const T& default_value) :
+    m_data(), m_active_indices(active_indices), m_default_value(default_value) {
+
+    m_data.reserve(active_indices.capacity());
+    m_data.resize(active_indices.size(), m_default_value);
+  }
+
+  virtual std::shared_ptr<Property_array_base<Index>> empty_clone(const std::vector<bool>& active_indices) override {
+    return std::make_shared<Property_array<Index, T>>(active_indices, m_default_value);
+  }
+
+  virtual std::shared_ptr<Property_array_base<Index>> clone(const std::vector<bool>& active_indices) override {
+    auto new_array = std::make_shared<Property_array<Index, T>>(active_indices, m_default_value);
+    new_array->m_data = m_data;
+    return new_array;
+  }
+
+  virtual void copy(const Property_array_base<Index>& other_base) override {
+    auto& other = dynamic_cast<const Property_array<Index, T>&>(other_base);
+    m_data = other.m_data;
+    CGAL_precondition(m_active_indices.size() == m_data.size());
+  }
+
+// deactived as MSVC 2017 as an issue with that but it is not currently used.
+#if 0
+  virtual void move(Property_array_base<Index>&& other_base) override {
+    auto&& other = static_cast<Property_array<Index, T>&&>(other_base);
+    m_data = std::move(other.m_data);
+    CGAL_precondition(m_active_indices.size() == m_data.size());
+  }
+#endif
+
+  virtual void append(const Property_array_base<Index>& other_base) override {
+    auto& other = dynamic_cast<const Property_array<Index, T>&>(other_base);
+    CGAL_precondition(m_data.size() + other.m_data.size() == m_active_indices.size());
+    m_data.insert(m_data.end(), other.m_data.begin(), other.m_data.end());
+  }
+
+  virtual void reserve(std::size_t n) override {
+    CGAL_precondition(m_active_indices.size() == n);
+    m_data.resize(n, m_default_value);
+  };
+
+  virtual void shrink_to_fit() override {
+    m_data.shrink_to_fit();
+  }
+
+  virtual void swap(Index a, Index b) override {
+    // todo: maybe cast to index, instead of casting index to size?
+    CGAL_precondition(std::size_t(a) < m_data.size() && std::size_t(b) < m_data.size());
+    std::iter_swap(m_data.begin() + a, m_data.begin() + b);
+  };
+
+  virtual void reset(Index i) override {
+    CGAL_precondition(std::size_t(i) < m_data.size());
+    m_data[std::size_t(i)] = m_default_value;
+  };
+
+  virtual const std::type_info& type() const override { return typeid(T); };
+
+  virtual void transfer_from(const Property_array_base<Index>& other_base,
+                             Index other_index, Index this_index) override {
+
+    CGAL_precondition(other_base.type() == type());
+    auto& other = dynamic_cast<const Property_array<Index, T>&>(other_base);
+    CGAL_precondition(std::size_t(other_index) < other.capacity() && std::size_t(this_index) < capacity());
+    m_data[this_index] = other.m_data[other_index];
+  }
+
+public:
+
+  // todo: there's not really a good reason to use these, maybe they should be removed
+
+  [[nodiscard]] std::size_t size() const { return std::count(m_active_indices.begin(), m_active_indices.end(), true); }
+
+  [[nodiscard]] std::size_t capacity() const { return m_data.size(); }
+
+  const_reference operator[](Index i) const {
+    CGAL_precondition(std::size_t(i) < m_data.size());
+    return m_data[std::size_t(i)];
+  }
+
+  reference operator[](Index i) {
+    CGAL_precondition(std::size_t(i) < m_data.size());
+    return m_data[std::size_t(i)];
+  }
+
+  iterator begin() { return m_data.begin(); }
+
+  iterator end() { return m_data.end(); }
+
+  const_iterator begin() const { return m_data.begin(); }
+
+  const_iterator end() const { return m_data.end(); }
+
+public:
+
+  bool operator==(const Property_array<Index, T>& other) const {
+    return &other == this;
+  }
+
+  bool operator!=(const Property_array<Index, T>& other) const { return !operator==(other); }
+
+};
+
+// todo: property maps/array handles should go in their own file
+
+// todo: add const/read-only handle
+template <typename Index, typename T>
+class Property_array_handle {
+
+  std::reference_wrapper<Property_array<Index, T>> m_array;
+
+public:
+
+  // Necessary for use as a boost::property_type
+  using key_type = Index;
+  using value_type = T;
+  using reference = typename std::vector<T>::reference;
+  using const_reference = typename std::vector<T>::const_reference;
+  using category = boost::lvalue_property_map_tag;
+
+  using iterator = typename std::vector<T>::iterator;
+  using const_iterator = typename std::vector<T>::const_iterator;
+
+  Property_array_handle(Property_array<Index, T>& array) : m_array(array) {}
+
+  [[nodiscard]] std::size_t size() const { return m_array.get().size(); }
+
+  [[nodiscard]] std::size_t capacity() const { return m_array.get().capacity(); }
+
+  Property_array<Index, T>& array() const { return m_array.get(); }
+
+  // todo: This might not be needed, if the other operator[] is made const
+  const_reference operator[](Index i) const { return m_array.get()[i]; }
+
+  reference operator[](Index i) { return m_array.get()[i]; }
+
+  // todo: maybe these can be const, in an lvalue property map?
+  iterator begin() { return m_array.get().begin(); }
+
+  iterator end() { return m_array.get().end(); }
+
+  const_iterator begin() const { return m_array.get().begin(); }
+
+  const_iterator end() const { return m_array.get().end(); }
+
+  bool operator==(const Property_array_handle<Index, T>& other) const { return other.m_array.get() == m_array.get(); }
+
+  bool operator!=(const Property_array_handle<Index, T>& other) const { return !operator==(other); }
+
+  inline friend reference get(Property_array_handle<Index, T> p, const Index& i) { return p[i]; }
+
+  inline friend void put(Property_array_handle<Index, T> p, const Index& i, const T& v) { p[i] = v; }
+
+};
+
+template <typename Index = std::size_t>
+class Property_container {
+
+  std::multimap<std::string, std::shared_ptr<Property_array_base<Index>>> m_properties;
+  std::vector<bool> m_active_indices{};
+
+public:
+
+  template<typename T>
+  using Array = Property_array<Index, T>;
+
+  Property_container() = default;
+
+  Property_container(const Property_container<Index>& other) {
+    m_active_indices = other.m_active_indices;
+
+    for (auto [name, array] : other.m_properties) {
+      // todo: this could probably be made faster using emplace_hint
+      m_properties.emplace(
+        name,
+        array->clone(m_active_indices)
+      );
+    }
+  }
+
+  Property_container(Property_container<Index>&& other) { *this = std::move(other); }
+
+  // This is not exactly an assignment as existing unique properties are kept.
+  Property_container<Index>& operator=(const Property_container<Index>& other) {
+    m_active_indices = other.m_active_indices;
+
+    for (auto [name, array] : other.m_properties) {
+      // search if property already exists
+      auto range = m_properties.equal_range(name);
+      auto it = range.first;
+      for (; it != range.second; it++) {
+        if (typeid(*array) == typeid((*it->second)))
+          break;
+      }
+
+      if (it != range.second)
+        it->second->copy(*array);
+      else
+        m_properties.emplace(name, array->clone(m_active_indices));
+    }
+
+    return *this;
+  }
+
+  // This is not exactly an assignment as existing unique properties are kept.
+  Property_container<Index>& operator=(Property_container<Index>&& other) {
+    m_active_indices = std::move(other.m_active_indices);
+
+    for (auto [name, array] : other.m_properties) {
+      // search if property already exists
+      auto range = m_properties.equal_range(name);
+      auto it = range.first;
+      for (; it != range.second; it++) {
+        if (typeid(*array) == typeid((*it->second)))
+          break;
+      }
+
+      if (it != range.second)
+        it->second->copy(std::move(*array));
+      else
+        m_properties.emplace(name, array->clone(m_active_indices));
+    }
+
+    // The moved-from property map should retain all of its properties, but contain 0 elements
+    other.reserve(0);
+    return *this;
+  }
+
+  template <typename T>
+  std::pair<std::reference_wrapper<Property_array<Index, T>>, bool>
+  get_or_add_property(const std::string& name, const T default_value = T()) {
+    auto range = m_properties.equal_range(name);
+    for (auto it = range.first; it != range.second; it++) {
+      Property_array<Index, T>* typed_array_ptr = dynamic_cast<Property_array<Index, T>*>(it->second.get());
+      if (typed_array_ptr != nullptr)
+        return { {*typed_array_ptr}, false };
+    }
+
+    auto it = m_properties.emplace(
+      name,
+      std::make_shared<Property_array<Index, T>>(
+        m_active_indices,
+        default_value
+      )
+    );
+
+    return {{*dynamic_cast<Property_array<Index, T>*>(it->second.get())}, true};
+  }
+
+  template <typename T>
+  Property_array<Index, T>& add_property(const std::string& name, const T default_value = T()) {
+    // todo: I'm not settled on the naming, but it's really convenient to have a function like this
+    auto [array, created] = get_or_add_property(name, default_value);
+    CGAL_precondition(created);
+    return array.get();
+  }
+
+/*
+  // todo: misleading name, maybe it could be add_same_properties?
+  void copy_properties(const Property_container<Index>& other) {
+    for (auto [name, other_array]: other.m_properties) {
+      // If this container doesn't have any property by this name, add it (with the same type as in other)
+      if (!property_exists(name))
+        m_property_arrays.emplace(name, other_array->empty_clone(m_active_indices));
+    }
+  }*/
+
+  template <typename T>
+  const Property_array<Index, T>& get_property(const std::string& name) const {
+    return *(get_property_if_exists<T>(name));
+  }
+
+  template <typename T>
+  Property_array<Index, T>& get_property(const std::string& name) {
+    return *(get_property_if_exists<T>(name));
+  }
+
+  template <typename T>
+  std::optional<std::reference_wrapper<Property_array<Index, T>>> get_property_if_exists(const std::string& name) const {
+    auto range = m_properties.equal_range(name);
+    for (auto it = range.first; it != range.second; it++) {
+      Property_array<Index, T>* typed_array_ptr = dynamic_cast<Property_array<Index, T>*>(it->second.get());
+      if (typed_array_ptr != nullptr)
+        return *typed_array_ptr;
+    }
+
+    return {};
+  }
+
+  template <typename T>
+  bool property_exists(const std::string& name) const {
+    auto range = m_properties.equal_range(name);
+
+    for (auto it = range.first; it != range.second; it++) {
+      Property_array<Index, T>* typed_array_ptr = dynamic_cast<Property_array<Index, T>*>(it->second.get());
+      if (typed_array_ptr != nullptr)
+        return true;
+    }
+
+    return false;
+  }
+
+  /*!
+   * Removes all properties with the name from the container.
+   *
+   * @param name
+   * @return number of removed properties.
+   */
+  std::size_t remove_properties(const std::string& name) { return m_properties.erase(name); }
+
+  template <typename T>
+  bool remove_property(const Property_array<Index, T>& arrayToRemove) {
+    const Property_array_base<Index>* ref = dynamic_cast<const Property_array_base<Index>*>(&arrayToRemove);
+    for (auto it = m_properties.begin(); it != m_properties.end(); it++) {
+      auto const& [name, array] = *it;
+      if (array.get() == ref) {
+        m_properties.erase(it);
+        return true;
+      }
+    }
+    return false;
+  }
+
+  void remove_all_properties_except(const std::vector<std::string>& preserved_names) {
+    // todo: if this is used often, it should take a parameter pack instead of a vector
+    // A fold expression could then be used in place of std::find for better performance
+    for (auto it = m_properties.begin(); it != m_properties.end();) {
+      auto const& [name, array] = *it;
+      if (std::find(preserved_names.begin(), preserved_names.end(), name) == preserved_names.end())
+        it = m_properties.erase(it);
+      else
+        it++;
+    }
+  }
+
+  std::vector<std::string> properties() const {
+    std::vector<std::string> property_names{};
+    for (auto const& [name, _]: m_properties)
+      property_names.emplace_back(name);
+    return property_names;
+  }
+
+  std::size_t num_properties() const { return m_properties.size(); }
+
+/* Deactivated as there may be several Property_maps with different types but the same name.
+  const std::type_info& property_type(const std::string& name) const {
+    if (auto it = m_property_arrays.find(name); it != m_property_arrays.end())
+      return it->second->type();
+    else
+      return typeid(void);
+  }*/
+
+public:
+
+  void reserve(std::size_t n) {
+    m_active_indices.resize(n);
+    for (auto [name, array]: m_properties)
+      array->reserve(n);
+  }
+
+  void resize(std::size_t n) {
+    reserve(n);
+    std::fill(m_active_indices.begin(), m_active_indices.end(), true);
+  }
+
+  [[nodiscard]] std::size_t size() const { return std::count(m_active_indices.begin(), m_active_indices.end(), true); }
+
+  [[nodiscard]] std::size_t capacity() const { return m_active_indices.size(); }
+
+  Index emplace_back() {
+
+    // Expand the storage and return the last element
+    reserve(capacity() + 1);
+    m_active_indices.back() = true;
+    auto first_new_index = Index(capacity() - 1);
+    reset(first_new_index);
+    return first_new_index;
+  }
+
+  Index emplace() {
+
+    // If there are empty slots, return the index of one of them and mark it as full
+    auto first_unused = std::find_if(m_active_indices.begin(), m_active_indices.end(), [](bool used) { return !used; });
+    if (first_unused != m_active_indices.end()) {
+      *first_unused = true;
+      auto index = Index(std::distance(m_active_indices.begin(), first_unused));
+      reset(index);
+      return index;
+    }
+
+    return emplace_back();
+  }
+
+  Index emplace_group_back(std::size_t n) {
+
+    // Expand the storage and return the start of the new region
+    reserve(capacity() + n);
+    for (auto it = m_active_indices.end() - n; it < m_active_indices.end(); ++it)
+      *it = true;
+    return Index(capacity() - n);
+  }
+
+  Index emplace_group(std::size_t n) {
+
+    auto search_start = m_active_indices.begin();
+    while (search_start != m_active_indices.end()) {
+
+      // Find the first unused cell
+      auto unused_begin = std::find_if(
+        search_start, m_active_indices.end(),
+        [](bool used) { return !used; }
+      );
+
+      auto unused_end = unused_begin;
+
+      // Determine if the group fits
+      if (std::distance(unused_begin, m_active_indices.end()) >= static_cast<typename std::iterator_traits<std::vector<bool>::iterator>::difference_type>(n))
+        unused_end = std::find_if(
+          unused_begin, (std::min)(unused_begin + n, m_active_indices.end()),
+          [](bool used) { return used; }
+      );
+
+      // If the discovered range was large enough
+      if (std::distance(unused_begin, unused_end) >= static_cast<typename std::iterator_traits<std::vector<bool>::iterator>::difference_type>(n)) {
+
+        // Mark the indices as used, and reset the properties of each of them
+        // todo: it would be better to provide a function to set a range
+        for (auto it = unused_begin; it < unused_end; ++it) {
+          *it = true;
+          reset(Index(std::distance(m_active_indices.begin(), it)));
+        }
+
+        // Return the first index of the range
+        return Index(std::distance(m_active_indices.begin(), unused_begin));
+      }
+
+      // If we didn't find a large enough region, continue our search after the end
+      search_start = unused_end;
+    }
+
+    // If no empty regions were found, expand the storage
+    return emplace_group_back(n);
+  }
+
+  void swap(Index a, Index b) {
+    for (auto [name, array]: m_properties)
+      array->swap(a, b);
+  }
+
+  void reset(Index i) {
+    for (auto [name, array]: m_properties)
+      array->reset(i);
+  }
+
+  void erase(Index i) {
+    m_active_indices[i] = false;
+    for (auto [name, array]: m_properties)
+      array->reset(i);
+  }
+
+  bool is_erased(Index i) const {
+    return !m_active_indices[i];
+  }
+
+  // todo: I'd prefer to eliminate this, if possible
+  void mark_active(Index i) {
+    return m_active_indices[i] = true;
+  }
+
+  void mark_inactive(Index i) {
+    return m_active_indices[i] = false;
+  }
+
+  std::vector<Index> active_list() const {
+    std::vector<Index> indices;
+    for (std::size_t i = 0; i < m_active_indices.size(); ++i)
+      if (m_active_indices[i]) indices.emplace_back(i);
+    return indices;
+  }
+
+  std::vector<Index> inactive_list() const {
+    std::vector<Index> indices;
+    for (std::size_t i = 0; i < m_active_indices.size(); ++i)
+      if (!m_active_indices[i]) indices.emplace_back(i);
+    return indices;
+  }
+
+  void shrink_to_fit() {
+    for (auto [name, array]: m_properties)
+      array->shrink_to_fit();
+  }
+
+  /*!
+   * Adds the elements of the other container to this container for each property which is present in this container.
+   *
+   * Gaps in both containers are preserved, and all elements of the other container are guaranteed
+   * to appear after the elements of this container.
+   * Properties in this container which don't appear in the other container are extended with default values.
+   * Properties in the other container which don't appear in this one are not included.
+   * todo: merge() would be useful as well, but could break contiguous regions in the other container
+   *
+   * @param other
+   */
+  void append(const Property_container<Index>& other) {
+
+    m_active_indices.insert(m_active_indices.end(), other.m_active_indices.begin(), other.m_active_indices.end());
+    for (auto [name, array]: m_properties) {
+
+      auto range = other.m_properties.equal_range(name);
+      auto it = range.first;
+      for (; it != range.second; it++) {
+        if (typeid(array.get()) == typeid((it->second.get())))
+          break;
+      }
+
+      if (it != range.second)
+        array->append(*it->second.get());
+      else
+        array->reserve(m_active_indices.size());
+    }
+  }
+
+/*
+  // todo: maybe should be renamed to transfer_from, but I'd rather remove this functionality entirely
+  void transfer(const Property_container<Index>& other, Index other_index, Index this_index) {
+    CGAL_precondition(other.m_property_arrays.size() == m_property_arrays.size());
+    for (auto [name, array]: m_property_arrays) {
+      auto other_array = other.m_property_arrays.at(name);
+      array->transfer_from(*other_array, other_index, this_index);
+    }
+  }*/
+
+  // todo: maybe a compress() method?
+};
+
+}
+
+#endif // DOXYGEN_RUNNING
+
+#endif //CGAL_PROPERTY_CONTAINTER_H

Property_map/test/Property_map/CMakeLists.txt

@@ -8,6 +8,7 @@ create_single_source_cgal_program("test_property_map.cpp")
 create_single_source_cgal_program("dynamic_property_map.cpp")
 create_single_source_cgal_program("dynamic_properties_test.cpp")
 create_single_source_cgal_program("kernel_converter_properties_test.cpp")
+create_single_source_cgal_program("test_Property_container.cpp")
 
 find_package(OpenMesh QUIET)
 if(OpenMesh_FOUND)

Property_map/test/Property_map/test_Property_container.cpp

@@ -0,0 +1,261 @@
+
+#include <CGAL/Property_container.h>
+#include <CGAL/use.h>
+
+using namespace CGAL::Properties::Experimental;
+
+void test_property_creation() {
+
+  Property_container properties;
+
+  // Should return an integer array which didn't previously exist
+  auto [integers, created] = properties.get_or_add_property("integer", 5);
+  static_assert(std::is_same_v<decltype(integers), std::reference_wrapper<Property_array<std::size_t, int>>>);
+  assert(created);
+  assert(properties.num_properties() == 1);
+
+  auto [floats, _] = properties.get_or_add_property<float>("float");
+  static_assert(std::is_same_v<decltype(floats), std::reference_wrapper<Property_array<std::size_t, float>>>);
+  assert(properties.num_properties() == 2);
+
+  // get() should retreive the same arrays
+  assert(integers.get() == properties.get_property<int>("integer"));
+  assert(floats.get() == properties.get_property<float>("float"));
+
+  // remove() should delete a property array & return if it existed
+  assert(!properties.remove_properties("not-a-real-property"));
+  auto removed = properties.remove_property<int>(integers);
+  assert(removed);
+  assert(properties.num_properties() == 1);
+
+  // Add a new property
+  auto [bools, bools_created] = properties.get_or_add_property<bool>("bools", false);
+  static_assert(std::is_same_v<decltype(bools), std::reference_wrapper<Property_array<std::size_t, bool>>>);
+  Property_array<std::size_t, bool>& b = bools.get();
+  CGAL_USE(b);
+}
+
+void test_element_access() {
+
+  Property_container properties;
+
+  auto& integers = properties.add_property("integers", 5);
+
+  // Reserve space for 100 elements
+  properties.reserve(100);
+  assert(properties.capacity() == 100);
+  assert(properties.size() == 0);
+
+  // Newly emplaced elements should go at the front
+  assert(properties.emplace() == 0);
+  assert(properties.emplace() == 1);
+  assert(properties.emplace() == 2);
+  assert(properties.size() == 3);
+
+  // Make sure that the new elements are equal to the default value
+  assert(integers[0] == 5);
+  assert(integers[1] == 5);
+  assert(integers[2] == 5);
+
+  // Add a new property
+  auto& floats = properties.add_property("floats", 6.0f);
+
+  // The new property array should already be of the right size
+  assert(floats.capacity() == 100);
+  assert(properties.size() == 3);
+
+  // Pre-existing elements should contain the default value
+  assert(floats[0] == 6.0f);
+  assert(floats[1] == 6.0f);
+  assert(floats[2] == 6.0f);
+
+  // Update values for a few elements
+  floats[0] = 1.0f;
+  floats[1] = 2.0f;
+  floats[2] = 3.0f;
+  integers[2] = -2;
+  assert(floats[0] == 1.0f);
+  assert(floats[1] == 2.0f);
+  assert(floats[2] == 3.0f);
+  assert(integers[2] == -2);
+
+  // Reset an element, and all of its properties should revert to the defaults
+  properties.reset(2);
+  assert(floats[2] == 6.0f);
+  assert(integers[2] == 5);
+
+  // Erase an element, and the size should be reduced
+  properties.erase(1);
+  assert(properties.size() == 2);
+  assert(properties.capacity() == 100);
+  assert(properties.active_list().size() == 2);
+  assert(properties.inactive_list().size() == 98);
+
+  // A newly emplaced element should take the empty slot
+  assert(properties.emplace() == 1);
+  assert(properties.size() == 3);
+  // todo: should the new element have default properties?
+  assert(properties.emplace() == 3);
+  assert(properties.size() == 4);
+
+  // Swapping a pair of elements swaps all of their properties
+  properties.swap(0, 3);
+  assert(integers[0] == 5);
+  assert(floats[0] == 6.0f);
+  assert(integers[3] == 5);
+  assert(floats[3] == 1.0f);
+
+}
+
+void test_emplace_group() {
+
+  Property_container properties;
+
+  auto& a = properties.add_property("a", 5);
+  CGAL_USE(a);
+  // Insert a group of 100 elements
+  properties.emplace_group(100);
+  assert(properties.size() == 100);
+
+  // Eliminate a few regions
+  properties.erase(3);
+  assert(properties.is_erased(3));
+  assert(properties.size() == 99);
+  for (int i = 20; i < 25; ++i)
+    properties.erase(i);
+  assert(properties.is_erased(23));
+  assert(properties.size() == 94);
+  for (int i = 50; i < 80; ++i)
+    properties.erase(i);
+  assert(properties.is_erased(53));
+  assert(properties.size() == 64);
+
+  // A group of size 4 should only fit in the empty region fo size 5
+  assert(properties.emplace_group(4) == 20);
+  assert(properties.size() == 68);
+  assert(properties.capacity() == 100);
+
+  // A group of size 16 should only fit in the empty region fo size 30
+  assert(properties.emplace_group(16) == 50);
+  assert(properties.size() == 84);
+  assert(properties.capacity() == 100);
+
+  // Another group of size 16 should require the storage to expand, because the largest empty region is mostly full now
+  assert(properties.emplace_group(16) == 100);
+  assert(properties.size() == 100);
+  assert(properties.capacity() == 116);
+
+}
+
+void test_append() {
+
+  // Create a pair of property containers with similar contents
+  Property_container properties_a, properties_b;
+  properties_a.add_property("ints", 1);
+  properties_b.add_property("ints", 2);
+  properties_a.add_property("floats", 3.0f);
+  properties_b.add_property("floats", 4.0f);
+
+  // One container will also contain an extra property
+  properties_a.add_property("bools", true);
+
+  // Add some values to both property sets
+  properties_a.emplace_group(10);
+  properties_b.emplace_group(5);
+  assert(properties_a.size() == 10);
+  assert(properties_b.size() == 5);
+
+  // Add the second group to the end of the first
+  properties_a.append(properties_b);
+  assert(properties_a.size() == 15);
+  assert(properties_b.size() == 5);
+
+  // Initialized values from the second group should appear after those of the first
+  assert(properties_a.get_property<int>("ints")[5] == 1);
+  assert(properties_a.get_property<int>("ints")[12] == 2);
+  assert(properties_a.get_property<float>("floats")[5] == 3.0f);
+  assert(properties_a.get_property<float>("floats")[12] == 4.0f);
+
+  // Additional properties in the first group should have expanded too, and been filled with defaults
+  // note: the property array must be const, because non const operator[] doesn't work for vector<bool>!
+  assert(std::as_const(properties_a).get_property<bool>("bools")[12] == true);
+}
+
+void test_constructors() {
+
+  // Default constructor should have no properties
+  Property_container<std::size_t> a{};
+  assert(a.num_properties() == 0);
+
+  // Copy constructor should duplicate all properties
+  auto& a_ints = a.add_property("ints", 0);
+  auto& a_floats = a.add_property("floats", 0.0f);
+  a.emplace_group(10);
+  a.get_property<int>("ints")[3] = 1;
+  a.get_property<float>("floats")[3] = 1.0f;
+  Property_container<std::size_t> b{a};
+  assert(b.num_properties() == a.num_properties() && b.num_properties() == 2);
+  assert(b.get_property<int>("ints")[3] == a.get_property<int>("ints")[3] && b.get_property<int>("ints")[3] == 1);
+  assert(b.get_property<float>("floats")[3] == a.get_property<float>("floats")[3] && b.get_property<float>("floats")[3] == 1.0f);
+
+  // Copy-assignment operator should do effectively the same thing as the copy constructor
+  Property_container<std::size_t> c;
+  c = a;
+  assert(c.num_properties() == a.num_properties() && c.num_properties() == 2);
+  assert(c.get_property<int>("ints")[3] == a.get_property<int>("ints")[3] && c.get_property<int>("ints")[3] == 1);
+  assert(c.get_property<float>("floats")[3] == a.get_property<float>("floats")[3] && c.get_property<float>("floats")[3] == 1.0f);
+
+  // Copied property containers should not be synced with the original
+  a.add_property("more_ints", 2);
+  assert(a.num_properties() == 3);
+  assert(b.num_properties() == 2);
+  assert(c.num_properties() == 2);
+  a.get_property<int>("ints")[4] = 2;
+  assert(a.get_property<int>("ints")[4] == 2);
+  assert(b.get_property<int>("ints")[4] == 0);
+  assert(c.get_property<int>("ints")[4] == 0);
+
+  // Copy assignment should not invalidate previously obtained array references,
+  // but it should update their values
+  auto &b_ints = b.get_property<int>("ints");
+  auto &b_floats = b.get_property<float>("floats");
+  assert(b_ints[4] == 0);
+  b = a;
+  assert(b.num_properties() == 3);
+  assert(b_ints[4] == 2);
+
+  // Move assignment shouldn't invalidate references either
+  Property_container<std::size_t> d{c};
+  auto &d_ints = d.get_property<int>("ints");
+  assert(d_ints[4] == 0);
+  d = std::move(a);
+  assert(d.num_properties() == 3);
+  assert(d_ints[4] == 2);
+
+  // Moved-from should be empty
+  // All properties are preserved, though
+  assert(a.num_properties() == 3);
+  assert(a.size() == 0);
+  assert(a_ints.capacity() == 0);
+  assert(a_floats.capacity() == 0);
+
+  // Move constructor should behave like move assignment
+  Property_container<std::size_t> e{std::move(b)};
+  assert(e.num_properties() == 3);
+  assert(b.num_properties() == 3);
+  assert(b.size() == 0);
+  assert(b_ints.capacity() == 0);
+  assert(b_floats.capacity() == 0);
+}
+
+
+int main() {
+
+  test_property_creation();
+  test_element_access();
+  test_emplace_group();
+  test_append();
+  test_constructors();
+
+  return 0;
+}

Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h

@@ -253,12 +253,12 @@ public:
   /*!
     Retrieves the point property map.
   */
-  const Point_map &point_map() const { return m_point_pmap; }
+  const Point_map &point_map() const { return *m_point_pmap; }
 
   /*!
     Retrieves the normal property map.
   */
-  const Normal_map &normal() const { return m_normal_pmap; }
+  const Normal_map &normal() const { return *m_normal_pmap; }
 
   Input_iterator input_iterator_first() const {
     return m_input_iterator_first;
@@ -361,13 +361,13 @@ public:
         m_direct_octrees[s] = new Direct_octree(
                 m_traits, last + 1,
                 last + subsetSize + 1,
-                m_point_pmap,
+                *m_point_pmap,
                 remainingPoints - subsetSize);
       } else
         m_direct_octrees[0] = new Direct_octree(
                 m_traits, m_input_iterator_first,
                 m_input_iterator_first + (subsetSize),
-                m_point_pmap,
+                *m_point_pmap,
                 0);
 
       m_available_octree_sizes[s] = subsetSize;
@@ -378,7 +378,7 @@ public:
 
     m_global_octree = new Indexed_octree(
             m_traits, m_input_iterator_first, m_input_iterator_beyond,
-            m_point_pmap
+            *m_point_pmap
     );
     m_global_octree->refine(m_options.cluster_epsilon);
 
@@ -496,7 +496,7 @@ public:
     }
 
     // Use bounding box diagonal as reference for default values
-    Bbox_3 bbox = m_global_octree->boundingBox();
+    auto bbox = m_global_octree->boundingBox();
     FT bbox_diagonal = (FT) CGAL::sqrt(
             (bbox.xmax() - bbox.xmin()) * (bbox.xmax() - bbox.xmin())
             + (bbox.ymax() - bbox.ymin()) * (bbox.ymax() - bbox.ymin())
@@ -574,14 +574,14 @@ public:
                     static_cast<unsigned int>(m_num_available_points));
                 while (m_shape_index[first_sample] != -1);
 
-                done = drawSamplesFromCellContainingPoint
-                  (m_global_octree,
-                   get(m_point_pmap,
-                       *(m_input_iterator_first + first_sample)),
-                   select_random_octree_level(),
-                   indices,
-                   m_shape_index,
-                   m_required_samples);
+                done = drawSamplesFromCellContainingPoint(
+                  m_global_octree,
+                  get(*m_point_pmap, *(m_input_iterator_first + first_sample)),
+                  select_random_octree_level(),
+                  indices,
+                  m_shape_index,
+                  m_required_samples
+                );
 
                 if (callback && !callback(num_invalid / double(m_num_total_points))) {
                   clear(num_invalid, candidates);
@@ -605,8 +605,8 @@ public:
                 p->compute(indices,
                            m_input_iterator_first,
                            m_traits,
-                           m_point_pmap,
-                           m_normal_pmap,
+                           *m_point_pmap,
+                           *m_normal_pmap,
                            m_options.epsilon,
                            m_options.normal_threshold);
 
@@ -1086,7 +1086,7 @@ private:
       Cell cell = stack.top();
       stack.pop();
 
-      FT width = octree->width() / (1 << (cell.depth()));
+      FT width = octree->width() / (1 << (octree->depth(cell)));
 
       FT diag = CGAL::sqrt(FT(3) * width * width) + epsilon;
 
@@ -1097,10 +1097,10 @@ private:
 
       // differ between full or partial overlap?
       // if full overlap further traversal of this branch is not necessary
-      if (cell.is_leaf()) {
+      if (octree->is_leaf(cell)) {
         std::vector<std::size_t> indices;
-        indices.reserve(cell.size());
-        for (std::size_t i = 0; i < cell.size(); i++) {
+        indices.reserve(octree->points(cell).size());
+        for (std::size_t i = 0; i < octree->points(cell).size(); i++) {
           if (shapeIndex[octree->index(cell, i)] == -1) {
             indices.push_back(octree->index(cell, i));
           }
@@ -1111,10 +1111,10 @@ private:
                                  indices);
       } else {
 
-        if (!cell.is_leaf()) {
+        if (!octree->is_leaf(cell)) {
           for (std::size_t i = 0; i < 8; i++) {
-            if (!cell[i].empty())
-              stack.push(cell[i]);
+            if (octree->points(octree->child(cell, i)).size() != 0)
+              stack.push(octree->child(cell, i));
           }
         }
       }
@@ -1129,30 +1129,11 @@ private:
   const typename Octree::Node node_containing_point(const Octree *octree, const Point &p, std::size_t level) {
 
     // Find the node containing the point
-    typename Octree::Node cur = octree->root();
-    while (!cur.is_null() && cur.depth() < level) {
+    typename Octree::Node n = octree->locate(p);
+    while (octree->depth(n) > level)
+      n = octree->parent(n);
 
-      // If cur is a leaf node, its child is null
-      if (cur.is_leaf())
-        return typename Octree::Node();
-
-      // If that child is empty, return null
-      if (cur.empty())
-        return typename Octree::Node();
-
-      // Determine the coordinate of the child
-      Point center = octree->barycenter(cur);
-      std::bitset<3> coordinate;
-      coordinate[0] = center.x() <= p.x();
-      coordinate[1] = center.y() <= p.y();
-      coordinate[2] = center.z() <= p.z();
-
-      // Otherwise, return the correct child of cur
-      cur = cur[coordinate.to_ulong()];
-
-    }
-
-    return cur;
+    return n;
   }
 
   template<class Octree>
@@ -1167,13 +1148,9 @@ private:
 
     const Cell cur = node_containing_point(octree, p, level);
 
-    // Stop if the node we need doesn't exist
-    if (cur.is_null())
-      return false;
-
     // Count point indices that map to -1 in the shape index
     std::size_t enough = 0;
-    for (auto j : cur) {
+    for (const auto j : octree->points(cur)) {
       if (shapeIndex[j] == -1)
         enough++;
       if (enough >= requiredSamples)
@@ -1186,7 +1163,7 @@ private:
 
     do {
       std::size_t p = CGAL::get_default_random().
-              uniform_int<std::size_t>(0, cur.size() - 1);
+              uniform_int<std::size_t>(0, octree->points(cur).size() - 1);
       std::size_t j = octree->index(cur, p);
 
       if (shapeIndex[j] == -1)
@@ -1225,8 +1202,8 @@ private:
   // iterators of input data
   bool m_valid_iterators;
   Input_iterator m_input_iterator_first, m_input_iterator_beyond;
-  Point_map m_point_pmap;
-  Normal_map m_normal_pmap;
+  std::optional<Point_map> m_point_pmap;
+  std::optional<Normal_map> m_normal_pmap;
 };
 
 

Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC_traits.h

@@ -39,6 +39,8 @@ namespace CGAL {
             class InputPointMap,
             class InputNormalMap>
   struct Efficient_RANSAC_traits {
+    ///
+    typedef Gt GeomTraits;
     ///
     typedef typename Gt::FT FT;
     ///

Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Octree.h

@@ -32,7 +32,7 @@ namespace CGAL {
 namespace Shape_detection {
 
 // Forward declaration needed for automatic traits detection without
-// including the deprecated header itself…
+// including the deprecated header itself
 template <typename Gt, typename IR, typename IPM, typename INM>
 struct Shape_detection_traits;
 
@@ -56,8 +56,9 @@ class RANSAC_octree {
   typedef std::vector<std::size_t> Input_range;
   typedef Random_index_access_property_map<Input_iterator, Point_map> Indexed_point_map;
 
-  typedef CGAL::Octree<typename Traits_base<Traits>::type,
-                       Input_range, Indexed_point_map> Octree;
+  typedef Orthtree_traits_point<typename Traits_base<Traits>::type, Input_range, Indexed_point_map, false, 3> OTraits;
+
+  typedef CGAL::Orthtree<OTraits> Octree;
 
   Traits m_traits;
   Input_range m_input_range;
@@ -70,7 +71,8 @@ class RANSAC_octree {
 
 public:
 
-  typedef typename Octree::Node Node;
+  typedef typename Octree::Node_index Node;
+  typedef typename OTraits::Node_data Node_data;
 
   RANSAC_octree(const Traits &traits,
                 Input_iterator begin,
@@ -81,18 +83,26 @@ public:
           m_input_range(boost::counting_iterator<std::size_t>(0),
                         boost::counting_iterator<std::size_t>(end - begin)),
           m_index_map(begin, point_map),
-          m_octree(m_input_range, m_index_map, 1.0),
+          m_octree(OTraits(m_input_range, m_index_map)),
           m_bBox (bbox_3(make_transform_iterator_from_property_map(begin, point_map),
                          make_transform_iterator_from_property_map(end, point_map))),
           m_offset(offset) {}
 
   std::size_t index (Node node, std::size_t i) const
   {
-    return m_offset + *(node.begin() + i);
+    return m_offset + *(m_octree.data(node).begin() + i);
+  }
+
+  std::size_t depth(const Node& node) const {
+    return m_octree.depth(node);
+  }
+
+  bool is_leaf(const Node& node) const {
+    return m_octree.is_leaf(node);
   }
 
   std::size_t size() const {
-    return m_octree.root().size();
+    return m_input_range.size();
   }
 
   std::size_t maxLevel() const {
@@ -127,17 +137,27 @@ public:
     return m_width;
   }
 
+  Node child(const Node& node, std::size_t i) const {
+    return m_octree.child(node, i);
+  }
+
+  Node parent(const Node& node) const {
+    return m_octree.parent(node);
+  }
+
   Node locate(const typename Traits::Point_3 &p) const {
     return m_octree.locate(p);
   }
 
   Node root() const { return m_octree.root(); }
 
+  Node_data points(const Node& n) const { return m_octree.data(n); }
+
   typename Traits::Point_3 barycenter(const Node &node) const {
     return m_octree.barycenter(node);
   }
 
-  Bbox_3 boundingBox() const {
+  typename Traits::GeomTraits::Iso_cuboid_3 boundingBox() const {
       return m_octree.bbox(m_octree.root());
   }
 };

Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Shape_base.h

@@ -374,7 +374,7 @@ namespace CGAL {
      */
     typename boost::property_traits< typename Traits::Point_map >::reference
     point(std::size_t i) const {
-      return get(this->m_point_pmap, *(this->m_first + i));
+      return get(this->m_point_pmap.value(), *(this->m_first + i));
     }
 
     /*!
@@ -382,7 +382,7 @@ namespace CGAL {
      */
     typename boost::property_traits< typename Traits::Normal_map >::reference
     normal(std::size_t i) const {
-      return get(this->m_normal_pmap, *(this->m_first + i));
+      return get(this->m_normal_pmap.value(), *(this->m_first + i));
     }
 
     /*!
@@ -694,8 +694,8 @@ namespace CGAL {
     Input_iterator m_first;
 
     Traits m_traits;
-    Point_map m_point_pmap;
-    Normal_map m_normal_pmap;
+    std::optional<Point_map> m_point_pmap;
+    std::optional<Normal_map> m_normal_pmap;
     /// \endcond
   };
 }

Shape_detection/package_info/Shape_detection/dependencies

@@ -6,6 +6,7 @@ Circulator
 Distance_2
 Distance_3
 Filtered_kernel
+Hash_map
 Homogeneous_kernel
 Installation
 Intersections_2
@@ -16,6 +17,7 @@ Kernel_d
 Modular_arithmetic
 Number_types
 Orthtree
+Point_set_2
 Principal_component_analysis
 Principal_component_analysis_LGPL
 Profiling_tools
@@ -25,4 +27,7 @@ STL_Extension
 Shape_detection
 Solver_interface
 Spatial_searching
+Spatial_sorting
 Stream_support
+TDS_2
+Triangulation_2

Shape_regularization/package_info/Shape_regularization/dependencies

@@ -1,5 +1,6 @@
 Algebraic_foundations
 BGL
+Cartesian_kernel
 Circulator
 Distance_2
 Distance_3
@@ -14,6 +15,7 @@ Kernel_d
 Modular_arithmetic
 Number_types
 Orthtree
+Point_set_2
 Principal_component_analysis_LGPL
 Profiling_tools
 Property_map