commit de1aab5971fc77f7d21c1b1892c80bbf350826f0
parent 1b56bcfd3e3951a89d3647ccd5defd9910a03c88
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 27 May 2020 15:53:42 +0200
Fix the spectral sampling in shortwave
Diffstat:
4 files changed, 300 insertions(+), 45 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -709,6 +709,39 @@ htrdr_fflush(struct htrdr* htrdr, FILE* stream)
/*******************************************************************************
* Local functions
******************************************************************************/
+double
+compute_sky_min_band_len
+ (struct htsky* sky,
+ const double range[2])
+{
+ double min_band_len = DBL_MAX;
+ size_t nbands;
+ ASSERT(sky && range && range[0] <= range[1]);
+
+ nbands = htsky_get_spectral_bands_count(sky);
+
+ if(eq_eps(range[0], range[1], 1.e-6)) {
+ ASSERT(nbands == 1);
+ min_band_len = 0;
+ } else {
+ size_t i = 0;
+
+ /* Compute the length of the current band clamped to the submitted range */
+ FOR_EACH(i, 0, nbands) {
+ const size_t iband = htsky_get_spectral_band_id(sky, i);
+ double wlens[2];
+ HTSKY(get_spectral_band_bounds(sky, iband, wlens));
+
+ /* Adjust band boundaries to the submitted range */
+ wlens[0] = MMAX(wlens[0], range[0]);
+ wlens[1] = MMIN(wlens[1], range[1]);
+
+ min_band_len = MMIN(wlens[1] - wlens[0], min_band_len);
+ }
+ }
+ return min_band_len;
+}
+
res_T
brightness_temperature
(struct htrdr* htrdr,
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -192,6 +192,13 @@ planck
}
}
+/* Return the minimum length in nanometer of the sky spectral bands
+ * clamped to in [range[0], range[1]]. */
+extern LOCAL_SYM double
+compute_sky_min_band_len
+ (struct htsky* sky,
+ const double range[2]);
+
extern LOCAL_SYM res_T
brightness_temperature
(struct htrdr* htrdr,
diff --git a/src/htrdr_cie_xyz.c b/src/htrdr_cie_xyz.c
@@ -19,6 +19,8 @@
#include "htrdr_c.h"
#include "htrdr_cie_xyz.h"
+#include <high_tune/htsky.h>
+
#include <rsys/algorithm.h>
#include <rsys/cstr.h>
#include <rsys/dynamic_array_double.h>
@@ -27,6 +29,8 @@
#include <math.h> /* nextafter */
+/*#define SKY_BANDS*/
+
struct htrdr_cie_xyz {
struct darray_double cdf_X;
struct darray_double cdf_Y;
@@ -58,7 +62,7 @@ trapezoidal_integration
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 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;
@@ -165,11 +169,41 @@ 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,
const char* func_name,
- const double range[2],
const size_t nbands)
{
enum { X, Y, Z }; /* Helper constant */
@@ -179,10 +213,10 @@ setup_cie_xyz
size_t i;
res_T res = RES_OK;
- ASSERT(cie && func_name && range && nbands);
- ASSERT(range[0] >= HTRDR_CIE_XYZ_RANGE_DEFAULT[0]);
- ASSERT(range[1] <= HTRDR_CIE_XYZ_RANGE_DEFAULT[1]);
- ASSERT(range[0] < range[1]);
+ 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) { \
@@ -203,13 +237,13 @@ setup_cie_xyz
#undef SETUP_STIMULUS
/* Compute the *unormalized* pdf of the tristimulus */
- cie->range[0] = range[0];
- cie->range[1] = range[1];
- cie->band_len = (range[1] - range[0]) / (double)nbands;
+
FOR_EACH(i, 0, nbands) {
- const double lambda_lo = range[0] + (double)i * cie->band_len;
- const double lambda_hi = lambda_lo + cie->band_len;
+ 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);
@@ -255,6 +289,172 @@ 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)
{
@@ -274,10 +474,12 @@ 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 discretisze the CIE tristimulus */
+ const size_t bands_count, /* # bands used to discretisze the CIE tristimulus */
struct htrdr_cie_xyz** out_cie)
{
struct htrdr_cie_xyz* cie = NULL;
+ double min_band_len = 0;
+ size_t nbands = bands_count;
res_T res = RES_OK;
ASSERT(htrdr && range && nbands && out_cie);
@@ -294,13 +496,36 @@ htrdr_cie_xyz_create
darray_double_init(htrdr->allocator, &cie->cdf_X);
darray_double_init(htrdr->allocator, &cie->cdf_Y);
darray_double_init(htrdr->allocator, &cie->cdf_Z);
+ cie->range[0] = range[0];
+ cie->range[1] = range[1];
+
+ min_band_len = compute_sky_min_band_len(cie->htrdr->sky, range);
+ cie->band_len = (range[1] - range[0]) / (double)nbands;
- res = setup_cie_xyz(cie, FUNC_NAME, range, nbands);
+ /* Adjust the band length to ensure that each sky spectral interval is
+ * overlapped by at least one band */
+ if(cie->band_len > min_band_len) {
+ cie->band_len = min_band_len;
+ nbands = (size_t)ceil((range[1] - range[0]) / cie->band_len);
+ 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;
@@ -329,8 +554,15 @@ 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
@@ -339,8 +571,15 @@ 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
@@ -349,7 +588,14 @@ 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
}
diff --git a/src/htrdr_ran_lw.c b/src/htrdr_ran_lw.c
@@ -122,37 +122,6 @@ error:
}
static double
-compute_sky_min_band_len(struct htrdr_ran_lw* ran_lw)
-{
- double min_band_len = DBL_MAX;
- size_t nbands;
- ASSERT(ran_lw && htsky_is_long_wave(ran_lw->htrdr->sky));
-
- nbands = htsky_get_spectral_bands_count(ran_lw->htrdr->sky);
-
- if(eq_eps(ran_lw->range[0], ran_lw->range[1], 1.e-6)) {
- ASSERT(nbands == 1);
- min_band_len = 0;
- } else {
- size_t i = 0;
-
- /* Compute the *unormalized* probability to sample a long wave band */
- FOR_EACH(i, 0, nbands) {
- const size_t iband = htsky_get_spectral_band_id(ran_lw->htrdr->sky, i);
- double wlens[2];
- HTSKY(get_spectral_band_bounds(ran_lw->htrdr->sky, iband, wlens));
-
- /* Adjust band boundaries to the submitted range */
- wlens[0] = MMAX(wlens[0], ran_lw->range[0]);
- wlens[1] = MMIN(wlens[1], ran_lw->range[1]);
-
- min_band_len = MMIN(wlens[1] - wlens[0], min_band_len);
- }
- }
- return min_band_len;
-}
-
-static double
ran_lw_sample_continue
(const struct htrdr_ran_lw* ran_lw,
const double r,
@@ -330,7 +299,7 @@ htrdr_ran_lw_create
ran_lw->ref_temperature = ref_temperature;
ran_lw->nbands = nbands;
- min_band_len = compute_sky_min_band_len(ran_lw);
+ min_band_len = compute_sky_min_band_len(ran_lw->htrdr->sky, ran_lw->range);
if(nbands == HTRDR_RAN_LW_CONTINUE) {
ran_lw->band_len = 0;
