GR 4D ray tracing — Phase 0 prototype

minkowski_sky is a deliberately small, CPU-only, single-frame Phase 0 benchmark. It renders point sources from a sky catalog through an analytic backend. It is not a sky texture: each source remains a direction, temperature, and amplitude until its sub-pixel Gaussian PSF is splatted.

Build and render the default 1280 x 720 image:

make run

The program first creates assets/sky_grid_5deg.csv when it is missing. The synthetic catalog places stars every 2 degrees on the union of longitude and latitude lines spaced 10 degrees apart; the two poles are stored only once. The eight octants (four 90-degree longitude sectors in each hemisphere) alternate red temperature_K = 3000 and blue temperature_K = 12000. Red stars use amplitude = 1; blue stars use amplitude = 0.00141095580387, which equalizes their CIE/linear-sRGB luminance under the renderer's blackbody integration. Longitude boundaries belong to the sector to their east and the equator to the northern hemisphere, so boundary stars have a deterministic color. The default output is PNG at output/imgs/minkowski_sky.png.

Optional HIP PSF backend

The default PSF_BACKEND=cpu uses the established OpenMP/private-HDR path. Build the direct-atomic HIP PSF backend explicitly with PSF_BACKEND=hip; it requires HIP/ROCm and a usable GPU agent, and writes a separate _hip binary so it cannot overwrite the CPU renderer:

make PSF_BACKEND=hip SPACETIME=minkowski backend
./build/Release/minkowski_sky_hip --catalog assets/sky_grid_5deg.csv \
  --output output/imgs/minkowski_sky_hip.png

The HIP backend accelerates only cache-eligible PSF events. Existing direct fallbacks remain CPU reference evaluations, with an ordered HDR transfer before and after each fallback. A HIP initialization, upload, kernel, or download error terminates the render; it never switches to the CPU backend silently. Use make PSF_BACKEND=hip SPACETIME=minkowski hip-psf-test for the small GPU-vs-CPU cache-HDR regression. HIP renders also report event, batch, H2D, kernel, and HDR-download timing after each frame. To bound diagnostic-event storage, the report times at most the first 64 batches and labels the timed batch count explicitly; the event and total-batch counts are not sampled.

Movie PNG sequence (Phase A)

Movie mode consumes a canonical observer-track CSV rather than a fixed camera. Each row stores coordinate time, proper time, Cartesian position, and the full four-by-four tetrad (21 columns total). Generate the first reproducible Minkowski benchmark—two coordinate seconds at 30 fps, accelerating from rest to about 0.95c—then render its numbered PNG frames:

make clean && make
mkdir -p output/imgs
./build/minkowski_sky --write-minkowski-accel-track output/minkowski_accel_2s.csv \
  --duration 2 --fps 30 --proper-acceleration 1.52
./build/minkowski_sky --observer-track output/minkowski_accel_2s.csv \
  --frames-dir output/imgs --frames-prefix minkowski_accel \
  --start-time 0 --duration 2 --fps 30 --exposure 1e-5

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. 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.

Reuse a completed lens map

Ray tracing and adaptive mesh refinement are independent of catalog lookup, PSF evaluation, exposure, and tone mapping. --lens-map-output FILE writes the finalized local inverse-lens mesh to one versioned .grlens file after all ray-trace/refinement generations finish; the same invocation still renders its ordinary image. The file stores every final vertex's image position, camera direction, infinity direction, frequency shift, terminal status, and the triangle topology. It also stores the frame dimensions, horizontal FOV, and, for a movie, all frame IDs and times.

For example, trace an analytic Schwarzschild frame once and retain the map:

make SPACETIME=schwarzschild backend
mkdir -p output/imgs output/maps
./build/Release/schwarzschild_sky --catalog assets/sky_grid_5deg.csv \
  --width 640 --height 360 --coarse-cell-pixels 8 --fov-deg 60 \
  --refine-max-level 2 --lens-map-output output/maps/schwarzschild.grlens \
  --output output/imgs/schwarzschild_trace.png

Later, render that saved map with a different catalog, PSF, or exposure without constructing an observer, spacetime source, or geodesic rays:

./build/Release/schwarzschild_sky --lens-map-input output/maps/schwarzschild.grlens \
  --catalog assets/sky_grid_5deg.csv --psf-fwhm-pixels 6 --psf-moffat-beta 3 \
  --exposure 0.002 --output output/imgs/schwarzschild_restyled.png

--lens-map-input and --lens-map-output are mutually exclusive. An imported map always retains its original pixel width, height, and horizontal FOV; an explicit conflicting --width, --height, or --fov-deg is rejected. This is intentional: the stored inverse map and PSF coordinates are in the original pixel geometry. The reader validates the format version, finite values, unit directions, triangle indices, and a per-frame CRC before rendering.

A movie export writes every final frame mesh into the same .grlens file. Export it during the usual movie render, then import it with --frames-dir and --frames-prefix; importing a multi-frame map does not require --observer-track:

./build/Release/minkowski_sky --lens-map-input output/maps/movie.grlens \
  --catalog assets/sky_grid_5deg.csv --frames-dir output/imgs \
  --frames-prefix restyled

PNG is the default output and the default build links libpng:

make clean && make
mkdir -p output/imgs
./build/minkowski_sky --output output/imgs/minkowski_sky.png

If libpng is unavailable, rebuild with make clean && make ENABLE_PNG=0. That intentionally selects the binary-PPM fallback, whose default path is output/imgs/minkowski_sky.ppm; pass a .ppm path for explicit output.

Run the flat-spacetime geodesic regression with:

make test

This also checks that the Kerr--Schild metric remains finite at r=2M and that the central ray from the default Schwarzschild camera is classified as captured.

Build an independent analytic Schwarzschild executable in Cartesian ingoing Kerr--Schild coordinates (regular at the horizon), then render the test catalog to PNG:

make clean && make SPACETIME=schwarzschild
mkdir -p output/imgs
./build/schwarzschild_sky --catalog assets/sky_grid_5deg.csv \
  --width 640 --height 360 --coarse-cell-pixels 8 --fov-deg 60 \
  --output output/imgs/schwarzschild_test_catalog.png

SPACETIME=minkowski (the default) and SPACETIME=schwarzschild select source files at compile time, so each executable contains exactly one metric provider. The Schwarzschild demonstration uses mass M=1 and places a static camera at coordinate radius 30 by default. --look-ra-deg and --look-dec-deg define the direction from the camera to the hole; the camera is placed at the opposite direction from the origin and its local forward axis points radially inward. Use --observer-radius R to select any R > 2; it is a coordinate radius in Cartesian Kerr--Schild coordinates. The backend escapes at r=256 and declares capture at r=1.5, safely inside the horizon at r=2. Those rendering thresholds are Phase-1 demonstration values, not settled production refinement or integration settings.

For a local radial boost relative to that static camera, pass --observer-inward-speed V, where 0 <= V < 1 is measured in the static observer's orthonormal frame and positive values point toward the hole. The default is 0, preserving the static camera.

For rays that asymptote to the future horizon in coordinate-time backward integration, the Schwarzschild demo also terminates at log(alpha p^0) = 8. This is the normalized-momentum horizon diagnostic already evolved by the integrator; it is disabled by default and does not replace the AH-calibrated spatial capture criterion planned for nmesh data.

Useful options:

./build/minkowski_sky --width 1920 --height 1080 --fov-deg 30 \
  --catalog assets/sky_grid_5deg.csv --output output/imgs/frame.png
./build/minkowski_sky --catalog assets/2mass/processed/2mass_psc_m31_0p5deg_stars.csv \
  --look-ra-deg 10.6847083 --look-dec-deg 41.26875 --fov-deg 1.8 \
  --exposure 1e15 --output output/imgs/2mass_m31.png
./build/minkowski_sky --catalog assets/2mass/processed/2mass_psc_m44_1p0deg_stars.csv \
  --look-ra-deg 129.99165 --look-dec-deg 19.54139 --fov-deg 2.0 \
  --exposure 1e15 --width 1920 --height 1920 --output output/imgs/2mass_m44.png
