Clear And Efficient 3D Range Tree Implementation

Question

I'm working on this project where I have to search for objects in 3d space, and efficiency is a huge concern, I think Range Tree is perfect for what I'm trying to do, Interval Tree would also work but I'm not going to delete anything from the tree, once I add every object in 3D space, I'll only use the structure to do a search.

Here's how I'm going to use the structure:

Let say that I have an array(let's call it queryArr) of objects (~ 10,000 objects) each objects have 3 parameter (x,y,z) I have another very large array(let's call it totalArr) of objects ( > 5,000,000 objects).

What I'm trying to do here is given element of queryArr, find the most similar(or the same element in totalArr) In some cases there will be an object with the same parameters in totalArr, but in most cases, there won't be an object with same parameters.

So, I'll search all the values in between (x+10,y+10,z+10) and (x-10,y-10,z-10). If it doesn't yield any result, I'll multiply x,y,z by 2 and try again until it yields some results.

The simplest way of doing this is a naive search method, which will have complexity of O(N*M) (N = size of queryArr, M = sie of totalArr) but this method is incredibly slow and dumb.

I think Range Tree's is the way to go, but I have never implemented one myself and I don't quite understand how the range tree works on dimensions bigger than 2. So does anyone know a good implementation of the range trees? I believe if I can have a source code, I'll be able to understand how they really work.

By the way if you think, there is a better structure than Range Tree for this task, let me know I'm open to suggestions. (I have already considered kd-Trees and Interval trees)

Thanks,

Answer 1

I just read the wikipedia article. Lets see if I can write an n dimensional range tree. Because anything worth doing in 3 dimensions is worth doing in n.

So the basic part of an n-dimensional range tree is that it can be recursively defined in terms of lower dimensional range trees.

Some property classes to work with relatively generic value types. Specialize element_properties<T> to set the scalar type of whatever your n-dimensional value is, and specialize get<i>(T const&) to get the ith dimension of your n-dimensional value.

#include <memory>
#include <cstddef>
#include <vector>
#include <iostream>
#include <algorithm>
#include <string>
#include <sstream>

void Assert(bool test) {
  if (!test)
  {
    std::cout << "Assert failed" << std::endl;
    exit(-1);
  }
}
template<typename... Args>
struct Print {
  static void Do(Args... args) {}
};
template<typename Arg, typename... Tail>
struct Print<Arg, Tail...> {
  static void Do(Arg arg, Tail... args) {
      std::cout << arg;
      Print<Tail...>::Do(args...);
  }
};
template<typename... Args>
void Debug(Args... args) {
    std::cout << "DEBUG:[";
    Print<Args...>::Do(args...);
    std::cout << "]\n";
}

template<typename T>
struct element_properties {
  typedef typename T::value_type value_type;
};
template<>
struct element_properties<int> {
  typedef int value_type;
};
template<size_t d, typename T>
typename element_properties<T>::value_type get( T const & t );

template<size_t d>
typename element_properties<int>::value_type get( int i ) { return i; }

template<size_t d, typename U, typename A>
typename element_properties<std::vector<U,A>>::value_type get( std::vector<U,A> const& v) {
  return v[d];
}

template<typename T, size_t dim, typename Order = std::less< typename element_properties<T>::value_type> >
struct range_tree {
  typedef typename element_properties<T>::value_type value_type;
  struct sorter {
    bool operator()( T const& left, T const& right ) const {
      return Order()( get<dim-1>(left), get<dim-1>(right) );
    }
  };
  struct printer {
    std::string operator()( T const& t ) const {
      std::string retval = "[ ";
      retval += print_elements( t );
      retval += "]";
      return retval;
    }
    std::string print_elements( T const& t ) const {
      std::stringstream ss;
      typedef typename range_tree<T, dim-1, Order>::printer next_printer;
      ss << next_printer().print_elements(t);
      ss << get<dim-1>(t) << " ";
      return ss.str();
    }
  };
  template<typename Iterator>
  range_tree( Iterator begin, Iterator end ) {
    std::sort( begin, end, sorter() );
    root.reset( new tree_node( begin, end ) );
  }

