Feat: Add CSV CIE 1931 blackbody LUT backend
Replace the in-tree Wyman analytic CIE fit with a repository CIE 1931 2-degree LUT derived from the 360-830 nm, 1 nm-linear CSV reference. The GRBBLUT3 table stores 1024 uniform log(T) XYZ nodes over [670.146556, 101408.88] K with four-point cubic Lagrange interpolation, and the same CSV reference supplies a three-term Rayleigh-Jeans form above the table and an inverse-temperature/log-XYZ crossover with channel-specific endpoint forms below it. The loader validates the header and payload checksum and rejects legacy formats and runtime generation. Initialize the backend in the renderer, frame test, and PSF capture tool so the new default is active wherever colors are produced.
This commit is contained in:
wyj committed
1 parent
2fc2b43e8d
commit
d2a901cab1
11 files changed
+866
-65
No files matched your search
@@ -0,0 +1,128 @@
|
||||
#include "blackbody_internal.h"
|
||||
|
||||
#include <math.h>
|
||||
|
||||
void blackbody_rayleigh_jeans_xyz(double temperature_K, int terms,
|
||||
double xyz[3]) {
|
||||
/* CSV-derived coefficients multiply T, 1, and T^-1. */
|
||||
static const double coefficients[3][3] = {
|
||||
{9.75239167492079072e3, -1.3438558214879372e8,
|
||||
6.31051886434512943e11},
|
||||
{9.51500076632225143e3, -1.25823411343233966e8,
|
||||
5.58059336995288075e11},
|
||||
{2.13884044303540753e4, -3.43402928472968325e8,
|
||||
1.84249991169126495e12}};
|
||||
xyz[0] = xyz[1] = xyz[2] = 0.0;
|
||||
if (!isfinite(temperature_K) || temperature_K <= 0.0 || terms <= 0)
|
||||
return;
|
||||
if (terms > 3) terms = 3;
|
||||
for (int channel = 0; channel < 3; ++channel) {
|
||||
xyz[channel] = coefficients[channel][0] * temperature_K;
|
||||
if (terms >= 2) xyz[channel] += coefficients[channel][1];
|
||||
if (terms >= 3) xyz[channel] += coefficients[channel][2] / temperature_K;
|
||||
}
|
||||
}
|
||||
|
||||
static double chebyshev_clenshaw(const double coefficients[9], double x) {
|
||||
double b1 = 0.0, b2 = 0.0;
|
||||
for (int degree = 8; degree >= 1; --degree) {
|
||||
const double b0 = 2.0 * x * b1 - b2 + coefficients[degree];
|
||||
b2 = b1;
|
||||
b1 = b0;
|
||||
}
|
||||
return coefficients[0] + x * b1 - b2;
|
||||
}
|
||||
|
||||
int blackbody_low_hybrid_log_xyz(double temperature_K, double log_xyz[3]) {
|
||||
static const double crossover_K = 4.09732109813541501;
|
||||
static const double maximum_K = 670.146556;
|
||||
static const double q_over_endpoint[3] = {
|
||||
17334.66117474619, 17334.66117474619, 22135.028884675903};
|
||||
static const double endpoint_delta[3] = {
|
||||
6.27453428003507278e-6, 6.27523604457500139e-6,
|
||||
9.72024116337713046e-5};
|
||||
/* The 1 nm-linear endpoint makes these polynomials exactly cubic. Z has
|
||||
* zero constant term because its 649--650 nm segment ends at zero. */
|
||||
static const double endpoint[3][4] = {
|
||||
{1.81136495194567529e-2, 6.61028930839653126e-5,
|
||||
1.48916636199141013e-8, 1.27817061855715528e-12},
|
||||
{6.54117161006302892e-3, 2.38707969128581198e-5,
|
||||
5.37761358258185981e-9, 4.61567454126977889e-13},
|
||||
{0.0, 1.18022769903543601e-3,
|
||||
2.13277824065096828e-7, 1.44529622149769957e-11}};
|
||||
static const double edges[7] = {
|
||||
4.09732109813541501, 9.58175285991856125, 22.4073207028704393,
|
||||
52.4004353296953668, 122.540559808645089, 286.566107776330601,
|
||||
670.146556};
|
||||
static const double coefficients[6][3][9] = {
|
||||
{
|
||||
{-2.19805342404174431,-0.427989370548201992,0.0457830231212280024,-0.00653024253368962181,0.00104792258851775496,-0.000179376495463268954,0.0000320194392828742251,-6.07285820395394865e-6,1.09768907792163728e-6},
|
||||
{-3.21660305217802182,-0.427989317567608284,0.0457830153418609421,-0.00653024167298543334,0.00104792254862916667,-0.000179376508222333492,0.0000320194445897039005,-6.07285961285998572e-6,1.09768936710106915e-6},
|
||||
{-3.15707121695679764,-0.842260091928597102,0.0900358213706993461,-0.0130901684838339702,0.00216325774923957118,-0.000378541811375747837,0.0000677130171872672182,-0.0000125764470117522098,2.20543446698214866e-6}
|
||||
},{
|
||||
{-1.31770949444839666,-0.441299064080369115,0.0486743388559153732,-0.00715733212557194619,0.00118377405373184402,-0.000208808548120484381,0.0000384098393642844184,-7.51938342275489625e-6,1.39706879054438949e-6},
|
||||
{-2.33625922712095605,-0.441299020893064677,0.0486743393223627113,-0.00715733306978699914,0.00118377427098721058,-0.000208808583754120261,0.0000384098442083096037,-7.51938398617148892e-6,1.39706884379154304e-6},
|
||||
{-1.41692623695155927,-0.878068580096058675,0.0982372095233778281,-0.0142171838603874451,0.00222410258776171213,-0.000361979109489253779,0.000060753034151493841,-0.0000107563963609356147,1.83055417864655942e-6}
|
||||
},{
|
||||
{-0.3921912858481657,-0.475886418801351195,0.0565482685596524812,-0.00894861708837536552,0.00159121925979185938,-0.000301479141093084493,0.0000595321929174356658,-0.0000125546605490984269,2.48503406235378882e-6},
|
||||
{-1.41074105985953395,-0.475886424111238694,0.0565482743959804094,-0.00894861785536701427,0.00159121932381919629,-0.000301479143653568406,0.0000595321922909604336,-0.0000125546602334654497,2.48503397406341601e-6},
|
||||
{0.405784497738078867,-0.915430672087484161,0.0971591394650644897,-0.0129050910771596898,0.00191184765152995853,-0.000308981746742567431,0.0000538107592092635484,-0.0000103283467419743846,1.93585445253475723e-6}
|
||||
},{
|
||||
{0.661994617272024548,-0.581350758628830698,0.0836597467165894418,-0.0158988753362009282,0.00336718167075529186,-0.000753662101271342622,0.000174585741711791468,-0.0000433048808785784171,9.75843572764750381e-6},
|
||||
{-0.356555122545496873,-0.581350782933150466,0.0836597486771623816,-0.015898875146527317,0.00336718160135198675,-0.000753662094742459443,0.000174585746021139991,-0.0000433048853568368614,9.75843776186480067e-6},
|
||||
{2.25724030901313405,-0.906149909784837897,0.0947967083167351412,-0.0147236528292987198,0.00294177345982300343,-0.000664178311762264752,0.000153812749558811817,-0.0000358188260289889306,7.26093914036807359e-6}
|
||||
},{
|
||||
{2.24517368371397051,-1.10137638803233863,0.260109466476170278,-0.0728062330999542422,0.0205487830596007146,-0.00547019170643925872,0.00127956756145661419,-0.000220004358600067039,0.0000138071701243527944},
|
||||
{1.22677135347875279,-1.10165295376793399,0.260338122908468991,-0.0729732322401341771,0.0206568819007716045,-0.00553248986531515299,0.00131198334055904886,-0.000236511389192968018,0.0000196226106014253749},
|
||||
{4.22881792094750399,-1.09273677796148737,0.167321326581346532,-0.0312233937525439919,0.00516074702621717587,-0.000696731180600597722,0.000080598581477044639,-0.000015356852062622951,5.09084046909121171e-6}
|
||||
},{
|
||||
{6.33932556236418281,-3.07589360106486531,0.456226630959950593,-0.00444896560449986224,-0.00803034195544535392,-0.000626248357061336987,0.000421665206117149112,0.0000469192966084690651,-0.0000234956663324828471},
|
||||
{5.3740651556541423,-3.15361482486290811,0.488370588919019434,-0.0119565427661776158,-0.0071239164324135692,-0.000725343801153405533,0.000461709147315171223,0.0000420002316640812232,-0.0000256745877581467219},
|
||||
{6.97759542644704118,-1.7029539661084809,0.314329784196094464,-0.0805402042898634484,0.0214126542337208624,-0.00404920923437737251,0.0000432174624337086157,0.000376147060473595442,-0.00015590604333367762}
|
||||
}};
|
||||
|
||||
if (!isfinite(temperature_K) || temperature_K <= 0.0 ||
|
||||
temperature_K > maximum_K)
|
||||
return -1;
|
||||
if (temperature_K <= crossover_K) {
|
||||
for (int channel = 0; channel < 3; ++channel) {
|
||||
double polynomial = endpoint[channel][3];
|
||||
for (int term = 2; term >= 0; --term)
|
||||
polynomial = endpoint[channel][term] + temperature_K * polynomial;
|
||||
log_xyz[channel] = -q_over_endpoint[channel] / temperature_K +
|
||||
log(temperature_K * polynomial) +
|
||||
endpoint_delta[channel] * temperature_K / crossover_K;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
int segment = 0;
|
||||
while (segment < 5 && temperature_K >= edges[segment + 1]) ++segment;
|
||||
const double inverse_center =
|
||||
0.5 * (1.0 / edges[segment] + 1.0 / edges[segment + 1]);
|
||||
const double inverse_half =
|
||||
0.5 * (1.0 / edges[segment] - 1.0 / edges[segment + 1]);
|
||||
const double x = (1.0 / temperature_K - inverse_center) / inverse_half;
|
||||
for (int channel = 0; channel < 3; ++channel)
|
||||
log_xyz[channel] = -q_over_endpoint[channel] / temperature_K +
|
||||
chebyshev_clenshaw(coefficients[segment][channel], x);
|
||||
return 0;
|
||||
}
|
||||
|
||||
int blackbody_low_hybrid_xyz(double temperature_K, double xyz[3]) {
|
||||
double log_xyz[3];
|
||||
if (blackbody_low_hybrid_log_xyz(temperature_K, log_xyz)) {
|
||||
xyz[0] = xyz[1] = xyz[2] = 0.0;
|
||||
return -1;
|
||||
}
|
||||
for (int channel = 0; channel < 3; ++channel)
|
||||
xyz[channel] = exp(log_xyz[channel]);
|
||||
return 0;
|
||||
}
|
||||
|
||||
LinearRgb blackbody_xyz_to_linear_rgb(const double xyz[3]) {
|
||||
return (LinearRgb){
|
||||
fmax(0.0, 3.24096994 * xyz[0] - 1.53738318 * xyz[1] - 0.49861076 * xyz[2]),
|
||||
fmax(0.0, -0.96924364 * xyz[0] + 1.87596750 * xyz[1] + 0.04155506 * xyz[2]),
|
||||
fmax(0.0, 0.05563008 * xyz[0] - 0.20397696 * xyz[1] + 1.05697151 * xyz[2])};
|
||||
}
|
||||
@@ -0,0 +1,12 @@
|
||||
#ifndef BLACKBODY_INTERNAL_H
|
||||
#define BLACKBODY_INTERNAL_H
|
||||
|
||||
#include "optics.h"
|
||||
|
||||
void blackbody_rayleigh_jeans_xyz(double temperature_K, int terms,
|
||||
double xyz[3]);
|
||||
int blackbody_low_hybrid_log_xyz(double temperature_K, double log_xyz[3]);
|
||||
int blackbody_low_hybrid_xyz(double temperature_K, double xyz[3]);
|
||||
LinearRgb blackbody_xyz_to_linear_rgb(const double xyz[3]);
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,175 @@
|
||||
#include "blackbody_internal.h"
|
||||
|
||||
#include <math.h>
|
||||
#include <stdint.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <time.h>
|
||||
|
||||
enum { BLACKBODY_LUT_VERSION_CIE1931 = 3,
|
||||
BLACKBODY_LUT_ELEMENT_F64 = 1,
|
||||
BLACKBODY_LUT_REFERENCE_CIE1931_1NM_LINEAR = 3 };
|
||||
|
||||
static const char default_table_path[] =
|
||||
"assets/blackbody/cie1931_2deg_xyz_1024.grbblut";
|
||||
|
||||
typedef struct {
|
||||
char magic[8];
|
||||
uint32_t version, header_bytes, endian_marker, element_type;
|
||||
uint32_t reference_type, wavelength_min_nm, wavelength_max_nm,
|
||||
wavelength_step_nm;
|
||||
uint64_t node_count;
|
||||
double log_temperature_min, log_temperature_max;
|
||||
double planck_h, planck_c, planck_k;
|
||||
char generator[64];
|
||||
uint64_t payload_checksum;
|
||||
} BlackbodyLutHeader;
|
||||
|
||||
_Static_assert(sizeof(BlackbodyLutHeader) == 160,
|
||||
"blackbody LUT header layout changed");
|
||||
|
||||
typedef struct {
|
||||
double *xyz;
|
||||
size_t nodes;
|
||||
double log_min, log_max, index_scale;
|
||||
int ready;
|
||||
} BlackbodyLut;
|
||||
|
||||
static BlackbodyLut table;
|
||||
|
||||
static uint64_t fnv1a64(const void *data, size_t bytes) {
|
||||
const unsigned char *p = data;
|
||||
uint64_t hash = UINT64_C(14695981039346656037);
|
||||
for (size_t i = 0; i < bytes; ++i) {
|
||||
hash ^= p[i];
|
||||
hash *= UINT64_C(1099511628211);
|
||||
}
|
||||
return hash;
|
||||
}
|
||||
|
||||
static int validate_header(const BlackbodyLutHeader *h) {
|
||||
return memcmp(h->magic, "GRBBLUT3", 8) == 0 &&
|
||||
h->version == BLACKBODY_LUT_VERSION_CIE1931 &&
|
||||
h->header_bytes == sizeof *h &&
|
||||
h->endian_marker == UINT32_C(0x01020304) &&
|
||||
h->element_type == BLACKBODY_LUT_ELEMENT_F64 &&
|
||||
h->reference_type == BLACKBODY_LUT_REFERENCE_CIE1931_1NM_LINEAR &&
|
||||
h->wavelength_min_nm == 360 && h->wavelength_max_nm == 830 &&
|
||||
h->wavelength_step_nm == 1 && h->node_count >= 4 &&
|
||||
h->node_count <= SIZE_MAX / (3 * sizeof(double)) &&
|
||||
isfinite(h->log_temperature_min) &&
|
||||
isfinite(h->log_temperature_max) &&
|
||||
h->log_temperature_max > h->log_temperature_min &&
|
||||
h->planck_h == 6.62607015e-34 && h->planck_c == 299792458.0 &&
|
||||
h->planck_k == 1.380649e-23;
|
||||
}
|
||||
|
||||
static int read_table(const char *path) {
|
||||
FILE *file = fopen(path, "rb");
|
||||
BlackbodyLutHeader h;
|
||||
if (file == NULL || fread(&h, sizeof h, 1, file) != 1 ||
|
||||
!validate_header(&h)) {
|
||||
if (file != NULL) fclose(file);
|
||||
return -1;
|
||||
}
|
||||
double *xyz = malloc(3 * (size_t)h.node_count * sizeof *xyz);
|
||||
const size_t count = 3 * (size_t)h.node_count;
|
||||
unsigned char extra;
|
||||
if (xyz == NULL || fread(xyz, sizeof *xyz, count, file) != count ||
|
||||
fread(&extra, 1, 1, file) != 0 || ferror(file) ||
|
||||
fnv1a64(xyz, count * sizeof *xyz) != h.payload_checksum) {
|
||||
free(xyz); fclose(file); return -1;
|
||||
}
|
||||
fclose(file);
|
||||
table = (BlackbodyLut){.xyz = xyz, .nodes = (size_t)h.node_count,
|
||||
.log_min = h.log_temperature_min,
|
||||
.log_max = h.log_temperature_max, .ready = 1};
|
||||
table.index_scale = (table.nodes - 1) / (table.log_max - table.log_min);
|
||||
return 0;
|
||||
}
|
||||
|
||||
int blackbody_backend_init(const char *table_path, size_t lut_nodes,
|
||||
double temperature_min_K,
|
||||
double temperature_max_K,
|
||||
const char *write_table_path, FILE *report) {
|
||||
struct timespec start, finish;
|
||||
timespec_get(&start, TIME_UTC);
|
||||
if (table.ready) return -1;
|
||||
if (lut_nodes != 0 || isfinite(temperature_min_K) ||
|
||||
isfinite(temperature_max_K) || write_table_path != NULL)
|
||||
return -1;
|
||||
const char *path = table_path != NULL ? table_path : default_table_path;
|
||||
if (read_table(path)) return -1;
|
||||
timespec_get(&finish, TIME_UTC);
|
||||
if (report != NULL)
|
||||
fprintf(report,
|
||||
"Blackbody LUT: %zu CIE 1931 2-deg 1 nm-linear XYZ nodes, "
|
||||
"T=[%.9g, %.9g] K, loaded %s in %.6f s; "
|
||||
"payload fnv1a64=%016llx\n",
|
||||
table.nodes,
|
||||
exp(table.log_min), exp(table.log_max),
|
||||
path,
|
||||
finish.tv_sec - start.tv_sec +
|
||||
1e-9 * (finish.tv_nsec - start.tv_nsec),
|
||||
(unsigned long long)fnv1a64(table.xyz,
|
||||
3 * table.nodes * sizeof *table.xyz));
|
||||
return 0;
|
||||
}
|
||||
|
||||
void blackbody_backend_destroy(void) {
|
||||
free(table.xyz);
|
||||
table = (BlackbodyLut){0};
|
||||
}
|
||||
|
||||
const char *blackbody_backend_name(void) { return "lut"; }
|
||||
|
||||
LinearRgb blackbody_to_linear_rgb(double temperature_K) {
|
||||
if (!table.ready || !isfinite(temperature_K) || temperature_K <= 0.0) {
|
||||
return (LinearRgb){0};
|
||||
}
|
||||
const double log_t = log(temperature_K);
|
||||
if (log_t < table.log_min) {
|
||||
double xyz[3];
|
||||
if (blackbody_low_hybrid_xyz(temperature_K, xyz))
|
||||
return (LinearRgb){0};
|
||||
return blackbody_xyz_to_linear_rgb(xyz);
|
||||
}
|
||||
if (log_t > table.log_max) {
|
||||
double xyz[3];
|
||||
blackbody_rayleigh_jeans_xyz(temperature_K, 3, xyz);
|
||||
return blackbody_xyz_to_linear_rgb(xyz);
|
||||
}
|
||||
const double position = (log_t - table.log_min) * table.index_scale;
|
||||
size_t low = (size_t)position;
|
||||
if (low >= table.nodes - 1) low = table.nodes - 2;
|
||||
const double fraction = position - low;
|
||||
double xyz[3];
|
||||
{
|
||||
size_t base;
|
||||
double x;
|
||||
if (low == 0) {
|
||||
base = 0;
|
||||
x = fraction;
|
||||
} else if (low >= table.nodes - 2) {
|
||||
base = table.nodes - 4;
|
||||
x = 2.0 + fraction;
|
||||
} else {
|
||||
base = low - 1;
|
||||
x = 1.0 + fraction;
|
||||
}
|
||||
const double weights[4] = {
|
||||
-(x - 1.0) * (x - 2.0) * (x - 3.0) / 6.0,
|
||||
x * (x - 2.0) * (x - 3.0) / 2.0,
|
||||
-x * (x - 1.0) * (x - 3.0) / 2.0,
|
||||
x * (x - 1.0) * (x - 2.0) / 6.0};
|
||||
for (int channel = 0; channel < 3; ++channel) {
|
||||
xyz[channel] = 0.0;
|
||||
for (int point = 0; point < 4; ++point)
|
||||
xyz[channel] += weights[point] *
|
||||
table.xyz[3 * (base + (size_t)point) + channel];
|
||||
xyz[channel] = fmax(0.0, xyz[channel]);
|
||||
}
|
||||
}
|
||||
return blackbody_xyz_to_linear_rgb(xyz);
|
||||
}
|
||||
+23
@@ -12,6 +12,7 @@
|
||||
#include <math.h>
|
||||
#include <omp.h>
|
||||
#include <stdio.h>
|
||||
#include <stdint.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <sys/stat.h>
|
||||
@@ -52,6 +53,7 @@ typedef struct {
|
||||
double slab_duration;
|
||||
double minkowski_proper_acceleration;
|
||||
int catalog_load_workers;
|
||||
const char *blackbody_table_path;
|
||||
RefinementConfig refinement;
|
||||
} Settings;
|
||||
|
||||
@@ -319,6 +321,8 @@ static int parse_args(int argc, char **argv, Settings *s,
|
||||
&s->minkowski_proper_acceleration)) {
|
||||
} else if (!strcmp(argv[i], "--catalog-load-workers") && i + 1 < argc &&
|
||||
!parse_int(argv[++i], &s->catalog_load_workers)) {
|
||||
} else if (!strcmp(argv[i], "--blackbody-table") && i + 1 < argc) {
|
||||
s->blackbody_table_path = argv[++i];
|
||||
} else if (!strcmp(argv[i], "--write-minkowski-accel-track") &&
|
||||
i + 1 < argc)
|
||||
s->write_minkowski_accel_track_path = argv[++i];
|
||||
@@ -377,6 +381,10 @@ static void print_help(const char *program) {
|
||||
" --psf-min-y Y Skip events below this linear HDR luminance (default: 0, disabled)\n"
|
||||
" --psf-direct Disable the PSF lookup cache (default: disabled)\n"
|
||||
" --catalog-load-workers N All-sky catalog loader workers (default: 4)\n"
|
||||
" --blackbody-table FILE Explicit GRBBLUT3 table\n"
|
||||
" (default: assets/blackbody/cie1931_2deg_xyz_1024.grbblut)\n",
|
||||
stdout);
|
||||
fputs(
|
||||
"\nAdaptive lens mesh:\n"
|
||||
" --coarse-cell-pixels N Initial mesh cell size in pixels (default: 16)\n"
|
||||
" --refine-max-level N Maximum refinement level (default: 0)\n"
|
||||
@@ -1111,6 +1119,18 @@ int main(int argc, char **argv) {
|
||||
}
|
||||
fprintf(stderr, "Created test catalog: %s\n", settings.catalog_path);
|
||||
}
|
||||
if (blackbody_backend_init(settings.blackbody_table_path,
|
||||
0, NAN, NAN, NULL, stderr)) {
|
||||
fprintf(stderr,
|
||||
"Blackbody backend '%s' initialization failed; the LUT loads the "
|
||||
"repository CIE GRBBLUT3 table by default or a valid explicit "
|
||||
"table.\n",
|
||||
blackbody_backend_name());
|
||||
catalog_destroy(&catalog);
|
||||
spacetime_destroy(&spacetime);
|
||||
return 1;
|
||||
}
|
||||
fprintf(stderr, "Blackbody backend: %s\n", blackbody_backend_name());
|
||||
if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf,
|
||||
settings.psf_relative_tail)) {
|
||||
#ifdef PSF_BACKEND_DUMMY
|
||||
@@ -1118,6 +1138,7 @@ int main(int argc, char **argv) {
|
||||
stderr);
|
||||
catalog_destroy(&catalog);
|
||||
spacetime_destroy(&spacetime);
|
||||
blackbody_backend_destroy();
|
||||
return 1;
|
||||
#else
|
||||
fputs("PSF cache construction failed; using direct evaluator.\n", stderr);
|
||||
@@ -1128,6 +1149,7 @@ int main(int argc, char **argv) {
|
||||
const int result = render_lens_map(&settings, &catalog);
|
||||
catalog_destroy(&catalog);
|
||||
psf_kernel_cache_destroy(&settings.psf_cache);
|
||||
blackbody_backend_destroy();
|
||||
return result == 0 ? 0 : 1;
|
||||
}
|
||||
int result = settings.frames_dir != NULL
|
||||
@@ -1137,5 +1159,6 @@ int main(int argc, char **argv) {
|
||||
spacetime_destroy(&spacetime);
|
||||
catalog_destroy(&catalog);
|
||||
psf_kernel_cache_destroy(&settings.psf_cache);
|
||||
blackbody_backend_destroy();
|
||||
return result == 0 ? 0 : 1;
|
||||
}
|
||||
@@ -21,68 +21,6 @@ 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,
|
||||
|
||||
@@ -50,6 +50,16 @@ extern "C" {
|
||||
/* 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);
|
||||
/* Backend initialization happens before renderer workers start. Integral and
|
||||
* fixed backends reject LUT-only options. The LUT backend loads the sole
|
||||
* repository CIE table by default, or validates an explicit equivalent file;
|
||||
* legacy runtime-generation arguments are rejected. */
|
||||
int blackbody_backend_init(const char *table_path, size_t lut_nodes,
|
||||
double temperature_min_K,
|
||||
double temperature_max_K,
|
||||
const char *write_table_path, FILE *report);
|
||||
void blackbody_backend_destroy(void);
|
||||
const char *blackbody_backend_name(void);
|
||||
int psf_kernel_cache_init(PsfKernelCache *cache,
|
||||
const PointSpreadFunction *psf,
|
||||
double relative_tail_fraction);
|
||||
|
||||
Reference in new issue
Block a user