Decouple movie tracing through metric slabs

src/geodesic.c

@@ -1,5 +1,6 @@
 #include "geodesic.h"
 #include <math.h>
+#include <stddef.h>
 
 typedef GeodesicRayState State;
 typedef struct {
@@ -29,11 +30,11 @@ static int invert(double a[3][3], double b[3][3]) {
 }
 
 /* Equation (4) and (5) of Bohn et al., arXiv:1410.7775. */
-static int rhs(const SpacetimeSource *source, double t, const State *s,
+static int rhs(const MetricSlab *slab, double t, const State *s,
                Derivative *out) {
   MetricData m;
   double inv[3][3], up[3] = {0}, da_pi = 0, k_pi_pi = 0;
-  if (spacetime_eval(source, t, s->x, &m) || m.alpha <= 0 ||
+  if (spacetime_slab_eval(slab, t, s->x, &m) || m.alpha <= 0 ||
       invert(m.gamma, inv))
     return -1;
   for (int i = 0; i < 3; i++)
@@ -79,19 +80,19 @@ static State add(const State *s, const Derivative *d, double h) {
   return r;
 }
 
-static int rk4(const SpacetimeSource *source, double t, double h, State *s) {
+static int rk4(const MetricSlab *slab, double t, double h, State *s) {
   Derivative a, b, c, d;
   State q;
-  if (rhs(source, t, s, &a))
+  if (rhs(slab, t, s, &a))
     return -1;
   q = add(s, &a, h / 2);
-  if (rhs(source, t + h / 2, &q, &b))
+  if (rhs(slab, t + h / 2, &q, &b))
     return -1;
   q = add(s, &b, h / 2);
-  if (rhs(source, t + h / 2, &q, &c))
+  if (rhs(slab, t + h / 2, &q, &c))
     return -1;
   q = add(s, &c, h);
-  if (rhs(source, t + h, &q, &d))
+  if (rhs(slab, t + h, &q, &d))
     return -1;
   for (int i = 0; i < 3; i++) {
     s->x[i] += h * (a.x[i] + 2 * b.x[i] + 2 * c.x[i] + d.x[i]) / 6;
@@ -104,13 +105,13 @@ static int rk4(const SpacetimeSource *source, double t, double h, State *s) {
   return 0;
 }
 
-int geodesic_initialize_past_ray(const SpacetimeSource *source,
+int geodesic_initialize_past_ray(const MetricSlab *slab,
                                  const ObserverState *o, const double n[3],
                                  State *s) {
   MetricData m;
   double k[4] = {o->tetrad[0][0], o->tetrad[0][1], o->tetrad[0][2],
                  o->tetrad[0][3]};
-  if (spacetime_eval(source, o->coordinate_time, o->coordinate_position, &m) ||
+  if (spacetime_slab_eval(slab, o->coordinate_time, o->coordinate_position, &m) ||
       m.alpha <= 0)
     return -1;
   for (int a = 0; a < 3; a++)
@@ -131,11 +132,11 @@ int geodesic_initialize_past_ray(const SpacetimeSource *source,
   return isfinite(s->log_alpha_p0) ? 0 : -1;
 }
 
-static int escaped_direction(const SpacetimeSource *source, double t,
+static int escaped_direction(const MetricSlab *slab, double t,
                              const State *s, double n[3]) {
   MetricData m;
   double inv[3][3], norm = 0;
-  if (spacetime_eval(source, t, s->x, &m) || invert(m.gamma, inv))
+  if (spacetime_slab_eval(slab, t, s->x, &m) || invert(m.gamma, inv))
     return -1;
   for (int i = 0; i < 3; i++) {
     n[i] = 0;
@@ -151,9 +152,9 @@ static int escaped_direction(const SpacetimeSource *source, double t,
 }
 
 GeodesicAdvanceResult geodesic_advance_past_ray(
-    const SpacetimeSource *source, State *s, double slab_left_time,
+    const MetricSlab *slab, State *s, double slab_left_time,
     const GeodesicTraceConfig *config, RayEndpoint *out) {
-  if (!source || !s || !config || !out || config->coordinate_time_step <= 0 ||
+  if (!slab || !s || !config || !out || config->coordinate_time_step <= 0 ||
       !config->max_steps || !isfinite(slab_left_time) ||
       slab_left_time > s->coordinate_time)
     return GEODESIC_ADVANCE_FAILED;
@@ -164,12 +165,12 @@ GeodesicAdvanceResult geodesic_advance_past_ray(
       return GEODESIC_ADVANCE_TERMINATED;
     }
     SpacetimeRayStatus status =
-        spacetime_classify(source, s->coordinate_time, s->x);
+        spacetime_slab_classify(slab, s->coordinate_time, s->x);
     if (status != SPACETIME_RAY_ACTIVE) {
       out->status = status == SPACETIME_RAY_ESCAPED ? RAY_ENDPOINT_ESCAPED
                                                     : RAY_ENDPOINT_CAPTURED;
       if (out->status == RAY_ENDPOINT_ESCAPED &&
-          escaped_direction(source, s->coordinate_time, s, out->n_infinity) == 0)
+          escaped_direction(slab, s->coordinate_time, s, out->n_infinity) == 0)
         out->frequency_ratio = exp(-s->log_alpha_p0);
       else if (out->status == RAY_ENDPOINT_ESCAPED)
         out->status = RAY_ENDPOINT_INTEGRATION_FAILURE;
@@ -183,7 +184,7 @@ GeodesicAdvanceResult geodesic_advance_past_ray(
     }
     const double h = -fmin(config->coordinate_time_step,
                            s->coordinate_time - slab_left_time);
-    if (rk4(source, s->coordinate_time, h, s))
+    if (rk4(slab, s->coordinate_time, h, s))
       return GEODESIC_ADVANCE_FAILED;
     s->coordinate_time += h;
     ++s->steps;
@@ -199,14 +200,23 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source,
                      .magnification = 1,
                      .status = RAY_ENDPOINT_INTEGRATION_FAILURE};
   State state;
+  MetricSlab *slab = NULL;
   if (!source || !observer || !config || config->coordinate_time_step <= 0 ||
-      !config->max_steps || fabs(dot(n, n) - 1) > 1e-10 ||
-      geodesic_initialize_past_ray(source, observer, n, &state))
+      !config->max_steps || fabs(dot(n, n) - 1) > 1e-10)
     return out;
+  if (spacetime_load_slab(source, observer->coordinate_time,
+                           observer->coordinate_time -
+                               config->coordinate_time_step * config->max_steps - 1.0,
+                           &slab) ||
+      geodesic_initialize_past_ray(slab, observer, n, &state)) {
+    spacetime_free_slab(slab);
+    return out;
+  }
   const double last_time =
       observer->coordinate_time - config->coordinate_time_step * config->max_steps;
-  if (geodesic_advance_past_ray(source, &state, last_time, config, &out) ==
+  if (geodesic_advance_past_ray(slab, &state, last_time, config, &out) ==
       GEODESIC_ADVANCE_ACTIVE)
     out.status = RAY_ENDPOINT_MAX_STEPS;
+  spacetime_free_slab(slab);
   return out;
 }

src/geodesic.h

@@ -48,12 +48,12 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source,
                                 const ObserverState *observer,
                                 const double camera_direction[3],
                                 const GeodesicTraceConfig *config);
-int geodesic_initialize_past_ray(const SpacetimeSource *source,
+int geodesic_initialize_past_ray(const MetricSlab *slab,
                                  const ObserverState *observer,
                                  const double camera_direction[3],
                                  GeodesicRayState *state);
 GeodesicAdvanceResult geodesic_advance_past_ray(
-    const SpacetimeSource *source, GeodesicRayState *state,
+    const MetricSlab *slab, GeodesicRayState *state,
     double slab_left_time, const GeodesicTraceConfig *config,
     RayEndpoint *endpoint);
 #endif

src/main.c

@@ -255,15 +255,19 @@ static int render_movie(const Settings *s, const StarCatalog *catalog,
     goto done;
   for (size_t f = 0; f < movie.frame_count; ++f)
     for (size_t v = 0; v < movie.frames[f].mesh.vertex_count; ++v)
-      if (ray_pool_append(&rays, spacetime, &movie.frames[f].observer,
+      if (ray_pool_append(&rays, &movie.frames[f].observer,
                           movie.frames[f].mesh.vertices[v].camera_direction,
                           f, v))
         goto done;
   double slab_hi = movie.frames[movie.frame_count - 1].coordinate_time;
   while (ray_pool_has_live(&rays)) {
     const double slab_lo = slab_hi - s->slab_duration;
-    ray_pool_activate_in_time_range(&rays, slab_hi, slab_lo);
-    ray_pool_advance_active(&rays, spacetime, slab_lo, &trace);
+    MetricSlab *slab = NULL;
+    if (spacetime_load_slab(spacetime, slab_hi, slab_lo, &slab))
+      goto done;
+    ray_pool_activate_in_time_range(&rays, slab);
+    ray_pool_advance_active(&rays, slab, &trace);
+    spacetime_free_slab(slab);
     slab_hi = slab_lo;
   }
   for (size_t i = 0; i < rays.count; ++i) {

src/ray.c

@@ -9,7 +9,8 @@ int ray_pool_init(RayPool *p, size_t capacity) {
   *p = (RayPool){.capacity = capacity};
 #define RAY_ALLOC(field) (p->field = calloc(capacity, sizeof *p->field))
   if (!(RAY_ALLOC(t) && RAY_ALLOC(x0) && RAY_ALLOC(x1) && RAY_ALLOC(x2) &&
-        RAY_ALLOC(p0) && RAY_ALLOC(p1) && RAY_ALLOC(p2) &&
+        RAY_ALLOC(p0) && RAY_ALLOC(p1) && RAY_ALLOC(p2) && RAY_ALLOC(observer) &&
+        RAY_ALLOC(direction0) && RAY_ALLOC(direction1) && RAY_ALLOC(direction2) &&
         RAY_ALLOC(log_alpha_p0) && RAY_ALLOC(steps) && RAY_ALLOC(frame_id) &&
         RAY_ALLOC(vertex_id) && RAY_ALLOC(status) && RAY_ALLOC(endpoint))) {
     ray_pool_destroy(p);
@@ -19,20 +20,19 @@ int ray_pool_init(RayPool *p, size_t capacity) {
   return 0;
 }
 
-int ray_pool_append(RayPool *p, const SpacetimeSource *source,
-                    const ObserverState *observer, const double direction[3],
+int ray_pool_append(RayPool *p, const ObserverState *observer,
+                    const double direction[3],
                     size_t frame_id, size_t vertex_id) {
   if (p == NULL || p->count == p->capacity)
     return -1;
-  GeodesicRayState s;
   const size_t i = p->count;
-  if (geodesic_initialize_past_ray(source, observer, direction, &s))
+  if (observer == NULL || direction == NULL)
     return -1;
-  p->t[i] = s.coordinate_time;
-  p->x0[i] = s.x[0]; p->x1[i] = s.x[1]; p->x2[i] = s.x[2];
-  p->p0[i] = s.Pi[0]; p->p1[i] = s.Pi[1]; p->p2[i] = s.Pi[2];
-  p->log_alpha_p0[i] = s.log_alpha_p0;
-  p->steps[i] = s.steps;
+  p->t[i] = observer->coordinate_time;
+  p->observer[i] = observer;
+  p->direction0[i] = direction[0];
+  p->direction1[i] = direction[1];
+  p->direction2[i] = direction[2];
   p->frame_id[i] = frame_id;
   p->vertex_id[i] = vertex_id;
   p->status[i] = RAY_POOL_PENDING;
@@ -42,14 +42,29 @@ int ray_pool_append(RayPool *p, const SpacetimeSource *source,
   return 0;
 }
 
-void ray_pool_activate_in_time_range(RayPool *p, double t_hi, double t_lo) {
-  for (size_t i = 0; i < p->count; ++i)
-    if (p->status[i] == RAY_POOL_PENDING && p->t[i] <= t_hi && p->t[i] > t_lo)
-      p->status[i] = RAY_POOL_ACTIVE;
+void ray_pool_activate_in_time_range(RayPool *p, const MetricSlab *slab) {
+  for (size_t i = 0; i < p->count; ++i) {
+    if (p->status[i] != RAY_POOL_PENDING || p->t[i] > slab->t_hi ||
+        p->t[i] <= slab->t_lo)
+      continue;
+    GeodesicRayState state;
+    if (geodesic_initialize_past_ray(
+            slab, p->observer[i],
+            (double[]){p->direction0[i], p->direction1[i], p->direction2[i]},
+            &state)) {
+      p->status[i] = RAY_POOL_FAILED;
+      continue;
+    }
+    p->x0[i] = state.x[0]; p->x1[i] = state.x[1]; p->x2[i] = state.x[2];
+    p->p0[i] = state.Pi[0]; p->p1[i] = state.Pi[1]; p->p2[i] = state.Pi[2];
+    p->log_alpha_p0[i] = state.log_alpha_p0;
+    p->steps[i] = state.steps;
+    p->status[i] = RAY_POOL_ACTIVE;
+  }
 }
 
-void ray_pool_advance_active(RayPool *p, const SpacetimeSource *source,
-                             double t_lo, const GeodesicTraceConfig *config) {
+void ray_pool_advance_active(RayPool *p, const MetricSlab *slab,
+                             const GeodesicTraceConfig *config) {
 #pragma omp parallel for schedule(static)
   for (size_t i = 0; i < p->count; ++i) {
     if (p->status[i] != RAY_POOL_ACTIVE)
@@ -60,7 +75,7 @@ void ray_pool_advance_active(RayPool *p, const SpacetimeSource *source,
                            .log_alpha_p0 = p->log_alpha_p0[i],
                            .steps = p->steps[i]};
     const GeodesicAdvanceResult result =
-        geodesic_advance_past_ray(source, &s, t_lo, config, &p->endpoint[i]);
+        geodesic_advance_past_ray(slab, &s, slab->t_lo, config, &p->endpoint[i]);
     p->t[i] = s.coordinate_time;
     p->x0[i] = s.x[0]; p->x1[i] = s.x[1]; p->x2[i] = s.x[2];
     p->p0[i] = s.Pi[0]; p->p1[i] = s.Pi[1]; p->p2[i] = s.Pi[2];
@@ -85,7 +100,8 @@ int ray_pool_has_live(const RayPool *p) {
 void ray_pool_destroy(RayPool *p) {
   if (p == NULL)
     return;
-  free(p->t); free(p->x0); free(p->x1); free(p->x2);
+  free(p->t); free(p->x0); free(p->x1); free(p->x2); free(p->observer);
+  free(p->direction0); free(p->direction1); free(p->direction2);
   free(p->p0); free(p->p1); free(p->p2); free(p->log_alpha_p0);
   free(p->steps); free(p->frame_id); free(p->vertex_id); free(p->status);
   free(p->endpoint);

src/ray.h

@@ -15,6 +15,8 @@ typedef enum {
 
 typedef struct {
   double *t, *x0, *x1, *x2, *p0, *p1, *p2, *log_alpha_p0;
+  const ObserverState **observer;
+  double *direction0, *direction1, *direction2;
   unsigned int *steps;
   size_t *frame_id, *vertex_id;
   uint8_t *status;
@@ -23,12 +25,12 @@ typedef struct {
 } RayPool;
 
 int ray_pool_init(RayPool *pool, size_t capacity);
-int ray_pool_append(RayPool *pool, const SpacetimeSource *source,
-                    const ObserverState *observer, const double direction[3],
+int ray_pool_append(RayPool *pool, const ObserverState *observer,
+                    const double direction[3],
                     size_t frame_id, size_t vertex_id);
-void ray_pool_activate_in_time_range(RayPool *pool, double t_hi, double t_lo);
-void ray_pool_advance_active(RayPool *pool, const SpacetimeSource *source,
-                             double t_lo, const GeodesicTraceConfig *config);
+void ray_pool_activate_in_time_range(RayPool *pool, const MetricSlab *slab);
+void ray_pool_advance_active(RayPool *pool, const MetricSlab *slab,
+                             const GeodesicTraceConfig *config);
 int ray_pool_has_live(const RayPool *pool);
 void ray_pool_destroy(RayPool *pool);
 

src/spacetime.h

@@ -18,12 +18,26 @@ typedef enum {
 } SpacetimeRayStatus;
 
 typedef struct SpacetimeSource SpacetimeSource;
+typedef struct MetricSlab MetricSlab;
+
+struct MetricSlab {
+  const SpacetimeSource *source;
+  double t_hi, t_lo;
+  void *context;
+};
 
 typedef struct {
   int (*eval)(const SpacetimeSource *source, double t, const double x[3],
               MetricData *metric);
   SpacetimeRayStatus (*classify)(const SpacetimeSource *source, double t,
                                  const double x[3]);
+  int (*load_slab)(const SpacetimeSource *source, double t_hi, double t_lo,
+                   MetricSlab **out);
+  void (*free_slab)(MetricSlab *slab);
+  int (*eval_slab)(const MetricSlab *slab, double t, const double x[3],
+                   MetricData *metric);
+  SpacetimeRayStatus (*classify_slab)(const MetricSlab *slab, double t,
+                                      const double x[3]);
   void (*destroy)(SpacetimeSource *source);
 } SpacetimeOps;
 
@@ -44,5 +58,12 @@ int spacetime_eval(const SpacetimeSource *source, double t, const double x[3],
                    MetricData *metric);
 SpacetimeRayStatus spacetime_classify(const SpacetimeSource *source, double t,
                                       const double x[3]);
+int spacetime_load_slab(const SpacetimeSource *source, double t_hi, double t_lo,
+                        MetricSlab **out);
+void spacetime_free_slab(MetricSlab *slab);
+int spacetime_slab_eval(const MetricSlab *slab, double t, const double x[3],
+                        MetricData *metric);
+SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t,
+                                           const double x[3]);
 
 #endif

src/spacetime_common.c

@@ -1,6 +1,7 @@
 #include "spacetime.h"
 
 #include <stddef.h>
+#include <stdlib.h>
 
 void spacetime_destroy(SpacetimeSource *source) {
   if (source != NULL && source->ops != NULL)
@@ -20,3 +21,45 @@ SpacetimeRayStatus spacetime_classify(const SpacetimeSource *source, double t,
              ? SPACETIME_RAY_CAPTURED
              : source->ops->classify(source, t, x);
 }
+
+int spacetime_load_slab(const SpacetimeSource *source, double t_hi, double t_lo,
+                        MetricSlab **out) {
+  if (source == NULL || source->ops == NULL || out == NULL || t_lo >= t_hi)
+    return -1;
+  if (source->ops->load_slab != NULL)
+    return source->ops->load_slab(source, t_hi, t_lo, out);
+  MetricSlab *slab = malloc(sizeof *slab);
+  if (slab == NULL)
+    return -1;
+  *slab = (MetricSlab){.source = source, .t_hi = t_hi, .t_lo = t_lo};
+  *out = slab;
+  return 0;
+}
+
+void spacetime_free_slab(MetricSlab *slab) {
+  if (slab == NULL)
+    return;
+  if (slab->source != NULL && slab->source->ops != NULL &&
+      slab->source->ops->free_slab != NULL)
+    slab->source->ops->free_slab(slab);
+  else
+    free(slab);
+}
+
+int spacetime_slab_eval(const MetricSlab *slab, double t, const double x[3],
+                        MetricData *metric) {
+  if (slab == NULL || t < slab->t_lo || t > slab->t_hi)
+    return -1;
+  if (slab->source->ops->eval_slab != NULL)
+    return slab->source->ops->eval_slab(slab, t, x, metric);
+  return spacetime_eval(slab->source, t, x, metric);
+}
+
+SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t,
+                                           const double x[3]) {
+  if (slab == NULL || t < slab->t_lo || t > slab->t_hi)
+    return SPACETIME_RAY_CAPTURED;
+  if (slab->source->ops->classify_slab != NULL)
+    return slab->source->ops->classify_slab(slab, t, x);
+  return spacetime_classify(slab->source, t, x);
+}