diff --git a/Orthtree/include/CGAL/Octree.h b/Orthtree/include/CGAL/Octree.h
index bd7e50b3b88..848b2731c73 100644
--- a/Orthtree/include/CGAL/Octree.h
+++ b/Orthtree/include/CGAL/Octree.h
@@ -35,11 +35,13 @@ namespace CGAL {
   \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`
  */
-template <typename GeomTraits, typename PointRange,
-          typename PointMap = Identity_property_map
-         <typename std::iterator_traits<typename PointRange::iterator>::value_type> >
+template <
+  typename GeomTraits,
+  typename PointRange,
+  typename PointMap = Identity_property_map<typename std::iterator_traits<typename PointRange::iterator>::value_type>
+>
 #ifdef DOXYGEN_RUNNING
-class Octree;
+  class Octree;
 #else
 using Octree = Orthtree<Orthtree_traits_3<GeomTraits>, PointRange, PointMap>;
 #endif
diff --git a/Orthtree/include/CGAL/Orthtree.h b/Orthtree/include/CGAL/Orthtree.h
index afc958b9754..43b42a53782 100644
--- a/Orthtree/include/CGAL/Orthtree.h
+++ b/Orthtree/include/CGAL/Orthtree.h
@@ -40,7 +40,8 @@
 #include <vector>
 #include <math.h>
 
-namespace CGAL {
+namespace CGAL
+{
 
 /*!
   \ingroup PkgOrthtreeClasses
@@ -63,8 +64,8 @@ namespace CGAL {
   \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 `Traits_::Point_d`
  */
-template<typename Traits_, typename PointRange_,
-         typename PointMap_ = Identity_property_map<typename Traits_::Point_d> >
+template <typename Traits_, typename PointRange_,
+  typename PointMap_ = Identity_property_map<typename Traits_::Point_d> >
 class Orthtree
 {
 
@@ -88,9 +89,9 @@ public:
   /// \cond SKIP_IN_MANUAL
   typedef typename Traits::Array Array; ///< Array type.
   typedef typename Traits::Construct_point_d_from_array
-  Construct_point_d_from_array;
+    Construct_point_d_from_array;
   typedef typename Traits::Construct_bbox_d
-  Construct_bbox_d;
+    Construct_bbox_d;
   /// \endcond
   /// @}
 
@@ -105,7 +106,7 @@ public:
   /*!
    * \brief Degree of the tree (number of children of non-leaf nodes).
    */
-  typedef Dimension_tag<(2 << (Dimension::value-1))> Degree;
+  typedef Dimension_tag<(2 << (Dimension::value - 1))> Degree;
 
   /*!
    * \brief The Sub-tree / Orthant type.
@@ -123,7 +124,7 @@ public:
 #ifdef DOXYGEN_RUNNING
   typedef unspecified_type Node_range;
 #else
-  typedef boost::iterator_range<Traversal_iterator<Node> > Node_range;
+  typedef boost::iterator_range<Traversal_iterator<const Node>> Node_range;
 #endif
 
   /// \cond SKIP_IN_MANUAL
@@ -131,7 +132,7 @@ public:
   /*!
    * \brief A function that determines the next node in a traversal given the current one.
    */
-  typedef std::function<Node(Node)> Node_traversal_method_const;
+  typedef std::function<const Node *(const Node *)> Node_traversal_method_const;
 
   /// \endcond
 
@@ -193,48 +194,48 @@ public:
            PointMap point_map = PointMap(),
            const FT enlarge_ratio = 1.2,
            Traits traits = Traits())
-    : m_traits (traits)
-    , m_range (point_range)
-    , m_point_map (point_map)
-    , m_root(Node(), 0)
+    : m_traits(traits)
+      , m_range(point_range)
+      , m_point_map(point_map)
+      , m_root() // todo: can this be default-constructed?
   {
     Array bbox_min;
     Array bbox_max;
 
     // init bbox with first values found
     {
-      const Point& p = get (m_point_map, *(point_range.begin()));
+      const Point& p = get(m_point_map, *(point_range.begin()));
       std::size_t i = 0;
-      for (const FT& x : cartesian_range(p))
+      for (const FT& x: cartesian_range(p))
       {
         bbox_min[i] = x;
         bbox_max[i] = x;
-        ++ i;
+        ++i;
       }
     }
 
-    for (const Range_type& r : point_range)
+    for (const Range_type& r: point_range)
     {
-      const Point& p = get (m_point_map, r);
+      const Point& p = get(m_point_map, r);
       std::size_t i = 0;
-      for (const FT& x : cartesian_range(p))
+      for (const FT& x: cartesian_range(p))
       {
         bbox_min[i] = (std::min)(x, bbox_min[i]);
         bbox_max[i] = (std::max)(x, bbox_max[i]);
-        ++ i;
+        ++i;
       }
     }
 
     Array bbox_centroid;
     FT max_length = FT(0);
-    for (std::size_t i = 0 ; i < Dimension::value; ++ i)
+    for (std::size_t i = 0; i < Dimension::value; ++i)
     {
       bbox_centroid[i] = (bbox_min[i] + bbox_max[i]) / FT(2);
       max_length = (std::max)(max_length, bbox_max[i] - bbox_min[i]);
     }
     max_length *= enlarge_ratio / FT(2);
 
-    for (std::size_t i = 0 ; i < Dimension::value; ++ i)
+    for (std::size_t i = 0; i < Dimension::value; ++i)
     {
       bbox_min[i] = bbox_centroid[i] - max_length;
       bbox_max[i] = bbox_centroid[i] + max_length;
@@ -254,47 +255,31 @@ public:
   /// \cond SKIP_IN_MANUAL
 
   // copy constructor
-  Orthtree (const Orthtree& other)
-    : m_traits (other.m_traits)
-    , m_range (other.m_range)
-    , m_point_map (other.m_point_map)
-    , m_root (other.m_root.deep_copy())
-    , m_bbox_min (other.m_bbox_min)
-    , m_side_per_depth(other.m_side_per_depth)
-  { }
+  Orthtree(const Orthtree& other)
+    : m_traits(other.m_traits)
+      , m_range(other.m_range)
+      , m_point_map(other.m_point_map)
+      , m_root(other.m_root.deep_copy())
+      , m_bbox_min(other.m_bbox_min)
+      , m_side_per_depth(other.m_side_per_depth) {}
 
   // move constructor
-  Orthtree (Orthtree&& other)
-    : m_traits (other.m_traits)
-    , m_range (other.m_range)
-    , m_point_map (other.m_point_map)
-    , m_root (other.m_root)
-    , m_bbox_min (other.m_bbox_min)
-    , m_side_per_depth(other.m_side_per_depth)
-  {
-    other.m_root = Node(Node(), 0);
+  Orthtree(Orthtree&& other)
+    : m_traits(other.m_traits)
+      , m_range(other.m_range)
+      , m_point_map(other.m_point_map)
+      , m_root(other.m_root)
+      , m_bbox_min(other.m_bbox_min)
+      , m_side_per_depth(other.m_side_per_depth) {
+    other.m_root = Node{};
   }
 
   // Non-necessary but just to be clear on the rule of 5:
 
   // assignment operators deleted (PointRange is a ref)
-  Orthtree& operator= (const Orthtree& other) = delete;
-  Orthtree& operator= (Orthtree&& other) = delete;
-  // Destructor
-  ~Orthtree()
-  {
-    std::queue<Node> nodes;
-    nodes.push(m_root);
-    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();
-    }
-  }
+  Orthtree& operator=(const Orthtree& other) = delete;
+
+  Orthtree& operator=(Orthtree&& other) = delete;
 
   // move constructor
   /// \endcond
@@ -320,52 +305,39 @@ public:
   void refine(const Split_predicate& split_predicate) {
 
     // If the tree has already been refined, reset it
-    if (!m_root.is_leaf()){
-      std::queue<Node> nodes;
-      for (std::size_t i = 0; i < Degree::value; ++ i)
-            nodes.push (m_root[i]);
-      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();
-      }
+    if (!m_root.is_leaf())
       m_root.unsplit();
-    }
 
     // Reset the side length map, too
     m_side_per_depth.resize(1);
 
     // Initialize a queue of nodes that need to be refined
-    std::queue<Node> todo;
-    todo.push(m_root);
+    std::queue<Node *> todo;
+    todo.push(&m_root);
 
     // Process items in the queue until it's consumed fully
     while (!todo.empty()) {
 
       // Get the next element
-      Node current = todo.front();
+      auto current = todo.front();
       todo.pop();
 
       // Check if this node needs to be processed
-      if (split_predicate(current)) {
+      if (split_predicate(*current)) {
 
         // Check if we've reached a new max depth
-        if (current.depth() == depth()) {
+        if (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(current[i]);
+          todo.push(&(*current)[i]);
 
       }
     }
@@ -401,52 +373,52 @@ public:
   void grade() {
 
     // Collect all the leaf nodes
-    std::queue<Node> leaf_nodes;
-    for (Node leaf : traverse(Orthtrees::Leaves_traversal())) {
+    std::queue<Node *> leaf_nodes;
+    for (const Node &leaf : traverse(Orthtrees::Leaves_traversal())) {
       // TODO: I'd like to find a better (safer) way of doing this
-      leaf_nodes.push(leaf);
+      leaf_nodes.push(const_cast<Node *>(&leaf));
     }
 
     // Iterate over the nodes
     while (!leaf_nodes.empty()) {
 
       // Get the next node
-      Node node = leaf_nodes.front();
+      Node *node = leaf_nodes.front();
       leaf_nodes.pop();
 
       // Skip this node if it isn't a leaf anymore
-      if (!node.is_leaf())
+      if (!node->is_leaf())
         continue;
 
       // Iterate over each of the neighbors
       for (int direction = 0; direction < 6; ++direction) {
 
         // Get the neighbor
-        Node neighbor = node.adjacent_node(direction);
+        auto *neighbor = node->adjacent_node(direction);
 
         // If it doesn't exist, skip it
-        if (neighbor.is_null())
+        if (!neighbor)
           continue;
 
         // Skip if this neighbor is a direct sibling (it's guaranteed to be the same depth)
         // TODO: This check might be redundant, if it doesn't affect performance maybe I could remove it
-        if (neighbor.parent() == node.parent())
+        if (neighbor->parent() == node->parent())
           continue;
 
         // If it's already been split, skip it
-        if (!neighbor.is_leaf())
+        if (!neighbor->is_leaf())
           continue;
 
         // Check if the neighbor breaks our grading rule
         // TODO: could the rule be parametrized?
-        if ((node.depth() - neighbor.depth()) > 1) {
+        if ((node->depth() - neighbor->depth()) > 1) {
 
           // Split the neighbor
-          split(neighbor);
+          split(*neighbor);
 
           // Add newly created children to the queue
           for (int i = 0; i < Degree::value; ++i) {
-            leaf_nodes.push(neighbor[i]);
+            leaf_nodes.push(&(*neighbor)[i]);
           }
         }
       }
@@ -459,9 +431,22 @@ public:
   /// @{
 
   /*!
-    \brief returns the root node.
+    \brief provides read-only access to the root node, and by
+    extension the rest of the tree.
+
+    \return a const reference to the root node of the tree.
    */
-  Node root() const { return m_root; }
+  const Node &root() const { return m_root; }
+
+  /*!
+    \brief provides read-write access to the root node, and by
+    extension the rest of the tree.
+
+    todo: why wasn't this provided previously?
+
+    \return a reference to the root node of the tree.
+   */
+  Node &root() { return m_root; }
 
   /*!
     \brief Convenience function to access the child nodes of the root
@@ -472,9 +457,9 @@ public:
     \sa `Node::operator[]()`
 
     \param index the index of the child node.
-    \return the accessed node.
+    \return a reference to the node.
    */
-  Node operator[](std::size_t index) const { return m_root[index]; }
+  const Node &operator[](std::size_t index) const { return m_root[index]; }
 
   /*!
     \brief returns the deepest level reached by a leaf node in this tree (root being level 0).
@@ -497,13 +482,13 @@ public:
   template<typename Traversal>
   Node_range traverse(const Traversal &traversal = Traversal()) const {
 
-    Node first = traversal.first(m_root);
+    const Node *first = traversal.first(&m_root);
 
     Node_traversal_method_const next
-      = [&](const Node& n) -> Node { return traversal.next(n); };
+      = [&](const Node* n) -> const Node * { return traversal.next(n); };
 
-    return boost::make_iterator_range(Traversal_iterator<Node>(first, next),
-                                      Traversal_iterator<Node>());
+    return boost::make_iterator_range(Traversal_iterator<const Node>(first, next),
+                                      Traversal_iterator<const Node>());
   }
 
   /*!
@@ -547,19 +532,19 @@ public:
     \param point query point.
     \return the node which contains the point.
    */
-  Node locate(const Point &point) const {
+  const Node& locate(const Point &point) const {
 
     // Make sure the point is enclosed by the orthtree
     CGAL_precondition (CGAL::do_intersect(point, bbox(m_root)));
 
     // Start at the root node
-    auto node_for_point = m_root;
+    auto *node_for_point = &m_root;
 
     // Descend the tree until reaching a leaf node
-    while (!node_for_point.is_leaf()) {
+    while (!node_for_point->is_leaf()) {
 
       // Find the point to split around
-      Point center = barycenter(node_for_point);
+      Point center = barycenter(*node_for_point);
 
       // Find the index of the correct sub-node
       typename Node::Local_coordinates index;
@@ -568,11 +553,11 @@ public:
         index[dimension ++] = (get<0>(r) < get<1>(r));
 
       // Find the correct sub-node of the current node
-      node_for_point = node_for_point[index.to_ulong()];
+      node_for_point = &(*node_for_point)[index.to_ulong()];
     }
 
     // Return the result
-    return node_for_point;
+    return *node_for_point;
   }
 
   /*!
@@ -582,15 +567,15 @@ public:
     `query`.
 
     \tparam OutputIterator a model of `OutputIterator` that accept `Point_d` objects.
-    \param query a query point.
-    \param k the number of neighbors.
-    \param output the output iterator.
+    \param query query point.
+    \param k number of neighbors.
+    \param output output iterator.
    */
-  template<typename OutputIterator>
-  OutputIterator nearest_neighbors (const Point& query,
-                                    std::size_t k,
-                                    OutputIterator output) const {
-    Sphere query_sphere (query, (std::numeric_limits<FT>::max)());
+  template <typename OutputIterator>
+  OutputIterator nearest_neighbors(const Point& query,
+                                   std::size_t k,
+                                   OutputIterator output) const {
+    Sphere query_sphere(query, (std::numeric_limits<FT>::max)());
     return nearest_k_neighbors_in_radius(query_sphere, k, output);
   }
 
@@ -604,8 +589,8 @@ public:
     \param query query sphere.
     \param output output iterator.
    */
-  template<typename OutputIterator>
-  OutputIterator nearest_neighbors (const Sphere& query, OutputIterator output) const {
+  template <typename OutputIterator>
+  OutputIterator nearest_neighbors(const Sphere& query, OutputIterator output) const {
     Sphere query_sphere = query;
     return nearest_k_neighbors_in_radius(query_sphere,
                                          (std::numeric_limits<std::size_t>::max)(), output);
@@ -625,7 +610,7 @@ public:
     \param output output iterator.
    */
   template<typename Query, typename OutputIterator>
-  OutputIterator intersected_nodes (const Query& query, OutputIterator output) const {
+  OutputIterator intersected_nodes(const Query& query, OutputIterator output) const {
     return intersected_nodes_recursive(query, root(), output);
   }
 
@@ -666,7 +651,7 @@ public:
 
   // TODO: Document this
   // TODO: Could this method name be reduced to just "center" ?
-  Point barycenter(const Node& node) const {
+  Point barycenter(const Node &node) const {
 
     // Determine the side length of this node
     FT size = m_side_per_depth[node.depth()];
@@ -677,7 +662,7 @@ public:
     for (const FT& f : cartesian_range(m_bbox_min))
     {
       bary[i] = FT(node.global_coordinates()[i]) * size + size / FT(2) + f;
-      ++ i;
+      ++i;
     }
 
     // Convert that location into a point
@@ -703,11 +688,12 @@ private: // functions :
     // Split the point collection around the center point on this dimension
     Range_iterator split_point = std::partition
       (begin, end,
-       [&](const Range_type &a) -> bool {
-        // This should be done with cartesian iterator but it seems
-        // complicated to do efficiently
+       [&](const Range_type& a) -> bool
+       {
+         // This should be done with cartesian iterator but it seems
+         // complicated to do efficiently
          return (get(m_point_map, a)[int(dimension)] < center[int(dimension)]);
-      });
+       });
 
     // Further subdivide the first side of the split
     std::bitset<Dimension::value> coord_left = coord;
@@ -721,7 +707,7 @@ private: // functions :
 
   }
 
-  void split(Node& node) {
+  void split(Node &node) {
 
     // Make sure the node hasn't already been split
     CGAL_precondition (node.is_leaf());
@@ -755,18 +741,19 @@ private: // functions :
     Point point;
     FT distance;
   };
-  struct Node_index_with_distance {
+
+  struct Node_index_with_distance
+  {
     typename Node::Local_coordinates index;
     FT distance;
 
-    Node_index_with_distance (const typename Node::Local_coordinates& index,
-                              const FT& distance)
-      : index(index), distance(distance)
-    { }
+    Node_index_with_distance(const typename Node::Local_coordinates& index,
+                             const FT& distance)
+      : index(index), distance(distance) {}
   };
 
-  void nearest_k_neighbors_recursive(Sphere& search_bounds, const Node &node,
-                                     std::vector<Point_with_distance> &results, FT epsilon = 0) const {
+  void nearest_k_neighbors_recursive(Sphere &search_bounds, const Node& node,
+                                     std::vector<Point_with_distance>& results, FT epsilon = 0) const {
 
     // Check whether the node has children
     if (node.is_leaf()) {
@@ -820,7 +807,7 @@ private: // functions :
 
       // Fill the list with child nodes
       for (int index = 0; index < Degree::value; ++index) {
-        Node child_node = node[index];
+        auto &child_node = node[index];
 
         // Add a child to the list, with its distance
         children_with_distances.emplace_back(typename Node::Local_coordinates(index),
@@ -834,7 +821,7 @@ private: // functions :
 
       // Loop over the children
       for (auto child_with_distance : children_with_distances) {
-        Node child_node = node[child_with_distance.index.to_ulong()];
+        auto &child_node = node[child_with_distance.index.to_ulong()];
 
         // Check whether the bounding box of the child intersects with the search bounds
         if (do_intersect(child_node, search_bounds)) {
@@ -853,9 +840,9 @@ private: // functions :
     // Check if the current node intersects with the query
     if (CGAL::do_intersect(query, bbox(node))) {
 
-      // if this node is a leaf, than it's considered an intersecting node
+      // if this node is a leaf, then it's considered an intersecting node
       if (node.is_leaf()) {
-        *output++ = node;
+        *output++ = &node;
         return output;
       }
 
@@ -884,9 +871,9 @@ private: // functions :
     \param k the number of points to find
     \param output the output iterator to add the found points to (in order of increasing distance)
    */
-  template<typename OutputIterator>
+  template <typename OutputIterator>
   OutputIterator nearest_k_neighbors_in_radius
-  (Sphere& query_sphere,
+  (Sphere &query_sphere,
    std::size_t k, OutputIterator output) const {
 
     // Create an empty list of points
@@ -907,18 +894,16 @@ private: // functions :
 public:
 
   /// \cond SKIP_IN_MANUAL
-  void dump_to_polylines (std::ostream& os) const
-  {
-    for (const Node& n : traverse<Orthtrees::Preorder_traversal>())
+  void dump_to_polylines(std::ostream& os) const {
+    for (const Node& n: traverse<Orthtrees::Preorder_traversal>())
       if (n.is_leaf())
       {
         Bbox box = bbox(n);
-        dump_box_to_polylines (box, os);
+        dump_box_to_polylines(box, os);
       }
   }
 
-  void dump_box_to_polylines (const Bbox_2& box, std::ostream& os) const
-  {
+  void dump_box_to_polylines(const Bbox_2& box, std::ostream& os) const {
     // dump in 3D for visualisation
     os << "5 "
        << box.xmin() << " " << box.ymin() << " 0 "
@@ -927,8 +912,8 @@ public:
        << box.xmax() << " " << box.ymin() << " 0 "
        << box.xmin() << " " << box.ymin() << " 0" << std::endl;
   }
-  void dump_box_to_polylines (const Bbox_3& box, std::ostream& os) const
-  {
+
+  void dump_box_to_polylines(const Bbox_3& box, std::ostream& os) const {
     // Back face
     os << "5 "
        << box.xmin() << " " << box.ymin() << " " << box.zmin() << " "
@@ -960,12 +945,12 @@ public:
        << box.xmax() << " " << box.ymax() << " " << box.zmax() << std::endl;
   }
 
-  friend std::ostream& operator<< (std::ostream& os, const Self& orthtree)
-  {
+  friend std::ostream& operator<<(std::ostream& os, const Self& orthtree) {
     // Create a range of nodes
     auto nodes = orthtree.traverse(Orthtrees::Preorder_traversal());
     // Iterate over the range
-    for (auto &n : nodes) {
+    for (auto& n: nodes)
+    {
       // Show the depth
       for (int i = 0; i < n.depth(); ++i)
         os << ". ";
diff --git a/Orthtree/include/CGAL/Orthtree/Node.h b/Orthtree/include/CGAL/Orthtree/Node.h
index a16d2fff192..42f439b5598 100644
--- a/Orthtree/include/CGAL/Orthtree/Node.h
+++ b/Orthtree/include/CGAL/Orthtree/Node.h
@@ -27,42 +27,20 @@
 namespace CGAL {
 
 /// \cond SKIP_IN_MANUAL
-namespace Orthtrees
-{
+namespace Orthtrees {
 
 // Non-documented, or testing purpose only
-struct Node_access
-{
+struct Node_access {
   template <typename Node, typename LC>
-  static Node create_node (Node parent, LC local_coordinates)
-  {
+  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(); }
+  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();
-    }
-  }
+  static void split(Node& node) { return node.split(); }
 
 };
 
@@ -79,9 +57,8 @@ struct Node_access
 
   \cgalModels `ConstRange`
  */
-template<typename Traits, typename PointRange, typename PointMap>
-class Orthtree<Traits, PointRange, PointMap>::Node
-{
+template <typename Traits, typename PointRange, typename PointMap>
+class Orthtree<Traits, PointRange, PointMap>::Node {
 
 public:
 
@@ -144,20 +121,25 @@ private:
   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;
+  // make Node trivially copiable
+//  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;
 
-    Data (Self parent)
-      : parent (parent), depth (0) { }
-  };
-
-  Data* m_data;
+  Point_range m_points;
+  Self* m_parent; // todo: use optional<reference_wrapper<Self>> instead of Self *
+  std::uint8_t m_depth;
+  Global_coordinates m_global_coordinates;
+  std::shared_ptr<Children> m_children;
 
 
   /// \cond SKIP_IN_MANUAL
@@ -173,8 +155,9 @@ private:
    * \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; }
+  Point_range& points() { return m_points; }
+
+  const Point_range& points() const { return m_points; }
 
   /// \name Construction
   /// @{
@@ -194,42 +177,38 @@ private:
     \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)) {
+  explicit Node(Self* parent, Local_coordinates local_coordinates)
+    : m_parent(parent) {
 
-    if (!parent.is_null()) {
-
-      m_data->depth = parent.m_data->depth + 1;
+    if (parent != nullptr) {
+      m_depth = parent->m_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];
+        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;
     }
-    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 deep_copy(Self parent = Node()) const {
 
     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>();
+    out.m_parent = m_parent;
+    out.m_points = m_points;
+    out.m_depth = m_depth;
+    out.m_global_coordinates = m_global_coordinates;
+
+    if (!is_leaf()) {
+      out.m_children = std::make_shared<Children>();
       for (int index = 0; index < Degree::value; index++)
-        (*out.m_data->children)[index] = (*this)[index].deep_copy(out);
+        (*out.m_children)[index] = (*m_children)[index].deep_copy(out);
     }
+
     return out;
   }
 
@@ -251,10 +230,10 @@ private:
 
     CGAL_precondition (is_leaf());
 
-    m_data->children = std::make_unique<Children>();
+    m_children = std::make_shared<Children>();
     for (int index = 0; index < Degree::value; index++) {
 
-      (*m_data->children)[index] = std::move(Self(*this, {Local_coordinates(index)}));
+      (*m_children)[index] = std::move(Self(this, {Local_coordinates(index)}));
 
     }
   }
@@ -267,7 +246,7 @@ private:
    */
   void unsplit() {
 
-    m_data->children.reset();
+    m_children.reset();
   }
 
   /// @}
@@ -279,10 +258,12 @@ public:
 
   /// \cond SKIP_IN_MANUAL
   // Default creates null node
-  Node() : m_data(nullptr) { }
+  // todo: default node is no longer null, but still isn't guaranteed to be valid
+  Node() = default;
 
   // Comparison operator
-  bool operator< (const Node& other) const { return m_data < other.m_data; }
+  // todo: where is this used
+  //bool operator<(const Node& other) const { return m_data < other.m_data; }
   /// \endcond
 
   /// \name Type & Location
@@ -291,36 +272,30 @@ public:
   /*!
     \brief returns `true` if the node is null, `false` otherwise.
   */
-  bool is_null() const { return (m_data == nullptr); }
+  //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();
+  bool is_root() const {
+    return m_parent == nullptr;
   }
 
   /*!
     \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);
+  bool is_leaf() const {
+    return (!m_children);
   }
 
   /*!
     \brief returns this node's depth.
     \pre `!is_null()`
   */
-  std::uint8_t depth() const
-  {
-    CGAL_precondition (!is_null());
-    return m_data->depth;
+  std::uint8_t depth() const {
+    return m_depth;
   }
 
   /*!
@@ -329,12 +304,9 @@ public:
   */
   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)
+    for (std::size_t i = 0; i < Dimension::value; ++i)
       result[i] = global_coordinates()[i] & 1;
 
     return result;
@@ -344,10 +316,8 @@ public:
     \brief returns this node's global coordinates.
     \pre `!is_null()`
   */
-  Global_coordinates global_coordinates() const
-  {
-    CGAL_precondition (!is_null());
-    return m_data->global_coordinates;
+  Global_coordinates global_coordinates() const {
+    return m_global_coordinates;
   }
 
 
@@ -356,18 +326,16 @@ public:
 
   /*!
     \brief returns this node's parent.
-    \pre `!is_null()`
-  */
-  Self parent() const
-  {
-    CGAL_precondition (!is_null());
-    return m_data->parent;
+    \pre `!is_root()`
+   */
+  const Self* parent() const {
+    CGAL_precondition (!is_root());
+    return m_parent;
   }
 
   /*!
     \brief returns the nth child of this node.
 
-    \pre `!is_null()`
     \pre `!is_leaf()`
     \pre `0 <= index && index < Degree::value`
 
@@ -418,13 +386,29 @@ public:
     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 {
+  Self& operator[](std::size_t index) {
 
-    CGAL_precondition (!is_null());
     CGAL_precondition (!is_leaf());
     CGAL_precondition (index < Degree::value);
 
-    return (*m_data->children)[index];
+    return (*m_children)[index];
+  }
+
+  /*!
+    \brief returns the nth child of this node.
+
+    \pre `!is_leaf()`
+    \pre `index < Degree::value`
+
+    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`.
+   */
+  const Self& operator[](std::size_t index) const {
+
+    CGAL_precondition (!is_leaf());
+    CGAL_precondition (index < Degree::value);
+
+    return (*m_children)[index];
   }
 
   /*!
@@ -478,9 +462,7 @@ public:
 
     \return the adjacent node if it exists, a null node otherwise.
   */
-  Self adjacent_node (Local_coordinates direction) const
-  {
-    CGAL_precondition(!is_null());
+  const Self* adjacent_node(Local_coordinates direction) const {
 
     // Direction:   LEFT  RIGHT  DOWN    UP  BACK FRONT
     // direction:    000    001   010   011   100   101
@@ -489,8 +471,7 @@ public:
     CGAL_precondition(direction.to_ulong() < Dimension::value * 2);
 
     // The root node has no adjacent nodes!
-    if (is_root())
-      return Self();
+    if (is_root()) return nullptr;
 
     // The least significant bit indicates the sign (which side of the node)
     bool sign = direction[0];
@@ -506,24 +487,22 @@ public:
 
     // 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];
+      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);
+    // Find the parent's neighbor in that direction, if it exists
+    const 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 (!adjacent_node_of_parent) return nullptr;
 
     // If the parent's adjacent node has no children, then it's this node's adjacent node
-    if (adjacent_node_of_parent.is_leaf())
+    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];
+    return &(*adjacent_node_of_parent)[local_coordinates().to_ulong() - offset];
 
   }
 
