Cache previous intersections in dmnsn_prtree_intersection().

Due to geometric locality of rays, this provides a very large speedup for
most scenes.

6 files changed, 75 insertions, 10 deletions

diff --git a/bench/libdimension/prtree.c b/bench/libdimension/prtree.c
index 6e793e4..b6f778a 100644
--- a/bench/libdimension/prtree.c
+++ b/bench/libdimension/prtree.c
@@ -96,7 +96,7 @@ main(void)
   dmnsn_intersection intersection;
 
   sandglass_bench_fine(&sandglass, {
-    dmnsn_prtree_intersection(tree, ray, &intersection);
+    dmnsn_prtree_intersection(tree, ray, &intersection, true);
   });
   printf("dmnsn_prtree_intersection(): %ld\n", sandglass.grains);
 
diff --git a/libdimension/csg.c b/libdimension/csg.c
index 9d02174..05f0b60 100644
--- a/libdimension/csg.c
+++ b/libdimension/csg.c
@@ -54,7 +54,7 @@ dmnsn_csg_union_intersection_fn(const dmnsn_object *csg,
                                 dmnsn_intersection *intersection)
 {
   dmnsn_prtree *prtree = csg->ptr;
-  return dmnsn_prtree_intersection(prtree, line, intersection);
+  return dmnsn_prtree_intersection(prtree, line, intersection, true);
 }
 
 /** CSG union inside callback. */
diff --git a/libdimension/prtree.c b/libdimension/prtree.c
index 4df8c63..9a329ce 100644
--- a/libdimension/prtree.c
+++ b/libdimension/prtree.c
@@ -25,6 +25,7 @@
  */
 
 #include "dimension-impl.h"
+#include <pthread.h>
 #include <stdlib.h>
 
 /** Number of children per PR-node. */
@@ -586,6 +587,9 @@ dmnsn_flatten_prtree(dmnsn_prtree_node *root)
   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)
@@ -611,6 +615,14 @@ dmnsn_new_prtree(const dmnsn_array *objects)
     prtree->bounding_box = dmnsn_zero_bounding_box();
   }
 
+  if (pthread_mutex_lock(&dmnsn_prtree_seq_mutex) != 0) {
+    dmnsn_error("Couldn't lock mutex.");
+  }
+  prtree->id = dmnsn_prtree_seq++;
+  if (pthread_mutex_unlock(&dmnsn_prtree_seq_mutex) != 0) {
+    dmnsn_error("Couldn't unlock mutex.");
+  }
+
   return prtree;
 }
 
