commit 3e60df19c701fd51aa43cdc939800aba6fd0030f
parent 689ceb1be730536399da4e5af57309a12534fcac
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 16 Apr 2020 22:20:01 +0200
Update the sampling of a long wave
First sample a wavelength. Then define its associated spectral band and
finally sample a quadrature point into it.
Diffstat:
6 files changed, 47 insertions(+), 140 deletions(-)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -24,6 +24,7 @@
#include "htrdr_camera.h"
#include "htrdr_ground.h"
#include "htrdr_mtl.h"
+#include "htrdr_ran_lw.h"
#include "htrdr_sun.h"
#include "htrdr_solve.h"
@@ -376,82 +377,6 @@ error:
goto exit;
}
-static res_T
-setup_lw_cdf(struct htrdr* htrdr)
-{
- /* Reference temperature used to compute the Long Wave cumulative */
- const double Tref = 290; /* In Kelvin */
- double* cdf = NULL;
- double* pdf = NULL;
- double sum = 0;
- size_t i;
- size_t nbands;
- res_T res = RES_OK;
- ASSERT(htrdr && htsky_is_long_wave(htrdr->sky));
-
- nbands = htsky_get_spectral_bands_count(htrdr->sky);
-
- res = darray_double_resize(&htrdr->lw_cdf, nbands);
- if(res != RES_OK) {
- htrdr_log_err(htrdr,
- "error allocating the CDF of the long wave spectral bands -- %s.\n",
- res_to_cstr(res));
- goto error;
- }
-
- res = darray_double_resize(&htrdr->lw_pdf, nbands);
- if(res != RES_OK) {
- htrdr_log_err(htrdr,
- "error allocating the PDF of the long wave spectral bands -- %s.\n",
- res_to_cstr(res));
- goto error;
- }
-
- pdf = darray_double_data_get(&htrdr->lw_pdf);
- cdf = darray_double_data_get(&htrdr->lw_cdf);
-
- /* 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(htrdr->sky, i);
- double wlens[2];
- HTSKY(get_spectral_band_bounds(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], Tref);
-
- /* Update the norm */
- sum += pdf[i];
- }
-
- /* Compute the cumulative of the previously computed probabilities */
- FOR_EACH(i, 0, nbands) {
- /* Normalize the probability */
- pdf[i] /= sum;
-
- /* Setup the cumulative */
- if(i == 0) {
- cdf[i] = pdf[i];
- } else {
- cdf[i] = pdf[i] + cdf[i-1];
- ASSERT(cdf[i] >= cdf[i-1]);
- }
- }
-
- cdf[nbands - 1] = 1.0;
-
-exit:
- return res;
-error:
- darray_double_clear(&htrdr->lw_cdf);
- darray_double_clear(&htrdr->lw_pdf);
- goto exit;
-}
-
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -480,12 +405,7 @@ htrdr_init
logger_set_stream(&htrdr->logger, LOG_OUTPUT, print_out, NULL);
logger_set_stream(&htrdr->logger, LOG_ERROR, print_err, NULL);
logger_set_stream(&htrdr->logger, LOG_WARNING, print_warn, NULL);
-
str_init(htrdr->allocator, &htrdr->output_name);
-
- darray_double_init(htrdr->allocator, &htrdr->lw_cdf);
- darray_double_init(htrdr->allocator, &htrdr->lw_pdf);
-
nthreads_max = MMAX(omp_get_max_threads(), omp_get_num_procs());
htrdr->dump_vtk = args->dump_vtk;
htrdr->verbose = args->verbose;
@@ -577,19 +497,38 @@ htrdr_init
iband1 = htsky_get_spectral_band_id(htrdr->sky, nbands-1);
HTSKY(get_spectral_band_bounds(htrdr->sky, iband0, wlen0));
HTSKY(get_spectral_band_bounds(htrdr->sky, iband1, wlen1));
- htrdr->wlen_range_m[0] = wlen0[0]*1e-9; /* Convert in meters */
- htrdr->wlen_range_m[1] = wlen1[1]*1e-9; /* Convert in meters */
- if(htsky_is_long_wave(htrdr->sky)) {
- /* 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);
+ /* Note that the bands are ranged in descending order wrt wavelength */
+ htrdr->wlen_range_m[1] = wlen0[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[0] = wlen1[1]*1e-9; /* Convert in meters */
+ ASSERT(htrdr->wlen_range_m[0] < htrdr->wlen_range_m[1]);
+
+ if(!htsky_is_long_wave(htrdr->sky)) { /* Short wave random variate */
+ const double* range = HTRDR_CIE_XYZ_RANGE_DEFAULT;
+ size_t n;
+
+ n = (size_t)(range[1] - range[0]);
+ res = htrdr_cie_xyz_create(htrdr, range, n, &htrdr->cie);
if(res != RES_OK) goto error;
+ } else { /* Long Wave random variate */
+ double range[2];
+ const double Tref = 290; /* In Kelvin */
+ size_t n;
+
+ range[0] = args->wlen_lw_range[0];
+ range[1] = args->wlen_lw_range[1];
+ n = (size_t)(range[1] - range[0]);
+
+ if(!n) {
+ if(eq_eps(range[0], range[1], 1.e-6)) {
+ range[0] -= 0.5;
+ range[1] += 0.5;
+ }
+ n = 1;
+ }
+ res = htrdr_ran_lw_create(htrdr, range, nbands, Tref, &htrdr->ran_lw);
+ if(res != RES_OK) goto error;
}
htrdr->lifo_allocators = MEM_CALLOC
@@ -656,6 +595,7 @@ htrdr_release(struct htrdr* htrdr)
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->ran_lw) htrdr_ran_lw_ref_put(htrdr->ran_lw);
if(htrdr->output && htrdr->output != stdout) fclose(htrdr->output);
if(htrdr->lifo_allocators) {
size_t i;
@@ -665,8 +605,6 @@ htrdr_release(struct htrdr* htrdr)
MEM_RM(htrdr->allocator, htrdr->lifo_allocators);
}
str_release(&htrdr->output_name);
- darray_double_release(&htrdr->lw_cdf);
- darray_double_release(&htrdr->lw_pdf);
logger_release(&htrdr->logger);
}
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -17,7 +17,6 @@
#ifndef HTRDR_H
#define HTRDR_H
-#include <rsys/dynamic_array_double.h>
#include <rsys/logger.h>
#include <rsys/ref_count.h>
#include <rsys/str.h>
@@ -51,6 +50,7 @@ struct htrdr {
struct htrdr_mtl* mtl;
struct htrdr_sun* sun;
struct htrdr_cie_xyz* cie;
+ struct htrdr_ran_lw* ran_lw;
struct htrdr_camera* cam;
struct htrdr_buffer* buf;
@@ -58,9 +58,6 @@ struct htrdr {
struct htsky* sky;
double wlen_range_m[2]; /* Integration range in *meters* */
- struct darray_double lw_cdf; /* CDF to sample a long wave bands */
- struct darray_double lw_pdf; /* PDF of the long wave bands*/
-
size_t spp; /* #samples per pixel */
size_t width; /* Image width */
size_t height; /* Image height */
diff --git a/src/htrdr_compute_radiance_lw.c b/src/htrdr_compute_radiance_lw.c
@@ -137,6 +137,7 @@ htrdr_compute_radiance_lw
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)
{
@@ -153,7 +154,6 @@ htrdr_compute_radiance_lw
double dir_next[3];
double band_bounds[2]; /* In nanometers */
double band_bounds_m[2]; /* In meters */
- double wlen;
double temperature;
double g;
double w = 0; /* Weight */
@@ -163,10 +163,6 @@ htrdr_compute_radiance_lw
/* Retrieve the band boundaries */
htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
- /* Arbitrarly use the wavelength at the center of the band as the wavelength
- * to use for data that depend on wavelength rather than spectral band as the
- * BRDF */
- wlen = (band_bounds[0] + band_bounds[1]) * 0.5;
band_bounds_m[0] = band_bounds[0] * 1e-9;
band_bounds_m[1] = band_bounds[1] * 1e-9;
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -21,6 +21,7 @@
#include "htrdr_buffer.h"
#include "htrdr_camera.h"
#include "htrdr_cie_xyz.h"
+#include "htrdr_ran_lw.h"
#include "htrdr_solve.h"
#include <high_tune/htsky.h>
@@ -479,36 +480,6 @@ error:
goto exit;
}
-static INLINE size_t
-sample_lw_spectral_interval
- (struct htrdr* htrdr,
- const double r,
- double* pdf)
-{
- const double* cdf = NULL;
- const double* find = NULL;
- double r_next = nextafter(r, DBL_MAX);
- size_t cdf_length = 0;
- size_t i;
- ASSERT(htrdr && r >= 0 && r < 1 && pdf);
- ASSERT(darray_double_size_get(&htrdr->lw_cdf)
- == darray_double_size_get(&htrdr->lw_pdf));
-
- cdf = darray_double_cdata_get(&htrdr->lw_cdf);
- cdf_length = darray_double_size_get(&htrdr->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);
-
- i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
- *pdf = darray_double_cdata_get(&htrdr->lw_pdf)[i];
-
- return htsky_get_spectral_band_id(htrdr->sky, i);
-}
-
static void
draw_pixel_sw
(struct htrdr* htrdr,
@@ -629,10 +600,11 @@ draw_pixel_lw
double ray_dir[3];
double weight;
double r0, r1;
+ double wlen;
size_t iband;
size_t iquad;
double usec;
- double band_pdf;
+ double pdf;
/* Begin the registration of the time spent to in the realisation */
time_current(&t0);
@@ -648,14 +620,17 @@ draw_pixel_lw
r0 = ssp_rng_canonical(rng);
r1 = ssp_rng_canonical(rng);
- /* Sample a spectral band and a quadrature point */
- iband = sample_lw_spectral_interval(htrdr, r0, &band_pdf);
+ /* Sample a wavelength */
+ wlen = htrdr_ran_lw_sample(htrdr->ran_lw, r0, &pdf);
+
+ /* 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);
/* Compute the luminance */
weight = htrdr_compute_radiance_lw
- (htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
- weight /= band_pdf;
+ (htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
+ weight /= 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
@@ -21,7 +21,7 @@
#include <rsys/algorithm.h>
#include <rsys/cstr.h>
-#include <rsys/dynamic_array.h>
+#include <rsys/dynamic_array_double.h>
#include <rsys/mem_allocator.h>
#include <rsys/ref_count.h>
@@ -148,7 +148,7 @@ htrdr_ran_lw_create
ASSERT(htrdr && range && nbands && out_ran_lw);
ran_lw = MEM_CALLOC(htrdr->allocator, 1, sizeof(*ran_lw));
- if(ran_lw) {
+ if(!ran_lw) {
res = RES_MEM_ERR;
htrdr_log_err(htrdr,
"%s: could not allocate long wave random variate data structure -- %s.\n",
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -63,7 +63,7 @@ htrdr_compute_radiance_sw
struct ssp_rng* rng,
const double pos[3],
const double dir[3],
- const double wlen,
+ const double wlen, /* In nanometer */
const size_t iband, /* Index of the spectral band */
const size_t iquad); /* Index of the quadrature point into the band */
@@ -74,6 +74,7 @@ htrdr_compute_radiance_lw
struct ssp_rng* rng,
const double pos[3],
const double dir[3],
+ const double wlen, /* In nanometer */
const size_t iband, /* Index of the spectral band */
const size_t iquad); /* Index of the quadrature point into the band */
