5
\$\begingroup\$

I've implemented a rather straightforward octree, but I don't feel that good about the coding style or some of the implementation. In particular, I use some nasty void pointers and extra news because I couldn't figure unions out, and I'd like to have that replaced (I'm thinking maybe templating can solve it, but any feedback is appreciated).

There are a few methods I'm not very interested in that are in the interface, but haven't been implemented. Feel free to ignore them, as they shouldn't have a substantive effect on the rest of the code.

Lastly, I'm not looking for a complete overhaul of the structure of the octree - I'm going to implement some other versions at a later date (i.e. pointerless duals, a few other fun ones), so please review this structure, where each node has a pointer to 8 children.

Octree.h

/*
    file - octree.h

    Templated implementation of an octree for a cpu

 */


#ifndef OCTREE_CPU_H
#define OCTREE_CPU_H

#include <algorithm>
#include <array>
#include <cstddef>
#include <iterator>
#include <limits>
#include <utility>
#include <type_traits>

#include "boundingbox.h"

using Point = std::array<double, 3>;

template <typename InputIterator, class PointExtractor, 
          size_t max_per_node = 16, size_t max_depth = 100>
class Octree {
 public:
  using tree_type = Octree<InputIterator, PointExtractor, max_per_node, max_depth>;

  Octree();

  Octree(InputIterator begin, InputIterator end);

  Octree(InputIterator begin, InputIterator end, PointExtractor f);

  Octree(const tree_type& rhs);

  template <size_t max_per_node_>
  Octree(const Octree<InputIterator, PointExtractor, max_per_node_, max_depth>& rhs);

  template <size_t max_depth_>
  Octree(const Octree<InputIterator, PointExtractor, max_per_node, max_depth_>& rhs);

  template <size_t max_per_node_, size_t max_depth_>
  Octree(const Octree<InputIterator, PointExtractor, max_per_node_, max_depth_>& rhs);

  Octree(tree_type&& rhs);

  void swap(tree_type& rhs);

  template <typename OutputIterator>
  bool search(const BoundingBox& box, OutputIterator& it) const;

  tree_type& operator=(tree_type rhs);

  tree_type& operator=(tree_type&& rhs);

  ~Octree();

  size_t size() const;

 private:  
  class Node;

  enum NodeContents {
    LEAF = 1,
    MAX_DEPTH_LEAF = 2,
    INTERNAL = 4
  };

  struct LeafNodeValues {
    std::array<std::pair<InputIterator, Point>, max_per_node> values_;
    size_t size_;
  };

  using childNodeArray = std::array<Node*, 8>;
  using maxItemNode = std::vector<std::pair<InputIterator, Point>>;

  class Node {
   public:    
    Node(const std::vector<std::pair<InputIterator, Point>>& input_values);

    Node(const std::vector<std::pair<InputIterator, Point>>& input_values, 
         const BoundingBox& box,
         size_t current_depth);

    ~Node();

    template <typename OutputIterator>
    bool search(const BoundingBox& box, OutputIterator& it) const;

   private:
    void* value_;
    BoundingBox extrema_;
    NodeContents tag_;

    void init_max_depth_leaf(const std::vector<std::pair<InputIterator, Point>>& input_values);

    void init_leaf(const std::vector<std::pair<InputIterator, Point>>& input_values);

    void init_internal(
        const std::vector<std::pair<InputIterator, Point>>& input_values,
        size_t current_depth);

    unsigned int getOctantIndex(const Point& p) const;

    struct InnerIterator {
      using wrapped_type = typename std::vector<std::pair<InputIterator, Point>>::const_iterator;
      wrapped_type it_;

      InnerIterator(wrapped_type it) : it_(it) {}

      Point operator*() const {
        return std::get<1>(*it_);
      }

      InnerIterator& operator++() {
        ++it_;
        return *this;
      }

      InnerIterator operator++(int) {
        InnerIterator other = *this;
        ++it_;
        return other;
      }

      bool operator==(const InnerIterator& rhs) const {
        return it_ == rhs.it_;
      }

      bool operator!=(const InnerIterator& rhs) const {
        return !operator==(rhs);
      }
    };
  };

  PointExtractor functor_;
  Node* head_;
  size_t size_;
};

// convenience macros to avoid typing so much
#define OCTREE Octree<InputIterator, PointExtractor, max_per_node, max_depth>
#define OCTREE_TEMPLATE typename InputIterator, class PointExtractor, size_t max_per_node, size_t max_depth

template <OCTREE_TEMPLATE>
OCTREE::Octree(): functor_(PointExtractor()), head_(nullptr), size_(0) {}

template <OCTREE_TEMPLATE>
OCTREE::Octree(InputIterator begin, InputIterator end)
  : Octree(begin, end, PointExtractor()) { }

template <OCTREE_TEMPLATE>
OCTREE::Octree(InputIterator begin, InputIterator end, PointExtractor f)
    : functor_(f), head_(nullptr), size_(0) {

  std::vector<std::pair<InputIterator, Point>> v;
  v.reserve(std::distance(begin, end));

  for (auto it = begin; it != end; ++it) {
    v.push_back(std::pair<InputIterator, Point>(it, functor_(*it)));
  }

  size_ = v.size();
  head_ = new Node(v);
}

template <OCTREE_TEMPLATE>
OCTREE::Octree(OCTREE::tree_type&& rhs) 
  : functor_(rhs.functor), head_(rhs.head_), size_(rhs.size_) { }

template <OCTREE_TEMPLATE>
void OCTREE::swap(OCTREE::tree_type& rhs) {
  std::swap(head_, rhs.head_);
  std::swap(functor_, rhs.functor_);
  std::swap(size_, rhs.size_);
}

template <OCTREE_TEMPLATE>
template <typename OutputIterator>
bool OCTREE::search(const BoundingBox& box, OutputIterator& it) const {
  return head_->search(box, it);
}

template <OCTREE_TEMPLATE>
typename OCTREE::tree_type& OCTREE::operator=(typename OCTREE::tree_type rhs) {
  swap(rhs);
  return *this;
}

template <OCTREE_TEMPLATE>
typename OCTREE::tree_type& OCTREE::operator=(typename OCTREE::tree_type&& rhs) {
  swap(rhs);
  return *this;
}

template <OCTREE_TEMPLATE>
OCTREE::~Octree() {
  delete head_;
}

template <OCTREE_TEMPLATE>
size_t OCTREE::size() const {
  return size_;
}

template <OCTREE_TEMPLATE>
OCTREE::Node::Node(const std::vector<std::pair<InputIterator, Point>>& input_values)
  : Node(input_values, 
         BoundingBox(InnerIterator(input_values.begin()), InnerIterator(input_values.end())),
         0) { }

template <OCTREE_TEMPLATE>
OCTREE::Node::Node(
    const std::vector<std::pair<InputIterator, Point>>& input_values, 
    const BoundingBox& box,
    size_t current_depth) : value_(nullptr), extrema_(box)  {
  if (current_depth > max_depth) {
    init_max_depth_leaf(input_values);
  } else if (input_values.size() <= max_per_node) {
    init_leaf(input_values);
  } else {
    init_internal(input_values, current_depth);
  }
}

template <OCTREE_TEMPLATE>
OCTREE::Node::~Node() {
  if (tag_ == NodeContents::INTERNAL) {
    childNodeArray* children = static_cast<childNodeArray*>(value_);
    for (size_t i = 0; i < 8; ++i) {
      delete children[0][i];
      children[0][i] = nullptr;
    }
    delete children;
  } else if (tag_ == NodeContents::LEAF) {
    delete static_cast<LeafNodeValues*>(value_);
  } else if (tag_ == NodeContents::MAX_DEPTH_LEAF) {
    delete static_cast<maxItemNode*>(value_);
  }

  value_ = nullptr;
}

template <OCTREE_TEMPLATE>
template <typename OutputIterator>
bool OCTREE::Node::search(const BoundingBox& p, OutputIterator& it) const {
  bool success = false;
  if (tag_ == NodeContents::INTERNAL) {
    childNodeArray& children = *static_cast<childNodeArray*>(value_);
    for (auto child : children) {
      if (child) {
        success = child->search(p, it) || success;
      }
    }
  } else if (tag_ == NodeContents::LEAF) {
    LeafNodeValues& children = *static_cast<LeafNodeValues*>(value_);
    for (size_t i = 0; i < children.size_; ++i) {
      Point& point = std::get<1>(children.values_[i]);
      if (p.contains(point)) {
        *it = std::get<0>(children.values_[i]);
        ++it;
        success = true;
      }
    }
  } else if (tag_ == NodeContents::MAX_DEPTH_LEAF) {
    maxItemNode& children = *static_cast<maxItemNode*>(value_);
    for (auto child : children) {
      Point& point = std::get<1>(child);
      if (p.contains(point)) {
        *it = std::get<0>(child);
        ++it;
        success = true;
      }
    }
  }
  return success;
}

template <OCTREE_TEMPLATE>
void OCTREE::Node::init_max_depth_leaf(
    const std::vector<std::pair<InputIterator, Point>>& input_values) {  
  value_ = new std::vector<std::pair<InputIterator, Point>>(input_values);
  tag_ = NodeContents::MAX_DEPTH_LEAF;
}

template <OCTREE_TEMPLATE>
void OCTREE::Node::init_leaf(
    const std::vector<std::pair<InputIterator, Point>>& input_values)  {
  std::array<std::pair<InputIterator, Point>, max_per_node> a;
  std::copy(input_values.begin(), input_values.end(), a.begin());
  value_ = new LeafNodeValues{a, input_values.size()};
  tag_ = NodeContents::LEAF;
}

template <OCTREE_TEMPLATE>
void OCTREE::Node::init_internal(
    const std::vector<std::pair<InputIterator, Point>>& input_values,
    size_t current_depth)  {
  std::array<std::vector<std::pair<InputIterator, Point>>, 8> childVectors;
  std::array<BoundingBox, 8> boxes = extrema_.partition();
  std::array<Node*, 8> children;

  for (unsigned child = 0; child < 8; ++child) {
    std::vector<std::pair<InputIterator, Point>>& childVector = childVectors[child];
    childVector.reserve(input_values.size() / 8);

    std::copy_if(
      input_values.begin(), 
      input_values.end(), 
      std::back_inserter(childVector),
      [&boxes, child](const std::pair<InputIterator, Point>& element) -> bool {
        Point p = std::get<1>(element);
        return boxes[child].contains(p);
      }
    );

    children[child] = childVector.empty()
        ? nullptr
        : new Node(childVector, boxes[child], ++current_depth);
  }

  value_ = new std::array<Node*, 8>(children);
  tag_ = NodeContents::INTERNAL;
}

template <OCTREE_TEMPLATE>
unsigned int OCTREE::Node::getOctantIndex(const Point& p) const {
  // children are ordered left to right, front to back, bottom to top.

  double xmid = (extrema_.xhi - extrema_.xlo) / 2.;
  double ymid = (extrema_.yhi - extrema_.ylo) / 2.;
  double zmid = (extrema_.zhi - extrema_.zlo) / 2.;
  bool left = p[0] < xmid && p[0] >= extrema_.xlo;
  bool front = p[1] < ymid && p[1] >= extrema_.ylo;
  bool bottom = p[2] < zmid && p[2] >= extrema_.zlo;

  if (bottom && left && front) {
    return 0;
  } else if (bottom && !left && front) {
    return 1;
  } else if (bottom && left && !front) {
    return 2;
  } else if (bottom && !left && !front) {
    return 3;
  } else if (!bottom && left && front) {
    return 4;
  } else if (!bottom && !left && front) {
    return 5;
  } else if (!bottom && left && !front) {
    return 6;
  } else {
    return 7;
  }
}

#endif // DEFINED OCTREE_CPU_H

boundingbox.h

#ifndef BOUNDINGBOX_H
#define BOUNDINGBOX_H

#include <array>
#include <limits>
#include <cmath>
#include <iostream>

// An axis aligned bounding boxed (AABB)
struct BoundingBox {
  BoundingBox() = default;
  template <typename InputIterator>
  BoundingBox(InputIterator begin, InputIterator end);
  BoundingBox(std::initializer_list<double> l);
  BoundingBox(BoundingBox&& rhs);
  BoundingBox(const BoundingBox&) = default;

  BoundingBox& operator=(BoundingBox&) = default;
  BoundingBox& operator=(BoundingBox&& rhs);
  BoundingBox& operator=(const BoundingBox& rhs);

  ~BoundingBox() = default;

  bool contains(const BoundingBox& other) const;
  bool contains(const std::array<double, 3>& point) const;

  bool overlap(const BoundingBox& other, BoundingBox* out) const;

  std::array<BoundingBox, 8> partition() const;

  bool operator==(const BoundingBox& rhs) const;

  bool operator!=(const BoundingBox& rhs) const;

  friend
  std::ostream& operator<<(std::ostream& stream, const BoundingBox& rhs);

  double xhi, xlo, yhi, ylo, zhi, zlo;
};

const BoundingBox initial = BoundingBox{
    std::numeric_limits<double>::min(), std::numeric_limits<double>::max(),
    std::numeric_limits<double>::min(), std::numeric_limits<double>::max(),
    std::numeric_limits<double>::min(), std::numeric_limits<double>::max()
};

const BoundingBox invalid = BoundingBox{
    std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), 
    std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), 
    std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()
};

