wyj/GR-raytracing
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Doc: reorganize bilingual READMEs and build guide with reference renders

Browse Source
This commit is contained in:
wyj committed 2026-09-06 02:03:50 -04:00
1 parent a18ebfbf1a
commit 94149d75e4
8 files changed
+612 -332
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

README.md
+122 -332
View File
@@ -1,360 +1,150 @@
# GR 4D ray tracing — Phase 0 prototype
# GR 4D ray tracing
`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.
**English** | [简体中文](README.zh-CN.md)
Build and render the default 1280 x 720 image:
An offline spacetime renderer focused on physical accuracy. The project aims
to turn time-dependent numerical-relativity simulations, including binary
black-hole mergers, into 4K movies by tracing light rays backward through the
four-dimensional spacetime.
[![Schwarzschild black-hole lensing of the 2MASS Galactic-center star field](assets/images/schwarzschild_galactic_center.png)](assets/images/schwarzschild_galactic_center.png)
*The 2MASS Galactic-center star field seen through Schwarzschild spacetime,
from a camera at radius 100 M. Click the image for the full 4K render;
[rendering command below](#example-galactic-center-field-through-schwarzschild-spacetime).*
Stars are individual catalog point sources with direction, temperature, and
amplitude. The renderer maps them into the camera image, accounting for
multiple images, gravitational lensing magnification, and frequency shifts,
then accumulates their sub-pixel point-spread functions (PSFs) into an HDR
image. Moving cameras are described by worldline and tetrad tracks.
## Current status
The current implementation supports analytic **Minkowski** and
**Schwarzschild** spacetimes, single images and observer-track image sequences,
adaptive lens meshes, and reusable lens-map files. It is written primarily in
C with OpenMP CPU parallelism; an optional HIP backend accelerates PSF
accumulation.
The [Nmesh](https://github.com/nmeshsource/nmesh) numerical-spacetime backend and BBH rendering are still planned.
The current scope is black-hole capture and distant stellar backgrounds;
local matter emission, accretion disks, and plasma are outside this stage.
See the [design document](nr_spacetime_movie_renderer_design.md) for the
architecture and development roadmap.
## Build
The default build requires a C11 compiler with OpenMP support, GNU Make, and
libpng development files. From the repository root:
```sh
make run
make -j
```
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`.
This builds both `build/Release/minkowski_sky` and
`build/Release/schwarzschild_sky`. For individual backends, Debug builds,
optional HDR/FITS output, HIP support, and regression checks, see
[build.md](build.md).
### Optional HIP PSF backend
## Prepare the stellar catalog
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:
The renderer accepts **any stellar catalog converted to the supported CSV
format**: four columns containing ICRS right ascension and declination in
degrees, temperature in kelvin, and amplitude, in that order. Use
`--catalog PATH` to load a single CSV.
**2MASS is the recommended survey catalog**, with download and processing
scripts provided by this project. Instructions are in
[assets/2mass/README.md](assets/2mass/README.md).
The full download takes approximately **three to four days**, so **contact me
for a compressed archive first** if possible. Place the processed tile files
under `assets/2mass/processed/all_sky/` for the commands below.
The included [assets/sky_grid_5deg.csv](assets/sky_grid_5deg.csv) is a synthetic
stellar grid for geometry and regression tests. It can be used directly with
`--catalog assets/sky_grid_5deg.csv`, without downloading survey data; see the
[test-grid example below](#example-synthetic-test-grid-with-mesh-overlay).
## Rendering examples
The examples below illustrate a few choices of sky field, spacetime, and
camera settings. After building and preparing the catalog, adapt these
commands to your own field of view, observer position, and rendering settings.
See [usage.md](usage.md) for the available controls and workflows.
### Example: Summer Triangle in flat spacetime
This example renders a field around the Summer Triangle:
```sh
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:
```sh
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
./build/Release/minkowski_sky \
--all-sky-catalog assets/2mass/processed/all_sky \
--look-ra-deg 296 --look-dec-deg 27 --fov-deg 72 \
--width 3840 --height 2160 --exposure 1e12 \
--output output/imgs/summer_triangle.png
```
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.
[![Summer Triangle in flat spacetime, exposure 1e12](assets/images/summer_triangle.png)](assets/images/summer_triangle.png)
## Reuse a completed lens map
*4K reference image with exposure `1e12`. Click to view at full resolution.*
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.
Angles are in degrees; `--fov-deg` is the horizontal field of view. Exposure
is an adjustable display multiplier. Use `--verbose` for progress during long
renders.
For example, trace an analytic Schwarzschild frame once and retain the map:
### Example: Galactic-center field through Schwarzschild spacetime
This example uses a camera at radius 100 M to render the lensed
Galactic-center field shown at the top of this page:
```sh
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:
```sh
./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`:
```sh
./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`:
```sh
make clean && make
mkdir -p output/imgs
./build/minkowski_sky --output output/imgs/minkowski_sky.png
./build/Release/schwarzschild_sky \
--all-sky-catalog assets/2mass/processed/all_sky \
--width 3840 --height 2160 \
--look-ra-deg 262.5 --look-dec-deg -30 --fov-deg 45 \
--observer-radius 100 --exposure 1e13 \
--coarse-cell-pixels 16 --refine-max-level 4 --refine-jacobian-min 0.2 \
--psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 \
--psf-relative-tail 1e-8 --psf-min-y 0 --max-cache-psf-flux 1e8 \
--catalog-load-workers 4 \
--output output/imgs/schwarzschild_galactic_center.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.
This uses the reference image's rendering settings, including its
`--max-cache-psf-flux 1e8` preview approximation, which clips bright PSF wings
to the cache radius. The [original benchmark record](benchmarks/2mass_galactic_center_blackhole.md)
preserves the command and terminal output; the command above omits the optional
HDR export.
Run the flat-spacetime geodesic regression with:
The Schwarzschild camera points toward the hole. Adaptive refinement should
be configured for the desired image accuracy; it is disabled by default.
Camera controls, movie sequences, lens-map reuse, PSF settings, and HDR output
are described in [usage.md](usage.md). Both binaries provide a complete option
list with `--help`.
### Example: Synthetic test grid with mesh overlay
This example uses `assets/sky_grid_5deg.csv` to inspect lensing and adaptive
mesh refinement in Schwarzschild spacetime. `--draw-mesh` overlays the final
image-plane triangles.
```sh
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:
```sh
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
./build/Release/schwarzschild_sky \
--max-cache-psf-flux 1e8 \
--catalog assets/sky_grid_5deg.csv \
--refine-max-level 3 --refine-jacobian-min 0.2 \
--width 3840 --height 2160 \
--look-ra-deg 0.1 --look-dec-deg 0.1 --fov-deg 45 \
--coarse-cell-pixels 32 \
--observer-radius 100 --exposure 0.2 \
--psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --draw-mesh \
--output output/imgs/schwarzschild_test_grid.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.
[![Synthetic stellar grid lensed by a Schwarzschild black hole, with adaptive mesh overlay](assets/images/schwarzschild_test_grid.png)](assets/images/schwarzschild_test_grid.png)
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:
```sh
./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:
```text
--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:
```sh
./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.
*4K test-grid reference image. Click to view at full resolution.*