Aggregate movie rays through time slabs

README.md

@@ -45,9 +45,12 @@ mkdir -p output/imgs
 This writes `minkowski_accel_000000.png` through
 `minkowski_accel_000060.png`.  The renderer treats the CSV as its observer
 input; the acceleration generator is only a reproducible flat-spacetime test
-fixture.  The current movie path still renders frames independently while the
+fixture.  Movie mode collects all current frame-mesh vertices into a single
-time-slab/RayPool scheduler is implemented next, so it must not yet be used
+SoA ray pool, activates rays as a newest-to-oldest coordinate-time scan reaches
-as a performance measurement for numerical-relativity data.  The current
+their observer event, and advances active rays to each slab boundary.  The
+analytic backends use logical slabs with no metric I/O; nmesh slab loading is
+the next backend step.  `--slab-duration` sets the coordinate-time width
+(default `64`) for this current fixed-mesh pass.  The current
 synthetic test catalog uses global default exposure `1e-3`; the accelerated
 benchmark explicitly uses `1e-5` because its physical Doppler blue shift
 otherwise clips the later frames.

src/geodesic.c

@@ -1,9 +1,7 @@
 #include "geodesic.h"
 #include <math.h>
 
-typedef struct {
+typedef GeodesicRayState State;
-  double x[3], Pi[3], log_alpha_p0;
-} State;
 typedef struct {
   double x[3], Pi[3], log_alpha_p0;
 } Derivative;
@@ -106,8 +104,9 @@ static int rk4(const SpacetimeSource *source, double t, double h, State *s) {
   return 0;
 }
 
-static int initialize(const SpacetimeSource *source, const ObserverState *o,
+int geodesic_initialize_past_ray(const SpacetimeSource *source,
-                      const double n[3], State *s) {
+                                 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]};
@@ -127,6 +126,8 @@ static int initialize(const SpacetimeSource *source, const ObserverState *o,
     s->Pi[i] /= m.alpha * k[0];
   }
   s->log_alpha_p0 = log(m.alpha * k[0]);
+  s->coordinate_time = o->coordinate_time;
+  s->steps = 0;
   return isfinite(s->log_alpha_p0) ? 0 : -1;
 }
 
@@ -149,6 +150,47 @@ static int escaped_direction(const SpacetimeSource *source, double t,
   return 0;
 }
 
+GeodesicAdvanceResult geodesic_advance_past_ray(
+    const SpacetimeSource *source, State *s, double slab_left_time,
+    const GeodesicTraceConfig *config, RayEndpoint *out) {
+  if (!source || !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;
+  while (s->coordinate_time > slab_left_time) {
+    if (config->capture_log_alpha_p0 > 0.0 &&
+        s->log_alpha_p0 >= config->capture_log_alpha_p0) {
+      out->status = RAY_ENDPOINT_CAPTURED;
+      return GEODESIC_ADVANCE_TERMINATED;
+    }
+    SpacetimeRayStatus status =
+        spacetime_classify(source, 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)
+        out->frequency_ratio = exp(-s->log_alpha_p0);
+      else if (out->status == RAY_ENDPOINT_ESCAPED)
+        out->status = RAY_ENDPOINT_INTEGRATION_FAILURE;
+      return out->status == RAY_ENDPOINT_INTEGRATION_FAILURE
+                 ? GEODESIC_ADVANCE_FAILED
+                 : GEODESIC_ADVANCE_TERMINATED;
+    }
+    if (s->steps >= config->max_steps) {
+      out->status = RAY_ENDPOINT_MAX_STEPS;
+      return GEODESIC_ADVANCE_TERMINATED;
+    }
+    const double h = -fmin(config->coordinate_time_step,
+                           s->coordinate_time - slab_left_time);
+    if (rk4(source, s->coordinate_time, h, s))
+      return GEODESIC_ADVANCE_FAILED;
+    s->coordinate_time += h;
+    ++s->steps;
+  }
+  return GEODESIC_ADVANCE_ACTIVE;
+}
+
 RayEndpoint geodesic_trace_past(const SpacetimeSource *source,
                                 const ObserverState *observer,
                                 const double n[3],
@@ -156,34 +198,15 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source,
   RayEndpoint out = {.frequency_ratio = 0,
                      .magnification = 1,
                      .status = RAY_ENDPOINT_INTEGRATION_FAILURE};
-  State s;
+  State state;
-  double t;
   if (!source || !observer || !config || config->coordinate_time_step <= 0 ||
       !config->max_steps || fabs(dot(n, n) - 1) > 1e-10 ||
-      initialize(source, observer, n, &s))
+      geodesic_initialize_past_ray(source, observer, n, &state))
     return out;