template <typename InputIterator>
BoundingBox::BoundingBox(InputIterator begin, InputIterator end) : BoundingBox(initial) {
  for (; begin != end; ++begin) {
    const std::array<double, 3>& point = *begin;

    if (point[0] < xlo) {
      xlo = point[0];
    } else if (point[0] > xhi) {
      xhi = point[0];
    }

    if (point[1] < ylo) {
      ylo = point[1];
    } else if (point[1] > yhi) {
      yhi = point[1];
    }

    if (point[2] < zlo) {
      zlo = point[2];
    } else if (point[2] > zhi) {
      zhi = point[2];
    }
  }
}

#endif // defined BOUNDINGBOX_H

boundingbox.cc

#include "boundingbox.h"

#include <algorithm>
#include <array>
#include <limits>
#include <cmath>
#include <iostream>
#include <utility>

using std::array;
using std::initializer_list;

BoundingBox::BoundingBox(BoundingBox&& rhs) {
  xhi = std::move(rhs.xhi);
  xlo = std::move(rhs.xlo);
  yhi = std::move(rhs.yhi);
  ylo = std::move(rhs.ylo);
  zhi = std::move(rhs.zhi);
  zlo = std::move(rhs.zlo);
}

BoundingBox::BoundingBox(initializer_list<double> l) {
  std::copy(l.begin(), l.end(), &xhi);
}

