commit 2d96123a255df15297e367eb1f221ae48e85da20
parent de1aab5971fc77f7d21c1b1892c80bbf350826f0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 27 May 2020 15:56:27 +0200
Rm the code used to track the spectral sampling bug
Diffstat:
| M | src/htrdr_cie_xyz.c | | | 230 | ------------------------------------------------------------------------------- |
1 file changed, 0 insertions(+), 230 deletions(-)
diff --git a/src/htrdr_cie_xyz.c b/src/htrdr_cie_xyz.c
@@ -29,8 +29,6 @@
#include <math.h> /* nextafter */
-/*#define SKY_BANDS*/
-
struct htrdr_cie_xyz {
struct darray_double cdf_X;
struct darray_double cdf_Y;
@@ -169,37 +167,6 @@ sample_cie_xyz
return lambda;
}
-static INLINE double
-sample_cie_xyz_from_sky_band
- (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[ */
-{
- double r0_next = nextafter(r0, DBL_MAX);
- double* find;
- double band_bounds[2];
- size_t i;
- size_t iband; /* Index of the sampled band */
- ASSERT(cie && cdf && cdf_length);
- ASSERT(0 <= r0 && r0 < 1);
- (void)f_bar;
-
- /* 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 */
- i = (size_t)(find - cdf);
- iband = htsky_get_spectral_band_id(cie->htrdr->sky, i);
- ASSERT(i< cdf_length);
- ASSERT(cdf[i] > r0 && (!iband || cdf[i-1] <= r0));
- HTSKY(get_spectral_band_bounds(cie->htrdr->sky, iband, band_bounds));
- return (band_bounds[0] + band_bounds[1]) * 0.5;
-}
-
static res_T
setup_cie_xyz
(struct htrdr_cie_xyz* cie,
@@ -289,172 +256,6 @@ error:
goto exit;
}
-static res_T
-setup_cie_xyz_from_sky_bands(struct htrdr_cie_xyz* cie, const char* func_name)
-{
- enum { X, Y, Z }; /* Helper constant */
- 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(cie && func_name);
- 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]);
-
- n = htsky_get_spectral_bands_count(cie->htrdr->sky);
- ASSERT(n);
-
- /* Allocate and reset the memory space for the tristimulus CDF */
- #define SETUP_STIMULUS(Stimulus) { \
- res = darray_double_resize(&cie->cdf_ ## Stimulus, n); \
- if(res != RES_OK) { \
- htrdr_log_err(cie->htrdr, \
- "%s: could not reserve the memory space for the CDF " \
- "of the "STR(X)" stimulus.\n", func_name); \
- goto error; \
- } \
- memset(darray_double_data_get(&cie->cdf_ ## Stimulus),0,n*sizeof(double)); \
- cdf[Stimulus] = darray_double_data_get(&cie->cdf_ ## Stimulus); \
- pdf[Stimulus] = cdf[Stimulus]; \
- memset(cdf[Stimulus], 0, n*sizeof(double)); \
- } (void)0
- SETUP_STIMULUS(X);
- SETUP_STIMULUS(Y);
- SETUP_STIMULUS(Z);
- #undef RESERVE
-
- /* Compute the *unormalized* pdf of the tristimulus */
- FOR_EACH(i, 0, n) {
- const size_t iband = htsky_get_spectral_band_id(cie->htrdr->sky, i);
- double band_bounds[2];
-
- HTSKY(get_spectral_band_bounds(cie->htrdr->sky, iband, band_bounds));
- ASSERT(band_bounds[0] < cie->range[1] && band_bounds[1] > cie->range[0]);
-
-
- ASSERT(band_bounds[0] <= band_bounds[1]);
- band_bounds[0] = MMAX(band_bounds[0], cie->range[0]);
- band_bounds[1] = MMIN(band_bounds[1], cie->range[1]);
-
- pdf[X][i] = trapezoidal_integration(band_bounds[0], band_bounds[1], fit_x_bar_1931);
- pdf[Y][i] = trapezoidal_integration(band_bounds[0], band_bounds[1], fit_y_bar_1931);
- pdf[Z][i] = trapezoidal_integration(band_bounds[0], band_bounds[1], 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, 0, 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(&cie->cdf_X);
- darray_double_clear(&cie->cdf_Y);
- darray_double_clear(&cie->cdf_Z);
- goto exit;
-}
-
-
-#include <star/ssp.h>
-#include <rsys/stretchy_array.h>
-
-static void
-test_integration(struct htrdr_cie_xyz* cie)
-{
- const size_t N = 500;
- const double dlambda = (cie->range[1] - cie->range[0]) / (double)N;
- double sum = 0;
- double ref = 0;
- size_t i;
- ASSERT(cie);
-
- FOR_EACH(i, 0, N) {
- const double lambda_lo = cie->range[0] + (double)i * dlambda;
- const double lambda_hi = lambda_lo + dlambda;
- sum += trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
- }
- ref = trapezoidal_integration(cie->range[0], cie->range[1], fit_x_bar_1931);
- printf("trapezoidal integration = %g; %g\n", ref, sum);
-}
-
-static void
-test_sample(struct htrdr_cie_xyz* cie)
-{
- const size_t N = 10000000;
- double* pdf = NULL;
- double* cdf = NULL;
- struct ssp_rng* rng = NULL;
- size_t i, n;
- ASSERT(cie);
-
- SSP(rng_create(cie->htrdr->allocator, &ssp_rng_mt19937_64, &rng));
- n = htsky_get_spectral_bands_count(cie->htrdr->sky);
- ASSERT(n);
-
- cdf = pdf = sa_add(pdf, n);
-
- FOR_EACH(i, 0, N) {
- const double r0 = ssp_rng_canonical(rng);
- const double r1 = ssp_rng_canonical(rng);
- const double lambda = htrdr_cie_xyz_sample_X(cie, r0, r1);
- const size_t iband = htsky_find_spectral_band(cie->htrdr->sky, lambda);
- const size_t j = iband - htsky_get_spectral_band_id(cie->htrdr->sky, 0);
- ASSERT(j < n);
- pdf[j] += 1;
- }
-
- FOR_EACH(i, 0, n) {
- pdf[i] /= (double)N;
- if(i == 0) {
- cdf[i] = pdf[i];
- } else {
- cdf[i] = pdf[i] + cdf[i-1];
- }
- printf("cdf[%lu] = %g\n", i, pdf[i]);
- }
-
- sa_release(pdf);
- SSP(rng_ref_put(rng));
-}
-
static void
release_cie_xyz(ref_T* ref)
{
@@ -510,22 +311,12 @@ htrdr_cie_xyz_create
printf("%lu\n", nbands);
}
-#ifdef SKY_BANDS
- (void)nbands;
- res = setup_cie_xyz_from_sky_bands(cie, FUNC_NAME);
-#else
res = setup_cie_xyz(cie, FUNC_NAME, nbands);
-#endif
if(res != RES_OK) goto error;
htrdr_log(htrdr, "Short wave interval defined on [%g, %g] nanometers.\n",
range[0], range[1]);
-#if 0
- test_integration(cie);
- test_sample(cie);
-#endif
-
exit:
*out_cie = cie;
return res;
@@ -554,15 +345,8 @@ htrdr_cie_xyz_sample_X
const double r0,
const double r1)
{
-#ifdef SKY_BANDS
- (void)r1;
- return sample_cie_xyz_from_sky_band
- (cie, darray_double_cdata_get(&cie->cdf_X),
- darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0);
-#else
return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
-#endif
}
double
@@ -571,15 +355,8 @@ htrdr_cie_xyz_sample_Y
const double r0,
const double r1)
{
-#ifdef SKY_BANDS
- (void)r1;
- return sample_cie_xyz_from_sky_band
- (cie, darray_double_cdata_get(&cie->cdf_Y),
- darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0);
-#else
return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
-#endif
}
double
@@ -588,14 +365,7 @@ htrdr_cie_xyz_sample_Z
const double r0,
const double r1)
{
-#ifdef SKY_BANDS
- (void)r1;
- return sample_cie_xyz_from_sky_band
- (cie, darray_double_cdata_get(&cie->cdf_Z),
- darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0);
-#else
return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
-#endif
}