./build/minkowski_sky --write-catalog assets/sky_grid_5deg.csv

The camera is a fixed inertial observer at coordinate position (0,0,0), with a tetrad whose forward direction is coordinate -Z and whose vertical direction is +Y. The frame first triangulates the image plane, then traces only its vertices backwards. Escaped endpoints form a triangulation on the source sky. For every locally invertible triangle, catalog stars inside its spherical source triangle are interpolated back to the image triangle and splatted as PSFs. Consequently multiple image triangles naturally create multiple images of the same star.

The ray state evolves (x^i, Pi_i, log(alpha p^0)) in coordinate time with RK4 using the 3+1 equations in Bohn et al. II.A, until the spacetime backend classifies the ray. spacetime.c is the only module containing the Minkowski metric or its infinity criterion; frame, observer, and integrator use only SpacetimeSource and MetricData. The initial regular mesh size is exposed as --coarse-cell-pixels; it is a Phase-0 sampling knob, not a settled production refinement threshold.

Adaptive image mesh refinement

Adaptive refinement is disabled by default (--refine-max-level 0), so the existing coarse-mesh renders remain unchanged. When enabled, its defaults are an absolute direction error of 1e-3 degrees, relative error 0.1, minimum long edge 0.5 pixels, minimum area 0.25 pixel-squared, and a provisional minimum discrete-Jacobian magnitude of 1e-3. Each value can be overridden independently:

--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

N caps the triangle refinement level. Let e be the angle between the traced longest-edge midpoint direction and the normalized endpoint interpolation, and let s be the angle between those two endpoint camera directions. Both are evaluated internally in radians; the absolute CLI threshold D is specified in degrees and converted before comparison. s is the angular geometric size of the image triangle's test edge, not a source-sky/lens-map length. A locally escaped triangle is split only when both e > D_rad (the converted --refine-angle-abs-deg D) and e / max(s, 1e-15) > --refine-angle-rel. P and A prevent selecting a leaf already at or below the requested image-plane long-edge and area scales. Triangles whose three vertices disagree between capture and escape are split independently of the direction-error thresholds, allowing the mesh to follow a shadow boundary.

Independently of the midpoint geometry test, an all-escaped triangle also computes the discrete lens Jacobian J = Omega_source / Omega_image. Both signed solid angles use 2 atan2(dot(a, cross(b,c)), 1 + dot(a,b) + dot(b,c) + dot(c,a)), with the ordered camera directions for Omega_image and their traced infinity directions for Omega_source. A J-driven split requires both a shared image edge whose incident triangles have opposite nonzero signs of J and min(abs(J_left), abs(J_right)) < --refine-jacobian-min. It then requests that shared edge on both leaves. Thus |J| bounds the fold selection instead of widening it as a standalone critical-curve band. A negative sign is physical parity and is retained. The 1e-3 default is deliberately provisional and should be tuned with the small Schwarzschild refinement diagnostic before being treated as a production threshold.

For a short Schwarzschild diagnostic that permits at most one actual split generation, for example:

./build/schwarzschild_sky --catalog assets/sky_grid_5deg.csv \
  --width 48 --height 48 --coarse-cell-pixels 24 --fov-deg 40 \
  --refine-max-level 1 --refine-angle-abs-deg 0.001 \
  --refine-angle-rel 0.001 --refine-jacobian-min 0.001 \
  --refine-min-edge-pixels 1 \
  --refine-min-area-pixels2 1 --draw-mesh \
  --output output/imgs/schwarzschild_refinement.png

For movies, each refinement generation completes the full newest-to-oldest time-slab sweep before any probe becomes a mesh vertex. Newly added vertices are therefore traced only by the next generation; the renderer never returns to a slab that has already been released. At the start of every generation, newly inserted vertices and geometry-only longest-edge probes for its new leaves are collected together, so both ray sets use the same parallel RayPool pass.

