commit ec918ad1f19de3c3aa36a992e3f35e413e5a91f2
parent 2279307e6c2b839fe1a7f5efc85a718a409ef9e3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 9 Jun 2020 15:14:02 +0200
Merge branch 'release_0.5'
Diffstat:
25 files changed, 1244 insertions(+), 1046 deletions(-)
diff --git a/README.md b/README.md
@@ -59,12 +59,39 @@ informations on CMake.
## Release notes
+### Version 0.5
+
+#### New feature
+
+Add support of shortwave integration with respect to the Planck function for a
+reference temperature whose default value is the blackbody temperature of the
+sun. Actually this is the counterpart of the longwave integration introduced by
+the "infrared rendering" in the 0.4 version. The main difference is that the
+source of radiation is the sun rather than the medium and its boundaries.
+
+The option `-l` that enabled the infrared rendering is now replaced by the new
+`-s` option that controls the spectral integration that can be CIE XYZ (i.e.
+regular image rendering), longwave or shortwave.
+
+#### Fixes
+
+- Fix the returned sun radiance: the precomputed per spectral band solar
+ incoming flux is removed and the sun radiance is now retrieved by directly
+ evaluating the monochromatic Planck for the blackbody temperature of the sun.
+- Fix CIE XYZ spectral integration: the pdf used to sample the CIE tristimulus
+ values was not correctly handled in the Monte-Carlo weight.
+- Fix the longwave spectral integration: the Monte-Carlo weight was wrong
+ leading to overestimated temperatures.
+
### Version 0.4
+#### New features
+
- Add support of infrared rendering: when defined, the new `-l` option setups
- the range of long waves into which the rendering is performed. In infrared
+ the range of longwave into which the rendering is performed. In infrared
rendering, each pixel stores the radiance per pixel and its associated
- brightness temperature.
+ brightness temperature. Spectral integration is done with respect to the
+ Planck function for a reference temperature of 290 K.
- The ground geometry can now have several materials whose data vary over the
spectrum. These materials are listed in a new
[htrdr-materials](https://gitlab.com/meso-star/htrdr/-/blob/master/doc/htrdr-materials.5.txt)
@@ -79,6 +106,9 @@ informations on CMake.
center of the sampled spectral band. Consequently, high
resolution data defined per wavelength (e.g. Mie's properties and the
reflectivity of the materials) are now fully taken into account.
+
+#### Fixes
+
- Fix a deadlock when `htrdr` was run through MPI.
- Fix a memory leak: the output file was not closed on exit.
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -25,11 +25,11 @@ option(NO_TEST "Do not build the tests" OFF)
# Check dependencies
################################################################################
find_package(AW 2.0 REQUIRED)
-find_package(HTSky 0.1 REQUIRED)
+find_package(HTSky 0.2 REQUIRED)
find_package(MruMtl 0.0 REQUIRED)
find_package(RCMake 0.3 REQUIRED)
find_package(RSys 0.7 REQUIRED)
-find_package(Star3D 0.6 REQUIRED)
+find_package(Star3D 0.7.1 REQUIRED)
find_package(StarSF 0.6 REQUIRED)
find_package(StarSP 0.8 REQUIRED)
find_package(StarVX 0.1 REQUIRED)
@@ -59,7 +59,7 @@ include_directories(
# Generate files
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 4)
+set(VERSION_MINOR 5)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -100,9 +100,10 @@ set(HTRDR_FILES_SRC
htrdr_interface.c
htrdr_main.c
htrdr_mtl.c
- htrdr_ran_lw.c
+ htrdr_ran_wlen.c
htrdr_rectangle.c
htrdr_slab.c
+ htrdr_spectral.c
htrdr_sun.c)
set(HTRDR_FILES_INC
htrdr.h
@@ -114,9 +115,10 @@ set(HTRDR_FILES_INC
htrdr_ground.h
htrdr_interface.h
htrdr_mtl.h
- htrdr_ran_lw.h
+ htrdr_ran_wlen.h
htrdr_rectangle.h
htrdr_slab.h
+ htrdr_spectral.h
htrdr_sun.h
htrdr_solve.h)
set(HTRDR_FILES_INC2 # Generated files
diff --git a/doc/htrdr-image.5.txt b/doc/htrdr-image.5.txt
@@ -30,25 +30,33 @@ thus ignored as well as empty lines. The first valid line stores 2 unsigned
integers that represent the image definition, i.e. the number of pixels per
line and per column. Then each line stores 8 pixel components.
-If the image is a regular rendering in the visible part of the spectrum, the
-pixel components are actually 4 pairs of floating points data representing the
-pixel color encoded in the CIE 1931 XYZ color space and the per radiative path
-computation time. The first, second and third pair encode the estimated
-integrated radiance in W.sr^-1.m^-2 of the X, Y and Z pixel components,
-respectively. The first value of each pair is the expected value of the
-estimated radiance while the second one is its associated standard deviation.
-The fourth pair saves the estimate in microseconds of the per radiative path
-computation time and its standard error.
-
-If the image is an infrared rendering, the first and second pixel components
-store the expected value and the standard error, respectively, of the
-estimated brightness temperature in Kelvin. The third and fourth component
-save the expected value and standard deviation of the pixel radiance in
-W.sr^-1.m^-2.nm^-1. The fifth and sixth pixel components are unused. Finally
-the last 2 pixel components save, as for a regular rendering, the estimate in
+If the image is a regular rendering in the visible part of the spectrum
+(*-s* _cie_xyz_ option in *htrdr*(1)), the pixel components are actually 4
+pairs of floating points data representing the pixel color encoded in the CIE
+1931 XYZ color space and the per radiative path computation time. The first,
+second and third pairs encode the estimated integrated radiance in W/sr/m^2 of
+the X, Y and Z pixel components, respectively. The first value of each pair is
+the expected value of the estimated radiance while the second one is its
+associated standard deviation. The fourth pair saves the estimate in
microseconds of the per radiative path computation time and its standard
error.
+If the image is an infrared rendering (*-s* *lw*=_wlen-min_,_wlen_max_ option
+in *htrdr*(1)), the first and second pixel components store the expected value
+and the standard error, respectively, of the estimated brightness temperature
+in Kelvin. The third and fourth components save the expected value and standard
+deviation of the pixel radiance that is either an integrated radiance in
+W/sr/m^2 or a spectral radiance in W/sr/m^2/nm whether this radiance is
+computed for a spectral range or for one wavelength. The fifth and sixth pixel
+components are unused. Finally the last 2 pixel components save, as for a
+regular rendering, the estimate in microseconds of the per radiative path
+computation time and its standard error.
+
+If it was generating from a shortwave rendering (*-s*
+*sw*=_wlen-min_,wlen-max_ option in *htrdr*(1)) the image is formatted as in
+longwave rendering mode exepted that the first and second pixel components are
+unused since no brightness temperature was evaluated in shortwave.
+
Pixels are sorted line by line, with the origin defined at the top left corner
of the image. With an image definition of N by M pixels, with N the number of
pixels per line and M the overall number of lines in the image, the first N
diff --git a/doc/htrdr-obj.5.txt b/doc/htrdr-obj.5.txt
@@ -60,7 +60,7 @@ f 1 2 3
f 3 2 4
-------
-Define a wooden cube whose faces Z-aligned faces are against a brick material.
+Define a wooden cube whose Z-aligned faces are against a brick material.
The remaining faces are in contact with the air.
[verse]
diff --git a/doc/htrdr.1.txt.in b/doc/htrdr.1.txt.in
@@ -30,7 +30,13 @@ SYNOPSIS
DESCRIPTION
-----------
*htrdr* is an image renderer of scenes composed of an atmospheric gas mixture,
-clouds, and a ground. Images can be rendered in the visible or the infrared
+clouds, and a ground. It evaluates the intensity incoming on each pixel of the
+sensor array. The underlying algorithm is based on a Monte-Carlo method: it
+consists in simulating a given number of optical paths originating from the
+camera, directed into the atmosphere, taking into account light absorption and
+scattering phenomena.
+
+Images can be rendered in the visible or the infrared
part of the spectrum. It uses spectral data that should be provided for the
pressure and temperature atmospheric vertical profile [1] (*-a* _atmosphere_),
the liquid water content in suspension within the clouds stored in a *htcp*(5)
@@ -44,27 +50,28 @@ whose materials are listed in the *htrdr-material*(5) file provided through the
*-M* option. Both, the clouds and the ground, can be infinitely repeated along
the X and Y axis by setting the *-r* and the *-R* options, respectively.
-*htrdr* evaluates the intensity incoming on each pixel of the sensor array. The
-underlying algorithm is based on a Monte-Carlo method: it consists in
-simulating a given number of optical paths originating from the camera,
-directed into the atmosphere, taking into account light absorption and
-scattering phenomena. When rendering in the visible part of the spectrum, the
-computation is performed over the visible part of the spectrum in [380, 780]
+Spectral dimension can be integrated in many ways (*-s* option). By default,
+the computation is performed for 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
-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.
-
-In *htrdr* the spatial unit 1.0 corresponds to one meter, the estimated
-radiance is given in W.sr^-1.m^-2 and the temperatures are expressed in Kelvin.
-The results are written to the output file if the *-o* option is defined and
-the standard output otherwise. The output image is a list of raw ASCII data
-formatted with respect to the *htrdr-image*(5) file format. Since *htrdr*
-relies on the Monte-Carlo method, each estimation is given with its numerical
-accuracy.
+observation position. The two other ways consist in explicitly defining the
+longwave or shortwave spectral range to handle and continuously sampling a
+wavelength in this range according to the Planck function for a reference
+temperature. Actually longwave and shortwave are keywords that mean that the
+source of radiation is whether external or internal to the medium,
+respectively. In shortwave, only the pixel radiance is evaluated and stored in
+the output image. In longwave this estimated radiance is then converted to its
+brightness temperature and both are saved in the image.
+
+In *htrdr* the spatial unit 1.0 corresponds to one meter and the temperatures
+are expressed in Kelvin. The estimated radiances are given in W/sr/m^2 excepted
+for monochromatic computations where the computed spectral radiance is defined
+in W/sr/m^2/nm. The results are written to the output file if the *-o* option
+is defined and the standard output otherwise. The output image is a list of raw
+ASCII data formatted with respect to the *htrdr-image*(5) file format. Since
+*htrdr* relies on the Monte-Carlo method, each estimation is given with its
+numerical accuracy.
During the simulation, *htrdr* dynamically loads/unloads cloud properties to
handle clouds whose data that do not feat in main memory. *htrdr* also supports
@@ -141,15 +148,11 @@ 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__,
- __wlen-max__] nanometers.
-
*-R*::
Infinitely repeat the _ground_ along the X and Y axis.
@@ -180,6 +183,44 @@ OPTIONS
File where *htrdr* writes its _output_ data. If not defined, write results to
standard output.
+*-s* <__spectral-parameter__:...>::
+ define the type and the range of the spectral integration. Available spectral
+ parameters are:
+
+ **cie_xyz**;;
+ the radiance is computed for the visible part of the spectrum in [380, 780]
+ nanometers with respect to the XYZ CIE 1931 tristimulus values. This is the
+ default comportment of *htrdr*.
+
+ **lw**=*_wlen-min_*,*_wlen-max_*;;
+ perform the spectral sampling continuously in the [_wlen-min_, _wlen-max_]
+ wavelength range (wavelength must be provided in nanometers) according to
+ the Planck function for a reference temperature. If _wlen-min_ and
+ _wlen-max_ are equals the computation is monochromatic. *lw* means for
+ longwave but is here a code word that really means "computation of radiance
+ using the internal source of radiation": in other words, radiation is
+ emitted by the medium and its boundaries. Because the application is for
+ the terrestrial atmosphere, internal radiation is emitted in the thermal
+ longwave part of the electromagnetic spectrum. Therefore the default value
+ of the reference temperature used in the spectral sampling is fixed at 290
+ K.
+
+ **sw**=*_wlen-min_*,*_wlen-max_*;;
+ perform the spectral sampling continuously in the [_wlen-min_, _wlen-max_]
+ wavelength range (wavelength must be provided in nanometers) according to
+ the Planck function for a reference temperature. If _wlen-min_ and
+ _wlen-max_ are equals the computation is monochromatic. In the present
+ case, *sw* means that the source of radiation is external to the medium:
+ because the application is for the terrestrial atmosphere, the value of the
+ reference temperature is by default fixed at 5778 K, i.e. the blackbody
+ temperature of the Sun.
+
+ **Tref**=**_temperature_**;;
+ reference temperature of the Planck function used to continuously sample the
+ longwave/shortwave spectral range. In longwave, it is set to 290 K by
+ default while in shortwave its default value is the blackbody temperature
+ of the sun (i.e. 5778 K).
+
*-T* _threshold_::
Optical thickness used as threshold criteria to partition the properties of
the clouds. By default its value is @HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@.
@@ -237,13 +278,14 @@ PPM image [4]:
$ htpp -o image.ppm output
Render the previous scene in infrared for the wavelengths in [9200, 10000]
-nanometers:
+nanometers with a reference temperature of 300 Kelvin:
- $ htrdr -l 9200,10000 -a gas.txt -m Mie.nc -g mountains.obj -R \
+ $ htrdr -a gas.txt -m Mie.nc -g mountains.obj -R \
-M materials.mtl \
-c clouds.htcp \
-C pos=0,0,400:tgt=0,1,0:up=0,0,1 \
-i def=800x600:spp=64 \
+ -s lw=9200,1000:Tref=300
-f -o output
Move the sun by setting its azimuthal and elevation angles to *120* and *40*
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -24,7 +24,7 @@
#include "htrdr_camera.h"
#include "htrdr_ground.h"
#include "htrdr_mtl.h"
-#include "htrdr_ran_lw.h"
+#include "htrdr_ran_wlen.h"
#include "htrdr_sun.h"
#include "htrdr_solve.h"
@@ -102,6 +102,23 @@ log_msg
}
}
+static enum htsky_spectral_type
+htrdr_to_sky_spectral_type(const enum htrdr_spectral_type type)
+{
+ enum htsky_spectral_type spectype;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ spectype = HTSKY_SPECTRAL_LW;
+ break;
+ case HTRDR_SPECTRAL_SW:
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ spectype = HTSKY_SPECTRAL_SW;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return spectype;
+}
+
static res_T
dump_buffer
(struct htrdr* htrdr,
@@ -117,13 +134,8 @@ dump_buffer
ASSERT(htrdr && buf && stream_name && stream);
(void)stream_name;
- if(htsky_is_long_wave(htrdr->sky)) {
- pixsz = sizeof(struct htrdr_pixel_lw);
- pixal = ALIGNOF(struct htrdr_pixel_lw);
- } else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
- }
+ pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
htrdr_buffer_get_layout(buf, &layout);
if(layout.elmt_size != pixsz || layout.alignment != pixal) {
@@ -140,8 +152,8 @@ dump_buffer
struct htrdr_estimate pix_time = HTRDR_ESTIMATE_NULL;
const struct htrdr_accum* pix_time_acc = NULL;
- if(htsky_is_long_wave(htrdr->sky)) {
- const struct htrdr_pixel_lw* pix = htrdr_buffer_at(buf, x, y);
+ if(htrdr->spectral_type != HTRDR_SPECTRAL_SW_CIE_XYZ){
+ const struct htrdr_pixel_xwave* pix = htrdr_buffer_at(buf, x, y);
fprintf(stream, "%g %g ",
pix->radiance_temperature.E, pix->radiance_temperature.SE);
fprintf(stream, "%g %g ", pix->radiance.E, pix->radiance.SE);
@@ -149,7 +161,7 @@ dump_buffer
pix_time_acc = &pix->time;
} else {
- const struct htrdr_pixel_sw* pix = htrdr_buffer_at(buf, x, y);
+ const struct htrdr_pixel_image* pix = htrdr_buffer_at(buf, x, y);
fprintf(stream, "%g %g ", pix->X.E, pix->X.SE);
fprintf(stream, "%g %g ", pix->Y.E, pix->Y.SE);
fprintf(stream, "%g %g ", pix->Z.E, pix->Z.SE);
@@ -389,6 +401,7 @@ htrdr_init
struct htsky_args htsky_args = HTSKY_ARGS_DEFAULT;
double proj_ratio;
double sun_dir[3];
+ double spectral_range[2];
const char* output_name = NULL;
size_t ithread;
int nthreads_max;
@@ -414,6 +427,8 @@ htrdr_init
htrdr->grid_max_definition[0] = args->grid_max_definition[0];
htrdr->grid_max_definition[1] = args->grid_max_definition[1];
htrdr->grid_max_definition[2] = args->grid_max_definition[2];
+ htrdr->spectral_type = args->spectral_type;
+ htrdr->ref_temperature = args->ref_temperature;
res = init_mpi(htrdr);
if(res != RES_OK) goto error;
@@ -484,32 +499,48 @@ htrdr_init
htsky_args.nthreads = htrdr->nthreads;
htsky_args.repeat_clouds = args->repeat_clouds;
htsky_args.verbose = htrdr->mpi_rank == 0 ? args->verbose : 0;
- htsky_args.wlen_lw_range[0] = args->wlen_lw_range[0];
- htsky_args.wlen_lw_range[1] = args->wlen_lw_range[1];
+ htsky_args.spectral_type = htrdr_to_sky_spectral_type(args->spectral_type);
+ htsky_args.wlen_range[0] = args->wlen_range[0];
+ htsky_args.wlen_range[1] = args->wlen_range[1];
res = htsky_create(&htrdr->logger, htrdr->allocator, &htsky_args, &htrdr->sky);
if(res != RES_OK) goto error;
- if(!htsky_is_long_wave(htrdr->sky)) { /* Short wave random variate */
- const double* range = HTRDR_CIE_XYZ_RANGE_DEFAULT;
- size_t n;
+ HTSKY(get_raw_spectral_bounds(htrdr->sky, spectral_range));
- n = (size_t)(range[1] - range[0]);
- res = htrdr_cie_xyz_create(htrdr, range, n, &htrdr->cie);
+ spectral_range[0] = MMAX(args->wlen_range[0], spectral_range[0]);
+ spectral_range[1] = MMIN(args->wlen_range[1], spectral_range[1]);
+ if(spectral_range[0] != args->wlen_range[0]
+ || spectral_range[1] != args->wlen_range[1]) {
+ htrdr_log_warn(htrdr,
+ "%s: the submitted spectral range overflowed the spectral data.\n", FUNC_NAME);
+ }
+
+ htrdr->wlen_range_m[0] = spectral_range[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[1] = spectral_range[1]*1e-9; /* Convert in meters */
+
+ if(htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) {
+ size_t n;
+ n = (size_t)(spectral_range[1] - spectral_range[0]);
+ res = htrdr_cie_xyz_create(htrdr, spectral_range, n, &htrdr->cie);
if(res != RES_OK) goto error;
- } else { /* Long Wave random variate */
- const double Tref = 290; /* In Kelvin */
+ } else {
size_t n;
- htrdr->wlen_range_m[0] = args->wlen_lw_range[0]*1e-9; /* Convert in meters */
- htrdr->wlen_range_m[1] = args->wlen_lw_range[1]*1e-9; /* Convert in meters */
+ if(htrdr->ref_temperature <= 0) {
+ htrdr_log_err(htrdr, "%s: invalid reference temperature %g K.\n",
+ FUNC_NAME, htrdr->ref_temperature);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
+ n = (size_t)(spectral_range[1] - spectral_range[0]);
- n = (size_t)(args->wlen_lw_range[1] - args->wlen_lw_range[0]);
- res = htrdr_ran_lw_create
- (htrdr, args->wlen_lw_range, n, Tref, &htrdr->ran_lw);
+ res = htrdr_ran_wlen_create
+ (htrdr, spectral_range, n, htrdr->ref_temperature, &htrdr->ran_wlen);
if(res != RES_OK) goto error;
- }
+ }
htrdr->lifo_allocators = MEM_CALLOC
(htrdr->allocator, htrdr->nthreads, sizeof(*htrdr->lifo_allocators));
@@ -535,16 +566,9 @@ htrdr_init
/* Create the image buffer only on the master process; the image parts
* rendered by the processes are gathered onto the master process. */
if(!htrdr->dump_vtk && htrdr->mpi_rank == 0) {
- size_t pixsz = 0; /* sizeof(pixel) */
- size_t pixal = 0; /* alignof(pixel) */
+ const size_t pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ const size_t pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
- if(htsky_is_long_wave(htrdr->sky)) {
- pixsz = sizeof(struct htrdr_pixel_lw);
- pixal = ALIGNOF(struct htrdr_pixel_lw);
- } else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
- }
res = htrdr_buffer_create(htrdr,
args->image.definition[0], /* Width */
args->image.definition[1], /* Height */
@@ -575,7 +599,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->ran_wlen) htrdr_ran_wlen_ref_put(htrdr->ran_wlen);
if(htrdr->output && htrdr->output != stdout) fclose(htrdr->output);
if(htrdr->lifo_allocators) {
size_t i;
@@ -709,64 +733,37 @@ htrdr_fflush(struct htrdr* htrdr, FILE* stream)
/*******************************************************************************
* Local functions
******************************************************************************/
-res_T
-brightness_temperature
- (struct htrdr* htrdr,
- const double lambda_min,
- const double lambda_max,
- const double radiance,
- double* temperature)
+double
+compute_sky_min_band_len
+ (struct htsky* sky,
+ const double range[2])
{
- const size_t MAX_ITER = 100;
- const double epsilon_T = 1e-4; /* In K */
- const double epsilon_B = radiance * 1e-8;
- double T, T0, T1, T2;
- double B, B0;
- size_t i;
- res_T res = RES_OK;
- ASSERT(temperature && lambda_min <= lambda_max);
+ double min_band_len = DBL_MAX;
+ size_t nbands;
+ ASSERT(sky && range && range[0] <= range[1]);
- /* Search for a brightness temperature whose radiance is greater than or
- * equal to the estimated radiance */
- T2 = 200;
- FOR_EACH(i, 0, MAX_ITER) {
- const double B2 = planck(lambda_min, lambda_max, T2);
- if(B2 >= radiance) break;
- T2 *= 2;
- }
- if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+ nbands = htsky_get_spectral_bands_count(sky);
- B0 = T0 = T1 = 0;
- FOR_EACH(i, 0, MAX_ITER) {
- T = (T1+T2)*0.5;
- B = planck(lambda_min, lambda_max, T);
+ if(eq_eps(range[0], range[1], 1.e-6)) {
+ ASSERT(nbands == 1);
+ min_band_len = 0;
+ } else {
+ size_t i = 0;
- if(B < radiance) {
- T1 = T;
- } else {
- T2 = T;
- }
+ /* 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));
- if(fabs(T-T0) < epsilon_T || fabs(B-B0) < epsilon_B)
- break;
+ /* Adjust band boundaries to the submitted range */
+ wlens[0] = MMAX(wlens[0], range[0]);
+ wlens[1] = MMIN(wlens[1], range[1]);
- T0 = T;
- B0 = B;
+ min_band_len = MMIN(wlens[1] - wlens[0], min_band_len);
+ }
}
- if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
-
- *temperature = T;
-
-exit:
- return res;
-error:
- htrdr_log_err(htrdr,
- "Could not compute the brightness temperature for the radiance %g "
- "estimated over [%g, %g] nanometers.\n",
- radiance,
- lambda_min*1e9,
- lambda_max*1e9);
- goto exit;
+ return min_band_len;
}
res_T
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -17,6 +17,8 @@
#ifndef HTRDR_H
#define HTRDR_H
+#include "htrdr_spectral.h"
+
#include <rsys/logger.h>
#include <rsys/ref_count.h>
#include <rsys/str.h>
@@ -50,13 +52,15 @@ struct htrdr {
struct htrdr_mtl* mtl;
struct htrdr_sun* sun;
struct htrdr_cie_xyz* cie;
- struct htrdr_ran_lw* ran_lw;
+ struct htrdr_ran_wlen* ran_wlen;
struct htrdr_camera* cam;
struct htrdr_buffer* buf;
struct htsky* sky;
+ enum htrdr_spectral_type spectral_type;
double wlen_range_m[2]; /* Integration range in *meters* */
+ double ref_temperature; /* Reference temperature in Kelvin */
size_t spp; /* #samples per pixel */
size_t width; /* Image width */
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -61,16 +61,6 @@ print_help(const char* cmd)
" -i <image> define the image to compute. Refer to the htrdr man\n"
" 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"
-" [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");
- printf(
-" -R infinitely repeat the ground along the X and Y axis.\n");
- printf(
-" -r infinitely repeat the clouds along the X and Y axis.\n");
- printf(
" -M MATERIALS file listing the ground materials.\n");
printf(
" -m MIE file of Mie's data.\n");
@@ -81,6 +71,14 @@ print_help(const char* cmd)
" -o OUTPUT file where data are written. If not defined, data are\n"
" written to standard output.\n");
printf(
+" -R infinitely repeat the ground along the X and Y axis.\n");
+ printf(
+" -r infinitely repeat the clouds along the X and Y axis.\n");
+ printf(
+" -s <spectral> define the type and range of the spectral\n"
+" integration. Refer to the htrdr man page for the list\n"
+" of spectral options\n");
+ printf(
" -T THRESHOLD optical thickness used as threshold during the octree\n"
" building. By default its value is `%g'.\n",
HTRDR_ARGS_DEFAULT.optical_thickness);
@@ -140,7 +138,7 @@ parse_fov(const char* str, double* out_fov)
fprintf(stderr, "Invalid field of view `%s'.\n", str);
return RES_BAD_ARG;
}
- if(fov < 30 || fov > 120) {
+ if(fov <= 0 || fov >= 180) {
fprintf(stderr, "The field of view %g is not in [30, 120].\n", fov);
return RES_BAD_ARG;
}
@@ -217,10 +215,11 @@ parse_camera_parameter(struct htrdr_args* args, const char* str)
char* val;
char* ctx;
res_T res = RES_OK;
+ ASSERT(args);
if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
fprintf(stderr,
- "Could not duplicate the rectangle option string `%s'.\n", str);
+ "Could not duplicate the camera option string `%s'.\n", str);
res = RES_MEM_ERR;
goto error;
}
@@ -262,6 +261,89 @@ error:
}
static res_T
+parse_spectral_range(const char* str, double wlen_range[2])
+{
+ double range[2];
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(wlen_range && str);
+
+ res = cstr_to_list_double(str, ',', range, &len, 2);
+ 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 spectral range `%s'.\n", str);
+ goto error;
+ }
+ wlen_range[0] = range[0];
+ wlen_range[1] = range[1];
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_spectral_parameter(struct htrdr_args* args, const char* str)
+{
+ char buf[128];
+ char* key;
+ char* val;
+ char* ctx;
+ res_T res = RES_OK;
+ ASSERT(args);
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr,
+ "Could not duplicate the spectral option string `%s'.\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &ctx);
+ val = strtok_r(NULL, "", &ctx);
+
+ if(!strcmp(key, "cie_xyz")) {
+ args->spectral_type = HTRDR_SPECTRAL_SW_CIE_XYZ;
+ args->wlen_range[0] = HTRDR_CIE_XYZ_RANGE_DEFAULT[0];
+ args->wlen_range[1] = HTRDR_CIE_XYZ_RANGE_DEFAULT[1];
+ } else {
+ if(!val) {
+ fprintf(stderr, "Missing value to the spectral option `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(!strcmp(key, "sw")) {
+ args->spectral_type = HTRDR_SPECTRAL_SW;
+ res = parse_spectral_range(val, args->wlen_range);
+ if(res != RES_OK) goto error;
+ } else if(!strcmp(key, "lw")) {
+ args->spectral_type = HTRDR_SPECTRAL_LW;
+ res = parse_spectral_range(val, args->wlen_range);
+ if(res != RES_OK) goto error;
+ } else if(!strcmp(key, "Tref")) {
+ res = cstr_to_double(val, &args->ref_temperature);
+ if(res == RES_OK && args->ref_temperature < 0) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ fprintf(stderr, "Invalid reference temperature Tref=%s.\n", val);
+ goto error;
+ }
+ } else {
+ fprintf(stderr, "Invalid spectral parameter `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
parse_multiple_parameters
(struct htrdr_args* args,
const char* str,
@@ -365,31 +447,6 @@ error:
goto exit;
}
-static res_T
-parse_lw_range(struct htrdr_args* args, const char* str)
-{
- double range[2];
- size_t len;
- res_T res = RES_OK;
- ASSERT(args && str);
-
- res = cstr_to_list_double(str, ',', range, &len, 2);
- 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);
- goto error;
- }
-
- args->wlen_lw_range[0] = range[0];
- args->wlen_lw_range[1] = range[1];
-
-exit:
- return res;
-error:
- goto exit;
-}
-
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -414,7 +471,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
}
}
- while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:l:M:m:O:o:RrT:t:V:v")) != -1) {
+ while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:M:m:O:o:Rrs:T:t:V:v")) != -1) {
switch(opt) {
case 'a': args->filename_gas = optarg; break;
case 'C':
@@ -435,15 +492,16 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
res = parse_multiple_parameters
(args, optarg, parse_image_parameter);
break;
- case 'l':
- res = parse_lw_range(args, optarg);
- break;
case 'M': args->filename_mtl = optarg; break;
case 'm': args->filename_mie = optarg; break;
case 'O': args->cache = optarg; break;
case 'o': args->output = optarg; break;
case 'r': args->repeat_clouds = 1; break;
case 'R': args->repeat_ground = 1; break;
+ case 's':
+ res = parse_multiple_parameters
+ (args, optarg, parse_spectral_parameter);
+ break;
case 'T':
res = cstr_to_double(optarg, &args->optical_thickness);
if(res == RES_OK && args->optical_thickness < 0) res = RES_BAD_ARG;
@@ -483,6 +541,22 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
goto error;
}
+ /* Setup default ref temperature if necessary */
+ if(args->ref_temperature <= 0) {
+ switch(args->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ args->ref_temperature = HTRDR_DEFAULT_LW_REF_TEMPERATURE;
+ break;
+ case HTRDR_SPECTRAL_SW:
+ args->ref_temperature = HTRDR_SUN_TEMPERATURE;
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ args->ref_temperature = -1; /* Unused */
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+
exit:
return res;
error:
diff --git a/src/htrdr_args.h.in b/src/htrdr_args.h.in
@@ -16,7 +16,8 @@
#ifndef HTRDR_ARGS_H
#define HTRDR_ARGS_H
-#include "htrdr_ground.h"
+#include "htrdr_spectral.h"
+#include "htrdr_cie_xyz.h"
#include <float.h>
#include <limits.h>
@@ -55,7 +56,9 @@ struct htrdr_args {
double optical_thickness; /* Threshold used during octree building */
unsigned grid_max_definition[3]; /* Maximum definition of the grid */
- double wlen_lw_range[2]; /* Long Wave range to handle in nm */
+ enum htrdr_spectral_type spectral_type;
+ double wlen_range[2]; /* Spectral range of integration in nm */
+ double ref_temperature; /* Planck reference temperature in Kelvin */
unsigned nthreads; /* Hint on the number of threads to use */
int force_overwriting;
@@ -92,7 +95,9 @@ struct htrdr_args {
90, /* Sun elevation */ \
@HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@, /* Optical thickness */ \
{UINT_MAX, UINT_MAX, UINT_MAX}, /* Maximum definition of the grid */ \
- {DBL_MAX, -DBL_MAX}, /* Long wave range. Degenerated <=> short wave */ \
+ HTRDR_SPECTRAL_SW_CIE_XYZ, /* Spectral type */ \
+ HTRDR_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
+ -1, /* Reference temperature */ \
(unsigned)~0, /* #threads */ \
0, /* Force overwriting */ \
0, /* dump VTK */ \
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -18,7 +18,6 @@
#define HTRDR_C_H
#include <rsys/rsys.h>
-#include <rsys/math.h> /* PI support */
#ifndef NDEBUG
#define MPI(Func) ASSERT(MPI_##Func == MPI_SUCCESS)
@@ -36,9 +35,6 @@ enum htrdr_mpi_message {
struct htrdr;
-#define SW_WAVELENGTH_MIN 380 /* In nanometer */
-#define SW_WAVELENGTH_MAX 780 /* In nanometer */
-
/* In nanometer */
static FINLINE double
wavenumber_to_wavelength(const double nu/*In cm^-1*/)
@@ -96,110 +92,12 @@ morton_xyz_decode_u21(const uint64_t code, uint32_t xyz[3])
xyz[2] = (uint32_t)morton3D_decode_u21(code >> 0);
}
-static INLINE double
-wiebelt(const double v)
-{
- int m;
- double w, v2, v4;
- /*.153989717364e+00;*/
- const double fifteen_over_pi_power_4 = 15.0/(PI*PI*PI*PI);
- const double z0 = 1.0/3.0;
- const double z1 = 1.0/8.0;
- const double z2 = 1.0/60.0;
- const double z4 = 1.0/5040.0;
- const double z6 = 1.0/272160.0;
- const double z8 = 1.0/13305600.0;
-
- if(v >= 2.) {
- w = 0;
- for(m=1; m<6 ;m++)
- w+=exp(-m*v)/(m*m*m*m) * (((m*v+3)*m*v+6)*m*v+6);
- w = w * fifteen_over_pi_power_4;
- } else {
- v2 = v*v;
- v4 = v2*v2;
- w = z0 - z1*v + z2*v2 - z4*v2*v2 + z6*v4*v2 - z8*v4*v4;
- w = 1. - fifteen_over_pi_power_4*v2*v*w;
- }
- ASSERT(w >= 0.0 && w <= 1.0);
- return w;
-}
-
-static INLINE double
-blackbody_fraction
- (const double lambda0, /* In meters */
- const double lambda1, /* In meters */
- const double temperature) /* In Kelvin */
-{
- const double C2 = 1.43877735e-2; /* m.K */
- double x0 = C2 / lambda0;
- double x1 = C2 / lambda1;
- double v0 = x0 / temperature;
- double v1 = x1 / temperature;
- double w0 = wiebelt(v0);
- double w1 = wiebelt(v1);
- return w1 - w0;
-}
-
-
-
-/* Return the Planck value in W/m^2/sr/m at a given wavelength */
-static INLINE double
-planck_monochromatic
- (const double lambda, /* In meters */
- const double temperature) /* In Kelvin */
-{
- const double c = 299792458; /* m/s */
- const double h = 6.62607015e-34; /* J.s */
- const double k = 1.380649e-23; /* J/K */
- const double lambda2 = lambda*lambda;
- const double lambda5 = lambda2*lambda2*lambda;
- const double B = ((2.0 * h * c*c) / lambda5) /* W/m^2/sr/m */
- / (exp(h*c/(lambda*k*temperature))-1.0);
- ASSERT(temperature > 0);
- return B;
-}
-
-/* Return the average Planck value in W/m^2/sr/m over the
- * [lambda_min, lambda_max] interval. */
-static INLINE double
-planck_interval
- (const double lambda_min, /* In meters */
- const double lambda_max, /* In meters */
- const double temperature) /* In Kelvin */
-{
- const double T2 = temperature*temperature;
- const double T4 = T2*T2;
- const double BOLTZMANN_CONSTANT = 5.6696e-8; /* W/m^2/K^4 */
- ASSERT(lambda_min < lambda_max && temperature > 0);
- return blackbody_fraction(lambda_min, lambda_max, temperature)
- * BOLTZMANN_CONSTANT * T4 / (PI * (lambda_max-lambda_min)); /* In W/m^2/sr/m */
-}
-
-/* Invoke planck_monochromatic or planck_interval whether the submitted
- * interval is null or not, respectively. The returned value is in W/m^2/sr/m */
-static FINLINE double
-planck
- (const double lambda_min, /* In meters */
- const double lambda_max, /* In meters */
- const double temperature) /* In Kelvin */
-{
- ASSERT(lambda_min <= lambda_max && temperature > 0);
- if(lambda_min == lambda_max) {
- return planck_monochromatic(lambda_min, temperature);
- } else {
- return planck_interval(lambda_min, lambda_max, temperature);
- }
-}
-
-extern LOCAL_SYM res_T
-brightness_temperature
- (struct htrdr* htrdr,
- const double lambda_min, /* In meters */
- const double lambda_max, /* In meters */
- /* Integrated over [lambda_min, lambda_max]. In W/m^2/sr/m */
- const double radiance,
- double* temperature);
+/* 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
open_output_stream
diff --git a/src/htrdr_camera.c b/src/htrdr_camera.c
@@ -73,7 +73,7 @@ htrdr_camera_create
ref_init(&cam->ref);
cam->htrdr = htrdr;
- if(fov <= 0) {
+ if(fov <= 0 || fov >= PI) {
htrdr_log_err(htrdr, "invalid horizontal camera field of view `%g'\n", fov);
res = RES_BAD_ARG;
goto error;
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>
@@ -31,6 +33,9 @@ struct htrdr_cie_xyz {
struct darray_double cdf_X;
struct darray_double cdf_Y;
struct darray_double cdf_Z;
+ double rcp_integral_X;
+ double rcp_integral_Y;
+ double rcp_integral_Z;
double range[2]; /* Boundaries of the handled CIE XYZ color space */
double band_len; /* Length in nanometers of a band */
@@ -58,7 +63,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;
@@ -159,7 +164,12 @@ sample_cie_xyz
} else if(lambda_min <= lambda_2 && lambda_2 < lambda_max) {
lambda = lambda_2;
} else {
- FATAL("Unexpected error.\n");
+ htrdr_log_warn(cie->htrdr,
+ "%s: cannot sample a wavelength in [%g, %g[. The possible wavelengths"
+ "were %g and %g.\n",
+ FUNC_NAME, lambda_min, lambda_max, lambda_1, lambda_2);
+ /* Arbitrarly choose the wavelength at the center of the sampled band */
+ lambda = (lambda_min + lambda_max)*0.5;
}
return lambda;
@@ -169,20 +179,19 @@ 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 */
double* pdf[3] = {NULL, NULL, NULL};
double* cdf[3] = {NULL, NULL, NULL};
- double sum[3] = {0,0,0};
+ double sum[3] = {0, 0, 0};
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 +212,12 @@ 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);
@@ -217,6 +225,15 @@ setup_cie_xyz
sum[Y] += pdf[Y][i];
sum[Z] += pdf[Z][i];
}
+ #define CHK_SUM(Sum, Range, Fit) \
+ ASSERT(eq_eps(Sum, trapezoidal_integration(Range[0], Range[1], Fit), 1.e-3))
+ CHK_SUM(sum[X], cie->range, fit_x_bar_1931);
+ CHK_SUM(sum[Y], cie->range, fit_y_bar_1931);
+ CHK_SUM(sum[Z], cie->range, fit_z_bar_1931);
+ #undef CHK_SUM
+ cie->rcp_integral_X = 1.0 / sum[X];
+ cie->rcp_integral_Y = 1.0 / sum[Y];
+ cie->rcp_integral_Z = 1.0 / sum[Z];
FOR_EACH(i, 0, nbands) {
/* Normalize the pdf */
@@ -274,10 +291,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,11 +313,24 @@ 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",
+ htrdr_log(htrdr, "CIE XYZ spectral interval defined on [%g, %g] nanometers.\n",
range[0], range[1]);
exit:
@@ -327,29 +359,38 @@ double
htrdr_cie_xyz_sample_X
(struct htrdr_cie_xyz* cie,
const double r0,
- const double r1)
+ const double r1,
+ double* pdf)
{
- return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_X;
+ return wlen;
}
double
htrdr_cie_xyz_sample_Y
(struct htrdr_cie_xyz* cie,
const double r0,
- const double r1)
+ const double r1,
+ double* pdf)
{
- return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_Y;
+ return wlen;
}
double
htrdr_cie_xyz_sample_Z
(struct htrdr_cie_xyz* cie,
const double r0,
- const double r1)
+ const double r1,
+ double* pdf)
{
- return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_Z;
+ return wlen;
}
diff --git a/src/htrdr_cie_xyz.h b/src/htrdr_cie_xyz.h
@@ -22,7 +22,9 @@ 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};
+#define HTRDR_CIE_XYZ_RANGE_DEFAULT__ {380, 780}
+static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] =
+ HTRDR_CIE_XYZ_RANGE_DEFAULT__;
extern LOCAL_SYM res_T
htrdr_cie_xyz_create
@@ -43,19 +45,22 @@ htrdr_cie_xyz_ref_put
extern LOCAL_SYM double
htrdr_cie_xyz_sample_X
(struct htrdr_cie_xyz* cie,
- const double r0, const double r1); /* Canonical numbers in [0, 1[ */
+ const double r0, const double r1, /* Canonical numbers in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
/* Return a wavelength in nanometer */
extern LOCAL_SYM double
htrdr_cie_xyz_sample_Y
(struct htrdr_cie_xyz* cie,
- const double r0, const double r1); /* Canonical number in [0, 1[ */
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
/* Return a wavelength in nanometer */
extern LOCAL_SYM double
htrdr_cie_xyz_sample_Z
(struct htrdr_cie_xyz* cie,
- const double r0, const double r1); /* Canonical number in [0, 1[ */
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
#endif /* HTRDR_cie_xyz_H */
diff --git a/src/htrdr_compute_radiance_lw.c b/src/htrdr_compute_radiance_lw.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+/*Spectralt (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
@@ -152,22 +152,13 @@ htrdr_compute_radiance_lw
double range[2];
double pos_next[3];
double dir_next[3];
- double band_bounds[2]; /* In nanometers */
- double band_bounds_m[2]; /* In meters */
double temperature;
+ double wlen_m = wlen * 1.e-9;
double g;
double w = 0; /* Weight */
ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads);
- /* Retrieve the band boundaries */
- htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
-
- /* Transform the band boundaries in meters and clamp them to the integration
- * domain */
- band_bounds_m[0] = MMAX(band_bounds[0] * 1e-9, htrdr->wlen_range_m[0]);
- band_bounds_m[1] = MMIN(band_bounds[1] * 1e-9, htrdr->wlen_range_m[1]);
-
/* Setup the phase function for this spectral band & quadrature point */
CHK(RES_OK == ssf_phase_create
(&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
@@ -220,7 +211,8 @@ htrdr_compute_radiance_lw
if(ctx.event_type == EVENT_ABSORPTION) {
ASSERT(!SVX_HIT_NONE(&svx_hit));
temperature = htsky_fetch_temperature(htrdr->sky, pos_next);
- w = planck(band_bounds_m[0], band_bounds_m[1], temperature);
+ /* weight is planck integrated over the spectral sub-interval */
+ w = planck_monochromatic(wlen_m, temperature);
break;
}
@@ -268,7 +260,8 @@ htrdr_compute_radiance_lw
if(temperature <= 0) {
w = 0;
} else {
- w = planck(band_bounds_m[0], band_bounds_m[1], temperature);
+ /* weight is planck integrated over the spectral sub-interval */
+ w = planck_monochromatic(wlen_m, temperature);
}
break;
}
@@ -279,7 +272,6 @@ htrdr_compute_radiance_lw
d3_set(dir, dir_next);
}
SSF(phase_ref_put(phase_hg));
-
return w;
}
diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c
@@ -48,7 +48,7 @@ static const struct scattering_context SCATTERING_CONTEXT_NULL = {
struct transmissivity_context {
struct ssp_rng* rng;
const struct htsky* sky;
- size_t iband; /* Index of the spectrald */
+ size_t iband; /* Index of the spectral */
size_t iquad; /* Index of the quadrature point into the band */
double Ts; /* Sampled optical thickness */
@@ -109,6 +109,13 @@ scattering_hit_filter
ASSERT(ctx->traversal_dst >= hit->distance[0]);
ASSERT(ctx->traversal_dst <= hit->distance[1]);
+ /* Stop the ray whenever the traversal distance without any scattering
+ * event is too high. It means the maximum scattering coefficient has a
+ * very small value, and the returned radiance is null. This can only
+ * happen when the voxel has a [quasi] infinite length in the propagation
+ * direction. */
+ if(ctx->traversal_dst > 1e9) break;
+
/* 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];
@@ -279,7 +286,8 @@ htrdr_compute_radiance_sw
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));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_SW
+ || htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ);
CHK(RES_OK == ssf_phase_create
(&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
@@ -298,7 +306,6 @@ htrdr_compute_radiance_sw
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);
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -21,7 +21,7 @@
#include "htrdr_buffer.h"
#include "htrdr_camera.h"
#include "htrdr_cie_xyz.h"
-#include "htrdr_ran_lw.h"
+#include "htrdr_ran_wlen.h"
#include "htrdr_solve.h"
#include <high_tune/htsky.h>
@@ -46,13 +46,13 @@
STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
enum pixel_format {
- PIXEL_LW,
- PIXEL_SW
+ PIXEL_XWAVE,
+ PIXEL_IMAGE
};
union pixel {
- struct htrdr_pixel_lw lw;
- struct htrdr_pixel_sw sw;
+ struct htrdr_pixel_xwave xwave;
+ struct htrdr_pixel_image image;
};
/* Tile of row ordered image pixels */
@@ -101,6 +101,19 @@ morton2D_encode(const uint16_t u16)
return u32;
}
+static INLINE enum pixel_format
+spectral_type_to_pixfmt(const enum htrdr_spectral_type spectral_type)
+{
+ enum pixel_format pixfmt;
+ switch(spectral_type) {
+ case HTRDR_SPECTRAL_LW: pixfmt = PIXEL_XWAVE; break;
+ case HTRDR_SPECTRAL_SW: pixfmt = PIXEL_XWAVE; break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ: pixfmt = PIXEL_IMAGE; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return pixfmt;
+}
+
static FINLINE struct tile*
tile_create(struct mem_allocator* allocator, const enum pixel_format fmt)
{
@@ -155,9 +168,11 @@ tile_at
ASSERT(tile && x < TILE_SIZE && y < TILE_SIZE);
return tile->data.pixels + (y*TILE_SIZE + x);
}
+
static void
write_tile_data
- (struct htrdr_buffer* buf,
+ (struct htrdr* htrdr,
+ struct htrdr_buffer* buf,
const struct tile_data* tile_data)
{
struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL;
@@ -165,13 +180,12 @@ write_tile_data
size_t irow_tile;
size_t ncols_tile, nrows_tile;
char* buf_mem;
- ASSERT(buf && tile_data);
+ ASSERT(htrdr && buf && tile_data);
+ (void)htrdr;
htrdr_buffer_get_layout(buf, &layout);
buf_mem = htrdr_buffer_get_data(buf);
- ASSERT(layout.elmt_size == (tile_data->format == PIXEL_LW
- ? sizeof(struct htrdr_pixel_lw)
- : sizeof(struct htrdr_pixel_sw)));
+ ASSERT(layout.elmt_size==htrdr_spectral_type_get_pixsz(htrdr->spectral_type));
/* Compute the row/column of the tile origin into the buffer */
icol = tile_data->x * (size_t)TILE_SIZE;
@@ -190,11 +204,11 @@ write_tile_data
FOR_EACH(x, 0, ncols_tile) {
switch(tile_data->format) {
- case PIXEL_LW:
- ((struct htrdr_pixel_lw*)buf_col)[x] = tile_row[x].lw;
+ case PIXEL_XWAVE:
+ ((struct htrdr_pixel_xwave*)buf_col)[x] = tile_row[x].xwave;
break;
- case PIXEL_SW:
- ((struct htrdr_pixel_sw*)buf_col)[x] = tile_row[x].sw;
+ case PIXEL_IMAGE:
+ ((struct htrdr_pixel_image*)buf_col)[x] = tile_row[x].image;
break;
default: FATAL("Unreachable code.\n"); break;
}
@@ -445,17 +459,18 @@ mpi_gather_tiles
LIST_FOR_EACH(node, tiles) {
struct tile* t = CONTAINER_OF(node, struct tile, node);
- write_tile_data(buf, &t->data);
+ write_tile_data(htrdr, buf, &t->data);
++itile;
}
if(itile != ntiles) {
+ enum pixel_format pixfmt;
ASSERT(htrdr->mpi_nprocs > 1);
/* Create a temporary tile to receive the tile data computed by the
* concurrent MPI processes */
- tile = tile_create
- (htrdr->allocator, htsky_is_long_wave(htrdr->sky) ? PIXEL_LW : PIXEL_SW);
+ pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
+ tile = tile_create(htrdr->allocator, pixfmt);
if(!tile) {
res = RES_MEM_ERR;
htrdr_log_err(htrdr,
@@ -468,7 +483,7 @@ mpi_gather_tiles
FOR_EACH(itile, itile, ntiles) {
MPI(Recv(&tile->data, (int)msg_sz, MPI_CHAR, MPI_ANY_SOURCE,
HTRDR_MPI_TILE_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE));
- write_tile_data(buf, &tile->data);
+ write_tile_data(htrdr, buf, &tile->data);
}
}
}
@@ -481,7 +496,7 @@ error:
}
static void
-draw_pixel_sw
+draw_pixel_image
(struct htrdr* htrdr,
const size_t ithread,
const size_t ipix[2],
@@ -489,13 +504,13 @@ draw_pixel_sw
const struct htrdr_camera* cam,
const size_t spp,
struct ssp_rng* rng,
- struct htrdr_pixel_sw* pixel)
+ struct htrdr_pixel_image* pixel)
{
struct htrdr_accum XYZ[3]; /* X, Y, and Z */
struct htrdr_accum time;
size_t ichannel;
ASSERT(ipix && ipix && pix_sz && cam && rng && pixel);
- ASSERT(!htsky_is_long_wave(htrdr->sky));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ);
/* Reset accumulators */
XYZ[0] = HTRDR_ACCUM_NULL;
@@ -514,6 +529,7 @@ draw_pixel_sw
double weight;
double r0, r1, r2;
double wlen; /* Sampled wavelength into the spectral band */
+ double pdf;
size_t iband; /* Sampled spectral band */
size_t iquad; /* Sampled quadrature point into the spectral band */
double usec;
@@ -535,20 +551,23 @@ draw_pixel_sw
/* Sample a spectral band and a quadrature point */
switch(ichannel) {
- 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;
+ case 0: wlen = htrdr_cie_xyz_sample_X(htrdr->cie, r0, r1, &pdf); break;
+ case 1: wlen = htrdr_cie_xyz_sample_Y(htrdr->cie, r0, r1, &pdf); break;
+ case 2: wlen = htrdr_cie_xyz_sample_Z(htrdr->cie, r0, r1, &pdf); break;
default: FATAL("Unreachable code.\n"); break;
}
iband = htsky_find_spectral_band(htrdr->sky, wlen);
iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r2, iband);
- /* Compute the luminance */
+ /* Compute the radiance in W/m^2/sr/m */
weight = htrdr_compute_radiance_sw
(htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
ASSERT(weight >= 0);
+ pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
+ weight /= pdf; /* In W/m^2/sr */
+
/* End the registration of the per realisation time */
time_sub(&t0, time_current(&t1), &t0);
usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
@@ -572,8 +591,39 @@ draw_pixel_sw
pixel->time = time;
}
+
+static INLINE double
+radiance_temperature
+ (struct htrdr* htrdr,
+ const double radiance) /* In W/m^2/sr */
+{
+ double temperature = 0;
+ double radiance_avg = radiance;
+ res_T res = RES_OK;
+ ASSERT(htrdr && radiance >= 0);
+
+ /* From integrated radiance to average radiance in W/m^2/sr/m */
+ if(htrdr->wlen_range_m[0] != htrdr->wlen_range_m[1]) { /* !monochromatic */
+ radiance_avg /= (htrdr->wlen_range_m[1] - htrdr->wlen_range_m[0]);
+ }
+
+ res = brightness_temperature
+ (htrdr,
+ htrdr->wlen_range_m[0],
+ htrdr->wlen_range_m[1],
+ radiance_avg,
+ &temperature);
+ if(res != RES_OK) {
+ htrdr_log_warn(htrdr,
+ "Could not compute the brightness temperature for the radiance %g.\n",
+ radiance_avg);
+ temperature = 0;
+ }
+ return temperature;
+}
+
static void
-draw_pixel_lw
+draw_pixel_xwave
(struct htrdr* htrdr,
const size_t ithread,
const size_t ipix[2],
@@ -581,14 +631,15 @@ draw_pixel_lw
const struct htrdr_camera* cam,
const size_t spp,
struct ssp_rng* rng,
- struct htrdr_pixel_lw* pixel)
+ struct htrdr_pixel_xwave* pixel)
{
struct htrdr_accum radiance;
struct htrdr_accum time;
size_t isamp;
double temp_min, temp_max;
ASSERT(ipix && ipix && pix_sz && cam && rng && pixel);
- ASSERT(htsky_is_long_wave(htrdr->sky));
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_LW
+ || htrdr->spectral_type == HTRDR_SPECTRAL_SW);
/* Reset the pixel accumulators */
radiance = HTRDR_ACCUM_NULL;
@@ -600,7 +651,7 @@ draw_pixel_lw
double ray_org[3];
double ray_dir[3];
double weight;
- double r0, r1;
+ double r0, r1, r2;
double wlen;
size_t iband;
size_t iquad;
@@ -620,23 +671,32 @@ draw_pixel_lw
r0 = ssp_rng_canonical(rng);
r1 = ssp_rng_canonical(rng);
+ r2 = ssp_rng_canonical(rng);
/* Sample a wavelength */
- wlen = htrdr_ran_lw_sample(htrdr->ran_lw, r0);
+ wlen = htrdr_ran_wlen_sample(htrdr->ran_wlen, r0, r1, &band_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);
-
- /* Retrieve the PDF to sample this sky band */
- band_pdf = htrdr_ran_lw_get_sky_band_pdf(htrdr->ran_lw, iband);
-
- /* Compute the luminance in W/m^2/sr/m */
- weight = htrdr_compute_radiance_lw
- (htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
- weight /= band_pdf;
+ iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r2, iband);
+
+ /* Compute the spectral radiance in W/m^2/sr/m */
+ switch(htrdr->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ weight = htrdr_compute_radiance_lw(htrdr, ithread, rng, ray_org,
+ ray_dir, wlen, iband, iquad);
+ break;
+ case HTRDR_SPECTRAL_SW:
+ weight = htrdr_compute_radiance_sw(htrdr, ithread, rng, ray_org,
+ ray_dir, wlen, iband, iquad);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
ASSERT(weight >= 0);
+ /* Importance sampling: correct weight with pdf */
+ weight /= band_pdf; /* In W/m^2/sr */
+
/* End the registration of the per realisation time */
time_sub(&t0, time_current(&t1), &t0);
usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
@@ -659,30 +719,13 @@ draw_pixel_lw
pixel->time = time;
/* Compute the brightness_temperature of the pixel and estimate its standard
- * error */
- #define BRIGHTNESS_TEMPERATURE(Radiance, Temperature) { \
- res_T res = brightness_temperature \
- (htrdr, \
- htrdr->wlen_range_m[0], \
- htrdr->wlen_range_m[1], \
- (Radiance), \
- &(Temperature)); \
- if(res != RES_OK) { \
- htrdr_log_warn(htrdr, \
- "Could not compute the brightness temperature for the radiance %g.\n", \
- (Radiance)); \
- (Temperature) = 0; \
- } \
- } (void)0
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E, pixel->radiance_temperature.E);
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E - pixel->radiance.SE, temp_min);
- BRIGHTNESS_TEMPERATURE(pixel->radiance.E + pixel->radiance.SE, temp_max);
- pixel->radiance_temperature.SE = temp_max - temp_min;
- #undef BRIGHTNESS_TEMPERATURE
-
- /* Transform the pixel radiance from W/m^2/sr/m to W/m^/sr/nm */
- pixel->radiance.E *= 1.0e-9;
- pixel->radiance.SE *= 1.0e-9;
+ * error if the sources were in the medium (<=> longwave) */
+ if(htrdr->spectral_type == HTRDR_SPECTRAL_LW) {
+ pixel->radiance_temperature.E = radiance_temperature(htrdr, pixel->radiance.E);
+ temp_min = radiance_temperature(htrdr, pixel->radiance.E - pixel->radiance.SE);
+ temp_max = radiance_temperature(htrdr, pixel->radiance.E + pixel->radiance.SE);
+ pixel->radiance_temperature.SE = temp_max - temp_min;
+ }
}
static res_T
@@ -724,10 +767,15 @@ draw_tile
ipix[1] = tile_org[1] + ipix_tile[1];
/* Draw the pixel */
- if(htsky_is_long_wave(htrdr->sky)) {
- draw_pixel_lw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->lw);
- } else {
- draw_pixel_sw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->sw);
+ switch(htrdr->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ draw_pixel_xwave(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ draw_pixel_image(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
}
}
return RES_OK;
@@ -773,7 +821,7 @@ draw_image
htrdr->mpi_working_procs[htrdr->mpi_rank] = 0;
--htrdr->mpi_nworking_procs;
- pixfmt = htsky_is_long_wave(htrdr->sky) ? PIXEL_LW : PIXEL_SW;
+ pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
omp_set_num_threads((int)nthreads);
#pragma omp parallel
@@ -927,19 +975,13 @@ htrdr_draw_radiance
proc_work_init(htrdr->allocator, &work);
if(htrdr->mpi_rank == 0) {
- size_t pixsz, pixal;
+ const size_t pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
+ const size_t pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
+
htrdr_buffer_get_layout(buf, &layout);
ASSERT(layout.width || layout.height || layout.elmt_size);
ASSERT(layout.width == width && layout.height == height);
- if(htsky_is_long_wave(htrdr->sky)) {
- pixsz = sizeof(struct htrdr_pixel_lw);
- pixal = ALIGNOF(struct htrdr_pixel_lw);
- } else {
- pixsz = sizeof(struct htrdr_pixel_sw);
- pixal = ALIGNOF(struct htrdr_pixel_sw);
- }
-
if(layout.elmt_size != pixsz || layout.alignment < pixal) {
htrdr_log_err(htrdr, "%s: invalid buffer layout.\n", FUNC_NAME);
res = RES_BAD_ARG;
diff --git a/src/htrdr_ran_lw.c b/src/htrdr_ran_lw.c
@@ -1,417 +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 /* nextafter */
-
-#include "htrdr.h"
-#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>
-#include <rsys/mem_allocator.h>
-#include <rsys/ref_count.h>
-
-#include <math.h> /* nextafter */
-
-struct htrdr_ran_lw {
- struct darray_double cdf;
- /* 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;
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static res_T
-setup_ran_lw_cdf
- (struct htrdr_ran_lw* ran_lw,
- 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 && ran_lw->nbands != HTRDR_RAN_LW_CONTINUE);
-
- 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",
- func_name, res_to_cstr(res));
- goto error;
- }
-
- cdf = darray_double_data_get(&ran_lw->cdf);
- pdf = cdf; /* Alias the CDF by the PDF since only one array is necessary */
-
- htrdr_log(ran_lw->htrdr,
- "Number of bands of the long wave cumulative: %lu\n",
- (unsigned long)ran_lw->nbands);
-
- /* Compute the *unormalized* probability to sample a long wave 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]);
- ASSERT(lambda_lo<= lambda_hi);
- ASSERT(lambda_lo > ran_lw->range[0] || eq_eps(lambda_lo, ran_lw->range[0], 1.e-6));
- ASSERT(lambda_lo < ran_lw->range[1] || eq_eps(lambda_lo, ran_lw->range[1], 1.e-6));
-
- /* Convert from nanometer to meter */
- lambda_lo *= 1.e-9;
- lambda_hi *= 1.e-9;
-
- /* Compute the probability of the current band */
- 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, ran_lw->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]);
- }
- }
- 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);
- goto exit;
-}
-
-static res_T
-compute_sky_bands_pdf
- (struct htrdr_ran_lw* ran_lw,
- const char* func_name,
- double* out_min_band_len)
-{
- double* pdf = NULL;
- double min_band_len = DBL_MAX;
- 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);
-
- if(eq_eps(ran_lw->range[0], ran_lw->range[1], 1.e-6)) {
- ASSERT(nbands == 1);
- pdf[0] = 1;
- min_band_len = 0;
- } else {
- double sum = 0;
- 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);
-
- /* Convert from nanometer to meter */
- wlens[0] *= 1.e-9;
- wlens[1] *= 1.e-9;
-
- /* Compute the probability of the current band */
- pdf[i] = blackbody_fraction(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:
- *out_min_band_len = min_band_len;
- 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_M = 1e-15;
- 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.0e-9;
- range_m[1] = ran_lw->range[1] * 1.0e-9;
-
- /* 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_M
- || 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)
-{
- struct htrdr_ran_lw* ran_lw = NULL;
- ASSERT(ref);
- ran_lw = CONTAINER_OF(ref, struct htrdr_ran_lw, ref);
- darray_double_release(&ran_lw->cdf);
- darray_double_release(&ran_lw->sky_bands_pdf);
- MEM_RM(ran_lw->htrdr->allocator, ran_lw);
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-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 discretized CDF */
- const double ref_temperature,
- struct htrdr_ran_lw** out_ran_lw)
-{
- struct htrdr_ran_lw* ran_lw = NULL;
- double min_band_len = 0;
- res_T res = RES_OK;
- 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) {
- res = RES_MEM_ERR;
- htrdr_log_err(htrdr,
- "%s: could not allocate long wave random variate data structure -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
- }
- ref_init(&ran_lw->ref);
- ran_lw->htrdr = htrdr;
- darray_double_init(htrdr->allocator, &ran_lw->cdf);
- darray_double_init(htrdr->allocator, &ran_lw->sky_bands_pdf);
-
- ran_lw->range[0] = range[0];
- ran_lw->range[1] = range[1];
- ran_lw->ref_temperature = ref_temperature;
- ran_lw->nbands = nbands;
-
- res = compute_sky_bands_pdf(ran_lw, FUNC_NAME, &min_band_len);
- if(res != RES_OK) goto error;
-
- if(nbands == HTRDR_RAN_LW_CONTINUE) {
- ran_lw->band_len = 0;
- } else {
- ran_lw->band_len = (range[1] - range[0]) / (double)ran_lw->nbands;
-
- /* Adjust the band length to ensure that each sky spectral interval is
- * overlapped by at least one band */
- if(ran_lw->band_len > min_band_len) {
- ran_lw->band_len = min_band_len;
- ran_lw->nbands = (size_t)ceil((range[1] - range[0]) / ran_lw->band_len);
- }
-
- res = setup_ran_lw_cdf(ran_lw, FUNC_NAME);
- if(res != RES_OK) goto error;
- }
-
- htrdr_log(htrdr, "Long wave interval defined on [%g, %g] nanometers.\n",
- range[0], range[1]);
-
-exit:
- *out_ran_lw = ran_lw;
- return res;
-error:
- if(ran_lw) htrdr_ran_lw_ref_put(ran_lw);
- goto exit;
-}
-
-void
-htrdr_ran_lw_ref_get(struct htrdr_ran_lw* ran_lw)
-{
- ASSERT(ran_lw);
- ref_get(&ran_lw->ref);
-}
-
-void
-htrdr_ran_lw_ref_put(struct htrdr_ran_lw* ran_lw)
-{
- ASSERT(ran_lw);
- ref_put(&ran_lw->ref, release_ran_lw);
-}
-
-double
-htrdr_ran_lw_sample
- (const struct htrdr_ran_lw* ran_lw,
- const double r)
-{
- ASSERT(ran_lw);
- if(ran_lw->nbands != HTRDR_RAN_LW_CONTINUE) {
- return ran_lw_sample_discreet(ran_lw, r, FUNC_NAME);
- } else if(eq_eps(ran_lw->range[0], ran_lw->range[1], 1.e-6)) {
- return ran_lw->range[0];
- } else {
- return ran_lw_sample_continue(ran_lw, r, FUNC_NAME);
- }
-}
-
-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);
-
- 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));
-
- /* Compute the index of the spectral band */
- i = iband - htsky_get_spectral_band_id(ran_lw->htrdr->sky, 0);
-
- /* 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
@@ -1,56 +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_RAN_LW_H
-#define HTRDR_RAN_LW_H
-
-#include <rsys/rsys.h>
-
-#define HTRDR_RAN_LW_CONTINUE 0
-
-struct htrdr;
-struct htrdr_ran_lw;
-
-extern LOCAL_SYM res_T
-htrdr_ran_lw_create
- (struct htrdr* htrdr,
- const double range[2], /* Must be included in [1000, 100000] nanometers */
- /* # bands used to discretisze the LW domain. HTRDR_RAN_LW_CONTINUE <=> no
- * discretisation */
- const size_t nbands, /* Hint on #bands used to discretised th CDF */
- const double ref_temperature, /* Reference temperature */
- struct htrdr_ran_lw** ran_lw);
-
-extern LOCAL_SYM void
-htrdr_ran_lw_ref_get
- (struct htrdr_ran_lw* ran_lw);
-
-extern LOCAL_SYM void
-htrdr_ran_lw_ref_put
- (struct htrdr_ran_lw* ran_lw);
-
-/* Return a wavelength in nanometer */
-extern LOCAL_SYM double
-htrdr_ran_lw_sample
- (const struct htrdr_ran_lw* ran_lw,
- 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 */
-
diff --git a/src/htrdr_ran_wlen.c b/src/htrdr_ran_wlen.c
@@ -0,0 +1,364 @@
+/* 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/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "htrdr.h"
+#include "htrdr_c.h"
+#include "htrdr_ran_wlen.h"
+
+#include <high_tune/htsky.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_ran_wlen {
+ struct darray_double pdf;
+ struct darray_double cdf;
+ double range[2]; /* Boundaries of the spectral integration interval */
+ double band_len; /* Length in nanometers of a band */
+ double ref_temperature; /* In Kelvin */
+ size_t nbands; /* # bands */
+
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+setup_wlen_ran_cdf
+ (struct htrdr_ran_wlen* wlen_ran,
+ const char* func_name)
+{
+ double* pdf = NULL;
+ double* cdf = NULL;
+ double sum = 0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(wlen_ran && func_name && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
+
+ res = darray_double_resize(&wlen_ran->cdf, wlen_ran->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(wlen_ran->htrdr,
+ "%s: Error allocating the CDF of the spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+ res = darray_double_resize(&wlen_ran->pdf, wlen_ran->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(wlen_ran->htrdr,
+ "%s: Error allocating the pdf of the spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+
+ cdf = darray_double_data_get(&wlen_ran->cdf);
+ pdf = darray_double_data_get(&wlen_ran->pdf); /* Now save pdf to correct weight */
+
+ htrdr_log(wlen_ran->htrdr,
+ "Number of bands of the spectrum cumulative: %lu\n",
+ (unsigned long)wlen_ran->nbands);
+
+ /* Compute the *unnormalized* probability to sample a small band */
+ FOR_EACH(i, 0, wlen_ran->nbands) {
+ double lambda_lo = wlen_ran->range[0] + (double)i * wlen_ran->band_len;
+ double lambda_hi = MMIN(lambda_lo + wlen_ran->band_len, wlen_ran->range[1]);
+ ASSERT(lambda_lo<= lambda_hi);
+ ASSERT(lambda_lo > wlen_ran->range[0] || eq_eps(lambda_lo, wlen_ran->range[0], 1.e-6));
+ ASSERT(lambda_lo < wlen_ran->range[1] || eq_eps(lambda_lo, wlen_ran->range[1], 1.e-6));
+
+ /* Convert from nanometer to meter */
+ lambda_lo *= 1.e-9;
+ lambda_hi *= 1.e-9;
+
+ /* Compute the probability of the current band */
+ pdf[i] = blackbody_fraction(lambda_lo, lambda_hi, wlen_ran->ref_temperature);
+
+ /* Update the norm */
+ sum += pdf[i];
+ }
+
+ /* Compute the cumulative of the previously computed probabilities */
+ FOR_EACH(i, 0, wlen_ran->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]);
+ }
+ }
+ ASSERT(eq_eps(cdf[wlen_ran->nbands-1], 1, 1.e-6));
+ cdf[wlen_ran->nbands - 1] = 1.0; /* Handle numerical issue */
+
+exit:
+ return res;
+error:
+ darray_double_clear(&wlen_ran->cdf);
+ darray_double_clear(&wlen_ran->pdf);
+ goto exit;
+}
+
+static double
+wlen_ran_sample_continue
+ (const struct htrdr_ran_wlen* wlen_ran,
+ const double r,
+ const double range[2], /* In nanometer */
+ const char* func_name,
+ double* pdf)
+{
+ /* Numerical parameters of the dichotomy search */
+ const size_t MAX_ITER = 100;
+ const double EPSILON_LAMBDA_M = 1e-15;
+ 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(wlen_ran && func_name);
+ ASSERT(range && range[0] < range[1]);
+ ASSERT(0 <= r && r < 1);
+
+ /* Convert the wavelength range in meters */
+ range_m[0] = range[0] * 1.0e-9;
+ range_m[1] = range[1] * 1.0e-9;
+
+ /* 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], wlen_ran->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, wlen_ran->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_M
+ || fabs(bf_prev - bf) < EPSILON_BF)
+ break;
+
+ lambda_m_prev = lambda_m;
+ bf_prev = bf;
+ }
+ if(i >= MAX_ITER) {
+ htrdr_log_warn(wlen_ran->htrdr,
+ "%s: could not sample a wavelength in the range [%g, %g] nanometers "
+ "for the reference temperature %g Kelvin.\n",
+ func_name, SPLIT2(range), wlen_ran->ref_temperature);
+ }
+
+ if(pdf) {
+ const double Tref = wlen_ran->ref_temperature;
+ const double B_lambda = planck(lambda_m, lambda_m, Tref);
+ const double B_mean = planck(range_m[0], range_m[1], Tref);
+ *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0]));
+ }
+
+ return lambda_m * 1.0e9; /* Convert in nanometers */
+}
+
+static double
+wlen_ran_sample_discrete
+ (const struct htrdr_ran_wlen* wlen_ran,
+ const double r0,
+ const double r1,
+ const char* func_name,
+ double* pdf)
+{
+ const double* cdf = NULL;
+ const double* find = NULL;
+ double r0_next = nextafter(r0, DBL_MAX);
+ double band_range[2];
+ double lambda = 0;
+ double pdf_continue = 0;
+ double pdf_band = 0;
+ size_t cdf_length = 0;
+ size_t i;
+ ASSERT(wlen_ran && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
+ ASSERT(0 <= r0 && r0 < 1);
+ ASSERT(0 <= r1 && r1 < 1);
+ (void)func_name;
+ (void)pdf_band;
+
+ cdf = darray_double_cdata_get(&wlen_ran->cdf);
+ cdf_length = darray_double_size_get(&wlen_ran->cdf);
+ ASSERT(cdf_length > 0);
+
+ /* Use r_next rather than r0 in order to find the first entry that is not less
+ * than *or equal* to r0 */
+ 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] > r0 && (!i || cdf[i-1] <= r0));
+
+ band_range[0] = wlen_ran->range[0] + (double)i*wlen_ran->band_len;
+ band_range[1] = band_range[0] + wlen_ran->band_len;
+
+ /* Fetch the pdf of the sampled band */
+ pdf_band = darray_double_cdata_get(&wlen_ran->pdf)[i];
+
+ /* Uniformly sample a wavelength in the sampled band */
+ lambda = band_range[0] + (band_range[1] - band_range[0]) * r1;
+ pdf_continue = 1.0 / ((band_range[1] - band_range[0])*1.e-9);
+
+ if(pdf) {
+ *pdf = pdf_band * pdf_continue;
+ }
+
+ return lambda;
+}
+
+static void
+release_wlen_ran(ref_T* ref)
+{
+ struct htrdr_ran_wlen* wlen_ran = NULL;
+ ASSERT(ref);
+ wlen_ran = CONTAINER_OF(ref, struct htrdr_ran_wlen, ref);
+ darray_double_release(&wlen_ran->cdf);
+ darray_double_release(&wlen_ran->pdf);
+ MEM_RM(wlen_ran->htrdr->allocator, wlen_ran);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_ran_wlen_create
+ (struct htrdr* htrdr,
+ /* range must be included in [200,1000] nm for solar or in [1000, 100000]
+ * nanometers for longwave (thermal)*/
+ const double range[2],
+ const size_t nbands, /* # bands used to discretized CDF */
+ const double ref_temperature,
+ struct htrdr_ran_wlen** out_wlen_ran)
+{
+ struct htrdr_ran_wlen* wlen_ran = NULL;
+ double min_band_len = 0;
+ res_T res = RES_OK;
+ ASSERT(htrdr && range && out_wlen_ran && ref_temperature > 0);
+ ASSERT(ref_temperature > 0);
+ ASSERT(range[0] <= range[1]);
+
+ wlen_ran = MEM_CALLOC(htrdr->allocator, 1, sizeof(*wlen_ran));
+ if(!wlen_ran) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: could not allocate longwave random variate data structure -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ ref_init(&wlen_ran->ref);
+ wlen_ran->htrdr = htrdr;
+ darray_double_init(htrdr->allocator, &wlen_ran->cdf);
+ darray_double_init(htrdr->allocator, &wlen_ran->pdf);
+
+ wlen_ran->range[0] = range[0];
+ wlen_ran->range[1] = range[1];
+ wlen_ran->ref_temperature = ref_temperature;
+ wlen_ran->nbands = nbands;
+
+ min_band_len = compute_sky_min_band_len(wlen_ran->htrdr->sky, wlen_ran->range);
+
+ if(nbands == HTRDR_WLEN_RAN_CONTINUE) {
+ wlen_ran->band_len = 0;
+ } else {
+ wlen_ran->band_len = (range[1] - range[0]) / (double)wlen_ran->nbands;
+
+ /* Adjust the band length to ensure that each sky spectral interval is
+ * overlapped by at least one band */
+ if(wlen_ran->band_len > min_band_len) {
+ wlen_ran->band_len = min_band_len;
+ wlen_ran->nbands = (size_t)ceil((range[1] - range[0]) / wlen_ran->band_len);
+ }
+
+ res = setup_wlen_ran_cdf(wlen_ran, FUNC_NAME);
+ if(res != RES_OK) goto error;
+ }
+
+ htrdr_log(htrdr, "Spectral interval defined on [%g, %g] nanometers.\n",
+ range[0], range[1]);
+
+exit:
+ *out_wlen_ran = wlen_ran;
+ return res;
+error:
+ if(wlen_ran) htrdr_ran_wlen_ref_put(wlen_ran);
+ goto exit;
+}
+
+void
+htrdr_ran_wlen_ref_get(struct htrdr_ran_wlen* wlen_ran)
+{
+ ASSERT(wlen_ran);
+ ref_get(&wlen_ran->ref);
+}
+
+void
+htrdr_ran_wlen_ref_put(struct htrdr_ran_wlen* wlen_ran)
+{
+ ASSERT(wlen_ran);
+ ref_put(&wlen_ran->ref, release_wlen_ran);
+}
+
+double
+htrdr_ran_wlen_sample
+ (const struct htrdr_ran_wlen* wlen_ran,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ ASSERT(wlen_ran);
+ if(wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE) { /* Discrete */
+ return wlen_ran_sample_discrete(wlen_ran, r0, r1, FUNC_NAME, pdf);
+ } else if(eq_eps(wlen_ran->range[0], wlen_ran->range[1], 1.e-6)) {
+ if(pdf) *pdf = 1;
+ return wlen_ran->range[0];
+ } else { /* Continue */
+ return wlen_ran_sample_continue
+ (wlen_ran, r0, wlen_ran->range, FUNC_NAME, pdf);
+ }
+}
+
diff --git a/src/htrdr_ran_wlen.h b/src/htrdr_ran_wlen.h
@@ -0,0 +1,54 @@
+/* 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/>. */
+
+#ifndef HTRDR_RAN_WLEN_H
+#define HTRDR_RAN_WLEN_H
+
+#include <rsys/rsys.h>
+
+#define HTRDR_WLEN_RAN_CONTINUE 0
+
+struct htrdr;
+struct htrdr_ran_wlen;
+
+extern LOCAL_SYM res_T
+htrdr_ran_wlen_create
+ (struct htrdr* htrdr,
+ const double range[2],
+ /* # bands used to discretisze the spectral domain. HTRDR_WLEN_RAN_CONTINUE
+ * <=> no discretisation */
+ const size_t nbands, /* Hint on #bands used to discretised th CDF */
+ const double ref_temperature, /* Reference temperature */
+ struct htrdr_ran_wlen** wlen_ran);
+
+extern LOCAL_SYM void
+htrdr_ran_wlen_ref_get
+ (struct htrdr_ran_wlen* wlen_ran);
+
+extern LOCAL_SYM void
+htrdr_ran_wlen_ref_put
+ (struct htrdr_ran_wlen* wlen_ran);
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_ran_wlen_sample
+ (const struct htrdr_ran_wlen* wlen_ran,
+ const double r0, /* Canonical number in [0, 1[ */
+ const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* May be NULL */
+
+#endif /* HTRDR_RAN_WLEN_H */
+
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -35,20 +35,20 @@ struct htrdr_estimate {
};
static const struct htrdr_estimate HTRDR_ESTIMATE_NULL;
-struct htrdr_pixel_sw {
- struct htrdr_estimate X; /* In W/m^2/sr */
- struct htrdr_estimate Y; /* In W/m^2/sr */
- struct htrdr_estimate Z; /* In W/m^2/sr */
+struct htrdr_pixel_xwave {
+ struct htrdr_estimate radiance; /* In W/m^2/sr */
+ struct htrdr_estimate radiance_temperature; /* In K */
struct htrdr_accum time; /* In microseconds */
};
-static const struct htrdr_pixel_sw HTRDR_PIXEL_SW_NULL;
+static const struct htrdr_pixel_xwave HTRDR_PIXEL_XWAVE_NULL;
-struct htrdr_pixel_lw {
- struct htrdr_estimate radiance; /* In W/m^2/sr */
- struct htrdr_estimate radiance_temperature; /* In K */
+struct htrdr_pixel_image {
+ struct htrdr_estimate X; /* In W/m^2/sr */
+ struct htrdr_estimate Y; /* In W/m^2/sr */
+ struct htrdr_estimate Z; /* In W/m^2/sr */
struct htrdr_accum time; /* In microseconds */
};
-static const struct htrdr_pixel_lw HTRDR_PIXEL_LW_NULL;
+static const struct htrdr_pixel_image HTRDR_PIXEL_IMAGE_NULL;
/* Forward declarations */
struct htrdr;
@@ -56,6 +56,7 @@ struct htrdr_camera;
struct s3d_hit;
struct ssp_rng;
+/* Return the shortwave radiance in W/m^2/sr/m */
extern LOCAL_SYM double
htrdr_compute_radiance_sw
(struct htrdr* htrdr,
@@ -67,6 +68,7 @@ htrdr_compute_radiance_sw
const size_t iband, /* Index of the spectral band */
const size_t iquad); /* Index of the quadrature point into the band */
+/* Return the longwave radiance in W/m^2/sr/m */
extern LOCAL_SYM double
htrdr_compute_radiance_lw
(struct htrdr* htrdr,
@@ -118,4 +120,38 @@ htrdr_accum_get_estimation
}
}
+static INLINE size_t
+htrdr_spectral_type_get_pixsz(const enum htrdr_spectral_type type)
+{
+ size_t sz = 0;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ sz = sizeof(struct htrdr_pixel_xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ sz = sizeof(struct htrdr_pixel_image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return sz;
+}
+
+static INLINE size_t
+htrdr_spectral_type_get_pixal(const enum htrdr_spectral_type type)
+{
+ size_t al = 0;
+ switch(type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ al = ALIGNOF(struct htrdr_pixel_xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ al = ALIGNOF(struct htrdr_pixel_image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return al;
+}
+
#endif /* HTRDR_SOLVE_H */
diff --git a/src/htrdr_spectral.c b/src/htrdr_spectral.c
@@ -0,0 +1,79 @@
+/* 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_spectral.h"
+
+res_T
+brightness_temperature
+ (struct htrdr* htrdr,
+ const double lambda_min,
+ const double lambda_max,
+ const double radiance, /* In W/m2/sr/m */
+ double* temperature)
+{
+ const size_t MAX_ITER = 100;
+ const double epsilon_T = 1e-4; /* In K */
+ const double epsilon_B = radiance * 1e-8;
+ double T, T0, T1, T2;
+ double B, B0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(temperature && lambda_min <= lambda_max);
+
+ /* Search for a brightness temperature whose radiance is greater than or
+ * equal to the estimated radiance */
+ T2 = 200;
+ FOR_EACH(i, 0, MAX_ITER) {
+ const double B2 = planck(lambda_min, lambda_max, T2);
+ if(B2 >= radiance) break;
+ T2 *= 2;
+ }
+ if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+
+ B0 = T0 = T1 = 0;
+ FOR_EACH(i, 0, MAX_ITER) {
+ T = (T1+T2)*0.5;
+ B = planck(lambda_min, lambda_max, T);
+
+ if(B < radiance) {
+ T1 = T;
+ } else {
+ T2 = T;
+ }
+
+ if(fabs(T-T0) < epsilon_T || fabs(B-B0) < epsilon_B)
+ break;
+
+ T0 = T;
+ B0 = B;
+ }
+ if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+
+ *temperature = T;
+
+exit:
+ return res;
+error:
+ htrdr_log_err(htrdr,
+ "Could not compute the brightness temperature for the estimated radiance %g "
+ "averaged over [%g, %g] nanometers.\n",
+ radiance,
+ lambda_min*1e9,
+ lambda_max*1e9);
+ goto exit;
+}
+
diff --git a/src/htrdr_spectral.h b/src/htrdr_spectral.h
@@ -0,0 +1,137 @@
+/* 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/>. */
+
+#ifndef HTRDR_SPECTRAL_H
+#define HTRDR_SPECTRAL_H
+
+#include <rsys/rsys.h>
+#include <rsys/math.h> /* PI support */
+
+#define HTRDR_SUN_TEMPERATURE 5778 /* In K */
+#define HTRDR_DEFAULT_LW_REF_TEMPERATURE 290 /* Default longwave temperature in K */
+
+enum htrdr_spectral_type {
+ HTRDR_SPECTRAL_LW, /* Longwave */
+ HTRDR_SPECTRAL_SW, /* Shortwave */
+ HTRDR_SPECTRAL_SW_CIE_XYZ /* Shortwave wrt the CIE XYZ tristimulus */
+};
+
+struct htrdr;
+
+static INLINE double
+wiebelt(const double v)
+{
+ int m;
+ double w, v2, v4;
+ /*.153989717364e+00;*/
+ const double fifteen_over_pi_power_4 = 15.0/(PI*PI*PI*PI);
+ const double z0 = 1.0/3.0;
+ const double z1 = 1.0/8.0;
+ const double z2 = 1.0/60.0;
+ const double z4 = 1.0/5040.0;
+ const double z6 = 1.0/272160.0;
+ const double z8 = 1.0/13305600.0;
+
+ if(v >= 2.) {
+ w = 0;
+ for(m=1; m<6 ;m++)
+ w+=exp(-m*v)/(m*m*m*m) * (((m*v+3)*m*v+6)*m*v+6);
+ w = w * fifteen_over_pi_power_4;
+ } else {
+ v2 = v*v;
+ v4 = v2*v2;
+ w = z0 - z1*v + z2*v2 - z4*v2*v2 + z6*v4*v2 - z8*v4*v4;
+ w = 1. - fifteen_over_pi_power_4*v2*v*w;
+ }
+ ASSERT(w >= 0.0 && w <= 1.0);
+ return w;
+}
+
+static INLINE double
+blackbody_fraction
+ (const double lambda0, /* In meters */
+ const double lambda1, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ const double C2 = 1.43877735e-2; /* m.K */
+ double x0 = C2 / lambda0;
+ double x1 = C2 / lambda1;
+ double v0 = x0 / temperature;
+ double v1 = x1 / temperature;
+ double w0 = wiebelt(v0);
+ double w1 = wiebelt(v1);
+ return w1 - w0;
+}
+
+/* Return the Planck value in W/m^2/sr/m at a given wavelength */
+static INLINE double
+planck_monochromatic
+ (const double lambda, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ const double c = 299792458; /* m/s */
+ const double h = 6.62607015e-34; /* J.s */
+ const double k = 1.380649e-23; /* J/K */
+ const double lambda2 = lambda*lambda;
+ const double lambda5 = lambda2*lambda2*lambda;
+ const double B = ((2.0 * h * c*c) / lambda5) /* W/m^2/sr/m */
+ / (exp(h*c/(lambda*k*temperature))-1.0);
+ ASSERT(temperature > 0);
+ return B;
+}
+
+/* Return the average Planck value in W/m^2/sr/m over the
+ * [lambda_min, lambda_max] interval. */
+static INLINE double
+planck_interval
+ (const double lambda_min, /* In meters */
+ const double lambda_max, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ const double T2 = temperature*temperature;
+ const double T4 = T2*T2;
+ const double BOLTZMANN_CONSTANT = 5.6696e-8; /* W/m^2/K^4 */
+ ASSERT(lambda_min < lambda_max && temperature > 0);
+ return blackbody_fraction(lambda_min, lambda_max, temperature)
+ * BOLTZMANN_CONSTANT * T4 / (PI * (lambda_max-lambda_min)); /* In W/m^2/sr/m */
+}
+
+/* Invoke planck_monochromatic or planck_interval whether the submitted
+ * interval is null or not, respectively. The returned value is in W/m^2/sr/m */
+static FINLINE double
+planck
+ (const double lambda_min, /* In meters */
+ const double lambda_max, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ ASSERT(lambda_min <= lambda_max && temperature > 0);
+ if(lambda_min == lambda_max) {
+ return planck_monochromatic(lambda_min, temperature);
+ } else {
+ return planck_interval(lambda_min, lambda_max, temperature);
+ }
+}
+
+extern LOCAL_SYM res_T
+brightness_temperature
+ (struct htrdr* htrdr,
+ const double lambda_min, /* In meters */
+ const double lambda_max, /* In meters */
+ /* Averaged over [lambda_min, lambda_max]. In W/m^2/sr/m */
+ const double radiance,
+ double* temperature);
+
+#endif /* HTRDR_SPECTRAL_H */
diff --git a/src/htrdr_sun.c b/src/htrdr_sun.c
@@ -20,23 +20,17 @@
#include <rsys/algorithm.h>
#include <rsys/double33.h>
-#include <rsys/dynamic_array_double.h>
#include <rsys/ref_count.h>
#include <rsys/math.h>
#include <star/ssp.h>
struct htrdr_sun {
- /* Short wave radiance in W.m^-2.sr^-1, for each spectral interval */
- struct darray_double radiances_sw;
-
- /* Short wave spectral interval boundaries, in cm^-1 */
- struct darray_double wavenumbers_sw;
-
double half_angle; /* In radian */
double cos_half_angle;
double solid_angle; /* In sr; solid_angle = 2*PI*(1 - cos(half_angle)) */
double frame[9];
+ double temperature; /* In K */
ref_T ref;
struct htrdr* htrdr;
@@ -51,8 +45,6 @@ release_sun(ref_T* ref)
struct htrdr_sun* sun;
ASSERT(ref);
sun = CONTAINER_OF(ref, struct htrdr_sun, ref);
- darray_double_release(&sun->radiances_sw);
- darray_double_release(&sun->wavenumbers_sw);
MEM_RM(sun->htrdr->allocator, sun);
}
@@ -62,63 +54,25 @@ release_sun(ref_T* ref)
res_T
htrdr_sun_create(struct htrdr* htrdr, struct htrdr_sun** out_sun)
{
- const double incoming_flux_sw[] = { /* In W.m^-2 */
- 12.793835026999544, 12.109561093845551, 20.365091338928245,
- 23.729742422870157, 22.427697221814142, 55.626612361454150,
- 102.93146523363953, 24.293596268358986, 345.73659325842243,
- 218.18441435866691, 347.18437832794524, 129.49426803812202,
- 50.146977730963876, 3.1197193425713365
- };
- const double wavenumbers_sw[] = { /* In cm^-1 */
- 820.000, 2600.00, 3250.00, 4000.00, 4650.00,
- 5150.00, 6150.00, 7700.00, 8050.00, 12850.0,
- 16000.0, 22650.0, 29000.0, 38000.0, 49999.0
- };
-
- const size_t nspectral_intervals = sizeof(incoming_flux_sw)/sizeof(double);
const double main_dir[3] = {0, 0, 1}; /* Default main sun direction */
struct htrdr_sun* sun = NULL;
- size_t i;
res_T res = RES_OK;
ASSERT(htrdr && out_sun);
- ASSERT(sizeof(wavenumbers_sw)/sizeof(double) == nspectral_intervals+1);
sun = MEM_CALLOC(htrdr->allocator, 1, sizeof(*sun));
if(!sun) {
htrdr_log_err(htrdr, "could not allocate the sun data structure.\n");
res = RES_MEM_ERR;
- goto error;
-}
+ goto error;
+ }
ref_init(&sun->ref);
sun->htrdr = htrdr;
- darray_double_init(htrdr->allocator, &sun->radiances_sw);
- darray_double_init(htrdr->allocator, &sun->wavenumbers_sw);
sun->half_angle = 4.6524e-3;
+ sun->temperature = 5778;
sun->cos_half_angle = cos(sun->half_angle);
sun->solid_angle = 2*PI*(1-sun->cos_half_angle);
d33_basis(sun->frame, main_dir);
- res = darray_double_resize(&sun->radiances_sw, nspectral_intervals);
- if(res != RES_OK) {
- htrdr_log_err(htrdr,
- "could not allocate the list of per spectral band radiance of the sun.\n");
- goto error;
- }
- res = darray_double_resize(&sun->wavenumbers_sw, nspectral_intervals+1);
- if(res != RES_OK) {
- htrdr_log_err(htrdr,
- "could not allocate the list of spectral band boundaries of the sun.\n");
- goto error;
- }
-
- FOR_EACH(i, 0, darray_double_size_get(&sun->radiances_sw)) {
- /* Convert the incoming flux in radiance */
- darray_double_data_get(&sun->radiances_sw)[i] = incoming_flux_sw[i] / PI;
- }
- FOR_EACH(i, 0, darray_double_size_get(&sun->wavenumbers_sw)) {
- darray_double_data_get(&sun->wavenumbers_sw)[i] = wavenumbers_sw[i];
- }
-
exit:
*out_sun = sun;
return res;
@@ -170,36 +124,9 @@ htrdr_sun_get_solid_angle(const struct htrdr_sun* sun)
}
double
-htrdr_sun_get_radiance(const struct htrdr_sun* sun, const double wavelength)
+htrdr_sun_get_radiance(const struct htrdr_sun* sun, const double wlen/*In nm*/)
{
- const double* wavenumbers;
- const double wnum = wavelength_to_wavenumber(wavelength);
- const double* wnum_upp;
- size_t nwavenumbers;
- size_t ispectral_band;
- ASSERT(sun && wavelength > 0);
-
- wavenumbers = darray_double_cdata_get(&sun->wavenumbers_sw);
- nwavenumbers = darray_double_size_get(&sun->wavenumbers_sw);
- ASSERT(nwavenumbers);
-
- if(wnum < wavenumbers[0] || wnum > wavenumbers[nwavenumbers-1]) {
- htrdr_log_warn(sun->htrdr,
- "the submitted wavelength is outside the sun spectrum.\n");
- }
-
- wnum_upp = search_lower_bound
- (&wnum, wavenumbers, nwavenumbers, sizeof(double), cmp_dbl);
-
- if(!wnum_upp) { /* Clamp to the upper spectral band */
- ispectral_band = nwavenumbers - 2;
- ASSERT(ispectral_band == darray_double_size_get(&sun->radiances_sw)-1);
- } else if(wnum_upp == wavenumbers) { /* Clamp to the lower spectral band */
- ispectral_band = 0;
- } else {
- ispectral_band = (size_t)(wnum_upp - wavenumbers - 1);
- }
- return darray_double_cdata_get(&sun->radiances_sw)[ispectral_band];
+ return planck_monochromatic(wlen*1.e-9/*From nm to m*/, sun->temperature);
}
int
@@ -214,59 +141,3 @@ htrdr_sun_is_dir_in_solar_cone(const struct htrdr_sun* sun, const double dir[3])
return dot >= sun->cos_half_angle;
}
-size_t
-htrdr_sun_get_spectral_bands_count(const struct htrdr_sun* sun)
-{
- ASSERT(sun);
- return darray_double_size_get(&sun->radiances_sw);
-}
-
-void
-htrdr_sun_get_spectral_band_bounds
- (const struct htrdr_sun* sun,
- const size_t ispectral_band,
- double bounds[2])
-{
- size_t iband_wlen;
- const double* wnums;
- ASSERT(sun && ispectral_band < htrdr_sun_get_spectral_bands_count(sun));
- ASSERT(ispectral_band + 1 < darray_double_size_get(&sun->wavenumbers_sw));
-
- /* The spectral bands are internally defined with wavenumbers, while the sun
- * API uses wavelengths. But since wavelengths are inversely proportionnal to
- * wavenumbers, one have to reversed the spectral band index */
- iband_wlen = htrdr_sun_get_spectral_bands_count(sun) - 1 - ispectral_band;
- wnums = darray_double_cdata_get(&sun->wavenumbers_sw);
- bounds[1] = wavenumber_to_wavelength(wnums[iband_wlen+0]);
- bounds[0] = wavenumber_to_wavelength(wnums[iband_wlen+1]);
- ASSERT(bounds[0] <= bounds[1]);
-}
-
-void
-htrdr_sun_get_CIE_XYZ_spectral_bands_range
- (const struct htrdr_sun* sun, size_t band_range[2])
-{
- /* Bounds of the X, Y and Z functions as defined by the CIE */
- const double nu_min = wavelength_to_wavenumber(SW_WAVELENGTH_MIN);/*In cm^-1*/
- const double nu_max = wavelength_to_wavenumber(SW_WAVELENGTH_MAX);/*In cm^-1*/
- const double* wnums = NULL;
- size_t iband_low_nu = SIZE_MAX;
- size_t iband_upp_nu = SIZE_MAX;
- size_t i;
- ASSERT(sun && band_range);
-
- wnums = darray_double_cdata_get(&sun->wavenumbers_sw);
-
- /* Perform a linear search since the number of spectral bands is small */
- FOR_EACH(i, 0, htrdr_sun_get_spectral_bands_count(sun)) {
- if(wnums[i] <= nu_min && nu_min < wnums[i+1]) iband_low_nu = i;
- if(wnums[i] <= nu_max && nu_max < wnums[i+1]) iband_upp_nu = i;
- }
- ASSERT(iband_low_nu <= iband_upp_nu);
-
- /* Transform the band index from "wavenumber space" to "wavelength space" */
- band_range[0] = htrdr_sun_get_spectral_bands_count(sun) - 1 - iband_upp_nu;
- band_range[1] = htrdr_sun_get_spectral_bands_count(sun) - 1 - iband_low_nu;
- ASSERT(band_range[0] <= band_range[1]);
-}
-
diff --git a/src/htrdr_sun.h b/src/htrdr_sun.h
@@ -54,7 +54,7 @@ extern LOCAL_SYM double
htrdr_sun_get_solid_angle
(const struct htrdr_sun* sun);
-extern LOCAL_SYM double /* W.sr^-1.m^-2 */
+extern LOCAL_SYM double /* W/m^2/sr/m */
htrdr_sun_get_radiance
(const struct htrdr_sun* sun,
const double wavelength);
@@ -64,21 +64,4 @@ htrdr_sun_is_dir_in_solar_cone
(const struct htrdr_sun* sun,
const double dir[3]);
-extern LOCAL_SYM size_t
-htrdr_sun_get_spectral_bands_count
- (const struct htrdr_sun* sun);
-
-extern LOCAL_SYM void
-htrdr_sun_get_spectral_band_bounds
- (const struct htrdr_sun* sun,
- const size_t ispectral_band,
- double bounds[2]); /* Lower and upper wavelength in nanometer */
-
-/* Return the ranges of the spectral bands where the CIE XYZ color space is
- * defined. CIE XYZ in [band_range[0], band_range[1]] */
-extern LOCAL_SYM void
-htrdr_sun_get_CIE_XYZ_spectral_bands_range
- (const struct htrdr_sun* sun,
- size_t band_range[2]);
-
#endif /* HTRDR_SUN_H */