@@ -677,9 +689,21 @@ dmnsn_ray_box_intersection(dmnsn_optimized_line optline,
   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 __thread dmnsn_array *dmnsn_prtree_cache = NULL;
+
 DMNSN_HOT bool
 dmnsn_prtree_intersection(const dmnsn_prtree *tree, dmnsn_line ray,
-                          dmnsn_intersection *intersection)
+                          dmnsn_intersection *intersection, bool reset)
 {
   double t = INFINITY;
 
@@ -697,18 +721,52 @@ dmnsn_prtree_intersection(const dmnsn_prtree *tree, dmnsn_line ray,
   /* 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 */
+  if (!dmnsn_prtree_cache) {
+    dmnsn_prtree_cache = dmnsn_new_array(sizeof(dmnsn_intersection_cache));
+  }
+  while (dmnsn_array_size(dmnsn_prtree_cache) <= tree->id) {
+    dmnsn_array_resize(dmnsn_prtree_cache,
+                       dmnsn_array_size(dmnsn_prtree_cache) + 1);
+    dmnsn_intersection_cache *cache = dmnsn_array_last(dmnsn_prtree_cache);
+    cache->i = 0;
+    for (size_t i = 0; i < DMNSN_PRTREE_CACHE_SIZE; ++i) {
+      cache->objects[i] = NULL;
+    }
+  }
+  dmnsn_intersection_cache *cache = dmnsn_array_first(dmnsn_prtree_cache);
+  cache += tree->id;
+  if (reset) {
+    cache->i = 0;
+  }
+  dmnsn_object *cached = NULL, *found = NULL;
+  if (cache->i < DMNSN_PRTREE_CACHE_SIZE) {
+    cached = cache->objects[cache->i];
+  }
+  if (cached && dmnsn_ray_box_intersection(optline, cached->bounding_box, t)) {
+    dmnsn_intersection local_intersection;
+    if (dmnsn_object_intersection(cached, ray, &local_intersection)) {
+      if (local_intersection.t < t) {
+        *intersection = local_intersection;
+        t = local_intersection.t;
+        found = cached;
+      }
+    }
+  }
+
   /* Search the bounded objects */
   dmnsn_flat_prnode *node = dmnsn_array_first(tree->bounded);
   while ((size_t)(node - (dmnsn_flat_prnode *)dmnsn_array_first(tree->bounded))
          < dmnsn_array_size(tree->bounded))
   {
     if (dmnsn_ray_box_intersection(optline, node->bounding_box, t)) {
-      if (node->object) {
+      if (node->object && node->object != cached) {
         dmnsn_intersection local_intersection;
         if (dmnsn_object_intersection(node->object, ray, &local_intersection)) {
           if (local_intersection.t < t) {
             *intersection = local_intersection;
             t = local_intersection.t;
+            found = node->object;
           }
         }
       }
@@ -719,6 +777,11 @@ dmnsn_prtree_intersection(const dmnsn_prtree *tree, dmnsn_line ray,
     }
   }
 
+  /* Update the cache */
+  if (cache->i < DMNSN_PRTREE_CACHE_SIZE) {
+    cache->objects[cache->i++] = found;
+  }
+
   return !isinf(t);
 }
 
diff --git a/libdimension/prtree.h b/libdimension/prtree.h
index 86707cf..26d6c1e 100644
--- a/libdimension/prtree.h
+++ b/libdimension/prtree.h
@@ -37,6 +37,7 @@ 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. */
@@ -46,7 +47,7 @@ void dmnsn_delete_prtree(dmnsn_prtree *tree);
 
 /** Find the closest ray-object intersection in the tree. */
 bool dmnsn_prtree_intersection(const dmnsn_prtree *tree, dmnsn_line ray,
-                               dmnsn_intersection *intersection);
+                               dmnsn_intersection *intersection, bool reset);
 /** Determine whether a point is inside any object in the tree. */
 bool dmnsn_prtree_inside(const dmnsn_prtree *tree, dmnsn_vector point);
 
diff --git a/libdimension/raytrace.c b/libdimension/raytrace.c
index 7a797bf..a35b1e6 100644
--- a/libdimension/raytrace.c
+++ b/libdimension/raytrace.c
@@ -253,10 +253,10 @@ dmnsn_raytrace_light_ray(const dmnsn_raytrace_state *state,
   dmnsn_color color = light->light_fn(light, state->r);
 
   unsigned int reclevel = state->reclevel;
-  while (reclevel) {
+  while (reclevel > 0) {
     dmnsn_intersection shadow_caster;
-    bool shadow_casted
-      = dmnsn_prtree_intersection(state->prtree, shadow_ray, &shadow_caster);
+    bool shadow_casted = dmnsn_prtree_intersection(state->prtree, shadow_ray,
+                                                   &shadow_caster, false);
 
     if (!shadow_casted || shadow_caster.t > 1.0) {
       break;
@@ -416,7 +416,8 @@ dmnsn_raytrace_shoot(dmnsn_raytrace_state *state, dmnsn_line ray)
   dmnsn_color color = dmnsn_raytrace_background(state, ray);
 
   dmnsn_intersection intersection;
-  if (dmnsn_prtree_intersection(state->prtree, ray, &intersection)) {
+  bool reset = state->reclevel == state->scene->reclimit - 1;
+  if (dmnsn_prtree_intersection(state->prtree, ray, &intersection, reset)) {
     state->intersection = &intersection;
     state->r = dmnsn_line_point(state->intersection->ray,
                                 state->intersection->t);
diff --git a/tests/libdimension/prtree.c b/tests/libdimension/prtree.c
index da69b97..cb36153 100644
--- a/tests/libdimension/prtree.c
+++ b/tests/libdimension/prtree.c
@@ -80,7 +80,7 @@ main(void)
     dmnsn_new_vector(0.0, 0.0, 1.0)
   );
 
-  if (!dmnsn_prtree_intersection(prtree, ray, &intersection)) {
+  if (!dmnsn_prtree_intersection(prtree, ray, &intersection, true)) {
     fprintf(stderr, "--- Didn't find intersection! ---\n");
     return EXIT_FAILURE;
   }