commit 4abbb2e19283261047025c0b9e48586da2bd3fc3
parent 6d9ea2c72a6d274fcf20463c13ba77f9221c0800
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 15 Apr 2020 10:13:51 +0200
Update the sampling of the spectral dimension in short wave
First sample a wavelength wrt the CIE 1931 XYZ color space. Then define
its associated spectral band and finally sample a quadrature point into
it.
Diffstat:
11 files changed, 831 insertions(+), 368 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -90,8 +90,8 @@ set(HTRDR_FILES_SRC
htrdr.c
htrdr_args.c
htrdr_buffer.c
- htrdr_CIE_XYZ.c
htrdr_camera.c
+ htrdr_cie_xyz.c
htrdr_compute_radiance_sw.c
htrdr_compute_radiance_lw.c
htrdr_draw_radiance.c
@@ -107,8 +107,8 @@ set(HTRDR_FILES_INC
htrdr.h
htrdr_c.h
htrdr_buffer.h
- htrdr_CIE_XYZ.h
htrdr_camera.h
+ htrdr_cie_xyz.h
htrdr_grid.h
htrdr_ground.h
htrdr_interface.h
diff --git a/src/:w b/src/:w
@@ -0,0 +1,457 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019 CNRS, 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/>. */
+
+#include "htrdr.h"
+#include "htrdr_c.h"
+#include "htrdr_interface.h"
+#include "htrdr_ground.h"
+#include "htrdr_solve.h"
+#include "htrdr_sun.h"
+
+#include <high_tune/htsky.h>
+
+#include <star/s3d.h>
+#include <star/ssf.h>
+#include <star/ssp.h>
+#include <star/svx.h>
+
+#include <rsys/double2.h>
+#include <rsys/float2.h>
+#include <rsys/float3.h>
+
+struct scattering_context {
+ struct ssp_rng* rng;
+ const struct htsky* sky;
+ size_t iband; /* Index of the spectral band */
+ size_t iquad; /* Index of the quadrature point into the band */
+
+ double Ts; /* Sampled optical thickness */
+ double traversal_dst; /* Distance traversed along the ray */
+};
+static const struct scattering_context SCATTERING_CONTEXT_NULL = {
+ NULL, NULL, 0, 0, 0, 0
+};
+
+struct transmissivity_context {
+ struct ssp_rng* rng;
+ const struct htsky* sky;
+ size_t iband; /* Index of the spectrald */
+ size_t iquad; /* Index of the quadrature point into the band */
+
+ double Ts; /* Sampled optical thickness */
+ double Tmin; /* Minimal optical thickness */
+ double traversal_dst; /* Distance traversed along the ray */
+
+ enum htsky_property prop;
+};
+static const struct transmissivity_context TRANSMISSION_CONTEXT_NULL = {
+ NULL, NULL, 0, 0, 0, 0, 0, 0
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static int
+scattering_hit_filter
+ (const struct svx_hit* hit,
+ const double org[3],
+ const double dir[3],
+ const double range[2],
+ void* context)
+{
+ struct scattering_context* ctx = context;
+ double ks_max;
+ int pursue_traversal = 1;
+ ASSERT(hit && ctx && !SVX_HIT_NONE(hit) && org && dir && range);
+ (void)range;
+
+ ks_max = htsky_fetch_svx_voxel_property(ctx->sky, HTSKY_Ks,
+ HTSKY_SVX_MAX, HTSKY_CPNT_MASK_ALL, ctx->iband, ctx->iquad, &hit->voxel);
+
+ ctx->traversal_dst = hit->distance[0];
+
+ /* Iterate until a collision occurs into the voxel or until the ray
+ * does not collide the voxel */
+ for(;;) {
+ /* Compute tau for the current leaf */
+ const double vox_dst = hit->distance[1] - ctx->traversal_dst;
+ const double T = vox_dst * ks_max;
+
+ /* A collision occurs behind `vox_dst' */
+ if(ctx->Ts > T) {
+ ctx->Ts -= T;
+ ctx->traversal_dst = hit->distance[1];
+ pursue_traversal = 1;
+ break;
+
+ /* A real/null collision occurs before `vox_dst' */
+ } else {
+ double pos[3];
+ double proba;
+ double ks;
+ const double collision_dst = ctx->Ts / ks_max;
+
+ /* Compute the traversed distance up to the challenged collision */
+ ctx->traversal_dst += collision_dst;
+ ASSERT(ctx->traversal_dst >= hit->distance[0]);
+ ASSERT(ctx->traversal_dst <= hit->distance[1]);
+
+ /* Compute the world space position where a collision may occur */
+ pos[0] = org[0] + ctx->traversal_dst * dir[0];
+ pos[1] = org[1] + ctx->traversal_dst * dir[1];
+ pos[2] = org[2] + ctx->traversal_dst * dir[2];
+
+ ks = htsky_fetch_raw_property(ctx->sky, HTSKY_Ks,
+ HTSKY_CPNT_MASK_ALL, ctx->iband, ctx->iquad, pos, -DBL_MAX, DBL_MAX);
+
+ /* Handle the case that ks_max is not *really* the max */
+ proba = ks / ks_max;
+
+ if(ssp_rng_canonical(ctx->rng) < proba) {/* Collide <=> real scattering */
+ pursue_traversal = 0;
+ break;
+ } else { /* Null collision */
+ ctx->Ts = ssp_ran_exp(ctx->rng, 1); /* Sample a new optical thickness */
+ }
+ }
+ }
+ return pursue_traversal;
+}
+
+static int
+transmissivity_hit_filter
+ (const struct svx_hit* hit,
+ const double org[3],
+ const double dir[3],
+ const double range[2],
+ void* context)
+{
+ struct transmissivity_context* ctx = context;
+ int comp_mask = HTSKY_CPNT_MASK_ALL;
+ double k_max;
+ double k_min;
+ int pursue_traversal = 1;
+ ASSERT(hit && ctx && !SVX_HIT_NONE(hit) && org && dir && range);
+ (void)range;
+
+ k_min = htsky_fetch_svx_voxel_property(ctx->sky, ctx->prop,
+ HTSKY_SVX_MIN, comp_mask, ctx->iband, ctx->iquad, &hit->voxel);
+ k_max = htsky_fetch_svx_voxel_property(ctx->sky, ctx->prop,
+ HTSKY_SVX_MAX, comp_mask, ctx->iband, ctx->iquad, &hit->voxel);
+ ASSERT(k_min <= k_max);
+
+ ctx->Tmin += (hit->distance[1] - hit->distance[0]) * k_min;
+ ctx->traversal_dst = hit->distance[0];
+
+ /* Iterate until a collision occurs into the voxel or until the ray
+ * does not collide the voxel */
+ for(;;) {
+ const double vox_dst = hit->distance[1] - ctx->traversal_dst;
+ const double Tdif = vox_dst * (k_max-k_min);
+
+ /* A collision occurs behind `vox_dst' */
+ if(ctx->Ts > Tdif) {
+ ctx->Ts -= Tdif;
+ ctx->traversal_dst = hit->distance[1];
+ pursue_traversal = 1;
+ break;
+
+ /* A real/null collision occurs before `vox_dst' */
+ } else {
+ double x[3];
+ double k;
+ double proba;
+ double collision_dst = ctx->Ts / (k_max - k_min);
+
+ /* Compute the traversed distance up to the challenged collision */
+ ctx->traversal_dst += collision_dst;
+ ASSERT(ctx->traversal_dst >= hit->distance[0]);
+ ASSERT(ctx->traversal_dst <= hit->distance[1]);
+
+ /* Compute the world space position where a collision may occur */
+ x[0] = org[0] + ctx->traversal_dst * dir[0];
+ x[1] = org[1] + ctx->traversal_dst * dir[1];
+ x[2] = org[2] + ctx->traversal_dst * dir[2];
+
+ k = htsky_fetch_raw_property(ctx->sky, ctx->prop,
+ comp_mask, ctx->iband, ctx->iquad, x, k_min, k_max);
+ ASSERT(k >= k_min && k <= k_max);
+
+ proba = (k - k_min) / (k_max - k_min);
+
+ if(ssp_rng_canonical(ctx->rng) < proba) { /* Collide */
+ pursue_traversal = 0;
+ break;
+ } else { /* Null collision */
+ ctx->Ts = ssp_ran_exp(ctx->rng, 1); /* Sample a new optical thickness */
+ }
+ }
+ }
+ return pursue_traversal;
+}
+
+static double
+transmissivity
+ (struct htrdr* htrdr,
+ struct ssp_rng* rng,
+ const enum htsky_property prop,
+ const size_t iband,
+ const size_t iquad,
+ const double pos[3],
+ const double dir[3],
+ const double range[2])
+{
+ struct svx_hit svx_hit;
+ struct transmissivity_context transmissivity_ctx = TRANSMISSION_CONTEXT_NULL;
+
+ ASSERT(htrdr && rng && pos && dir && range);
+
+ transmissivity_ctx.rng = rng;
+ transmissivity_ctx.sky = htrdr->sky;
+ transmissivity_ctx.iband = iband;
+ transmissivity_ctx.iquad = iquad;
+ transmissivity_ctx.Ts = ssp_ran_exp(rng, 1); /* Sample an optical thickness */
+ transmissivity_ctx.prop = prop;
+
+ /* Compute the transmissivity */
+ HTSKY(trace_ray(htrdr->sky, pos, dir, range, NULL,
+ transmissivity_hit_filter, &transmissivity_ctx, iband, iquad, &svx_hit));
+
+ if(SVX_HIT_NONE(&svx_hit)) {
+ return transmissivity_ctx.Tmin ? exp(-transmissivity_ctx.Tmin) : 1.0;
+ } else {
+ return 0;
+ }
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+double
+htrdr_compute_radiance_sw
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ struct ssp_rng* rng,
+ const double pos_in[3],
+ const double dir_in[3],
+ const double wlen, /* In nanometer */
+ const size_t iband,
+ const size_t iquad)
+{
+ struct s3d_hit s3d_hit = S3D_HIT_NULL;
+ struct s3d_hit s3d_hit_tmp = S3D_HIT_NULL;
+ struct s3d_hit s3d_hit_prev = S3D_HIT_NULL;
+ struct svx_hit svx_hit = SVX_HIT_NULL;
+ struct ssf_phase* phase_hg = NULL;
+ struct ssf_phase* phase_rayleigh = NULL;
+
+ double pos[3];
+ double dir[3];
+ double range[2];
+ double pos_next[3];
+ double dir_next[3];
+ double band_bounds[2]; /* In nanometers */
+
+ double wlen;
+ double R;
+ double r; /* Random number */
+ double wo[3]; /* -dir */
+ double pdf;
+ double sun_solid_angle;
+ double Tr; /* Overall transmissivity */
+ double Tr_abs; /* Absorption transmissivity */
+ double L_sun; /* Sun radiance in W.m^-2.sr^-1 */
+ double sun_dir[3];
+ double ksi = 1; /* Throughput */
+ double w = 0; /* MC weight */
+ double g = 0; /* Asymmetry parameter of the HG phase function */
+
+ ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads);
+ ASSERT(!htsky_is_long_wave(htrdr->sky));
+
+ CHK(RES_OK == ssf_phase_create
+ (&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
+ CHK(RES_OK == ssf_phase_create
+ (&htrdr->lifo_allocators[ithread], &ssf_phase_rayleigh, &phase_rayleigh));
+
+ /* Setup the phase function for this spectral band & quadrature point */
+ g = htsky_fetch_particle_phase_function_asymmetry_parameter
+ (htrdr->sky, iband, iquad);
+ SSF(phase_hg_setup(phase_hg, g));
+
+ /* Fetch sun properties. Note that the sun spectral data are defined by bands
+ * that, actually are the same of the SW spectral bands defined in the
+ * default "ecrad_opt_prot.txt" file provided by the HTGOP project. */
+ htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
+ ASSERT(band_bounds[0] <= wlen && wlen <= band_bounds[1]);
+ sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
+ L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
+
+ d3_set(pos, pos_in);
+ d3_set(dir, dir_in);
+
+ /* Add the directly contribution of the sun */
+ if(htrdr_sun_is_dir_in_solar_cone(htrdr->sun, dir)) {
+ /* Add the direct contribution of the sun */
+ d2(range, 0, FLT_MAX);
+
+ /* Check that the ray is not occlude along the submitted range */
+ HTRDR(ground_trace_ray(htrdr->ground, pos, dir, range, NULL, &s3d_hit_tmp));
+ if(!S3D_HIT_NONE(&s3d_hit_tmp)) {
+ Tr = 0;
+ } else {
+ Tr = transmissivity
+ (htrdr, rng, HTSKY_Kext, iband, iquad , pos, dir, range);
+ w = L_sun * Tr;
+ }
+ }
+
+ /* Radiative random walk */
+ for(;;) {
+ struct scattering_context scattering_ctx = SCATTERING_CONTEXT_NULL;
+ double bounce_reflectivity = 1;
+
+ /* Find the first intersection with a surface */
+ d2(range, 0, DBL_MAX);
+ HTRDR(ground_trace_ray
+ (htrdr->ground, pos, dir, range, &s3d_hit_prev, &s3d_hit));
+
+ /* Sample an optical thickness */
+ scattering_ctx.Ts = ssp_ran_exp(rng, 1);
+
+ /* Setup the remaining scattering context fields */
+ scattering_ctx.rng = rng;
+ scattering_ctx.sky = htrdr->sky;
+ scattering_ctx.iband = iband;
+ scattering_ctx.iquad = iquad;
+
+ /* Define if a scattering event occurs */
+ d2(range, 0, s3d_hit.distance);
+ HTSKY(trace_ray(htrdr->sky, pos, dir, range, NULL,
+ scattering_hit_filter, &scattering_ctx, iband, iquad, &svx_hit));
+
+ /* No scattering and no surface reflection. Stop the radiative random walk */
+ if(S3D_HIT_NONE(&s3d_hit) && SVX_HIT_NONE(&svx_hit)) {
+ break;
+ }
+ ASSERT(SVX_HIT_NONE(&svx_hit)
+ || ( svx_hit.distance[0] <= scattering_ctx.traversal_dst
+ && svx_hit.distance[1] >= scattering_ctx.traversal_dst));
+
+ /* Negate the incoming dir to match the convention of the SSF library */
+ d3_minus(wo, dir);
+
+ /* Compute the new position */
+ pos_next[0] = pos[0] + dir[0]*scattering_ctx.traversal_dst;
+ pos_next[1] = pos[1] + dir[1]*scattering_ctx.traversal_dst;
+ pos_next[2] = pos[2] + dir[2]*scattering_ctx.traversal_dst;
+
+ /* Define the absorption transmissivity from the current position to the
+ * next position */
+ d2(range, 0, scattering_ctx.traversal_dst);
+ Tr_abs = transmissivity
+ (htrdr, rng, HTSKY_Ka, iband, iquad, pos, dir, range);
+ if(Tr_abs <= 0) break;
+
+ /* Sample a sun direction */
+ htrdr_sun_sample_direction(htrdr->sun, rng, sun_dir);
+ d2(range, 0, FLT_MAX);
+
+ s3d_hit_prev = SVX_HIT_NONE(&svx_hit) ? s3d_hit : S3D_HIT_NULL;
+
+ /* Check that the sun is visible from the new position */
+ HTRDR(ground_trace_ray
+ (htrdr->ground, pos_next, sun_dir, range, &s3d_hit_prev, &s3d_hit_tmp));
+
+ /* Compute the sun transmissivity */
+ if(!S3D_HIT_NONE(&s3d_hit_tmp)) {
+ Tr = 0;
+ } else {
+ Tr = transmissivity
+ (htrdr, rng, HTSKY_Kext, iband, iquad, pos_next, sun_dir, range);
+ }
+
+ /* Scattering at a surface */
+ if(SVX_HIT_NONE(&svx_hit)) {
+ struct htrdr_interface interf = HTRDR_INTERFACE_NULL;
+ struct ssf_bsdf* bsdf = NULL;
+ double N[3];
+ int type;
+
+ /* Fetch the hit interface and build its BSDF */
+ htrdr_ground_get_interface(htrdr->ground, &s3d_hit, &interf);
+ HTRDR(interface_create_bsdf
+ (htrdr, &interf, ithread, wlen, pos_next, dir, rng, &s3d_hit, &bsdf));
+
+ d3_normalize(N, d3_set_f3(N, s3d_hit.normal));
+ if(d3_dot(N, wo) < 0) d3_minus(N, N);
+
+ bounce_reflectivity = ssf_bsdf_sample
+ (bsdf, rng, wo, N, dir_next, &type, &pdf);
+ if(!(type & SSF_REFLECTION)) { /* Handle only reflections */
+ bounce_reflectivity = 0;
+ }
+
+ if(d3_dot(N, sun_dir) < 0) { /* Below the ground */
+ R = 0;
+ } else {
+ R = ssf_bsdf_eval(bsdf, wo, N, sun_dir) * d3_dot(N, sun_dir);
+ }
+
+ /* Release the BSDF */
+ SSF(bsdf_ref_put(bsdf));
+
+ /* Scattering in the medium */
+ } else {
+ struct ssf_phase* phase;
+ double ks_particle; /* Scattering coefficient of the particles */
+ double ks_gas; /* Scattering coefficient of the gaz */
+ double ks; /* Overall scattering coefficient */
+
+ ks_gas = htsky_fetch_raw_property(htrdr->sky, HTSKY_Ks,
+ HTSKY_CPNT_FLAG_GAS, iband, iquad, pos_next, -DBL_MAX, DBL_MAX);
+ ks_particle = htsky_fetch_raw_property(htrdr->sky, HTSKY_Ks,
+ HTSKY_CPNT_FLAG_PARTICLES, iband, iquad, pos_next, -DBL_MAX, DBL_MAX);
+ ks = ks_particle + ks_gas;
+
+ r = ssp_rng_canonical(rng);
+ if(r < ks_gas / ks) { /* Gas scattering */
+ phase = phase_rayleigh;
+ } else { /* Cloud scattering */
+ phase = phase_hg;
+ }
+
+ ssf_phase_sample(phase, rng, wo, dir_next, NULL);
+ R = ssf_phase_eval(phase, wo, sun_dir);
+ }
+
+ /* Update the MC weight */
+ ksi *= Tr_abs;
+ w += ksi * L_sun * sun_solid_angle * Tr * R;
+
+ /* Russian roulette wrt surface scattering */
+ if(SVX_HIT_NONE(&svx_hit) && ssp_rng_canonical(rng) >= bounce_reflectivity)
+ break;
+
+ /* Setup the next random walk state */
+ d3_set(pos, pos_next);
+ d3_set(dir, dir_next);
+ }
+ SSF(phase_ref_put(phase_hg));
+ SSF(phase_ref_put(phase_rayleigh));
+ return w;
+}
+
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -20,6 +20,7 @@
#include "htrdr_c.h"
#include "htrdr_args.h"
#include "htrdr_buffer.h"
+#include "htrdr_cie_xyz.h"
#include "htrdr_camera.h"
#include "htrdr_ground.h"
#include "htrdr_mtl.h"
@@ -583,6 +584,12 @@ htrdr_init
/* Define the CDF used to sample a long wave band */
res = setup_lw_cdf(htrdr);
if(res != RES_OK) goto error;
+ } else {
+ /* Setup the CDF used to sample the short wave */
+ res = htrdr_cie_xyz_create
+ (htrdr, HTRDR_CIE_XYZ_RANGE_DEFAULT, 400, &htrdr->cie);
+ if(res != RES_OK) goto error;
+
}
htrdr->lifo_allocators = MEM_CALLOC
@@ -648,6 +655,7 @@ htrdr_release(struct htrdr* htrdr)
if(htrdr->cam) htrdr_camera_ref_put(htrdr->cam);
if(htrdr->buf) htrdr_buffer_ref_put(htrdr->buf);
if(htrdr->mtl) htrdr_mtl_ref_put(htrdr->mtl);
+ if(htrdr->cie) htrdr_cie_xyz_ref_put(htrdr->cie);
if(htrdr->output && htrdr->output != stdout) fclose(htrdr->output);
if(htrdr->lifo_allocators) {
size_t i;
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -35,6 +35,7 @@
struct htsky;
struct htrdr_args;
struct htrdr_buffer;
+struct htrdr_cie_xyz;
struct htrdr_rectangle;
struct mem_allocator;
struct mutext;
@@ -49,6 +50,7 @@ struct htrdr {
struct htrdr_ground* ground;
struct htrdr_mtl* mtl;
struct htrdr_sun* sun;
+ struct htrdr_cie_xyz* cie;
struct htrdr_camera* cam;
struct htrdr_buffer* buf;
diff --git a/src/htrdr_CIE_XYZ.c b/src/htrdr_CIE_XYZ.c
@@ -1,283 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |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
-
-#include "htrdr.h"
-#include "htrdr_c.h"
-#include "htrdr_CIE_XYZ.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 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_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 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,
- const double r) /* Canonical number in [0, 1[ */
-{
- double r_next = nextafter(r, DBL_MAX);
- double* find;
- double wlen;
- size_t i;
- ASSERT(cie && cdf && cdf_length && r >= 0 && r < 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(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
- ASSERT(find);
-
- i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
-
- /* Return the wavelength at the center of the sampled band */
- wlen = (double)i * cie->band_len + 0.5*cie->band_len;
- ASSERT(cie->range[0] < wlen && wlen < cie->range[1]);
- return wlen;
-}
-
-
-static void
-release_CIE_XYZ(ref_T* ref)
-{
- struct htrdr_CIE_XYZ* cie = 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);
-}
-
-/*******************************************************************************
- * 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 nbands, /* # bands used to discretisze the CIE tristimulus */
- struct htrdr_CIE_XYZ** out_cie)
-{
-
- enum { X, Y, Z }; /* Helper constant */
- struct htrdr_CIE_XYZ* cie = NULL;
- 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(htrdr && range && nbands && out_cie);
- ASSERT(range[0] >= HTRDR_CIE_XYZ_RANGE_DEFAULT[0]);
- ASSERT(range[1] <= HTRDR_CIE_XYZ_RANGE_DEFAULT[1]);
- ASSERT(range[0] < range[1]);
-
- cie = MEM_CALLOC(htrdr->allocator, 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);
- cie->htrdr = htrdr;
- cie->range[0] = range[0];
- cie->range[1] = range[1];
- darray_double_init(htrdr->allocator, &cie->cdf_X);
- darray_double_init(htrdr->allocator, &cie->cdf_Y);
- darray_double_init(htrdr->allocator, &cie->cdf_Z);
-
- /* 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(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, 0, nbands*sizeof(double)); \
- } (void)0
- SETUP_STIMULUS(X);
- SETUP_STIMULUS(Y);
- SETUP_STIMULUS(Z);
- #undef RESERVE
-
- /* Compute the *unormalized* pdf of the tristimulus */
- 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;
- 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];
- }
-
- 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(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));
-
-#ifndef NDEBUG
- FOR_EACH(i, 1, nbands) {
- 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][nbands-1] = 1.0;
- cdf[Y][nbands-1] = 1.0;
- cdf[Z][nbands-1] = 1.0;
-
-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 r)
-{
- return sample_CIE_XYZ(cie, darray_double_cdata_get(&cie->cdf_X),
- darray_double_size_get(&cie->cdf_X), r);
-}
-
-double
-htrdr_CIE_XYZ_sample_Y(struct htrdr_CIE_XYZ* cie, const double r)
-{
- return sample_CIE_XYZ(cie, darray_double_cdata_get(&cie->cdf_Y),
- darray_double_size_get(&cie->cdf_Y), r);
-}
-
-double
-htrdr_CIE_XYZ_sample_Z(struct htrdr_CIE_XYZ* cie, const double r)
-{
- return sample_CIE_XYZ(cie, darray_double_cdata_get(&cie->cdf_Z),
- darray_double_size_get(&cie->cdf_Z), r);
-}
-
diff --git a/src/htrdr_CIE_XYZ.h b/src/htrdr_CIE_XYZ.h
@@ -1,61 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |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 <rsys/rsys.h>
-
-struct htrdr;
-struct htrdr_CIE_XYZ;
-
-/* Wavelength boundaries of the CIE XYZ color space in nanometers */
-static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] = {380, 780};
-
-extern LOCAL_SYM 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 s*/
- struct htrdr_CIE_XYZ** cie);
-
-extern LOCAL_SYM void
-htrdr_CIE_XYZ_ref_get
- (struct htrdr_CIE_XYZ* cie);
-
-extern LOCAL_SYM void
-htrdr_CIE_XYZ_ref_put
- (struct htrdr_CIE_XYZ* cie);
-
-/* Return a wavelength in nanometer */
-extern LOCAL_SYM double
-htrdr_CIE_XYZ_sample_X
- (struct htrdr_CIE_XYZ* cie,
- const double r); /* Canonical number in [0, 1[ */
-
-/* Return a wavelength in nanometer */
-extern LOCAL_SYM double
-htrdr_CIE_XYZ_sample_Y
- (struct htrdr_CIE_XYZ* cie,
- const double r); /* Canonical number in [0, 1[ */
-
-/* Return a wavelength in nanometer */
-extern LOCAL_SYM double
-htrdr_CIE_XYZ_sample_Z
- (struct htrdr_CIE_XYZ* cie,
- const double r); /* Canonical number in [0, 1[ */
-
-#endif /* HTRDR_CIE_XYZ_H */
-
diff --git a/src/htrdr_cie_xyz.c b/src/htrdr_cie_xyz.c
@@ -0,0 +1,284 @@
+/* Copyright (C) 2018, 2019, 2020 |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
+
+#include "htrdr.h"
+#include "htrdr_c.h"
+#include "htrdr_cie_xyz.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 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_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 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,
+ const double r) /* Canonical number in [0, 1[ */
+{
+ double r_next = nextafter(r, DBL_MAX);
+ double* find;
+ double wlen;
+ size_t i;
+ ASSERT(cie && cdf && cdf_length && r >= 0 && r < 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(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ ASSERT(find);
+
+ i = (size_t)(find - cdf);
+ ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
+
+ /* Return the wavelength at the center of the sampled band */
+ wlen = cie->range[0] + cie->band_len * ((double)i + 0.5);
+ ASSERT(cie->range[0] < wlen && wlen < cie->range[1]);
+ return wlen;
+}
+
+
+static void
+release_cie_xyz(ref_T* ref)
+{
+ struct htrdr_cie_xyz* cie = 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);
+ MEM_RM(cie->htrdr->allocator, cie);
+}
+
+/*******************************************************************************
+ * 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 nbands, /* # bands used to discretisze the CIE tristimulus */
+ struct htrdr_cie_xyz** out_cie)
+{
+
+ enum { X, Y, Z }; /* Helper constant */
+ struct htrdr_cie_xyz* cie = NULL;
+ 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(htrdr && range && nbands && out_cie);
+ ASSERT(range[0] >= HTRDR_CIE_XYZ_RANGE_DEFAULT[0]);
+ ASSERT(range[1] <= HTRDR_CIE_XYZ_RANGE_DEFAULT[1]);
+ ASSERT(range[0] < range[1]);
+
+ cie = MEM_CALLOC(htrdr->allocator, 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);
+ cie->htrdr = htrdr;
+ cie->range[0] = range[0];
+ cie->range[1] = range[1];
+ darray_double_init(htrdr->allocator, &cie->cdf_X);
+ darray_double_init(htrdr->allocator, &cie->cdf_Y);
+ darray_double_init(htrdr->allocator, &cie->cdf_Z);
+
+ /* 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(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 */
+ 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;
+ 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];
+ }
+
+ 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(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));
+
+#ifndef NDEBUG
+ FOR_EACH(i, 1, nbands) {
+ 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][nbands-1] = 1.0;
+ cdf[Y][nbands-1] = 1.0;
+ cdf[Z][nbands-1] = 1.0;
+
+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 r)
+{
+ return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
+ darray_double_size_get(&cie->cdf_X), r);
+}
+
+double
+htrdr_cie_xyz_sample_Y(struct htrdr_cie_xyz* cie, const double r)
+{
+ return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
+ darray_double_size_get(&cie->cdf_Y), r);
+}
+
+double
+htrdr_cie_xyz_sample_Z(struct htrdr_cie_xyz* cie, const double r)
+{
+ return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
+ darray_double_size_get(&cie->cdf_Z), r);
+}
+
diff --git a/src/htrdr_cie_xyz.h b/src/htrdr_cie_xyz.h
@@ -0,0 +1,61 @@
+/* Copyright (C) 2018, 2019, 2020 |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 <rsys/rsys.h>
+
+struct htrdr;
+struct htrdr_cie_xyz;
+
+/* Wavelength boundaries of the CIE XYZ color space in nanometers */
+static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] = {380, 780};
+
+extern LOCAL_SYM 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 s*/
+ struct htrdr_cie_xyz** cie);
+
+extern LOCAL_SYM void
+htrdr_cie_xyz_ref_get
+ (struct htrdr_cie_xyz* cie);
+
+extern LOCAL_SYM void
+htrdr_cie_xyz_ref_put
+ (struct htrdr_cie_xyz* cie);
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_cie_xyz_sample_X
+ (struct htrdr_cie_xyz* cie,
+ const double r); /* Canonical number in [0, 1[ */
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_cie_xyz_sample_Y
+ (struct htrdr_cie_xyz* cie,
+ const double r); /* Canonical number in [0, 1[ */
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_cie_xyz_sample_Z
+ (struct htrdr_cie_xyz* cie,
+ const double r); /* Canonical number in [0, 1[ */
+
+#endif /* HTRDR_cie_xyz_H */
+
diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c
@@ -247,6 +247,7 @@ htrdr_compute_radiance_sw
struct ssp_rng* rng,
const double pos_in[3],
const double dir_in[3],
+ const double wlen, /* In nanometer */
const size_t iband,
const size_t iquad)
{
@@ -264,7 +265,6 @@ htrdr_compute_radiance_sw
double dir_next[3];
double band_bounds[2]; /* In nanometers */
- double wlen;
double R;
double r; /* Random number */
double wo[3]; /* -dir */
@@ -291,13 +291,11 @@ htrdr_compute_radiance_sw
(htrdr->sky, iband, iquad);
SSF(phase_hg_setup(phase_hg, g));
- /* Fetch sun properties. Arbitrarily use the wavelength at the center of the
- * band to retrieve the sun radiance of the current band. Note that the sun
- * spectral data are defined by bands that, actually are the same of the SW
- * spectral bands defined in the default "ecrad_opt_prot.txt" file provided
- * by the HTGOP project. */
+ /* Fetch sun properties. Note that the sun spectral data are defined by bands
+ * that, actually are the same of the SW spectral bands defined in the
+ * default "ecrad_opt_prot.txt" file provided by the HTGOP project. */
htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
- wlen = (band_bounds[0] + band_bounds[1]) * 0.5;
+ ASSERT(band_bounds[0] <= wlen && wlen <= band_bounds[1]);
sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -20,6 +20,7 @@
#include "htrdr_c.h"
#include "htrdr_buffer.h"
#include "htrdr_camera.h"
+#include "htrdr_cie_xyz.h"
#include "htrdr_solve.h"
#include <high_tune/htsky.h>
@@ -541,8 +542,9 @@ draw_pixel_sw
double ray_dir[3];
double weight;
double r0, r1;
- size_t iband;
- size_t iquad;
+ double wlen; /* Sampled wavelength into the spectral band */
+ size_t iband; /* Sampled spectral band */
+ size_t iquad; /* Sampled quadrature point into the spectral band */
double usec;
/* Begin the registration of the time spent to in the realisation */
@@ -561,24 +563,18 @@ draw_pixel_sw
/* Sample a spectral band and a quadrature point */
switch(ichannel) {
- case 0:
- htsky_sample_sw_spectral_data_CIE_1931_X
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
- case 1:
- htsky_sample_sw_spectral_data_CIE_1931_Y
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
- case 2:
- htsky_sample_sw_spectral_data_CIE_1931_Z
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
+ case 0: wlen = htrdr_cie_xyz_sample_X(htrdr->cie, r0); break;
+ case 1: wlen = htrdr_cie_xyz_sample_Y(htrdr->cie, r0); break;
+ case 2: wlen = htrdr_cie_xyz_sample_Z(htrdr->cie, r0); break;
default: FATAL("Unreachable code.\n"); break;
}
+ iband = htsky_find_spectral_band(htrdr->sky, wlen);
+ iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r1, iband);
+
/* Compute the luminance */
weight = htrdr_compute_radiance_sw
- (htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
+ (htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
ASSERT(weight >= 0);
/* End the registration of the per realisation time */
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -63,6 +63,7 @@ htrdr_compute_radiance_sw
struct ssp_rng* rng,
const double pos[3],
const double dir[3],
+ const double wlen,
const size_t iband, /* Index of the spectral band */
const size_t iquad); /* Index of the quadrature point into the band */