BoundingBox& BoundingBox::operator=(const BoundingBox& rhs) {
  std::copy(&rhs.xhi, &rhs.xhi + 6, &xhi);
  return *this;
}

BoundingBox& BoundingBox::operator=(BoundingBox&& rhs) {
  xhi = std::move(rhs.xhi);
  xlo = std::move(rhs.xlo);
  yhi = std::move(rhs.yhi);
  ylo = std::move(rhs.ylo);
  zhi = std::move(rhs.zhi);
  zlo = std::move(rhs.zlo);
  return *this;
}

bool BoundingBox::contains(const BoundingBox& other) const {
  return xlo <= other.xlo && xhi >= other.xhi &&
         ylo <= other.ylo && yhi >= other.yhi &&
         zlo <= other.zlo && zhi >= other.zhi;
}

bool BoundingBox::contains(const std::array<double, 3>& point) const {
  return xlo <= point[0] && xhi > point[0] &&
         ylo <= point[1] && yhi > point[1] &&
         zlo <= point[2] && zhi > point[2];
}

bool BoundingBox::overlap(const BoundingBox& other, BoundingBox* out) const {
  // trivial cases
  if (contains(other)) {
    *out = other;
    return true;
  } else if (other.contains(*this)) {
    *out = *this;
    return true;
  } 

  // Check if there is no intersection
  if (xhi < other.xlo || xlo > other.xhi ||
      yhi < other.ylo || ylo > other.yhi ||
      zhi < other.zlo || zlo > other.zhi) {
    *out = invalid;
    return false;
  }

  // Actually calculate the bounds
  double upperX = std::min(xhi, other.xhi);
  double upperY = std::min(yhi, other.yhi);
  double upperZ = std::min(zhi, other.zhi);

  double lowerX = std::max(xlo, other.xlo);
  double lowerY = std::max(ylo, other.ylo);
  double lowerZ = std::max(zlo, other.zlo);

  *out = BoundingBox{upperX, lowerX, upperY, lowerY, upperZ, lowerZ};
  return true;
}

