Allow other BVH implementations to be used.

dmnsn_bvh is now a generic API, which could potentially support
octrees, etc, in addition to PR-trees.

10 files changed, 600 insertions, 472 deletions

diff --git a/libdimension/Makefile.am b/libdimension/Makefile.am
index aa5921e..a8f06fc 100644
--- a/libdimension/Makefile.am
+++ b/libdimension/Makefile.am
@@ -1,5 +1,5 @@
 ###########################################################################
-## Copyright (C) 2009-2011 Tavian Barnes <tavianator@tavianator.com>     ##
+## Copyright (C) 2009-2012 Tavian Barnes <tavianator@tavianator.com>     ##
 ##                                                                       ##
 ## This file is part of The Dimension Build Suite.                       ##
 ##                                                                       ##
@@ -63,6 +63,8 @@ nobase_include_HEADERS = dimension.h                                           \
 lib_LTLIBRARIES = libdimension.la
 
 libdimension_la_SOURCES = $(nobase_include_HEADERS)                            \
+                          bvh.c                                                \
+                          bvh.h                                                \
                           camera.c                                             \
                           canvas.c                                             \
                           canvas_pigment.c                                     \
diff --git a/libdimension/bench/prtree.c b/libdimension/bench/prtree.c
index b168890..953c7b4 100644
--- a/libdimension/bench/prtree.c
+++ b/libdimension/bench/prtree.c
@@ -1,5 +1,5 @@
 /*************************************************************************
- * Copyright (C) 2009-2011 Tavian Barnes <tavianator@tavianator.com>     *
+ * Copyright (C) 2009-2012 Tavian Barnes <tavianator@tavianator.com>     *
  *                                                                       *
  * This file is part of The Dimension Benchmark Suite.                   *
  *                                                                       *
@@ -17,6 +17,7 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>. *
  *************************************************************************/
 
+#include "../bvh.c"
 #include "../prtree.c"
 #include "../threads.c"
 #include "../future.c"
@@ -88,13 +89,13 @@ main(void)
     dmnsn_array_push(objects, &object);
   }
 
-  dmnsn_prtree *tree;
+  dmnsn_bvh *bvh;
   sandglass_bench_noprecache(&sandglass, {
-    tree = dmnsn_new_prtree(objects);
+    bvh = dmnsn_new_bvh(objects, DMNSN_BVH_PRTREE);
   });
-  printf("dmnsn_new_prtree(): %ld\n", sandglass.grains);
+  printf("dmnsn_new_bvh(DMNSN_BVH_PRTREE): %ld\n", sandglass.grains);
 
