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

Optics: add PSF luminance cutoff

Browse Source
This commit is contained in:
wyj committed 2026-08-30 15:10:14 -04:00
1 parent 9e23c3fcdb
commit e8deb1ff19
7 files changed
+123 -38
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

src/frame.c
+29 -13
View File
@@ -851,15 +851,18 @@ typedef struct {
const PsfKernelCache *psf_cache;
double max_cache_psf_flux;
double psf_relative_tail;
double psf_min_y;
size_t images;
size_t direct_fallbacks;
size_t cached_wing_clipped;
size_t discarded_below_min_y;
} TriangleSplatContext;
typedef struct {
size_t images;
size_t direct_fallbacks;
size_t cached_wing_clipped;
size_t discarded_below_min_y;
#ifdef GR_DEBUG
double max_raw_magnification;
size_t magnification_clamped_triangles;
@@ -872,9 +875,11 @@ static void copy_psf_splat_stats(PsfSplatStats *destination,
if (destination == NULL)
return;
*destination = (PsfSplatStats){.cached_splats =
source.images - source.direct_fallbacks,
source.images - source.direct_fallbacks -
source.discarded_below_min_y,
.cached_wing_clipped = source.cached_wing_clipped,
.direct_fallbacks = source.direct_fallbacks,
.discarded_below_min_y = source.discarded_below_min_y,
#ifdef GR_DEBUG
.max_raw_magnification = source.max_raw_magnification,
.magnification_clamped_triangles =
@@ -919,9 +924,10 @@ static int splat_catalog_tile(const Star *stars, size_t count,
context->hdr, context->width, context->height, image_x, image_y, color,
flux,
context->psf, context->psf_cache, context->max_cache_psf_flux,
context->psf_relative_tail);
context->psf_relative_tail, context->psf_min_y);
context->direct_fallbacks += direct_fallback == 1;
context->cached_wing_clipped += direct_fallback == 2;
context->discarded_below_min_y += direct_fallback == 3;
#ifdef GR_DEBUG
if (direct_fallback == 1)
/* This is deliberately emitted by the active splat worker: a direct
@@ -944,6 +950,7 @@ static CatalogSplatStats splat_catalog_triangles(
int height, double exposure, const PointSpreadFunction *psf,
const PsfKernelCache *psf_cache, double max_magnification,
double max_cache_psf_flux, double psf_relative_tail,
double psf_min_y,
size_t first_triangle, size_t last_triangle) {
CatalogSplatStats stats = {0};
for (size_t t = first_triangle; t < last_triangle; ++t) {
@@ -978,12 +985,14 @@ static CatalogSplatStats splat_catalog_triangles(
.exposure = exposure, .magnification = magnification,
.psf = psf, .psf_cache = psf_cache,
.max_cache_psf_flux = max_cache_psf_flux,
.psf_relative_tail = psf_relative_tail};
.psf_relative_tail = psf_relative_tail,
.psf_min_y = psf_min_y};
if (catalog_visit_source_triangle(catalog, direction, 0, splat_catalog_tile,
&context) == 0)
stats.images += context.images;
stats.direct_fallbacks += context.direct_fallbacks;
stats.cached_wing_clipped += context.cached_wing_clipped;
stats.discarded_below_min_y += context.discarded_below_min_y;
}
return stats;
}
@@ -1016,6 +1025,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
double max_magnification,
double max_cache_psf_flux,
double psf_relative_tail,
double psf_min_y,
int limit_workers_by_memory,
int catalog_load_workers,
CatalogPrefetchStats *prefetch_stats,
@@ -1026,7 +1036,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
isnan(max_magnification) || max_magnification <= 0.0 ||
!isfinite(max_cache_psf_flux) || max_cache_psf_flux < 1.0 ||
!isfinite(psf_relative_tail) || psf_relative_tail <= 0.0 ||
psf_relative_tail >= 1.0)
psf_relative_tail >= 1.0 || !isfinite(psf_min_y) || psf_min_y < 0.0)
return 0;
/* A bounded parallel read phase completes before splatting. Its serial cache
@@ -1050,7 +1060,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
{
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
0, mesh->triangle_count);
copy_psf_splat_stats(psf_stats, stats);
return stats.images;
@@ -1067,7 +1077,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
{
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
0, mesh->triangle_count);
copy_psf_splat_stats(psf_stats, stats);
return stats.images;
@@ -1078,7 +1088,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
{
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
0, mesh->triangle_count);
copy_psf_splat_stats(psf_stats, stats);
return stats.images;
@@ -1095,13 +1105,14 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
free(private_hdr);
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, hdr, width, height, exposure, psf, psf_cache,
max_magnification, max_cache_psf_flux, psf_relative_tail,
max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y,
0, mesh->triangle_count);
copy_psf_splat_stats(psf_stats, stats);
return stats.images;
}
size_t images = 0, direct_fallbacks = 0, cached_wing_clipped = 0;
size_t images = 0, direct_fallbacks = 0, cached_wing_clipped = 0,
discarded_below_min_y = 0;
#ifdef GR_DEBUG
double max_raw_magnification = 0.0;
size_t magnification_clamped_triangles = 0;
@@ -1109,9 +1120,9 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
/* Keep the ordinary render loop byte-for-byte free of progress checks. */
if (progress != NULL && progress->worker_callback != NULL) {
#ifdef GR_DEBUG
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, magnification_clamped_triangles) reduction(max : max_raw_magnification)
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y, magnification_clamped_triangles) reduction(max : max_raw_magnification)
#else
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped)
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y)
#endif
{
const size_t worker = (size_t)omp_get_thread_num();
@@ -1126,10 +1137,12 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, private_hdr[worker], width, height, exposure, psf,
psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail,
psf_min_y,
triangle, triangle + 1);
images += stats.images;
direct_fallbacks += stats.direct_fallbacks;
cached_wing_clipped += stats.cached_wing_clipped;
discarded_below_min_y += stats.discarded_below_min_y;
#ifdef GR_DEBUG
max_raw_magnification = fmax(max_raw_magnification, stats.max_raw_magnification);
magnification_clamped_triangles += stats.magnification_clamped_triangles;
@@ -1147,9 +1160,9 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
}
} else {
#ifdef GR_DEBUG
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, magnification_clamped_triangles) reduction(max : max_raw_magnification)
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y, magnification_clamped_triangles) reduction(max : max_raw_magnification)
#else
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped)
#pragma omp parallel num_threads((int)worker_count) reduction(+ : images, direct_fallbacks, cached_wing_clipped, discarded_below_min_y)
#endif
{
const size_t worker = (size_t)omp_get_thread_num();
@@ -1158,10 +1171,12 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
const CatalogSplatStats stats = splat_catalog_triangles(
mesh, catalog, private_hdr[worker], width, height, exposure, psf,
psf_cache, max_magnification, max_cache_psf_flux, psf_relative_tail,
psf_min_y,
triangle, triangle + 1);
images += stats.images;
direct_fallbacks += stats.direct_fallbacks;
cached_wing_clipped += stats.cached_wing_clipped;
discarded_below_min_y += stats.discarded_below_min_y;
#ifdef GR_DEBUG
max_raw_magnification = fmax(max_raw_magnification, stats.max_raw_magnification);
magnification_clamped_triangles += stats.magnification_clamped_triangles;
@@ -1180,6 +1195,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
.images = images,
.direct_fallbacks = direct_fallbacks,
.cached_wing_clipped = cached_wing_clipped,
.discarded_below_min_y = discarded_below_min_y,
#ifdef GR_DEBUG
.max_raw_magnification = max_raw_magnification,
.magnification_clamped_triangles =
src/frame.h
+1
View File
@@ -127,6 +127,7 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh,
double max_magnification,
double max_cache_psf_flux,
double psf_relative_tail,
double psf_min_y,
int limit_workers_by_memory,
int catalog_load_workers,
CatalogPrefetchStats *prefetch_stats,
src/main.c
+7 -1
View File
@@ -25,6 +25,7 @@ typedef struct {
double max_magnification;
double max_cache_psf_flux;
double psf_relative_tail;
double psf_min_y;
double observer_radius;
double observer_inward_speed;
PointSpreadFunction psf;
@@ -188,6 +189,7 @@ static int parse_args(int argc, char **argv, Settings *s,
.max_magnification = INFINITY,
.max_cache_psf_flux = 1.0,
.psf_relative_tail = 1e-8,
.psf_min_y = 0.0,
.observer_radius = 30.0,
.psf = {2.7, 4.5},
.catalog_path = "assets/sky_grid_5deg.csv",
@@ -260,6 +262,8 @@ static int parse_args(int argc, char **argv, Settings *s,
} else if (!strcmp(argv[i], "--psf-relative-tail") && i + 1 < argc &&
!parse_finite_positive(argv[++i], &s->psf_relative_tail) &&
s->psf_relative_tail < 1.0) {
} else if (!strcmp(argv[i], "--psf-min-y") && i + 1 < argc &&
!parse_nonnegative(argv[++i], &s->psf_min_y)) {
} else if (!strcmp(argv[i], "--psf-direct")) {
s->psf_direct = 1;
} else if (!strcmp(argv[i], "--verbose")) {
@@ -469,6 +473,7 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog,
&mesh, catalog, hdr, s->width, s->height, s->exposure, &s->psf,
&s->psf_cache, s->max_magnification, s->max_cache_psf_flux,
s->psf_relative_tail,
s->psf_min_y,
spacetime_limits_render_workers_by_memory(spacetime),
s->catalog_load_workers, &prefetch, &psf_stats,
&(FrameSplatProgress){report_splat_progress,
@@ -648,6 +653,7 @@ static int render_movie(const Settings *s, StarCatalog *catalog,
&movie.frames[i].mesh, catalog, hdr, s->width, s->height, s->exposure,
&s->psf, &s->psf_cache, s->max_magnification, s->max_cache_psf_flux,
s->psf_relative_tail,
s->psf_min_y,
spacetime_limits_render_workers_by_memory(spacetime),
s->catalog_load_workers, &prefetch, &psf_stats,
&(FrameSplatProgress){report_splat_progress,
@@ -693,7 +699,7 @@ int main(int argc, char **argv) {
"[--exposure E] [--observer-radius R] [--observer-inward-speed V] "
"[--psf-fwhm-pixels N] [--psf-moffat-beta N] "
"[--max-magnification M] [--max-cache-psf-flux F] "
"[--psf-relative-tail R] "
"[--psf-relative-tail R] [--psf-min-y Y] "
"[--psf-direct] [--verbose] "
#ifdef ENABLE_HDR_OUTPUT
"[--hdr-output] "
src/optics.c
+50 -11
View File
@@ -295,10 +295,14 @@ void psf_kernel_cache_report(const PsfKernelCache *cache,
if (stream == NULL)
return;
(void)cache;
fprintf(stream, "PSF splats: cached %zu, cached wing-clipped %zu, direct fallbacks %zu\n",
fprintf(stream, "PSF splats: cached %zu, cached wing-clipped %zu, direct fallbacks %zu, "
"discarded below min-Y %zu\n",
stats == NULL ? 0u : stats->cached_splats,
stats == NULL ? 0u : stats->cached_wing_clipped,
stats == NULL ? 0u : stats->direct_fallbacks);
stats == NULL ? 0u : stats->direct_fallbacks,
stats == NULL ? 0u : stats->discarded_below_min_y);
if (stats != NULL && stats->discarded_below_min_y != 0)
fputs("Warning: --psf-min-y discarded one or more PSF events.\n", stream);
#ifdef GR_DEBUG
if (stats != NULL)
fprintf(stream, "Debug: max raw magnification %.6g; magnification-clamped "
@@ -325,18 +329,49 @@ static int moffat_row_offset_range(double support_radius, double fx, double fy,
return *first_offset_x <= *last_offset_x;
}
static double linear_rgb_luminance(LinearRgb color)
{
return 0.2126729 * color.r + 0.7151522 * color.g + 0.0721750 * color.b;
}
/* `min_y` is a deliberate rendering cutoff. It compares against the
* continuous Moffat luminance profile rather than adding a second pixel-level
* cache lookup at the boundary. */
static double moffat_min_y_radius(double alpha, double beta, LinearRgb color,
double flux, double min_y)
{
if (min_y <= 0.0)
return INFINITY;
const double peak_y = linear_rgb_luminance(color) * flux *
(beta - 1.0) / (pi * alpha * alpha);
if (!(peak_y > min_y) || !isfinite(peak_y))
return 0.0;
return alpha * sqrt(pow(min_y / peak_y, -1.0 / beta) - 1.0);
}
static double moffat_effective_support_radius(double alpha, double beta,
LinearRgb color, double flux,
double relative_tail_fraction,
double min_y)
{
const double strict_radius = moffat_support_radius(
alpha, beta, flux, relative_tail_fraction);
const double min_y_radius = moffat_min_y_radius(alpha, beta, color, flux, min_y);
return fmin(strict_radius, min_y_radius);
}
void splat_moffat_direct(double *hdr, int width, int height, double x, double y,
LinearRgb color, double flux,
const PointSpreadFunction *psf,
double relative_tail_fraction)
double relative_tail_fraction, double min_y)
{
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, relative_tail_fraction);
const double support_radius = moffat_effective_support_radius(
alpha, psf->moffat_beta, color, flux, relative_tail_fraction, min_y);
if (!isfinite(support_radius))
return;
const int support = (int)ceil(support_radius);
@@ -371,20 +406,24 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
const PointSpreadFunction *psf,
const PsfKernelCache *cache,
double max_cache_psf_flux,
double relative_tail_fraction)
double relative_tail_fraction, double min_y)
{
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, relative_tail_fraction);
if (!isfinite(min_y) || min_y < 0.0)
return 1;
const double support_radius = moffat_effective_support_radius(
alpha, psf->moffat_beta, color, flux, relative_tail_fraction, min_y);
if (support_radius == 0.0)
return 3;
if (cache == NULL || !cache->ready || cache->fwhm_pixels != psf->fwhm_pixels ||
cache->moffat_beta != psf->moffat_beta ||
cache->relative_tail_fraction != relative_tail_fraction ||
!isfinite(support_radius) ||
(support_radius > cache->max_radius_pixels && flux > max_cache_psf_flux)) {
splat_moffat_direct(hdr, width, height, x, y, color, flux, psf,
relative_tail_fraction);
relative_tail_fraction, min_y);
return 1;
}
const double base_x = floor(x), base_y = floor(y);
@@ -427,10 +466,10 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
void splat_moffat(double *hdr, int width, int height, double x, double y,
LinearRgb color, double flux,
const PointSpreadFunction *psf,
double relative_tail_fraction)
double relative_tail_fraction, double min_y)
{
splat_moffat_direct(hdr, width, height, x, y, color, flux, psf,
relative_tail_fraction);
relative_tail_fraction, min_y);
}
static unsigned char tonemap_channel(double hdr_value)
src/optics.h
+4 -3
View File
@@ -25,6 +25,7 @@ typedef struct {
size_t cached_splats;
size_t cached_wing_clipped;
size_t direct_fallbacks;
size_t discarded_below_min_y;
#ifdef GR_DEBUG
double max_raw_magnification;
size_t magnification_clamped_triangles;
@@ -46,11 +47,11 @@ void psf_kernel_cache_report(const PsfKernelCache *cache,
void splat_moffat_direct(double *hdr, int width, int height, double x, double y,
LinearRgb color, double flux,
const PointSpreadFunction *psf,
double relative_tail_fraction);
double relative_tail_fraction, double min_y);
void splat_moffat(double *hdr, int width, int height, double x, double y,
LinearRgb color, double flux,
const PointSpreadFunction *psf,
double relative_tail_fraction);
double relative_tail_fraction, double min_y);
/* Returns 0 for a complete cached splat, 1 for the direct reference fallback,
* and 2 when the cached core was used with its outer wing intentionally
* clipped. max_cache_psf_flux == 1 preserves the historical behavior. */
@@ -59,7 +60,7 @@ int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
const PointSpreadFunction *psf,
const PsfKernelCache *cache,
double max_cache_psf_flux,
double relative_tail_fraction);
double relative_tail_fraction, double min_y);
/* 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);