array<BoundingBox, 8> BoundingBox::partition() const {
  double xmid = (xhi - xlo) / 2.;
  double ymid = (yhi - ylo) / 2.;
  double zmid = (zhi - zlo) / 2.;

  std::array<BoundingBox, 8> ret{{
    BoundingBox{xmid, xlo, ymid, ylo, zmid, zlo}, // bottom left front
    BoundingBox{xhi, xmid, ymid, ylo, zmid, zlo}, // bottom right front
    BoundingBox{xmid, xlo, yhi, ymid, zmid, zlo}, // bottom left back
    BoundingBox{xhi, xmid, yhi, ymid, zmid, zlo}, // bottom right back
    BoundingBox{xmid, xlo, ymid, ylo, zhi, zmid}, // top left front
    BoundingBox{xhi, xmid, ymid, ylo, zhi, zmid}, // top right front
    BoundingBox{xmid, xlo, yhi, ymid, zhi, zmid}, // top left back
    BoundingBox{xhi, xmid, yhi, ymid, zhi, zmid}  // top right back
  }};
  return ret;
}

bool BoundingBox::operator==(const BoundingBox& rhs) const {
  // They're all equal, or they're all NaNs
  return (xhi == rhs.xhi && xlo == rhs.xlo &&
          yhi == rhs.yhi && ylo == rhs.ylo &&
          zhi == rhs.zhi && zlo == rhs.zlo) ||
         (std::isnan(xhi) && std::isnan(rhs.xhi) &&
          std::isnan(xlo) && std::isnan(rhs.xlo) &&
          std::isnan(yhi) && std::isnan(rhs.yhi) &&
          std::isnan(ylo) && std::isnan(rhs.ylo) &&
          std::isnan(zhi) && std::isnan(rhs.zhi) &&
          std::isnan(zlo) && std::isnan(rhs.zlo));
}

