wyj
2a0197229a
--draw-mesh no longer modifies the primary render. The clean HDR FITS and tone-mapped image are written first, then the final image-plane mesh is overlaid in place on the already-consumed HDR buffer and saved as a <stem>_mesh sibling. This applies uniformly to single-frame, movie-frame, and imported lens-map renders, with derived paths validated before catalog, spacetime, and render initialization. Dummy PSF backends still write nothing and report ignoring --draw-mesh. The reference-image target now uses the configured image extension so the ENABLE_PNG=0 test path passes, and the CLI tests cover clean-main identity, sibling naming, PATH_MAX overflow, and zero-tolerance HDR comparison. README, README.zh-CN, and usage.md describe the new naming; the mesh reference asset is renamed to match.
1346 lines
57 KiB
C
1346 lines
57 KiB
C
#include "catalog.h"
|
|
#include "frame.h"
|
|
#include "lens_map.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 <stdint.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 max_magnification;
|
|
double max_cache_psf_flux;
|
|
double psf_relative_tail;
|
|
double psf_min_y;
|
|
double observer_radius;
|
|
double observer_position[3], observer_velocity[3], camera_roll_deg;
|
|
int position_specified, velocity_specified, look_specified;
|
|
int radius_specified, roll_specified;
|
|
PointSpreadFunction psf;
|
|
PsfKernelCache psf_cache;
|
|
int fast_mode;
|
|
int fast_supersample;
|
|
FastPsfDeposit fast_deposit;
|
|
FastPsfAccumulator *fast_psf;
|
|
const char *catalog_path;
|
|
const char *all_sky_catalog_path;
|
|
const char *output_path;
|
|
const char *lens_map_input_path;
|
|
const char *lens_map_output_path;
|
|
int width_specified, height_specified, fov_specified;
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
int write_hdr_output;
|
|
char hdr_output_path[PATH_MAX];
|
|
#endif
|
|
const char *observer_track_path;
|
|
int movie_track_samples;
|
|
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;
|
|
const char *blackbody_table_path;
|
|
RefinementConfig refinement;
|
|
} 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_nonnegative_int(const char *text, unsigned int *value) {
|
|
char *end;
|
|
errno = 0;
|
|
unsigned long parsed = strtoul(text, &end, 10);
|
|
if (errno || *end || parsed > 1024)
|
|
return -1;
|
|
*value = (unsigned 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 == text || *end || !isfinite(*value) || *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 == text || *end || !isfinite(*value) || *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_finite(const char *text, double *value) {
|
|
char *end;
|
|
errno = 0;
|
|
*value = strtod(text, &end);
|
|
return errno || end == text || *end || !isfinite(*value) ? -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 parse_at_least_one(const char *text, double *value) {
|
|
char *end;
|
|
errno = 0;
|
|
*value = strtod(text, &end);
|
|
return errno || *end || !isfinite(*value) || *value < 1.0 ? -1 : 0;
|
|
}
|
|
|
|
static int parse_finite_positive(const char *text, double *value) {
|
|
char *end;
|
|
errno = 0;
|
|
*value = strtod(text, &end);
|
|
return errno || end == text || *end || !isfinite(*value) || *value <= 0.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;
|
|
}
|
|
|
|
/* Length of the path up to (but excluding) its final extension. Only the
|
|
* last path component is inspected, so dots inside directory names do not
|
|
* count, and a leading-dot basename such as ".png" has no extension. */
|
|
static size_t final_stem_length(const char *path) {
|
|
const char *filename = strrchr(path, '/');
|
|
const char *extension;
|
|
|
|
filename = filename == NULL ? path : filename + 1;
|
|
extension = strrchr(filename, '.');
|
|
return extension != NULL && extension != filename ? (size_t)(extension - path)
|
|
: strlen(path);
|
|
}
|
|
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
static int make_hdr_output_path(char path[PATH_MAX], const char *output_path) {
|
|
static const char hdr_suffix[] = "_HDR.fits";
|
|
const size_t stem_length = final_stem_length(output_path);
|
|
|
|
if (stem_length + sizeof hdr_suffix > PATH_MAX)
|
|
return -1;
|
|
memcpy(path, output_path, stem_length);
|
|
memcpy(path + stem_length, hdr_suffix, sizeof hdr_suffix);
|
|
return 0;
|
|
}
|
|
#endif
|
|
|
|
/* Derives the mesh-overlay sibling of a tone-mapped output path by inserting
|
|
* "_mesh" before the final extension: image.png -> image_mesh.png. */
|
|
static int make_mesh_output_path(char path[PATH_MAX], const char *output_path) {
|
|
static const char mesh_infix[] = "_mesh";
|
|
const size_t stem_length = final_stem_length(output_path);
|
|
const char *suffix = output_path + stem_length;
|
|
const size_t suffix_length = strlen(suffix);
|
|
const size_t infix_length = sizeof mesh_infix - 1;
|
|
|
|
if (stem_length + infix_length + suffix_length + 1 > PATH_MAX)
|
|
return -1;
|
|
memcpy(path, output_path, stem_length);
|
|
memcpy(path + stem_length, mesh_infix, infix_length);
|
|
memcpy(path + stem_length + infix_length, suffix, suffix_length + 1);
|
|
return 0;
|
|
}
|
|
|
|
typedef struct {
|
|
const char *output_path;
|
|
char mesh_path[PATH_MAX];
|
|
int draw_mesh;
|
|
} FrameOutputPaths;
|
|
|
|
/* Fills in the plain and mesh-overlay output paths for one frame, returning
|
|
* -1 if a requested mesh sibling would overflow PATH_MAX. Callers choose the
|
|
* timing; main() separately pre-validates the single-frame path before
|
|
* expensive initialization. */
|
|
static int build_frame_output_paths(const Settings *s, const char *output_path,
|
|
FrameOutputPaths *paths) {
|
|
paths->output_path = output_path;
|
|
paths->draw_mesh = s->draw_mesh;
|
|
if (paths->draw_mesh && make_mesh_output_path(paths->mesh_path, output_path)) {
|
|
fprintf(stderr, "Mesh overlay output path is too long: %s\n", output_path);
|
|
return -1;
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
/* Canonical output order for every frame: clean HDR, clean tone-mapped image,
|
|
* then the mesh overlay. The overlay reuses the already-consumed HDR buffer,
|
|
* so no second full-size framebuffer is allocated and nothing is re-rendered. */
|
|
static int write_frame_outputs(const Settings *s, const FrameLensMesh *mesh,
|
|
double *hdr, int width, int height, double fov_deg,
|
|
const FrameOutputPaths *paths, size_t images,
|
|
size_t stars, const char *note) {
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
if (s->write_hdr_output &&
|
|
write_hdr_fits(s->hdr_output_path, hdr, width, height, fov_deg)) {
|
|
fprintf(stderr, "Failed to write HDR FITS image: %s\n", s->hdr_output_path);
|
|
return -1;
|
|
}
|
|
#else
|
|
(void)s;
|
|
(void)fov_deg;
|
|
#endif
|
|
const int write_result =
|
|
write_tonemapped_image(paths->output_path, hdr, width, height);
|
|
fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s%s)\n",
|
|
images, stars, paths->output_path,
|
|
write_result == 0 ? "ok" : "write failed", note);
|
|
if (write_result)
|
|
return -1;
|
|
if (paths->draw_mesh) {
|
|
frame_draw_mesh(mesh, hdr, width, height, 0.5, 0.5);
|
|
if (write_tonemapped_image(paths->mesh_path, hdr, width, height)) {
|
|
fprintf(stderr, "Failed to write mesh overlay image: %s\n",
|
|
paths->mesh_path);
|
|
return -1;
|
|
}
|
|
fprintf(stderr, "Wrote mesh overlay image: %s\n", paths->mesh_path);
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
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,
|
|
#ifdef GR_DEBUG
|
|
/* Debug builds make existing low-frequency render progress
|
|
* visible by default. Their worker heartbeats are local
|
|
* 8/16/32... milestones, not shared progress counters. */
|
|
.verbose = 1,
|
|
#endif
|
|
.horizontal_fov_deg = 30.0,
|
|
/* Preserve the legacy -Z view after switching to standard
|
|
* right-handed ICRS Cartesian axes. */
|
|
.look_ra_deg = 90.0,
|
|
.look_dec_deg = -90.0,
|
|
.exposure = 1e-3,
|
|
.max_magnification = INFINITY,
|
|
.max_cache_psf_flux = 1.0,
|
|
.psf_relative_tail = 1e-8,
|
|
.psf_min_y = 0.0,
|
|
.observer_radius = 30.0,
|
|
.psf = {2.7, 4.5},
|
|
.fast_supersample = 2,
|
|
.fast_deposit = FAST_PSF_DEPOSIT_NEAREST,
|
|
.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,
|
|
.refinement = {.angle_absolute_rad =
|
|
1e-3 * 3.14159265358979323846 / 180.0,
|
|
.angle_relative = 0.1,
|
|
.jacobian_minimum = 1e-3,
|
|
.min_edge_pixels = 0.5,
|
|
.min_area_pixels2 = 0.25}};
|
|
*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];
|
|
else if (!strcmp(argv[i], "--lens-map-input") && i + 1 < argc)
|
|
s->lens_map_input_path = argv[++i];
|
|
else if (!strcmp(argv[i], "--lens-map-output") && i + 1 < argc)
|
|
s->lens_map_output_path = argv[++i];
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
else if (!strcmp(argv[i], "--hdr-output"))
|
|
s->write_hdr_output = 1;
|
|
#endif
|
|
else if (!strcmp(argv[i], "--width") && i + 1 < argc &&
|
|
!parse_int(argv[++i], &s->width)) {
|
|
s->width_specified = 1;
|
|
} else if (!strcmp(argv[i], "--height") && i + 1 < argc &&
|
|
!parse_int(argv[++i], &s->height)) {
|
|
s->height_specified = 1;
|
|
} else if (!strcmp(argv[i], "--coarse-cell-pixels") && i + 1 < argc &&
|
|
!parse_int(argv[++i], &s->coarse_cell_pixels)) {
|
|
} else if (!strcmp(argv[i], "--refine-max-level") && i + 1 < argc &&
|
|
!parse_nonnegative_int(argv[++i], &s->refinement.max_level)) {
|
|
} else if (!strcmp(argv[i], "--refine-angle-abs-deg") && i + 1 < argc &&
|
|
!parse_positive(argv[++i], &s->refinement.angle_absolute_rad)) {
|
|
s->refinement.angle_absolute_rad *= 3.14159265358979323846 / 180.0;
|
|
} else if (!strcmp(argv[i], "--refine-angle-rel") && i + 1 < argc &&
|
|
!parse_positive(argv[++i], &s->refinement.angle_relative)) {
|
|
} else if (!strcmp(argv[i], "--refine-jacobian-min") && i + 1 < argc &&
|
|
!parse_positive(argv[++i], &s->refinement.jacobian_minimum)) {
|
|
} else if (!strcmp(argv[i], "--refine-min-edge-pixels") && i + 1 < argc &&
|
|
!parse_positive(argv[++i], &s->refinement.min_edge_pixels)) {
|
|
} else if (!strcmp(argv[i], "--refine-min-area-pixels2") && i + 1 < argc &&
|
|
!parse_positive(argv[++i], &s->refinement.min_area_pixels2)) {
|
|
} 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)) {
|
|
s->fov_specified = 1;
|
|
} else if (!strcmp(argv[i], "--look-ra-deg") && i + 1 < argc &&
|
|
!parse_ra_deg(argv[++i], &s->look_ra_deg)) {
|
|
s->look_specified = 1;
|
|
} else if (!strcmp(argv[i], "--look-dec-deg") && i + 1 < argc &&
|
|
!parse_dec_deg(argv[++i], &s->look_dec_deg)) {
|
|
s->look_specified = 1;
|
|
} else if (!strcmp(argv[i], "--exposure") && i + 1 < argc &&
|
|
!parse_positive(argv[++i], &s->exposure)) {
|
|
} else if ((!strcmp(argv[i], "--observer-position") ||
|
|
!strcmp(argv[i], "--observer-velocity")) && i + 3 < argc) {
|
|
const int position = !strcmp(argv[i], "--observer-position");
|
|
double *v = position ? s->observer_position : s->observer_velocity;
|
|
for (int component = 0; component < 3; ++component)
|
|
if (parse_finite(argv[++i], &v[component])) return -1;
|
|
if (position) s->position_specified = 1;
|
|
else s->velocity_specified = 1;
|
|
} else if (!strcmp(argv[i], "--camera-roll-deg") && i + 1 < argc &&
|
|
!parse_finite(argv[++i], &s->camera_roll_deg)) {
|
|
s->roll_specified = 1;
|
|
} else if (!strcmp(argv[i], "--observer-radius") && i + 1 < argc &&
|
|
!parse_finite_positive(argv[++i], &s->observer_radius)) {
|
|
s->radius_specified = 1;
|
|
} 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], "--max-magnification") && i + 1 < argc &&
|
|
!parse_finite_positive(argv[++i], &s->max_magnification)) {
|
|
} else if (!strcmp(argv[i], "--max-cache-psf-flux") && i + 1 < argc &&
|
|
!parse_at_least_one(argv[++i], &s->max_cache_psf_flux)) {
|
|
} else if (!strcmp(argv[i], "--psf-relative-tail") && i + 1 < argc &&
|
|
!parse_finite_positive(argv[++i], &s->psf_relative_tail) &&
|
|
s->psf_relative_tail < 1.0) {
|
|
} else if (!strcmp(argv[i], "--psf-min-y") && i + 1 < argc &&
|
|
!parse_nonnegative(argv[++i], &s->psf_min_y)) {
|
|
} else if (!strcmp(argv[i], "--psf-direct")) {
|
|
s->psf_direct = 1;
|
|
} else if (!strcmp(argv[i], "--fast-mode")) {
|
|
s->fast_mode = 1;
|
|
} else if (!strcmp(argv[i], "--fast-supersample") && i + 1 < argc &&
|
|
!parse_int(argv[++i], &s->fast_supersample)) {
|
|
} else if (!strcmp(argv[i], "--fast-deposit") && i + 1 < argc) {
|
|
const char *mode = argv[++i];
|
|
if (!strcmp(mode, "nearest"))
|
|
s->fast_deposit = FAST_PSF_DEPOSIT_NEAREST;
|
|
else if (!strcmp(mode, "bilinear"))
|
|
s->fast_deposit = FAST_PSF_DEPOSIT_BILINEAR;
|
|
else
|
|
return -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], "--movie-track-samples"))
|
|
s->movie_track_samples = 1;
|
|
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], "--blackbody-table") && i + 1 < argc) {
|
|
s->blackbody_table_path = argv[++i];
|
|
} 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;
|
|
}
|
|
|
|
static void print_help(const char *program) {
|
|
#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
|
|
|
|
printf("Usage: %s [options]\n\n", program);
|
|
fputs("Input/output:\n"
|
|
" --catalog PATH Single CSV catalog (default: assets/sky_grid_5deg.csv)\n"
|
|
" --all-sky-catalog DIR Lazy-loaded all-sky catalog directory (default: disabled)\n"
|
|
" --output PATH Tonemapped image path (default: ", stdout);
|
|
printf("%s)\n", default_output_path);
|
|
fputs(" --lens-map-output FILE Save finalized ray-traced lens map; keeps rendering normally\n"
|
|
" --lens-map-input FILE Re-render a saved map without observer, spacetime, or ray tracing\n"
|
|
" (map width, height, and FOV are fixed; mutually exclusive with --lens-map-output)\n",
|
|
stdout);
|
|
fputs(" --write-catalog PATH Write the generated test catalog and exit (default: disabled)\n"
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
" --hdr-output Also write a single-frame linear HDR FITS file (default: disabled)\n"
|
|
#endif
|
|
"\nImage and camera:\n"
|
|
" --width N Image width in pixels (default: 1280)\n"
|
|
" --height N Image height in pixels (default: 720)\n"
|
|
" --fov-deg D Horizontal field of view in degrees (default: 30)\n"
|
|
" --look-ra-deg D ICRS look direction right ascension in degrees (default: 90)\n"
|
|
" --look-dec-deg D ICRS look direction declination in degrees (default: -90)\n"
|
|
" --exposure E Linear exposure multiplier (default: 1e-3)\n"
|
|
" --observer-position X Y Z Coordinate position; alone implies looking at the origin\n"
|
|
" --observer-radius R Infer position = -R * look direction (default R: 30); conflicts with position\n"
|
|
" --observer-velocity VX VY VZ Coordinate dx/dt, dy/dt, dz/dt (default: 0 0 0); must be timelike\n"
|
|
" --camera-roll-deg ANGLE Rotate up toward right about forward (default: 0)\n"
|
|
" Explicit look alone implies position = -R * look direction.\n"
|
|
" Look is projected into the moving camera rest space.\n"
|
|
#ifdef SPACETIME_SCHWARZSCHILD
|
|
" Default position: (0,0,30); look RA=90, Dec=-90.\n"
|
|
#else
|
|
" Default position: (0,0,0); look RA=90, Dec=-90.\n"
|
|
#endif
|
|
" Single-frame camera options conflict with track/map input.\n"
|
|
"\nPSF and catalog splatting:\n"
|
|
" --psf-fwhm-pixels N Moffat PSF FWHM in pixels (default: 2.7)\n"
|
|
" --psf-moffat-beta N Moffat PSF beta, greater than 1 (default: 4.5)\n"
|
|
" --max-magnification M Preview magnification cap (default: unlimited)\n"
|
|
" --max-cache-psf-flux F Cache flux limit before direct fallback (default: 1)\n"
|
|
" --psf-relative-tail R Maximum omitted relative PSF tail fraction (default: 1e-8)\n"
|
|
" --psf-min-y Y Skip events below this linear HDR luminance (default: 0, disabled)\n"
|
|
" --psf-direct Disable the PSF lookup cache (default: disabled)\n"
|
|
" --fast-mode Deposit each image as a supersampled delta and run one\n"
|
|
" global PSF convolution + downsample (CPU only; preview)\n"
|
|
" --fast-supersample N Fast-mode supersample factor, 1..8 (default: 2)\n"
|
|
" --fast-deposit MODE nearest (exact FWHM/beta, 0.5/N px quantization) or\n"
|
|
" bilinear (exact centroid, broadens FWHM); default: nearest\n"
|
|
" --catalog-load-workers N All-sky catalog loader workers (default: 4)\n"
|
|
" --blackbody-table FILE Explicit GRBBLUT3 table\n"
|
|
" (default: assets/blackbody/cie1931_2deg_xyz_1024.grbblut)\n",
|
|
stdout);
|
|
fputs(
|
|
"\nAdaptive lens mesh:\n"
|
|
" --coarse-cell-pixels N Initial mesh cell size in pixels (default: 16)\n"
|
|
" --refine-max-level N Maximum refinement level (default: 0)\n"
|
|
" --refine-angle-abs-deg D Absolute angular interpolation error limit (default: 0.001)\n"
|
|
" --refine-angle-rel R Relative angular interpolation error limit (default: 0.1)\n"
|
|
" --refine-jacobian-min J Fold-refinement Jacobian threshold (default: 1e-3)\n"
|
|
" --refine-min-edge-pixels P Stop refinement below this edge length (default: 0.5)\n"
|
|
" --refine-min-area-pixels2 A Stop refinement below this triangle area (default: 0.25)\n",
|
|
stdout);
|
|
#ifdef ENABLE_PNG
|
|
fputs(" --draw-mesh Also write the final lens-mesh overlay as <output-stem>_mesh.png\n",
|
|
stdout);
|
|
#else
|
|
fputs(" --draw-mesh Also write the final lens-mesh overlay as <output-stem>_mesh.ppm\n",
|
|
stdout);
|
|
#endif
|
|
fputs(
|
|
"\nMovie and observer track:\n"
|
|
" --observer-track PATH Observer worldline/tetrad CSV for movie rendering (default: disabled)\n"
|
|
" --movie-track-samples One frame per CSV row; ignores start-time/duration/fps (default: disabled)\n"
|
|
" --frames-dir DIR Write a movie image sequence to this directory (default: disabled)\n"
|
|
" --frames-prefix NAME Movie frame filename prefix (default: frame)\n"
|
|
" --start-time T Movie start coordinate time (default: 0)\n"
|
|
" --duration T Movie duration (default: 2)\n"
|
|
" --fps N Movie frame rate (default: 30)\n"
|
|
" --slab-duration T Metric time-slab duration (default: 64)\n"
|
|
" --proper-acceleration A Minkowski accelerated-track proper acceleration (default: 1.52)\n"
|
|
" --write-minkowski-accel-track PATH\n"
|
|
" Write a generated Minkowski acceleration track and exit (default: disabled)\n"
|
|
"\nOther:\n"
|
|
" --verbose Print low-frequency rendering progress (default: disabled)\n"
|
|
" --help Print this help and exit\n",
|
|
stdout);
|
|
}
|
|
|
|
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:
|
|
#ifdef PSF_BACKEND_DUMMY
|
|
fprintf(stderr, "Frame %zu: catalog classification finished in %.1f s; reporting chunk statistics...\n",
|
|
progress->frame_id, omp_get_wtime() - progress->splat_start);
|
|
#else
|
|
fprintf(stderr, "Frame %zu: catalog splatting finished in %.1f s; writing image...\n",
|
|
progress->frame_id, omp_get_wtime() - progress->splat_start);
|
|
#endif
|
|
break;
|
|
}
|
|
}
|
|
|
|
static void report_splat_worker_progress(void *context, size_t worker_id,
|
|
size_t worker_count,
|
|
size_t triangle_count, int finished) {
|
|
const RenderProgress *progress = context;
|
|
if (progress == NULL || !progress->verbose)
|
|
return;
|
|
if (finished)
|
|
fprintf(stderr, "Frame %zu: splat worker %zu/%zu finished after %zu local triangles.\n",
|
|
progress->frame_id, worker_id + 1, worker_count, triangle_count);
|
|
else if (triangle_count == 0)
|
|
fprintf(stderr, "Frame %zu: splat worker %zu/%zu started.\n",
|
|
progress->frame_id, worker_id + 1, worker_count);
|
|
else
|
|
fprintf(stderr, "Frame %zu: splat worker %zu/%zu reached %zu local triangles.\n",
|
|
progress->frame_id, worker_id + 1, worker_count, triangle_count);
|
|
}
|
|
|
|
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 void report_frame_refinement(void *context, size_t generation,
|
|
size_t sample_count, size_t vertex_count,
|
|
size_t triangle_count, int added_vertices,
|
|
int finished) {
|
|
const Settings *settings = context;
|
|
if (settings == NULL || !settings->verbose)
|
|
return;
|
|
if (!finished)
|
|
fprintf(stderr,
|
|
"Frame 0: refinement generation %zu tracing %zu samples "
|
|
"from %zu vertices and %zu triangles.\n",
|
|
generation, sample_count, vertex_count, triangle_count);
|
|
else
|
|
fprintf(stderr,
|
|
"Frame 0: refinement generation %zu finished; added %d vertices, "
|
|
"now %zu vertices and %zu triangles.\n",
|
|
generation, added_vertices, vertex_count, triangle_count);
|
|
}
|
|
|
|
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 resolve_camera(Settings *s) {
|
|
if (s->movie_track_samples && (!s->observer_track_path || !s->frames_dir ||
|
|
s->lens_map_input_path)) {
|
|
fputs("--movie-track-samples requires --observer-track and --frames-dir, without --lens-map-input.\n", stderr);
|
|
return -1;
|
|
}
|
|
const int camera_specified = s->position_specified || s->look_specified ||
|
|
s->radius_specified || s->velocity_specified || s->roll_specified;
|
|
if (s->position_specified && s->radius_specified) {
|
|
fputs("--observer-position and --observer-radius are mutually exclusive.\n", stderr);
|
|
return -1;
|
|
}
|
|
if (camera_specified && (s->observer_track_path || s->frames_dir ||
|
|
s->lens_map_input_path)) {
|
|
fputs("Single-frame camera options cannot be combined with movie/observer-track or --lens-map-input.\n", stderr);
|
|
return -1;
|
|
}
|
|
if (s->position_specified && !s->look_specified) {
|
|
const double *x = s->observer_position;
|
|
const double radius = hypot(hypot(x[0], x[1]), x[2]);
|
|
if (!isfinite(radius) || radius == 0.0) {
|
|
fputs("Cannot infer a look direction from this position; specify --look-ra-deg and/or --look-dec-deg.\n", stderr);
|
|
return -1;
|
|
}
|
|
const double degrees = 180.0 / 3.14159265358979323846;
|
|
s->look_ra_deg = (x[0] == 0.0 && x[1] == 0.0)
|
|
? 0.0 : atan2(-x[1], -x[0]) * degrees;
|
|
if (s->look_ra_deg < 0.0) s->look_ra_deg += 360.0;
|
|
if (s->look_ra_deg >= 360.0) s->look_ra_deg = 0.0;
|
|
s->look_dec_deg = atan2(-x[2], hypot(x[0], x[1])) * degrees;
|
|
}
|
|
int infer_position = s->look_specified || s->radius_specified;
|
|
#ifdef SPACETIME_SCHWARZSCHILD
|
|
infer_position = 1;
|
|
#endif
|
|
if (!s->position_specified && infer_position) {
|
|
const ObserverState pointing = observer_fixed_at_origin_look_at(
|
|
s->look_ra_deg, s->look_dec_deg);
|
|
for (int i = 0; i < 3; ++i)
|
|
s->observer_position[i] = -s->observer_radius * pointing.tetrad[1][i + 1];
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
static int build_observer(const Settings *s, const SpacetimeSource *spacetime,
|
|
ObserverState *observer) {
|
|
ObserverCamera camera = {.look_ra_deg = s->look_ra_deg,
|
|
.look_dec_deg = s->look_dec_deg,
|
|
.roll_deg = s->camera_roll_deg};
|
|
for (int i = 0; i < 3; ++i) {
|
|
camera.position[i] = s->observer_position[i];
|
|
camera.velocity[i] = s->observer_velocity[i];
|
|
}
|
|
if (spacetime_classify(spacetime, camera.coordinate_time, camera.position) ==
|
|
SPACETIME_RAY_CAPTURED) {
|
|
fputs("Camera position is inside the backend capture cutoff or invalid.\n", stderr);
|
|
return -1;
|
|
}
|
|
MetricData metric;
|
|
if (spacetime_eval(spacetime, camera.coordinate_time, camera.position, &metric)) {
|
|
fputs("Could not evaluate metric at the camera event.\n", stderr);
|
|
return -1;
|
|
}
|
|
double q;
|
|
const ObserverBuildResult result = observer_from_coordinate_camera(
|
|
&metric, &camera, observer, &q);
|
|
if (result != OBSERVER_BUILD_OK) {
|
|
fprintf(stderr, "Camera construction failed (%s): position=(%.17g, %.17g, %.17g), "
|
|
"coordinate velocity=(%.17g, %.17g, %.17g), Q=%.17g.\n",
|
|
result == OBSERVER_BUILD_NON_TIMELIKE ? "velocity is not timelike" :
|
|
result == OBSERVER_BUILD_INVALID_INPUT ? "invalid input/metric" :
|
|
"invalid tetrad",
|
|
camera.position[0], camera.position[1], camera.position[2],
|
|
camera.velocity[0], camera.velocity[1], camera.velocity[2], q);
|
|
return -1;
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
static void report_psf_splat(const Settings *s, const PsfSplatStats *stats) {
|
|
if (s->fast_mode) {
|
|
fprintf(stderr,
|
|
"Fast PSF splats: deposited %zu, wing-clipped %zu, discarded "
|
|
"below min-Y %zu\n",
|
|
stats->cached_splats, stats->cached_wing_clipped,
|
|
stats->discarded_below_min_y);
|
|
return;
|
|
}
|
|
psf_kernel_cache_report(&s->psf_cache, stats, stderr);
|
|
}
|
|
|
|
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_destroy(&mesh);
|
|
free(hdr);
|
|
return -1;
|
|
}
|
|
if (s->verbose)
|
|
fprintf(stderr, "Frame 0: tracing %zu initial rays from %zu mesh triangles...\n",
|
|
mesh.vertex_count, mesh.triangle_count);
|
|
const double initial_trace_start = omp_get_wtime();
|
|
if (frame_lens_mesh_trace(&mesh, spacetime, observer, &trace)) {
|
|
frame_lens_mesh_destroy(&mesh);
|
|
free(hdr);
|
|
return -1;
|
|
}
|
|
if (s->verbose)
|
|
fprintf(stderr, "Frame 0: initial ray trace finished in %.3f s.\n",
|
|
omp_get_wtime() - initial_trace_start);
|
|
if (s->refinement.max_level > 0 && s->verbose)
|
|
fprintf(stderr, "Frame 0: starting adaptive ray-trace refinement (max level %u)...\n",
|
|
s->refinement.max_level);
|
|
const double refinement_start = omp_get_wtime();
|
|
if (frame_lens_mesh_refine_with_progress(
|
|
&mesh, spacetime, observer, &trace, &s->refinement,
|
|
s->verbose ? report_frame_refinement : NULL, (void *)s)) {
|
|
frame_lens_mesh_destroy(&mesh);
|
|
free(hdr);
|
|
return -1;
|
|
}
|
|
if (s->refinement.max_level > 0 && s->verbose)
|
|
fprintf(stderr,
|
|
"Frame 0: adaptive ray-trace refinement finished in %.3f s; "
|
|
"%zu vertices, %zu triangles.\n",
|
|
omp_get_wtime() - refinement_start, mesh.vertex_count,
|
|
mesh.triangle_count);
|
|
if (s->lens_map_output_path != NULL) {
|
|
const LensMapFrame map_frame = {.frame_id = 0,
|
|
.coordinate_time = 0.0,
|
|
.proper_time = 0.0,
|
|
.mesh = mesh};
|
|
if (lens_map_write(s->lens_map_output_path, s->width, s->height,
|
|
s->horizontal_fov_deg, &map_frame, 1)) {
|
|
fprintf(stderr, "Failed to write lens map: %s\n", s->lens_map_output_path);
|
|
frame_lens_mesh_destroy(&mesh);
|
|
free(hdr);
|
|
return -1;
|
|
}
|
|
fprintf(stderr, "Wrote lens map: %s (%zu vertices, %zu triangles)\n",
|
|
s->lens_map_output_path, mesh.vertex_count, mesh.triangle_count);
|
|
}
|
|
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, s->max_magnification, s->max_cache_psf_flux,
|
|
s->psf_relative_tail,
|
|
s->psf_min_y,
|
|
spacetime_limits_render_workers_by_memory(spacetime),
|
|
s->catalog_load_workers, &prefetch, &psf_stats,
|
|
&(FrameSplatProgress){report_splat_progress,
|
|
s->verbose ? report_splat_worker_progress : NULL,
|
|
&progress},
|
|
s->fast_psf);
|
|
if (images == SIZE_MAX) {
|
|
frame_lens_mesh_destroy(&mesh);
|
|
free(hdr);
|
|
return -1;
|
|
}
|
|
FrameOutputPaths output_paths;
|
|
if (build_frame_output_paths(s, output_path, &output_paths)) {
|
|
frame_lens_mesh_destroy(&mesh);
|
|
free(hdr);
|
|
return -1;
|
|
}
|
|
const int result = write_frame_outputs(
|
|
s, &mesh, hdr, s->width, s->height, s->horizontal_fov_deg, &output_paths,
|
|
images, catalog->count, "");
|
|
report_psf_splat(s, &psf_stats);
|
|
frame_lens_mesh_destroy(&mesh);
|
|
free(hdr);
|
|
return result;
|
|
}
|
|
|
|
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 trace_movie_generation(Movie *movie, const Settings *s,
|
|
const SpacetimeSource *spacetime,
|
|
const GeodesicTraceConfig *trace,
|
|
size_t generation) {
|
|
RayPool rays = {0};
|
|
size_t ray_count = 0;
|
|
size_t total_added = 0;
|
|
for (size_t f = 0; f < movie->frame_count; ++f) {
|
|
const int prepared = frame_lens_mesh_prepare_generation(
|
|
&movie->frames[f].mesh, &s->refinement);
|
|
if (prepared < 0) return -1;
|
|
ray_count += (size_t)prepared;
|
|
}
|
|
if (ray_count == 0) return 0;
|
|
if (s->verbose)
|
|
fprintf(stderr,
|
|
"Ray trace generation %zu: collected %zu new samples across %zu frames.\n",
|
|
generation, ray_count, movie->frame_count);
|
|
if (ray_pool_init(&rays, ray_count)) return -1;
|
|
for (size_t f = 0; f < movie->frame_count; ++f) {
|
|
size_t count = 0;
|
|
const FrameSample *samples = frame_lens_mesh_samples(&movie->frames[f].mesh, &count);
|
|
for (size_t sample = 0; sample < count; ++sample)
|
|
if (ray_pool_append(&rays, &movie->frames[f].observer,
|
|
samples[sample].vertex.camera_direction, f, sample)) {
|
|
ray_pool_destroy(&rays);
|
|
return -1;
|
|
}
|
|
}
|
|
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 (spacetime_load_slab(spacetime, slab_hi, slab_lo, &slab)) {
|
|
ray_pool_destroy(&rays);
|
|
return -1;
|
|
}
|
|
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);
|
|
if (s->verbose)
|
|
fprintf(stderr,
|
|
"Ray trace generation %zu, slab %zu [%.6g, %.6g]: activated %zu; "
|
|
"live %zu -> %zu, terminated %zu, failed %zu.\n",
|
|
generation, ++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)
|
|
if (frame_lens_mesh_install_sample(&movie->frames[rays.frame_id[i]].mesh,
|
|
rays.vertex_id[i], &rays.endpoint[i])) {
|
|
ray_pool_destroy(&rays);
|
|
return -1;
|
|
}
|
|
ray_pool_destroy(&rays);
|
|
if (s->verbose)
|
|
fprintf(stderr, "Ray trace generation %zu: installing endpoints and refining meshes.\n",
|
|
generation);
|
|
for (size_t f = 0; f < movie->frame_count; ++f) {
|
|
/* Frames converge independently; only finish those traced this pass. */
|
|
if (movie->frames[f].mesh.sample_count == 0)
|
|
continue;
|
|
const int added = frame_lens_mesh_finish_generation(&movie->frames[f].mesh,
|
|
&s->refinement);
|
|
if (added < 0) {
|
|
fprintf(stderr, "Ray trace generation %zu: frame %zu refinement failed.\n",
|
|
generation, movie->frames[f].frame_id);
|
|
return -1;
|
|
}
|
|
total_added += (size_t)added;
|
|
if (s->verbose && added > 0)
|
|
fprintf(stderr, "Ray trace generation %zu: frame %zu added %d vertices.\n",
|
|
generation, movie->frames[f].frame_id, added);
|
|
}
|
|
if (s->verbose)
|
|
fprintf(stderr,
|
|
"Ray trace generation %zu: refinement finished; added %zu vertices "
|
|
"across %zu frames.\n",
|
|
generation, total_added, movie->frame_count);
|
|
return 1;
|
|
}
|
|
|
|
static int render_movie(const Settings *s, StarCatalog *catalog,
|
|
const SpacetimeSource *spacetime) {
|
|
ObserverTrack track = {0};
|
|
Movie movie = {0};
|
|
const GeodesicTraceConfig trace = trace_config();
|
|
int result = -1;
|
|
if (s->observer_track_path == NULL ||
|
|
observer_track_load_csv(&track, s->observer_track_path) ||
|
|
(s->movie_track_samples ? movie_init_track_samples(&movie, &track) :
|
|
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;
|
|
for (size_t generation = 0;; ++generation) {
|
|
const int traced = trace_movie_generation(&movie, s, spacetime, &trace,
|
|
generation);
|
|
if (traced < 0)
|
|
goto done;
|
|
if (traced == 0)
|
|
break;
|
|
}
|
|
if (s->lens_map_output_path != NULL) {
|
|
LensMapFrame *map_frames = calloc(movie.frame_count, sizeof *map_frames);
|
|
if (map_frames == NULL) goto done;
|
|
for (size_t i = 0; i < movie.frame_count; ++i)
|
|
map_frames[i] = (LensMapFrame){.frame_id = movie.frames[i].frame_id,
|
|
.coordinate_time = movie.frames[i].coordinate_time,
|
|
.proper_time = movie.frames[i].proper_time,
|
|
.mesh = movie.frames[i].mesh};
|
|
const int write_failed = lens_map_write(s->lens_map_output_path, s->width,
|
|
s->height, s->horizontal_fov_deg,
|
|
map_frames, movie.frame_count);
|
|
free(map_frames);
|
|
if (write_failed) {
|
|
fprintf(stderr, "Failed to write lens map: %s\n", s->lens_map_output_path);
|
|
goto done;
|
|
}
|
|
fprintf(stderr, "Wrote lens-map movie: %s (%zu frames)\n",
|
|
s->lens_map_output_path, movie.frame_count);
|
|
}
|
|
for (size_t i = 0; i < movie.frame_count; ++i) {
|
|
char output_path[PATH_MAX];
|
|
FrameOutputPaths output_paths;
|
|
if (frame_output_path(output_path, s, movie.frames[i].frame_id) ||
|
|
build_frame_output_paths(s, output_path, &output_paths))
|
|
goto done;
|
|
double *hdr = calloc((size_t)s->width * s->height * 3, sizeof *hdr);
|
|
if (hdr == NULL)
|
|
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, s->max_magnification, s->max_cache_psf_flux,
|
|
s->psf_relative_tail,
|
|
s->psf_min_y,
|
|
spacetime_limits_render_workers_by_memory(spacetime),
|
|
s->catalog_load_workers, &prefetch, &psf_stats,
|
|
&(FrameSplatProgress){report_splat_progress,
|
|
s->verbose ? report_splat_worker_progress : NULL,
|
|
&progress},
|
|
s->fast_psf);
|
|
if (images == SIZE_MAX) {
|
|
free(hdr);
|
|
goto done;
|
|
}
|
|
const int write_result = write_frame_outputs(
|
|
s, &movie.frames[i].mesh, hdr, s->width, s->height,
|
|
s->horizontal_fov_deg, &output_paths, images, catalog->count, "");
|
|
free(hdr);
|
|
report_psf_splat(s, &psf_stats);
|
|
if (write_result)
|
|
goto done;
|
|
}
|
|
result = 0;
|
|
done:
|
|
movie_destroy(&movie);
|
|
observer_track_destroy(&track);
|
|
return result;
|
|
}
|
|
|
|
static int render_lens_map(const Settings *s, StarCatalog *catalog) {
|
|
LensMap map = {0};
|
|
if (lens_map_read(s->lens_map_input_path, &map)) {
|
|
fprintf(stderr, "Failed to read or validate lens map: %s\n", s->lens_map_input_path);
|
|
return -1;
|
|
}
|
|
if ((s->width_specified && s->width != map.width) ||
|
|
(s->height_specified && s->height != map.height) ||
|
|
(s->fov_specified && fabs(s->horizontal_fov_deg - map.horizontal_fov_deg) > 1e-12)) {
|
|
fputs("Imported lens-map geometry conflicts with --width, --height, or --fov-deg.\n", stderr);
|
|
lens_map_destroy(&map);
|
|
return -1;
|
|
}
|
|
if ((map.frame_count == 1 && s->frames_dir != NULL) ||
|
|
(map.frame_count > 1 && s->frames_dir == NULL)) {
|
|
fputs("A single-frame lens map uses --output; a multi-frame lens map requires --frames-dir.\n", stderr);
|
|
lens_map_destroy(&map);
|
|
return -1;
|
|
}
|
|
int result = 0;
|
|
FastPsfAccumulator local_fast = {0};
|
|
FastPsfAccumulator *fast = NULL;
|
|
if (s->fast_mode &&
|
|
fast_psf_accumulator_init(&local_fast, map.width, map.height,
|
|
s->fast_supersample, s->fast_deposit, &s->psf,
|
|
s->psf_relative_tail, s->psf_min_y,
|
|
s->psf_direct)) {
|
|
fputs("Fast-mode accumulator construction failed for the imported map.\n",
|
|
stderr);
|
|
lens_map_destroy(&map);
|
|
return -1;
|
|
}
|
|
if (s->fast_mode) {
|
|
fast = &local_fast;
|
|
fast_psf_accumulator_report(&local_fast, stderr);
|
|
fast_psf_accumulator_set_verbose(&local_fast, s->verbose);
|
|
}
|
|
for (size_t i = 0; i < map.frame_count; ++i) {
|
|
const char *output_path = s->output_path;
|
|
char movie_path[PATH_MAX];
|
|
#ifndef PSF_BACKEND_DUMMY
|
|
FrameOutputPaths output_paths;
|
|
#endif
|
|
if (map.frame_count > 1) {
|
|
if (frame_output_path(movie_path, s, (size_t)map.frames[i].frame_id)) {
|
|
fputs("Could not construct imported-map movie output path.\n", stderr);
|
|
result = -1; break;
|
|
}
|
|
output_path = movie_path;
|
|
}
|
|
#ifndef PSF_BACKEND_DUMMY
|
|
if (build_frame_output_paths(s, output_path, &output_paths)) {
|
|
result = -1; break;
|
|
}
|
|
#else
|
|
(void)output_path;
|
|
#endif
|
|
#ifdef PSF_BACKEND_DUMMY
|
|
/* The dummy backend classifies real events but never touches pixels. */
|
|
double *hdr = calloc(1, sizeof *hdr);
|
|
#else
|
|
double *hdr = calloc((size_t)map.width * map.height * 3, sizeof *hdr);
|
|
#endif
|
|
if (hdr == NULL) { result = -1; break; }
|
|
CatalogPrefetchStats prefetch = {0};
|
|
PsfSplatStats psf_stats = {0};
|
|
RenderProgress progress = {.verbose = s->verbose,
|
|
.frame_id = (size_t)map.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(
|
|
&map.frames[i].mesh, catalog, hdr, map.width, map.height, s->exposure,
|
|
&s->psf, &s->psf_cache, s->max_magnification, s->max_cache_psf_flux,
|
|
s->psf_relative_tail, s->psf_min_y, 0, s->catalog_load_workers,
|
|
&prefetch, &psf_stats,
|
|
&(FrameSplatProgress){report_splat_progress,
|
|
s->verbose ? report_splat_worker_progress : NULL,
|
|
&progress},
|
|
fast);
|
|
if (images == SIZE_MAX) {
|
|
free(hdr);
|
|
result = -1;
|
|
break;
|
|
}
|
|
#ifndef PSF_BACKEND_DUMMY
|
|
const int write_result = write_frame_outputs(
|
|
s, &map.frames[i].mesh, hdr, map.width, map.height,
|
|
map.horizontal_fov_deg, &output_paths, images, catalog->count,
|
|
"; imported lens map");
|
|
free(hdr);
|
|
report_psf_splat(s, &psf_stats);
|
|
if (write_result) { result = -1; break; }
|
|
#else
|
|
if (s->draw_mesh)
|
|
fputs("Dummy PSF backend ignores --draw-mesh.\n", stderr);
|
|
free(hdr);
|
|
fprintf(stderr,
|
|
"Dummy PSF classified %zu images from %zu catalog stars; no HDR, PNG, or PPM was written.\n",
|
|
images, catalog->count);
|
|
report_psf_splat(s, &psf_stats);
|
|
#endif
|
|
}
|
|
fast_psf_accumulator_destroy(&local_fast);
|
|
lens_map_destroy(&map);
|
|
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 (argc == 2 && !strcmp(argv[1], "--help")) {
|
|
print_help(argv[0]);
|
|
return 0;
|
|
}
|
|
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] "
|
|
"[--lens-map-input FILE | --lens-map-output FILE] "
|
|
"[--exposure E] [--observer-radius R | --observer-position X Y Z] "
|
|
"[--observer-velocity VX VY VZ] [--camera-roll-deg ANGLE] "
|
|
"[--psf-fwhm-pixels N] [--psf-moffat-beta N] "
|
|
"[--max-magnification M] [--max-cache-psf-flux F] "
|
|
"[--psf-relative-tail R] [--psf-min-y Y] "
|
|
"[--psf-direct] [--fast-mode --fast-supersample N "
|
|
"--fast-deposit nearest|bilinear] [--verbose] "
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
"[--hdr-output] "
|
|
#endif
|
|
"[--coarse-cell-pixels N] [--refine-max-level N "
|
|
"--refine-angle-abs-deg D --refine-angle-rel R "
|
|
"--refine-jacobian-min J "
|
|
"--refine-min-edge-pixels P --refine-min-area-pixels2 A] "
|
|
"[--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 | --movie-track-samples] "
|
|
"[--proper-acceleration A --write-minkowski-accel-track PATH]\n",
|
|
argv[0]);
|
|
return 2;
|
|
}
|
|
#ifdef GR_DEBUG
|
|
fputs("Debug build: low-frequency render progress is enabled; direct PSF "
|
|
"fallbacks log immediately from their splat worker.\n", stderr);
|
|
#endif
|
|
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);
|
|
if (settings.lens_map_input_path != NULL && settings.lens_map_output_path != NULL) {
|
|
fputs("--lens-map-input and --lens-map-output are mutually exclusive.\n", stderr);
|
|
return 2;
|
|
}
|
|
#ifdef PSF_BACKEND_DUMMY
|
|
if (settings.lens_map_input_path == NULL) {
|
|
fputs("The dummy PSF backend requires --lens-map-input and never traces or writes images.\n",
|
|
stderr);
|
|
return 2;
|
|
}
|
|
if (settings.psf_direct) {
|
|
fputs("The dummy PSF backend measures cache-event chunks and does not support --psf-direct.\n",
|
|
stderr);
|
|
return 2;
|
|
}
|
|
fputs("Dummy PSF backend: --output and --hdr-output are accepted for command parity but no image files will be written.\n",
|
|
stderr);
|
|
#endif
|
|
#if defined(PSF_BACKEND_HIP) || defined(PSF_BACKEND_DUMMY)
|
|
if (settings.fast_mode) {
|
|
fputs("--fast-mode is available only in the CPU PSF backend build.\n",
|
|
stderr);
|
|
return 2;
|
|
}
|
|
#endif
|
|
if (settings.fast_mode &&
|
|
(settings.fast_supersample < 1 || settings.fast_supersample > 8)) {
|
|
fputs("--fast-supersample must be between 1 and 8.\n", stderr);
|
|
return 2;
|
|
}
|
|
if (settings.fast_mode)
|
|
fputs("Fast mode is a preview approximation; --max-cache-psf-flux is "
|
|
"ignored and one global kernel is used for every event.\n",
|
|
stderr);
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
if (settings.frames_dir != NULL && settings.write_hdr_output) {
|
|
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;
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
if (settings.write_hdr_output &&
|
|
make_hdr_output_path(settings.hdr_output_path, settings.output_path)) {
|
|
fputs("HDR output path is too long.\n", stderr);
|
|
return 2;
|
|
}
|
|
#endif
|
|
if (settings.draw_mesh && settings.frames_dir == NULL) {
|
|
char mesh_path[PATH_MAX];
|
|
if (make_mesh_output_path(mesh_path, settings.output_path)) {
|
|
fputs("Mesh overlay output path is too long.\n", stderr);
|
|
return 2;
|
|
}
|
|
}
|
|
if (resolve_camera(&settings)) return 2;
|
|
SpacetimeSource spacetime = {0};
|
|
ObserverState observer;
|
|
if (settings.lens_map_input_path == NULL) {
|
|
if (spacetime_create_default(&spacetime)) {
|
|
fputs("Could not create spacetime source\n", stderr);
|
|
return 1;
|
|
}
|
|
if (settings.frames_dir == NULL &&
|
|
build_observer(&settings, &spacetime, &observer)) {
|
|
spacetime_destroy(&spacetime);
|
|
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);
|
|
spacetime_destroy(&spacetime);
|
|
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);
|
|
spacetime_destroy(&spacetime);
|
|
return 1;
|
|
}
|
|
fprintf(stderr, "Created test catalog: %s\n", settings.catalog_path);
|
|
}
|
|
if (blackbody_backend_init(settings.blackbody_table_path,
|
|
0, NAN, NAN, NULL, stderr)) {
|
|
fprintf(stderr,
|
|
"Blackbody backend '%s' initialization failed; the LUT loads the "
|
|
"repository CIE GRBBLUT3 table by default or a valid explicit "
|
|
"table.\n",
|
|
blackbody_backend_name());
|
|
catalog_destroy(&catalog);
|
|
spacetime_destroy(&spacetime);
|
|
return 1;
|
|
}
|
|
fprintf(stderr, "Blackbody backend: %s\n", blackbody_backend_name());
|
|
FastPsfAccumulator fast_accumulator = {0};
|
|
if (settings.fast_mode && settings.lens_map_input_path == NULL) {
|
|
if (fast_psf_accumulator_init(&fast_accumulator, settings.width,
|
|
settings.height, settings.fast_supersample,
|
|
settings.fast_deposit, &settings.psf,
|
|
settings.psf_relative_tail, settings.psf_min_y,
|
|
settings.psf_direct)) {
|
|
fputs("Fast-mode accumulator construction failed.\n", stderr);
|
|
catalog_destroy(&catalog);
|
|
spacetime_destroy(&spacetime);
|
|
blackbody_backend_destroy();
|
|
return 1;
|
|
}
|
|
settings.fast_psf = &fast_accumulator;
|
|
fast_psf_accumulator_report(&fast_accumulator, stderr);
|
|
fast_psf_accumulator_set_verbose(&fast_accumulator, settings.verbose);
|
|
} else if (settings.fast_mode) {
|
|
fputs("Fast mode on an imported lens map uses the map's own dimensions.\n",
|
|
stderr);
|
|
}
|
|
if (!settings.fast_mode && !settings.psf_direct &&
|
|
psf_kernel_cache_init(&settings.psf_cache, &settings.psf,
|
|
settings.psf_relative_tail)) {
|
|
#ifdef PSF_BACKEND_DUMMY
|
|
fputs("PSF cache construction failed; dummy chunk statistics are unavailable.\n",
|
|
stderr);
|
|
catalog_destroy(&catalog);
|
|
spacetime_destroy(&spacetime);
|
|
blackbody_backend_destroy();
|
|
return 1;
|
|
#else
|
|
fputs("PSF cache construction failed; using direct evaluator.\n", stderr);
|
|
#endif
|
|
}
|
|
if (!settings.fast_mode)
|
|
psf_kernel_cache_report_ready(&settings.psf_cache, stderr);
|
|
if (settings.lens_map_input_path != NULL) {
|
|
const int result = render_lens_map(&settings, &catalog);
|
|
catalog_destroy(&catalog);
|
|
psf_kernel_cache_destroy(&settings.psf_cache);
|
|
fast_psf_accumulator_destroy(&fast_accumulator);
|
|
blackbody_backend_destroy();
|
|
return result == 0 ? 0 : 1;
|
|
}
|
|
int result = settings.frames_dir != NULL
|
|
? render_movie(&settings, &catalog, &spacetime)
|
|
: render_observer_frame(&settings, &catalog, &spacetime,
|
|
&observer, settings.output_path);
|
|
spacetime_destroy(&spacetime);
|
|
catalog_destroy(&catalog);
|
|
psf_kernel_cache_destroy(&settings.psf_cache);
|
|
fast_psf_accumulator_destroy(&fast_accumulator);
|
|
blackbody_backend_destroy();
|
|
return result == 0 ? 0 : 1;
|
|
}
|