diff options
Diffstat (limited to 'libdimension/prtree.c')
| -rw-r--r-- | libdimension/prtree.c | 439 |
1 files changed, 30 insertions, 409 deletions
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; -} |