  template<size_t n, typename Func>
  void walk(Func f) const {
      if (root) root->walk<n>(f);
  }
  template<size_t n, typename Func>
  void walk(Func f) {
      if (root) root->walk<n>(f);
  }
  struct tree_node {
    std::unique_ptr< range_tree<T, dim-1, Order> > subtree;
    T value;
    template<size_t n, typename Func>
    void walk(Func f) const {
      if (n==dim && !left && !right)
        f(value);
      if (left)
        left->walk<n>(f);
      if (right)
        right->walk<n>(f);
      if (subtree)
        subtree->walk<n>(f);
    }
    template<size_t n, typename Func>
    void walk(Func f) {
      if (n==dim && !left && !right)
        f(value);
      if (left)
        left->walk<n>(f);
      if (right)
        right->walk<n>(f);
      if (subtree)
        subtree->walk<n>(f);
    }
    void find_path( T const& t, std::vector< tree_node const* >& vec ) {
      vec.push_back(this);
      if ( sorter()(t, value) ) {
        if (left)
          left->find_path(t, vec);
      } else if (sorter()(value, t)) {
        if (right)
          right->find_path(t, vec);
      } else {
        // found it!
        return;
      }
    }
    std::vector< tree_node const* > range_search( T const& left, T const& right )
    {
      std::vector<tree_node const*> left_path;
      std::vector<tree_node const*> right_path;
      find_path( left, left_path );
      find_path( right, right_path );
      // erase common path:
      {
        auto it1 = left_path.begin();
        auto it2 = right_path.begin();
        for( ; it1 != left_path.end() && it2 != right_path.end(); ++it1, ++it2) {
          if (*it1 != *it2)
          {
            Debug( "Different: ", printer()( (*it1)->value ), ", ", printer()( (*it2)->value ) );
            break;
          }

          Debug( "Identical: ", printer()( (*it1)->value ), ", ", printer()( (*it2)->value ) );
        }
        // remove identical prefixes:
        if (it2 == right_path.end() && it2 != right_path.begin())
            --it2;
        if (it1 == left_path.end() && it1 != left_path.begin())
            --it1;
        right_path.erase( right_path.begin(), it2 );
        left_path.erase( left_path.begin(), it1 );
      }
      for (auto it = left_path.begin(); it != left_path.end(); ++it) {
        if (*it && (*it)->right) {
          Debug( "Has right child: ", printer()( (*it)->value ) );
          *it = (*it)->right.get();
          Debug( "It is: ", printer()( (*it)->value ) );
        } else {
          Debug( "Has no right child: ", printer()( (*it)->value ) );
          if ( sorter()( (*it)->value, left) || sorter()( right, (*it)->value) ) {
            Debug( printer()( (*it)->value ), "<", printer()( left ), " so erased" );
            *it = 0;
          }
        }
      }
      for (auto it = right_path.begin(); it != right_path.end(); ++it) {
        if (*it && (*it)->left) {
          Debug( "Has left child: ", printer()( (*it)->value ) );
          *it = (*it)->left.get();
          Debug( "It is: ", printer()( (*it)->value ) );
        } else {
          Debug( "Has no left child: ", printer()( (*it)->value ) );
          if ( sorter()( (*it)->value, left) || sorter()( right, (*it)->value) ) {
            Debug( printer()( right ), "<", printer()( (*it)->value ), " so erased" );
            *it = 0;
          }
        }
      }
      left_path.insert( left_path.end(), right_path.begin(), right_path.end() );
      // remove duds and duplicates:
      auto highwater = std::remove_if( left_path.begin(), left_path.end(), []( tree_node const* n) { return n==0; } );
      std::sort( left_path.begin(), highwater );
      left_path.erase( std::unique( left_path.begin(), highwater ), left_path.end() );
      return left_path;
    }

    std::unique_ptr<tree_node> left;
    std::unique_ptr<tree_node> right;
    // rounds down:
    template<typename Iterator>
    static Iterator middle( Iterator begin, Iterator end ) {
      return (end-begin-1)/2 + begin ;
    }
    template<typename Iterator>
    tree_node( Iterator begin, Iterator end ):value(*middle(begin,end)) {
      Debug( "Inserted ", get<dim-1>(value), " at level ", dim );
      Iterator mid = middle(begin,end);
      Assert( begin != end );
      if (begin +1 != end) { // not a leaf
        Debug( "Not a leaf at level ", dim );
        ++mid; // so *mid was the last element in the left sub tree 
        Assert(mid!=begin);
        Assert(mid!=end);
        left.reset( new tree_node( begin, mid ) );
        right.reset( new tree_node( mid, end ) );
      } else {
        Debug( "Leaf at level ", dim );
      }
      if (dim > 0) {
        subtree.reset( new range_tree<T, dim-1, Order>( begin, end ) );
      }
    }
  };
  std::unique_ptr<tree_node> root;
};
// makes the code above a tad easier:
template<typename T, typename Order >
struct range_tree< T, 0, Order > {
  typedef typename element_properties<T>::value_type value_type;
  struct printer { template<typename Unused>std::string print_elements(Unused const&) {return std::string();} };
  range_tree(...) {};
  struct tree_node {}; // maybe some stub functions in here
  template<size_t n, typename Func>
  void walk(Func f) {}
};

int main() {
  typedef std::vector<int> vector_type;
  std::vector<vector_type> test;
  test.push_back( vector_type{5,2} );
  test.push_back( vector_type{2,3} );
  range_tree< vector_type, 2 > tree( test.begin(), test.end() );
  std::cout << "Walking dim 2:";
  auto print_node = [](vector_type const& v){ std::cout << "(" << v[0] << "," << v[1] << ")"; };
  tree.walk<2>( print_node );
  std::cout << "\nWalking dim 1:";
  tree.walk<1>( print_node );
  std::cout << "\n";

  std::cout << "Range search from {3,3} to {10,10}\n";
  auto nodes = tree.root->range_search( vector_type{3,3}, vector_type{10,10} );
  for (auto it = nodes.begin(); it != nodes.end(); ++it)
  {
    (*it)->walk<2>( print_node );
  }
}

which is pretty damn close to an n-dimensional range tree. The 0 dimension tree naturally contains nothing.

Basic facilities to search (in one dimension at a time) have been added now. You can manually do the recursions into lower dimensions, or it up so that the range_search always returns level 1 tree_node*s.