From 80b891d2c8709a242b4fb4c68955195c16dad320 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Fri, 28 Aug 2026 13:53:06 -0400 Subject: [PATCH] Optics: cache pixel-area Moffat kernels --- Makefile | 9 +- README.md | 13 + .../2mass_all_sky_psf_lookup_2026-08-28.md | 65 ++++ psf_lookup_optimization_plan.md | 129 +++++++ src/frame.c | 73 ++-- src/frame.h | 4 +- src/main.c | 41 ++- src/optics.c | 320 ++++++++++++++++-- src/optics.h | 37 ++ tests/test_frame.c | 58 +++- 10 files changed, 693 insertions(+), 56 deletions(-) create mode 100644 benchmarks/2mass_all_sky_psf_lookup_2026-08-28.md create mode 100644 psf_lookup_optimization_plan.md diff --git a/Makefile b/Makefile index 9408958..5f52ada 100644 --- a/Makefile +++ b/Makefile @@ -22,8 +22,9 @@ FRAME_TEST_TARGET := build/test_frame SCHWARZSCHILD_TEST_TARGET := build/test_schwarzschild OBSERVER_TRACK_TEST_TARGET := build/test_observer_track CATALOG_PREFETCH_TEST_TARGET := build/test_catalog_prefetch +PSF_HDR_TEST_TARGET := build/minkowski_psf_hdr_test -.PHONY: all clean run test minkowski schwarzschild +.PHONY: all clean run test minkowski schwarzschild psf-hdr-test ifeq ($(SPACETIME),minkowski) BACKEND_CPPFLAGS := -DSPACETIME_MINKOWSKI @@ -51,6 +52,12 @@ minkowski: schwarzschild: $(MAKE) SPACETIME=schwarzschild all +# Deliberately separate from the production binary: enables --hdr-output PFM. +$(PSF_HDR_TEST_TARGET): $(COMMON_SOURCES) src/spacetime_minkowski.c src/main.c | build + $(CC) $(CPPFLAGS) -DENABLE_HDR_DEBUG -DSPACETIME_MINKOWSKI $(CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ + +psf-hdr-test: $(PSF_HDR_TEST_TARGET) + $(TEST_TARGET): tests/test_geodesic.c $(CORE_MINKOWSKI_SOURCES) | build $(CC) $(CPPFLAGS) $(CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ diff --git a/README.md b/README.md index 760ad84..8d4d6ab 100644 --- a/README.md +++ b/README.md @@ -158,6 +158,19 @@ Point sources use a flux-normalized circular Moffat PSF by default 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. + +For PSF validation only, `make psf-hdr-test` builds +`build/minkowski_psf_hdr_test`, a separate binary with a `--hdr-output PATH` +option. It writes the pre-tone-mapping RGB framebuffer as a 32-bit float PFM +image; the ordinary binaries do not contain this option or writer. + 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. diff --git a/benchmarks/2mass_all_sky_psf_lookup_2026-08-28.md b/benchmarks/2mass_all_sky_psf_lookup_2026-08-28.md new file mode 100644 index 0000000..9362cba --- /dev/null +++ b/benchmarks/2mass_all_sky_psf_lookup_2026-08-28.md @@ -0,0 +1,65 @@ +# 2MASS 全天 1080p Moffat lookup 基准（2026-08-28） + +## 目的 + +测量 `PsfKernelCache` 在真实全天 2MASS catalog、高星像数单帧中的端到端 +收益，并与查表实现之前的直接 Moffat 路径比较。两个运行使用相同的相机、 +catalog、输出尺寸与曝光；唯一意图中的算法差异是 PSF 计算。 + +## 可复现环境 + +- Git base revision: `a71180706f3abfb9fe7669a2366944ea3141d78c`。查表结果来自 + 其上的 PSF lookup 工作树，随后由本次提交固化；基线为改动前二进制。 +- Host CPU: 12th Gen Intel Core i7-12700K，16 个 online logical CPUs。 +- Build defaults: C11、`-march=native -O2 -pipe`、OpenMP、PNG output。 +- Catalog loading: `--catalog-load-workers 4`（默认值）。 +- Render workers: 1920 x 1080 的 RGB `double` private HDR buffer 约为 + 47.5 MiB；在 `FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES = 512 MiB` 限制下， + 最多使用 **10 个** private-HDR render workers。 +- Catalog: `assets/2mass/processed/all_sky`；本视场请求并加载 9,328 个 + tiles，共 18,148,775 颗 catalog stars。 + +## 控制变量与命令 + +相机朝向南天极 `(RA, Dec) = (0 deg, -90 deg)`，60-degree 水平视场， +1920 x 1080，`--coarse-cell-pixels 16`，exposure `1e12`。两次均产生 +16,555,241 个星像。 + +```sh +time ./build/minkowski_sky \ + --all-sky-catalog assets/2mass/processed/all_sky \ + --look-ra-deg 0 --look-dec-deg -90 \ + --width 1920 --height 1080 \ + --coarse-cell-pixels 16 --fov-deg 60 \ + --exposure 1e12 \ + --output output/imgs/2mass_1080p_60deg_south_pole_1e12_psf_cache.png +``` + +用 `--psf-direct` 可在同一查表版本中运行 8-point 直接积分参考；历史基线 +在该选项加入前使用逐像素 `pow()` 直接 Moffat 求值。 + +## 结果 + +| 路径 | images | cache build | cached / direct fallback | prefetch | user | system | wall | +|---|---:|---:|---:|---:|---:|---:|---:| +| 历史直接 Moffat | 16,555,241 | — | — | 3.913 s | 1617.84 s | 2.33 s | 165.59 s | +| 64x64 lookup | 16,555,241 | 6.520 s | 16,555,241 / 0 | 3.981 s | 645.24 s | 4.02 s | 74.48 s | + +查表版本端到端 wall time 加速为 `165.59 / 74.48 = 2.22x`，节省 91.11 s +（55.0%）；user CPU time 加速为 `2.51x`。预取差异仅 0.068 s，不能解释 +该提升。扣除一次性 kernel build 后，查表运行剩余 67.96 s；该数仅用于说明 +多帧或更多星像时的摊销潜力，不应与没有建表步骤的历史单帧总时间混为同一 +端到端指标。 + +这次 lookup 运行的 `cached splats` 等于 image count，且 direct fallback 为 +零，因此该速度结果没有把亮星回退的成本隐藏在统计之外。 + +## 输出一致性与范围 + +本基准确认 catalog 预取、source-triangle 查询和星像计数在优化前后不变。 +PSF 的 float-HDR 像素精度验收另见单星直接积分比较：默认 `FWHM=2.7, +beta=4.5` 的最大误差为直接参考峰值的 `3.18e-5`，RMS 为 `8.25e-7`。 + +本记录比较的是历史点采样直接路径与新的像素面积积分 lookup 路径的性能。 +在需要逐像素数值回归时，使用新二进制的 `--psf-direct`（8-point +pixel-area reference）而非历史路径。 diff --git a/psf_lookup_optimization_plan.md b/psf_lookup_optimization_plan.md new file mode 100644 index 0000000..978a58f --- /dev/null +++ b/psf_lookup_optimization_plan.md @@ -0,0 +1,129 @@ +# CPU PSF lookup-table optimization plan + +## Goal + +Speed up the current CPU point-source rendering path when the physically +correct lens map produces very many star images. This plan changes only how a +known image's Moffat PSF is evaluated and accumulated. It does **not** change +the catalog representation, inverse lens map, number of images, frequency +shift, magnification, adaptive-mesh work, or use of the GPU. + +The target is to replace the per-covered-pixel Moffat `pow()` evaluation with +a lookup of the same continuously translated PSF. It must preserve continuous +sub-pixel image positions and have an explicit, testable bound on PSF-tail +error, including for bright sources. + +## Interface and ownership + +Add a `PsfKernelCache` owned by the rendering settings for the lifetime of one +program invocation. Exactly one cache is built after command-line parsing, +using the process-wide values of: + +- `--psf-fwhm-pixels`; +- `--psf-moffat-beta`; +- a named PSF-tail/display-error configuration; +- a fixed sub-pixel phase resolution. + +It is immutable after construction, shared read-only by splat workers, and is +destroyed after rendering. It is not serialized to disk and no multi-parameter +cache or LRU is required: the current program has one PSF parameter pair for +the entire image or movie. + +`splat_moffat()` gains a cache-aware path (or is replaced by an equivalently +named cache-aware function). A direct-evaluation implementation remains +available as a reference and as a fallback if cache construction fails. + +## Kernel construction + +For each sub-pixel phase `(fx, fy)`, where + +`fx = x - floor(x)` and `fy = y - floor(y)`, build a scalar kernel containing +the **pixel-area integral** of the normalized circular Moffat profile over each +integer pixel footprint. Kernel weights are independent of color and flux. +Use a fixed, documented numerical quadrature whose convergence is verified +against a higher-accuracy reference. + +At runtime, retain the image event's full-precision `(x, y)`. Select the four +neighboring phase tables and bilinearly interpolate their weights. Do not +round to a nearest phase table. Normalize every interior phase kernel so its +weight sum is one; an image clipped by the framebuffer boundary retains the +current behavior of losing the off-frame contribution rather than being +renormalized. + +Start with a 64 by 64 phase grid. This is an initial implementation value, +not an accuracy claim: increase it if the phase-continuity and reference-image +tests below do not meet their thresholds. + +## Tail support and bright sources + +A Moffat profile has infinite support, so a finite renderer needs an explicit +finite-support approximation. Do not rely solely on the current fixed +relative omitted-flux fraction for every source. For each image event, choose +the support radius from its final post-lens, post-frequency-shift, post-exposure +flux so that both of these hold: + +- omitted integrated flux is below a configured absolute frame error budget; +- the PSF value at the support boundary is below a configured per-pixel HDR + error budget. + +The cache contains weights through the largest radius required by the selected +configuration. A dim image uses a cropped subset; an unusually bright image +uses a larger subset. If an event requires a radius beyond the cache maximum, +use the direct Moffat evaluator for that event and record it in diagnostics. +This prevents a visible hard cutoff in bright-star wings without silently +discarding their flux. + +The selected tail budgets, maximum cache radius, and fallback count must be +reported with the render diagnostics. Their final numerical values are to be +chosen from the validation experiments, not inferred from a performance goal. + +## Accumulation and parallel behavior + +For every image event, compute color and final flux exactly as in the existing +path. Then multiply the scalar cached PSF weights by that color and flux and +add them to the worker's existing private HDR buffer. Preserve the current +private-HDR ownership, memory cap, reduction ordering, and serial fallback; +this plan does not introduce atomics, locks, or a different parallel +partitioning. + +## Validation and benchmark + +Add tests and a reproducible benchmark before making the lookup path default: + +1. Compare cached and high-accuracy direct pixel-area integration for a grid + of FWHM, beta, phase, flux, and support-radius cases. Check kernel sum, + per-pixel error, and omitted-flux bound. +2. Translate a star in sub-pixel increments across pixel boundaries. Verify + continuous image values and absence of phase quantization jumps. +3. Test bright stars whose required support is larger than the ordinary cache + crop, including the direct-fallback case. Verify the configured tail-error + bound and no visible-radius discontinuity relative to the reference. +4. Compare serial and OpenMP outputs under the existing private-HDR contract. + Preserve the current expected identity/tolerance policy explicitly. +5. Benchmark the Schwarzschild 1080p, 60-degree, FWHM 1.0 all-sky case and + separately report cache-build time, number of image events, cached-splat + time, direct-fallback count, and total render time. Compare output against + the direct evaluator with the same error budgets. + +The lookup path becomes the default only after it meets the agreed visual and +numerical tolerances and demonstrates a material speedup in the PSF-dominated +case. The direct evaluator remains a selectable regression reference. + +## Implementation record (2026-08-28) + +Implemented as one immutable `PsfKernelCache` per invocation, built after CLI +parsing and passed read-only through the existing frame splat workers. The +cache uses 64 by 64 phase intervals (65 nodes on each axis), 4-point +Gauss-Legendre pixel-area quadrature, and bilinear phase interpolation. The +`--psf-direct` regression mode uses an independent 8-point quadrature +reference; events whose requested support exceeds the cache use that path +automatically. + +The selected HDR budgets are omitted flux `1e-6` and boundary contribution +`1e-7`; the former `1e-8` relative tail remains part of the support selection. +The cache reports build time, maximum radius, cached splat count, and direct +fallback count. `tests/test_frame.c` covers cache construction, phases near +pixel boundaries, comparison with the direct reference, and bright-event +fallback. The float-HDR default-PSF acceptance result and the 1080p full-sky +performance result are recorded in +`benchmarks/2mass_all_sky_psf_lookup_2026-08-28.md`. diff --git a/src/frame.c b/src/frame.c index d5611de..2f1b2dc 100644 --- a/src/frame.c +++ b/src/frame.c @@ -184,7 +184,9 @@ typedef struct { int width, height; double exposure, magnification; const PointSpreadFunction *psf; + const PsfKernelCache *psf_cache; size_t images; + size_t direct_fallbacks; } TriangleSplatContext; static int splat_catalog_tile(const Star *stars, size_t count, @@ -218,9 +220,10 @@ static int splat_catalog_tile(const Star *stars, size_t count, weights[1] * context->vertex[1]->log_frequency_ratio + weights[2] * context->vertex[2]->log_frequency_ratio; const LinearRgb color = blackbody_to_linear_rgb(star->temperature_K * exp(log_g)); - splat_moffat(context->hdr, context->width, context->height, image_x, image_y, - color, context->exposure * star->amplitude * context->magnification, - context->psf); + context->direct_fallbacks += splat_moffat_cached( + context->hdr, context->width, context->height, image_x, image_y, color, + context->exposure * star->amplitude * context->magnification, + context->psf, context->psf_cache); ++context->images; } return 0; @@ -230,8 +233,10 @@ static size_t splat_catalog_triangles(const FrameLensMesh *mesh, StarCatalog *catalog, double *hdr, int width, int height, double exposure, const PointSpreadFunction *psf, + const PsfKernelCache *psf_cache, size_t first_triangle, - size_t last_triangle) { + size_t last_triangle, + size_t *direct_fallbacks) { size_t images = 0; for (size_t t = first_triangle; t < last_triangle; ++t) { const LensVertex *vertex[3]; @@ -251,10 +256,11 @@ static size_t splat_catalog_triangles(const FrameLensMesh *mesh, .triangle = &mesh->triangles[t], .hdr = hdr, .width = width, .height = height, .exposure = exposure, .magnification = magnification, - .psf = psf}; + .psf = psf, .psf_cache = psf_cache}; if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile, &context) == 0) images += context.images; + *direct_fallbacks += context.direct_fallbacks; } return images; } @@ -283,8 +289,10 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, StarCatalog *catalog, double *hdr, int width, int height, double exposure, const PointSpreadFunction *psf, + const PsfKernelCache *psf_cache, int catalog_load_workers, - CatalogPrefetchStats *prefetch_stats) { + CatalogPrefetchStats *prefetch_stats, + PsfSplatStats *psf_stats) { if (mesh == NULL || catalog == NULL || hdr == NULL || exposure <= 0.0 || psf == NULL || width <= 0 || height <= 0 || catalog_load_workers <= 0) return 0; @@ -292,25 +300,45 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, /* A bounded parallel read phase completes before splatting. Its serial cache * commit leaves immutable tile data for the OpenMP splat workers. */ prefetch_catalog_for_mesh(mesh, catalog, catalog_load_workers, prefetch_stats); + if (psf_stats != NULL) + *psf_stats = (PsfSplatStats){0}; const size_t pixel_count = (size_t)width * height * 3; if (pixel_count > SIZE_MAX / sizeof(double) || pixel_count * sizeof(double) > FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES / 2) - return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure, - psf, 0, mesh->triangle_count); + { + size_t fallbacks = 0; + const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height, + exposure, psf, psf_cache, 0, + mesh->triangle_count, &fallbacks); + if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks}; + return images; + } const size_t buffer_bytes = pixel_count * sizeof(double); size_t worker_count = FRAME_SPLAT_MAX_PRIVATE_HDR_BYTES / buffer_bytes; const int max_threads = omp_get_max_threads(); if (worker_count > (size_t)max_threads) worker_count = (size_t)max_threads; if (worker_count < 2 || worker_count > INT_MAX) - return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure, - psf, 0, mesh->triangle_count); + { + size_t fallbacks = 0; + const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height, + exposure, psf, psf_cache, 0, + mesh->triangle_count, &fallbacks); + if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks}; + return images; + } double **private_hdr = calloc(worker_count, sizeof *private_hdr); if (private_hdr == NULL) - return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure, - psf, 0, mesh->triangle_count); + { + size_t fallbacks = 0; + const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height, + exposure, psf, psf_cache, 0, + mesh->triangle_count, &fallbacks); + if (psf_stats != NULL) *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks}; + return images; + } size_t allocated = 0; for (; allocated < worker_count; ++allocated) { private_hdr[allocated] = calloc(pixel_count, sizeof **private_hdr); @@ -321,12 +349,17 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, while (allocated > 0) free(private_hdr[--allocated]); free(private_hdr); - return splat_catalog_triangles(mesh, catalog, hdr, width, height, exposure, - psf, 0, mesh->triangle_count); + size_t fallbacks = 0; + const size_t images = splat_catalog_triangles(mesh, catalog, hdr, width, height, + exposure, psf, psf_cache, 0, + mesh->triangle_count, &fallbacks); + if (psf_stats != NULL) + *psf_stats = (PsfSplatStats){images - fallbacks, fallbacks}; + return images; } - size_t images = 0; -#pragma omp parallel num_threads((int)worker_count) reduction(+ : images) + size_t images = 0, direct_fallbacks = 0; +#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks) { const size_t worker = (size_t)omp_get_thread_num(); /* Source density and lens magnification can vary by orders of magnitude @@ -335,9 +368,9 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, * idle. Each worker still owns its HDR buffer exclusively. */ #pragma omp for schedule(dynamic, 1) for (size_t triangle = 0; triangle < mesh->triangle_count; ++triangle) - images += splat_catalog_triangles(mesh, catalog, private_hdr[worker], - width, height, exposure, psf, triangle, - triangle + 1); + images += splat_catalog_triangles(mesh, catalog, private_hdr[worker], width, + height, exposure, psf, psf_cache, triangle, + triangle + 1, &direct_fallbacks); } #pragma omp parallel for schedule(static) for (size_t pixel = 0; pixel < pixel_count; ++pixel) @@ -346,6 +379,8 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, for (size_t worker = 0; worker < worker_count; ++worker) free(private_hdr[worker]); free(private_hdr); + if (psf_stats != NULL) + *psf_stats = (PsfSplatStats){images - direct_fallbacks, direct_fallbacks}; return images; } diff --git a/src/frame.h b/src/frame.h index e1f5f7d..13f6e9a 100644 --- a/src/frame.h +++ b/src/frame.h @@ -39,8 +39,10 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, StarCatalog *catalog, double *hdr, int width, int height, double exposure, const PointSpreadFunction *psf, + const PsfKernelCache *psf_cache, int catalog_load_workers, - CatalogPrefetchStats *prefetch_stats); + CatalogPrefetchStats *prefetch_stats, + PsfSplatStats *psf_stats); void frame_draw_mesh(const FrameLensMesh *mesh, double *hdr, int width, int height, double gray, double opacity); void frame_lens_mesh_destroy(FrameLensMesh *mesh); diff --git a/src/main.c b/src/main.c index d35b261..7223a5c 100644 --- a/src/main.c +++ b/src/main.c @@ -18,13 +18,18 @@ typedef struct { int width, height; int coarse_cell_pixels; int draw_mesh; + int psf_direct; double horizontal_fov_deg, look_ra_deg, look_dec_deg, exposure; double observer_radius; double observer_inward_speed; PointSpreadFunction psf; + PsfKernelCache psf_cache; const char *catalog_path; const char *all_sky_catalog_path; const char *output_path; +#ifdef ENABLE_HDR_DEBUG + const char *hdr_output_path; +#endif const char *observer_track_path; const char *frames_dir; const char *frames_prefix; @@ -126,6 +131,10 @@ static int parse_args(int argc, char **argv, Settings *s, s->all_sky_catalog_path = argv[++i]; else if (!strcmp(argv[i], "--output") && i + 1 < argc) s->output_path = argv[++i]; +#ifdef ENABLE_HDR_DEBUG + else if (!strcmp(argv[i], "--hdr-output") && i + 1 < argc) + s->hdr_output_path = argv[++i]; +#endif else if (!strcmp(argv[i], "--width") && i + 1 < argc && !parse_int(argv[++i], &s->width)) { } else if (!strcmp(argv[i], "--height") && i + 1 < argc && @@ -150,6 +159,8 @@ static int parse_args(int argc, char **argv, Settings *s, !parse_positive(argv[++i], &s->psf.fwhm_pixels)) { } else if (!strcmp(argv[i], "--psf-moffat-beta") && i + 1 < argc && !parse_moffat_beta(argv[++i], &s->psf.moffat_beta)) { + } else if (!strcmp(argv[i], "--psf-direct")) { + s->psf_direct = 1; } else if (!strcmp(argv[i], "--write-catalog") && i + 1 < argc) *write_path = argv[++i]; else if (!strcmp(argv[i], "--observer-track") && i + 1 < argc) @@ -219,15 +230,26 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog, return -1; } CatalogPrefetchStats prefetch = {0}; + PsfSplatStats psf_stats = {0}; size_t images = frame_splat_catalog( &mesh, catalog, hdr, s->width, s->height, s->exposure, &s->psf, - s->catalog_load_workers, &prefetch); + &s->psf_cache, s->catalog_load_workers, &prefetch, &psf_stats); if (s->draw_mesh) frame_draw_mesh(&mesh, hdr, s->width, s->height, 0.5, 0.5); +#ifdef ENABLE_HDR_DEBUG + if (s->hdr_output_path != NULL && + write_hdr_pfm(s->hdr_output_path, hdr, s->width, s->height)) { + perror(s->hdr_output_path); + frame_lens_mesh_destroy(&mesh); + free(hdr); + return -1; + } +#endif int result = write_tonemapped_image(output_path, hdr, s->width, s->height); fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n", images, catalog->count, output_path, result == 0 ? "ok" : "write failed"); + psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr); if (catalog->kind == STAR_CATALOG_ALL_SKY) fprintf(stderr, "Catalog prefetch: %zu requested, %zu newly loaded (%zu stars), " @@ -319,15 +341,17 @@ static int render_movie(const Settings *s, StarCatalog *catalog, goto done; } CatalogPrefetchStats prefetch = {0}; + PsfSplatStats psf_stats = {0}; const size_t images = frame_splat_catalog( &movie.frames[i].mesh, catalog, hdr, s->width, s->height, s->exposure, - &s->psf, s->catalog_load_workers, &prefetch); + &s->psf, &s->psf_cache, s->catalog_load_workers, &prefetch, &psf_stats); if (s->draw_mesh) frame_draw_mesh(&movie.frames[i].mesh, hdr, s->width, s->height, 0.5, 0.5); const int write_result = write_tonemapped_image(output_path, hdr, s->width, s->height); free(hdr); fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n", images, catalog->count, output_path, write_result == 0 ? "ok" : "write failed"); + psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr); if (catalog->kind == STAR_CATALOG_ALL_SKY) fprintf(stderr, "Catalog prefetch: %zu requested, %zu newly loaded (%zu stars), " @@ -369,6 +393,10 @@ int main(int argc, char **argv) { "N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] " "[--exposure E] [--observer-radius R] [--observer-inward-speed V] " "[--psf-fwhm-pixels N] [--psf-moffat-beta N] " + "[--psf-direct] " +#ifdef ENABLE_HDR_DEBUG + "[--hdr-output PATH] " +#endif "[--coarse-cell-pixels N] [--draw-mesh] [--write-catalog PATH] " "[--catalog-load-workers N] " "[--observer-track PATH --frames-dir DIR --frames-prefix NAME " @@ -384,6 +412,12 @@ int main(int argc, char **argv) { return write_minkowski_accel_track(&settings) == 0 ? 0 : (perror(settings.write_minkowski_accel_track_path), 1); +#ifdef ENABLE_HDR_DEBUG + if (settings.frames_dir != NULL && settings.hdr_output_path != NULL) { + fputs("--hdr-output is available only for a single-frame render.\n", stderr); + return 2; + } +#endif StarCatalog catalog = {0}; if (settings.all_sky_catalog_path != NULL) { if (catalog_load_all_sky(&catalog, settings.all_sky_catalog_path)) { @@ -404,10 +438,13 @@ int main(int argc, char **argv) { catalog_destroy(&catalog); return 1; } + if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf)) + fputs("PSF cache construction failed; using direct evaluator.\n", stderr); int result = settings.frames_dir != NULL ? render_movie(&settings, &catalog, &spacetime) : render_frame(&settings, &catalog, &spacetime); spacetime_destroy(&spacetime); catalog_destroy(&catalog); + psf_kernel_cache_destroy(&settings.psf_cache); return result == 0 ? 0 : 1; } diff --git a/src/optics.c b/src/optics.c index ad59ebc..77bb4c9 100644 --- a/src/optics.c +++ b/src/optics.c @@ -1,9 +1,11 @@ #include "optics.h" #include +#include #include #include #include +#include #ifdef ENABLE_PNG #include @@ -76,39 +78,272 @@ LinearRgb blackbody_to_linear_rgb(double temperature_K) 0.0, INFINITY)}; } +enum { + PSF_PHASE_RESOLUTION = 64, + PSF_MAX_CACHE_RADIUS_PIXELS = 48, + PSF_QUADRATURE_ORDER = 4, +}; + +/* These are display-space HDR error budgets, deliberately independent of a + * relative flux cut. The historical 1e-8 relative tail remains a floor for + * ordinary images; brighter images grow their support or use the reference + * fallback rather than acquiring a visible clipped wing. */ +static const double psf_relative_tail_fraction = 1e-8; +static const double psf_tail_absolute_hdr_budget = 1e-6; +static const double psf_boundary_hdr_budget = 1e-7; +static const double pi = 3.14159265358979323846; +static const double gauss4_x[PSF_QUADRATURE_ORDER] = { + -0.8611363115940526, -0.3399810435848563, + 0.3399810435848563, 0.8611363115940526}; +static const double gauss4_w[PSF_QUADRATURE_ORDER] = { + 0.3478548451374539, 0.6521451548625461, + 0.6521451548625461, 0.3478548451374539}; +static const double gauss8_x[8] = { + -0.9602898564975363, -0.7966664774136267, -0.5255324099163290, + -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, + 0.7966664774136267, 0.9602898564975363}; +static const double gauss8_w[8] = { + 0.1012285362903763, 0.2223810344533745, 0.3137066458778873, + 0.3626837833783620, 0.3626837833783620, 0.3137066458778873, + 0.2223810344533745, 0.1012285362903763}; + +static int valid_psf(const PointSpreadFunction *psf) +{ + return psf != NULL && isfinite(psf->fwhm_pixels) && + isfinite(psf->moffat_beta) && psf->fwhm_pixels > 0.0 && + psf->moffat_beta > 1.0; +} + +static double moffat_alpha(const PointSpreadFunction *psf) +{ + return psf->fwhm_pixels / + (2.0 * sqrt(pow(2.0, 1.0 / psf->moffat_beta) - 1.0)); +} + +static double moffat_support_radius(double alpha, double beta, + double scaled_flux) +{ + const double normalization = scaled_flux * (beta - 1.0) / (pi * alpha * alpha); + const double relative_tail = fmin(psf_relative_tail_fraction, + psf_tail_absolute_hdr_budget / scaled_flux); + const double tail_radius = alpha * sqrt(pow(relative_tail, + 1.0 / (1.0 - beta)) - 1.0); + const double boundary_ratio = psf_boundary_hdr_budget / normalization; + const double boundary_radius = boundary_ratio >= 1.0 ? 0.0 : alpha * sqrt( + pow(boundary_ratio, -1.0 / beta) - 1.0); + return fmax(tail_radius, boundary_radius); +} + +static double moffat_pixel_integral_quadrature(double alpha, double beta, + double pixel_x, double pixel_y, + double star_x, double star_y, + const double *nodes, + const double *weights, int order) +{ + const double normalization = (beta - 1.0) / (pi * alpha * alpha); + double sum = 0.0; + for (int iy = 0; iy < order; ++iy) + for (int ix = 0; ix < order; ++ix) { + const double sx = pixel_x + 0.5 + 0.5 * nodes[ix] - star_x; + const double sy = pixel_y + 0.5 + 0.5 * nodes[iy] - star_y; + sum += 0.25 * weights[ix] * weights[iy] * normalization * + pow(1.0 + (sx * sx + sy * sy) / (alpha * alpha), -beta); + } + return sum; +} + +static double moffat_pixel_integral_cached(double alpha, double beta, + double pixel_x, double pixel_y, + double star_x, double star_y) +{ + return moffat_pixel_integral_quadrature(alpha, beta, pixel_x, pixel_y, + star_x, star_y, gauss4_x, gauss4_w, 4); +} + +static double moffat_pixel_integral_reference(double alpha, double beta, + double pixel_x, double pixel_y, + double star_x, double star_y) +{ + return moffat_pixel_integral_quadrature(alpha, beta, pixel_x, pixel_y, + star_x, star_y, gauss8_x, gauss8_w, 8); +} + +static size_t kernel_index(const PsfKernelCache *cache, int phase_x, + int phase_y, int offset_x, int offset_y) +{ + const size_t nodes = (size_t)cache->phase_resolution + 1; + const size_t side = (size_t)cache->radius_pixels * 2 + 1; + return (((size_t)phase_y * nodes + phase_x) * side + + (size_t)(offset_y + cache->radius_pixels)) * side + + (size_t)(offset_x + cache->radius_pixels); +} + +void psf_kernel_cache_destroy(PsfKernelCache *cache) +{ + if (cache == NULL) + return; + free(cache->weights); + *cache = (PsfKernelCache){0}; +} + +int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf) +{ + if (cache == NULL || !valid_psf(psf)) + return -1; + psf_kernel_cache_destroy(cache); + const double alpha = moffat_alpha(psf); + const double requested_radius = moffat_support_radius(alpha, psf->moffat_beta, 1.0); + if (!isfinite(requested_radius) || requested_radius <= 0.0) + return -1; + const int radius = (int)fmin((double)PSF_MAX_CACHE_RADIUS_PIXELS, + ceil(requested_radius) + 1.0); + const size_t nodes = PSF_PHASE_RESOLUTION + 1u; + const size_t side = (size_t)radius * 2 + 1u; + if (nodes > SIZE_MAX / nodes || nodes * nodes > SIZE_MAX / side || + nodes * nodes * side > SIZE_MAX / side || + nodes * nodes * side * side > SIZE_MAX / sizeof(float)) + return -1; + const size_t count = nodes * nodes * side * side; + float *weights = malloc(count * sizeof *weights); + if (weights == NULL) + return -1; + const clock_t start = clock(); + PsfKernelCache building = {.fwhm_pixels = psf->fwhm_pixels, + .moffat_beta = psf->moffat_beta, + .alpha_pixels = alpha, + .max_radius_pixels = radius, + .weights = weights, + .phase_resolution = PSF_PHASE_RESOLUTION, + .radius_pixels = radius}; + for (int phase_y = 0; phase_y <= PSF_PHASE_RESOLUTION; ++phase_y) + for (int phase_x = 0; phase_x <= PSF_PHASE_RESOLUTION; ++phase_x) { + const double star_x = (double)phase_x / PSF_PHASE_RESOLUTION; + const double star_y = (double)phase_y / PSF_PHASE_RESOLUTION; + double sum = 0.0; + for (int offset_y = -radius; offset_y <= radius; ++offset_y) + for (int offset_x = -radius; offset_x <= radius; ++offset_x) { + const size_t index = kernel_index(&building, phase_x, phase_y, + offset_x, offset_y); + const double weight = moffat_pixel_integral_cached( + alpha, psf->moffat_beta, offset_x, offset_y, star_x, star_y); + weights[index] = (float)weight; + sum += weight; + } + if (!(sum > 0.0) || !isfinite(sum)) { + free(weights); + return -1; + } + for (int offset_y = -radius; offset_y <= radius; ++offset_y) + for (int offset_x = -radius; offset_x <= radius; ++offset_x) { + const size_t index = kernel_index(&building, phase_x, phase_y, + offset_x, offset_y); + weights[index] = (float)(weights[index] / sum); + } + } + building.build_seconds = (double)(clock() - start) / CLOCKS_PER_SEC; + building.ready = 1; + *cache = building; + return 0; +} + +void psf_kernel_cache_report(const PsfKernelCache *cache, + const PsfSplatStats *stats, FILE *stream) +{ + if (stream == NULL) + return; + if (cache == NULL || !cache->ready) { + fprintf(stream, "PSF cache: disabled; cached splats 0, direct fallbacks %zu\n", + stats == NULL ? 0u : stats->direct_fallbacks); + return; + } + fprintf(stream, "PSF cache: 64x64 phases, radius %.0f px, tail abs %.0e, " + "boundary %.0e, build %.3f s; cached splats %zu, direct fallbacks %zu\n", + cache->max_radius_pixels, psf_tail_absolute_hdr_budget, + psf_boundary_hdr_budget, cache->build_seconds, + stats == NULL ? 0u : stats->cached_splats, + stats == NULL ? 0u : stats->direct_fallbacks); +} + +void splat_moffat_direct(double *hdr, int width, int height, double x, double y, + LinearRgb color, double flux, + const PointSpreadFunction *psf) +{ + if (hdr == NULL || width <= 0 || height <= 0 || flux <= 0.0 || !valid_psf(psf)) + return; + const double alpha = moffat_alpha(psf); + if (!isfinite(flux)) + return; + const double support_radius = moffat_support_radius(alpha, psf->moffat_beta, flux); + const int min_x = fmax(0.0, floor(x - support_radius)); + const int max_x = fmin((double)width - 1.0, ceil(x + support_radius)); + const int min_y = fmax(0.0, floor(y - support_radius)); + const int max_y = fmin((double)height - 1.0, ceil(y + support_radius)); + for (int py = min_y; py <= max_y; ++py) for (int px = min_x; px <= max_x; ++px) { + if (px < 0 || px >= width || py < 0 || py >= height) + continue; + const double dx = (px + 0.5) - x, dy = (py + 0.5) - y; + if (dx * dx + dy * dy > support_radius * support_radius) + continue; + const double weight = moffat_pixel_integral_reference(alpha, psf->moffat_beta, + px, py, x, y); + double *pixel = &hdr[3 * (py * width + px)]; + pixel[0] += color.r * flux * weight; + pixel[1] += color.g * flux * weight; + pixel[2] += color.b * flux * weight; + } +} + +int splat_moffat_cached(double *hdr, int width, int height, double x, double y, + LinearRgb color, double flux, + const PointSpreadFunction *psf, + const PsfKernelCache *cache) +{ + if (hdr == NULL || width <= 0 || height <= 0 || flux <= 0.0 || !valid_psf(psf)) + return 1; + const double alpha = moffat_alpha(psf); + const double support_radius = moffat_support_radius(alpha, psf->moffat_beta, + flux); + if (cache == NULL || !cache->ready || cache->fwhm_pixels != psf->fwhm_pixels || + cache->moffat_beta != psf->moffat_beta || !isfinite(support_radius) || + support_radius > cache->max_radius_pixels) { + splat_moffat_direct(hdr, width, height, x, y, color, flux, psf); + return 1; + } + const double base_x = floor(x), base_y = floor(y); + const double fx = x - base_x, fy = y - base_y; + const double phase_x = fx * cache->phase_resolution; + const double phase_y = fy * cache->phase_resolution; + const int x0 = (int)floor(phase_x), y0 = (int)floor(phase_y); + const int x1 = x0 + 1, y1 = y0 + 1; + const double tx = phase_x - x0, ty = phase_y - y0; + const int support = (int)ceil(support_radius); + for (int offset_y = -support; offset_y <= support; ++offset_y) + for (int offset_x = -support; offset_x <= support; ++offset_x) { + const int px = (int)base_x + offset_x, py = (int)base_y + offset_y; + if (px < 0 || px >= width || py < 0 || py >= height) + continue; + const double dx = offset_x + 0.5 - fx, dy = offset_y + 0.5 - fy; + if (dx * dx + dy * dy > support_radius * support_radius) + continue; + const double w00 = cache->weights[kernel_index(cache, x0, y0, offset_x, offset_y)]; + const double w10 = cache->weights[kernel_index(cache, x1, y0, offset_x, offset_y)]; + const double w01 = cache->weights[kernel_index(cache, x0, y1, offset_x, offset_y)]; + const double w11 = cache->weights[kernel_index(cache, x1, y1, offset_x, offset_y)]; + const double weight = (1.0 - ty) * ((1.0 - tx) * w00 + tx * w10) + + ty * ((1.0 - tx) * w01 + tx * w11); + double *pixel = &hdr[3 * (py * width + px)]; + pixel[0] += color.r * flux * weight; + pixel[1] += color.g * flux * weight; + pixel[2] += color.b * flux * weight; + } + return 0; +} + void splat_moffat(double *hdr, int width, int height, double x, double y, LinearRgb color, double flux, const PointSpreadFunction *psf) { - /* The tail omitted outside this radius contains 1e-8 of the Moffat's - * total flux. Unlike the old fixed 3-sigma box, this is both circular and - * far below the displayed HDR precision for the chosen beta. */ - const double tail_fraction = 1e-8; - if (hdr == NULL || psf == NULL || flux <= 0.0 || - psf->fwhm_pixels <= 0.0 || psf->moffat_beta <= 1.0) - return; - const double beta = psf->moffat_beta; - const double alpha = psf->fwhm_pixels / - (2.0 * sqrt(pow(2.0, 1.0 / beta) - 1.0)); - const double support_radius = alpha * sqrt( - pow(tail_fraction, 1.0 / (1.0 - beta)) - 1.0); - const double support_radius_squared = support_radius * support_radius; - const int min_x = (int)floor(x - support_radius); - const int max_x = (int)ceil(x + support_radius); - const int min_y = (int)floor(y - support_radius); - const int max_y = (int)ceil(y + support_radius); - const double normalization = flux * (beta - 1.0) / - (3.14159265358979323846 * alpha * alpha); - for (int py = min_y; py <= max_y; ++py) for (int px = min_x; px <= max_x; ++px) { - if (px < 0 || px >= width || py < 0 || py >= height) continue; - double dx = (px + 0.5) - x, dy = (py + 0.5) - y; - const double radius_squared = dx * dx + dy * dy; - if (radius_squared > support_radius_squared) continue; - const double w = normalization * - pow(1.0 + radius_squared / (alpha * alpha), -beta); - double *pixel = &hdr[3 * (py * width + px)]; - pixel[0] += color.r * w; pixel[1] += color.g * w; pixel[2] += color.b * w; - } + splat_moffat_direct(hdr, width, height, x, y, color, flux, psf); } static unsigned char tonemap_channel(double hdr_value) @@ -186,3 +421,30 @@ int write_tonemapped_image(const char *path, const double *hdr, int width, int h return -1; #endif } + +#ifdef ENABLE_HDR_DEBUG +int write_hdr_pfm(const char *path, const double *hdr, int width, int height) +{ + if (path == NULL || hdr == NULL || width <= 0 || height <= 0) + return -1; + FILE *file = fopen(path, "wb"); + if (file == NULL) + return -1; + const uint16_t endian_probe = 1; + const char *scale = *(const unsigned char *)&endian_probe == 1 ? "-1.0" : "1.0"; + int result = fprintf(file, "PF\n%d %d\n%s\n", width, height, scale) < 0 ? -1 : 0; + /* PFM rows are stored bottom-to-top. Its negative scale declares little-endian + * float samples, avoiding an unnecessary byte swap on the normal test host. */ + for (int row = height - 1; result == 0 && row >= 0; --row) + for (int column = 0; column < width * 3; ++column) { + const float sample = (float)hdr[(size_t)row * width * 3 + column]; + if (fwrite(&sample, sizeof sample, 1, file) != 1) { + result = -1; + break; + } + } + if (fclose(file) != 0) + result = -1; + return result; +} +#endif diff --git a/src/optics.h b/src/optics.h index 40ea112..3b809dc 100644 --- a/src/optics.h +++ b/src/optics.h @@ -1,20 +1,57 @@ #ifndef OPTICS_H #define OPTICS_H +#include +#include + typedef struct { double r, g, b; } LinearRgb; typedef struct { double fwhm_pixels; double moffat_beta; } PointSpreadFunction; +/* One immutable process-wide kernel for the one PSF parameter pair accepted + * by the current renderer. Its storage remains private to optics.c. */ +typedef struct { + double fwhm_pixels, moffat_beta, alpha_pixels; + double max_radius_pixels, build_seconds; + float *weights; + int phase_resolution, radius_pixels; + int ready; +} PsfKernelCache; + +typedef struct { + size_t cached_splats; + size_t direct_fallbacks; +} PsfSplatStats; + /* Integrate a Planck spectrum into absolute linear-sRGB spectral radiance * (W m^-2 sr^-1), before catalog amplitude and display exposure. */ LinearRgb blackbody_to_linear_rgb(double temperature_K); +int psf_kernel_cache_init(PsfKernelCache *cache, + const PointSpreadFunction *psf); +void psf_kernel_cache_destroy(PsfKernelCache *cache); +void psf_kernel_cache_report(const PsfKernelCache *cache, + const PsfSplatStats *stats, FILE *stream); +/* Reference implementation: pixel-area-integrated Moffat with the same tail + * budgets as the cache-aware renderer. */ +void splat_moffat_direct(double *hdr, int width, int height, double x, double y, + LinearRgb color, double flux, + const PointSpreadFunction *psf); void splat_moffat(double *hdr, int width, int height, double x, double y, LinearRgb color, double flux, const PointSpreadFunction *psf); +/* Returns nonzero when this event used the direct reference fallback. */ +int splat_moffat_cached(double *hdr, int width, int height, double x, double y, + LinearRgb color, double flux, + const PointSpreadFunction *psf, + const PsfKernelCache *cache); /* Writes PNG when built with libpng; non-libpng builds use PPM fallback. */ int write_tonemapped_image(const char *path, const double *hdr, int width, int height); +#ifdef ENABLE_HDR_DEBUG +/* Test-build-only: writes the pre-tone-mapping framebuffer as RGB float PFM. */ +int write_hdr_pfm(const char *path, const double *hdr, int width, int height); +#endif #endif diff --git a/tests/test_frame.c b/tests/test_frame.c index a194de4..2169b88 100644 --- a/tests/test_frame.c +++ b/tests/test_frame.c @@ -27,7 +27,7 @@ int main(void) { goto done; const size_t images = frame_splat_catalog(&mesh, &catalog, hdr, width, height, test_exposure, - &psf, 1, NULL); + &psf, NULL, 1, NULL, NULL); if (images != 1 || hdr[3 * (50 * width + 50)] <= 0.0) { fputs("flat-space inverse lens-map regression failed\n", stderr); goto done; @@ -44,10 +44,10 @@ int main(void) { omp_set_dynamic(0); omp_set_num_threads(1); const size_t serial_images = frame_splat_catalog( - &mesh, &catalog, serial_hdr, width, height, test_exposure, &psf, 1, NULL); + &mesh, &catalog, serial_hdr, width, height, test_exposure, &psf, NULL, 1, NULL, NULL); omp_set_num_threads(4); const size_t parallel_images = frame_splat_catalog( - &mesh, &catalog, parallel_hdr, width, height, test_exposure, &psf, 1, NULL); + &mesh, &catalog, parallel_hdr, width, height, test_exposure, &psf, NULL, 1, NULL, NULL); omp_set_num_threads(original_threads); for (int value = 0; value < width * height * 3; ++value) if (fabs(serial_hdr[value] - parallel_hdr[value]) > @@ -83,6 +83,56 @@ int main(void) { fputs("Moffat normalization or wing regression failed\n", stderr); goto done; } + /* The cache stores 4-point pixel-area integrals over a 64x64 sub-pixel + * lattice. Compare its bilinear interpolation with the independent 8-point + * direct reference at phases on both sides of a pixel boundary. */ + PsfKernelCache cache = {0}; + double *cached_hdr = calloc((size_t)width * height * 3, sizeof *cached_hdr); + double *reference_hdr = calloc((size_t)width * height * 3, sizeof *reference_hdr); + if (cached_hdr == NULL || reference_hdr == NULL || + psf_kernel_cache_init(&cache, &psf)) { + free(cached_hdr); + free(reference_hdr); + psf_kernel_cache_destroy(&cache); + fputs("PSF cache construction regression failed\n", stderr); + goto done; + } + const double phases[][2] = {{0.01, 0.99}, {0.499, 0.501}, {0.999, 0.001}}; + for (size_t phase = 0; phase < sizeof phases / sizeof *phases; ++phase) { + memset(cached_hdr, 0, (size_t)width * height * 3 * sizeof *cached_hdr); + memset(reference_hdr, 0, (size_t)width * height * 3 * sizeof *reference_hdr); + if (splat_moffat_cached(cached_hdr, width, height, 50.0 + phases[phase][0], + 50.0 + phases[phase][1], (LinearRgb){1.0, 1.0, 1.0}, + 1.0, &psf, &cache) != 0) { + fputs("ordinary PSF cache unexpectedly fell back\n", stderr); + free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); + goto done; + } + splat_moffat_direct(reference_hdr, width, height, + 50.0 + phases[phase][0], 50.0 + phases[phase][1], + (LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf); + double peak = 0.0, max_error = 0.0; + for (int value = 0; value < width * height * 3; ++value) { + peak = fmax(peak, reference_hdr[value]); + max_error = fmax(max_error, fabs(cached_hdr[value] - reference_hdr[value])); + } + if (peak <= 0.0 || max_error > 4e-5 * peak) { + fputs("PSF cache interpolation accuracy regression failed\n", stderr); + free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); + goto done; + } + } + /* A bright event must avoid a cached hard cutoff by selecting the direct + * reference path when the requested support exceeds the cache. */ + if (splat_moffat_cached(cached_hdr, width, height, 50.5, 50.5, + (LinearRgb){1.0, 1.0, 1.0}, 1000.0, &psf, &cache) != 1) { + fputs("bright PSF direct-fallback regression failed\n", stderr); + free(cached_hdr); free(reference_hdr); psf_kernel_cache_destroy(&cache); + goto done; + } + free(cached_hdr); + free(reference_hdr); + psf_kernel_cache_destroy(&cache); frame_draw_mesh(&mesh, hdr, width, height, 0.5, 0.5); if (hdr[3 * (10 * width + 20)] != 0.25) { fputs("mesh diagnostic overlay regression failed\n", stderr); @@ -102,7 +152,7 @@ int main(void) { if (frame_lens_mesh_build_coarse(&fine_mesh, width, height, 1, 0.1) || frame_lens_mesh_trace(&fine_mesh, &spacetime, &observer, &trace) || frame_splat_catalog(&fine_mesh, &fine_catalog, hdr, width, height, - test_exposure, &psf, 1, NULL) != 1) { + test_exposure, &psf, NULL, 1, NULL, NULL) != 1) { fputs("fine source-triangle containment regression failed\n", stderr); frame_lens_mesh_destroy(&fine_mesh); goto done;