bool BoundingBox::operator!=(const BoundingBox& rhs) const {
  return !operator==(rhs);
}

std::ostream& operator<<(std::ostream& stream, const BoundingBox& rhs) {
  stream << "{" << rhs.xhi << ", " << rhs.xlo << ", "
                << rhs.yhi << ", " << rhs.ylo << ", "
                << rhs.zhi << ", " << rhs.zlo << ", "
         << "}";
  return stream;
}
\$\endgroup\$
5
\$\begingroup\$

I'm looking only at the BoundingBox class. There are a lot of interesting cleanups you can apply to it:

Ditch the constructors

BoundingBox is a Plain Old Data (POD) type. It only holds a handful of doubles, so the compiler will supply all the necessary default constructors and operators, including:

  • Copy constructor
  • Assignment operator
  • No-op destructor

Move operator and move-assign are pointless. You cannot move a double, it fits in a machine register, so the ones you've implemented basically do the same as a copy constructor. Your copy/assignment code is essentially doing unsafe memcopies and might even be less efficient than the defaults the compiler would generate had you omitted them.

But we still got these two:

BoundingBox(InputIterator begin, InputIterator end);
BoundingBox(std::initializer_list<double> l);

The initializer list constructor is cool, it lets you write code like this:

BoundingBox bb = { 1.0, 2.0, 3.0,  4.0, 5.0, 6.0 };

