commit c2cb907802e65c023d683e735398dd7f998b8870
parent 9d002482a618421a40eaf90931850c85b7b6903b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 28 Oct 2022 16:38:30 +0200
Rename htrdr_cie_xyz and htrdr_ran_wlen
Both are wavelength distributions. Therefore, htrdr_cie_xyz is renamed
to htrdr_ran_wlen_cie_xyz and htrdr_ran_wlen is renamed to
htrdr_ran_wlen_planck
Diffstat:
16 files changed, 931 insertions(+), 925 deletions(-)
diff --git a/cmake/core/CMakeLists.txt b/cmake/core/CMakeLists.txt
@@ -84,12 +84,12 @@ set(HTRDR_CORE_FILES_SRC
htrdr.c
htrdr_args.c
htrdr_buffer.c
- htrdr_cie_xyz.c
htrdr_draw_map.c
htrdr_geometry.c
htrdr_log.c
htrdr_materials.c
- htrdr_ran_wlen.c
+ htrdr_ran_wlen_cie_xyz.c
+ htrdr_ran_wlen_planck.c
htrdr_rectangle.c
htrdr_slab.c
htrdr_spectral.c)
@@ -98,13 +98,13 @@ set(HTRDR_CORE_FILES_INC
htrdr_accum.h
htrdr_buffer.h
htrdr_c.h
- htrdr_cie_xyz.h
htrdr_draw_map.h
htrdr_geometry.c
htrdr_interface.h
htrdr_log.h
htrdr_materials.h
- htrdr_ran_wlen.h
+ htrdr_ran_wlen_cie_xyz.h
+ htrdr_ran_wlen_planck.h
htrdr_rectangle.h
htrdr_sensor.h
htrdr_slab.h
diff --git a/src/atmosphere/htrdr_atmosphere.c b/src/atmosphere/htrdr_atmosphere.c
@@ -24,10 +24,10 @@
#include "atmosphere/htrdr_atmosphere_sun.h"
#include "core/htrdr_buffer.h"
-#include "core/htrdr_cie_xyz.h"
#include "core/htrdr_log.h"
#include "core/htrdr_materials.h"
-#include "core/htrdr_ran_wlen.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_ran_wlen_planck.h"
#include "core/htrdr_rectangle.h"
#include <high_tune/htsky.h>
@@ -306,8 +306,8 @@ atmosphere_release(ref_T* ref)
if(cmd->ground) htrdr_atmosphere_ground_ref_put(cmd->ground);
if(cmd->mats) htrdr_materials_ref_put(cmd->mats);
if(cmd->sun) htrdr_atmosphere_sun_ref_put(cmd->sun);
- if(cmd->cie) htrdr_cie_xyz_ref_put(cmd->cie);
- if(cmd->ran_wlen) htrdr_ran_wlen_ref_put(cmd->ran_wlen);
+ if(cmd->cie) htrdr_ran_wlen_cie_xyz_ref_put(cmd->cie);
+ if(cmd->planck) htrdr_ran_wlen_planck_ref_put(cmd->planck);
if(cmd->camera) SCAM(ref_put(cmd->camera));
if(cmd->flux_map) htrdr_rectangle_ref_put(cmd->flux_map);
if(cmd->buf) htrdr_buffer_ref_put(cmd->buf);
@@ -439,7 +439,8 @@ htrdr_atmosphere_create
nintervals = compute_spectral_bands_count(cmd);
if(cmd->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) {
- res = htrdr_cie_xyz_create(htrdr, spectral_range, nintervals, &cmd->cie);
+ res = htrdr_ran_wlen_cie_xyz_create
+ (htrdr, spectral_range, nintervals, &cmd->cie);
if(res != RES_OK) goto error;
} else {
if(cmd->ref_temperature <= 0) {
@@ -448,8 +449,8 @@ htrdr_atmosphere_create
res = RES_BAD_ARG;
goto error;
}
- res = htrdr_ran_wlen_create
- (htrdr, spectral_range, nintervals, cmd->ref_temperature, &cmd->ran_wlen);
+ res = htrdr_ran_wlen_planck_create
+ (htrdr, spectral_range, nintervals, cmd->ref_temperature, &cmd->planck);
if(res != RES_OK) goto error;
}
diff --git a/src/atmosphere/htrdr_atmosphere_c.h b/src/atmosphere/htrdr_atmosphere_c.h
@@ -81,17 +81,17 @@ struct htsky;
struct htrdr;
struct htrdr_atmosphere_args;
struct htrdr_buffer;
-struct htrdr_cie_xyz;
struct htrdr_materials;
-struct htrdr_ran_wlen;
+struct htrdr_ran_wlen_cie_xyz;
+struct htrdr_ran_wlen_planck;
struct ssp_rng;
struct htrdr_atmosphere {
struct htrdr_atmosphere_ground* ground;
struct htrdr_atmosphere_sun* sun;
struct htrdr_materials* mats;
- struct htrdr_cie_xyz* cie;
- struct htrdr_ran_wlen* ran_wlen;
+ struct htrdr_ran_wlen_cie_xyz* cie;
+ struct htrdr_ran_wlen_planck* planck;
struct scam* camera; /* Camera */
struct htrdr_rectangle* flux_map; /* Flux map */
diff --git a/src/atmosphere/htrdr_atmosphere_draw_map.c b/src/atmosphere/htrdr_atmosphere_draw_map.c
@@ -21,11 +21,11 @@
#include "core/htrdr.h"
#include "core/htrdr_buffer.h"
-#include "core/htrdr_cie_xyz.h"
#include "core/htrdr_draw_map.h"
#include "core/htrdr_interface.h"
#include "core/htrdr_log.h"
-#include "core/htrdr_ran_wlen.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_ran_wlen_planck.h"
#include "core/htrdr_rectangle.h"
#include <high_tune/htsky.h>
@@ -160,9 +160,9 @@ draw_pixel_image
/* Sample a spectral band and a quadrature point */
switch(ichannel) {
- case 0: wlen = htrdr_cie_xyz_sample_X(cmd->cie, r0, r1, &pdf); break;
- case 1: wlen = htrdr_cie_xyz_sample_Y(cmd->cie, r0, r1, &pdf); break;
- case 2: wlen = htrdr_cie_xyz_sample_Z(cmd->cie, r0, r1, &pdf); break;
+ case 0: wlen = htrdr_ran_wlen_cie_xyz_sample_X(cmd->cie, r0, r1, &pdf); break;
+ case 1: wlen = htrdr_ran_wlen_cie_xyz_sample_Y(cmd->cie, r0, r1, &pdf); break;
+ case 2: wlen = htrdr_ran_wlen_cie_xyz_sample_Z(cmd->cie, r0, r1, &pdf); break;
default: FATAL("Unreachable code.\n"); break;
}
pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
@@ -249,7 +249,7 @@ draw_pixel_flux
r2 = ssp_rng_canonical(args->rng);
/* Sample a wavelength */
- wlen = htrdr_ran_wlen_sample(cmd->ran_wlen, r0, r1, &band_pdf);
+ wlen = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &band_pdf);
band_pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
/* Select the associated band and sample a quadrature point */
@@ -370,7 +370,7 @@ draw_pixel_xwave
r2 = ssp_rng_canonical(args->rng);
/* Sample a wavelength */
- wlen = htrdr_ran_wlen_sample(cmd->ran_wlen, r0, r1, &band_pdf);
+ wlen = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &band_pdf);
band_pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
/* Select the associated band and sample a quadrature point */
diff --git a/src/core/htrdr_args.c b/src/core/htrdr_args.c
@@ -444,8 +444,8 @@ parse_spectral_parameter(const char* str, void* ptr)
if(!strcmp(key, "cie_xyz")) {
args->spectral_type = HTRDR_SPECTRAL_SW_CIE_XYZ;
- args->wlen_range[0] = HTRDR_CIE_XYZ_RANGE_DEFAULT[0];
- args->wlen_range[1] = HTRDR_CIE_XYZ_RANGE_DEFAULT[1];
+ args->wlen_range[0] = HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[0];
+ args->wlen_range[1] = HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[1];
} else {
if(!val) {
fprintf(stderr, "Missing value to the spectral option `%s'.\n", key);
diff --git a/src/core/htrdr_args.h.in b/src/core/htrdr_args.h.in
@@ -18,7 +18,7 @@
#ifndef HTRDR_ARGS_H
#define HTRDR_ARGS_H
-#include "core/htrdr_cie_xyz.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
#include "core/htrdr_spectral.h"
#include <rsys/rsys.h>
@@ -115,7 +115,7 @@ struct htrdr_args_spectral {
enum htrdr_spectral_type spectral_type;
};
#define HTRDR_ARGS_SPECTRAL_DEFAULT__ { \
- HTRDR_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
+ HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
-1, /* Reference temperature */ \
HTRDR_SPECTRAL_SW_CIE_XYZ, /* Spectral type */ \
}
diff --git a/src/core/htrdr_cie_xyz.c b/src/core/htrdr_cie_xyz.c
@@ -1,389 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#define _POSIX_C_SOURCE 200112L /* nextafter */
-
-#include "core/htrdr.h"
-#include "core/htrdr_c.h"
-#include "core/htrdr_cie_xyz.h"
-#include "core/htrdr_log.h"
-
-#include <rsys/algorithm.h>
-#include <rsys/cstr.h>
-#include <rsys/dynamic_array_double.h>
-#include <rsys/mem_allocator.h>
-#include <rsys/ref_count.h>
-
-#include <math.h> /* nextafter */
-
-struct htrdr_cie_xyz {
- struct darray_double cdf_X;
- struct darray_double cdf_Y;
- struct darray_double cdf_Z;
- double rcp_integral_X;
- double rcp_integral_Y;
- double rcp_integral_Z;
- double range[2]; /* Boundaries of the handled CIE XYZ color space */
- double band_len; /* Length in nanometers of a band */
-
- struct htrdr* htrdr;
- ref_T ref;
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static INLINE double
-trapezoidal_integration
- (const double lambda_lo, /* Integral lower bound. In nanometer */
- const double lambda_hi, /* Integral upper bound. In nanometer */
- double (*f_bar)(const double lambda)) /* Function to integrate */
-{
- double dlambda;
- size_t i, n;
- double integral = 0;
- ASSERT(lambda_lo <= lambda_hi);
- ASSERT(lambda_lo > 0);
-
- n = (size_t)(lambda_hi - lambda_lo) + 1;
- dlambda = (lambda_hi - lambda_lo) / (double)n;
-
- FOR_EACH(i, 0, n) {
- const double lambda1 = lambda_lo + dlambda*(double)(i+0);
- const double lambda2 = lambda_lo + dlambda*(double)(i+1);
- const double f1 = f_bar(lambda1);
- const double f2 = f_bar(lambda2);
- integral += (f1 + f2)*dlambda*0.5;
- }
- return integral;
-}
-
-/* The following 3 functions are used to fit the CIE Xbar, Ybar and Zbar curved
- * has defined by the 1931 standard. These analytical fits are propsed by C.
- * Wyman, P. P. Sloan & P. Shirley in "Simple Analytic Approximations to the
- * CIE XYZ Color Matching Functions" - JCGT 2013. */
-static INLINE double
-fit_x_bar_1931(const double lambda)
-{
- const double a = (lambda - 442.0) * (lambda < 442.0 ? 0.0624 : 0.0374);
- const double b = (lambda - 599.8) * (lambda < 599.8 ? 0.0264 : 0.0323);
- const double c = (lambda - 501.1) * (lambda < 501.1 ? 0.0490 : 0.0382);
- return 0.362*exp(-0.5*a*a) + 1.056*exp(-0.5f*b*b) - 0.065*exp(-0.5*c*c);
-}
-
-static FINLINE double
-fit_y_bar_1931(const double lambda)
-{
- const double a = (lambda - 568.8) * (lambda < 568.8 ? 0.0213 : 0.0247);
- const double b = (lambda - 530.9) * (lambda < 530.9 ? 0.0613 : 0.0322);
- return 0.821*exp(-0.5*a*a) + 0.286*exp(-0.5*b*b);
-}
-
-static FINLINE double
-fit_z_bar_1931(const double lambda)
-{
- const double a = (lambda - 437.0) * (lambda < 437.0 ? 0.0845 : 0.0278);
- const double b = (lambda - 459.0) * (lambda < 459.0 ? 0.0385 : 0.0725);
- return 1.217*exp(-0.5*a*a) + 0.681*exp(-0.5*b*b);
-}
-
-static INLINE double
-sample_cie_xyz
- (const struct htrdr_cie_xyz* cie,
- const double* cdf,
- const size_t cdf_length,
- double (*f_bar)(const double lambda), /* Function to integrate */
- const double r0, /* Canonical number in [0, 1[ */
- const double r1) /* Canonical number in [0, 1[ */
-{
- double r0_next = nextafter(r0, DBL_MAX);
- double* find;
- double f_min, f_max; /* CIE 1931 value for the band boundaries */
- double lambda; /* Sampled wavelength */
- double lambda_min, lambda_max; /* Boundaries of the sampled band */
- double lambda_1, lambda_2; /* Solutions if the equation to solve */
- double a, b, c, d; /* Equation parameters */
- double delta, sqrt_delta;
- size_t iband; /* Index of the sampled band */
- ASSERT(cie && cdf && cdf_length);
- ASSERT(0 <= r0 && r0 < 1);
- ASSERT(0 <= r1 && r1 < 1);
-
- /* Use r_next rather than r in order to find the first entry that is not less
- * than *or equal* to r */
- find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
- ASSERT(find);
-
- /* Define and check the sampled band */
- iband = (size_t)(find - cdf);
- ASSERT(iband < cdf_length);
- ASSERT(cdf[iband] > r0 && (!iband || cdf[iband-1] <= r0));
-
- /* Define the boundaries of the sampled band */
- lambda_min = cie->range[0] + cie->band_len * (double)iband;
- lambda_max = lambda_min + cie->band_len;
-
- /* Define the value of the CIE 1931 function for the boudaries of the sampled
- * band */
- f_min = f_bar(lambda_min);
- f_max = f_bar(lambda_max);
-
- /* Compute the equation constants */
- a = 0.5 * (f_max - f_min) / cie->band_len;
- b = (lambda_max * f_min - lambda_min * f_max) / cie->band_len;
- c = -lambda_min * f_min + lambda_min*lambda_min * a;
- d = 0.5 * (f_max + f_min) * cie->band_len;
-
- delta = b*b - 4*a*(c-d*r1);
- if(delta < 0 && eq_eps(delta, 0, 1.e-6)) {
- delta = 0;
- }
- ASSERT(delta > 0);
- sqrt_delta = sqrt(delta);
-
- /* Compute the roots that solve the equation */
- lambda_1 = (-b - sqrt_delta) / (2*a);
- lambda_2 = (-b + sqrt_delta) / (2*a);
-
- /* Select the solution */
- if(lambda_min <= lambda_1 && lambda_1 < lambda_max) {
- lambda = lambda_1;
- } else if(lambda_min <= lambda_2 && lambda_2 < lambda_max) {
- lambda = lambda_2;
- } else {
- htrdr_log_warn(cie->htrdr,
- "%s: cannot sample a wavelength in [%g, %g[. The possible wavelengths"
- "were %g and %g.\n",
- FUNC_NAME, lambda_min, lambda_max, lambda_1, lambda_2);
- /* Arbitrarly choose the wavelength at the center of the sampled band */
- lambda = (lambda_min + lambda_max)*0.5;
- }
-
- return lambda;
-}
-
-static res_T
-setup_cie_xyz
- (struct htrdr_cie_xyz* cie,
- const char* func_name,
- const size_t nbands)
-{
- enum { X, Y, Z }; /* Helper constant */
- double* pdf[3] = {NULL, NULL, NULL};
- double* cdf[3] = {NULL, NULL, NULL};
- double sum[3] = {0, 0, 0};
- size_t i;
- res_T res = RES_OK;
-
- ASSERT(cie && func_name && nbands);
- ASSERT(cie->range[0] >= HTRDR_CIE_XYZ_RANGE_DEFAULT[0]);
- ASSERT(cie->range[1] <= HTRDR_CIE_XYZ_RANGE_DEFAULT[1]);
- ASSERT(cie->range[0] < cie->range[1]);
-
- /* Allocate and reset the memory space for the tristimulus CDF */
- #define SETUP_STIMULUS(Stimulus) { \
- res = darray_double_resize(&cie->cdf_ ## Stimulus, nbands); \
- if(res != RES_OK) { \
- htrdr_log_err(cie->htrdr, \
- "%s: Could not reserve the memory space for the CDF " \
- "of the "STR(X)" stimulus -- %s.\n", func_name, res_to_cstr(res)); \
- goto error; \
- } \
- cdf[Stimulus] = darray_double_data_get(&cie->cdf_ ## Stimulus); \
- pdf[Stimulus] = cdf[Stimulus]; \
- memset(cdf[Stimulus], 0, nbands*sizeof(double)); \
- } (void)0
- SETUP_STIMULUS(X);
- SETUP_STIMULUS(Y);
- SETUP_STIMULUS(Z);
- #undef SETUP_STIMULUS
-
- /* Compute the *unormalized* pdf of the tristimulus */
- FOR_EACH(i, 0, nbands) {
- const double lambda_lo = cie->range[0] + (double)i * cie->band_len;
- const double lambda_hi = MMIN(lambda_lo + cie->band_len, cie->range[1]);
- ASSERT(lambda_lo <= lambda_hi);
- ASSERT(lambda_lo >= cie->range[0]);
- ASSERT(lambda_hi <= cie->range[1]);
- pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
- pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931);
- pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931);
- sum[X] += pdf[X][i];
- sum[Y] += pdf[Y][i];
- sum[Z] += pdf[Z][i];
- }
- #define CHK_SUM(Sum, Range, Fit) \
- ASSERT(eq_eps(Sum, trapezoidal_integration(Range[0], Range[1], Fit), 1.e-3))
- CHK_SUM(sum[X], cie->range, fit_x_bar_1931);
- CHK_SUM(sum[Y], cie->range, fit_y_bar_1931);
- CHK_SUM(sum[Z], cie->range, fit_z_bar_1931);
- #undef CHK_SUM
- cie->rcp_integral_X = 1.0 / sum[X];
- cie->rcp_integral_Y = 1.0 / sum[Y];
- cie->rcp_integral_Z = 1.0 / sum[Z];
-
- FOR_EACH(i, 0, nbands) {
- /* Normalize the pdf */
- pdf[X][i] /= sum[X];
- pdf[Y][i] /= sum[Y];
- pdf[Z][i] /= sum[Z];
- /* Setup the cumulative */
- if(i == 0) {
- cdf[X][i] = pdf[X][i];
- cdf[Y][i] = pdf[Y][i];
- cdf[Z][i] = pdf[Z][i];
- } else {
- cdf[X][i] = pdf[X][i] + cdf[X][i-1];
- cdf[Y][i] = pdf[Y][i] + cdf[Y][i-1];
- cdf[Z][i] = pdf[Z][i] + cdf[Z][i-1];
- ASSERT(cdf[X][i] >= cdf[X][i-1]);
- ASSERT(cdf[Y][i] >= cdf[Y][i-1]);
- ASSERT(cdf[Z][i] >= cdf[Z][i-1]);
- }
- }
- ASSERT(eq_eps(cdf[X][nbands-1], 1, 1.e-6));
- ASSERT(eq_eps(cdf[Y][nbands-1], 1, 1.e-6));
- ASSERT(eq_eps(cdf[Z][nbands-1], 1, 1.e-6));
-
- /* Handle numerical issue */
- cdf[X][nbands-1] = 1.0;
- cdf[Y][nbands-1] = 1.0;
- cdf[Z][nbands-1] = 1.0;
-
-exit:
- return res;
-error:
- darray_double_clear(&cie->cdf_X);
- darray_double_clear(&cie->cdf_Y);
- darray_double_clear(&cie->cdf_Z);
- goto exit;
-}
-
-static void
-release_cie_xyz(ref_T* ref)
-{
- struct htrdr_cie_xyz* cie = NULL;
- struct htrdr* htrdr = NULL;
- ASSERT(ref);
- cie = CONTAINER_OF(ref, struct htrdr_cie_xyz, ref);
- darray_double_release(&cie->cdf_X);
- darray_double_release(&cie->cdf_Y);
- darray_double_release(&cie->cdf_Z);
- htrdr = cie->htrdr;
- MEM_RM(htrdr_get_allocator(cie->htrdr), cie);
- htrdr_ref_put(htrdr);
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_cie_xyz_create
- (struct htrdr* htrdr,
- const double range[2], /* Must be included in [380, 780] nanometers */
- const size_t bands_count, /* # bands used to discretize the CIE tristimulus */
- struct htrdr_cie_xyz** out_cie)
-{
- struct htrdr_cie_xyz* cie = NULL;
- size_t nbands = bands_count;
- res_T res = RES_OK;
- ASSERT(htrdr && range && nbands && out_cie);
-
- cie = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cie));
- if(!cie) {
- res = RES_MEM_ERR;
- htrdr_log_err(htrdr,
- "%s: could not allocate the CIE XYZ data structure -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
- }
- ref_init(&cie->ref);
- darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_X);
- darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Y);
- darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Z);
- cie->range[0] = range[0];
- cie->range[1] = range[1];
- htrdr_ref_get(htrdr);
- cie->htrdr = htrdr;
-
- cie->band_len = (range[1] - range[0]) / (double)nbands;
-
- res = setup_cie_xyz(cie, FUNC_NAME, nbands);
- if(res != RES_OK) goto error;
-
- htrdr_log(htrdr, "CIE XYZ spectral interval defined on [%g, %g] nanometers.\n",
- range[0], range[1]);
-
-exit:
- *out_cie = cie;
- return res;
-error:
- if(cie) htrdr_cie_xyz_ref_put(cie);
- goto exit;
-}
-
-void
-htrdr_cie_xyz_ref_get(struct htrdr_cie_xyz* cie)
-{
- ASSERT(cie);
- ref_get(&cie->ref);
-}
-
-void
-htrdr_cie_xyz_ref_put(struct htrdr_cie_xyz* cie)
-{
- ASSERT(cie);
- ref_put(&cie->ref, release_cie_xyz);
-}
-
-double
-htrdr_cie_xyz_sample_X
- (struct htrdr_cie_xyz* cie,
- const double r0,
- const double r1,
- double* pdf)
-{
- const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
- darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
- if(pdf) *pdf = cie->rcp_integral_X;
- return wlen;
-}
-
-double
-htrdr_cie_xyz_sample_Y
- (struct htrdr_cie_xyz* cie,
- const double r0,
- const double r1,
- double* pdf)
-{
- const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
- darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
- if(pdf) *pdf = cie->rcp_integral_Y;
- return wlen;
-}
-
-double
-htrdr_cie_xyz_sample_Z
- (struct htrdr_cie_xyz* cie,
- const double r0,
- const double r1,
- double* pdf)
-{
- const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
- darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
- if(pdf) *pdf = cie->rcp_integral_Z;
- return wlen;
-}
-
diff --git a/src/core/htrdr_cie_xyz.h b/src/core/htrdr_cie_xyz.h
@@ -1,72 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTRDR_CIE_XYZ_H
-#define HTRDR_CIE_XYZ_H
-
-#include "core/htrdr.h"
-#include <rsys/rsys.h>
-
-/* Wavelength boundaries of the CIE XYZ color space in nanometers */
-#define HTRDR_CIE_XYZ_RANGE_DEFAULT__ {380, 780}
-static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] =
- HTRDR_CIE_XYZ_RANGE_DEFAULT__;
-
-/* Forward declarations */
-struct htrdr;
-struct htrdr_cie_xyz;
-
-BEGIN_DECLS
-
-HTRDR_CORE_API res_T
-htrdr_cie_xyz_create
- (struct htrdr* htrdr,
- const double range[2], /* Must be included in [380, 780] nanometers */
- const size_t nbands, /* # bands used to discretize the CIE tristimulus */
- struct htrdr_cie_xyz** cie);
-
-HTRDR_CORE_API void
-htrdr_cie_xyz_ref_get
- (struct htrdr_cie_xyz* cie);
-
-HTRDR_CORE_API void
-htrdr_cie_xyz_ref_put
- (struct htrdr_cie_xyz* cie);
-
-/* Return a wavelength in nanometer */
-HTRDR_CORE_API double
-htrdr_cie_xyz_sample_X
- (struct htrdr_cie_xyz* cie,
- const double r0, const double r1, /* Canonical numbers in [0, 1[ */
- double* pdf); /* In nm^-1. May be NULL */
-
-/* Return a wavelength in nanometer */
-HTRDR_CORE_API double
-htrdr_cie_xyz_sample_Y
- (struct htrdr_cie_xyz* cie,
- const double r0, const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* In nm^-1. May be NULL */
-
-/* Return a wavelength in nanometer */
-HTRDR_CORE_API double
-htrdr_cie_xyz_sample_Z
- (struct htrdr_cie_xyz* cie,
- const double r0, const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* In nm^-1. May be NULL */
-
-END_DECLS
-
-#endif /* HTRDR_cie_xyz_H */
-
diff --git a/src/core/htrdr_ran_wlen.c b/src/core/htrdr_ran_wlen.c
@@ -1,364 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019, 2021 CNRS
- * Copyright (C) 2018, 2019 Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#define _POSIX_C_SOURCE 200112L /* nextafter */
-
-#include "core/htrdr.h"
-#include "core/htrdr_c.h"
-#include "core/htrdr_log.h"
-#include "core/htrdr_ran_wlen.h"
-#include "core/htrdr_spectral.h"
-
-#include <rsys/algorithm.h>
-#include <rsys/cstr.h>
-#include <rsys/dynamic_array_double.h>
-#include <rsys/mem_allocator.h>
-#include <rsys/ref_count.h>
-
-#include <math.h> /* nextafter */
-
-struct htrdr_ran_wlen {
- struct darray_double pdf;
- struct darray_double cdf;
- double range[2]; /* Boundaries of the spectral integration interval */
- double band_len; /* Length in nanometers of a band */
- double ref_temperature; /* In Kelvin */
- size_t nbands; /* # bands */
-
- struct htrdr* htrdr;
- ref_T ref;
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static res_T
-setup_wlen_ran_cdf
- (struct htrdr_ran_wlen* wlen_ran,
- const char* func_name)
-{
- double* pdf = NULL;
- double* cdf = NULL;
- double sum = 0;
- size_t i;
- res_T res = RES_OK;
- ASSERT(wlen_ran && func_name && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
-
- res = darray_double_resize(&wlen_ran->cdf, wlen_ran->nbands);
- if(res != RES_OK) {
- htrdr_log_err(wlen_ran->htrdr,
- "%s: Error allocating the CDF of the spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
- res = darray_double_resize(&wlen_ran->pdf, wlen_ran->nbands);
- if(res != RES_OK) {
- htrdr_log_err(wlen_ran->htrdr,
- "%s: Error allocating the pdf of the spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
-
- cdf = darray_double_data_get(&wlen_ran->cdf);
- pdf = darray_double_data_get(&wlen_ran->pdf); /* Now save pdf to correct weight */
-
- htrdr_log(wlen_ran->htrdr,
- "Number of bands of the spectrum cumulative: %lu\n",
- (unsigned long)wlen_ran->nbands);
-
- /* Compute the *unnormalized* probability to sample a small band */
- FOR_EACH(i, 0, wlen_ran->nbands) {
- double lambda_lo = wlen_ran->range[0] + (double)i * wlen_ran->band_len;
- double lambda_hi = MMIN(lambda_lo + wlen_ran->band_len, wlen_ran->range[1]);
- ASSERT(lambda_lo<= lambda_hi);
- ASSERT(lambda_lo > wlen_ran->range[0]
- || eq_eps(lambda_lo, wlen_ran->range[0], 1.e-6));
- ASSERT(lambda_lo < wlen_ran->range[1]
- || eq_eps(lambda_lo, wlen_ran->range[1], 1.e-6));
-
- /* Convert from nanometer to meter */
- lambda_lo *= 1.e-9;
- lambda_hi *= 1.e-9;
-
- /* Compute the probability of the current band */
- pdf[i] = htrdr_blackbody_fraction
- (lambda_lo, lambda_hi, wlen_ran->ref_temperature);
-
- /* Update the norm */
- sum += pdf[i];
- }
-
- /* Compute the cumulative of the previously computed probabilities */
- FOR_EACH(i, 0, wlen_ran->nbands) {
- /* Normalize the probability */
- pdf[i] /= sum;
-
- /* Setup the cumulative */
- if(i == 0) {
- cdf[i] = pdf[i];
- } else {
- cdf[i] = pdf[i] + cdf[i-1];
- ASSERT(cdf[i] >= cdf[i-1]);
- }
- }
- ASSERT(eq_eps(cdf[wlen_ran->nbands-1], 1, 1.e-6));
- cdf[wlen_ran->nbands - 1] = 1.0; /* Handle numerical issue */
-
-exit:
- return res;
-error:
- darray_double_clear(&wlen_ran->cdf);
- darray_double_clear(&wlen_ran->pdf);
- goto exit;
-}
-
-static double
-wlen_ran_sample_continue
- (const struct htrdr_ran_wlen* wlen_ran,
- const double r,
- const double range[2], /* In nanometer */
- const char* func_name,
- double* pdf)
-{
- /* Numerical parameters of the dichotomy search */
- const size_t MAX_ITER = 100;
- const double EPSILON_LAMBDA_M = 1e-15;
- const double EPSILON_BF = 1e-6;
-
- /* Local variables */
- double bf = 0;
- double bf_prev = 0;
- double bf_min_max = 0;
- double lambda_m = 0;
- double lambda_m_prev = 0;
- double lambda_m_min = 0;
- double lambda_m_max = 0;
- double range_m[2] = {0, 0};
- size_t i;
-
- /* Check precondition */
- ASSERT(wlen_ran && func_name);
- ASSERT(range && range[0] < range[1]);
- ASSERT(0 <= r && r < 1);
-
- /* Convert the wavelength range in meters */
- range_m[0] = range[0] * 1.0e-9;
- range_m[1] = range[1] * 1.0e-9;
-
- /* Setup the dichotomy search */
- lambda_m_min = range_m[0];
- lambda_m_max = range_m[1];
- bf_min_max = htrdr_blackbody_fraction
- (range_m[0], range_m[1], wlen_ran->ref_temperature);
-
- /* Numerically search the lambda corresponding to the submitted canonical
- * number */
- FOR_EACH(i, 0, MAX_ITER) {
- double r_test;
- lambda_m = (lambda_m_min + lambda_m_max) * 0.5;
- bf = htrdr_blackbody_fraction
- (range_m[0], lambda_m, wlen_ran->ref_temperature);
-
- r_test = bf / bf_min_max;
- if(r_test < r) {
- lambda_m_min = lambda_m;
- } else {
- lambda_m_max = lambda_m;
- }
-
- if(fabs(lambda_m_prev - lambda_m) < EPSILON_LAMBDA_M
- || fabs(bf_prev - bf) < EPSILON_BF)
- break;
-
- lambda_m_prev = lambda_m;
- bf_prev = bf;
- }
- if(i >= MAX_ITER) {
- htrdr_log_warn(wlen_ran->htrdr,
- "%s: could not sample a wavelength in the range [%g, %g] nanometers "
- "for the reference temperature %g Kelvin.\n",
- func_name, SPLIT2(range), wlen_ran->ref_temperature);
- }
-
- if(pdf) {
- const double Tref = wlen_ran->ref_temperature;
- const double B_lambda = htrdr_planck(lambda_m, lambda_m, Tref);
- const double B_mean = htrdr_planck(range_m[0], range_m[1], Tref);
- *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0]));
- *pdf *= 1.e-9; /* Transform from m^-1 to nm^-1 */
- }
-
- return lambda_m * 1.e9; /* Convert in nanometers */
-}
-
-static double
-wlen_ran_sample_discrete
- (const struct htrdr_ran_wlen* wlen_ran,
- const double r0,
- const double r1,
- const char* func_name,
- double* pdf)
-{
- const double* cdf = NULL;
- const double* find = NULL;
- double r0_next = nextafter(r0, DBL_MAX);
- double band_range[2];
- double lambda = 0;
- double pdf_continue = 0;
- double pdf_band = 0;
- size_t cdf_length = 0;
- size_t i;
- ASSERT(wlen_ran && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
- ASSERT(0 <= r0 && r0 < 1);
- ASSERT(0 <= r1 && r1 < 1);
- (void)func_name;
- (void)pdf_band;
-
- cdf = darray_double_cdata_get(&wlen_ran->cdf);
- cdf_length = darray_double_size_get(&wlen_ran->cdf);
- ASSERT(cdf_length > 0);
-
- /* Use r_next rather than r0 in order to find the first entry that is not less
- * than *or equal* to r0 */
- find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
- ASSERT(find);
-
- i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r0 && (!i || cdf[i-1] <= r0));
-
- band_range[0] = wlen_ran->range[0] + (double)i*wlen_ran->band_len;
- band_range[1] = band_range[0] + wlen_ran->band_len;
-
- /* Fetch the pdf of the sampled band */
- pdf_band = darray_double_cdata_get(&wlen_ran->pdf)[i];
-
- /* Uniformly sample a wavelength in the sampled band */
- lambda = band_range[0] + (band_range[1] - band_range[0]) * r1;
- pdf_continue = 1.0 / (band_range[1] - band_range[0]);
-
- if(pdf) {
- *pdf = pdf_band * pdf_continue;
- }
-
- return lambda;
-}
-
-static void
-release_wlen_ran(ref_T* ref)
-{
- struct htrdr_ran_wlen* wlen_ran = NULL;
- struct htrdr* htrdr = NULL;
- ASSERT(ref);
- wlen_ran = CONTAINER_OF(ref, struct htrdr_ran_wlen, ref);
- darray_double_release(&wlen_ran->cdf);
- darray_double_release(&wlen_ran->pdf);
- htrdr = wlen_ran->htrdr;
- MEM_RM(htrdr_get_allocator(htrdr), wlen_ran);
- htrdr_ref_put(wlen_ran->htrdr);
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_ran_wlen_create
- (struct htrdr* htrdr,
- /* range must be included in [200,1000] nm for shortwave or in [1000,100000]
- * nanometers for longwave (thermal) */
- const double range[2],
- const size_t nbands, /* # bands used to discretized CDF */
- const double ref_temperature,
- struct htrdr_ran_wlen** out_wlen_ran)
-{
- struct htrdr_ran_wlen* wlen_ran = NULL;
- res_T res = RES_OK;
- ASSERT(htrdr && range && out_wlen_ran && ref_temperature > 0);
- ASSERT(ref_temperature > 0);
- ASSERT(range[0] <= range[1]);
-
- wlen_ran = MEM_CALLOC(htrdr->allocator, 1, sizeof(*wlen_ran));
- if(!wlen_ran) {
- res = RES_MEM_ERR;
- htrdr_log_err(htrdr,
- "%s: could not allocate longwave random variate data structure -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
- }
- ref_init(&wlen_ran->ref);
- darray_double_init(htrdr->allocator, &wlen_ran->cdf);
- darray_double_init(htrdr->allocator, &wlen_ran->pdf);
- htrdr_ref_get(htrdr);
- wlen_ran->htrdr = htrdr;
-
- wlen_ran->range[0] = range[0];
- wlen_ran->range[1] = range[1];
- wlen_ran->ref_temperature = ref_temperature;
- wlen_ran->nbands = nbands;
-
- if(nbands == HTRDR_WLEN_RAN_CONTINUE) {
- wlen_ran->band_len = 0;
- } else {
- wlen_ran->band_len = (range[1] - range[0]) / (double)wlen_ran->nbands;
-
- res = setup_wlen_ran_cdf(wlen_ran, FUNC_NAME);
- if(res != RES_OK) goto error;
- }
-
- htrdr_log(htrdr, "Spectral interval defined on [%g, %g] nanometers.\n",
- range[0], range[1]);
-
-exit:
- *out_wlen_ran = wlen_ran;
- return res;
-error:
- if(wlen_ran) htrdr_ran_wlen_ref_put(wlen_ran);
- goto exit;
-}
-
-void
-htrdr_ran_wlen_ref_get(struct htrdr_ran_wlen* wlen_ran)
-{
- ASSERT(wlen_ran);
- ref_get(&wlen_ran->ref);
-}
-
-void
-htrdr_ran_wlen_ref_put(struct htrdr_ran_wlen* wlen_ran)
-{
- ASSERT(wlen_ran);
- ref_put(&wlen_ran->ref, release_wlen_ran);
-}
-
-double
-htrdr_ran_wlen_sample
- (const struct htrdr_ran_wlen* wlen_ran,
- const double r0,
- const double r1,
- double* pdf)
-{
- ASSERT(wlen_ran);
- if(wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE) { /* Discrete */
- return wlen_ran_sample_discrete(wlen_ran, r0, r1, FUNC_NAME, pdf);
- } else if(eq_eps(wlen_ran->range[0], wlen_ran->range[1], 1.e-6)) {
- if(pdf) *pdf = 1;
- return wlen_ran->range[0];
- } else { /* Continue */
- return wlen_ran_sample_continue
- (wlen_ran, r0, wlen_ran->range, FUNC_NAME, pdf);
- }
-}
-
diff --git a/src/core/htrdr_ran_wlen.h b/src/core/htrdr_ran_wlen.h
@@ -1,61 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019, 2021 CNRS
- * Copyright (C) 2018, 2019 Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTRDR_RAN_WLEN_H
-#define HTRDR_RAN_WLEN_H
-
-#include "core/htrdr.h"
-#include <rsys/rsys.h>
-
-#define HTRDR_WLEN_RAN_CONTINUE 0
-
-/* Forward declarations */
-struct htrdr;
-struct htrdr_ran_wlen;
-
-BEGIN_DECLS
-
-HTRDR_CORE_API res_T
-htrdr_ran_wlen_create
- (struct htrdr* htrdr,
- const double range[2],
- /* # bands used to discretisze the spectral domain. HTRDR_WLEN_RAN_CONTINUE
- * <=> no discretisation */
- const size_t nbands, /* Hint on #bands used to discretised th CDF */
- const double ref_temperature, /* Reference temperature */
- struct htrdr_ran_wlen** wlen_ran);
-
-HTRDR_CORE_API void
-htrdr_ran_wlen_ref_get
- (struct htrdr_ran_wlen* wlen_ran);
-
-HTRDR_CORE_API void
-htrdr_ran_wlen_ref_put
- (struct htrdr_ran_wlen* wlen_ran);
-
-/* Return a wavelength in nanometer */
-HTRDR_CORE_API double
-htrdr_ran_wlen_sample
- (const struct htrdr_ran_wlen* wlen_ran,
- const double r0, /* Canonical number in [0, 1[ */
- const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* In nm^-1. May be NULL */
-
-END_DECLS
-
-#endif /* HTRDR_RAN_WLEN_H */
-
diff --git a/src/core/htrdr_ran_wlen_cie_xyz.c b/src/core/htrdr_ran_wlen_cie_xyz.c
@@ -0,0 +1,391 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "core/htrdr.h"
+#include "core/htrdr_c.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_log.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <math.h> /* nextafter */
+
+struct htrdr_ran_wlen_cie_xyz {
+ struct darray_double cdf_X;
+ struct darray_double cdf_Y;
+ struct darray_double cdf_Z;
+ double rcp_integral_X;
+ double rcp_integral_Y;
+ double rcp_integral_Z;
+ double range[2]; /* Boundaries of the handled CIE XYZ color space */
+ double band_len; /* Length in nanometers of a band */
+
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE double
+trapezoidal_integration
+ (const double lambda_lo, /* Integral lower bound. In nanometer */
+ const double lambda_hi, /* Integral upper bound. In nanometer */
+ double (*f_bar)(const double lambda)) /* Function to integrate */
+{
+ double dlambda;
+ size_t i, n;
+ double integral = 0;
+ ASSERT(lambda_lo <= lambda_hi);
+ ASSERT(lambda_lo > 0);
+
+ n = (size_t)(lambda_hi - lambda_lo) + 1;
+ dlambda = (lambda_hi - lambda_lo) / (double)n;
+
+ FOR_EACH(i, 0, n) {
+ const double lambda1 = lambda_lo + dlambda*(double)(i+0);
+ const double lambda2 = lambda_lo + dlambda*(double)(i+1);
+ const double f1 = f_bar(lambda1);
+ const double f2 = f_bar(lambda2);
+ integral += (f1 + f2)*dlambda*0.5;
+ }
+ return integral;
+}
+
+/* The following 3 functions are used to fit the CIE Xbar, Ybar and Zbar curved
+ * has defined by the 1931 standard. These analytical fits are propsed by C.
+ * Wyman, P. P. Sloan & P. Shirley in "Simple Analytic Approximations to the
+ * CIE XYZ Color Matching Functions" - JCGT 2013. */
+static INLINE double
+fit_x_bar_1931(const double lambda)
+{
+ const double a = (lambda - 442.0) * (lambda < 442.0 ? 0.0624 : 0.0374);
+ const double b = (lambda - 599.8) * (lambda < 599.8 ? 0.0264 : 0.0323);
+ const double c = (lambda - 501.1) * (lambda < 501.1 ? 0.0490 : 0.0382);
+ return 0.362*exp(-0.5*a*a) + 1.056*exp(-0.5f*b*b) - 0.065*exp(-0.5*c*c);
+}
+
+static FINLINE double
+fit_y_bar_1931(const double lambda)
+{
+ const double a = (lambda - 568.8) * (lambda < 568.8 ? 0.0213 : 0.0247);
+ const double b = (lambda - 530.9) * (lambda < 530.9 ? 0.0613 : 0.0322);
+ return 0.821*exp(-0.5*a*a) + 0.286*exp(-0.5*b*b);
+}
+
+static FINLINE double
+fit_z_bar_1931(const double lambda)
+{
+ const double a = (lambda - 437.0) * (lambda < 437.0 ? 0.0845 : 0.0278);
+ const double b = (lambda - 459.0) * (lambda < 459.0 ? 0.0385 : 0.0725);
+ return 1.217*exp(-0.5*a*a) + 0.681*exp(-0.5*b*b);
+}
+
+static INLINE double
+sample_cie_xyz
+ (const struct htrdr_ran_wlen_cie_xyz* cie,
+ const double* cdf,
+ const size_t cdf_length,
+ double (*f_bar)(const double lambda), /* Function to integrate */
+ const double r0, /* Canonical number in [0, 1[ */
+ const double r1) /* Canonical number in [0, 1[ */
+{
+ double r0_next = nextafter(r0, DBL_MAX);
+ double* find;
+ double f_min, f_max; /* CIE 1931 value for the band boundaries */
+ double lambda; /* Sampled wavelength */
+ double lambda_min, lambda_max; /* Boundaries of the sampled band */
+ double lambda_1, lambda_2; /* Solutions if the equation to solve */
+ double a, b, c, d; /* Equation parameters */
+ double delta, sqrt_delta;
+ size_t iband; /* Index of the sampled band */
+ ASSERT(cie && cdf && cdf_length);
+ ASSERT(0 <= r0 && r0 < 1);
+ ASSERT(0 <= r1 && r1 < 1);
+
+ /* Use r_next rather than r in order to find the first entry that is not less
+ * than *or equal* to r */
+ find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ ASSERT(find);
+
+ /* Define and check the sampled band */
+ iband = (size_t)(find - cdf);
+ ASSERT(iband < cdf_length);
+ ASSERT(cdf[iband] > r0 && (!iband || cdf[iband-1] <= r0));
+
+ /* Define the boundaries of the sampled band */
+ lambda_min = cie->range[0] + cie->band_len * (double)iband;
+ lambda_max = lambda_min + cie->band_len;
+
+ /* Define the value of the CIE 1931 function for the boudaries of the sampled
+ * band */
+ f_min = f_bar(lambda_min);
+ f_max = f_bar(lambda_max);
+
+ /* Compute the equation constants */
+ a = 0.5 * (f_max - f_min) / cie->band_len;
+ b = (lambda_max * f_min - lambda_min * f_max) / cie->band_len;
+ c = -lambda_min * f_min + lambda_min*lambda_min * a;
+ d = 0.5 * (f_max + f_min) * cie->band_len;
+
+ delta = b*b - 4*a*(c-d*r1);
+ if(delta < 0 && eq_eps(delta, 0, 1.e-6)) {
+ delta = 0;
+ }
+ ASSERT(delta > 0);
+ sqrt_delta = sqrt(delta);
+
+ /* Compute the roots that solve the equation */
+ lambda_1 = (-b - sqrt_delta) / (2*a);
+ lambda_2 = (-b + sqrt_delta) / (2*a);
+
+ /* Select the solution */
+ if(lambda_min <= lambda_1 && lambda_1 < lambda_max) {
+ lambda = lambda_1;
+ } else if(lambda_min <= lambda_2 && lambda_2 < lambda_max) {
+ lambda = lambda_2;
+ } else {
+ htrdr_log_warn(cie->htrdr,
+ "%s: cannot sample a wavelength in [%g, %g[. The possible wavelengths"
+ "were %g and %g.\n",
+ FUNC_NAME, lambda_min, lambda_max, lambda_1, lambda_2);
+ /* Arbitrarly choose the wavelength at the center of the sampled band */
+ lambda = (lambda_min + lambda_max)*0.5;
+ }
+
+ return lambda;
+}
+
+static res_T
+setup_cie_xyz
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const char* func_name,
+ const size_t nbands)
+{
+ enum { X, Y, Z }; /* Helper constant */
+ double* pdf[3] = {NULL, NULL, NULL};
+ double* cdf[3] = {NULL, NULL, NULL};
+ double sum[3] = {0, 0, 0};
+ size_t i;
+ res_T res = RES_OK;
+
+ ASSERT(cie && func_name && nbands);
+ ASSERT(cie->range[0] >= HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[0]);
+ ASSERT(cie->range[1] <= HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[1]);
+ ASSERT(cie->range[0] < cie->range[1]);
+
+ /* Allocate and reset the memory space for the tristimulus CDF */
+ #define SETUP_STIMULUS(Stimulus) { \
+ res = darray_double_resize(&cie->cdf_ ## Stimulus, nbands); \
+ if(res != RES_OK) { \
+ htrdr_log_err(cie->htrdr, \
+ "%s: Could not reserve the memory space for the CDF " \
+ "of the "STR(X)" stimulus -- %s.\n", func_name, res_to_cstr(res)); \
+ goto error; \
+ } \
+ cdf[Stimulus] = darray_double_data_get(&cie->cdf_ ## Stimulus); \
+ pdf[Stimulus] = cdf[Stimulus]; \
+ memset(cdf[Stimulus], 0, nbands*sizeof(double)); \
+ } (void)0
+ SETUP_STIMULUS(X);
+ SETUP_STIMULUS(Y);
+ SETUP_STIMULUS(Z);
+ #undef SETUP_STIMULUS
+
+ /* Compute the *unormalized* pdf of the tristimulus */
+ FOR_EACH(i, 0, nbands) {
+ const double lambda_lo = cie->range[0] + (double)i * cie->band_len;
+ const double lambda_hi = MMIN(lambda_lo + cie->band_len, cie->range[1]);
+ ASSERT(lambda_lo <= lambda_hi);
+ ASSERT(lambda_lo >= cie->range[0]);
+ ASSERT(lambda_hi <= cie->range[1]);
+ pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
+ pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931);
+ pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931);
+ sum[X] += pdf[X][i];
+ sum[Y] += pdf[Y][i];
+ sum[Z] += pdf[Z][i];
+ }
+ #define CHK_SUM(Sum, Range, Fit) \
+ ASSERT(eq_eps(Sum, trapezoidal_integration(Range[0], Range[1], Fit), 1.e-3))
+ CHK_SUM(sum[X], cie->range, fit_x_bar_1931);
+ CHK_SUM(sum[Y], cie->range, fit_y_bar_1931);
+ CHK_SUM(sum[Z], cie->range, fit_z_bar_1931);
+ #undef CHK_SUM
+ cie->rcp_integral_X = 1.0 / sum[X];
+ cie->rcp_integral_Y = 1.0 / sum[Y];
+ cie->rcp_integral_Z = 1.0 / sum[Z];
+
+ FOR_EACH(i, 0, nbands) {
+ /* Normalize the pdf */
+ pdf[X][i] /= sum[X];
+ pdf[Y][i] /= sum[Y];
+ pdf[Z][i] /= sum[Z];
+ /* Setup the cumulative */
+ if(i == 0) {
+ cdf[X][i] = pdf[X][i];
+ cdf[Y][i] = pdf[Y][i];
+ cdf[Z][i] = pdf[Z][i];
+ } else {
+ cdf[X][i] = pdf[X][i] + cdf[X][i-1];
+ cdf[Y][i] = pdf[Y][i] + cdf[Y][i-1];
+ cdf[Z][i] = pdf[Z][i] + cdf[Z][i-1];
+ ASSERT(cdf[X][i] >= cdf[X][i-1]);
+ ASSERT(cdf[Y][i] >= cdf[Y][i-1]);
+ ASSERT(cdf[Z][i] >= cdf[Z][i-1]);
+ }
+ }
+ ASSERT(eq_eps(cdf[X][nbands-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Y][nbands-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Z][nbands-1], 1, 1.e-6));
+
+ /* Handle numerical issue */
+ cdf[X][nbands-1] = 1.0;
+ cdf[Y][nbands-1] = 1.0;
+ cdf[Z][nbands-1] = 1.0;
+
+exit:
+ return res;
+error:
+ darray_double_clear(&cie->cdf_X);
+ darray_double_clear(&cie->cdf_Y);
+ darray_double_clear(&cie->cdf_Z);
+ goto exit;
+}
+
+static void
+release_cie_xyz(ref_T* ref)
+{
+ struct htrdr_ran_wlen_cie_xyz* cie = NULL;
+ struct htrdr* htrdr = NULL;
+ ASSERT(ref);
+ cie = CONTAINER_OF(ref, struct htrdr_ran_wlen_cie_xyz, ref);
+ darray_double_release(&cie->cdf_X);
+ darray_double_release(&cie->cdf_Y);
+ darray_double_release(&cie->cdf_Z);
+ htrdr = cie->htrdr;
+ MEM_RM(htrdr_get_allocator(cie->htrdr), cie);
+ htrdr_ref_put(htrdr);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_ran_wlen_cie_xyz_create
+ (struct htrdr* htrdr,
+ const double range[2], /* Must be included in [380, 780] nanometers */
+ const size_t bands_count, /* # bands used to discretize the CIE tristimulus */
+ struct htrdr_ran_wlen_cie_xyz** out_cie)
+{
+ struct htrdr_ran_wlen_cie_xyz* cie = NULL;
+ size_t nbands = bands_count;
+ res_T res = RES_OK;
+ ASSERT(htrdr && range && nbands && out_cie);
+
+ cie = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cie));
+ if(!cie) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: could not allocate the CIE XYZ data structure -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ ref_init(&cie->ref);
+ darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_X);
+ darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Y);
+ darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Z);
+ cie->range[0] = range[0];
+ cie->range[1] = range[1];
+ htrdr_ref_get(htrdr);
+ cie->htrdr = htrdr;
+
+ cie->band_len = (range[1] - range[0]) / (double)nbands;
+
+ res = setup_cie_xyz(cie, FUNC_NAME, nbands);
+ if(res != RES_OK) goto error;
+
+ htrdr_log(htrdr, "CIE XYZ spectral interval defined on [%g, %g] nanometers.\n",
+ range[0], range[1]);
+
+exit:
+ *out_cie = cie;
+ return res;
+error:
+ if(cie) htrdr_ran_wlen_cie_xyz_ref_put(cie);
+ goto exit;
+}
+
+void
+htrdr_ran_wlen_cie_xyz_ref_get(struct htrdr_ran_wlen_cie_xyz* cie)
+{
+ ASSERT(cie);
+ ref_get(&cie->ref);
+}
+
+void
+htrdr_ran_wlen_cie_xyz_ref_put(struct htrdr_ran_wlen_cie_xyz* cie)
+{
+ ASSERT(cie);
+ ref_put(&cie->ref, release_cie_xyz);
+}
+
+double
+htrdr_ran_wlen_cie_xyz_sample_X
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
+ darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_X;
+ return wlen;
+}
+
+double
+htrdr_ran_wlen_cie_xyz_sample_Y
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
+ darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_Y;
+ return wlen;
+}
+
+double
+htrdr_ran_wlen_cie_xyz_sample_Z
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
+ darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_Z;
+ return wlen;
+}
+
diff --git a/src/core/htrdr_ran_wlen_cie_xyz.h b/src/core/htrdr_ran_wlen_cie_xyz.h
@@ -0,0 +1,74 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_RAN_WLEN_CIE_XYZ_H
+#define HTRDR_RAN_WLEN_CIE_XYZ_H
+
+#include "core/htrdr.h"
+#include <rsys/rsys.h>
+
+/* Wavelength boundaries of the CIE XYZ color space in nanometers */
+#define HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT__ {380, 780}
+static const double HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[2] =
+ HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT__;
+
+/* Forward declarations */
+struct htrdr;
+struct htrdr_ran_wlen_cie_xyz;
+
+BEGIN_DECLS
+
+HTRDR_CORE_API res_T
+htrdr_ran_wlen_cie_xyz_create
+ (struct htrdr* htrdr,
+ const double range[2], /* Must be included in [380, 780] nanometers */
+ const size_t nbands, /* # bands used to discretize the CIE tristimulus */
+ struct htrdr_ran_wlen_cie_xyz** cie);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_cie_xyz_ref_get
+ (struct htrdr_ran_wlen_cie_xyz* cie);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_cie_xyz_ref_put
+ (struct htrdr_ran_wlen_cie_xyz* cie);
+
+/* Return a wavelength in nanometer */
+HTRDR_CORE_API double
+htrdr_ran_wlen_cie_xyz_sample_X
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0, const double r1, /* Canonical numbers in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
+
+/* Return a wavelength in nanometer */
+HTRDR_CORE_API double
+htrdr_ran_wlen_cie_xyz_sample_Y
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
+
+/* Return a wavelength in nanometer */
+HTRDR_CORE_API double
+htrdr_ran_wlen_cie_xyz_sample_Z
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
+
+END_DECLS
+
+#endif /* HTRDR_RAN_WLEN_CIE_XYZ_H */
+
diff --git a/src/core/htrdr_ran_wlen_planck.c b/src/core/htrdr_ran_wlen_planck.c
@@ -0,0 +1,365 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019 Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "core/htrdr.h"
+#include "core/htrdr_c.h"
+#include "core/htrdr_log.h"
+#include "core/htrdr_ran_wlen_planck.h"
+#include "core/htrdr_spectral.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <math.h> /* nextafter */
+
+struct htrdr_ran_wlen_planck {
+ struct darray_double pdf;
+ struct darray_double cdf;
+ double range[2]; /* Boundaries of the spectral integration interval */
+ double band_len; /* Length in nanometers of a band */
+ double ref_temperature; /* In Kelvin */
+ size_t nbands; /* # bands */
+
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+setup_wlen_planck_ran_cdf
+ (struct htrdr_ran_wlen_planck* planck,
+ const char* func_name)
+{
+ double* pdf = NULL;
+ double* cdf = NULL;
+ double sum = 0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(planck && func_name);
+ ASSERT(planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE);
+
+ res = darray_double_resize(&planck->cdf, planck->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(planck->htrdr,
+ "%s: Error allocating the CDF of the spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+ res = darray_double_resize(&planck->pdf, planck->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(planck->htrdr,
+ "%s: Error allocating the pdf of the spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+
+ cdf = darray_double_data_get(&planck->cdf);
+ pdf = darray_double_data_get(&planck->pdf); /* Now save pdf to correct weight */
+
+ htrdr_log(planck->htrdr,
+ "Number of bands of the spectrum cumulative: %lu\n",
+ (unsigned long)planck->nbands);
+
+ /* Compute the *unnormalized* probability to sample a small band */
+ FOR_EACH(i, 0, planck->nbands) {
+ double lambda_lo = planck->range[0] + (double)i * planck->band_len;
+ double lambda_hi = MMIN(lambda_lo + planck->band_len, planck->range[1]);
+ ASSERT(lambda_lo<= lambda_hi);
+ ASSERT(lambda_lo > planck->range[0]
+ || eq_eps(lambda_lo, planck->range[0], 1.e-6));
+ ASSERT(lambda_lo < planck->range[1]
+ || eq_eps(lambda_lo, planck->range[1], 1.e-6));
+
+ /* Convert from nanometer to meter */
+ lambda_lo *= 1.e-9;
+ lambda_hi *= 1.e-9;
+
+ /* Compute the probability of the current band */
+ pdf[i] = htrdr_blackbody_fraction
+ (lambda_lo, lambda_hi, planck->ref_temperature);
+
+ /* Update the norm */
+ sum += pdf[i];
+ }
+
+ /* Compute the cumulative of the previously computed probabilities */
+ FOR_EACH(i, 0, planck->nbands) {
+ /* Normalize the probability */
+ pdf[i] /= sum;
+
+ /* Setup the cumulative */
+ if(i == 0) {
+ cdf[i] = pdf[i];
+ } else {
+ cdf[i] = pdf[i] + cdf[i-1];
+ ASSERT(cdf[i] >= cdf[i-1]);
+ }
+ }
+ ASSERT(eq_eps(cdf[planck->nbands-1], 1, 1.e-6));
+ cdf[planck->nbands - 1] = 1.0; /* Handle numerical issue */
+
+exit:
+ return res;
+error:
+ darray_double_clear(&planck->cdf);
+ darray_double_clear(&planck->pdf);
+ goto exit;
+}
+
+static double
+wlen_ran_sample_continue
+ (const struct htrdr_ran_wlen_planck* planck,
+ const double r,
+ const double range[2], /* In nanometer */
+ const char* func_name,
+ double* pdf)
+{
+ /* Numerical parameters of the dichotomy search */
+ const size_t MAX_ITER = 100;
+ const double EPSILON_LAMBDA_M = 1e-15;
+ const double EPSILON_BF = 1e-6;
+
+ /* Local variables */
+ double bf = 0;
+ double bf_prev = 0;
+ double bf_min_max = 0;
+ double lambda_m = 0;
+ double lambda_m_prev = 0;
+ double lambda_m_min = 0;
+ double lambda_m_max = 0;
+ double range_m[2] = {0, 0};
+ size_t i;
+
+ /* Check precondition */
+ ASSERT(planck && func_name);
+ ASSERT(range && range[0] < range[1]);
+ ASSERT(0 <= r && r < 1);
+
+ /* Convert the wavelength range in meters */
+ range_m[0] = range[0] * 1.0e-9;
+ range_m[1] = range[1] * 1.0e-9;
+
+ /* Setup the dichotomy search */
+ lambda_m_min = range_m[0];
+ lambda_m_max = range_m[1];
+ bf_min_max = htrdr_blackbody_fraction
+ (range_m[0], range_m[1], planck->ref_temperature);
+
+ /* Numerically search the lambda corresponding to the submitted canonical
+ * number */
+ FOR_EACH(i, 0, MAX_ITER) {
+ double r_test;
+ lambda_m = (lambda_m_min + lambda_m_max) * 0.5;
+ bf = htrdr_blackbody_fraction
+ (range_m[0], lambda_m, planck->ref_temperature);
+
+ r_test = bf / bf_min_max;
+ if(r_test < r) {
+ lambda_m_min = lambda_m;
+ } else {
+ lambda_m_max = lambda_m;
+ }
+
+ if(fabs(lambda_m_prev - lambda_m) < EPSILON_LAMBDA_M
+ || fabs(bf_prev - bf) < EPSILON_BF)
+ break;
+
+ lambda_m_prev = lambda_m;
+ bf_prev = bf;
+ }
+ if(i >= MAX_ITER) {
+ htrdr_log_warn(planck->htrdr,
+ "%s: could not sample a wavelength in the range [%g, %g] nanometers "
+ "for the reference temperature %g Kelvin.\n",
+ func_name, SPLIT2(range), planck->ref_temperature);
+ }
+
+ if(pdf) {
+ const double Tref = planck->ref_temperature;
+ const double B_lambda = htrdr_planck(lambda_m, lambda_m, Tref);
+ const double B_mean = htrdr_planck(range_m[0], range_m[1], Tref);
+ *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0]));
+ *pdf *= 1.e-9; /* Transform from m^-1 to nm^-1 */
+ }
+
+ return lambda_m * 1.e9; /* Convert in nanometers */
+}
+
+static double
+wlen_ran_sample_discrete
+ (const struct htrdr_ran_wlen_planck* planck,
+ const double r0,
+ const double r1,
+ const char* func_name,
+ double* pdf)
+{
+ const double* cdf = NULL;
+ const double* find = NULL;
+ double r0_next = nextafter(r0, DBL_MAX);
+ double band_range[2];
+ double lambda = 0;
+ double pdf_continue = 0;
+ double pdf_band = 0;
+ size_t cdf_length = 0;
+ size_t i;
+ ASSERT(planck && planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE);
+ ASSERT(0 <= r0 && r0 < 1);
+ ASSERT(0 <= r1 && r1 < 1);
+ (void)func_name;
+ (void)pdf_band;
+
+ cdf = darray_double_cdata_get(&planck->cdf);
+ cdf_length = darray_double_size_get(&planck->cdf);
+ ASSERT(cdf_length > 0);
+
+ /* Use r_next rather than r0 in order to find the first entry that is not less
+ * than *or equal* to r0 */
+ find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ ASSERT(find);
+
+ i = (size_t)(find - cdf);
+ ASSERT(i < cdf_length && cdf[i] > r0 && (!i || cdf[i-1] <= r0));
+
+ band_range[0] = planck->range[0] + (double)i*planck->band_len;
+ band_range[1] = band_range[0] + planck->band_len;
+
+ /* Fetch the pdf of the sampled band */
+ pdf_band = darray_double_cdata_get(&planck->pdf)[i];
+
+ /* Uniformly sample a wavelength in the sampled band */
+ lambda = band_range[0] + (band_range[1] - band_range[0]) * r1;
+ pdf_continue = 1.0 / (band_range[1] - band_range[0]);
+
+ if(pdf) {
+ *pdf = pdf_band * pdf_continue;
+ }
+
+ return lambda;
+}
+
+static void
+release_wlen_planck_ran(ref_T* ref)
+{
+ struct htrdr_ran_wlen_planck* planck = NULL;
+ struct htrdr* htrdr = NULL;
+ ASSERT(ref);
+ planck = CONTAINER_OF(ref, struct htrdr_ran_wlen_planck, ref);
+ darray_double_release(&planck->cdf);
+ darray_double_release(&planck->pdf);
+ htrdr = planck->htrdr;
+ MEM_RM(htrdr_get_allocator(htrdr), planck);
+ htrdr_ref_put(planck->htrdr);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_ran_wlen_planck_create
+ (struct htrdr* htrdr,
+ /* range must be included in [200,1000] nm for shortwave or in [1000,100000]
+ * nanometers for longwave (thermal) */
+ const double range[2],
+ const size_t nbands, /* # bands used to discretized CDF */
+ const double ref_temperature,
+ struct htrdr_ran_wlen_planck** out_wlen_planck_ran)
+{
+ struct htrdr_ran_wlen_planck* planck = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && range && out_wlen_planck_ran && ref_temperature > 0);
+ ASSERT(ref_temperature > 0);
+ ASSERT(range[0] <= range[1]);
+
+ planck = MEM_CALLOC(htrdr->allocator, 1, sizeof(*planck));
+ if(!planck) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: could not allocate longwave random variate data structure -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ ref_init(&planck->ref);
+ darray_double_init(htrdr->allocator, &planck->cdf);
+ darray_double_init(htrdr->allocator, &planck->pdf);
+ htrdr_ref_get(htrdr);
+ planck->htrdr = htrdr;
+
+ planck->range[0] = range[0];
+ planck->range[1] = range[1];
+ planck->ref_temperature = ref_temperature;
+ planck->nbands = nbands;
+
+ if(nbands == HTRDR_WLEN_RAN_PLANCK_CONTINUE) {
+ planck->band_len = 0;
+ } else {
+ planck->band_len = (range[1] - range[0]) / (double)planck->nbands;
+
+ res = setup_wlen_planck_ran_cdf(planck, FUNC_NAME);
+ if(res != RES_OK) goto error;
+ }
+
+ htrdr_log(htrdr, "Spectral interval defined on [%g, %g] nanometers.\n",
+ range[0], range[1]);
+
+exit:
+ *out_wlen_planck_ran = planck;
+ return res;
+error:
+ if(planck) htrdr_ran_wlen_planck_ref_put(planck);
+ goto exit;
+}
+
+void
+htrdr_ran_wlen_planck_ref_get(struct htrdr_ran_wlen_planck* planck)
+{
+ ASSERT(planck);
+ ref_get(&planck->ref);
+}
+
+void
+htrdr_ran_wlen_planck_ref_put(struct htrdr_ran_wlen_planck* planck)
+{
+ ASSERT(planck);
+ ref_put(&planck->ref, release_wlen_planck_ran);
+}
+
+double
+htrdr_ran_wlen_planck_sample
+ (const struct htrdr_ran_wlen_planck* planck,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ ASSERT(planck);
+ if(planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE) { /* Discrete */
+ return wlen_ran_sample_discrete(planck, r0, r1, FUNC_NAME, pdf);
+ } else if(eq_eps(planck->range[0], planck->range[1], 1.e-6)) {
+ if(pdf) *pdf = 1;
+ return planck->range[0];
+ } else { /* Continue */
+ return wlen_ran_sample_continue
+ (planck, r0, planck->range, FUNC_NAME, pdf);
+ }
+}
+
diff --git a/src/core/htrdr_ran_wlen_planck.h b/src/core/htrdr_ran_wlen_planck.h
@@ -0,0 +1,61 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019 Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_RAN_WLEN_PLANCK_H
+#define HTRDR_RAN_WLEN_PLANCK_H
+
+#include "core/htrdr.h"
+#include <rsys/rsys.h>
+
+#define HTRDR_WLEN_RAN_PLANCK_CONTINUE 0
+
+/* Forward declarations */
+struct htrdr;
+struct htrdr_ran_wlen_planck;
+
+BEGIN_DECLS
+
+HTRDR_CORE_API res_T
+htrdr_ran_wlen_planck_create
+ (struct htrdr* htrdr,
+ const double range[2],
+ /* # bands used to discretisze the spectral domain. HTRDR_WLEN_RAN_CONTINUE
+ * <=> no discretisation */
+ const size_t nbands, /* Hint on #bands used to discretised th CDF */
+ const double ref_temperature, /* Reference temperature */
+ struct htrdr_ran_wlen_planck** planck);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_planck_ref_get
+ (struct htrdr_ran_wlen_planck* planck);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_planck_ref_put
+ (struct htrdr_ran_wlen_planck* planck);
+
+/* Return a wavelength in nanometer */
+HTRDR_CORE_API double
+htrdr_ran_wlen_planck_sample
+ (const struct htrdr_ran_wlen_planck* planck,
+ const double r0, /* Canonical number in [0, 1[ */
+ const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
+
+END_DECLS
+
+#endif /* HTRDR_RAN_WLEN_PLANCK_H */
+
diff --git a/src/planeto/htrdr_planeto.c b/src/planeto/htrdr_planeto.c
@@ -18,8 +18,8 @@
#define _POSIX_C_SOURCE 200112L /* fdopen, nextafter, rint */
#include "core/htrdr.h"
-#include "core/htrdr_cie_xyz.h"
-#include "core/htrdr_ran_wlen.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_ran_wlen_planck.h"
#include "core/htrdr_log.h"
#include "planeto/htrdr_planeto.h"
@@ -247,14 +247,14 @@ setup_spectral_domain
/* Planck distribution */
case HTRDR_SPECTRAL_LW:
case HTRDR_SPECTRAL_SW:
- res = htrdr_ran_wlen_create(cmd->htrdr, cmd->spectral_domain.wlen_range,
- nintervals, cmd->spectral_domain.ref_temperature,
- &cmd->ran_wlen_planck);
+ res = htrdr_ran_wlen_planck_create(cmd->htrdr,
+ cmd->spectral_domain.wlen_range, nintervals,
+ cmd->spectral_domain.ref_temperature, &cmd->planck);
break;
/* CIE XYZ distribution */
case HTRDR_SPECTRAL_SW_CIE_XYZ:
- res = htrdr_cie_xyz_create(cmd->htrdr, cmd->spectral_domain.wlen_range,
- nintervals, &cmd->ran_wlen_cie);
+ res = htrdr_ran_wlen_cie_xyz_create(cmd->htrdr,
+ cmd->spectral_domain.wlen_range, nintervals, &cmd->cie);
break;
default: FATAL("Unreachable code\n"); break;
@@ -443,8 +443,8 @@ planeto_release(ref_T* ref)
if(cmd->atmosphere) RNATM(ref_put(cmd->atmosphere));
if(cmd->ground) RNGRD(ref_put(cmd->ground));
if(cmd->source) htrdr_planeto_source_ref_put(cmd->source);
- if(cmd->ran_wlen_cie) htrdr_cie_xyz_ref_put(cmd->ran_wlen_cie);
- if(cmd->ran_wlen_planck) htrdr_ran_wlen_ref_put(cmd->ran_wlen_planck);
+ if(cmd->cie) htrdr_ran_wlen_cie_xyz_ref_put(cmd->cie);
+ if(cmd->planck) htrdr_ran_wlen_planck_ref_put(cmd->planck);
if(cmd->octrees_storage) CHK(fclose(cmd->octrees_storage) == 0);
if(cmd->output && cmd->output != stdout) CHK(fclose(cmd->output) == 0);
if(cmd->buf) htrdr_buffer_ref_put(cmd->buf);
diff --git a/src/planeto/htrdr_planeto_c.h b/src/planeto/htrdr_planeto_c.h
@@ -29,9 +29,9 @@
/* Forward declarations */
struct htrdr;
-struct htrdr_cie;
struct htrdr_pixel_format;
-struct htrdr_ran_wlen;
+struct htrdr_ran_wlen_cie_xyz;
+struct htrdr_ran_wlen_planck;
struct rnatm;
struct rngrd;
struct scam;
@@ -70,8 +70,8 @@ struct htrdr_planeto {
struct htrdr_planeto_source* source;
struct htrdr_args_spectral spectral_domain;
- struct htrdr_cie_xyz* ran_wlen_cie; /* HTRDR_SPECTRAL_SW_CIE_XYZ */
- struct htrdr_ran_wlen* ran_wlen_planck; /* HTRDR_SPECTRAL_<LW|SW> */
+ struct htrdr_ran_wlen_cie_xyz* cie; /* HTRDR_SPECTRAL_SW_CIE_XYZ */
+ struct htrdr_ran_wlen_planck* planck; /* HTRDR_SPECTRAL_<LW|SW> */
FILE* octrees_storage;