@@ -531,7 +510,21 @@ public:
     \brief equivalent to `adjacent_node()`, with an adjacency direction
     rather than a bitset.
    */
-  Self adjacent_node(Adjacency adjacency) const {
+  const Self* adjacent_node(Adjacency adjacency) const {
+    return adjacent_node(std::bitset<Dimension::value>(static_cast<int>(adjacency)));
+  }
+
+  /*!
+   * \brief equivalent to adjacent_node, except non-const
+   */
+  Self* adjacent_node(std::bitset<Dimension::value> direction) {
+    return const_cast<Self*>(const_cast<const Self*>(this)->adjacent_node(direction));
+  }
+
+  /*!
+   * \brief equivalent to adjacent_node, with a Direction rather than a bitset and non-const
+   */
+  Self* adjacent_node(Adjacency adjacency) {
     return adjacent_node(std::bitset<Dimension::value>(static_cast<int>(adjacency)));
   }
 
@@ -544,27 +537,27 @@ public:
     \brief checks whether the node is empty of points or not.
    */
   bool empty() const {
-    return m_data->points.empty();
+    return m_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()));
+    return std::size_t(std::distance(m_points.begin(), m_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(); }
+  const_iterator begin() const { return m_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(); }
+  const_iterator end() const { return m_points.end(); }
 
   /// @}
 
@@ -575,18 +568,23 @@ public:
   /*!
    * \brief compares the topology of this node to another node.
    *
-   * \todo
+   * \todo This seems out of date, the implementation I see compares for direct equality
    *
    * \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;
+  bool operator==(const Self& rhs) const {
+
+    // todo: This is a trivial implementation, maybe it can be set to =default in c++17?
+    return rhs.m_parent == m_parent &&
+           rhs.m_children == m_children &&
+           rhs.m_points == m_points &&
+           rhs.m_depth == m_depth &&
+           rhs.m_global_coordinates == m_global_coordinates;
   }
 
-  static bool is_topology_equal (const Self& a, const Self& b)
-  {
-    CGAL_assertion (!a.is_null() && !b.is_null());
+  // todo: this does what the documentation for operator== claims to do!
+  static bool is_topology_equal(const Self& a, const Self& b) {
 
     // If one node is a leaf, and the other isn't, they're not the same
     if (a.is_leaf() != b.is_leaf())
@@ -607,8 +605,7 @@ public:
     return (a.global_coordinates() == b.global_coordinates());
   }
 
-  friend std::ostream& operator<< (std::ostream& os, const Self& node)
-  {
+  friend std::ostream& operator<<(std::ostream& os, const Self& node) {
     return internal::print_orthtree_node(os, node);
   }
   /// \endcond
diff --git a/Orthtree/include/CGAL/Orthtree/Traversal_iterator.h b/Orthtree/include/CGAL/Orthtree/Traversal_iterator.h
index d8a594e4ef0..6670b480cd6 100644
--- a/Orthtree/include/CGAL/Orthtree/Traversal_iterator.h
+++ b/Orthtree/include/CGAL/Orthtree/Traversal_iterator.h
@@ -44,7 +44,7 @@ public:
    *
    * \todo
    */
-  typedef std::function<Value(Value)> Traversal_function;
+  typedef std::function<Value *(Value *)> Traversal_function;
 
   /// @}
 
@@ -58,7 +58,7 @@ public:
    *
    * \todo
    */
-  Traversal_iterator() : m_value(), m_next() {}
+  Traversal_iterator() : m_value(nullptr), m_next() {}
 
   /*!
    * \brief
@@ -68,7 +68,7 @@ public:
    * \param first
    * \param next
    */
-  Traversal_iterator(Value first, const Traversal_function &next) : m_value(first), m_next(next) {}
+  Traversal_iterator(Value *first, const Traversal_function &next) : m_value(first), m_next(next) {}
 
   /// @}
 
@@ -84,12 +84,12 @@ private:
   }
 
   Value &dereference() const {
-    return const_cast<Value&>(m_value);
+    return *m_value;
   }
 
 private:
 
