diff --git a/README.md b/README.md index 159c339..bae7cd7 100644 --- a/README.md +++ b/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 -time-slab/RayPool scheduler is implemented next, so it must not yet be used -as a performance measurement for numerical-relativity data. The current +fixture. Movie mode collects all current frame-mesh vertices into a single +SoA ray pool, activates rays as a newest-to-oldest coordinate-time scan reaches +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. diff --git a/src/geodesic.c b/src/geodesic.c index 81bd245..1a554ea 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -1,9 +1,7 @@ #include "geodesic.h" #include -typedef struct { - double x[3], Pi[3], log_alpha_p0; -} State; +typedef GeodesicRayState 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, - const double n[3], State *s) { +int geodesic_initialize_past_ray(const SpacetimeSource *source, + 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; - double t; + State state; 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; - for (unsigned int i = 0; i < config->max_steps; i++) { - if (config->capture_log_alpha_p0 > 0.0 && - s.log_alpha_p0 >= config->capture_log_alpha_p0) { - out.status = RAY_ENDPOINT_CAPTURED; - 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; + 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) == + GEODESIC_ADVANCE_ACTIVE) + out.status = RAY_ENDPOINT_MAX_STEPS; return out; } diff --git a/src/geodesic.h b/src/geodesic.h index 8a31273..5631315 100644 --- a/src/geodesic.h +++ b/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 diff --git a/src/main.c b/src/main.c index e14e1b6..e59f500 100644 --- a/src/main.c +++ b/src/main.c @@ -3,6 +3,7 @@ #include "movie.h" #include "observer_track.h" #include "optics.h" +#include "ray.h" #include "spacetime.h" #include @@ -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) || - render_observer_frame(s, catalog, spacetime, &movie.frames[i].observer, - output_path)) + double *hdr = calloc((size_t)s->width * s->height * 3, sizeof *hdr); + if (hdr == NULL || frame_output_path(output_path, s, movie.frames[i].frame_id)) { + 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; diff --git a/src/movie.c b/src/movie.c index 76c4472..eeb2c43 100644 --- a/src/movie.c +++ b/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}; } diff --git a/src/movie.h b/src/movie.h index 327ad22..4cc17d8 100644 --- a/src/movie.h +++ b/src/movie.h @@ -2,6 +2,7 @@ #define MOVIE_H #include "observer_track.h" +#include "frame.h" #include @@ -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 diff --git a/src/ray.c b/src/ray.c new file mode 100644 index 0000000..7a4aeb4 --- /dev/null +++ b/src/ray.c @@ -0,0 +1,93 @@ +#include "ray.h" + +#include +#include + +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}; +} diff --git a/src/ray.h b/src/ray.h new file mode 100644 index 0000000..0248ef1 --- /dev/null +++ b/src/ray.h @@ -0,0 +1,35 @@ +#ifndef RAY_H +#define RAY_H + +#include "geodesic.h" + +#include +#include + +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