Catalog directions and --look-ra-deg/--look-dec-deg use standard right-handed ICRS Cartesian axes: +X is RA 0 degrees/Dec 0 degrees, +Y is RA 90 degrees/Dec 0 degrees, and +Z is the north celestial pole. The local camera axes are forward, celestial north, and celestial west, so an image with north up has decreasing RA to the right. The defaults preserve the original -Z view. --exposure converts a catalog's physical flux normalization to the prototype HDR scale. The current synthetic catalog is calibrated for default exposure 1e-3; a 2MASS blackbody normalization in steradians requires a much larger display exposure such as the example above. The optics path integrates each fitted Planck spectrum through CIE 1931 color-matching functions and converts the resulting radiance to linear sRGB; it does not use an empirical color-temperature RGB approximation.

Point sources use a flux-normalized circular Moffat PSF by default (--psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5). The FWHM matches the former 1.15-pixel Gaussian core while the Moffat wings remain continuous; both values are display/optics calibration parameters.

The default renderer builds one immutable, process-wide 64-by-64 sub-pixel Moffat lookup kernel. Its weights are pixel-area integrals and are bilinearly interpolated between phase tables. The renderer reports its build time and cached/direct-fallback image counts. Pass --psf-direct to use the slower 8-point quadrature reference evaluator for regression comparisons; an image whose required HDR-tail support exceeds the cache radius selects that reference path automatically.

--psf-relative-tail R controls the maximum omitted PSF tail fraction per image; it defaults to 1e-8. Larger values intentionally shorten the Moffat support and rebuild the immutable cache at the corresponding radius, which is useful when a faster, lower-fidelity render is acceptable.

--psf-min-y Y defaults to 0 (disabled). A positive value is a linear-HDR luminance cutoff: a PSF stops where its continuous Moffat Y profile falls below Y, and an image whose central value is already below Y is omitted. The final PSF report counts such omitted events and emits a warning when any occur.

For bounded preview renders, --max-cache-psf-flux F (default 1) allows images with 1 < flux <= F to use the existing cache instead of the direct evaluator. This does not rebuild or enlarge the cache: the image retains its true core flux, color, and sub-pixel position, while its Moffat wing is clipped at the existing cache radius. Images above F retain the direct fallback. --max-magnification M (default unlimited) caps the per-triangle rendering magnification before flux is formed; it is an explicit preview approximation.

The PSF-cache completion line is printed before tracing and catalog splatting begin. For long renders, pass --verbose to print catalog-prefetch state, splat-worker local heartbeats (8, 16, 32, ... completed triangles per worker), and image-write boundaries. The worker heartbeats use neither global progress accounting nor cross-worker synchronization. Verbose ray-trace output reports the initial mesh trace and refinement stages for single frames; movie mode additionally reports each generation's sample count and each time slab's activation and terminal-ray summary. Movie renders always print one summary per time slab; --verbose also prints the ray counts before each slab is loaded.

To preserve a render for later exposure and tone-mapping work, build the desired backend with ENABLE_HDR=1, for example make SPACETIME=schwarzschild ENABLE_HDR=1. This produces build/Release/schwarzschild_sky, which accepts --hdr-output. This switch writes the HDR file next to the ordinary output, replacing its extension with _HDR.fits; for example, --output output/imgs/ring.png --hdr-output writes output/imgs/ring_HDR.fits. It writes the pre-tone-mapping RGB framebuffer as a three-plane, 32-bit float FITS image. Values remain linear HDR at the renderer's arbitrary scale; no tone mapping or per-frame normalization is applied. The ordinary binaries do not contain this option or writer.

The FITS header describes a synthetic 8640-by-5760, 36-by-24 mm full-frame sensor with 4.1667 um pixels. Each render records its active centered crop and derives FOCALLEN from that render's width and horizontal --fov-deg; these camera fields support plate-solving workflows but do not calibrate flux.

Pass --draw-mesh to alpha-composite image-plane triangle edges as one-pixel-wide 0.5 linear-gray diagnostic lines at 0.5 opacity. The line rasterizer uses coverage-based antialiasing.

S
Description
No description provided
Readme
73 MiB
0 Stars 1 Watchers 0 Forks
Languages
C 74%
Python 15.7%
HIP 5.2%
Assembly 3.6%
Makefile 1.1%
Other 0.4%