-  t = observer->coordinate_time;
+  const double last_time =
-  for (unsigned int i = 0; i < config->max_steps; i++) {
+      observer->coordinate_time - config->coordinate_time_step * config->max_steps;
-    if (config->capture_log_alpha_p0 > 0.0 &&
+  if (geodesic_advance_past_ray(source, &state, last_time, config, &out) ==
-        s.log_alpha_p0 >= config->capture_log_alpha_p0) {
+      GEODESIC_ADVANCE_ACTIVE)
-      out.status = RAY_ENDPOINT_CAPTURED;
+    out.status = RAY_ENDPOINT_MAX_STEPS;
-      return out;
-    }
-    SpacetimeRayStatus status = spacetime_classify(source, t, 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, t, &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;
-      return out;
-    }
-    if (rk4(source, t, -config->coordinate_time_step, &s))
-      return out;
-    t -= config->coordinate_time_step;
-  }
-  out.status = RAY_ENDPOINT_MAX_STEPS;
   return out;
 }

src/geodesic.h

@@ -28,10 +28,32 @@ typedef struct {
   double capture_log_alpha_p0;
 } GeodesicTraceConfig;
 
+typedef struct {
+  double coordinate_time;
+  double x[3];
+  double Pi[3];
+  double log_alpha_p0;
+  unsigned int steps;
+} GeodesicRayState;
+
+typedef enum {
+  GEODESIC_ADVANCE_ACTIVE = 0,
+  GEODESIC_ADVANCE_TERMINATED = 1,
+  GEODESIC_ADVANCE_FAILED = -1
+} GeodesicAdvanceResult;
+
 /* camera_direction is a unit vector in the observer's (forward, up, right)
  * tetrad. */
 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,
+                                 const ObserverState *observer,
+                                 const double camera_direction[3],
+                                 GeodesicRayState *state);
+GeodesicAdvanceResult geodesic_advance_past_ray(
+    const SpacetimeSource *source, GeodesicRayState *state,
+    double slab_left_time, const GeodesicTraceConfig *config,
+    RayEndpoint *endpoint);
 #endif

src/main.c

@@ -3,6 +3,7 @@
 #include "movie.h"
 #include "observer_track.h"
 #include "optics.h"
+#include "ray.h"
 #include "spacetime.h"
 
 #include <errno.h>
@@ -27,6 +28,7 @@ typedef struct {
   const char *frames_prefix;
   const char *write_minkowski_accel_track_path;
   double movie_start_time, movie_duration, movie_fps;
+  double slab_duration;
   double minkowski_proper_acceleration;
 } Settings;
 