-  Value m_value;
+  Value *m_value;
   Traversal_function m_next;
 };
 }
diff --git a/Orthtree/include/CGAL/Orthtree/Traversals.h b/Orthtree/include/CGAL/Orthtree/Traversals.h
index 885af4f9b4d..b7de300b3c9 100644
--- a/Orthtree/include/CGAL/Orthtree/Traversals.h
+++ b/Orthtree/include/CGAL/Orthtree/Traversals.h
@@ -24,7 +24,7 @@ namespace CGAL {
 
 /// \cond SKIP_IN_MANUAL
 // Forward declaration
-template<typename T, typename PR, typename PM>
+template <typename T, typename PR, typename PM>
 class Orthtree;
 /// \endcond
 
@@ -33,86 +33,85 @@ namespace Orthtrees {
 /// \cond SKIP_IN_MANUAL
 
 template <typename Node>
-Node next_sibling(Node n) {
+const Node* next_sibling(const Node* n) {
 
   // Passing null returns the first node
-  if (n.is_null())
-    return Node();
+  if (nullptr == n)
+    return nullptr;
 
   // If this node has no parent, it has no siblings
-  if (n.parent().is_null())
-    return Node();
+  if (n->is_root())
+    return nullptr;
 
   // Find out which child this is
-  std::size_t index = n.local_coordinates().to_ulong();
+  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();
+    return nullptr;
 
   // Otherwise, return the next child
-  return n.parent()[index + 1];
+  return &((*n->parent())[index + 1]);
 }
 
 template <typename Node>
-Node next_sibling_up(Node n) {
+const Node* next_sibling_up(const Node* n) {
 
-  if (n.is_null())
-    return Node();
+  if (!n || n->is_root()) return nullptr;
 
-  Node up = n.parent();
+  auto up = n->parent();
+  while (nullptr != up) {
 
-  while (!up.is_null()) {
-
-    if (!next_sibling(up).is_null())
+    if (nullptr != next_sibling(up))
       return next_sibling(up);
 
-    up = up.parent();
+    if (up->is_root()) return nullptr;
+
+    up = up->parent();
   }
 
-  return Node();
+  return nullptr;
 }
 
 template <typename Node>
-Node deepest_first_child(Node n) {
+const Node* deepest_first_child(const Node* n) {
 
-  if (n.is_null())
-    return Node();
+  if (!n)
+    return nullptr;
 
   // Find the deepest child on the left
-  Node first = n;
-  while (!first.is_leaf())
-    first = first[0];
+  auto 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) {
+const Node& first_child_at_depth(const Node* n, std::size_t depth) {
 
-  if (n.is_null())
-    return Node();
+  if (!n)
+    return nullptr;
 
-  std::stack<Node> todo;
+  std::stack<const Node*> todo;
   todo.push(n);
 
-  if (n.depth() == depth)
+  if (n->depth() == depth)
     return n;
 
-  while (!todo.empty())
-  {
-    Node node = todo.top();
+  while (!todo.empty()) {
+    const Node* node = todo.top();
     todo.pop();
 
-    if (node.depth() == depth)
+    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)]);
+    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();
+  return nullptr;
 }
 
 /// \endcond
@@ -128,25 +127,30 @@ Node first_child_at_depth(Node n, std::size_t depth) {
 struct Preorder_traversal {
 
   template <typename Node>
-  Node first(Node root) const {
+  const Node* first(const Node* root) const {
     return root;
   }
 
   template <typename Node>
-  Node next(Node n) const {
+  const Node* next(const Node* n) const {
 
-    if (n.is_leaf()) {
+    if (n->is_leaf()) {
 
-      Node next = next_sibling(n);
+      auto next = next_sibling(n);
+
+      if (nullptr == next) {
 
-      if (next.is_null())
         return next_sibling_up(n);
+      }
 
       return next;
 
+    } else {
+
+      // Return the first child of this node
+      return &(*n)[0];
     }
-    else // Return the first child of this node
-      return n[0];
+
   }
 };
 
@@ -161,18 +165,18 @@ struct Preorder_traversal {
 struct Postorder_traversal {
 
   template <typename Node>
-  Node first(Node root) const {
+  const Node* first(const Node* root) const {
 
     return deepest_first_child(root);
   }
 
   template <typename Node>
-  Node next(Node n) const {
+  const Node* next(const Node* n) const {
 
-    Node next = deepest_first_child(next_sibling(n));
+    auto next = deepest_first_child(next_sibling(n));
 
-    if (next.is_null())
-      next = n.parent();
+    if (!next)
+      next = n->parent();
 
     return next;
   }
@@ -189,17 +193,17 @@ struct Postorder_traversal {
 struct Leaves_traversal {
 
   template <typename Node>
-  Node first(Node root) const {
+  const Node* first(const Node* root) const {
 
     return deepest_first_child(root);
   }
 
   template <typename Node>
-  Node next(Node n) const {
+  const Node* next(const Node* n) const {
 
-    Node next = deepest_first_child(next_sibling(n));
+    auto next = deepest_first_child(next_sibling(n));
 
-    if (next.is_null())
+    if (!next)
       next = deepest_first_child(next_sibling_up(n));
 
     return next;
@@ -226,29 +230,28 @@ public:
   /*!
     constructs a `depth`-level traversal.
   */
-  Level_traversal (std::size_t depth) : depth(depth) { }
+  Level_traversal(std::size_t depth) : depth(depth) {}
 
   template <typename Node>
-  Node first(Node root) const {
+  const Node* first(const Node* root) const {
     return first_child_at_depth(root, depth);
   }
 
   template <typename Node>
-  Node next(Node n) const {
+  const Node* next(const Node* n) const {
+    // fixme: leftover from debugging?
     std::cerr << depth << " ";
-    Node next = next_sibling(n);
+    const Node* next = next_sibling(n);
 
-    if (next.is_null())
-    {
-      Node up = n;
-      do
-      {
+    if (!next) {
+      const Node* up = n;
+      do {
         up = next_sibling_up(up);
-        if (up.is_null())
-          return Node();
+        if (!up)
+          return nullptr;
+
         next = first_child_at_depth(up, depth);
-      }
-      while (next.is_null());
+      } while (!next);
     }
 
     return next;
diff --git a/Orthtree/test/Orthtree/test_node_adjacent.cpp b/Orthtree/test/Orthtree/test_node_adjacent.cpp
index 084d72dc517..1d4e55535d0 100644
--- a/Orthtree/test/Orthtree/test_node_adjacent.cpp
+++ b/Orthtree/test/Orthtree/test_node_adjacent.cpp
@@ -10,7 +10,7 @@ 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;
+  Octree;
 typedef Octree::Node Node;
 typedef Octree::Traits Traits;
 
@@ -42,33 +42,36 @@ 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.root().adjacent_node(0) == nullptr);
+  assert(octree.root().adjacent_node(1) == nullptr);
+  assert(octree.root().adjacent_node(2) == nullptr);
+  assert(octree.root().adjacent_node(3) == nullptr);
+  assert(octree.root().adjacent_node(4) == nullptr);
+  assert(octree.root().adjacent_node(5) == nullptr);
 
   // Left Top Front node should have siblings to the Right, Down, and Back
   auto left_top_back = octree.root()[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.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) == nullptr);
+  assert(left_top_back.adjacent_node(Traits::UP) == nullptr);
+  assert(left_top_back.adjacent_node(Traits::BACK) == nullptr);
 
   std::cout << octree.root()[Traits::LEFT_BOTTOM_BACK] << std::endl;
 
   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.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) != nullptr);
+  assert(right_top_back_of_left_bottom_back.adjacent_node(Traits::UP) != nullptr);
+  assert(right_top_back_of_left_bottom_back.adjacent_node(Traits::DOWN) != nullptr);
+  assert(right_top_back_of_left_bottom_back.adjacent_node(Traits::FRONT) != nullptr);
+
+  // A node at the back of the tree should have no neighbor to its back
+  assert(right_top_back_of_left_bottom_back.adjacent_node(Traits::BACK) == nullptr);
 
   return 0;
 }
diff --git a/Orthtree/test/Orthtree/test_octree_grade.cpp b/Orthtree/test/Orthtree/test_octree_grade.cpp
index a0a574bd9ef..686068028fc 100644
--- a/Orthtree/test/Orthtree/test_octree_grade.cpp
+++ b/Orthtree/test/Orthtree/test_octree_grade.cpp
@@ -27,10 +27,10 @@ std::size_t count_jumps(Octree &octree) {
 
       auto adjacent_node = node.adjacent_node(direction);
 
-      if (adjacent_node.is_null())
+      if (adjacent_node == nullptr)
         continue;
 
-      if ((node.depth() - adjacent_node.depth()) > 1)
+      if ((node.depth() - adjacent_node->depth()) > 1)
         jumps++;
     }
   }
diff --git a/Orthtree/test/Orthtree/test_octree_intersecting.cpp b/Orthtree/test/Orthtree/test_octree_intersecting.cpp
index 59ce353fd28..4d83929c2e0 100644
--- a/Orthtree/test/Orthtree/test_octree_intersecting.cpp
+++ b/Orthtree/test/Orthtree/test_octree_intersecting.cpp
@@ -46,14 +46,14 @@ int main(void) {
     auto query = Point{1, 1, 1};
 
     // Get a list of nodes intersected
-    std::vector<Octree::Node> nodes{};
+    std::vector<const Octree::Node *> nodes{};
     octree.intersected_nodes(query, std::back_inserter(nodes));
 
     // A point should only intersect one node
     assert(1 == nodes.size());
 
     // That node should be the node leaf that contains the point
-    assert(octree.locate(Point(1, 1, 1)) == nodes[0]);
+    assert(octree.locate(Point(1, 1, 1)) == *nodes[0]);
   }
 
   // Intersection with a sphere
@@ -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<const Octree::Node *> 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[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]);
   }
 
   // Intersection with a ray
@@ -81,19 +81,19 @@ 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<const Octree::Node *> 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[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]);
   }
 
   return EXIT_SUCCESS;
diff --git a/Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp b/Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp
index 8e87eaab152..bec700dbc46 100644
--- a/Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp
+++ b/Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp
@@ -20,7 +20,7 @@ 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;
+  Octree;
 
 typedef CGAL::Search_traits_3<Kernel> Kd_tree_traits;
 typedef CGAL::Orthogonal_k_neighbor_search<Kd_tree_traits> Kd_tree_search;
@@ -47,7 +47,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,7 +72,6 @@ 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();
@@ -109,9 +108,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 --> "
@@ -143,6 +142,7 @@ void kdtree_vs_octree(std::size_t dataset_size, std::size_t K) {
 
 int main(void) {
 
+  naive_vs_octree(21);
   naive_vs_octree(500);
   naive_vs_octree(1000);
   naive_vs_octree(10000);
diff --git a/Orthtree/test/Orthtree/test_octree_refine.cpp b/Orthtree/test/Orthtree/test_octree_refine.cpp
index ff517782179..82b0862e327 100644
--- a/Orthtree/test/Orthtree/test_octree_refine.cpp
+++ b/Orthtree/test/Orthtree/test_octree_refine.cpp
@@ -25,12 +25,11 @@ void test_1_point() {
   octree.refine(10, 1);
 
   // Check that the topology matches
-  Node single_node = CGAL::Orthtrees::Node_access::create_node(Node(), 0);
+  Node single_node = CGAL::Orthtrees::Node_access::create_node(static_cast<Node *>(nullptr), 0);
   CGAL::Orthtrees::Node_access::points(single_node)
     = CGAL::Orthtrees::Node_access::points(octree.root());
   assert(Node::is_topology_equal(single_node, octree.root()));
   assert(0 == octree.depth());
-  CGAL::Orthtrees::Node_access::free(single_node);
 }
 
 void test_2_points() {
@@ -45,12 +44,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);
+  Node other = CGAL::Orthtrees::Node_access::create_node(static_cast<Node *>(nullptr), 0);
   CGAL::Orthtrees::Node_access::split(other);
   assert(Node::is_topology_equal(other, octree.root()));
   assert(1 == octree.depth());
-  CGAL::Orthtrees::Node_access::free(other);
-
 }
 
 void test_4_points() {
@@ -66,13 +63,12 @@ void test_4_points() {
   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);
+  Node other = CGAL::Orthtrees::Node_access::create_node(static_cast<Node *>(nullptr), 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()));
   assert(2 == octree.depth());
-  CGAL::Orthtrees::Node_access::free(other);
 }
 
 int main(void) {