wyj/GR-raytracing
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
4aa5f6666a1709cba3f77ab5ff5d4111ff1d7f00
GR-raytracing/src/main.c
T

557 lines
22 KiB
C
Raw Blame History

#include "catalog.h"
#include "frame.h"
#include "movie.h"
#include "observer_track.h"
#include "optics.h"
#include "ray.h"
#include "spacetime.h"
#include <errno.h>
#include <limits.h>
#include <math.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
typedef struct {
int width, height;
int coarse_cell_pixels;
int draw_mesh;
int psf_direct;
int verbose;
double horizontal_fov_deg, look_ra_deg, look_dec_deg, exposure;
double observer_radius;
double observer_inward_speed;
PointSpreadFunction psf;
PsfKernelCache psf_cache;
const char *catalog_path;
const char *all_sky_catalog_path;
const char *output_path;
#ifdef ENABLE_HDR_DEBUG
const char *hdr_output_path;
#endif
const char *observer_track_path;
const char *frames_dir;
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;
int catalog_load_workers;
} Settings;
static int parse_int(const char *text, int *value) {
char *end;
errno = 0;
long parsed = strtol(text, &end, 10);
if (errno || *end || parsed <= 0 || parsed > 16384)
return -1;
*value = (int)parsed;
return 0;
}
static int parse_double(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value <= 0.0 || *value >= 179.0 ? -1 : 0;
}
static int parse_ra_deg(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < 0.0 || *value >= 360.0 ? -1 : 0;
}
static int parse_dec_deg(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < -90.0 || *value > 90.0 ? -1 : 0;
}
static int parse_positive(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value <= 0.0 ? -1 : 0;
}
static int parse_nonnegative(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || !isfinite(*value) || *value < 0.0 ? -1 : 0;
}
static int parse_speed(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value < 0.0 || *value >= 1.0 ? -1 : 0;
}
static int parse_moffat_beta(const char *text, double *value) {
char *end;
errno = 0;
*value = strtod(text, &end);
return errno || *end || *value <= 1.0 ? -1 : 0;
}
static int validate_tonemapped_output_path(const char *path) {
const size_t path_length = strlen(path);
#ifdef ENABLE_PNG
if (path_length >= 4 && strcmp(path + path_length - 4, ".png") == 0)
return 0;
fputs("PNG output is enabled; use a .png output path.\n", stderr);
#else
if (path_length >= 4 && strcmp(path + path_length - 4, ".ppm") == 0)
return 0;
fputs("PNG output is unavailable; use a .ppm output path or rebuild with libpng.\n",
stderr);
#endif
return -1;
}
static int parse_args(int argc, char **argv, Settings *s,
const char **write_path) {
#ifdef ENABLE_PNG
const char *default_output_path = "output/imgs/minkowski_sky.png";
#else
const char *default_output_path = "output/imgs/minkowski_sky.ppm";
#endif
*s = (Settings){.width = 1280,
.height = 720,
.coarse_cell_pixels = 16,
.horizontal_fov_deg = 30.0,
.look_ra_deg = 270.0,
.look_dec_deg = 0.0,
.exposure = 1e-3,
.observer_radius = 30.0,
.psf = {2.7, 4.5},
.catalog_path = "assets/sky_grid_5deg.csv",
.output_path = default_output_path,
.frames_prefix = "frame",
.movie_duration = 2.0,
.movie_fps = 30.0,
.slab_duration = 64.0,
.minkowski_proper_acceleration = 1.52,
.catalog_load_workers = 4};
*write_path = NULL;
for (int i = 1; i < argc; ++i) {
if (!strcmp(argv[i], "--catalog") && i + 1 < argc)
s->catalog_path = argv[++i];
else if (!strcmp(argv[i], "--all-sky-catalog") && i + 1 < argc)
s->all_sky_catalog_path = argv[++i];
else if (!strcmp(argv[i], "--output") && i + 1 < argc)
s->output_path = argv[++i];
#ifdef ENABLE_HDR_DEBUG
else if (!strcmp(argv[i], "--hdr-output") && i + 1 < argc)
s->hdr_output_path = argv[++i];
#endif
else if (!strcmp(argv[i], "--width") && i + 1 < argc &&
!parse_int(argv[++i], &s->width)) {
} else if (!strcmp(argv[i], "--height") && i + 1 < argc &&
!parse_int(argv[++i], &s->height)) {
} else if (!strcmp(argv[i], "--coarse-cell-pixels") && i + 1 < argc &&
!parse_int(argv[++i], &s->coarse_cell_pixels)) {
} else if (!strcmp(argv[i], "--draw-mesh")) {
s->draw_mesh = 1;
} else if (!strcmp(argv[i], "--fov-deg") && i + 1 < argc &&
!parse_double(argv[++i], &s->horizontal_fov_deg)) {
} else if (!strcmp(argv[i], "--look-ra-deg") && i + 1 < argc &&
!parse_ra_deg(argv[++i], &s->look_ra_deg)) {
} else if (!strcmp(argv[i], "--look-dec-deg") && i + 1 < argc &&
!parse_dec_deg(argv[++i], &s->look_dec_deg)) {
} else if (!strcmp(argv[i], "--exposure") && i + 1 < argc &&
!parse_positive(argv[++i], &s->exposure)) {
} else if (!strcmp(argv[i], "--observer-inward-speed") && i + 1 < argc &&
!parse_speed(argv[++i], &s->observer_inward_speed)) {
} else if (!strcmp(argv[i], "--observer-radius") && i + 1 < argc &&
!parse_positive(argv[++i], &s->observer_radius)) {
} else if (!strcmp(argv[i], "--psf-fwhm-pixels") && i + 1 < argc &&
!parse_positive(argv[++i], &s->psf.fwhm_pixels)) {
} else if (!strcmp(argv[i], "--psf-moffat-beta") && i + 1 < argc &&
!parse_moffat_beta(argv[++i], &s->psf.moffat_beta)) {
} else if (!strcmp(argv[i], "--psf-direct")) {
s->psf_direct = 1;
} else if (!strcmp(argv[i], "--verbose")) {
s->verbose = 1;
} else if (!strcmp(argv[i], "--write-catalog") && i + 1 < argc)
*write_path = argv[++i];
else if (!strcmp(argv[i], "--observer-track") && i + 1 < argc)
s->observer_track_path = argv[++i];
else if (!strcmp(argv[i], "--frames-dir") && i + 1 < argc)
s->frames_dir = argv[++i];
else if (!strcmp(argv[i], "--frames-prefix") && i + 1 < argc)
s->frames_prefix = argv[++i];
else if (!strcmp(argv[i], "--start-time") && i + 1 < argc &&
!parse_nonnegative(argv[++i], &s->movie_start_time)) {
} else if (!strcmp(argv[i], "--duration") && i + 1 < argc &&
!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)) {
} else if (!strcmp(argv[i], "--catalog-load-workers") && i + 1 < argc &&
!parse_int(argv[++i], &s->catalog_load_workers)) {
} else if (!strcmp(argv[i], "--write-minkowski-accel-track") &&
i + 1 < argc)
s->write_minkowski_accel_track_path = argv[++i];
else
return -1;
}
return 0;
}
typedef struct {
int verbose;
size_t frame_id;
double splat_start;
const CatalogPrefetchStats *prefetch;
int catalog_load_workers;
int all_sky_catalog;
} RenderProgress;
static void report_splat_progress(void *context, FrameSplatProgressStage stage,
size_t completed, size_t total) {
RenderProgress *progress = context;
(void)completed;
if (progress == NULL)
return;
if (stage == FRAME_SPLAT_PROGRESS_PREFETCH_END && progress->all_sky_catalog) {
const CatalogPrefetchStats *prefetch = progress->prefetch;
fprintf(stderr,
"Catalog prefetch: %zu requested, %zu newly loaded (%zu stars), "
"%zu unavailable in %.3f s; %d loader workers\n",
prefetch->requested_tiles, prefetch->newly_loaded_tiles,
prefetch->newly_loaded_stars, prefetch->unavailable_tiles,
prefetch->load_seconds, progress->catalog_load_workers);
}
if (!progress->verbose)
return;
switch (stage) {
case FRAME_SPLAT_PROGRESS_PREFETCH_BEGIN:
fprintf(stderr, "Frame %zu: finding and prefetching catalog tiles...\n",
progress->frame_id);
break;
case FRAME_SPLAT_PROGRESS_PREFETCH_END:
fprintf(stderr, "Frame %zu: catalog prefetch finished (%zu candidate tiles).\n",
progress->frame_id, total);
break;
case FRAME_SPLAT_PROGRESS_BEGIN:
progress->splat_start = omp_get_wtime();
fprintf(stderr, "Frame %zu: splatting %zu lens triangles...\n",
progress->frame_id, total);
break;
case FRAME_SPLAT_PROGRESS_END:
fprintf(stderr, "Frame %zu: catalog splatting finished in %.1f s; writing image...\n",
progress->frame_id, omp_get_wtime() - progress->splat_start);
break;
}
}
static void ray_pool_status_counts(const RayPool *rays, size_t *pending,
size_t *active, size_t *terminated,
size_t *failed) {
*pending = *active = *terminated = *failed = 0;
for (size_t i = 0; i < rays->count; ++i)
switch (rays->status[i]) {
case RAY_POOL_PENDING: ++*pending; break;
case RAY_POOL_ACTIVE: ++*active; break;
case RAY_POOL_TERMINATED: ++*terminated; break;
case RAY_POOL_FAILED: ++*failed; break;
}
}
static GeodesicTraceConfig trace_config(void) {
#ifdef SPACETIME_SCHWARZSCHILD
return (GeodesicTraceConfig){.coordinate_time_step = 0.1,
.max_steps = 4096,
.capture_log_alpha_p0 = 8.0};
#else
return (GeodesicTraceConfig){.coordinate_time_step = 1.0,
.max_steps = 2048};
#endif
}
static int default_observer(const Settings *s, ObserverState *observer) {
#ifdef SPACETIME_SCHWARZSCHILD
return observer_inward_schwarzschild_ks_look_at(
1.0, s->observer_radius, s->look_ra_deg, s->look_dec_deg,
s->observer_inward_speed, observer);
#else
*observer = observer_fixed_at_origin_look_at(s->look_ra_deg, s->look_dec_deg);
return 0;
#endif
}
static int render_observer_frame(const Settings *s, StarCatalog *catalog,
const SpacetimeSource *spacetime,
const ObserverState *observer,
const char *output_path) {
const GeodesicTraceConfig trace = trace_config();
FrameLensMesh mesh = {0};
double *hdr = calloc((size_t)s->width * s->height * 3, sizeof *hdr);
if (hdr == NULL ||
frame_lens_mesh_build_coarse(&mesh, s->width, s->height,
s->coarse_cell_pixels,
s->horizontal_fov_deg) ||
frame_lens_mesh_trace(&mesh, spacetime, observer, &trace)) {
frame_lens_mesh_destroy(&mesh);
free(hdr);
return -1;
}
if (s->verbose)
fprintf(stderr, "Frame 0: traced %zu lens vertices; starting catalog render.\n",
mesh.vertex_count);
CatalogPrefetchStats prefetch = {0};
PsfSplatStats psf_stats = {0};
RenderProgress progress = {.verbose = s->verbose,
.frame_id = 0,
.prefetch = &prefetch,
.catalog_load_workers = s->catalog_load_workers,
.all_sky_catalog = catalog->kind == STAR_CATALOG_ALL_SKY};
size_t images = frame_splat_catalog(
&mesh, catalog, hdr, s->width, s->height, s->exposure, &s->psf,
&s->psf_cache, spacetime_limits_render_workers_by_memory(spacetime),
s->catalog_load_workers, &prefetch, &psf_stats,
&(FrameSplatProgress){report_splat_progress, &progress});
if (s->draw_mesh)
frame_draw_mesh(&mesh, hdr, s->width, s->height, 0.5, 0.5);
#ifdef ENABLE_HDR_DEBUG
if (s->hdr_output_path != NULL &&
write_hdr_fits(s->hdr_output_path, hdr, s->width, s->height,
s->horizontal_fov_deg)) {
fprintf(stderr, "Failed to write HDR FITS image: %s\n", s->hdr_output_path);
frame_lens_mesh_destroy(&mesh);
free(hdr);
return -1;
}
#endif
int result = write_tonemapped_image(output_path, hdr, s->width, s->height);
fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n",
images, catalog->count, output_path,
result == 0 ? "ok" : "write failed");
psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr);
frame_lens_mesh_destroy(&mesh);
free(hdr);
return result;
}
static int render_frame(const Settings *s, StarCatalog *catalog,
const SpacetimeSource *spacetime) {
ObserverState observer;
return default_observer(s, &observer) ? -1
: render_observer_frame(s, catalog, spacetime,
&observer, s->output_path);
}
static int frame_output_path(char path[PATH_MAX], const Settings *s,
size_t frame_id) {
#ifdef ENABLE_PNG
const char *extension = "png";
#else
const char *extension = "ppm";
#endif
struct stat st;
if (s->frames_dir == NULL || s->frames_prefix == NULL ||
strchr(s->frames_prefix, '/') != NULL || stat(s->frames_dir, &st) ||
!S_ISDIR(st.st_mode))
return -1;
const int written = snprintf(path, PATH_MAX, "%s/%s_%06zu.%s", s->frames_dir,
s->frames_prefix, frame_id, extension);
return written < 0 || written >= PATH_MAX ? -1 : 0;
}
static int render_movie(const Settings *s, 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) ||
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, &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;
size_t slab_id = 0;
while (ray_pool_has_live(&rays)) {
const double slab_lo = slab_hi - s->slab_duration;
MetricSlab *slab = NULL;
size_t pending_before, active_before, terminated_before, failed_before;
ray_pool_status_counts(&rays, &pending_before, &active_before,
&terminated_before, &failed_before);
if (s->verbose)
fprintf(stderr, "Time slab %zu: loading [%.6g, %.6g] with %zu pending and %zu active rays.\n",
slab_id + 1, slab_hi, slab_lo, pending_before, active_before);
if (spacetime_load_slab(spacetime, slab_hi, slab_lo, &slab))
goto done;
ray_pool_activate_in_time_range(&rays, slab);
size_t pending_active, active_active, terminated_active, failed_active;
ray_pool_status_counts(&rays, &pending_active, &active_active,
&terminated_active, &failed_active);
ray_pool_advance_active(&rays, slab, &trace);
spacetime_free_slab(slab);
size_t pending_after, active_after, terminated_after, failed_after;
ray_pool_status_counts(&rays, &pending_after, &active_after,
&terminated_after, &failed_after);
fprintf(stderr,
"Time slab %zu [%.6g, %.6g]: activated %zu; live %zu -> %zu, "
"terminated %zu, failed %zu\n",
++slab_id, slab_hi, slab_lo, active_active - active_before,
pending_before + active_before, pending_after + active_after,
terminated_after, failed_after);
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];
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;
}
CatalogPrefetchStats prefetch = {0};
PsfSplatStats psf_stats = {0};
RenderProgress progress = {
.verbose = s->verbose,
.frame_id = movie.frames[i].frame_id,
.prefetch = &prefetch,
.catalog_load_workers = s->catalog_load_workers,
.all_sky_catalog = catalog->kind == STAR_CATALOG_ALL_SKY};
const size_t images = frame_splat_catalog(
&movie.frames[i].mesh, catalog, hdr, s->width, s->height, s->exposure,
&s->psf, &s->psf_cache,
spacetime_limits_render_workers_by_memory(spacetime),
s->catalog_load_workers, &prefetch, &psf_stats,
&(FrameSplatProgress){report_splat_progress, &progress});
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");
psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr);
if (write_result)
goto done;
}
result = 0;
done:
ray_pool_destroy(&rays);
movie_destroy(&movie);
observer_track_destroy(&track);
return result;
}
static int write_minkowski_accel_track(const Settings *s) {
ObserverTrack track = {0};
const int result = observer_track_generate_minkowski_acceleration(
&track, s->minkowski_proper_acceleration, s->movie_duration,
1.0 / s->movie_fps) ||
observer_track_write_csv(&track,
s->write_minkowski_accel_track_path)
? -1
: 0;
observer_track_destroy(&track);
return result;
}
int main(int argc, char **argv) {
Settings settings;
const char *write_path;
if (parse_args(argc, argv, &settings, &write_path)) {
fprintf(stderr,
"Usage: %s [--catalog PATH | --all-sky-catalog DIR] [--output PATH] [--width N] [--height "
"N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] "
"[--exposure E] [--observer-radius R] [--observer-inward-speed V] "
"[--psf-fwhm-pixels N] [--psf-moffat-beta N] "
"[--psf-direct] [--verbose] "
#ifdef ENABLE_HDR_DEBUG
"[--hdr-output PATH] "
#endif
"[--coarse-cell-pixels N] [--draw-mesh] [--write-catalog PATH] "
"[--catalog-load-workers N] "
"[--observer-track PATH --frames-dir DIR --frames-prefix NAME "
"--start-time T --duration T --fps N] "
"[--proper-acceleration A --write-minkowski-accel-track PATH]\n",
argv[0]);
return 2;
}
if (write_path != NULL)
return catalog_write_octant_grid(write_path) == 0 ? 0
: (perror(write_path), 1);
if (settings.write_minkowski_accel_track_path != NULL)
return write_minkowski_accel_track(&settings) == 0
? 0
: (perror(settings.write_minkowski_accel_track_path), 1);
#ifdef ENABLE_HDR_DEBUG
if (settings.frames_dir != NULL && settings.hdr_output_path != NULL) {
fputs("--hdr-output is available only for a single-frame render.\n", stderr);
return 2;
}
#endif
if (settings.frames_dir == NULL &&
validate_tonemapped_output_path(settings.output_path))
return 2;
StarCatalog catalog = {0};
if (settings.all_sky_catalog_path != NULL) {
if (catalog_load_all_sky(&catalog, settings.all_sky_catalog_path)) {
perror(settings.all_sky_catalog_path);
return 1;
}
} else if (catalog_load_csv(&catalog, settings.catalog_path)) {
if (catalog_write_octant_grid(settings.catalog_path) ||
catalog_load_csv(&catalog, settings.catalog_path)) {
perror(settings.catalog_path);
return 1;
}
fprintf(stderr, "Created test catalog: %s\n", settings.catalog_path);
}
SpacetimeSource spacetime = {0};
if (spacetime_create_default(&spacetime)) {
fputs("Could not create spacetime source\n", stderr);
catalog_destroy(&catalog);
return 1;
}
if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf))
fputs("PSF cache construction failed; using direct evaluator.\n", stderr);
psf_kernel_cache_report_ready(&settings.psf_cache, stderr);
int result = settings.frames_dir != NULL
? render_movie(&settings, &catalog, &spacetime)
: render_frame(&settings, &catalog, &spacetime);
spacetime_destroy(&spacetime);
catalog_destroy(&catalog);
psf_kernel_cache_destroy(&settings.psf_cache);
return result == 0 ? 0 : 1;
}
Reference in new issue View Git Blame Copy Permalink