@@ -104,6 +106,7 @@ static int parse_args(int argc, char **argv, Settings *s,
                   .frames_prefix = "frame",
                   .movie_duration = 2.0,
                   .movie_fps = 30.0,
+                  .slab_duration = 64.0,
                   .minkowski_proper_acceleration = 1.52};
   *write_path = NULL;
   for (int i = 1; i < argc; ++i) {
@@ -147,6 +150,8 @@ static int parse_args(int argc, char **argv, Settings *s,
                !parse_nonnegative(argv[++i], &s->movie_duration)) {
     } else if (!strcmp(argv[i], "--fps") && i + 1 < argc &&
                !parse_positive(argv[++i], &s->movie_fps)) {
+    } else if (!strcmp(argv[i], "--slab-duration") && i + 1 < argc &&
+               !parse_positive(argv[++i], &s->slab_duration)) {
     } else if (!strcmp(argv[i], "--proper-acceleration") && i + 1 < argc &&
                !parse_nonnegative(argv[++i],
                                   &s->minkowski_proper_acceleration)) {
@@ -233,21 +238,64 @@ static int render_movie(const Settings *s, const StarCatalog *catalog,
                         const SpacetimeSource *spacetime) {
   ObserverTrack track = {0};
   Movie movie = {0};
+  RayPool rays = {0};
+  const GeodesicTraceConfig trace = trace_config();
   int result = -1;
   if (s->observer_track_path == NULL ||
       observer_track_load_csv(&track, s->observer_track_path) ||
       movie_init(&movie, &track, s->movie_start_time, s->movie_duration,
-                 s->movie_fps))
+                 s->movie_fps) ||
+      movie_build_coarse_meshes(&movie, s->width, s->height,
+                                s->coarse_cell_pixels, s->horizontal_fov_deg))
     goto done;
+  size_t ray_count = 0;
+  for (size_t i = 0; i < movie.frame_count; ++i)
+    ray_count += movie.frames[i].mesh.vertex_count;
+  if (ray_count == 0 || ray_pool_init(&rays, ray_count))
+    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,
+                          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);
+    slab_hi = slab_lo;
+  }
+  for (size_t i = 0; i < rays.count; ++i) {
+    LensVertex *vertex = &movie.frames[rays.frame_id[i]].mesh.vertices[rays.vertex_id[i]];
+    vertex->status = rays.endpoint[i].status;
+    if (vertex->status == RAY_ENDPOINT_ESCAPED) {
+      for (int axis = 0; axis < 3; ++axis)
+        vertex->n_infinity[axis] = rays.endpoint[i].n_infinity[axis];
+      vertex->log_frequency_ratio = log(rays.endpoint[i].frequency_ratio);
+    }
+  }
   for (size_t i = 0; i < movie.frame_count; ++i) {
     char output_path[PATH_MAX];
-    if (frame_output_path(output_path, s, movie.frames[i].frame_id) ||
+    double *hdr = calloc((size_t)s->width * s->height * 3, sizeof *hdr);
-        render_observer_frame(s, catalog, spacetime, &movie.frames[i].observer,
+    if (hdr == NULL || frame_output_path(output_path, s, movie.frames[i].frame_id)) {
-                              output_path))
+      free(hdr);
+      goto done;
+    }
+    const size_t images = frame_splat_catalog(&movie.frames[i].mesh, catalog, hdr,
+                                              s->width, s->height, s->exposure, &s->psf);
+    if (s->draw_mesh)
+      frame_draw_mesh(&movie.frames[i].mesh, hdr, s->width, s->height, 0.5, 0.5);
+    const int write_result = write_tonemapped_image(output_path, hdr, s->width, s->height);
+    free(hdr);
+    fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n",
+            images, catalog->count, output_path, write_result == 0 ? "ok" : "write failed");
+    if (write_result)
       goto done;
   }
   result = 0;
 done:
+  ray_pool_destroy(&rays);
   movie_destroy(&movie);
   observer_track_destroy(&track);
   return result;

src/movie.c

@@ -33,9 +33,25 @@ int movie_init(Movie *movie, const ObserverTrack *track, double start_time,
   return 0;
 }
 
+int movie_build_coarse_meshes(Movie *movie, int width, int height,
+                              int cell_pixels, double horizontal_fov_deg) {
+  if (movie == NULL)
+    return -1;
+  for (size_t i = 0; i < movie->frame_count; ++i)
+    if (frame_lens_mesh_build_coarse(&movie->frames[i].mesh, width, height,
+                                     cell_pixels, horizontal_fov_deg)) {
+      for (size_t j = 0; j <= i; ++j)
+        frame_lens_mesh_destroy(&movie->frames[j].mesh);
+      return -1;
+    }
+  return 0;
+}
+
 void movie_destroy(Movie *movie) {
   if (movie == NULL)
     return;
+  for (size_t i = 0; i < movie->frame_count; ++i)
+    frame_lens_mesh_destroy(&movie->frames[i].mesh);
   free(movie->frames);
   *movie = (Movie){0};
 }

