commit 30e8f59edd9b7d188f9bb14f369659abead4494d
parent 1b56bcfd3e3951a89d3647ccd5defd9910a03c88
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 27 May 2020 16:05:28 +0200
Merge branch 'hotfix_shortwave' into develop
Diffstat:
6 files changed, 80 insertions(+), 55 deletions(-)
diff --git a/doc/htrdr.1.txt.in b/doc/htrdr.1.txt.in
@@ -53,7 +53,7 @@ computation is performed over the visible part of the spectrum in [380, 780]
nanometers, for the three components of the CIE 1931 XYZ colorimetric space
that are subsequently recombined in order to obtain the final color for each
pixel, and finally the whole image of the scene as seen from the set
-observation position. In long waves, the rendering is performed for the range
+observation position. In longwave, the rendering is performed for the range
of wavelengths defined by the *-l* option. The estimated radiance per pixel is
then converted to its brightness temperature and both are saved in the output
image.
@@ -141,13 +141,13 @@ OPTIONS
Number of samples per pixel estimation. In regular image rendering, a pixel
will use "3 * _samples-count_" Monte-Carlo realisations, one set of
_samples-count_ realisations for each X, Y and Z component of the CIE 1931
- XYZ color space. In long wave rendering (*-l* option) only one set of
+ XYZ color space. In longwave rendering (*-l* option) only one set of
_samples-count_ is used to estimate the pixel radiance and the resulting
brightness temperature for the submitted range of wavelengths. By default,
*spp* is set to @HTRDR_ARGS_DEFAULT_IMG_SPP@.
*-l* __wlen-min__,__wlen-max__::
- Switch in infrared rendering for the long wave interval in [__wlen-min__,
+ Switch in infrared rendering for the longwave interval in [__wlen-min__,
__wlen-max__] nanometers.
*-R*::
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_args.c b/src/htrdr_args.c
@@ -62,7 +62,7 @@ print_help(const char* cmd)
" page for the list of image options\n");
printf(
" -l WLEN_MIN,WLEN_MAX\n"
-" switch in infrared rendering for the long waves in\n"
+" switch in infrared rendering for longwave in\n"
" [WLEN_MIN, WLEN_MAX], in nanometers. By default, the\n"
" rendering is performed for the visible part of the\n"
" spectrum in [380, 780] nanometers.\n");
@@ -377,7 +377,7 @@ parse_lw_range(struct htrdr_args* args, const char* str)
if(res == RES_OK && len != 2) res = RES_BAD_ARG;
if(res == RES_OK && range[0] > range[1]) res = RES_BAD_ARG;
if(res != RES_OK) {
- fprintf(stderr, "Invalid long wave range `%s'.\n", str);
+ fprintf(stderr, "Invalid longwave range `%s'.\n", str);
goto error;
}
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>
@@ -58,7 +60,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;
@@ -169,7 +171,6 @@ 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 +180,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 +204,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);
@@ -274,10 +275,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,8 +297,21 @@ 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;
+
+ /* 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);
+ }
- res = setup_cie_xyz(cie, FUNC_NAME, range, nbands);
+ res = setup_cie_xyz(cie, FUNC_NAME, nbands);
if(res != RES_OK) goto error;
htrdr_log(htrdr, "Short wave interval defined on [%g, %g] nanometers.\n",
diff --git a/src/htrdr_ran_lw.c b/src/htrdr_ran_lw.c
@@ -59,14 +59,14 @@ setup_ran_lw_cdf
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",
+ "%s: Error allocating the CDF of the longwave spectral bands -- %s.\n",
func_name, res_to_cstr(res));
goto error;
}
res = darray_double_resize(&ran_lw->pdf, ran_lw->nbands);
if(res != RES_OK) {
htrdr_log_err(ran_lw->htrdr,
- "%s: Error allocating the pdf of the long wave spectral bands -- %s.\n",
+ "%s: Error allocating the pdf of the longwave spectral bands -- %s.\n",
func_name, res_to_cstr(res));
goto error;
}
@@ -75,10 +75,10 @@ setup_ran_lw_cdf
pdf = darray_double_data_get(&ran_lw->pdf); /* Now save pdf to correct weight */
htrdr_log(ran_lw->htrdr,
- "Number of bands of the long wave cumulative: %lu\n",
+ "Number of bands of the longwave cumulative: %lu\n",
(unsigned long)ran_lw->nbands);
- /* Compute the *unormalized* probability to sample a long wave band */
+ /* Compute the *unormalized* probability to sample a longwave band */
FOR_EACH(i, 0, ran_lw->nbands) {
double lambda_lo = ran_lw->range[0] + (double)i * ran_lw->band_len;
double lambda_hi = MMIN(lambda_lo + ran_lw->band_len, ran_lw->range[1]);
@@ -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,
@@ -316,7 +285,7 @@ htrdr_ran_lw_create
if(!ran_lw) {
res = RES_MEM_ERR;
htrdr_log_err(htrdr,
- "%s: could not allocate long wave random variate data structure -- %s.\n",
+ "%s: could not allocate longwave random variate data structure -- %s.\n",
FUNC_NAME, res_to_cstr(res));
goto error;
}
@@ -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;
