202 lines
9.0 KiB
C
202 lines
9.0 KiB
C
#include "frame.h"
|
|
#include "optics.h"
|
|
|
|
#include <math.h>
|
|
#include <omp.h>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
|
|
int main(void) {
|
|
const int width = 100, height = 100;
|
|
const double test_exposure = 1e-3;
|
|
const PointSpreadFunction psf = {.fwhm_pixels = 2.7, .moffat_beta = 4.5};
|
|
const GeodesicTraceConfig trace = {.coordinate_time_step = 0.25,
|
|
.max_steps = 100};
|
|
const ObserverState observer = observer_fixed_at_origin();
|
|
Star star = {
|
|
.direction = {0.0, 0.0, -1.0}, .temperature_K = 7000.0, .amplitude = 1.0};
|
|
StarCatalog catalog = {.stars = &star, .count = 1};
|
|
SpacetimeSource spacetime = {0};
|
|
FrameLensMesh mesh = {0};
|
|
double *hdr = calloc((size_t)width * height * 3, sizeof *hdr);
|
|
int result = 1;
|
|
if (hdr == NULL || spacetime_create_minkowski(&spacetime, 10.0) ||
|
|
frame_lens_mesh_build_coarse(&mesh, width, height, 20, 30.0) ||
|
|
frame_lens_mesh_trace(&mesh, &spacetime, &observer, &trace))
|
|
goto done;
|
|
const size_t images =
|
|
frame_splat_catalog(&mesh, &catalog, hdr, width, height, test_exposure,
|
|
&psf, NULL, 0, 1, NULL, NULL, NULL);
|
|
if (images != 1 || hdr[3 * (50 * width + 50)] <= 0.0) {
|
|
fputs("flat-space inverse lens-map regression failed\n", stderr);
|
|
goto done;
|
|
}
|
|
/* Private HDR accumulation must preserve the serial splat result. */
|
|
double *serial_hdr = calloc((size_t)width * height * 3, sizeof *serial_hdr);
|
|
double *parallel_hdr = calloc((size_t)width * height * 3, sizeof *parallel_hdr);
|
|
if (serial_hdr == NULL || parallel_hdr == NULL) {
|
|
free(serial_hdr);
|
|
free(parallel_hdr);
|
|
goto done;
|
|
}
|
|
const int original_threads = omp_get_max_threads();
|
|
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, NULL, 0,
|
|
1, NULL, NULL, NULL);
|
|
omp_set_num_threads(4);
|
|
const size_t parallel_images = frame_splat_catalog(
|
|
&mesh, &catalog, parallel_hdr, width, height, test_exposure, &psf, NULL, 0,
|
|
1, NULL, 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]) >
|
|
1e-12 * fmax(1.0, fabs(serial_hdr[value]))) {
|
|
fputs("parallel catalog splat regression failed\n", stderr);
|
|
free(serial_hdr);
|
|
free(parallel_hdr);
|
|
goto done;
|
|
}
|
|
free(serial_hdr);
|
|
free(parallel_hdr);
|
|
if (serial_images != 1 || parallel_images != serial_images) {
|
|
fputs("parallel catalog image-count regression failed\n", stderr);
|
|
goto done;
|
|
}
|
|
const LinearRgb cool = blackbody_to_linear_rgb(3000.0);
|
|
const LinearRgb hot = blackbody_to_linear_rgb(10000.0);
|
|
if (!(cool.r > cool.b && hot.b > hot.r &&
|
|
hot.r + hot.g + hot.b > cool.r + cool.g + cool.b)) {
|
|
fputs("blackbody spectral-color regression failed\n", stderr);
|
|
goto done;
|
|
}
|
|
/* The Moffat is flux-normalized and retains a measurable, continuous wing
|
|
* beyond the former Gaussian's 3-sigma raster box. */
|
|
memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr);
|
|
splat_moffat(hdr, width, height, 50.5, 50.5,
|
|
(LinearRgb){1.0, 1.0, 1.0}, 1.0, &psf);
|
|
double moffat_flux = 0.0;
|
|
for (int pixel = 0; pixel < width * height; ++pixel)
|
|
moffat_flux += hdr[3 * pixel];
|
|
if (fabs(moffat_flux - 1.0) > 0.01 ||
|
|
hdr[3 * (50 * width + 62)] <= 0.0) {
|
|
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);
|
|
goto done;
|
|
}
|
|
/* A fixed absolute edge tolerance used to make tiny source triangles claim
|
|
* sources far outside their field. */
|
|
FrameLensMesh fine_mesh = {0};
|
|
Star fine_stars[2] = {{.direction = {0.0, 0.0, -1.0},
|
|
.temperature_K = 7000.0,
|
|
.amplitude = 1.0},
|
|
{.direction = {0.01, 0.0, -0.9999499987499375},
|
|
.temperature_K = 7000.0,
|
|
.amplitude = 1.0}};
|
|
StarCatalog fine_catalog = {.stars = fine_stars, .count = 2};
|
|
memset(hdr, 0, (size_t)width * height * 3 * sizeof *hdr);
|
|
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, NULL, 0, 1, NULL, NULL, NULL) != 1) {
|
|
fputs("fine source-triangle containment regression failed\n", stderr);
|
|
frame_lens_mesh_destroy(&fine_mesh);
|
|
goto done;
|
|
}
|
|
frame_lens_mesh_destroy(&fine_mesh);
|
|
/* 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. */
|
|
FrameLensMesh adaptive_mesh = {0};
|
|
const RefinementConfig refine = {.max_level = 1,
|
|
.angle_absolute_rad = 1e-4,
|
|
.angle_relative = 1e-4,
|
|
.min_edge_pixels = 1.0,
|
|
.min_area_pixels2 = 1.0};
|
|
if (frame_lens_mesh_build_coarse(&adaptive_mesh, width, height, 100, 30.0))
|
|
goto done;
|
|
for (size_t i = 0; i < adaptive_mesh.vertex_count; ++i) {
|
|
adaptive_mesh.vertices[i].traced = 1;
|
|
adaptive_mesh.vertices[i].status = RAY_ENDPOINT_ESCAPED;
|
|
adaptive_mesh.vertices[i].n_infinity[0] = 1.0;
|
|
}
|
|
if (frame_lens_mesh_prepare_generation(&adaptive_mesh, &refine) != 1) {
|
|
fputs("adaptive shared-edge probe setup regression failed\n", stderr);
|
|
frame_lens_mesh_destroy(&adaptive_mesh);
|
|
goto done;
|
|
}
|
|
const RayEndpoint bent_probe = {.n_infinity = {0.0, 1.0, 0.0},
|
|
.frequency_ratio = 1.0,
|
|
.status = RAY_ENDPOINT_ESCAPED};
|
|
if (frame_lens_mesh_install_sample(&adaptive_mesh, 0, &bent_probe) ||
|
|
frame_lens_mesh_finish_generation(&adaptive_mesh, &refine) != 1 ||
|
|
adaptive_mesh.vertex_count != 5 || adaptive_mesh.triangle_count != 4) {
|
|
fputs("adaptive shared-edge split regression failed\n", stderr);
|
|
frame_lens_mesh_destroy(&adaptive_mesh);
|
|
goto done;
|
|
}
|
|
frame_lens_mesh_destroy(&adaptive_mesh);
|
|
result = 0;
|
|
done:
|
|
frame_lens_mesh_destroy(&mesh);
|
|
spacetime_destroy(&spacetime);
|
|
free(hdr);
|
|
return result;
|
|
}
|