src/movie.h

@@ -2,6 +2,7 @@
 #define MOVIE_H
 
 #include "observer_track.h"
+#include "frame.h"
 
 #include <stddef.h>
 
@@ -10,6 +11,7 @@ typedef struct {
   double coordinate_time;
   double proper_time;
   ObserverState observer;
+  FrameLensMesh mesh;
 } MovieFrame;
 
 typedef struct {
@@ -19,6 +21,8 @@ typedef struct {
 
 int movie_init(Movie *movie, const ObserverTrack *track, double start_time,
                double duration, double frames_per_second);
+int movie_build_coarse_meshes(Movie *movie, int width, int height,
+                              int cell_pixels, double horizontal_fov_deg);
 void movie_destroy(Movie *movie);
 
 #endif

src/ray.c

@@ -0,0 +1,93 @@
+#include "ray.h"
+
+#include <omp.h>
+#include <stdlib.h>
+
+int ray_pool_init(RayPool *p, size_t capacity) {
+  if (p == NULL || capacity == 0)
+    return -1;
+  *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(log_alpha_p0) && RAY_ALLOC(steps) && RAY_ALLOC(frame_id) &&
+        RAY_ALLOC(vertex_id) && RAY_ALLOC(status) && RAY_ALLOC(endpoint))) {
+    ray_pool_destroy(p);
+    return -1;
+  }
+#undef RAY_ALLOC
+  return 0;
+}
+
+int ray_pool_append(RayPool *p, const SpacetimeSource *source,
+                    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))
+    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->frame_id[i] = frame_id;
+  p->vertex_id[i] = vertex_id;
+  p->status[i] = RAY_POOL_PENDING;
+  p->endpoint[i] = (RayEndpoint){.magnification = 1.0,
+                                  .status = RAY_ENDPOINT_INTEGRATION_FAILURE};
+  ++p->count;
+  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_advance_active(RayPool *p, const SpacetimeSource *source,
+                             double t_lo, 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)
+      continue;
+    GeodesicRayState s = {.coordinate_time = p->t[i],
+                           .x = {p->x0[i], p->x1[i], p->x2[i]},
+                           .Pi = {p->p0[i], p->p1[i], p->p2[i]},
+                           .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]);
+    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;
+    if (result == GEODESIC_ADVANCE_TERMINATED)
+      p->status[i] = RAY_POOL_TERMINATED;
+    else if (result == GEODESIC_ADVANCE_FAILED)
+      p->status[i] = RAY_POOL_FAILED;
+  }
+}
+
+int ray_pool_has_live(const RayPool *p) {
+  if (p == NULL)
+    return 0;
+  for (size_t i = 0; i < p->count; ++i)
+    if (p->status[i] == RAY_POOL_PENDING || p->status[i] == RAY_POOL_ACTIVE)
+      return 1;
+  return 0;
+}
+
+void ray_pool_destroy(RayPool *p) {
+  if (p == NULL)
+    return;
+  free(p->t); free(p->x0); free(p->x1); free(p->x2);
+  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);
+  *p = (RayPool){0};
+}

src/ray.h

@@ -0,0 +1,35 @@
+#ifndef RAY_H
+#define RAY_H
+
+#include "geodesic.h"
+
+#include <stddef.h>
+#include <stdint.h>
+
+typedef enum {
+  RAY_POOL_PENDING,
+  RAY_POOL_ACTIVE,
+  RAY_POOL_TERMINATED,
+  RAY_POOL_FAILED
+} RayPoolStatus;
+
+typedef struct {
+  double *t, *x0, *x1, *x2, *p0, *p1, *p2, *log_alpha_p0;
+  unsigned int *steps;
+  size_t *frame_id, *vertex_id;
+  uint8_t *status;
+  RayEndpoint *endpoint;
+  size_t count, capacity;
+} 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],
+                    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);
+int ray_pool_has_live(const RayPool *pool);
+void ray_pool_destroy(RayPool *pool);
+
+#endif