diff --git a/mk/reference_images.mk b/mk/reference_images.mk index 9aa408d..40f26d3 100644 --- a/mk/reference_images.mk +++ b/mk/reference_images.mk @@ -13,7 +13,4 @@ test-reference-images: $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/minkowski_ra1_dec1_6 python3 $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/minkowski_ra1_dec1_640x360_HDR.fits $(REFERENCE_TMP_DIR)/minkowski_ra1_dec1_640x360_HDR.fits $(MAKE) SPACETIME=schwarzschild ENABLE_HDR=1 backend OMP_NUM_THREADS=16 $(BUILD_DIR)/schwarzschild_sky --catalog assets/sky_grid_5deg.csv --output $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360.png --hdr-output --width 640 --height 360 --fov-deg 60 --look-ra-deg 1 --look-dec-deg 1 --exposure 0.1 --observer-radius 30 --observer-velocity 0 0 0 --psf-fwhm-pixels 2.7 --psf-moffat-beta 4.5 --max-magnification 1e300 --max-cache-psf-flux 1 --psf-relative-tail 1e-8 --psf-min-y 0 --coarse-cell-pixels 16 --refine-max-level 0 --refine-angle-abs-deg 0.001 --refine-angle-rel 0.1 --refine-jacobian-min 1e-3 --refine-min-edge-pixels 0.5 --refine-min-area-pixels2 0.25 --catalog-load-workers 4 - # The general metric-based tetrad changes four float32 samples by one ULP - # versus the legacy analytic static tetrad. Keep the original fixture and - # allow at most float32 relative rounding (2^-23), with zero abs tolerance. - python3 $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits --rel-tolerance 1.1920928955078125e-7 + python3 $(FLOATDIFF_SCRIPT) $(REFERENCE_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits $(REFERENCE_TMP_DIR)/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index 9cfffe5..135fc90 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -250,6 +250,21 @@ typedef struct { # 8. 局部 adaptive refinement +局部 inverse map 的球面面积权重必须满足凸组合约束。当前实现以无符号 +子面积除以总面积得到权重,在插值位置和频移前检查权重和有限、为正,且 +`abs(sum(weights) - 1) <= 1e-8`;不满足时只跳过当前恒星在当前 triangle +中的像。通过检查后,将权重除以其总和,消除小的归一化误差。整 tile +包含的 catalog 快速路径也必须执行这项检查。 + +这是反演有效性保护,不是 refinement 阈值:细长 source triangle 的带容差 +边测试可能误接纳外部点,无符号子面积之和便大于总面积,直接插值会把像 +放到 image triangle 外。不能先归一化明显无效的权重来掩盖误接纳。 +固定的无量纲容差 `1e-8` 依据 Schwarzschild 3840×2160、45° FOV、R=100、 +coarse cell=32、level=4、J=0.2 网格上的诊断选定:正对中心及偏离 +0.00161°、方位为 1°/44°/23° 的单星测试中,严格包含样本的 double +权重和误差不超过约 `2.5e-13`,正对中心的越界样本最小误差约 `5.9e-4`。 +这不保证临界曲线附近的离散映射已经收敛,也不保证数值亮环无缺口。 + 每个 triangle 只关心自身局部映射。 不需要先把整个像平面映射到天球、再全局分类“第几阶像”。 diff --git a/src/frame.c b/src/frame.c index 1612d71..56901e1 100644 --- a/src/frame.c +++ b/src/frame.c @@ -930,6 +930,17 @@ static int spherical_barycentric_weights(const double point[3], const double a[3 weights[0] = spherical_area(point, b, c) / area; weights[1] = spherical_area(point, c, a) / area; weights[2] = spherical_area(point, a, b) / area; + /* A thin source triangle can admit an exterior point through the tolerant + * edge test. Unsigned subareas then do not partition the total area, and + * interpolating with their ratios can move an image outside its triangle. + * Reject that case before normalizing roundoff in valid convex weights. + * The 1e-8 bound is dimensionless; see the inverse-map design notes. */ + const double weight_sum = weights[0] + weights[1] + weights[2]; + if (!isfinite(weight_sum) || weight_sum <= 0.0 || + fabs(weight_sum - 1.0) > 1e-8) + return -1; + for (int i = 0; i < 3; ++i) + weights[i] /= weight_sum; return 0; } @@ -1055,7 +1066,7 @@ static int splat_catalog_tile(const Star *stars, size_t count, context->vertex[0]->n_infinity, context->vertex[1]->n_infinity, context->vertex[2]->n_infinity, weights)) - return -1; + continue; } if (!owns_source_boundary(context->triangle, weights)) continue; diff --git a/tests/data/psf_event_sink_reference/minkowski_ra1_dec1_640x360_HDR.fits b/tests/data/psf_event_sink_reference/minkowski_ra1_dec1_640x360_HDR.fits index 678196a..064bfd2 100644 Binary files a/tests/data/psf_event_sink_reference/minkowski_ra1_dec1_640x360_HDR.fits and b/tests/data/psf_event_sink_reference/minkowski_ra1_dec1_640x360_HDR.fits differ diff --git a/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits b/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits index 0aca09c..8b12615 100644 Binary files a/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits and b/tests/data/psf_event_sink_reference/schwarzschild_ra1_dec1_fov60_640x360_HDR.fits differ diff --git a/tests/test_frame.c b/tests/test_frame.c index e472bd0..ccd08b3 100644 --- a/tests/test_frame.c +++ b/tests/test_frame.c @@ -330,6 +330,50 @@ int main(void) { goto done; } frame_lens_mesh_destroy(&fine_mesh); + /* Schwarzschild level-4/J=0.2 ring triangle: the old edge tolerance accepts + * (1,0,0) although it is outside, producing unsigned weights summing to + * 1.09608. Keep a genuine interior source after it to check that rejecting + * one source does not discard subsequent stars. Exercise both parities. */ + const double thin_directions[3][3] = { + {0.99999998891116071, 0.00013431364716360775, 6.4323579270168807e-05}, + {0.99999997228355786, -0.00021216644782556991, -0.00010206998544149922}, + {0.9999927016945267, -0.0034455730513416835, -0.001650631403158936}}; + LensVertex thin_vertices[3] = { + {.image_x = 40, .image_y = 40, .status = RAY_ENDPOINT_ESCAPED}, + {.image_x = 48, .image_y = 40, .status = RAY_ENDPOINT_ESCAPED}, + {.image_x = 40, .image_y = 48, .status = RAY_ENDPOINT_ESCAPED}}; + LensTriangle thin_triangle = {.vertex = {0, 1, 2}}; + FrameLensMesh thin_mesh = {.vertices = thin_vertices, .vertex_count = 3, + .triangles = &thin_triangle, .triangle_count = 1}; + Star thin_stars[2] = { + {.direction = {1, 0, 0}, .temperature_K = 7000, .amplitude = 1}, + {.temperature_K = 7000, .amplitude = 1}}; + for (int i = 0; i < 3; ++i) { + memcpy(thin_vertices[i].n_infinity, thin_directions[i], + sizeof thin_directions[i]); + memcpy(thin_vertices[i].camera_direction, thin_directions[i], + sizeof thin_directions[i]); + for (int axis = 0; axis < 3; ++axis) + thin_stars[1].direction[axis] += thin_directions[i][axis]; + } + const double thin_norm = hypot(hypot(thin_stars[1].direction[0], + thin_stars[1].direction[1]), + thin_stars[1].direction[2]); + for (int axis = 0; axis < 3; ++axis) + thin_stars[1].direction[axis] /= thin_norm; + StarCatalog thin_catalog = {.stars = thin_stars, .count = 2}; + for (int parity = 0; parity < 2; ++parity) { + thin_triangle.vertex[1] = parity ? 2 : 1; + thin_triangle.vertex[2] = parity ? 1 : 2; + memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr); + if (frame_splat_catalog(&thin_mesh, &thin_catalog, hdr, width, height, + test_exposure, &psf, NULL, 1.0, 1.0, + psf_relative_tail, 0.0, 0, 1, NULL, NULL, NULL) != 1 || + hdr[3 * (43 * width + 43)] <= 0.0) { + fputs("thin source-triangle inverse-map regression failed\n", stderr); + goto done; + } + } /* Refinement probes are temporary until their generation is complete. A * shared diagonal probe must produce one stable midpoint and conforming * children only after its endpoint has been installed. */