Splitting is done by index (internally)

Pre-allocating nodes is no longer necessary, since refine() isn't broken by pointer invalidation

Orthtree/include/CGAL/Orthtree.h

@@ -198,7 +198,7 @@ public:
       , m_point_map(point_map) {
 
     // fixme: this can be removed once traversal doesn't require pointer stability
-    m_nodes.reserve(1'000'000);
+    //m_nodes.reserve(1'000'000);
     m_nodes.emplace_back();
 
     Array bbox_min;
@@ -311,8 +311,8 @@ public:
     m_side_per_depth.resize(1);
 
     // Initialize a queue of nodes that need to be refined
-    std::queue<Node*> todo;
-    todo.push(&root());
+    std::queue<std::size_t> todo;
+    todo.push(0);
 
     // Process items in the queue until it's consumed fully
     while (!todo.empty()) {
@@ -322,21 +322,21 @@ public:
       todo.pop();
 
       // Check if this node needs to be processed
-      if (split_predicate(*current)) {
+      if (split_predicate(m_nodes[current])) {
 
         // Check if we've reached a new max depth
-        if (current->depth() == depth()) {
+        if (m_nodes[current].depth() == depth()) {
 
           // Update the side length map
           m_side_per_depth.push_back(*(m_side_per_depth.end() - 1) / 2);
         }
 
         // Split the node, redistributing its points to its children
-        split((*current));
+        split(current);
 
         // Process each of its children
         for (int i = 0; i < Degree::value; ++i)
-          todo.push(&children(*current)[i]);
+          todo.push(m_nodes[current].m_children_index.get() + i);
 
       }
     }
@@ -780,23 +780,26 @@ public:
   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(Node& node) {
+  void split(Node& node) { split(index(node)); }
+
+  void split(std::size_t n) {
 
     // Make sure the node hasn't already been split
-    CGAL_precondition (node.is_leaf());
+    CGAL_precondition (m_nodes[n].is_leaf());
 
     // Split the node to create children
     using Local_coordinates = typename Node::Local_coordinates;
     for (int i = 0; i < Degree::value; i++) {
-      m_nodes.emplace_back(&node, index(node), Local_coordinates{i});
+      m_nodes.emplace_back(n, m_nodes[n].global_coordinates(), m_nodes[n].depth() + 1, Local_coordinates{i});
     }
-    node.m_children_index = m_nodes.size() - Degree::value;
+    // todo: this assumes that the new nodes always are allocated at the end
+    m_nodes[n].m_children_index = m_nodes.size() - Degree::value;
 
     // Find the point to around which the node is split
-    Point center = barycenter(node);
+    Point center = barycenter(m_nodes[n]);
 
     // Add the node's points to its children
-    reassign_points(node, node.points().begin(), node.points().end(), center);
+    reassign_points(m_nodes[n], m_nodes[n].points().begin(), m_nodes[n].points().end(), center);
   }
 
   /*!
@@ -810,6 +813,10 @@ public:
     node.m_children_index.reset();
   }
 
+  void unsplit(std::size_t n) {
+    unsplit(m_nodes[n]);
+  }
+
   // todo: documentation
   Point barycenter(const Node& node) const {
 

Orthtree/include/CGAL/Orthtree/Node.h

@@ -120,7 +120,7 @@ public:
   /// @{
 
   /// \cond SKIP_IN_MANUAL
-  Node() = default;
+  Node() = default; // constructs a root node
   /// \endcond
 
   /*!
@@ -138,23 +138,13 @@ public:
     \param parent the node containing this one
     \param index this node's relationship to its parent
   */
-  explicit Node(Self* parent, boost::optional<std::size_t> parent_index, Local_coordinates local_coordinates) :
-    m_parent_index(parent_index) {
+  explicit Node(std::size_t parent_index, Global_coordinates parent_coordinates,
+                std::size_t depth, Local_coordinates local_coordinates) :
+    m_parent_index(parent_index), m_depth(depth) {
 
-    // todo: this can be cleaned up; it probably doesn't need to take a reference to the parent
+    for (int i = 0; i < Dimension::value; i++)
+      m_global_coordinates[i] = (2 * parent_coordinates[i]) + local_coordinates[i];
 
-    if (parent != nullptr) {
-      m_depth = parent->m_depth + 1;
-
-      for (int i = 0; i < Dimension::value; i++)
-        m_global_coordinates[i] = (2 * parent->m_global_coordinates[i]) + local_coordinates[i];
-
-    } else {
-      m_depth = 0;
-
-      for (int i = 0; i < Dimension::value; i++)
-        m_global_coordinates[i] = 0;
-    }
   }
 
   /// @}
@@ -177,11 +167,6 @@ public:
   /// \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()`