htrdr_ran_wlen_cie_xyz.c (13385B)
1 /* Copyright (C) 2018-2019, 2022-2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2020-2022 Institut Mines Télécom Albi-Carmaux 3 * Copyright (C) 2022-2025 Institut Pierre-Simon Laplace 4 * Copyright (C) 2022-2025 Institut de Physique du Globe de Paris 5 * Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com) 6 * Copyright (C) 2022-2025 Observatoire de Paris 7 * Copyright (C) 2022-2025 Université de Reims Champagne-Ardenne 8 * Copyright (C) 2022-2025 Université de Versaille Saint-Quentin 9 * Copyright (C) 2018-2019, 2022-2025 Université Paul Sabatier 10 * 11 * This program is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU General Public License as published by 13 * the Free Software Foundation, either version 3 of the License, or 14 * (at your option) any later version. 15 * 16 * This program is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 * GNU General Public License for more details. 20 * 21 * You should have received a copy of the GNU General Public License 22 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 23 24 #define _POSIX_C_SOURCE 200112L /* nextafter */ 25 26 #include "core/htrdr.h" 27 #include "core/htrdr_c.h" 28 #include "core/htrdr_ran_wlen_cie_xyz.h" 29 #include "core/htrdr_log.h" 30 31 #include <rsys/algorithm.h> 32 #include <rsys/cstr.h> 33 #include <rsys/dynamic_array_double.h> 34 #include <rsys/mem_allocator.h> 35 #include <rsys/ref_count.h> 36 37 #include <math.h> /* nextafter */ 38 39 struct htrdr_ran_wlen_cie_xyz { 40 struct darray_double cdf_X; 41 struct darray_double cdf_Y; 42 struct darray_double cdf_Z; 43 double rcp_integral_X; 44 double rcp_integral_Y; 45 double rcp_integral_Z; 46 double range[2]; /* Boundaries of the handled CIE XYZ color space */ 47 double band_len; /* Length in nanometers of a band */ 48 49 struct htrdr* htrdr; 50 ref_T ref; 51 }; 52 53 /******************************************************************************* 54 * Helper functions 55 ******************************************************************************/ 56 static INLINE double 57 trapezoidal_integration 58 (const double lambda_lo, /* Integral lower bound. In nanometer */ 59 const double lambda_hi, /* Integral upper bound. In nanometer */ 60 double (*f_bar)(const double lambda)) /* Function to integrate */ 61 { 62 double dlambda; 63 size_t i, n; 64 double integral = 0; 65 ASSERT(lambda_lo <= lambda_hi); 66 ASSERT(lambda_lo > 0); 67 68 n = (size_t)(lambda_hi - lambda_lo) + 1; 69 dlambda = (lambda_hi - lambda_lo) / (double)n; 70 71 FOR_EACH(i, 0, n) { 72 const double lambda1 = lambda_lo + dlambda*(double)(i+0); 73 const double lambda2 = lambda_lo + dlambda*(double)(i+1); 74 const double f1 = f_bar(lambda1); 75 const double f2 = f_bar(lambda2); 76 integral += (f1 + f2)*dlambda*0.5; 77 } 78 return integral; 79 } 80 81 /* The following 3 functions are used to fit the CIE Xbar, Ybar and Zbar curved 82 * has defined by the 1931 standard. These analytical fits are propsed by C. 83 * Wyman, P. P. Sloan & P. Shirley in "Simple Analytic Approximations to the 84 * CIE XYZ Color Matching Functions" - JCGT 2013. */ 85 static INLINE double 86 fit_x_bar_1931(const double lambda) 87 { 88 const double a = (lambda - 442.0) * (lambda < 442.0 ? 0.0624 : 0.0374); 89 const double b = (lambda - 599.8) * (lambda < 599.8 ? 0.0264 : 0.0323); 90 const double c = (lambda - 501.1) * (lambda < 501.1 ? 0.0490 : 0.0382); 91 return 0.362*exp(-0.5*a*a) + 1.056*exp(-0.5f*b*b) - 0.065*exp(-0.5*c*c); 92 } 93 94 static FINLINE double 95 fit_y_bar_1931(const double lambda) 96 { 97 const double a = (lambda - 568.8) * (lambda < 568.8 ? 0.0213 : 0.0247); 98 const double b = (lambda - 530.9) * (lambda < 530.9 ? 0.0613 : 0.0322); 99 return 0.821*exp(-0.5*a*a) + 0.286*exp(-0.5*b*b); 100 } 101 102 static FINLINE double 103 fit_z_bar_1931(const double lambda) 104 { 105 const double a = (lambda - 437.0) * (lambda < 437.0 ? 0.0845 : 0.0278); 106 const double b = (lambda - 459.0) * (lambda < 459.0 ? 0.0385 : 0.0725); 107 return 1.217*exp(-0.5*a*a) + 0.681*exp(-0.5*b*b); 108 } 109 110 static INLINE double 111 sample_cie_xyz 112 (const struct htrdr_ran_wlen_cie_xyz* cie, 113 const double* cdf, 114 const size_t cdf_length, 115 double (*f_bar)(const double lambda), /* Function to integrate */ 116 const double r0, /* Canonical number in [0, 1[ */ 117 const double r1) /* Canonical number in [0, 1[ */ 118 { 119 double r0_next = nextafter(r0, DBL_MAX); 120 double* find; 121 double f_min, f_max; /* CIE 1931 value for the band boundaries */ 122 double lambda; /* Sampled wavelength */ 123 double lambda_min, lambda_max; /* Boundaries of the sampled band */ 124 double lambda_1, lambda_2; /* Solutions if the equation to solve */ 125 double a, b, c, d; /* Equation parameters */ 126 double delta, sqrt_delta; 127 size_t iband; /* Index of the sampled band */ 128 ASSERT(cie && cdf && cdf_length); 129 ASSERT(0 <= r0 && r0 < 1); 130 ASSERT(0 <= r1 && r1 < 1); 131 132 /* Use r_next rather than r in order to find the first entry that is not less 133 * than *or equal* to r */ 134 find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl); 135 ASSERT(find); 136 137 /* Define and check the sampled band */ 138 iband = (size_t)(find - cdf); 139 ASSERT(iband < cdf_length); 140 ASSERT(cdf[iband] > r0 && (!iband || cdf[iband-1] <= r0)); 141 142 /* Define the boundaries of the sampled band */ 143 lambda_min = cie->range[0] + cie->band_len * (double)iband; 144 lambda_max = lambda_min + cie->band_len; 145 146 /* Define the value of the CIE 1931 function for the boudaries of the sampled 147 * band */ 148 f_min = f_bar(lambda_min); 149 f_max = f_bar(lambda_max); 150 151 /* Compute the equation constants */ 152 a = 0.5 * (f_max - f_min) / cie->band_len; 153 b = (lambda_max * f_min - lambda_min * f_max) / cie->band_len; 154 c = -lambda_min * f_min + lambda_min*lambda_min * a; 155 d = 0.5 * (f_max + f_min) * cie->band_len; 156 157 delta = b*b - 4*a*(c-d*r1); 158 if(delta < 0 && eq_eps(delta, 0, 1.e-6)) { 159 delta = 0; 160 } 161 ASSERT(delta > 0); 162 sqrt_delta = sqrt(delta); 163 164 /* Compute the roots that solve the equation */ 165 lambda_1 = (-b - sqrt_delta) / (2*a); 166 lambda_2 = (-b + sqrt_delta) / (2*a); 167 168 /* Select the solution */ 169 if(lambda_min <= lambda_1 && lambda_1 < lambda_max) { 170 lambda = lambda_1; 171 } else if(lambda_min <= lambda_2 && lambda_2 < lambda_max) { 172 lambda = lambda_2; 173 } else { 174 htrdr_log_warn(cie->htrdr, 175 "%s: cannot sample a wavelength in [%g, %g[. The possible wavelengths" 176 "were %g and %g.\n", 177 FUNC_NAME, lambda_min, lambda_max, lambda_1, lambda_2); 178 /* Arbitrarly choose the wavelength at the center of the sampled band */ 179 lambda = (lambda_min + lambda_max)*0.5; 180 } 181 182 return lambda; 183 } 184 185 static res_T 186 setup_cie_xyz 187 (struct htrdr_ran_wlen_cie_xyz* cie, 188 const char* func_name, 189 const size_t nbands) 190 { 191 enum { X, Y, Z }; /* Helper constant */ 192 double* pdf[3] = {NULL, NULL, NULL}; 193 double* cdf[3] = {NULL, NULL, NULL}; 194 double sum[3] = {0, 0, 0}; 195 size_t i; 196 res_T res = RES_OK; 197 198 ASSERT(cie && func_name && nbands); 199 ASSERT(cie->range[0] >= HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[0]); 200 ASSERT(cie->range[1] <= HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[1]); 201 ASSERT(cie->range[0] < cie->range[1]); 202 203 /* Allocate and reset the memory space for the tristimulus CDF */ 204 #define SETUP_STIMULUS(Stimulus) { \ 205 res = darray_double_resize(&cie->cdf_ ## Stimulus, nbands); \ 206 if(res != RES_OK) { \ 207 htrdr_log_err(cie->htrdr, \ 208 "%s: Could not reserve the memory space for the CDF " \ 209 "of the "STR(X)" stimulus -- %s.\n", func_name, res_to_cstr(res)); \ 210 goto error; \ 211 } \ 212 cdf[Stimulus] = darray_double_data_get(&cie->cdf_ ## Stimulus); \ 213 pdf[Stimulus] = cdf[Stimulus]; \ 214 memset(cdf[Stimulus], 0, nbands*sizeof(double)); \ 215 } (void)0 216 SETUP_STIMULUS(X); 217 SETUP_STIMULUS(Y); 218 SETUP_STIMULUS(Z); 219 #undef SETUP_STIMULUS 220 221 /* Compute the *unormalized* pdf of the tristimulus */ 222 FOR_EACH(i, 0, nbands) { 223 const double lambda_lo = cie->range[0] + (double)i * cie->band_len; 224 const double lambda_hi = MMIN(lambda_lo + cie->band_len, cie->range[1]); 225 ASSERT(lambda_lo <= lambda_hi); 226 ASSERT(lambda_lo >= cie->range[0]); 227 ASSERT(lambda_hi <= cie->range[1]); 228 pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931); 229 pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931); 230 pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931); 231 sum[X] += pdf[X][i]; 232 sum[Y] += pdf[Y][i]; 233 sum[Z] += pdf[Z][i]; 234 } 235 #define CHK_SUM(Sum, Range, Fit) \ 236 ASSERT(eq_eps(Sum, trapezoidal_integration(Range[0], Range[1], Fit), 1.e-3)) 237 CHK_SUM(sum[X], cie->range, fit_x_bar_1931); 238 CHK_SUM(sum[Y], cie->range, fit_y_bar_1931); 239 CHK_SUM(sum[Z], cie->range, fit_z_bar_1931); 240 #undef CHK_SUM 241 cie->rcp_integral_X = 1.0 / sum[X]; 242 cie->rcp_integral_Y = 1.0 / sum[Y]; 243 cie->rcp_integral_Z = 1.0 / sum[Z]; 244 245 FOR_EACH(i, 0, nbands) { 246 /* Normalize the pdf */ 247 pdf[X][i] /= sum[X]; 248 pdf[Y][i] /= sum[Y]; 249 pdf[Z][i] /= sum[Z]; 250 /* Setup the cumulative */ 251 if(i == 0) { 252 cdf[X][i] = pdf[X][i]; 253 cdf[Y][i] = pdf[Y][i]; 254 cdf[Z][i] = pdf[Z][i]; 255 } else { 256 cdf[X][i] = pdf[X][i] + cdf[X][i-1]; 257 cdf[Y][i] = pdf[Y][i] + cdf[Y][i-1]; 258 cdf[Z][i] = pdf[Z][i] + cdf[Z][i-1]; 259 ASSERT(cdf[X][i] >= cdf[X][i-1]); 260 ASSERT(cdf[Y][i] >= cdf[Y][i-1]); 261 ASSERT(cdf[Z][i] >= cdf[Z][i-1]); 262 } 263 } 264 ASSERT(eq_eps(cdf[X][nbands-1], 1, 1.e-6)); 265 ASSERT(eq_eps(cdf[Y][nbands-1], 1, 1.e-6)); 266 ASSERT(eq_eps(cdf[Z][nbands-1], 1, 1.e-6)); 267 268 /* Handle numerical issue */ 269 cdf[X][nbands-1] = 1.0; 270 cdf[Y][nbands-1] = 1.0; 271 cdf[Z][nbands-1] = 1.0; 272 273 exit: 274 return res; 275 error: 276 darray_double_clear(&cie->cdf_X); 277 darray_double_clear(&cie->cdf_Y); 278 darray_double_clear(&cie->cdf_Z); 279 goto exit; 280 } 281 282 static void 283 release_cie_xyz(ref_T* ref) 284 { 285 struct htrdr_ran_wlen_cie_xyz* cie = NULL; 286 struct htrdr* htrdr = NULL; 287 ASSERT(ref); 288 cie = CONTAINER_OF(ref, struct htrdr_ran_wlen_cie_xyz, ref); 289 darray_double_release(&cie->cdf_X); 290 darray_double_release(&cie->cdf_Y); 291 darray_double_release(&cie->cdf_Z); 292 htrdr = cie->htrdr; 293 MEM_RM(htrdr_get_allocator(cie->htrdr), cie); 294 htrdr_ref_put(htrdr); 295 } 296 297 /******************************************************************************* 298 * Local functions 299 ******************************************************************************/ 300 res_T 301 htrdr_ran_wlen_cie_xyz_create 302 (struct htrdr* htrdr, 303 const double range[2], /* Must be included in [380, 780] nanometers */ 304 const size_t bands_count, /* # bands used to discretize the CIE tristimulus */ 305 struct htrdr_ran_wlen_cie_xyz** out_cie) 306 { 307 struct htrdr_ran_wlen_cie_xyz* cie = NULL; 308 size_t nbands = bands_count; 309 res_T res = RES_OK; 310 ASSERT(htrdr && range && nbands && out_cie); 311 312 cie = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cie)); 313 if(!cie) { 314 res = RES_MEM_ERR; 315 htrdr_log_err(htrdr, 316 "%s: could not allocate the CIE XYZ data structure -- %s.\n", 317 FUNC_NAME, res_to_cstr(res)); 318 goto error; 319 } 320 ref_init(&cie->ref); 321 darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_X); 322 darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Y); 323 darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Z); 324 cie->range[0] = range[0]; 325 cie->range[1] = range[1]; 326 htrdr_ref_get(htrdr); 327 cie->htrdr = htrdr; 328 329 cie->band_len = (range[1] - range[0]) / (double)nbands; 330 331 res = setup_cie_xyz(cie, FUNC_NAME, nbands); 332 if(res != RES_OK) goto error; 333 334 htrdr_log(htrdr, "CIE XYZ spectral interval defined on [%g, %g] nanometers.\n", 335 range[0], range[1]); 336 337 exit: 338 *out_cie = cie; 339 return res; 340 error: 341 if(cie) htrdr_ran_wlen_cie_xyz_ref_put(cie); 342 goto exit; 343 } 344 345 void 346 htrdr_ran_wlen_cie_xyz_ref_get(struct htrdr_ran_wlen_cie_xyz* cie) 347 { 348 ASSERT(cie); 349 ref_get(&cie->ref); 350 } 351 352 void 353 htrdr_ran_wlen_cie_xyz_ref_put(struct htrdr_ran_wlen_cie_xyz* cie) 354 { 355 ASSERT(cie); 356 ref_put(&cie->ref, release_cie_xyz); 357 } 358 359 double 360 htrdr_ran_wlen_cie_xyz_sample_X 361 (struct htrdr_ran_wlen_cie_xyz* cie, 362 const double r0, 363 const double r1, 364 double* pdf) 365 { 366 const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X), 367 darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1); 368 if(pdf) *pdf = cie->rcp_integral_X; 369 return wlen; 370 } 371 372 double 373 htrdr_ran_wlen_cie_xyz_sample_Y 374 (struct htrdr_ran_wlen_cie_xyz* cie, 375 const double r0, 376 const double r1, 377 double* pdf) 378 { 379 const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y), 380 darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1); 381 if(pdf) *pdf = cie->rcp_integral_Y; 382 return wlen; 383 } 384 385 double 386 htrdr_ran_wlen_cie_xyz_sample_Z 387 (struct htrdr_ran_wlen_cie_xyz* cie, 388 const double r0, 389 const double r1, 390 double* pdf) 391 { 392 const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z), 393 darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1); 394 if(pdf) *pdf = cie->rcp_integral_Z; 395 return wlen; 396 } 397