668 lines
29 KiB
C
668 lines
29 KiB
C
#include "optics.h"
|
|
|
|
#include <limits.h>
|
|
#include <math.h>
|
|
#include <omp.h>
|
|
#include <stdint.h>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
#include <fitsio.h>
|
|
#endif
|
|
|
|
#ifdef ENABLE_PNG
|
|
#include <png.h>
|
|
#endif
|
|
|
|
static double clamp(double value, double low, double high)
|
|
{
|
|
return value < low ? low : value > high ? high : value;
|
|
}
|
|
|
|
/* Wyman, Sloan, and Shirley (2013), Eq. 4: analytic fits to the CIE 1931
|
|
* 2-degree color-matching functions. Wavelength is in nanometres. */
|
|
static void cie_1931_2deg(double wavelength_nm, double xyz[3])
|
|
{
|
|
const double x1 = (wavelength_nm - 442.0) *
|
|
(wavelength_nm < 442.0 ? 0.0624 : 0.0374);
|
|
const double x2 = (wavelength_nm - 599.8) *
|
|
(wavelength_nm < 599.8 ? 0.0264 : 0.0323);
|
|
const double x3 = (wavelength_nm - 501.1) *
|
|
(wavelength_nm < 501.1 ? 0.0490 : 0.0382);
|
|
const double y1 = (wavelength_nm - 568.8) *
|
|
(wavelength_nm < 568.8 ? 0.0213 : 0.0247);
|
|
const double y2 = (wavelength_nm - 530.9) *
|
|
(wavelength_nm < 530.9 ? 0.0613 : 0.0322);
|
|
const double z1 = (wavelength_nm - 437.0) *
|
|
(wavelength_nm < 437.0 ? 0.0845 : 0.0278);
|
|
const double z2 = (wavelength_nm - 459.0) *
|
|
(wavelength_nm < 459.0 ? 0.0385 : 0.0725);
|
|
xyz[0] = 0.362 * exp(-0.5 * x1 * x1) +
|
|
1.056 * exp(-0.5 * x2 * x2) - 0.065 * exp(-0.5 * x3 * x3);
|
|
xyz[1] = 0.821 * exp(-0.5 * y1 * y1) +
|
|
0.286 * exp(-0.5 * y2 * y2);
|
|
xyz[2] = 1.217 * exp(-0.5 * z1 * z1) +
|
|
0.681 * exp(-0.5 * z2 * z2);
|
|
}
|
|
|
|
static double planck_radiance_lambda(double wavelength_m, double temperature_K)
|
|
{
|
|
const double h = 6.62607015e-34;
|
|
const double c = 299792458.0;
|
|
const double k = 1.380649e-23;
|
|
const double exponent = h * c / (wavelength_m * k * temperature_K);
|
|
return 2.0 * h * c * c /
|
|
(pow(wavelength_m, 5.0) * expm1(exponent));
|
|
}
|
|
|
|
LinearRgb blackbody_to_linear_rgb(double temperature_K)
|
|
{
|
|
/* Integrate Planck spectral radiance from 380 to 780 nm into CIE XYZ,
|
|
* then transform XYZ to linear sRGB. Results are W m^-2 sr^-1 in each
|
|
* linear-primary channel, before the catalog amplitude and exposure. */
|
|
double xyz[3] = {0.0, 0.0, 0.0};
|
|
const double wavelength_step_m = 5e-9;
|
|
if (!isfinite(temperature_K) || temperature_K <= 0.0)
|
|
return (LinearRgb){0.0, 0.0, 0.0};
|
|
for (int wavelength_nm = 380; wavelength_nm <= 780; wavelength_nm += 5) {
|
|
double matching[3];
|
|
const double radiance =
|
|
planck_radiance_lambda(wavelength_nm * 1e-9, temperature_K);
|
|
cie_1931_2deg(wavelength_nm, matching);
|
|
for (int channel = 0; channel < 3; ++channel)
|
|
xyz[channel] += radiance * matching[channel] * wavelength_step_m;
|
|
}
|
|
return (LinearRgb){
|
|
clamp(3.24096994 * xyz[0] - 1.53738318 * xyz[1] - 0.49861076 * xyz[2],
|
|
0.0, INFINITY),
|
|
clamp(-0.96924364 * xyz[0] + 1.87596750 * xyz[1] + 0.04155506 * xyz[2],
|
|
0.0, INFINITY),
|
|
clamp(0.05563008 * xyz[0] - 0.20397696 * xyz[1] + 1.05697151 * xyz[2],
|
|
0.0, INFINITY)};
|
|
}
|
|
|
|
enum {
|
|
PSF_PHASE_RESOLUTION = 64,
|
|
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_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,
|
|
double relative_tail_fraction)
|
|
{
|
|
if (!isfinite(relative_tail_fraction) || relative_tail_fraction <= 0.0 ||
|
|
relative_tail_fraction >= 1.0)
|
|
return NAN;
|
|
const double normalization = scaled_flux * (beta - 1.0) / (pi * alpha * alpha);
|
|
const double relative_tail = fmin(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,
|
|
double relative_tail_fraction)
|
|
{
|
|
if (cache == NULL || !valid_psf(psf) ||
|
|
!isfinite(relative_tail_fraction) || relative_tail_fraction <= 0.0 ||
|
|
relative_tail_fraction >= 1.0)
|
|
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, relative_tail_fraction);
|
|
if (!isfinite(requested_radius) || requested_radius <= 0.0)
|
|
return -1;
|
|
/* The ordinary cache must cover the support implied by this invocation's
|
|
* PSF parameters. In particular, a low beta has broad Moffat wings even
|
|
* when its FWHM is small; a fixed radius would make every ordinary event
|
|
* fall back to the much slower direct evaluator. Brighter events can
|
|
* still exceed this flux=1 reference support and use that fallback. */
|
|
if (requested_radius > (double)INT_MAX - 1.0)
|
|
return -1;
|
|
const int radius = (int)ceil(requested_radius) + 1;
|
|
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 double start = omp_get_wtime();
|
|
PsfKernelCache building = {.fwhm_pixels = psf->fwhm_pixels,
|
|
.moffat_beta = psf->moffat_beta,
|
|
.alpha_pixels = alpha,
|
|
.relative_tail_fraction = relative_tail_fraction,
|
|
.max_radius_pixels = radius,
|
|
.weights = weights,
|
|
.phase_resolution = PSF_PHASE_RESOLUTION,
|
|
.radius_pixels = radius};
|
|
int invalid_kernel = 0;
|
|
/* Each sub-pixel phase owns a disjoint kernel plane. Parallelizing this
|
|
* regular construction needs no locks, and the finished kernel remains
|
|
* immutable to splat workers. */
|
|
#pragma omp parallel for collapse(2) schedule(static)
|
|
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)) {
|
|
/* An error cannot return from an OpenMP structured block. */
|
|
#pragma omp atomic write
|
|
invalid_kernel = 1;
|
|
continue;
|
|
}
|
|
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);
|
|
}
|
|
}
|
|
if (invalid_kernel) {
|
|
free(weights);
|
|
return -1;
|
|
}
|
|
building.build_seconds = omp_get_wtime() - start;
|
|
building.ready = 1;
|
|
*cache = building;
|
|
return 0;
|
|
}
|
|
|
|
void psf_kernel_cache_report_ready(const PsfKernelCache *cache, FILE *stream) {
|
|
if (stream == NULL)
|
|
return;
|
|
if (cache == NULL || !cache->ready) {
|
|
fputs("PSF cache: disabled; using the direct evaluator.\n", stream);
|
|
return;
|
|
}
|
|
fprintf(stream, "PSF cache ready: 64x64 phases, radius %.0f px, relative tail %.0e, "
|
|
"tail abs %.0e, boundary %.0e, build %.3f s\n",
|
|
cache->max_radius_pixels, cache->relative_tail_fraction,
|
|
psf_tail_absolute_hdr_budget,
|
|
psf_boundary_hdr_budget, cache->build_seconds);
|
|
}
|
|
|
|
void psf_kernel_cache_report(const PsfKernelCache *cache,
|
|
const PsfSplatStats *stats, FILE *stream)
|
|
{
|
|
if (stream == NULL)
|
|
return;
|
|
(void)cache;
|
|
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->discarded_below_min_y);
|
|
if (stats != NULL && stats->gpu_batch_count != 0)
|
|
fprintf(stream, "HIP PSF: %zu events in %zu batches (%zu timed); upload %.6f s, "
|
|
"kernel %.6f s, download %.6f s\n",
|
|
stats->gpu_event_count, stats->gpu_batch_count,
|
|
stats->gpu_timed_batch_count,
|
|
stats->gpu_upload_seconds, stats->gpu_kernel_seconds,
|
|
stats->gpu_download_seconds);
|
|
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 "
|
|
"triangles %zu\n",
|
|
stats->max_raw_magnification, stats->magnification_clamped_triangles);
|
|
#endif
|
|
}
|
|
|
|
/* Return the exact integer x-offset interval lying inside the circular support
|
|
* on one y-offset row. This avoids visiting the square box's corner pixels. */
|
|
static int moffat_row_offset_range(double support_radius, double fx, double fy,
|
|
int support, int offset_y,
|
|
int *first_offset_x, int *last_offset_x)
|
|
{
|
|
const double dy = offset_y + 0.5 - fy;
|
|
const double remaining = support_radius * support_radius - dy * dy;
|
|
if (remaining < 0.0)
|
|
return 0;
|
|
const double half_span = sqrt(remaining);
|
|
const int first = (int)ceil(fx - 0.5 - half_span);
|
|
const int last = (int)floor(fx - 0.5 + half_span);
|
|
*first_offset_x = first > -support ? first : -support;
|
|
*last_offset_x = last < support ? last : support;
|
|
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 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_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);
|
|
const int base_x = (int)floor(x), base_y = (int)floor(y);
|
|
const double fx = x - base_x, fy = y - base_y;
|
|
for (int offset_y = -support; offset_y <= support; ++offset_y) {
|
|
const int py = base_y + offset_y;
|
|
if (py < 0 || py >= height)
|
|
continue;
|
|
int first_offset_x, last_offset_x;
|
|
if (!moffat_row_offset_range(support_radius, fx, fy, support, offset_y,
|
|
&first_offset_x, &last_offset_x))
|
|
continue;
|
|
if (first_offset_x < -base_x)
|
|
first_offset_x = -base_x;
|
|
if (last_offset_x >= width - base_x)
|
|
last_offset_x = width - base_x - 1;
|
|
for (int offset_x = first_offset_x; offset_x <= last_offset_x; ++offset_x) {
|
|
const int px = base_x + offset_x;
|
|
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 psf_prepare_cached_event(PsfCachedEvent *event, double x, double y,
|
|
LinearRgb color, double flux,
|
|
const PointSpreadFunction *psf,
|
|
const PsfKernelCache *cache,
|
|
double max_cache_psf_flux,
|
|
double relative_tail_fraction, double min_y)
|
|
{
|
|
if (event == NULL || flux <= 0.0 || !valid_psf(psf))
|
|
return 1;
|
|
const double alpha = moffat_alpha(psf);
|
|
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)) {
|
|
return 1;
|
|
}
|
|
*event = (PsfCachedEvent){.x = x, .y = y, .color = color, .flux = flux,
|
|
.support_radius = support_radius};
|
|
return support_radius > cache->max_radius_pixels ? 2 : 0;
|
|
}
|
|
|
|
int splat_moffat_cached(double *hdr, int width, int height, double x, double y,
|
|
LinearRgb color, double flux,
|
|
const PointSpreadFunction *psf,
|
|
const PsfKernelCache *cache,
|
|
double max_cache_psf_flux,
|
|
double relative_tail_fraction, double min_y)
|
|
{
|
|
if (hdr == NULL || width <= 0 || height <= 0)
|
|
return 1;
|
|
PsfCachedEvent event;
|
|
const int status = psf_prepare_cached_event(
|
|
&event, x, y, color, flux, psf, cache, max_cache_psf_flux,
|
|
relative_tail_fraction, min_y);
|
|
if (status == 1) {
|
|
splat_moffat_direct(hdr, width, height, x, y, color, flux, psf,
|
|
relative_tail_fraction, min_y);
|
|
return status;
|
|
}
|
|
if (status == 3)
|
|
return status;
|
|
splat_prepared_cached_event(hdr, width, height, &event, cache);
|
|
return status;
|
|
}
|
|
|
|
void splat_prepared_cached_event(double *hdr, int width, int height,
|
|
const PsfCachedEvent *event,
|
|
const PsfKernelCache *cache)
|
|
{
|
|
if (hdr == NULL || width <= 0 || height <= 0 || event == NULL || cache == NULL ||
|
|
!cache->ready)
|
|
return;
|
|
const double base_x = floor(event->x), base_y = floor(event->y);
|
|
const double fx = event->x - base_x, fy = event->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)fmin(ceil(event->support_radius), cache->max_radius_pixels);
|
|
for (int offset_y = -support; offset_y <= support; ++offset_y) {
|
|
const int py = (int)base_y + offset_y;
|
|
if (py < 0 || py >= height)
|
|
continue;
|
|
int first_offset_x, last_offset_x;
|
|
if (!moffat_row_offset_range(event->support_radius, fx, fy, support, offset_y,
|
|
&first_offset_x, &last_offset_x))
|
|
continue;
|
|
if (first_offset_x < -(int)base_x)
|
|
first_offset_x = -(int)base_x;
|
|
if (last_offset_x >= width - (int)base_x)
|
|
last_offset_x = width - (int)base_x - 1;
|
|
for (int offset_x = first_offset_x; offset_x <= last_offset_x; ++offset_x) {
|
|
const int px = (int)base_x + offset_x;
|
|
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] += event->color.r * event->flux * weight;
|
|
pixel[1] += event->color.g * event->flux * weight;
|
|
pixel[2] += event->color.b * event->flux * weight;
|
|
}
|
|
}
|
|
}
|
|
|
|
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 min_y)
|
|
{
|
|
splat_moffat_direct(hdr, width, height, x, y, color, flux, psf,
|
|
relative_tail_fraction, min_y);
|
|
}
|
|
|
|
static unsigned char tonemap_channel(double hdr_value)
|
|
{
|
|
/* Reinhard tone mapping followed by the sRGB display transfer curve. */
|
|
const double linear = hdr_value / (1.0 + hdr_value);
|
|
const double display = linear <= 0.0031308 ? 12.92 * linear
|
|
: 1.055 * pow(linear, 1.0 / 2.4) - 0.055;
|
|
return (unsigned char)lround(255.0 * clamp(display, 0.0, 1.0));
|
|
}
|
|
|
|
#ifndef ENABLE_PNG
|
|
static int write_tonemapped_ppm(const char *path, const double *hdr, int width,
|
|
int height)
|
|
{
|
|
FILE *file = fopen(path, "wb");
|
|
if (file == NULL) return -1;
|
|
fprintf(file, "P6\n%d %d\n255\n", width, height);
|
|
for (int i = 0; i < width * height * 3; ++i) {
|
|
const unsigned char value = tonemap_channel(hdr[i]);
|
|
if (fwrite(&value, 1, 1, file) != 1) { fclose(file); return -1; }
|
|
}
|
|
return fclose(file) == 0 ? 0 : -1;
|
|
}
|
|
#endif
|
|
|
|
#ifdef ENABLE_PNG
|
|
static int write_tonemapped_png(const char *path, const double *hdr, int width,
|
|
int height)
|
|
{
|
|
FILE *file = fopen(path, "wb");
|
|
png_structp png = NULL;
|
|
png_infop info = NULL;
|
|
unsigned char *pixels = NULL;
|
|
int result = -1;
|
|
if (file == NULL) return -1;
|
|
png = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
|
|
if (png == NULL) goto done;
|
|
info = png_create_info_struct(png);
|
|
if (info == NULL || setjmp(png_jmpbuf(png))) goto done;
|
|
pixels = malloc((size_t)width * height * 3);
|
|
if (pixels == NULL) goto done;
|
|
for (int i = 0; i < width * height * 3; ++i)
|
|
pixels[i] = tonemap_channel(hdr[i]);
|
|
png_init_io(png, file);
|
|
png_set_IHDR(png, info, (png_uint_32)width, (png_uint_32)height, 8,
|
|
PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
|
|
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
|
|
png_write_info(png, info);
|
|
for (int row = 0; row < height; ++row)
|
|
png_write_row(png, &pixels[(size_t)row * width * 3]);
|
|
png_write_end(png, info);
|
|
result = 0;
|
|
done:
|
|
free(pixels);
|
|
png_destroy_write_struct(&png, &info);
|
|
if (fclose(file) != 0) result = -1;
|
|
return result;
|
|
}
|
|
#endif
|
|
|
|
int write_tonemapped_image(const char *path, const double *hdr, int width, int height)
|
|
{
|
|
const size_t path_length = strlen(path);
|
|
#ifdef ENABLE_PNG
|
|
if (path_length >= 4 && strcmp(path + path_length - 4, ".png") == 0)
|
|
return write_tonemapped_png(path, hdr, width, height);
|
|
fputs("PNG output is enabled; use a .png output path.\n", stderr);
|
|
return -1;
|
|
#else
|
|
if (path_length >= 4 && strcmp(path + path_length - 4, ".ppm") == 0)
|
|
return write_tonemapped_ppm(path, hdr, width, height);
|
|
fputs("PNG output is unavailable; use a .ppm output path or rebuild with libpng.\n",
|
|
stderr);
|
|
return -1;
|
|
#endif
|
|
}
|
|
|
|
#ifdef ENABLE_HDR_OUTPUT
|
|
int write_hdr_fits(const char *path, const double *hdr, int width, int height,
|
|
double horizontal_fov_deg)
|
|
{
|
|
if (path == NULL || hdr == NULL || width <= 0 || height <= 0 ||
|
|
!isfinite(horizontal_fov_deg) || horizontal_fov_deg <= 0.0 ||
|
|
horizontal_fov_deg >= 179.0)
|
|
return -1;
|
|
const size_t path_length = strlen(path);
|
|
char *overwrite_path = malloc(path_length + 2);
|
|
float *scanline = malloc((size_t)width * sizeof *scanline);
|
|
fitsfile *file = NULL;
|
|
int status = 0;
|
|
int result = -1;
|
|
long dimensions[3] = {width, height, 3};
|
|
char creator[] = "GR_4d_raytracing";
|
|
char unit[] = "linear HDR; arbitrary renderer scale";
|
|
char color_axis[] = "RGB";
|
|
char instrument[] = "GR4D virtual full-frame 50MP";
|
|
char detector_size[] = "[1:8640,1:5760]";
|
|
char crop_size[FLEN_VALUE];
|
|
double pixel_size_um = 36e3 / 8640.0;
|
|
const double crop_width_mm = width * pixel_size_um / 1000.0;
|
|
const double crop_height_mm = height * pixel_size_um / 1000.0;
|
|
double focal_length_mm =
|
|
crop_width_mm / (2.0 * tan(horizontal_fov_deg * pi / 360.0));
|
|
if (overwrite_path == NULL || scanline == NULL)
|
|
goto done;
|
|
snprintf(crop_size, sizeof crop_size, "%.6g x %.6g mm centered crop",
|
|
crop_width_mm, crop_height_mm);
|
|
/* CFITSIO uses a leading ! to give this output the same overwrite behavior
|
|
* as the former PFM writer. */
|
|
overwrite_path[0] = '!';
|
|
memcpy(overwrite_path + 1, path, path_length + 1);
|
|
fits_create_file(&file, overwrite_path, &status);
|
|
fits_create_img(file, FLOAT_IMG, 3, dimensions, &status);
|
|
fits_update_key(file, TSTRING, "CREATOR", creator, NULL, &status);
|
|
fits_update_key(file, TSTRING, "BUNIT", unit, NULL, &status);
|
|
fits_update_key(file, TSTRING, "CTYPE3", color_axis, NULL, &status);
|
|
fits_update_key(file, TSTRING, "INSTRUME", instrument,
|
|
"synthetic camera metadata", &status);
|
|
fits_update_key(file, TSTRING, "DETSIZE", detector_size,
|
|
"full-frame detector pixels", &status);
|
|
fits_update_key(file, TDOUBLE, "FOCALLEN", &focal_length_mm,
|
|
"mm; derived from horizontal field of view", &status);
|
|
fits_update_key(file, TDOUBLE, "XPIXSZ", &pixel_size_um,
|
|
"um; full-frame detector pixel size", &status);
|
|
fits_update_key(file, TDOUBLE, "YPIXSZ", &pixel_size_um,
|
|
"um; full-frame detector pixel size", &status);
|
|
fits_update_key(file, TSTRING, "CROPSIZE", crop_size,
|
|
"active sensor area for this render", &status);
|
|
/* FITS image coordinates start at the lower-left. The HDR buffer is
|
|
* top-down like PNG, so reverse rows while preserving its visual orientation. */
|
|
for (int channel = 0; status == 0 && channel < 3; ++channel)
|
|
for (int row = 0; status == 0 && row < height; ++row) {
|
|
const size_t offset = (size_t)row * width * 3;
|
|
for (int column = 0; column < width; ++column)
|
|
scanline[column] = (float)hdr[offset + 3 * column + channel];
|
|
long first_pixel[3] = {1, height - row, channel + 1};
|
|
fits_write_pix(file, TFLOAT, first_pixel, width, scanline, &status);
|
|
}
|
|
if (status == 0)
|
|
result = 0;
|
|
done:
|
|
if (file != NULL) {
|
|
int close_status = 0;
|
|
fits_close_file(file, &close_status);
|
|
if (status == 0 && close_status != 0)
|
|
status = close_status;
|
|
}
|
|
if (status != 0)
|
|
fits_report_error(stderr, status);
|
|
free(scanline);
|
|
free(overwrite_path);
|
|
return status == 0 ? result : -1;
|
|
}
|
|
#endif
|