-  /* dmnsn_prtree_intersection() */
+  /* dmnsn_bvh_intersection() */
   dmnsn_line ray = dmnsn_new_line(
     dmnsn_new_vector( 1.0,  1.0, -2.0),
     dmnsn_new_vector(-0.5, -0.5,  1.0)
@@ -102,23 +103,23 @@ main(void)
   dmnsn_intersection intersection;
 
   sandglass_bench_fine(&sandglass, {
-    dmnsn_prtree_intersection(tree, ray, &intersection, true);
+    dmnsn_bvh_intersection(bvh, ray, &intersection, true);
   });
-  printf("dmnsn_prtree_intersection(): %ld\n", sandglass.grains);
+  printf("dmnsn_bvh_intersection(): %ld\n", sandglass.grains);
 
   sandglass_bench_fine(&sandglass, {
-    dmnsn_prtree_intersection(tree, ray, &intersection, false);
+    dmnsn_bvh_intersection(bvh, ray, &intersection, false);
   });
-  printf("dmnsn_prtree_intersection(nocache): %ld\n", sandglass.grains);
+  printf("dmnsn_bvh_intersection(nocache): %ld\n", sandglass.grains);
 
-  /* dmnsn_prtree_inside() */
+  /* dmnsn_bvh_inside() */
   sandglass_bench_fine(&sandglass, {
-    dmnsn_prtree_inside(tree, dmnsn_zero);
+    dmnsn_bvh_inside(bvh, dmnsn_zero);
   });
-  printf("dmnsn_prtree_inside(): %ld\n", sandglass.grains);
+  printf("dmnsn_bvh_inside(): %ld\n", sandglass.grains);
 
   /* Cleanup */
-  dmnsn_delete_prtree(tree);
+  dmnsn_delete_bvh(bvh);
   DMNSN_ARRAY_FOREACH (dmnsn_object **, object, objects) {
     dmnsn_delete_object(*object);
   }
diff --git a/libdimension/bvh.c b/libdimension/bvh.c
new file mode 100644
index 0000000..a038d11
--- /dev/null
+++ b/libdimension/bvh.c
@@ -0,0 +1,448 @@
+/*************************************************************************
+ * Copyright (C) 2012 Tavian Barnes <tavianator@tavianator.com>          *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * BVH implementation.  These are the hottest code paths in libdimension.
+ */
+
+#include "dimension-internal.h"
+
+/** Implementation for DMNSN_BVH_NONE: just stick all objects in one node. */
+static dmnsn_bvh_node *
+dmnsn_new_stupid_bvh(const dmnsn_array *objects)
+{
+  dmnsn_bvh_node *root = dmnsn_new_bvh_node(dmnsn_array_size(objects));
+
+  DMNSN_ARRAY_FOREACH (dmnsn_object **, object, objects) {
+    dmnsn_bvh_node *leaf = dmnsn_new_bvh_leaf_node(*object);
+    dmnsn_bvh_node_add(root, leaf);
+  }
+
+  return root;
+}
+
+/* Implementation of opaque dmnsn_bvh type. */
+struct dmnsn_bvh {
+  dmnsn_array *unbounded;          /**< The unbounded objects. */
+  dmnsn_array *bounded;            /**< The BVH of the bounded objects. */
+  size_t id;                       /**< A unique ID for this BVH. */
+};
+
+/** A flat BVH node for storing in an array for fast pre-order traversal. */
+typedef struct dmnsn_flat_bvh_node {
+  dmnsn_bounding_box bounding_box; /* The bounding box of this node. */
+  dmnsn_object *object;            /* The referenced object, for leaf nodes. */
+  ptrdiff_t skip;                  /* Displacement to the next sibling. */
+} dmnsn_flat_bvh_node;
+
+/** Add an object or its children, if any, to an array. */
+static void
+dmnsn_split_add_object(dmnsn_array *objects, const dmnsn_object *object)
+{
+  if (object->split_children) {
+    DMNSN_ARRAY_FOREACH (const dmnsn_object **, child, object->children) {
+      dmnsn_split_add_object(objects, *child);
+    }
+  } else {
+    dmnsn_array_push(objects, &object);
+  }
+}
+
+/** Split unions to create the input for the BVH. */
+static dmnsn_array *
+dmnsn_split_objects(const dmnsn_array *objects)
+{
+  dmnsn_array *split = dmnsn_new_array(sizeof(dmnsn_object *));
+  DMNSN_ARRAY_FOREACH (const dmnsn_object **, object, objects) {
+    dmnsn_split_add_object(split, *object);
+  }
+  return split;
+}
+
+/** Split unbounded objects into a new array. */
+static dmnsn_array *
+dmnsn_split_unbounded(dmnsn_array *objects)
+{
+  dmnsn_array *unbounded = dmnsn_new_array(sizeof(dmnsn_object *));
+
+  dmnsn_object **array = dmnsn_array_first(objects);
+  size_t i, skip;
+  for (i = 0, skip = 0; i < dmnsn_array_size(objects); ++i) {
+    if (dmnsn_bounding_box_is_infinite(array[i]->bounding_box)) {
+      dmnsn_array_push(unbounded, &array[i]);
+      ++skip;
+    } else {
+      array[i - skip] = array[i];
+    }
+  }
+  dmnsn_array_resize(objects, i - skip);
+
+  return unbounded;
+}
+
+/** Recursively flatten a BVH into an array of flat nodes. */
+static void
+dmnsn_flatten_bvh_recursive(dmnsn_bvh_node *node, dmnsn_array *flat)
+{
+  size_t currenti = dmnsn_array_size(flat);
+  dmnsn_array_resize(flat, currenti + 1);
+  dmnsn_flat_bvh_node *flatnode = dmnsn_array_at(flat, currenti);
+
+  flatnode->bounding_box = node->bounding_box;
+  flatnode->object       = node->object;
+
+  for (size_t i = 0; i < node->nchildren && node->children[i]; ++i) {
+    dmnsn_flatten_bvh_recursive(node->children[i], flat);
+  }
+
+  /* Array could have been realloc()'d somewhere else above */
+  flatnode = dmnsn_array_at(flat, currenti);
+  flatnode->skip = dmnsn_array_size(flat) - currenti;
+}
+
+/** Flatten a BVH into an array of flat nodes. */
+static dmnsn_array *
+dmnsn_flatten_bvh(dmnsn_bvh_node *root)
+{
+  dmnsn_array *flat = dmnsn_new_array(sizeof(dmnsn_flat_bvh_node));
+  if (root) {
+    dmnsn_flatten_bvh_recursive(root, flat);
+  }
+  return flat;
+}
+
+static size_t dmnsn_bvh_seq = 0;
+static pthread_mutex_t dmnsn_bvh_seq_mutex = PTHREAD_MUTEX_INITIALIZER;
+
+dmnsn_bvh *dmnsn_new_bvh(const dmnsn_array *objects, dmnsn_bvh_kind kind)
+{
+  dmnsn_bvh *bvh = dmnsn_malloc(sizeof(dmnsn_bvh));
+
+  dmnsn_array *bounded = dmnsn_split_objects(objects);
+  bvh->unbounded = dmnsn_split_unbounded(bounded);
+
+  dmnsn_bvh_node *root = NULL;
+  if (dmnsn_array_size(bounded) > 0) {
+    switch (kind) {
+    case DMNSN_BVH_NONE:
+      root = dmnsn_new_stupid_bvh(bounded);
+      break;
+    case DMNSN_BVH_PRTREE:
+      root = dmnsn_new_prtree(bounded);
+      break;
+    default:
+      dmnsn_unreachable("Invalid BVH kind.");
+    }
+  }
+  bvh->bounded = dmnsn_flatten_bvh(root);
+
+  dmnsn_delete_bvh_node(root);
+  dmnsn_delete_array(bounded);
+
+  dmnsn_lock_mutex(&dmnsn_bvh_seq_mutex);
+    bvh->id = dmnsn_bvh_seq++;
+  dmnsn_unlock_mutex(&dmnsn_bvh_seq_mutex);
+
+  return bvh;
+}
+
+void
+dmnsn_delete_bvh(dmnsn_bvh *bvh)
+{
+  if (bvh) {
+    dmnsn_delete_array(bvh->bounded);
+    dmnsn_delete_array(bvh->unbounded);
+    dmnsn_free(bvh);
+  }
+}
+
+/** A line with pre-calculated reciprocals to avoid divisions. */
+typedef struct dmnsn_optimized_line {
+  dmnsn_vector x0;    /**< The origin of the line. */
+  dmnsn_vector n_inv; /**< The inverse of each component of the line's slope .*/
+} dmnsn_optimized_line;
+
+/** Precompute inverses for faster ray-box intersection tests. */
+static inline dmnsn_optimized_line
+dmnsn_optimize_line(dmnsn_line line)
+{
+  dmnsn_optimized_line optline = {
+    .x0    = line.x0,
+    .n_inv = dmnsn_new_vector(1.0/line.n.x, 1.0/line.n.y, 1.0/line.n.z)
+  };
+  return optline;
+}
+
+/** Ray-AABB intersection test, by the slab method.  Highly optimized. */
+static inline bool
+dmnsn_ray_box_intersection(dmnsn_optimized_line optline,
+                           dmnsn_bounding_box box, double t)
+{
+  /*
+   * This is actually correct, even though it appears not to handle edge cases
+   * (line.n.{x,y,z} == 0).  It works because the infinities that result from
+   * dividing by zero will still behave correctly in the comparisons.  Lines
+   * which are parallel to an axis and outside the box will have tmin == inf
+   * or tmax == -inf, while lines inside the box will have tmin and tmax
+   * unchanged.
+   */
+
+  double tx1 = (box.min.x - optline.x0.x)*optline.n_inv.x;
+  double tx2 = (box.max.x - optline.x0.x)*optline.n_inv.x;
+
+  double tmin = dmnsn_min(tx1, tx2);
+  double tmax = dmnsn_max(tx1, tx2);
+
+  double ty1 = (box.min.y - optline.x0.y)*optline.n_inv.y;
+  double ty2 = (box.max.y - optline.x0.y)*optline.n_inv.y;
+
+  tmin = dmnsn_max(tmin, dmnsn_min(ty1, ty2));
+  tmax = dmnsn_min(tmax, dmnsn_max(ty1, ty2));
+
+  double tz1 = (box.min.z - optline.x0.z)*optline.n_inv.z;
+  double tz2 = (box.max.z - optline.x0.z)*optline.n_inv.z;
+
+  tmin = dmnsn_max(tmin, dmnsn_min(tz1, tz2));
+  tmax = dmnsn_min(tmax, dmnsn_max(tz1, tz2));
+
+  return tmax >= dmnsn_max(0.0, tmin) && tmin < t;
+}
+
+/** The number of intersections to cache. */
+#define DMNSN_INTERSECTION_CACHE_SIZE 32
+
+/** An array of cached intersections. */
+typedef struct dmnsn_intersection_cache {
+  size_t i;
+  dmnsn_object *objects[DMNSN_INTERSECTION_CACHE_SIZE];
+} dmnsn_intersection_cache;
+
+/** The thread-specific intersection cache. */
+static pthread_key_t dmnsn_intersection_caches;
+/** Initialize the thread-specific pointer exactly once. */
+static pthread_once_t dmnsn_intersection_caches_once = PTHREAD_ONCE_INIT;
+
+static void
+dmnsn_delete_intersection_caches(void *caches)
+{
+  dmnsn_delete_array(caches);
+}
+
+static void
+dmnsn_initialize_intersection_caches(void)
+{
+  dmnsn_key_create(&dmnsn_intersection_caches,
+                   dmnsn_delete_intersection_caches);
+}
+
+static dmnsn_array *
+dmnsn_get_intersection_caches(void)
+{
+  dmnsn_once(&dmnsn_intersection_caches_once,
+             dmnsn_initialize_intersection_caches);
+  return pthread_getspecific(dmnsn_intersection_caches);
+}
+
+/** Needed because pthreads doesn't destroy data from the main thread unless
+    it exits with pthread_exit(). */
+DMNSN_DESTRUCTOR static void
+dmnsn_delete_main_intersection_caches(void)
+{
+  dmnsn_delete_intersection_caches(dmnsn_get_intersection_caches());
+  dmnsn_key_delete(dmnsn_intersection_caches);
+}
+
+static dmnsn_intersection_cache *
+dmnsn_get_intersection_cache(size_t id)
+{
+  dmnsn_array *caches = dmnsn_get_intersection_caches();
+  if (!caches) {
+    caches = dmnsn_new_array(sizeof(dmnsn_intersection_cache));
+    dmnsn_setspecific(dmnsn_intersection_caches, caches);
+  }
+
+  while (dmnsn_array_size(caches) <= id) {
+    dmnsn_array_resize(caches, dmnsn_array_size(caches) + 1);
+    dmnsn_intersection_cache *cache = dmnsn_array_last(caches);
+    cache->i = 0;
+    for (size_t i = 0; i < DMNSN_INTERSECTION_CACHE_SIZE; ++i) {
+      cache->objects[i] = NULL;
+    }
+  }
+
+  return dmnsn_array_at(caches, id);
+}
+
+/** Test for a closer object intersection than we've found so far. */
+static inline bool
+dmnsn_closer_intersection(dmnsn_object *object, dmnsn_line ray,
+                          dmnsn_intersection *intersection, double *t)
+{
+  dmnsn_intersection local_intersection;
+  if (dmnsn_object_intersection(object, ray, &local_intersection)) {
+    if (local_intersection.t < *t) {
+      *intersection = local_intersection;
+      *t = local_intersection.t;
+      return true;
+    }
+  }
+  return false;
+}
+
+DMNSN_HOT bool
+dmnsn_bvh_intersection(const dmnsn_bvh *bvh, dmnsn_line ray,
+                       dmnsn_intersection *intersection, bool reset)
+{
+  double t = INFINITY;
+
+  /* Search the unbounded objects */
+  DMNSN_ARRAY_FOREACH (dmnsn_object **, object, bvh->unbounded) {
+    dmnsn_closer_intersection(*object, ray, intersection, &t);
+  }
+
+  /* Precalculate 1.0/ray.n.{x,y,z} to save time in intersection tests */
+  dmnsn_optimized_line optline = dmnsn_optimize_line(ray);
+
+  /* Search the intersection cache */
+  dmnsn_intersection_cache *cache = dmnsn_get_intersection_cache(bvh->id);
+  if (dmnsn_unlikely(reset)) {
+    cache->i = 0;
+  }
+  dmnsn_object *cached = NULL, *found = NULL;
+  if (dmnsn_likely(cache->i < DMNSN_INTERSECTION_CACHE_SIZE)) {
+    cached = cache->objects[cache->i];
+  }
+  if (cached && dmnsn_ray_box_intersection(optline, cached->bounding_box, t)) {
+    if (dmnsn_closer_intersection(cached, ray, intersection, &t)) {
+      found = cached;
+    }
+  }
+
+  /* Search the bounded objects */
+  dmnsn_flat_bvh_node *node = dmnsn_array_first(bvh->bounded);
+  dmnsn_flat_bvh_node *last = dmnsn_array_last(bvh->bounded);
+  while (node <= last) {
+    if (dmnsn_ray_box_intersection(optline, node->bounding_box, t)) {
+      if (node->object && node->object != cached) {
+        if (dmnsn_closer_intersection(node->object, ray, intersection, &t)) {
+          found = node->object;
+        }
+      }
+      ++node;
+    } else {
+      node += node->skip;
+    }
+  }
+
+  /* Update the cache */
+  if (dmnsn_likely(cache->i < DMNSN_INTERSECTION_CACHE_SIZE)) {
+    cache->objects[cache->i] = found;
+    ++cache->i;
+  }
+
+  return !isinf(t);
+}
+
+DMNSN_HOT bool
+dmnsn_bvh_inside(const dmnsn_bvh *bvh, dmnsn_vector point)
+{
+  /* Search the unbounded objects */
+  DMNSN_ARRAY_FOREACH (dmnsn_object **, object, bvh->unbounded) {
+    if (dmnsn_object_inside(*object, point))
+      return true;
+  }
+
+  /* Search the bounded objects */
+  dmnsn_flat_bvh_node *node = dmnsn_array_first(bvh->bounded);
+  dmnsn_flat_bvh_node *last = dmnsn_array_last(bvh->bounded);
+  while (node <= last) {
+    if (dmnsn_bounding_box_contains(node->bounding_box, point)) {
+      if (node->object && dmnsn_object_inside(node->object, point)) {
+        return true;
+      }
+      ++node;
+    } else {
+      node += node->skip;
+    }
+  }
+
+  return false;
+}
+
+dmnsn_bounding_box
+dmnsn_bvh_bounding_box(const dmnsn_bvh *bvh)
+{
+  if (dmnsn_array_size(bvh->unbounded) > 0) {
+    return dmnsn_infinite_bounding_box();
+  } else if (dmnsn_array_size(bvh->bounded) > 0) {
+    dmnsn_flat_bvh_node *root = dmnsn_array_first(bvh->bounded);
+    return root->bounding_box;
+  } else {
+    return dmnsn_zero_bounding_box();
+  }
+}
+
+dmnsn_bvh_node *
+dmnsn_new_bvh_node(size_t max_children)
+{
+  dmnsn_bvh_node *node = dmnsn_malloc(sizeof(dmnsn_bvh_node)
+                                      + max_children*sizeof(dmnsn_bvh_node *));
+  node->bounding_box = dmnsn_zero_bounding_box();
+  node->object = NULL;
+  node->nchildren = 0;
+  node->max_children = max_children;
+  return node;
+}
+
+dmnsn_bvh_node *
+dmnsn_new_bvh_leaf_node(dmnsn_object *object)
+{
+  dmnsn_bvh_node *node = dmnsn_malloc(sizeof(dmnsn_bvh_node));
+  node->bounding_box = object->bounding_box;
+  node->object = object;
+  node->nchildren = 0;
+  node->max_children = 0;
+  return node;
+}
+
+void
+dmnsn_delete_bvh_node(dmnsn_bvh_node *node)
+{
+  if (node) {
+    for (size_t i = 0; i < node->nchildren; ++i) {
+      dmnsn_delete_bvh_node(node->children[i]);
+    }
+    dmnsn_free(node);
+  }
+}
+
+void
+dmnsn_bvh_node_add(dmnsn_bvh_node *parent, dmnsn_bvh_node *child)
+{
+  dmnsn_assert(parent->nchildren < parent->max_children,
+               "Too many BVH children inserted.");
+
+  parent->bounding_box.min = dmnsn_vector_min(parent->bounding_box.min,
+                                              child->bounding_box.min);
+  parent->bounding_box.max = dmnsn_vector_max(parent->bounding_box.max,
+                                              child->bounding_box.max);
+  parent->children[parent->nchildren++] = child;
+}
diff --git a/libdimension/bvh.h b/libdimension/bvh.h
new file mode 100644
index 0000000..2c3a8a4
--- /dev/null
+++ b/libdimension/bvh.h
@@ -0,0 +1,76 @@
+/*************************************************************************
+ * Copyright (C) 2012 Tavian Barnes <tavianator@tavianator.com>          *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file.
+ * Bounding volume hierarchy.  This generic interface allows different data
+ * structures to be represented in the same way, thus sharing code for their
+ * traversal algorithm.
+ */
+
+#include <stdbool.h>
+
+/** A bounding volume hierarchy. */
+typedef struct dmnsn_bvh dmnsn_bvh;
+
+/** Available BVH implementations. */
+typedef enum dmnsn_bvh_kind {
+  DMNSN_BVH_NONE,
+  DMNSN_BVH_PRTREE,
+} dmnsn_bvh_kind;
+
+/** Create a BVH. */
+DMNSN_INTERNAL dmnsn_bvh *dmnsn_new_bvh(const dmnsn_array *objects,
+                                        dmnsn_bvh_kind kind);
+/** Delete a BVH. */
+DMNSN_INTERNAL void dmnsn_delete_bvh(dmnsn_bvh *bvh);
+
+/** Find the closest ray-object intersection in the tree. */
+DMNSN_INTERNAL bool dmnsn_bvh_intersection(const dmnsn_bvh *bvh,
+                                           dmnsn_line ray,
+                                           dmnsn_intersection *intersection,
+                                           bool reset);
+/** Determine whether a point is inside any object in the tree. */
+DMNSN_INTERNAL bool dmnsn_bvh_inside(const dmnsn_bvh *bvh, dmnsn_vector point);
+/** Return the bounding box of the whole hierarchy. */
+DMNSN_INTERNAL dmnsn_bounding_box dmnsn_bvh_bounding_box(const dmnsn_bvh *bvh);
+
+/** A non-flat BVH representation, used by BVH implementations. */
+typedef struct dmnsn_bvh_node {
+  dmnsn_bounding_box bounding_box;   /** The bounding box of this node. */
+  dmnsn_object *object;              /** The object, for leaf nodes. */
+  int data;                          /** Extra field for implementation use. */
+  size_t nchildren;                  /** How many children this node has. */
+  size_t max_children;               /** Maximum number of children. */
+  struct dmnsn_bvh_node *children[]; /** Flexible array of children. */
+} dmnsn_bvh_node;
+
+/** Create a BVH node. */
+DMNSN_INTERNAL dmnsn_bvh_node *dmnsn_new_bvh_node(size_t max_children);
+
+/** Create a BVH leaf node. */
+DMNSN_INTERNAL dmnsn_bvh_node *dmnsn_new_bvh_leaf_node(dmnsn_object *object);
+
+/** Delete a BVH node. */
+DMNSN_INTERNAL void dmnsn_delete_bvh_node(dmnsn_bvh_node *node);
+
+/** Add a child to a BVH node. */
+DMNSN_INTERNAL void dmnsn_bvh_node_add(dmnsn_bvh_node *parent,
+                                       dmnsn_bvh_node *child);
diff --git a/libdimension/csg.c b/libdimension/csg.c
index d65a300..2595f7d 100644
--- a/libdimension/csg.c
+++ b/libdimension/csg.c
@@ -1,5 +1,5 @@
 /*************************************************************************
- * Copyright (C) 2010-2011 Tavian Barnes <tavianator@tavianator.com>     *
+ * Copyright (C) 2010-2012 Tavian Barnes <tavianator@tavianator.com>     *
  *                                                                       *
  * This file is part of The Dimension Library.                           *
  *                                                                       *
@@ -36,16 +36,16 @@ dmnsn_csg_union_intersection_fn(const dmnsn_object *csg,
                                 dmnsn_line line,
                                 dmnsn_intersection *intersection)
 {
-  dmnsn_prtree *prtree = csg->ptr;
-  return dmnsn_prtree_intersection(prtree, line, intersection, true);
+  dmnsn_bvh *bvh = csg->ptr;
+  return dmnsn_bvh_intersection(bvh, line, intersection, true);
 }
 
 /** CSG union inside callback. */
 static bool
 dmnsn_csg_union_inside_fn(const dmnsn_object *csg, dmnsn_vector point)
 {
-  dmnsn_prtree *prtree = csg->ptr;
-  return dmnsn_prtree_inside(prtree, point);
+  dmnsn_bvh *bvh = csg->ptr;
+  return dmnsn_bvh_inside(bvh, point);
 }
 
 /** CSG union initialization callback. */
@@ -54,16 +54,16 @@ dmnsn_csg_union_initialize_fn(dmnsn_object *csg)
 {
   csg->trans = dmnsn_identity_matrix();
 
-  dmnsn_prtree *prtree = dmnsn_new_prtree(csg->children);
-  csg->ptr = prtree;
-  csg->bounding_box = prtree->bounding_box;
+  dmnsn_bvh *bvh = dmnsn_new_bvh(csg->children, DMNSN_BVH_PRTREE);
+  csg->ptr = bvh;
+  csg->bounding_box = dmnsn_bvh_bounding_box(bvh);
 }
 
 /** CSG union destruction callback. */
 static void
 dmnsn_csg_union_free_fn(void *ptr)
 {
-  dmnsn_delete_prtree(ptr);
+  dmnsn_delete_bvh(ptr);
 }
 
 /* Bulk-load a union */
diff --git a/libdimension/dimension-internal.h b/libdimension/dimension-internal.h
index eae585d..3f71fb1 100644
--- a/libdimension/dimension-internal.h
+++ b/libdimension/dimension-internal.h
@@ -1,5 +1,5 @@
 /*************************************************************************
- * Copyright (C) 2009-2011 Tavian Barnes <tavianator@tavianator.com>     *
+ * Copyright (C) 2009-2012 Tavian Barnes <tavianator@tavianator.com>     *
  *                                                                       *
  * This file is part of The Dimension Library.                           *
  *                                                                       *
@@ -35,6 +35,7 @@
 #include "platform.h"
 #include "future-impl.h"
 #include "threads.h"
+#include "bvh.h"
 #include "prtree.h"
 #include "rgba16.h"
 
diff --git a/libdimension/prtree.c b/libdimension/prtree.c
index 2451f6b..ba26d30 100644
--- a/libdimension/prtree.c
+++ b/libdimension/prtree.c
@@ -1,5 +1,5 @@
 /*************************************************************************
- * Copyright (C) 2010-2011 Tavian Barnes <tavianator@tavianator.com>     *
+ * Copyright (C) 2010-2012 Tavian Barnes <tavianator@tavianator.com>     *
  *                                                                       *
  * This file is part of The Dimension Library.                           *
  *                                                                       *
@@ -20,12 +20,10 @@
 
 /**
  * @file
- * Priority R-tree implementation.  These are the hottest code paths in
- * libdimension.
+ * Priority R-tree implementation.
  */
 
 #include "dimension-internal.h"
-#include <pthread.h>
 #include <stdlib.h>
 
 /** Number of children per PR-node. */
@@ -33,13 +31,6 @@
 /** Number of priority leaves per pseudo-PR-node (must be 2*ndimensions). */
 #define DMNSN_PSEUDO_B 6
 
-/** A flat node for storing in an array for fast pre-order traversal. */
-typedef struct dmnsn_flat_prnode {
-  dmnsn_bounding_box bounding_box;
-  dmnsn_object *object;
-  size_t skip;
-} dmnsn_flat_prnode;
-
 /** The side of the split that a node ended up on. */
 typedef enum dmnsn_prnode_location {
   DMNSN_PRTREE_LEAF, /**< Priority leaf. */
@@ -47,50 +38,15 @@ typedef enum dmnsn_prnode_location {
   DMNSN_PRTREE_RIGHT /**< Right child. */
 } dmnsn_prnode_location;
 
-/** Pseudo PR-tree node. */
-typedef struct dmnsn_prnode {
-  dmnsn_bounding_box bounding_box;
-
-  dmnsn_object *object;
-  struct dmnsn_prnode *children[DMNSN_PRTREE_B];
-
-  dmnsn_prnode_location location;
-} dmnsn_prnode;
-
 /** Construct an empty PR-node. */
-static inline dmnsn_prnode *
+static inline dmnsn_bvh_node *
 dmnsn_new_prnode(void)
 {
-  dmnsn_prnode *node = dmnsn_malloc(sizeof(dmnsn_prnode));
-  node->bounding_box = dmnsn_zero_bounding_box();
-  node->object       = NULL;
-  node->location     = DMNSN_PRTREE_LEFT; /* Mustn't be _LEAF */
-  for (size_t i = 0; i < DMNSN_PRTREE_B; ++i) {
-    node->children[i] = NULL;
-  }
+  dmnsn_bvh_node *node = dmnsn_new_bvh_node(DMNSN_PRTREE_B);
+  node->data = DMNSN_PRTREE_LEFT; /* Mustn't be _LEAF */
   return node;
 }
 
-/** Free a non-flat PR-tree. */
-static void
-dmnsn_delete_prnode(dmnsn_prnode *node)
-{
-  if (node) {
-    for (size_t i = 0; i < DMNSN_PRTREE_B && node->children[i]; ++i) {
-      dmnsn_delete_prnode(node->children[i]);
-    }
-    dmnsn_free(node);
-  }
-}
-
-/** Expand a node to contain the bounding box \p box. */
-static void
-dmnsn_prnode_swallow(dmnsn_prnode *node, dmnsn_bounding_box box)
-{
-  node->bounding_box.min = dmnsn_vector_min(node->bounding_box.min, box.min);
-  node->bounding_box.max = dmnsn_vector_max(node->bounding_box.max, box.max);
-}
-
 /** Comparator types. */
 enum {
   DMNSN_XMIN,
@@ -103,7 +59,7 @@ enum {
 
 /** Get a coordinate of the bounding box of a node. */
 static inline double
-dmnsn_get_coordinate(const dmnsn_prnode * const *node, int comparator)
+dmnsn_get_coordinate(const dmnsn_bvh_node * const *node, int comparator)
 {
   switch (comparator) {
   case DMNSN_XMIN:
@@ -191,23 +147,22 @@ dmnsn_add_priority_leaves(dmnsn_array *sorted_leaves[DMNSN_PSEUDO_B],
                           dmnsn_array *new_leaves)
 {
   for (size_t i = 0; i < DMNSN_PSEUDO_B; ++i) {
-    dmnsn_prnode *leaf = NULL;
-    dmnsn_prnode **leaves = dmnsn_array_first(sorted_leaves[i]);
-    for (size_t j = 0, count = 0, size = dmnsn_array_size(sorted_leaves[i]);
-         j < size && count < DMNSN_PRTREE_B;
+    dmnsn_bvh_node *leaf = NULL;
+    dmnsn_bvh_node **leaves = dmnsn_array_first(sorted_leaves[i]);
+    for (size_t j = 0, size = dmnsn_array_size(sorted_leaves[i]);
+         j < size && (!leaf || leaf->nchildren < DMNSN_PRTREE_B);
          ++j)
     {
       /* Skip all the previously found extreme nodes */
-      if (leaves[j]->location == DMNSN_PRTREE_LEAF) {
+      if (leaves[j]->data == DMNSN_PRTREE_LEAF) {
         continue;
       }
 
       if (!leaf) {
         leaf = dmnsn_new_prnode();
       }
-      leaves[j]->location = DMNSN_PRTREE_LEAF;
-      leaf->children[count++] = leaves[j];
-      dmnsn_prnode_swallow(leaf, leaves[j]->bounding_box);
+      leaves[j]->data = DMNSN_PRTREE_LEAF;
+      dmnsn_bvh_node_add(leaf, leaves[j]);
     }
 
     if (leaf) {
@@ -225,10 +180,10 @@ dmnsn_split_sorted_leaves(dmnsn_array *sorted_leaves[DMNSN_PSEUDO_B],
                           size_t i)
 {
   /* Get rid of the extreme nodes */
-  dmnsn_prnode **leaves = dmnsn_array_first(sorted_leaves[i]);
+  dmnsn_bvh_node **leaves = dmnsn_array_first(sorted_leaves[i]);
   size_t j, skip;
   for (j = 0, skip = 0; j < dmnsn_array_size(sorted_leaves[i]); ++j) {
-    if (leaves[j]->location == DMNSN_PRTREE_LEAF) {
+    if (leaves[j]->data == DMNSN_PRTREE_LEAF) {
       ++skip;
     } else {
       leaves[j - skip] = leaves[j];
@@ -242,24 +197,24 @@ dmnsn_split_sorted_leaves(dmnsn_array *sorted_leaves[DMNSN_PSEUDO_B],
 
   /* Split the appropriate list and mark the left and right child nodes */
   right_sorted_leaves[i] = dmnsn_array_split(sorted_leaves[i]);
-  DMNSN_ARRAY_FOREACH (dmnsn_prnode **, node, sorted_leaves[i]) {
-    (*node)->location = DMNSN_PRTREE_LEFT;
+  DMNSN_ARRAY_FOREACH (dmnsn_bvh_node **, node, sorted_leaves[i]) {
+    (*node)->data = DMNSN_PRTREE_LEFT;
   }
-  DMNSN_ARRAY_FOREACH (dmnsn_prnode **, node, right_sorted_leaves[i]) {
-    (*node)->location = DMNSN_PRTREE_RIGHT;
+  DMNSN_ARRAY_FOREACH (dmnsn_bvh_node **, node, right_sorted_leaves[i]) {
+    (*node)->data = DMNSN_PRTREE_RIGHT;
   }
 
   /* Split the rest of the lists */
   for (size_t j = 0; j < DMNSN_PSEUDO_B; ++j) {
     if (j != i) {
-      right_sorted_leaves[j] = dmnsn_new_array(sizeof(dmnsn_prnode *));
+      right_sorted_leaves[j] = dmnsn_new_array(sizeof(dmnsn_bvh_node *));
 
-      dmnsn_prnode **leaves = dmnsn_array_first(sorted_leaves[j]);
+      dmnsn_bvh_node **leaves = dmnsn_array_first(sorted_leaves[j]);
       size_t k, skip;
       for (k = 0, skip = 0; k < dmnsn_array_size(sorted_leaves[j]); ++k) {
-        if (leaves[k]->location == DMNSN_PRTREE_LEAF) {
+        if (leaves[k]->data == DMNSN_PRTREE_LEAF) {
           ++skip;
-        } else if (leaves[k]->location == DMNSN_PRTREE_RIGHT) {
+        } else if (leaves[k]->data == DMNSN_PRTREE_RIGHT) {
           dmnsn_array_push(right_sorted_leaves[j], &leaves[k]);
           ++skip;
         } else {
@@ -306,7 +261,7 @@ dmnsn_priority_leaves(const dmnsn_array *leaves)
     dmnsn_array_sort(sorted_leaves[i], dmnsn_comparators[i]);
   }
 
-  dmnsn_array *new_leaves = dmnsn_new_array(sizeof(dmnsn_prnode *));
+  dmnsn_array *new_leaves = dmnsn_new_array(sizeof(dmnsn_bvh_node *));
   dmnsn_priority_leaves_recursive(sorted_leaves, new_leaves, 0);
 
   for (size_t i = 0; i < DMNSN_PSEUDO_B; ++i) {
@@ -316,20 +271,18 @@ dmnsn_priority_leaves(const dmnsn_array *leaves)
   return new_leaves;
 }
 
-/** Construct a non-flat PR-tree. */
-static dmnsn_prnode *
-dmnsn_make_prtree(const dmnsn_array *objects)
+dmnsn_bvh_node *
+dmnsn_new_prtree(const dmnsn_array *objects)
 {
   if (dmnsn_array_size(objects) == 0) {
     return NULL;
   }
 
   /* Make the initial array of leaves */
-  dmnsn_array *leaves = dmnsn_new_array(sizeof(dmnsn_prnode *));
+  dmnsn_array *leaves = dmnsn_new_array(sizeof(dmnsn_bvh_node *));
   DMNSN_ARRAY_FOREACH (dmnsn_object **, object, objects) {
-    dmnsn_prnode *node = dmnsn_new_prnode();
-    node->bounding_box = (*object)->bounding_box;
-    node->object       = *object;
+    dmnsn_bvh_node *node = dmnsn_new_bvh_leaf_node(*object);
+    node->data = DMNSN_PRTREE_LEFT; /* Mustn't be _LEAF */
     dmnsn_array_push(leaves, &node);
   }
 
@@ -339,339 +292,7 @@ dmnsn_make_prtree(const dmnsn_array *objects)
     leaves = new_leaves;
   }
 
-  dmnsn_prnode *root = *(dmnsn_prnode **)dmnsn_array_first(leaves);
+  dmnsn_bvh_node *root = *(dmnsn_bvh_node **)dmnsn_array_first(leaves);
   dmnsn_delete_array(leaves);
   return root;
 }
-
-/** Add an object or its children, if any, to an array. */
-static void
-dmnsn_split_add_object(dmnsn_array *objects, const dmnsn_object *object)
-{
-  if (object->split_children) {
-    DMNSN_ARRAY_FOREACH (const dmnsn_object **, child, object->children) {
-      dmnsn_split_add_object(objects, *child);
-    }
-  } else {
-    dmnsn_array_push(objects, &object);
-  }
-}
-
-/** Split unions to create the input for the PR-tree. */
-static dmnsn_array *
-dmnsn_split_objects(const dmnsn_array *objects)
-{
-  dmnsn_array *split = dmnsn_new_array(sizeof(dmnsn_object *));
-  DMNSN_ARRAY_FOREACH (const dmnsn_object **, object, objects) {
-    dmnsn_split_add_object(split, *object);
-  }
-  return split;
-}
-
-/** Split unbounded objects into a new array. */
-static dmnsn_array *
-dmnsn_split_unbounded(dmnsn_array *objects)
-{
-  dmnsn_array *unbounded = dmnsn_new_array(sizeof(dmnsn_object *));
-
-  dmnsn_object **array = dmnsn_array_first(objects);
-  size_t i, skip;
-  for (i = 0, skip = 0; i < dmnsn_array_size(objects); ++i) {
-    if (dmnsn_bounding_box_is_infinite(array[i]->bounding_box)) {
-      dmnsn_array_push(unbounded, &array[i]);
-      ++skip;
-    } else {
-      array[i - skip] = array[i];
-    }
-  }
-  dmnsn_array_resize(objects, i - skip);
-
-  return unbounded;
-}
-
-/** Recursively flatten a PR-tree into an array of flat nodes. */
-static void
-dmnsn_flatten_prtree_recursive(dmnsn_prnode *node, dmnsn_array *flat)
-{
-  size_t currenti = dmnsn_array_size(flat);
-  dmnsn_array_resize(flat, currenti + 1);
-  dmnsn_flat_prnode *flatnode = dmnsn_array_at(flat, currenti);
-
-  flatnode->bounding_box = node->bounding_box;
-  flatnode->object       = node->object;
-
-  for (size_t i = 0; i < DMNSN_PRTREE_B && node->children[i]; ++i) {
-    dmnsn_flatten_prtree_recursive(node->children[i], flat);
-  }
-
-  /* Array could have been realloc()'d somewhere else above */
-  flatnode = dmnsn_array_at(flat, currenti);
-  flatnode->skip = dmnsn_array_size(flat) - currenti;
-}
-
-/** Flatten a PR-tree into an array of flat nodes. */
-static dmnsn_array *
-dmnsn_flatten_prtree(dmnsn_prnode *root)
-{
-  dmnsn_array *flat = dmnsn_new_array(sizeof(dmnsn_flat_prnode));
-  if (root) {
-    dmnsn_flatten_prtree_recursive(root, flat);
-  }
-  return flat;
-}
-
-static size_t dmnsn_prtree_seq = 0;
-static pthread_mutex_t dmnsn_prtree_seq_mutex = PTHREAD_MUTEX_INITIALIZER;
-
-/* Construct a PR-tree from a bulk of objects */
-dmnsn_prtree *
-dmnsn_new_prtree(const dmnsn_array *objects)
-{
-  dmnsn_prtree *prtree = dmnsn_malloc(sizeof(dmnsn_prtree));
-
-  dmnsn_array *bounded = dmnsn_split_objects(objects);
-  prtree->unbounded = dmnsn_split_unbounded(bounded);
-
-  dmnsn_prnode *root = dmnsn_make_prtree(bounded);
-  prtree->bounded = dmnsn_flatten_prtree(root);
-
-  dmnsn_delete_prnode(root);
-  dmnsn_delete_array(bounded);
-
-  if (dmnsn_array_size(prtree->unbounded) > 0) {
-    prtree->bounding_box = dmnsn_infinite_bounding_box();
-  } else if (dmnsn_array_size(prtree->bounded) > 0) {
-    dmnsn_flat_prnode *root = dmnsn_array_first(prtree->bounded);
-    prtree->bounding_box = root->bounding_box;
-  } else {
-    prtree->bounding_box = dmnsn_zero_bounding_box();
-  }
-
-  dmnsn_lock_mutex(&dmnsn_prtree_seq_mutex);
-    prtree->id = dmnsn_prtree_seq++;
-  dmnsn_unlock_mutex(&dmnsn_prtree_seq_mutex);
-
-  return prtree;
-}
-
-/** Free a PR-tree. */
-void
-dmnsn_delete_prtree(dmnsn_prtree *tree)
-{
-  if (tree) {
-    dmnsn_delete_array(tree->bounded);
-    dmnsn_delete_array(tree->unbounded);
-    dmnsn_free(tree);
-  }
-}
-
-/** A line with pre-calculated reciprocals to avoid divisions. */
-typedef struct dmnsn_optimized_line {
-  dmnsn_vector x0;    /**< The origin of the line. */
-  dmnsn_vector n_inv; /**< The inverse of each component of the line's slope .*/
-} dmnsn_optimized_line;
-
-/** Precompute inverses for faster ray-box intersection tests. */
-static inline dmnsn_optimized_line
-dmnsn_optimize_line(dmnsn_line line)
-{
-  dmnsn_optimized_line optline = {
-    .x0    = line.x0,
-    .n_inv = dmnsn_new_vector(1.0/line.n.x, 1.0/line.n.y, 1.0/line.n.z)
-  };
-  return optline;
-}
-
-/** Ray-AABB intersection test, by the slab method.  Highly optimized. */
-static inline bool
-dmnsn_ray_box_intersection(dmnsn_optimized_line optline,
-                           dmnsn_bounding_box box, double t)
-{
-  /*
-   * This is actually correct, even though it appears not to handle edge cases
-   * (line.n.{x,y,z} == 0).  It works because the infinities that result from
-   * dividing by zero will still behave correctly in the comparisons.  Lines
-   * which are parallel to an axis and outside the box will have tmin == inf
-   * or tmax == -inf, while lines inside the box will have tmin and tmax
-   * unchanged.
-   */
-
-  double tx1 = (box.min.x - optline.x0.x)*optline.n_inv.x;
-  double tx2 = (box.max.x - optline.x0.x)*optline.n_inv.x;
-
-  double tmin = dmnsn_min(tx1, tx2);
-  double tmax = dmnsn_max(tx1, tx2);
-
-  double ty1 = (box.min.y - optline.x0.y)*optline.n_inv.y;
-  double ty2 = (box.max.y - optline.x0.y)*optline.n_inv.y;
-
-  tmin = dmnsn_max(tmin, dmnsn_min(ty1, ty2));
-  tmax = dmnsn_min(tmax, dmnsn_max(ty1, ty2));
-
-  double tz1 = (box.min.z - optline.x0.z)*optline.n_inv.z;
-  double tz2 = (box.max.z - optline.x0.z)*optline.n_inv.z;
-
-  tmin = dmnsn_max(tmin, dmnsn_min(tz1, tz2));
-  tmax = dmnsn_min(tmax, dmnsn_max(tz1, tz2));
-
-  return tmax >= dmnsn_max(0.0, tmin) && tmin < t;
-}
-
-/** The number of intersections to cache. */
-#define DMNSN_PRTREE_CACHE_SIZE 32
-
-/** An array of cached intersections. */
-typedef struct dmnsn_intersection_cache {
-  size_t i;
-  dmnsn_object *objects[DMNSN_PRTREE_CACHE_SIZE];
-} dmnsn_intersection_cache;
-
-/** The thread-specific intersection cache. */
-static pthread_key_t dmnsn_prtree_caches;
-/** Initialize the thread-specific pointer exactly once. */
-static pthread_once_t dmnsn_prtree_caches_once = PTHREAD_ONCE_INIT;
-
-static void
-dmnsn_delete_prtree_caches(void *caches)
-{
-  dmnsn_delete_array(caches);
-}
-
-static void
-dmnsn_initialize_prtree_caches(void)
-{
-  dmnsn_key_create(&dmnsn_prtree_caches, dmnsn_delete_prtree_caches);
-}
-
-static dmnsn_array *
-dmnsn_get_prtree_caches(void)
-{
-  dmnsn_once(&dmnsn_prtree_caches_once, dmnsn_initialize_prtree_caches);
-  return pthread_getspecific(dmnsn_prtree_caches);
-}
-
-/** Needed because pthreads doesn't destroy data from the main thread unless
-    it exits with pthread_exit(). */
-DMNSN_DESTRUCTOR static void
-dmnsn_delete_main_prtree_caches(void)
-{
-  dmnsn_delete_prtree_caches(dmnsn_get_prtree_caches());
-  dmnsn_key_delete(dmnsn_prtree_caches);
-}
-
-static dmnsn_intersection_cache *
-dmnsn_get_intersection_cache(size_t id)
-{
-  dmnsn_array *caches = dmnsn_get_prtree_caches();
-  if (!caches) {
-    caches = dmnsn_new_array(sizeof(dmnsn_intersection_cache));
-    dmnsn_setspecific(dmnsn_prtree_caches, caches);
-  }
-
-  while (dmnsn_array_size(caches) <= id) {
-    dmnsn_array_resize(caches, dmnsn_array_size(caches) + 1);
-    dmnsn_intersection_cache *cache = dmnsn_array_last(caches);
-    cache->i = 0;
-    for (size_t i = 0; i < DMNSN_PRTREE_CACHE_SIZE; ++i) {
-      cache->objects[i] = NULL;
-    }
-  }
-
-  return dmnsn_array_at(caches, id);
-}
-
-/** Test for a closer object intersection than we've found so far. */
-static inline bool
-dmnsn_closer_intersection(dmnsn_object *object, dmnsn_line ray,
-                          dmnsn_intersection *intersection, double *t)
-{
-  dmnsn_intersection local_intersection;
-  if (dmnsn_object_intersection(object, ray, &local_intersection)) {
-    if (local_intersection.t < *t) {
-      *intersection = local_intersection;
-      *t = local_intersection.t;
-      return true;
-    }
-  }
-  return false;
-}
-
-DMNSN_HOT bool
-dmnsn_prtree_intersection(const dmnsn_prtree *tree, dmnsn_line ray,
-                          dmnsn_intersection *intersection, bool reset)
-{
-  double t = INFINITY;
-
-  /* Search the unbounded objects */
-  DMNSN_ARRAY_FOREACH (dmnsn_object **, object, tree->unbounded) {
-    dmnsn_closer_intersection(*object, ray, intersection, &t);
-  }
-
-  /* Precalculate 1.0/ray.n.{x,y,z} to save time in intersection tests */
-  dmnsn_optimized_line optline = dmnsn_optimize_line(ray);
-
-  /* Search the intersection cache */
-  dmnsn_intersection_cache *cache = dmnsn_get_intersection_cache(tree->id);
-  if (dmnsn_unlikely(reset)) {
-    cache->i = 0;
-  }
-  dmnsn_object *cached = NULL, *found = NULL;
-  if (dmnsn_likely(cache->i < DMNSN_PRTREE_CACHE_SIZE)) {
-    cached = cache->objects[cache->i];
-  }
-  if (cached && dmnsn_ray_box_intersection(optline, cached->bounding_box, t)) {
-    if (dmnsn_closer_intersection(cached, ray, intersection, &t)) {
-      found = cached;
-    }
-  }
-
-  /* Search the bounded objects */
-  dmnsn_flat_prnode *node = dmnsn_array_first(tree->bounded);
-  dmnsn_flat_prnode *last = dmnsn_array_last(tree->bounded);
-  while (node <= last) {
-    if (dmnsn_ray_box_intersection(optline, node->bounding_box, t)) {
-      if (node->object && node->object != cached) {
-        if (dmnsn_closer_intersection(node->object, ray, intersection, &t)) {
-          found = node->object;
-        }
-      }
-      ++node;
-    } else {
-      node += node->skip;
-    }
-  }
-
-  /* Update the cache */
-  if (dmnsn_likely(cache->i < DMNSN_PRTREE_CACHE_SIZE)) {
-    cache->objects[cache->i] = found;
-    ++cache->i;
-  }
-
-  return !isinf(t);
-}
-
-DMNSN_HOT bool
-dmnsn_prtree_inside(const dmnsn_prtree *tree, dmnsn_vector point)
-{
-  /* Search the unbounded objects */
-  DMNSN_ARRAY_FOREACH (dmnsn_object **, object, tree->unbounded) {
-    if (dmnsn_object_inside(*object, point))
-      return true;
-  }
-
-  /* Search the bounded objects */
-  dmnsn_flat_prnode *node = dmnsn_array_first(tree->bounded);
-  dmnsn_flat_prnode *last = dmnsn_array_last(tree->bounded);
-  while (node <= last) {
-    if (dmnsn_bounding_box_contains(node->bounding_box, point)) {
-      if (node->object && dmnsn_object_inside(node->object, point)) {
-        return true;
-      }
-      ++node;
-    } else {
-      node += node->skip;
-    }
-  }
-
-  return false;
-}
diff --git a/libdimension/prtree.h b/libdimension/prtree.h
index 733f854..0a8224a 100644
--- a/libdimension/prtree.h
+++ b/libdimension/prtree.h
@@ -1,5 +1,5 @@
 /*************************************************************************
- * Copyright (C) 2010-2011 Tavian Barnes <tavianator@tavianator.com>     *
+ * Copyright (C) 2010-2012 Tavian Barnes <tavianator@tavianator.com>     *
  *                                                                       *
  * This file is part of The Dimension Library.                           *
  *                                                                       *
@@ -26,27 +26,5 @@
  * from B-trees.
  */
 
-#include <stdbool.h>
-
-/** A priority R-tree; the spatial subdivision method used for intersection
-    queries against the scene graph. */
-typedef struct dmnsn_prtree {
-  dmnsn_bounding_box bounding_box; /**< The bounding box for the whole scene. */
-  dmnsn_array *unbounded;          /**< The unbounded objects. */
-  dmnsn_array *bounded;            /**< A PR-tree of the bounded objects. */
-  size_t id;                       /**< A unique ID for the PR-tree. */
-} dmnsn_prtree;
-
 /** Create a PR-tree. */
-DMNSN_INTERNAL dmnsn_prtree *dmnsn_new_prtree(const dmnsn_array *objects);
-/** Delete a PR-tree. */
-DMNSN_INTERNAL void dmnsn_delete_prtree(dmnsn_prtree *tree);
-
-/** Find the closest ray-object intersection in the tree. */
-DMNSN_INTERNAL bool dmnsn_prtree_intersection(const dmnsn_prtree *tree,
-                                              dmnsn_line ray,
-                                              dmnsn_intersection *intersection,
-                                              bool reset);
-/** Determine whether a point is inside any object in the tree. */
-DMNSN_INTERNAL bool dmnsn_prtree_inside(const dmnsn_prtree *tree,
-                                        dmnsn_vector point);
+DMNSN_INTERNAL dmnsn_bvh_node *dmnsn_new_prtree(const dmnsn_array *objects);
diff --git a/libdimension/ray_trace.c b/libdimension/ray_trace.c
index 127eb1e..7587e2c 100644
--- a/libdimension/ray_trace.c
+++ b/libdimension/ray_trace.c
@@ -1,5 +1,5 @@
 /*************************************************************************
- * Copyright (C) 2010-2011 Tavian Barnes <tavianator@tavianator.com>     *
+ * Copyright (C) 2010-2012 Tavian Barnes <tavianator@tavianator.com>     *
  *                                                                       *
  * This file is part of The Dimension Library.                           *
  *                                                                       *
@@ -34,7 +34,7 @@
 typedef struct {
   dmnsn_future *future;
   dmnsn_scene *scene;
-  dmnsn_prtree *prtree;
+  dmnsn_bvh *bvh;
 } dmnsn_ray_trace_payload;
 
 /* Ray-trace a scene */
@@ -81,7 +81,7 @@ dmnsn_ray_trace_scene_thread(void *ptr)
 
   /* Time the bounding tree construction */
   dmnsn_timer_start(&payload->scene->bounding_timer);
-    payload->prtree = dmnsn_new_prtree(payload->scene->objects);
+    payload->bvh = dmnsn_new_bvh(payload->scene->objects, DMNSN_BVH_PRTREE);
   dmnsn_timer_stop(&payload->scene->bounding_timer);
 
   /* Set up the future object */
@@ -93,7 +93,7 @@ dmnsn_ray_trace_scene_thread(void *ptr)
                                          payload, payload->scene->nthreads);
   dmnsn_timer_stop(&payload->scene->render_timer);
 
-  dmnsn_delete_prtree(payload->prtree);
+  dmnsn_delete_bvh(payload->bvh);
   dmnsn_free(payload);
 
   return ret;
@@ -111,7 +111,7 @@ typedef struct dmnsn_rtstate {
   const dmnsn_intersection *intersection;
   const dmnsn_texture *texture;
   const dmnsn_interior *interior;
-  const dmnsn_prtree *prtree;
+  const dmnsn_bvh *bvh;
   unsigned int reclevel;
 
   dmnsn_vector r;
@@ -146,12 +146,12 @@ dmnsn_ray_trace_scene_concurrent(void *ptr, unsigned int thread,
   const dmnsn_ray_trace_payload *payload = ptr;
   dmnsn_future *future = payload->future;
   dmnsn_scene *scene = payload->scene;
-  dmnsn_prtree *prtree = payload->prtree;
+  dmnsn_bvh *bvh = payload->bvh;
 
   dmnsn_rtstate state = {
     .parent = NULL,
     .scene  = scene,
-    .prtree = prtree,
+    .bvh = bvh,
   };
 
   /* Iterate through each pixel */
@@ -230,8 +230,8 @@ dmnsn_ray_shoot(dmnsn_rtstate *state, dmnsn_line ray)
 
   dmnsn_intersection intersection;
   bool reset = state->reclevel == state->scene->reclimit - 1;
-  dmnsn_prtree_intersection(state->prtree, ray, &intersection, reset);
-  if (dmnsn_prtree_intersection(state->prtree, ray, &intersection, reset)) {
+  dmnsn_bvh_intersection(state->bvh, ray, &intersection, reset);
+  if (dmnsn_bvh_intersection(state->bvh, ray, &intersection, reset)) {
     /* Found an intersection */
     dmnsn_rtstate_initialize(state, &intersection);
 
@@ -360,8 +360,8 @@ dmnsn_trace_light_ray(dmnsn_rtstate *state, const dmnsn_light *light)
 
   /* Test for shadow ray intersections */
   dmnsn_intersection shadow_caster;
-  bool in_shadow = dmnsn_prtree_intersection(state->prtree, shadow_ray,
-                                             &shadow_caster, false);
+  bool in_shadow = dmnsn_bvh_intersection(state->bvh, shadow_ray,
+                                          &shadow_caster, false);
   if (!in_shadow || !light->shadow_fn(light, shadow_caster.t)) {
     return true;
   }
diff --git a/libdimension/tests/prtree.c b/libdimension/tests/prtree.c
index d89a18d..7e69acd 100644
--- a/libdimension/tests/prtree.c
+++ b/libdimension/tests/prtree.c
@@ -1,5 +1,5 @@
 /*************************************************************************
- * Copyright (C) 2010-2011 Tavian Barnes <tavianator@tavianator.com>     *
+ * Copyright (C) 2010-2012 Tavian Barnes <tavianator@tavianator.com>     *
  *                                                                       *
  * This file is part of The Dimension Test Suite.                        *
  *                                                                       *
@@ -21,6 +21,7 @@
  * Basic tests of PR-trees
  */
 
+#include "../bvh.c"
 #include "../prtree.c"
 #include "../threads.c"
 #include "../future.c"
@@ -73,7 +74,7 @@ main(void)
     dmnsn_array_push(objects, &object);
   }
 
-  dmnsn_prtree *prtree = dmnsn_new_prtree(objects);
+  dmnsn_bvh *bvh = dmnsn_new_bvh(objects, DMNSN_BVH_PRTREE);
 
   dmnsn_intersection intersection;
   dmnsn_line ray = dmnsn_new_line(
@@ -81,7 +82,7 @@ main(void)
     dmnsn_new_vector(0.0, 0.0, 1.0)
   );
 
-  if (!dmnsn_prtree_intersection(prtree, ray, &intersection, true)) {
+  if (!dmnsn_bvh_intersection(bvh, ray, &intersection, true)) {
     fprintf(stderr, "--- Didn't find intersection! ---\n");
     return EXIT_FAILURE;
   }
@@ -93,7 +94,7 @@ main(void)
     return EXIT_FAILURE;
   }
 
-  dmnsn_delete_prtree(prtree);
+  dmnsn_delete_bvh(bvh);
   DMNSN_ARRAY_FOREACH (dmnsn_object **, object, objects) {
     dmnsn_delete_object(*object);
   }