commit 90b51beaeb6807ef7958ec656c08b808e61dd633
parent 3f539eecb55a386a6de92a68178d7009fba9d01d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 20 Apr 2020 14:54:47 +0200
Add support of continue LW/SW sampling
Diffstat:
6 files changed, 313 insertions(+), 553 deletions(-)
diff --git a/src/:w b/src/:w
@@ -1,457 +0,0 @@
-/* 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_cie_xyz.c b/src/htrdr_cie_xyz.c
@@ -100,26 +100,66 @@ 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 (*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 r_next = nextafter(r, DBL_MAX);
+ double r0_next = nextafter(r0, DBL_MAX);
double* find;
- double wlen;
- size_t i;
- ASSERT(cie && cdf && cdf_length && r >= 0 && r < 1);
+ 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(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ 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] > r && (!i || cdf[i-1] <= r));
+ /* 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);
+ 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 {
+ FATAL("Unexpected error.\n");
+ }
- /* 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;
+ return lambda;
}
static res_T
@@ -278,23 +318,32 @@ htrdr_cie_xyz_ref_put(struct htrdr_cie_xyz* cie)
}
double
-htrdr_cie_xyz_sample_X(struct htrdr_cie_xyz* cie, const double r)
+htrdr_cie_xyz_sample_X
+ (struct htrdr_cie_xyz* cie,
+ const double r0,
+ const double r1)
{
return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
- darray_double_size_get(&cie->cdf_X), r);
+ darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
}
double
-htrdr_cie_xyz_sample_Y(struct htrdr_cie_xyz* cie, const double r)
+htrdr_cie_xyz_sample_Y
+ (struct htrdr_cie_xyz* cie,
+ const double r0,
+ const double r1)
{
return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
- darray_double_size_get(&cie->cdf_Y), r);
+ darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
}
double
-htrdr_cie_xyz_sample_Z(struct htrdr_cie_xyz* cie, const double r)
+htrdr_cie_xyz_sample_Z
+ (struct htrdr_cie_xyz* cie,
+ const double r0,
+ const double r1)
{
return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
- darray_double_size_get(&cie->cdf_Z), r);
+ darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
}
diff --git a/src/htrdr_cie_xyz.h b/src/htrdr_cie_xyz.h
@@ -43,19 +43,19 @@ htrdr_cie_xyz_ref_put
extern LOCAL_SYM double
htrdr_cie_xyz_sample_X
(struct htrdr_cie_xyz* cie,
- const double r); /* Canonical number in [0, 1[ */
+ const double r0, const double r1); /* Canonical numbers 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[ */
+ const double r0, const double r1); /* 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[ */
+ const double r0, const double r1); /* Canonical number in [0, 1[ */
#endif /* HTRDR_cie_xyz_H */
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -512,7 +512,7 @@ draw_pixel_sw
double ray_org[3];
double ray_dir[3];
double weight;
- double r0, r1;
+ double r0, r1, r2;
double wlen; /* Sampled wavelength into the spectral band */
size_t iband; /* Sampled spectral band */
size_t iquad; /* Sampled quadrature point into the spectral band */
@@ -531,17 +531,18 @@ draw_pixel_sw
r0 = ssp_rng_canonical(rng);
r1 = ssp_rng_canonical(rng);
+ r2 = ssp_rng_canonical(rng);
/* Sample a spectral band and a quadrature point */
switch(ichannel) {
- 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;
+ case 0: wlen = htrdr_cie_xyz_sample_X(htrdr->cie, r0, r1); break;
+ case 1: wlen = htrdr_cie_xyz_sample_Y(htrdr->cie, r0, r1); break;
+ case 2: wlen = htrdr_cie_xyz_sample_Z(htrdr->cie, r0, r1); 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);
+ iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r2, iband);
/* Compute the luminance */
weight = htrdr_compute_radiance_sw
@@ -604,9 +605,9 @@ draw_pixel_lw
size_t iband;
size_t iquad;
double usec;
- double pdf;
+ double band_pdf;
- /* Begin the registration of the time spent to in the realisation */
+ /* Begin the registration of the time spent in the realisation */
time_current(&t0);
/* Sample a position into the pixel, in the normalized image plane */
@@ -621,16 +622,19 @@ draw_pixel_lw
r1 = ssp_rng_canonical(rng);
/* Sample a wavelength */
- wlen = htrdr_ran_lw_sample(htrdr->ran_lw, r0, &pdf);
+ wlen = htrdr_ran_lw_sample(htrdr->ran_lw, r0);
/* Select the associated band and sample a quadrature point */
iband = htsky_find_spectral_band(htrdr->sky, wlen);
iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r1, iband);
+ /* Retrieve the PDF to sample this sky band */
+ band_pdf = htrdr_ran_lw_get_sky_band_pdf(htrdr->ran_lw, iband);
+
/* Compute the luminance */
weight = htrdr_compute_radiance_lw
(htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
- weight /= pdf;
+ weight /= band_pdf;
ASSERT(weight >= 0);
/* End the registration of the per realisation time */
diff --git a/src/htrdr_ran_lw.c b/src/htrdr_ran_lw.c
@@ -19,6 +19,8 @@
#include "htrdr_c.h"
#include "htrdr_ran_lw.h"
+#include <high_tune/htsky.h>
+
#include <rsys/algorithm.h>
#include <rsys/cstr.h>
#include <rsys/dynamic_array_double.h>
@@ -29,9 +31,12 @@
struct htrdr_ran_lw {
struct darray_double cdf;
- struct darray_double pdf;
+ /* Probas to sample a sky band overlapped by the range */
+ struct darray_double sky_bands_pdf;
double range[2]; /* Boundaries of the handled CIE XYZ color space */
double band_len; /* Length in nanometers of a band */
+ double ref_temperature; /* In Kelvin */
+ size_t nbands; /* # bands */
struct htrdr* htrdr;
ref_T ref;
@@ -41,24 +46,18 @@ struct htrdr_ran_lw {
* Helper functions
******************************************************************************/
static res_T
-setup_ran_lw
+setup_ran_lw_cdf
(struct htrdr_ran_lw* ran_lw,
- const char* func_name,
- const double range[2],
- const size_t nbands,
- const double ref_temperature)
+ const char* func_name)
{
double* pdf = NULL;
double* cdf = NULL;
double sum = 0;
size_t i;
res_T res = RES_OK;
- ASSERT(ran_lw && func_name && range && range && ref_temperature);
- ASSERT(range[0] >= 1000);
- ASSERT(range[1] <= 100000);
- ASSERT(range[0] < range[1]);
+ ASSERT(ran_lw && func_name && ran_lw->nbands != HTRDR_RAN_LW_CONTINUE);
- res = darray_double_resize(&ran_lw->cdf, nbands);
+ res = darray_double_resize(&ran_lw->cdf, ran_lw->nbands);
if(res != RES_OK) {
htrdr_log_err(ran_lw->htrdr,
"%s: Error allocating the CDF of the long wave spectral bands -- %s.\n",
@@ -66,23 +65,12 @@ setup_ran_lw
goto error;
}
- res = darray_double_resize(&ran_lw->pdf, nbands);
- if(res != RES_OK) {
- htrdr_log_err(ran_lw->htrdr,
- "%s: Error allocating the PDF of the long wave spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
-
cdf = darray_double_data_get(&ran_lw->cdf);
- pdf = darray_double_data_get(&ran_lw->pdf);
+ pdf = cdf; /* Alias the CDF by the PDF since only one array is necessary */
/* Compute the *unormalized* probability to sample a long wave band */
- ran_lw->range[0] = range[0];
- ran_lw->range[1] = range[1];
- ran_lw->band_len = (range[1] - range[0]) / (double)nbands;
- FOR_EACH(i, 0, nbands) {
- double lambda_lo = range[0] + (double)i * ran_lw->band_len;
+ FOR_EACH(i, 0, ran_lw->nbands) {
+ double lambda_lo = ran_lw->range[0] + (double)i * ran_lw->band_len;
double lambda_hi = lambda_lo + ran_lw->band_len;
ASSERT(lambda_lo <= lambda_hi);
@@ -91,14 +79,14 @@ setup_ran_lw
lambda_hi *= 1.e-9;
/* Compute the probability of the current band */
- pdf[i] = planck(lambda_lo, lambda_hi, ref_temperature);
+ pdf[i] = planck(lambda_lo, lambda_hi, ran_lw->ref_temperature);
/* Update the norm */
sum += pdf[i];
}
/* Compute the cumulative of the previously computed probabilities */
- FOR_EACH(i, 0, nbands) {
+ FOR_EACH(i, 0, ran_lw->nbands) {
/* Normalize the probability */
pdf[i] /= sum;
@@ -110,17 +98,170 @@ setup_ran_lw
ASSERT(cdf[i] >= cdf[i-1]);
}
}
- ASSERT(eq_eps(cdf[nbands-1], 1, 1.e-6));
- cdf[nbands - 1] = 1.0; /* Handle numerical issue */
+ ASSERT(eq_eps(cdf[ran_lw->nbands-1], 1, 1.e-6));
+ cdf[ran_lw->nbands - 1] = 1.0; /* Handle numerical issue */
exit:
return res;
error:
darray_double_clear(&ran_lw->cdf);
- darray_double_clear(&ran_lw->pdf);
goto exit;
}
+static res_T
+compute_sky_bands_pdf(struct htrdr_ran_lw* ran_lw, const char* func_name)
+{
+ double* pdf = NULL;
+ double sum = 0;
+ size_t i;
+ size_t nbands;
+ res_T res = RES_OK;
+ ASSERT(ran_lw && htsky_is_long_wave(ran_lw->htrdr->sky) && func_name);
+
+ nbands = htsky_get_spectral_bands_count(ran_lw->htrdr->sky);
+
+ res = darray_double_resize(&ran_lw->sky_bands_pdf, nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(ran_lw->htrdr,
+ "%s: error allocating the PDF of the long wave spectral bands "
+ "of the sky -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+
+ pdf = darray_double_data_get(&ran_lw->sky_bands_pdf);
+
+ /* Compute the *unormalized* probability to sample a long wave band */
+ sum = 0;
+ 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));
+
+ /* Convert from nanometer to meter */
+ wlens[0] = wlens[0] * 1.e-9;
+ wlens[1] = wlens[1] * 1.e-9;
+
+ /* Compute the probability of the current band */
+ pdf[i] = planck(wlens[0], wlens[1], ran_lw->ref_temperature);
+
+ /* Update the norm */
+ sum += pdf[i];
+ }
+
+ /* Normalize the probabilities */
+ FOR_EACH(i, 0, nbands) pdf[i] /= sum;
+
+exit:
+ return res;
+error:
+ darray_double_clear(&ran_lw->sky_bands_pdf);
+ goto exit;
+
+}
+
+static double
+ran_lw_sample_discreet
+ (const struct htrdr_ran_lw* ran_lw,
+ const double r,
+ const char* func_name)
+{
+ const double* cdf = NULL;
+ const double* find = NULL;
+ double r_next = nextafter(r, DBL_MAX);
+ double wlen = 0;
+ size_t cdf_length = 0;
+ size_t i;
+ ASSERT(ran_lw && ran_lw->nbands != HTRDR_RAN_LW_CONTINUE);
+ ASSERT(0 <= r && r < 1);
+ (void)func_name;
+
+ cdf = darray_double_cdata_get(&ran_lw->cdf);
+ cdf_length = darray_double_size_get(&ran_lw->cdf);
+ ASSERT(cdf_length > 0);
+
+ /* 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 = ran_lw->range[0] + ran_lw->band_len * ((double)i + 0.5);
+ ASSERT(ran_lw->range[0] < wlen && wlen < ran_lw->range[1]);
+
+ return wlen;
+}
+
+static double
+ran_lw_sample_continue
+ (const struct htrdr_ran_lw* ran_lw,
+ const double r,
+ const char* func_name)
+{
+ /* Numerical parameters of the dichotomy search */
+ const size_t MAX_ITER = 100;
+ const double EPSILON_LAMBDA = 1e-6;
+ 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(ran_lw && ran_lw->nbands == HTRDR_RAN_LW_CONTINUE && func_name);
+ ASSERT(0 <= r && r < 1);
+
+ /* Convert the wavelength range in meters */
+ range_m[0] = ran_lw->range[0] / 1.0e9;
+ range_m[1] = ran_lw->range[1] / 1.0e9;
+
+ /* Setup the dichotomy search */
+ lambda_m_min = range_m[0];
+ lambda_m_max = range_m[1];
+ bf_min_max = blackbody_fraction
+ (range_m[0], range_m[1], ran_lw->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 = blackbody_fraction(range_m[0], lambda_m, ran_lw->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
+ || fabs(bf_prev - bf) < EPSILON_BF)
+ break;
+
+ lambda_m_prev = lambda_m;
+ bf_prev = bf;
+ }
+ if(i >= MAX_ITER) {
+ htrdr_log_warn(ran_lw->htrdr,
+ "%s: could not sample a wavelength in the range [%g, %g] nanometers "
+ "for the reference temperature %g Kelvin.\n",
+ func_name, SPLIT2(ran_lw->range), ran_lw->ref_temperature);
+ }
+
+ return lambda_m * 1.0e9; /* Convert in nanometers */
+}
+
static void
release_ran_lw(ref_T* ref)
{
@@ -128,7 +269,7 @@ release_ran_lw(ref_T* ref)
ASSERT(ref);
ran_lw = CONTAINER_OF(ref, struct htrdr_ran_lw, ref);
darray_double_release(&ran_lw->cdf);
- darray_double_release(&ran_lw->pdf);
+ darray_double_release(&ran_lw->sky_bands_pdf);
MEM_RM(ran_lw->htrdr->allocator, ran_lw);
}
@@ -139,13 +280,17 @@ res_T
htrdr_ran_lw_create
(struct htrdr* htrdr,
const double range[2], /* Must be included in [1000, 100000] nanometers */
- const size_t nbands, /* # bands used to discretisze the CIE tristimulus */
+ const size_t nbands, /* # bands used to discretized the CIE tristimulus */
const double ref_temperature,
struct htrdr_ran_lw** out_ran_lw)
{
struct htrdr_ran_lw* ran_lw = NULL;
res_T res = RES_OK;
- ASSERT(htrdr && range && nbands && out_ran_lw);
+ ASSERT(htrdr && range && out_ran_lw && ref_temperature > 0);
+ ASSERT(ref_temperature > 0);
+ ASSERT(range[0] >= 1000);
+ ASSERT(range[1] <= 100000);
+ ASSERT(range[0] < range[1]);
ran_lw = MEM_CALLOC(htrdr->allocator, 1, sizeof(*ran_lw));
if(!ran_lw) {
@@ -158,9 +303,22 @@ htrdr_ran_lw_create
ref_init(&ran_lw->ref);
ran_lw->htrdr = htrdr;
darray_double_init(htrdr->allocator, &ran_lw->cdf);
- darray_double_init(htrdr->allocator, &ran_lw->pdf);
+ darray_double_init(htrdr->allocator, &ran_lw->sky_bands_pdf);
- res = setup_ran_lw(ran_lw, FUNC_NAME, range, nbands, ref_temperature);
+ ran_lw->range[0] = range[0];
+ ran_lw->range[1] = range[1];
+ ran_lw->ref_temperature = ref_temperature;
+ ran_lw->nbands = nbands;
+
+ if(nbands == HTRDR_RAN_LW_CONTINUE) {
+ ran_lw->band_len = 0;
+ } else {
+ ran_lw->band_len = (range[1] - range[0]) / (double)nbands;
+ res = setup_ran_lw_cdf(ran_lw, FUNC_NAME);
+ if(res != RES_OK) goto error;
+ }
+
+ res = compute_sky_bands_pdf(ran_lw, FUNC_NAME);
if(res != RES_OK) goto error;
exit:
@@ -188,38 +346,36 @@ htrdr_ran_lw_ref_put(struct htrdr_ran_lw* ran_lw)
double
htrdr_ran_lw_sample
(const struct htrdr_ran_lw* ran_lw,
- const double r,
- double* pdf)
+ const double r)
{
- const double* cdf = NULL;
- const double* find = NULL;
- double r_next = nextafter(r, DBL_MAX);
- double wlen = 0;
- size_t cdf_length = 0;
- size_t i;
- ASSERT(ran_lw && 0 <= r && r < 1);
- ASSERT(darray_double_size_get(&ran_lw->cdf)
- == darray_double_size_get(&ran_lw->pdf));
-
- cdf = darray_double_cdata_get(&ran_lw->cdf);
- cdf_length = darray_double_size_get(&ran_lw->cdf);
-
- /* 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);
+ ASSERT(ran_lw);
+ if(ran_lw->band_len == 0) {
+ return ran_lw_sample_continue(ran_lw, r, FUNC_NAME);
+ } else {
+ return ran_lw_sample_discreet(ran_lw, r, FUNC_NAME);
+ }
+}
- i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
+double
+htrdr_ran_lw_get_sky_band_pdf
+ (const struct htrdr_ran_lw* ran_lw,
+ const size_t iband)
+{
+ size_t i;
+ size_t nbands;
+ (void)nbands;
+ ASSERT(ran_lw);
- /* Return the wavelength at the center of the sampled band */
- wlen = ran_lw->range[0] + ran_lw->band_len * ((double)i + 0.5);
- ASSERT(ran_lw->range[0] < wlen && wlen < ran_lw->range[1]);
+ nbands = htsky_get_spectral_bands_count(ran_lw->htrdr->sky);
+ ASSERT(nbands);
+ ASSERT(iband >= htsky_get_spectral_band_id(ran_lw->htrdr->sky, 0));
+ ASSERT(iband <= htsky_get_spectral_band_id(ran_lw->htrdr->sky, nbands-1));
- if(pdf) {
- *pdf = darray_double_cdata_get(&ran_lw->pdf)[i];
- }
+ /* Compute the index of the spectral band */
+ i = iband - htsky_get_spectral_band_id(ran_lw->htrdr->sky, 0);
- return wlen;
+ /* Fetch its PDF */
+ ASSERT(i < darray_double_size_get(&ran_lw->sky_bands_pdf));
+ return darray_double_cdata_get(&ran_lw->sky_bands_pdf)[i];
}
diff --git a/src/htrdr_ran_lw.h b/src/htrdr_ran_lw.h
@@ -18,6 +18,8 @@
#include <rsys/rsys.h>
+#define HTRDR_RAN_LW_CONTINUE 0
+
struct htrdr;
struct htrdr_ran_lw;
@@ -25,7 +27,9 @@ extern LOCAL_SYM res_T
htrdr_ran_lw_create
(struct htrdr* htrdr,
const double range[2], /* Must be included in [1000, 100000] nanometers */
- const size_t nbands, /* # bands used to discretisze the LW domain */
+ /* # bands used to discretisze the LW domain. HTRDR_RAN_LW_CONTINUE <=> no
+ * discretisation */
+ const size_t nbands,
const double ref_temperature, /* Reference temperature */
struct htrdr_ran_lw** ran_lw);
@@ -41,8 +45,12 @@ htrdr_ran_lw_ref_put
extern LOCAL_SYM double
htrdr_ran_lw_sample
(const struct htrdr_ran_lw* ran_lw,
- const double r,
- double* pdf); /* May be NULL */
+ const double r); /* Canonical number in [0, 1[ */
+
+extern LOCAL_SYM double
+htrdr_ran_lw_get_sky_band_pdf
+ (const struct htrdr_ran_lw* ran_lw,
+ const size_t iband);
#endif /* HTRDR_RAN_LW_H */
