commit 7ddf06d6bfa49fb5715422c86822c10f6cff4a62
parent 782c5a211aaf4a79c94bd5dee2935d328670daa3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 6 Aug 2018 18:06:40 +0200
Implement and test htgop_get_sw_spectral_intervals_CIE_XYZ
Diffstat:
5 files changed, 391 insertions(+), 286 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -42,7 +42,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(HTGOP_FILES_SRC
htgop.c
- htgop_sample_sw_spectral_interval_CIE_1931_XYZ.c)
+ htgop_sw_spectral_interval_CIE_XYZ.c)
set(HTGOP_FILES_INC
htgop_c.h
htgop_fetch_radiative_properties.h
diff --git a/src/htgop.h b/src/htgop.h
@@ -363,6 +363,12 @@ htgop_spectral_interval_sample_quadrature
const double r, /* Canonical random number in [0, 1[ */
size_t* iquad_point); /* Id of the sample quadrature point */
+/* Return the inclusive range of spectral interval that overlaps the CIE XYZ
+ * color space */
+HTGOP_API res_T
+htgop_get_sw_spectral_intervals_CIE_XYZ
+ (const struct htgop* htgop,
+ size_t specint_raneg[2]);
END_DECLS
diff --git a/src/htgop_sample_sw_spectral_interval_CIE_1931_XYZ.c b/src/htgop_sample_sw_spectral_interval_CIE_1931_XYZ.c
@@ -1,285 +0,0 @@
-/* Copyright (C) |Meso|Star> 2018 (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 support */
-
-#include "htgop.h"
-#include "htgop_c.h"
-
-#include <rsys/algorithm.h>
-
-/* Helper macros to convert from wave number to wavelength and vice versa */
-#define WNUM_TO_WLEN(Num) (1.0e7 / ((double)Num))
-#define WLEN_TO_WNUM(Len) WNUM_TO_WLEN(Len)
-
-/* Boundaries of the CIE XYZ color space in wavelength and wave number */
-#define CIE_XYZ_WAVELENGTH_MIN 380.0 /* In nanometer */
-#define CIE_XYZ_WAVELENGTH_MAX 780.0 /* In nanometer */
-#define CIE_XYZ_WAVENUMBER_MIN WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX) /*In cm^-1*/
-#define CIE_XYZ_WAVENUMBER_MAX WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN) /*In cm^-1*/
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static FINLINE 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_hi + 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 FINLINE 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 res_T
-sample_sw_spectral_interval_CIE
- (const struct htgop* htgop,
- const double* cdf,
- const size_t cdf_length,
- const double r, /* Canonical number in [0, 1[ */
- size_t* id)
-{
- double r_next = nextafter(r, DBL_MAX);
- double* find;
- size_t i;
- ASSERT(htgop && cdf && cdf_length);
-
- if(!id) return RES_BAD_ARG;
- if(r < 0 || r >= 1) {
- log_err(htgop, "Could nont sample a spectral interval in the CIE XYZ color "
- "space: invalid canonical random number `%g'.\n", r);
- return RES_BAD_ARG;
- }
-
- /* 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(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
- if(!find) { /* <=> CIE XYZ band does not intersect the loaded intervals */
- log_err(htgop, "Could not sample a spectral interval in the CIE XYZ color "
- "space : the loaded data does not overlap the CIE XYZ color space.\n");
- return RES_BAD_OP;
- }
-
- i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
- *id = i;
- return RES_OK;
-}
-
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-res_T
-htgop_sample_sw_spectral_interval_CIE_1931_X
- (const struct htgop* htgop,
- const double r, /* Canonical number in [0, 1[ */
- size_t* ispecint)
-{
- if(!htgop) return RES_BAD_ARG;
- return sample_sw_spectral_interval_CIE
- (htgop, darray_double_cdata_get(&htgop->sw_X_cdf),
- darray_double_size_get(&htgop->sw_X_cdf), r, ispecint);
-}
-
-res_T
-htgop_sample_sw_spectral_interval_CIE_1931_Y
- (const struct htgop* htgop,
- const double r, /* Canonical number in [0, 1[ */
- size_t* ispecint)
-{
- if(!htgop) return RES_BAD_ARG;
- return sample_sw_spectral_interval_CIE
- (htgop, darray_double_cdata_get(&htgop->sw_Y_cdf),
- darray_double_size_get(&htgop->sw_Y_cdf), r, ispecint);
-}
-
-res_T
-htgop_sample_sw_spectral_interval_CIE_1931_Z
- (const struct htgop* htgop,
- const double r, /* Canonical number in [0, 1[ */
- size_t* ispecint)
-{
- if(!htgop) return RES_BAD_ARG;
- return sample_sw_spectral_interval_CIE
- (htgop, darray_double_cdata_get(&htgop->sw_Z_cdf),
- darray_double_size_get(&htgop->sw_Z_cdf), r, ispecint);
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-compute_CIE_1931_XYZ_cdf(struct htgop* htgop)
-{
- enum { X, Y, Z }; /* Helper constant */
- size_t iband_min, iband_max;
- const double* wnums = NULL;
- double sum[3] = {0,0,0};
- double* pdf[3] = {NULL, NULL, NULL};
- double* cdf[3] = {NULL, NULL, NULL};
- size_t i, n;
- res_T res = RES_OK;
- ASSERT(htgop);
-
- HTGOP(get_sw_spectral_intervals_wave_numbers(htgop, &wnums));
- HTGOP(get_sw_spectral_intervals_count(htgop, &n));
- ASSERT(n);
-
- /* Allocate and reset the memory space for the tristimulus CDF */
- #define SETUP_STIMULUS(Stimulus) { \
- res = darray_double_resize(&htgop->sw_ ## Stimulus ## _cdf, n); \
- if(res != RES_OK) { \
- log_err(htgop, \
- "Could not reserve the memory space for the CDF " \
- "of the "STR(X)" stimulus.\n"); \
- goto error; \
- } \
- memset(darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf), 0, \
- n*sizeof(double)); \
- cdf[Stimulus] = darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf); \
- pdf[Stimulus] = cdf[Stimulus]; \
- } (void)0
- SETUP_STIMULUS(X);
- SETUP_STIMULUS(Y);
- SETUP_STIMULUS(Z);
- #undef RESERVE
-
- ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MAX, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN), 1.e-6));
- ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MIN, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX), 1.e-6));
- ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MAX, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MIN), 1.e-6));
- ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MIN, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MAX), 1.e-6));
-
- /* The sw spectral intervals does not intersect the CIE XYZ interval */
- if(wnums[0] > CIE_XYZ_WAVENUMBER_MAX || wnums[n] < CIE_XYZ_WAVENUMBER_MIN) {
- log_warn(htgop,
- "The loaded SW spectral intervals does not intersect "
- "the CIE XYZ interval.\n");
- goto exit;
- }
-
- /* Search the lower and upper bound of the bands that overlap the XYZ
- * spectral interval. Use a simple linear search since the overall number of
- * bands should be small */
- iband_min = 0, iband_max = n-1;
- FOR_EACH(i, 0, n) {
- if(wnums[i]<=CIE_XYZ_WAVENUMBER_MIN && CIE_XYZ_WAVENUMBER_MIN<=wnums[i+1])
- iband_min = i;
- if(wnums[i]<=CIE_XYZ_WAVENUMBER_MAX && CIE_XYZ_WAVENUMBER_MAX<=wnums[i+1])
- iband_max = i;
- }
- ASSERT(iband_min <= iband_max);
-
- /* Compute the *unormalized* pdf of the tristimulus */
- FOR_EACH(i, iband_min, iband_max+1) {
- const double lambda_lo = WNUM_TO_WLEN(wnums[i+1]);
- const double lambda_hi = WNUM_TO_WLEN(wnums[i+0]);
- ASSERT(lambda_lo <= lambda_hi);
- 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];
- }
-
- /* Compute the cdf of the tristimulus. Note that the upper bound is 'n' and not
- * 'iband_max+1' in order to ensure that the CDF is strictly increasing. */
- FOR_EACH(i, iband_min, n) {
- /* 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(eq_eps(cdf[X][n-1], 1, 1.e-6));
- ASSERT(eq_eps(cdf[Y][n-1], 1, 1.e-6));
- ASSERT(eq_eps(cdf[Z][n-1], 1, 1.e-6));
-
-#ifndef NDEBUG
- FOR_EACH(i, 1, n) {
- ASSERT(cdf[X][i] >= cdf[X][i-1]);
- ASSERT(cdf[Y][i] >= cdf[Y][i-1]);
- ASSERT(cdf[Z][i] >= cdf[Z][i-1]);
- }
-#endif
-
- /* Handle numerical issue */
- cdf[X][n-1] = 1.0;
- cdf[Y][n-1] = 1.0;
- cdf[Z][n-1] = 1.0;
-
-exit:
- return res;
-error:
- darray_double_clear(&htgop->sw_X_cdf);
- darray_double_clear(&htgop->sw_Y_cdf);
- darray_double_clear(&htgop->sw_Z_cdf);
- goto exit;
-}
-
-
diff --git a/src/htgop_sw_spectral_interval_CIE_XYZ.c b/src/htgop_sw_spectral_interval_CIE_XYZ.c
@@ -0,0 +1,336 @@
+/* Copyright (C) |Meso|Star> 2018 (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 support */
+
+#include "htgop.h"
+#include "htgop_c.h"
+
+#include <rsys/algorithm.h>
+
+/* Helper macros to convert from wave number to wavelength and vice versa */
+#define WNUM_TO_WLEN(Num) (1.0e7 / ((double)Num))
+#define WLEN_TO_WNUM(Len) WNUM_TO_WLEN(Len)
+
+/* Boundaries of the CIE XYZ color space in wavelength and wave number */
+#define CIE_XYZ_WAVELENGTH_MIN 380.0 /* In nanometer */
+#define CIE_XYZ_WAVELENGTH_MAX 780.0 /* In nanometer */
+#define CIE_XYZ_WAVENUMBER_MIN WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX) /*In cm^-1*/
+#define CIE_XYZ_WAVENUMBER_MAX WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN) /*In cm^-1*/
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static FINLINE 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_hi + 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 FINLINE 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 res_T
+sample_sw_spectral_interval_CIE
+ (const struct htgop* htgop,
+ const double* cdf,
+ const size_t cdf_length,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* id)
+{
+ double r_next = nextafter(r, DBL_MAX);
+ double* find;
+ size_t i;
+ ASSERT(htgop && cdf && cdf_length);
+
+ if(!id) return RES_BAD_ARG;
+ if(r < 0 || r >= 1) {
+ log_err(htgop, "Could nont sample a spectral interval in the CIE XYZ color "
+ "space: invalid canonical random number `%g'.\n", r);
+ return RES_BAD_ARG;
+ }
+
+ /* 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(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ if(!find) { /* <=> CIE XYZ band does not intersect the loaded intervals */
+ log_err(htgop, "Could not sample a spectral interval in the CIE XYZ color "
+ "space : the loaded data do not overlap the CIE XYZ color space.\n");
+ return RES_BAD_OP;
+ }
+
+ i = (size_t)(find - cdf);
+ ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
+ *id = i;
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+htgop_sample_sw_spectral_interval_CIE_1931_X
+ (const struct htgop* htgop,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* ispecint)
+{
+ if(!htgop) return RES_BAD_ARG;
+ return sample_sw_spectral_interval_CIE
+ (htgop, darray_double_cdata_get(&htgop->sw_X_cdf),
+ darray_double_size_get(&htgop->sw_X_cdf), r, ispecint);
+}
+
+res_T
+htgop_sample_sw_spectral_interval_CIE_1931_Y
+ (const struct htgop* htgop,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* ispecint)
+{
+ if(!htgop) return RES_BAD_ARG;
+ return sample_sw_spectral_interval_CIE
+ (htgop, darray_double_cdata_get(&htgop->sw_Y_cdf),
+ darray_double_size_get(&htgop->sw_Y_cdf), r, ispecint);
+}
+
+res_T
+htgop_sample_sw_spectral_interval_CIE_1931_Z
+ (const struct htgop* htgop,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* ispecint)
+{
+ if(!htgop) return RES_BAD_ARG;
+ return sample_sw_spectral_interval_CIE
+ (htgop, darray_double_cdata_get(&htgop->sw_Z_cdf),
+ darray_double_size_get(&htgop->sw_Z_cdf), r, ispecint);
+}
+
+res_T
+htgop_get_sw_spectral_intervals_CIE_XYZ
+ (const struct htgop* htgop,
+ size_t specint_range[2])
+{
+ const double* wnums = NULL;
+ const double wnum_min = nextafter(CIE_XYZ_WAVENUMBER_MIN, DBL_MAX);
+ const double wnum_max = CIE_XYZ_WAVENUMBER_MAX;
+ size_t nwnums;
+ res_T res = RES_OK;
+ double* low = NULL;
+ double* upp = NULL;
+ if(!htgop || !specint_range) return RES_BAD_ARG;
+
+ wnums = darray_double_cdata_get(&htgop->sw_specints.wave_numbers);
+ nwnums = darray_double_size_get(&htgop->sw_specints.wave_numbers);
+
+ low = search_lower_bound(&wnum_min, wnums, nwnums, sizeof(double), cmp_dbl);
+ upp = search_lower_bound(&wnum_max, wnums, nwnums, sizeof(double), cmp_dbl);
+
+ if(!low || upp == wnums) {
+ log_err(htgop,
+ "%s: the loaded data do not overlap the CIE XYZ color space.\n",
+ FUNC_NAME);
+ res = RES_BAD_OP;
+ goto error;
+ }
+ ASSERT(low[0] >= wnum_min && (low == wnums || low[-1] < wnum_min));
+ ASSERT(!upp || (upp[0] >= wnum_max && upp[-1] < wnum_max));
+
+ specint_range[0] = low == wnums ? 0 : (size_t)(low - wnums) - 1;
+ specint_range[1] = !upp ? nwnums - 1 : (size_t)(upp - wnums);
+ /* Transform the upper wnum bound in spectral interval id */
+ specint_range[1] = specint_range[1] - 1;
+
+ ASSERT(specint_range[0] < specint_range[1]);
+ ASSERT(wnums[specint_range[0]+0] <= CIE_XYZ_WAVENUMBER_MIN);
+ ASSERT(wnums[specint_range[0]+1] > CIE_XYZ_WAVENUMBER_MIN);
+ ASSERT(wnums[specint_range[1]+0] < CIE_XYZ_WAVENUMBER_MAX);
+ ASSERT(wnums[specint_range[1]+1] >= CIE_XYZ_WAVENUMBER_MAX);
+
+exit:
+ return res;
+error:
+ if(specint_range) {
+ specint_range[0] = SIZE_MAX;
+ specint_range[1] = 0;
+ }
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+compute_CIE_1931_XYZ_cdf(struct htgop* htgop)
+{
+ enum { X, Y, Z }; /* Helper constant */
+ size_t iband_min, iband_max;
+ const double* wnums = NULL;
+ double sum[3] = {0,0,0};
+ double* pdf[3] = {NULL, NULL, NULL};
+ double* cdf[3] = {NULL, NULL, NULL};
+ size_t i, n;
+ res_T res = RES_OK;
+ ASSERT(htgop);
+
+ HTGOP(get_sw_spectral_intervals_wave_numbers(htgop, &wnums));
+ HTGOP(get_sw_spectral_intervals_count(htgop, &n));
+ ASSERT(n);
+
+ /* Allocate and reset the memory space for the tristimulus CDF */
+ #define SETUP_STIMULUS(Stimulus) { \
+ res = darray_double_resize(&htgop->sw_ ## Stimulus ## _cdf, n); \
+ if(res != RES_OK) { \
+ log_err(htgop, \
+ "Could not reserve the memory space for the CDF " \
+ "of the "STR(X)" stimulus.\n"); \
+ goto error; \
+ } \
+ memset(darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf), 0, \
+ n*sizeof(double)); \
+ cdf[Stimulus] = darray_double_data_get(&htgop->sw_ ## Stimulus ## _cdf); \
+ pdf[Stimulus] = cdf[Stimulus]; \
+ } (void)0
+ SETUP_STIMULUS(X);
+ SETUP_STIMULUS(Y);
+ SETUP_STIMULUS(Z);
+ #undef RESERVE
+
+ ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MAX, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MIN), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVENUMBER_MIN, WLEN_TO_WNUM(CIE_XYZ_WAVELENGTH_MAX), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MAX, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MIN), 1.e-6));
+ ASSERT(eq_eps(CIE_XYZ_WAVELENGTH_MIN, WNUM_TO_WLEN(CIE_XYZ_WAVENUMBER_MAX), 1.e-6));
+
+ /* The sw spectral intervals does not intersect the CIE XYZ interval */
+ if(wnums[0] > CIE_XYZ_WAVENUMBER_MAX || wnums[n] < CIE_XYZ_WAVENUMBER_MIN) {
+ log_warn(htgop,
+ "The loaded SW spectral intervals does not intersect "
+ "the CIE XYZ interval.\n");
+ goto exit;
+ }
+
+ /* Search the lower and upper bound of the bands that overlap the XYZ
+ * spectral interval. Use a simple linear search since the overall number of
+ * bands should be small */
+ iband_min = 0, iband_max = n-1;
+ FOR_EACH(i, 0, n) {
+ if(wnums[i]<=CIE_XYZ_WAVENUMBER_MIN && CIE_XYZ_WAVENUMBER_MIN<=wnums[i+1])
+ iband_min = i;
+ if(wnums[i]<=CIE_XYZ_WAVENUMBER_MAX && CIE_XYZ_WAVENUMBER_MAX<=wnums[i+1])
+ iband_max = i;
+ }
+ ASSERT(iband_min <= iband_max);
+
+ /* Compute the *unormalized* pdf of the tristimulus */
+ FOR_EACH(i, iband_min, iband_max+1) {
+ const double lambda_lo = WNUM_TO_WLEN(wnums[i+1]);
+ const double lambda_hi = WNUM_TO_WLEN(wnums[i+0]);
+ ASSERT(lambda_lo <= lambda_hi);
+ 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];
+ }
+
+ /* Compute the cdf of the tristimulus. Note that the upper bound is 'n' and not
+ * 'iband_max+1' in order to ensure that the CDF is strictly increasing. */
+ FOR_EACH(i, iband_min, n) {
+ /* 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(eq_eps(cdf[X][n-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Y][n-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Z][n-1], 1, 1.e-6));
+
+#ifndef NDEBUG
+ FOR_EACH(i, 1, n) {
+ ASSERT(cdf[X][i] >= cdf[X][i-1]);
+ ASSERT(cdf[Y][i] >= cdf[Y][i-1]);
+ ASSERT(cdf[Z][i] >= cdf[Z][i-1]);
+ }
+#endif
+
+ /* Handle numerical issue */
+ cdf[X][n-1] = 1.0;
+ cdf[Y][n-1] = 1.0;
+ cdf[Z][n-1] = 1.0;
+
+exit:
+ return res;
+error:
+ darray_double_clear(&htgop->sw_X_cdf);
+ darray_double_clear(&htgop->sw_Y_cdf);
+ darray_double_clear(&htgop->sw_Z_cdf);
+ goto exit;
+}
+
+
diff --git a/src/test_htgop_load.c b/src/test_htgop_load.c
@@ -64,6 +64,53 @@ check_position_to_layer_id(const struct htgop* htgop)
}
}
+static void
+check_sw_specints_CIE_XYZ(const struct htgop* htgop)
+{
+ struct htgop_spectral_interval specint_low;
+ struct htgop_spectral_interval specint_upp;
+ const double wnum_min_CIE_XYZ = 1.e7/780.0;
+ const double wnum_max_CIE_XYZ = 1.e7/380.0;
+ size_t range[2];
+ size_t range2[2];
+ size_t i;
+ size_t nspecints;
+
+ CHK(htgop_get_sw_spectral_intervals_CIE_XYZ(NULL, range) == RES_BAD_ARG);
+ CHK(htgop_get_sw_spectral_intervals_CIE_XYZ(htgop, NULL) == RES_BAD_ARG);
+ CHK(htgop_get_sw_spectral_intervals_CIE_XYZ(htgop, range) == RES_OK);
+ CHK(range[0] <= range[1]);
+
+ CHK(htgop_get_sw_spectral_interval(htgop, range[0], &specint_low) == RES_OK);
+ CHK(htgop_get_sw_spectral_interval(htgop, range[1], &specint_upp) == RES_OK);
+
+ CHK(specint_low.wave_numbers[0] < specint_upp.wave_numbers[1]);
+ CHK(specint_low.wave_numbers[0] <= wnum_min_CIE_XYZ);
+ CHK(specint_upp.wave_numbers[1] >= wnum_max_CIE_XYZ);
+
+ CHK(htgop_get_sw_spectral_intervals_count(htgop, &nspecints) == RES_OK);
+ range2[0] = SIZE_MAX;
+ range2[1] = SIZE_MAX;
+ FOR_EACH(i, 0, nspecints) {
+ struct htgop_spectral_interval specint;
+ CHK(htgop_get_sw_spectral_interval(htgop, i, &specint) == RES_OK);
+
+ if(specint.wave_numbers[0] <= wnum_min_CIE_XYZ
+ && specint.wave_numbers[1] > wnum_min_CIE_XYZ) {
+ range2[0] = i;
+ }
+ if(specint.wave_numbers[0] < wnum_max_CIE_XYZ
+ && specint.wave_numbers[1] >= wnum_max_CIE_XYZ) {
+ range2[1] = i;
+ }
+ }
+ CHK(range2[0] != SIZE_MAX);
+ CHK(range2[1] != SIZE_MAX);
+ CHK(range2[0] <= range2[1]);
+ CHK(range2[0] == range[0]);
+ CHK(range2[1] == range[1]);
+}
+
int
main(int argc, char** argv)
{
@@ -359,6 +406,7 @@ main(int argc, char** argv)
}
check_position_to_layer_id(htgop);
+ check_sw_specints_CIE_XYZ(htgop);
CHK(fclose(fp) == 0);
reader_release(&rdr);