But it turns out that you don't even need to define one. If you instead remove all the other constructors and leave all data public, as it is now, you can already initialize a struct instance with that syntax.

But to get the aggregate-style initialization working we need to remove all user defined constructors, so the begin/end range constructor has to go. I would consider this a fair trade, given how much it would simplify the code, so instead, consider turning that constructor into a helper function:

template <typename InputIterator>
BoundingBox makeBoundingBox(InputIterator begin, InputIterator end);

BoundingBox is cheap to copy, so that function should be just as fast as the constructor, assuming it isn't inlined (which is likely, since it is a template).

Types are your friends

So don't be afraid to make new friends (i.e. create other helper types) ;)

This member data:

double xhi, xlo, yhi, ylo, zhi, zlo;

Should be a type, like a Point3D:

struct Point3D {
    double x;
    double y;
    double z;
};

I can see that you have a quasi-Point3D already by using a std::array of 3 elements, but I'd advise making it a struct in its own right, so it can have methods and also the nicer v.x, v.y v.z syntax over just v[0], v[1], v[2].

Miscellaneous

  • The stream output operator doesn't have to be a friend function. The struct has no private data. You only use friend when you want to allow access to the private data of the type in question. Just make it a free-standing function.

  • The comparison operator is way too big. Assuming we went with the Point3D idea, it could be rewriting like so:

    bool BoundingBox::operator == (const BoundingBox& rhs) const {
        const bool allEqual = (mins == rhs.mins && maxs == rhs.maxs);
        if (allEqual) return true;
    
        const bool allNaN = (mins.isNaN() && rhs.mins.isNaN() && 
                             maxs.isNaN() && rhs.maxs.isNaN());
        if (allNaN) return true;
    
        return false;
    }
    

    A couple notes on the above:

    • I've used the names mins and maxs for our bounding box extents. That's a very common notation for this particular data structure, so I'd suggest sticking to that.

    • Also note that comparing floating point numbers with == is not very reliable. They are prone to floating point rounding errors, so you might end up with two very close numbers that for all purposes would be equal but end up comparing not equal. Take a look at Stack Overflow on how to compare floats for equality within a given epsilon.


Summing up, this is how your axis-aligned bounding-box could look like after the suggested changes:

struct Point3D {
    double x;
    double y;
    double z;

    bool isNaN() const;
    bool operator == (const Point3D& other) const;
};

struct BoundingBox {
    Point3D mins;
    Point3D maxs;

    bool contains(const BoundingBox& other) const;
    bool contains(const Point3D& point) const;

    bool overlap(const BoundingBox& other, BoundingBox* out) const;
    std::array<BoundingBox, 8> partition() const;

    bool operator == (const BoundingBox& rhs) const;
    bool operator != (const BoundingBox& rhs) const;
};

template <typename InputIterator>
BoundingBox makeBoundingBox(InputIterator begin, InputIterator end);
std::ostream& operator << (std::ostream& stream, const BoundingBox& rhs);
\$\endgroup\$
3
\$\begingroup\$

By no means a complete review (there is a lot of code), just few notes on what jumped out at the first reading.

  • getOctantCode computes the result in a very long way. Consider

        return (!bottom << 2) | (!left << 1) | (!front);
    
  • I don't think that special-casing "trivial cases" simplifies BoundingBox::overlap in any respect.

  • I am not sure that there is any guarantee of how xhi, xlo, yhi, ylo, zhi, zlo; are going to be lay out in memory. The code

        std::copy(&rhs.xhi, &rhs.xhi + 6, &xhi);
    

seems to assume that they are sequential, and it is safe to access them as an array. Maybe it is all right and mandated by Standard, but a reviewer is definitely stumbled.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.