commit 6faf2caafe28c02ec18717002952d874deda22d3
parent a1b7875594d4c9d04852a54e71b1b9a45e58c8de
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 19 Nov 2020 11:56:02 +0100
Merge branch 'feature_flux_map' into develop
Diffstat:
27 files changed, 2584 insertions(+), 1836 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -28,7 +28,7 @@ find_package(AW 2.0 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(RSys 0.11 REQUIRED)
find_package(Star3D 0.7.1 REQUIRED)
find_package(StarSF 0.6 REQUIRED)
find_package(StarSP 0.8 REQUIRED)
@@ -67,10 +67,15 @@ set(HTRDR_ARGS_DEFAULT_CAMERA_POS "0,0,0")
set(HTRDR_ARGS_DEFAULT_CAMERA_TGT "0,1,0")
set(HTRDR_ARGS_DEFAULT_CAMERA_UP "0,0,1")
set(HTRDR_ARGS_DEFAULT_CAMERA_FOV "70")
+set(HTRDR_ARGS_DEFAULT_RECTANGLE_POS "0,0,0")
+set(HTRDR_ARGS_DEFAULT_RECTANGLE_TGT "0,0,1")
+set(HTRDR_ARGS_DEFAULT_RECTANGLE_UP "0,1,0")
+set(HTRDR_ARGS_DEFAULT_RECTANGLE_SZ "1,1")
set(HTRDR_ARGS_DEFAULT_IMG_WIDTH "320")
set(HTRDR_ARGS_DEFAULT_IMG_HEIGHT "240")
set(HTRDR_ARGS_DEFAULT_IMG_SPP "1")
set(HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD "1")
+set(HTRDR_ARGS_DEFAULT_SKY_MTL_NAME "\"air\"")
configure_file(${HTRDR_SOURCE_DIR}/htrdr_version.h.in
${CMAKE_CURRENT_BINARY_DIR}/htrdr_version.h @ONLY)
@@ -94,14 +99,14 @@ set(HTRDR_FILES_SRC
htrdr_cie_xyz.c
htrdr_compute_radiance_sw.c
htrdr_compute_radiance_lw.c
- htrdr_draw_radiance.c
+ htrdr_draw_map.c
htrdr_grid.c
htrdr_ground.c
- htrdr_interface.c
htrdr_main.c
- htrdr_mtl.c
+ htrdr_materials.c
htrdr_ran_wlen.c
htrdr_rectangle.c
+ htrdr_sensor.c
htrdr_slab.c
htrdr_spectral.c
htrdr_sun.c)
@@ -114,9 +119,10 @@ set(HTRDR_FILES_INC
htrdr_grid.h
htrdr_ground.h
htrdr_interface.h
- htrdr_mtl.h
+ htrdr_materials.h
htrdr_ran_wlen.h
htrdr_rectangle.h
+ htrdr_sensor.h
htrdr_slab.h
htrdr_spectral.h
htrdr_sun.h
diff --git a/doc/htrdr-image.5.txt b/doc/htrdr-image.5.txt
@@ -30,33 +30,41 @@ 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
-(*-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
+If the image is a regular rendering in the visible part of the spectrum (*-C*
+_camera_ and *-s* _cie_xyz_ options 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
+If the image is an infrared rendering (*-C* _camera_ and *-s*
+*lw*=_wlen-min_,_wlen_max_ options 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 (*-C* _camera_ and *-s*
+*sw*=_wlen-min_,wlen-max_ options 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.
+For flux computations (*-p* _rectangle_ option in *htrdr*(1)), the first and
+second pixel component stores the expected value and the standard error of the
+pixel flux in W/m^2 for the part of the pixel that is outside any geometry. As
+previously, the seventh and eighth pixel components store the estimate of the
+radiative path computation time in microseconds and its standard error. The
+remaining components, i.e. the components 3 to 6, are unused.
+
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-materials.5.txt b/doc/htrdr-materials.5.txt
@@ -26,11 +26,14 @@ DESCRIPTION
-----------
A *htrdr-materials* file lists the materials that can be used by the ground
geometry provided through a *htrdr-obj*(5) file to the *htrdr*(1) program.
-Each line of the file gives the name of the material and the path toward the
-*mrumtl*(5) file storing the spectral properties of its associated
-Bidirectional Reflectance Distribution Function. The material name can be
-composed of any characters excepted for spaces and tabulations. The path toward
-the *mrumtl*(5) file must be a valid path relative to the file system.
+Each line of the file gives the name of the material. For opaque materials the
+material name is followed by the path toward the *mrumtl*(5) file storing the
+spectral properties of its associated Bidirectional Reflectance Distribution
+Function. Furthermore, the temperature of the material must be defined too.
+
+The material name can be composed of any characters excepted for spaces and
+tabulations. The path toward the *mrumtl*(5) file must be a valid path
+relative to the file system.
Characters behind the hash mark (#) are considered as comments and are thus
ignored. Empty lines, i.e. lines without any characters or composed of spaces
@@ -43,22 +46,28 @@ GRAMMAR
-------
<htrdr-materials> ::= <material>
[ <material> ... ]
-<material> ::= <name> <mrumtl-path>
+<material> ::= <name> <properties>
<name> ::= STRING
+<properties> ::= none | <mrumtl-path> <temperature>
<mrumtl-path> ::= PATH
+<temperature> ::= REAL # In Kelvin
-------
EXAMPLE
-------
-The following file lists 2 materials. The first one named *grass* has its
+The following file lists 3 materials. The first one named *grass* has its
spectral BRDF defines in the *A001.mrumtl* file. The second one is named
-*sand* and has its spectral properties saved in the *B002.mrumtl* file.
+*sand* and has its spectral properties saved in the *B002.mrumtl* file. Both
+materials have a temperature of 300 Kelvin. The last material is a
+semi-transparent material named *air* with no additionnal properties defined
+in this file.
[verse]
-------
-grass /opt/materials/A001.mrumtl
-sand /opt/materials/B002.mrumtl
+grass /opt/materials/A001.mrumtl 300
+sand /opt/materials/B002.mrumtl 300
+air none
-------
SEE ALSO
diff --git a/doc/htrdr.1.txt.in b/doc/htrdr.1.txt.in
@@ -20,7 +20,7 @@ htrdr(1)
NAME
----
-htrdr - image renderer of cloudy atmospheres
+htrdr - simulate radiative transfert in cloudy atmospheres
SYNOPSIS
--------
@@ -29,48 +29,55 @@ SYNOPSIS
DESCRIPTION
-----------
-*htrdr* is an image renderer of scenes composed of an atmospheric gas mixture,
-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)
-file (*-c* _clouds_), and the optical properties of water droplets that have
-been obtained from a Mie code and formatted according to the *htmie*(5) file
-format (*-m* _mie_). The user also has to provide: the characteristics of the
-simulated camera (*-C* _camera_), the sensor definition (*-i* _image_), and the
-position of the sun (*-D* _azimuth_,_elevation_). It is also possible to
-provide an *htrdr-obj*(5) file representing the ground geometry (*-g* _ground_)
-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.
-
-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. The two other ways consist in explicitly defining the
-longwave or shortwave spectral range to handle and continuously sampling a
-wavelength in this range. 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.
+*htrdr* simulate radiative transfert in scenes composed of an atmospheric gas
+mixture, 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 sensor, directed into the atmosphere, taking into account light
+absorption and scattering phenomena.
+
+Radiative transfert can be evaluated 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) file (*-c*
+_clouds_), and the optical properties of water droplets that have been obtained
+from a Mie code and formatted according to the *htmie*(5) file format (*-m*
+_mie_). The user also has to set the position of the sun (*-D*
+_azimuth_,_elevation_), the sensor type (*-C* _camera_ or *-p* _rectangle_) and
+its definition (*-i* _image_). It is also possible to provide an *htrdr-obj*(5)
+file representing the ground geometry (*-g* _ground_) 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.
+
+Two types of sensor are supported by *htrdr*. The camera (*-C* _camera_) is
+used to render an image of the scene from the given point of view while the
+rectangle sensor (*-p* _rectangle_) is used to compute a flux map.
+
+Spectral dimension can be integrated in many ways (*-s* option). When rendering
+an image (*-C* _camera_), the computation is by default 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. The two other ways consist
+in explicitly defining the longwave or shortwave spectral range to handle and
+continuously sampling a wavelength in this range. Actually longwave and
+shortwave are keywords that mean that the source of radiation is whether
+external or internal to the medium, respectively. In shortwave rendering, only
+the pixel radiance is evaluated and stored in the output image. For longwave
+rendering this estimated radiance is then converted to its brightness
+temperature and both are saved in the image. When computing a flux map (*-p*
+_rectangle_), the per pixel flux is saved into the output map whether spectral
+domain is longwave or shortwave.
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.
+in W/sr/m^2/nm. The fluxes are saved in W/m^2. 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
@@ -137,7 +144,7 @@ OPTIONS
List short help and exit.
*-i* <__image-parameter__:...>::
- Define the image to render. Available image parameters are:
+ Define the sensor array. Available image parameters are:
**def**=**_width_**x**_height_**;;
Definition of the image. By default the image definition is
@@ -147,10 +154,9 @@ 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 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@.
+ XYZ color space. In shortwave/longwave rendering or flux computation, only
+ one set of _samples-count_ is used. By default, *spp* is set to
+ @HTRDR_ARGS_DEFAULT_IMG_SPP@.
*-R*::
Infinitely repeat the _ground_ along the X and Y axis.
@@ -164,6 +170,10 @@ OPTIONS
*-m* _mie_::
Path toward a *htmie*(5) file defining the optical properties of water
droplets.
+
+*-n* _sky-mtl_::
+ Name of the material representing the sky in the *htrdr-materials*(5) file.
+ By default, _sky-mtl_ is @HTRDR_ARGS_DEFAULT_SKY_MTL_NAME@.
*-O* _cache_::
File used to cache the sky data. If the _cache_ file does not exists, it is
@@ -182,8 +192,29 @@ OPTIONS
File where *htrdr* writes its _output_ data. If not defined, write results to
standard output.
+*-p* <__rectangle-parameter__:...>::
+ Switch in flux map computation. The flux is computed for the part of the
+ sensor that is outside any geometry. The rectangular sensor onto which the
+ flux is integrated is defined by the following parameters:
+
+ **pos**=**_x_**,**_y_**,**_z_**;;
+ Position of the center of the rectangle. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_POS@}.
+
+ **tgt**=**_x_**,**_y_**,**_z_**;;
+ Position targeted by the rectangle, i.e. *tgt* - *pos* is the rectangle
+ normal. By default it is set to {@HTRDR_ARGS_DEFAULT_RECTANGLE_TGT@}.
+
+ **up**=**_x_**,**_y_**,**_z_**;;
+ Up vector of the rectangle. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_UP@}.
+
+ **sz**=**_width_**,**_height_**;;
+ Size of the rectangle. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_SZ@}.
+
*-s* <__spectral-parameter__:...>::
- define the type and the range of the spectral integration. Available spectral
+ Define the type and the range of the spectral integration. Available spectral
parameters are:
**cie_xyz**;;
@@ -276,8 +307,8 @@ PPM image [4]:
-f -o output
$ htpp -o image.ppm output
-Render the previous scene in infrared for the wavelengths in [9200, 10000]
-nanometers with a reference temperature of 300 Kelvin:
+Render the previous scene in infrared for the wavelengths in [*9200*, *10000*]
+nanometers with a reference temperature of *300* Kelvin:
$ htrdr -a gas.txt -m Mie.nc -g mountains.obj -R \
-M materials.mtl \
@@ -301,6 +332,21 @@ output and convert the resulting image in PPM before visualising it through the
-C pos=0,0,400:tgt=0,1,0:up=0,0,1 \
-i def=1280x720:spp=1024 | htpp | feh -
+Compute the downward flux for the shortwave interval [*350*, *4000*] nanometers
+on a square of *100* meters side positionned at the origin at *1* meter height.
+The resolution of the flux map is *500* by *500* pixels and *1000* realisations
+is used to estimate the flux per pixel. It is saved in the *flux_map.txt* file
+even though this file already exists:
+
+ $ htrdr -D0,90 -a gas.txt -m Mie.nc -g plane.obj -R \
+ -M materials.mtl \
+ -c clouds.htcp \
+ -O my_cache \
+ -p pos=0,0,1:tgt=0,0,2:up=0,1,0:sz=100,100 \
+ -i def=500x500:spp=1000 \
+ -s sw=350,4000 \
+ -f -o flux_map.txt
+
Write into *output* the data structures used to partition the clouds
properties. Use the *csplit*(1) Unix command to extract from *output* the list
of the generated grids and save each of them in a specific VTK file whose name
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -23,8 +23,9 @@
#include "htrdr_cie_xyz.h"
#include "htrdr_camera.h"
#include "htrdr_ground.h"
-#include "htrdr_mtl.h"
+#include "htrdr_materials.h"
#include "htrdr_ran_wlen.h"
+#include "htrdr_rectangle.h"
#include "htrdr_sun.h"
#include "htrdr_solve.h"
#include "htrdr_version.h"
@@ -125,6 +126,7 @@ dump_buffer
(struct htrdr* htrdr,
struct htrdr_buffer* buf,
struct htrdr_accum* time_acc, /* May be NULL */
+ struct htrdr_accum* flux_acc, /* May be NULL */
const char* stream_name,
FILE* stream)
{
@@ -135,8 +137,10 @@ dump_buffer
ASSERT(htrdr && buf && stream_name && stream);
(void)stream_name;
- pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
- pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
+ pixsz = htrdr_spectral_type_get_pixsz
+ (htrdr->spectral_type, htrdr->sensor.type);
+ pixal = htrdr_spectral_type_get_pixal
+ (htrdr->spectral_type, htrdr->sensor.type);
htrdr_buffer_get_layout(buf, &layout);
if(layout.elmt_size != pixsz || layout.alignment != pixal) {
@@ -148,25 +152,46 @@ dump_buffer
fprintf(stream, "%lu %lu\n", layout.width, layout.height);
if(time_acc) *time_acc = HTRDR_ACCUM_NULL;
+ if(flux_acc) *flux_acc = HTRDR_ACCUM_NULL;
FOR_EACH(y, 0, layout.height) {
FOR_EACH(x, 0, layout.width) {
struct htrdr_estimate pix_time = HTRDR_ESTIMATE_NULL;
const struct htrdr_accum* pix_time_acc = NULL;
- 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);
- fprintf(stream, "0 0 ");
+ if(htrdr->sensor.type == HTRDR_SENSOR_RECTANGLE) {
+ const struct htrdr_pixel_flux* pix = htrdr_buffer_at(buf, x, y);
+ struct htrdr_estimate flux = HTRDR_ESTIMATE_NULL;
+
+ if(pix->flux.nweights == 0) {
+ fprintf(stream, "0 0 0 0 0 0 ");
+ } else {
+ htrdr_accum_get_estimation(&pix->flux, &flux);
+ fprintf(stream, "%g %g 0 0 0 0 ", flux.E, flux.SE);
+
+ if(flux_acc) {
+ flux_acc->sum_weights += pix->flux.sum_weights;
+ flux_acc->sum_weights_sqr += pix->flux.sum_weights_sqr;
+ flux_acc->nweights += pix->flux.nweights;
+ }
+ }
pix_time_acc = &pix->time;
} else {
- 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);
- pix_time_acc = &pix->time;
+ 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);
+ fprintf(stream, "0 0 ");
+ pix_time_acc = &pix->time;
+
+ } else {
+ 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);
+ pix_time_acc = &pix->time;
+ }
}
htrdr_accum_get_estimation(pix_time_acc, &pix_time);
@@ -390,6 +415,47 @@ error:
goto exit;
}
+static res_T
+setup_sensor(struct htrdr* htrdr, const struct htrdr_args* args)
+{
+ double proj_ratio;
+ res_T res = RES_OK;
+ ASSERT(htrdr && args);
+
+ htrdr->sensor.type = args->sensor_type;
+
+ if(args->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ
+ && args->sensor_type != HTRDR_SENSOR_CAMERA) {
+ htrdr_log_err(htrdr, "the CIE 1931 XYZ spectral integration can be used "
+ "only with a camera sensor.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ switch(args->sensor_type) {
+ case HTRDR_SENSOR_CAMERA:
+ proj_ratio =
+ (double)args->image.definition[0]
+ / (double)args->image.definition[1];
+ res = htrdr_camera_create(htrdr, args->camera.pos, args->camera.tgt,
+ args->camera.up, proj_ratio, MDEG2RAD(args->camera.fov_y),
+ &htrdr->sensor.camera);
+ break;
+ case HTRDR_SENSOR_RECTANGLE:
+ res = htrdr_rectangle_create(htrdr, args->rectangle.sz,
+ args->rectangle.pos, args->rectangle.tgt, args->rectangle.up,
+ &htrdr->sensor.rectangle);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -400,7 +466,6 @@ htrdr_init
struct htrdr* htrdr)
{
struct htsky_args htsky_args = HTSKY_ARGS_DEFAULT;
- double proj_ratio;
double sun_dir[3];
double spectral_range[2];
const char* output_name = NULL;
@@ -430,6 +495,7 @@ htrdr_init
htrdr->grid_max_definition[2] = args->grid_max_definition[2];
htrdr->spectral_type = args->spectral_type;
htrdr->ref_temperature = args->ref_temperature;
+ htrdr->sky_mtl_name = args->sky_mtl_name;
res = init_mpi(htrdr);
if(res != RES_OK) goto error;
@@ -471,7 +537,7 @@ htrdr_init
/* Materials are necessary only if a ground geometry is defined */
if(args->filename_obj) {
- res = htrdr_mtl_create(htrdr, args->filename_mtl, &htrdr->mtl);
+ res = htrdr_materials_create(htrdr, args->filename_mtl, &htrdr->mats);
if(res != RES_OK) goto error;
}
@@ -479,11 +545,7 @@ htrdr_init
&htrdr->ground);
if(res != RES_OK) goto error;
- proj_ratio =
- (double)args->image.definition[0]
- / (double)args->image.definition[1];
- res = htrdr_camera_create(htrdr, args->camera.pos, args->camera.tgt,
- args->camera.up, proj_ratio, MDEG2RAD(args->camera.fov_y), &htrdr->cam);
+ res = setup_sensor(htrdr, args);
if(res != RES_OK) goto error;
res = htrdr_sun_create(htrdr, &htrdr->sun);
@@ -570,8 +632,10 @@ 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) {
- const size_t pixsz = htrdr_spectral_type_get_pixsz(htrdr->spectral_type);
- const size_t pixal = htrdr_spectral_type_get_pixal(htrdr->spectral_type);
+ const size_t pixsz = htrdr_spectral_type_get_pixsz
+ (htrdr->spectral_type, htrdr->sensor.type);
+ const size_t pixal = htrdr_spectral_type_get_pixal
+ (htrdr->spectral_type, htrdr->sensor.type);
res = htrdr_buffer_create(htrdr,
args->image.definition[0], /* Width */
@@ -599,9 +663,10 @@ htrdr_release(struct htrdr* htrdr)
if(htrdr->ground) htrdr_ground_ref_put(htrdr->ground);
if(htrdr->sky) HTSKY(ref_put(htrdr->sky));
if(htrdr->sun) htrdr_sun_ref_put(htrdr->sun);
- if(htrdr->cam) htrdr_camera_ref_put(htrdr->cam);
+ if(htrdr->sensor.camera) htrdr_camera_ref_put(htrdr->sensor.camera);
+ if(htrdr->sensor.rectangle) htrdr_rectangle_ref_put(htrdr->sensor.rectangle);
if(htrdr->buf) htrdr_buffer_ref_put(htrdr->buf);
- if(htrdr->mtl) htrdr_mtl_ref_put(htrdr->mtl);
+ if(htrdr->mats) htrdr_materials_ref_put(htrdr->mats);
if(htrdr->cie) htrdr_cie_xyz_ref_put(htrdr->cie);
if(htrdr->ran_wlen) htrdr_ran_wlen_ref_put(htrdr->ran_wlen);
if(htrdr->output && htrdr->output != stdout) fclose(htrdr->output);
@@ -639,14 +704,16 @@ htrdr_run(struct htrdr* htrdr)
}
}
} else {
- res = htrdr_draw_radiance(htrdr, htrdr->cam, htrdr->width,
- htrdr->height, htrdr->spp, htrdr->buf);
+ res = htrdr_draw_map(htrdr, &htrdr->sensor, htrdr->width, htrdr->height,
+ htrdr->spp, htrdr->buf);
if(res != RES_OK) goto error;
if(htrdr->mpi_rank == 0) {
struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL;
+ struct htrdr_accum flux_acc = HTRDR_ACCUM_NULL;
struct htrdr_estimate path_time;
+ struct htrdr_estimate flux;
- res = dump_buffer(htrdr, htrdr->buf, &path_time_acc,
+ res = dump_buffer(htrdr, htrdr->buf, &path_time_acc, &flux_acc,
str_cget(&htrdr->output_name), htrdr->output);
if(res != RES_OK) goto error;
@@ -655,6 +722,14 @@ htrdr_run(struct htrdr* htrdr)
"Time per radiative path (in micro seconds): %g +/- %g\n",
path_time.E,
path_time.SE);
+
+ if(htrdr->sensor.type == HTRDR_SENSOR_RECTANGLE) {
+ htrdr_accum_get_estimation(&flux_acc, &flux);
+ htrdr_log(htrdr,
+ "Radiative flux density (in W/(external m^2)): %g +/- %g\n",
+ flux.E,
+ flux.SE);
+ }
}
}
exit:
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -17,6 +17,7 @@
#ifndef HTRDR_H
#define HTRDR_H
+#include "htrdr_sensor.h"
#include "htrdr_spectral.h"
#include <rsys/logger.h>
@@ -37,6 +38,7 @@ struct htsky;
struct htrdr_args;
struct htrdr_buffer;
struct htrdr_cie_xyz;
+struct htrdr_materials;
struct htrdr_rectangle;
struct mem_allocator;
struct mutext;
@@ -49,15 +51,16 @@ struct htrdr {
struct s3d_device* s3d;
struct htrdr_ground* ground;
- struct htrdr_mtl* mtl;
+ struct htrdr_materials* mats;
struct htrdr_sun* sun;
struct htrdr_cie_xyz* cie;
struct htrdr_ran_wlen* ran_wlen;
- struct htrdr_camera* cam;
+ struct htrdr_sensor sensor;
struct htrdr_buffer* buf;
struct htsky* sky;
+ const char* sky_mtl_name;
enum htrdr_spectral_type spectral_type;
double wlen_range_m[2]; /* Integration range in *meters* */
double ref_temperature; /* Reference temperature in Kelvin */
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -34,8 +34,8 @@ print_help(const char* cmd)
ASSERT(cmd);
printf("Usage: %s [OPION]... -a ATMOSPHERE\n", cmd);
printf(
-"Render an image for scenes composed of an atmospheric gas mixture,\n"
-"clouds and a ground.\n\n");
+"Render an image or compute a flux map for scenes composed of an\n"
+"atmospheric gas mixture, clouds and a ground.\n\n");
printf(
" -a ATMOSPHERE gas optical properties of the atmosphere.\n");
printf(
@@ -65,12 +65,20 @@ print_help(const char* cmd)
printf(
" -m MIE file of Mie's data.\n");
printf(
+" -n SKY-NAME name used to identify the sky in the MATERIALS file.\n"
+" Its default value is `%s'.\n",
+ HTRDR_ARGS_DEFAULT.sky_mtl_name);
+ printf(
" -O CACHE name of the cache file used to store/restore the sky\n"
" data.\n");
printf(
" -o OUTPUT file where data are written. If not defined, data are\n"
" written to standard output.\n");
printf(
+" -p <rectangle> switch in flux computation by defining the rectangular\n"
+" sensor onto wich the flux is computed. Refer to the\n"
+" htrdr man page for the list of rectangle options.\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");
@@ -99,7 +107,7 @@ print_help(const char* cmd)
"Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier. htrdr is free\n"
"software released under the GNU GPL license, version 3 or later. You\n"
"are free to change or redistribute it under certain conditions\n"
-"<http://gnu.org/licenses/gpl.html>\n");
+"<http://gnu.org/licenses/gpl.html>.\n");
}
static INLINE res_T
@@ -261,6 +269,62 @@ error:
}
static res_T
+parse_rectangle_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 rectangle option string `%s'.\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ /* pos=0,0,10.1; key <- pos, val <- 0,0,10 */
+ key = strtok_r(buf, "=", &ctx);
+ val = strtok_r(NULL, "", &ctx);
+
+ if(!val) {
+ fprintf(stderr, "Missing value to the rectangle option `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ #define PARSE(Name, Func) { \
+ if(RES_OK != (res = Func)) { \
+ fprintf(stderr, "Invalid rectangle "Name" `%s'.\n", val); \
+ goto error; \
+ } \
+ } (void)0
+ if(!strcmp(key, "pos")) {
+ PARSE("position", parse_doubleX(val, args->rectangle.pos, 3));
+ } else if(!strcmp(key, "tgt")) {
+ PARSE("target", parse_doubleX(val, args->rectangle.tgt, 3));
+ } else if(!strcmp(key, "up")) {
+ PARSE("up vector", parse_doubleX(val, args->rectangle.up, 3));
+ } else if(!strcmp(key, "sz")) {
+ PARSE("size", parse_doubleX(val, args->rectangle.sz, 2));
+ } else {
+ fprintf(stderr, "Invalid rectangle parameter `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ #undef PARSE
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+
+
+static res_T
parse_spectral_range(const char* str, double wlen_range[2])
{
double range[2];
@@ -471,10 +535,11 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
}
}
- while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:M:m:O:o:Rrs:T:t:V:v")) != -1) {
+ while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:M:m:n:O:o:p:Rrs:T:t:V:v")) != -1) {
switch(opt) {
case 'a': args->filename_gas = optarg; break;
case 'C':
+ args->sensor_type = HTRDR_SENSOR_CAMERA;
res = parse_multiple_parameters
(args, optarg, parse_camera_parameter);
break;
@@ -494,8 +559,14 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
break;
case 'M': args->filename_mtl = optarg; break;
case 'm': args->filename_mie = optarg; break;
+ case 'n': args->sky_mtl_name = optarg; break;
case 'O': args->cache = optarg; break;
case 'o': args->output = optarg; break;
+ case 'p':
+ args->sensor_type = HTRDR_SENSOR_RECTANGLE;
+ res = parse_multiple_parameters
+ (args, optarg, parse_rectangle_parameter);
+ break;
case 'r': args->repeat_clouds = 1; break;
case 'R': args->repeat_ground = 1; break;
case 's':
@@ -544,7 +615,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
/* Setup default ref temperature if necessary */
if(args->ref_temperature <= 0) {
switch(args->spectral_type) {
- case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_LW:
args->ref_temperature = HTRDR_DEFAULT_LW_REF_TEMPERATURE;
break;
case HTRDR_SPECTRAL_SW:
diff --git a/src/htrdr_args.h.in b/src/htrdr_args.h.in
@@ -16,6 +16,7 @@
#ifndef HTRDR_ARGS_H
#define HTRDR_ARGS_H
+#include "htrdr_sensor.h"
#include "htrdr_spectral.h"
#include "htrdr_cie_xyz.h"
@@ -31,9 +32,10 @@ struct htrdr_args {
const char* filename_mtl; /* Path of the materials */
const char* cache;
const char* output;
+ const char* sky_mtl_name;
struct {
- double pos[3]; /* Position */
+ double pos[3]; /* Center of the renctangle */
double tgt[3]; /* Target */
double up[3]; /* Up vector */
double sz[2]; /* Plane size in world space */
@@ -60,6 +62,7 @@ struct htrdr_args {
double wlen_range[2]; /* Spectral range of integration in nm */
double ref_temperature; /* Planck reference temperature in Kelvin */
+ enum htrdr_sensor_type sensor_type;
unsigned nthreads; /* Hint on the number of threads to use */
int force_overwriting;
int dump_vtk; /* Dump the loaded cloud properties in a VTK file */
@@ -77,11 +80,12 @@ struct htrdr_args {
NULL, /* Mtl filename */ \
NULL, /* Cache filename */ \
NULL, /* Output filename */ \
+ @HTRDR_ARGS_DEFAULT_SKY_MTL_NAME@, /* Sky mtl name */ \
{ \
- {0, 0, 0}, /* plane position */ \
- {0, 0, 1}, /* plane target */ \
- {0, 1, 0}, /* plane up */ \
- {1, 1}, /* plane size */ \
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_POS@}, /* Rectangle center */ \
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_TGT@}, /* Rectangle target */ \
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_UP@}, /* Rectangle up */ \
+ {@HTRDR_ARGS_DEFAULT_RECTANGLE_SZ@}, /* Rectangle size */ \
}, { \
{@HTRDR_ARGS_DEFAULT_CAMERA_POS@}, /* Camera position */ \
{@HTRDR_ARGS_DEFAULT_CAMERA_TGT@}, /* Camera target */ \
@@ -98,6 +102,7 @@ struct htrdr_args {
HTRDR_SPECTRAL_SW_CIE_XYZ, /* Spectral type */ \
HTRDR_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
-1, /* Reference temperature */ \
+ HTRDR_SENSOR_CAMERA, /* sensor type */ \
(unsigned)~0, /* #threads */ \
0, /* Force overwriting */ \
0, /* dump VTK */ \
diff --git a/src/htrdr_compute_radiance_lw.c b/src/htrdr_compute_radiance_lw.c
@@ -228,6 +228,7 @@ htrdr_compute_radiance_lw
/* Scattering at a surface */
} else {
struct htrdr_interface interf = HTRDR_INTERFACE_NULL;
+ const struct htrdr_mtl* mtl = NULL;
struct ssf_bsdf* bsdf = NULL;
double bounce_reflectivity = 0;
double N[3];
@@ -235,10 +236,10 @@ htrdr_compute_radiance_lw
ASSERT(ctx.event_type == EVENT_NONE);
ASSERT(!S3D_HIT_NONE(&s3d_hit));
- /* Fetch the hit interface and build its BSDF */
+ /* Fetch the hit interface materal and build its BSDF */
htrdr_ground_get_interface(htrdr->ground, &s3d_hit, &interf);
- HTRDR(interface_create_bsdf
- (htrdr, &interf, ithread, wlen, pos_next, dir, rng, &s3d_hit, &bsdf));
+ mtl = htrdr_interface_fetch_hit_mtl(&interf, dir, &s3d_hit);
+ HTRDR(mtl_create_bsdf(htrdr, mtl, ithread, wlen, rng, &bsdf));
d3_normalize(N, d3_set_f3(N, s3d_hit.normal));
if(d3_dot(N, wo) < 0) d3_minus(N, N);
@@ -253,19 +254,11 @@ htrdr_compute_radiance_lw
SSF(bsdf_ref_put(bsdf));
if(ssp_rng_canonical(rng) >= bounce_reflectivity) { /* Absorbed at boundary */
- /* Retrieve the temperature of the interface. Anyway, since we do not
- * have this data yet, we use the temperature of the sky at the current
- * position as the temperature of the surface. */
- temperature = htsky_fetch_temperature(htrdr->sky, pos_next);
- if(temperature <= 0) {
- w = 0;
- } else {
- /* weight is planck integrated over the spectral sub-interval */
- w = planck_monochromatic(wlen_m, temperature);
- }
+ temperature = mtl->temperature; /* Fetch mtl temperature */
+ /* Weight is planck integrated over the spectral sub-interval */
+ w = temperature > 0 ? planck_monochromatic(wlen_m, temperature) : 0;
break;
}
-
s3d_hit_prev = s3d_hit;
}
d3_set(pos, pos_next);
diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c
@@ -252,6 +252,7 @@ htrdr_compute_radiance_sw
(struct htrdr* htrdr,
const size_t ithread,
struct ssp_rng* rng,
+ const int cpnt_mask, /* Combination of enum htrdr_radiance_cpnt_flag */
const double pos_in[3],
const double dir_in[3],
const double wlen, /* In nanometer */
@@ -276,7 +277,6 @@ htrdr_compute_radiance_sw
double r; /* Random number */
double wo[3]; /* -dir */
double pdf;
- double sun_solid_angle;
double Tr; /* Overall transmissivity */
double Tr_abs; /* Absorption transmissivity */
double L_sun; /* Sun radiance in W.m^-2.sr^-1 */
@@ -304,17 +304,14 @@ htrdr_compute_radiance_sw
* default "ecrad_opt_prot.txt" file provided by the HTGOP project. */
htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
ASSERT(band_bounds[0] <= wlen && wlen <= band_bounds[1]);
- sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
d3_set(pos, pos_in);
d3_set(dir, dir_in);
- /* Add the directly contribution of the sun */
- if(htrdr_sun_is_dir_in_solar_cone(htrdr->sun, dir)) {
- /* Add the direct contribution of the sun */
- d2(range, 0, FLT_MAX);
-
+ if((cpnt_mask & HTRDR_RADIANCE_DIRECT) /* Handle direct contribuation */
+ && htrdr_sun_is_dir_in_solar_cone(htrdr->sun, dir)) {
/* Check that the ray is not occlude along the submitted range */
+ d2(range, 0, FLT_MAX);
HTRDR(ground_trace_ray(htrdr->ground, pos, dir, range, NULL, &s3d_hit_tmp));
if(!S3D_HIT_NONE(&s3d_hit_tmp)) {
Tr = 0;
@@ -325,10 +322,19 @@ htrdr_compute_radiance_sw
}
}
+ if((cpnt_mask & HTRDR_RADIANCE_DIFFUSE) == 0)
+ goto exit; /* Discard diffuse contribution */
+
/* Radiative random walk */
for(;;) {
struct scattering_context scattering_ctx = SCATTERING_CONTEXT_NULL;
+ struct ssf_bsdf* bsdf = NULL;
+ struct ssf_phase* phase;
+ double N[3];
double bounce_reflectivity = 1;
+ double sun_dir_pdf;
+ int surface_scattering = 0; /* Define if hit a surface */
+ int bsdf_type = 0;
/* Find the first intersection with a surface */
d2(range, 0, DBL_MAX);
@@ -372,57 +378,30 @@ htrdr_compute_radiance_sw
(htrdr, rng, HTSKY_Ka, iband, iquad, pos, dir, range);
if(Tr_abs <= 0) break;
- /* Sample a sun direction */
- htrdr_sun_sample_direction(htrdr->sun, rng, sun_dir);
- d2(range, 0, FLT_MAX);
-
- s3d_hit_prev = SVX_HIT_NONE(&svx_hit) ? s3d_hit : S3D_HIT_NULL;
-
- /* Check that the sun is visible from the new position */
- HTRDR(ground_trace_ray
- (htrdr->ground, pos_next, sun_dir, range, &s3d_hit_prev, &s3d_hit_tmp));
+ surface_scattering = SVX_HIT_NONE(&svx_hit);
- /* Compute the sun transmissivity */
- if(!S3D_HIT_NONE(&s3d_hit_tmp)) {
- Tr = 0;
- } else {
- Tr = transmissivity
- (htrdr, rng, HTSKY_Kext, iband, iquad, pos_next, sun_dir, range);
- }
-
- /* Scattering at a surface */
- if(SVX_HIT_NONE(&svx_hit)) {
+ /* Sample the scattering direction */
+ if(surface_scattering) { /* Scattering at a surface */
struct htrdr_interface interf = HTRDR_INTERFACE_NULL;
- struct ssf_bsdf* bsdf = NULL;
- double N[3];
- int type;
+ const struct htrdr_mtl* mtl = NULL;
- /* Fetch the hit interface and build its BSDF */
+ /* Fetch the hit interface materal and build its BSDF */
htrdr_ground_get_interface(htrdr->ground, &s3d_hit, &interf);
- HTRDR(interface_create_bsdf
- (htrdr, &interf, ithread, wlen, pos_next, dir, rng, &s3d_hit, &bsdf));
+ mtl = htrdr_interface_fetch_hit_mtl(&interf, dir, &s3d_hit);
+ HTRDR(mtl_create_bsdf(htrdr, mtl, ithread, wlen, rng, &bsdf));
+ /* Revert the normal if necessary to match the SSF convention */
d3_normalize(N, d3_set_f3(N, s3d_hit.normal));
if(d3_dot(N, wo) < 0) d3_minus(N, N);
+ /* Sample scattering direction */
bounce_reflectivity = ssf_bsdf_sample
- (bsdf, rng, wo, N, dir_next, &type, &pdf);
- if(!(type & SSF_REFLECTION)) { /* Handle only reflections */
+ (bsdf, rng, wo, N, dir_next, &bsdf_type, &pdf);
+ if(!(bsdf_type & SSF_REFLECTION)) { /* Handle only reflections */
bounce_reflectivity = 0;
}
- if(d3_dot(N, sun_dir) < 0) { /* Below the ground */
- R = 0;
- } else {
- R = ssf_bsdf_eval(bsdf, wo, N, sun_dir) * d3_dot(N, sun_dir);
- }
-
- /* Release the BSDF */
- SSF(bsdf_ref_put(bsdf));
-
- /* Scattering in the medium */
- } else {
- struct ssf_phase* phase;
+ } else { /* Scattering in a volume */
double ks_particle; /* Scattering coefficient of the particles */
double ks_gas; /* Scattering coefficient of the gaz */
double ks; /* Overall scattering coefficient */
@@ -440,22 +419,74 @@ htrdr_compute_radiance_sw
phase = phase_hg;
}
+ /* Sample scattering direction */
ssf_phase_sample(phase, rng, wo, dir_next, NULL);
- R = ssf_phase_eval(phase, wo, sun_dir);
+ ssf_phase_ref_get(phase);
+ }
+
+ /* Sample the direction of the direct contribution */
+ if(surface_scattering && (bsdf_type & SSF_SPECULAR)) {
+ if(!htrdr_sun_is_dir_in_solar_cone(htrdr->sun, dir_next)) {
+ R = 0; /* No direct lightning */
+ } else {
+ sun_dir[0] = dir_next[0];
+ sun_dir[1] = dir_next[1];
+ sun_dir[2] = dir_next[2];
+ R = d3_dot(N, sun_dir)<0/* Below the ground*/ ? 0 : bounce_reflectivity;
+ }
+ sun_dir_pdf = 1.0;
+ } else {
+ /* Sample a sun direction */
+ sun_dir_pdf = htrdr_sun_sample_direction(htrdr->sun, rng, sun_dir);
+ if(surface_scattering) {
+ R = d3_dot(N, sun_dir) < 0/* Below the ground */
+ ? 0 : ssf_bsdf_eval(bsdf, wo, N, sun_dir) * d3_dot(N, sun_dir);
+ } else {
+ R = ssf_phase_eval(phase, wo, sun_dir);
+ }
+ }
+
+ /* The direct contribution to the scattering point is not null so we need
+ * to compute the transmissivity from sun to scatt pt */
+ if(R <= 0) {
+ Tr = 0;
+ } else {
+ /* Check that the sun is visible from the new position */
+ d2(range, 0, FLT_MAX);
+ s3d_hit_prev = SVX_HIT_NONE(&svx_hit) ? s3d_hit : S3D_HIT_NULL;
+ HTRDR(ground_trace_ray
+ (htrdr->ground, pos_next, sun_dir, range, &s3d_hit_prev, &s3d_hit_tmp));
+
+ /* Compute the sun transmissivity */
+ if(!S3D_HIT_NONE(&s3d_hit_tmp)) {
+ Tr = 0;
+ } else {
+ Tr = transmissivity
+ (htrdr, rng, HTSKY_Kext, iband, iquad, pos_next, sun_dir, range);
+ }
+ }
+
+ /* Release the scattering function */
+ if(surface_scattering) {
+ SSF(bsdf_ref_put(bsdf));
+ } else {
+ SSF(phase_ref_put(phase));
}
/* Update the MC weight */
ksi *= Tr_abs;
- w += ksi * L_sun * sun_solid_angle * Tr * R;
+ w += ksi * L_sun * Tr * R / sun_dir_pdf;
/* Russian roulette wrt surface scattering */
- if(SVX_HIT_NONE(&svx_hit) && ssp_rng_canonical(rng) >= bounce_reflectivity)
+ if(surface_scattering && ssp_rng_canonical(rng) >= bounce_reflectivity)
break;
/* Setup the next random walk state */
d3_set(pos, pos_next);
d3_set(dir, dir_next);
}
+
+exit:
SSF(phase_ref_put(phase_hg));
SSF(phase_ref_put(phase_rayleigh));
return w;
diff --git a/src/htrdr_draw_map.c b/src/htrdr_draw_map.c
@@ -0,0 +1,1207 @@
+/* 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 /* nanosleep && nextafter */
+
+#include "htrdr.h"
+#include "htrdr_c.h"
+#include "htrdr_buffer.h"
+#include "htrdr_camera.h"
+#include "htrdr_cie_xyz.h"
+#include "htrdr_ran_wlen.h"
+#include "htrdr_rectangle.h"
+#include "htrdr_solve.h"
+#include "htrdr_sun.h"
+
+#include <high_tune/htsky.h>
+
+#include <rsys/algorithm.h>
+#include <rsys/clock_time.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_u32.h>
+#include <rsys/math.h>
+#include <rsys/mutex.h>
+#include <star/ssp.h>
+
+#include <omp.h>
+#include <mpi.h>
+#include <time.h>
+#include <unistd.h>
+
+#define RNG_SEQUENCE_SIZE 10000
+
+#define TILE_MCODE_NULL UINT32_MAX
+#define TILE_SIZE 32 /* Definition in X & Y of a tile */
+STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
+
+enum pixel_format {
+ PIXEL_XWAVE,
+ PIXEL_IMAGE
+};
+
+union pixel {
+ struct htrdr_pixel_flux flux;
+ struct htrdr_pixel_xwave xwave;
+ struct htrdr_pixel_image image;
+};
+
+/* Tile of row ordered image pixels */
+struct tile {
+ struct list_node node;
+ struct mem_allocator* allocator;
+ ref_T ref;
+
+ struct tile_data {
+ uint16_t x, y; /* 2D coordinates of the tile in tile space */
+ enum pixel_format format;
+ /* Simulate the flexible array member of the C99 standard. */
+ union pixel pixels[1/*Dummy element*/];
+ } data;
+};
+
+/* List of tile to compute onto the MPI process. */
+struct proc_work {
+ struct mutex* mutex;
+ struct darray_u32 tiles; /* #tiles to render */
+ size_t itile; /* Next tile to render in the above list of tiles */
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static FINLINE uint16_t
+morton2D_decode(const uint32_t u32)
+{
+ uint32_t x = u32 & 0x55555555;
+ x = (x | (x >> 1)) & 0x33333333;
+ x = (x | (x >> 2)) & 0x0F0F0F0F;
+ x = (x | (x >> 4)) & 0x00FF00FF;
+ x = (x | (x >> 8)) & 0x0000FFFF;
+ return (uint16_t)x;
+}
+
+static FINLINE uint32_t
+morton2D_encode(const uint16_t u16)
+{
+ uint32_t u32 = u16;
+ u32 = (u32 | (u32 << 8)) & 0x00FF00FF;
+ u32 = (u32 | (u32 << 4)) & 0X0F0F0F0F;
+ u32 = (u32 | (u32 << 2)) & 0x33333333;
+ u32 = (u32 | (u32 << 1)) & 0x55555555;
+ 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)
+{
+ struct tile* tile;
+ const size_t tile_sz =
+ sizeof(struct tile) - sizeof(union pixel)/*rm dummy pixel*/;
+ const size_t buf_sz = /* Flexiblbe array element */
+ TILE_SIZE*TILE_SIZE*sizeof(union pixel);
+ ASSERT(allocator);
+
+ tile = MEM_ALLOC(allocator, tile_sz+buf_sz);
+ if(!tile) return NULL;
+
+ tile->data.format = fmt;
+
+ ref_init(&tile->ref);
+ list_init(&tile->node);
+ tile->allocator = allocator;
+ ASSERT(IS_ALIGNED(&tile->data.pixels, ALIGNOF(union pixel)));
+
+ return tile;
+}
+
+static INLINE void
+tile_ref_get(struct tile* tile)
+{
+ ASSERT(tile);
+ tile_ref_get(tile);
+}
+
+static INLINE void
+release_tile(ref_T* ref)
+{
+ struct tile* tile = CONTAINER_OF(ref, struct tile, ref);
+ ASSERT(ref);
+ MEM_RM(tile->allocator, tile);
+}
+
+static INLINE void
+tile_ref_put(struct tile* tile)
+{
+ ASSERT(tile);
+ ref_put(&tile->ref, release_tile);
+}
+
+static FINLINE union pixel*
+tile_at
+ (struct tile* tile,
+ const size_t x, /* In tile space */
+ const size_t y) /* In tile space */
+{
+ ASSERT(tile && x < TILE_SIZE && y < TILE_SIZE);
+ return tile->data.pixels + (y*TILE_SIZE + x);
+}
+
+static void
+write_tile_data
+ (struct htrdr* htrdr,
+ const struct htrdr_sensor* sensor,
+ struct htrdr_buffer* buf,
+ const struct tile_data* tile_data)
+{
+ struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL;
+ size_t icol, irow;
+ size_t irow_tile;
+ size_t ncols_tile, nrows_tile;
+ char* buf_mem;
+ ASSERT(htrdr && sensor && buf && tile_data);
+ (void)htrdr, (void)sensor;
+
+ htrdr_buffer_get_layout(buf, &layout);
+ buf_mem = htrdr_buffer_get_data(buf);
+ ASSERT(layout.elmt_size
+ == htrdr_spectral_type_get_pixsz(htrdr->spectral_type, sensor->type));
+
+ /* Compute the row/column of the tile origin into the buffer */
+ icol = tile_data->x * (size_t)TILE_SIZE;
+ irow = tile_data->y * (size_t)TILE_SIZE;
+
+ /* Define the number of tile row/columns to write into the buffer */
+ ncols_tile = MMIN(icol + TILE_SIZE, layout.width) - icol;
+ nrows_tile = MMIN(irow + TILE_SIZE, layout.height) - irow;
+
+ /* Copy the row ordered tile data */
+ FOR_EACH(irow_tile, 0, nrows_tile) {
+ char* buf_row = buf_mem + (irow + irow_tile) * layout.pitch;
+ char* buf_col = buf_row + icol * layout.elmt_size;
+ const union pixel* tile_row = tile_data->pixels + irow_tile*TILE_SIZE;
+ size_t x;
+
+ FOR_EACH(x, 0, ncols_tile) {
+ switch(tile_data->format) {
+ case PIXEL_XWAVE:
+ ((struct htrdr_pixel_xwave*)buf_col)[x] = tile_row[x].xwave;
+ break;
+ case PIXEL_IMAGE:
+ ((struct htrdr_pixel_image*)buf_col)[x] = tile_row[x].image;
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+ }
+}
+
+static INLINE void
+proc_work_init(struct mem_allocator* allocator, struct proc_work* work)
+{
+ ASSERT(work);
+ darray_u32_init(allocator, &work->tiles);
+ work->itile = 0;
+ CHK(work->mutex = mutex_create());
+}
+
+static INLINE void
+proc_work_release(struct proc_work* work)
+{
+ darray_u32_release(&work->tiles);
+ mutex_destroy(work->mutex);
+}
+
+static INLINE void
+proc_work_reset(struct proc_work* work)
+{
+ ASSERT(work);
+ mutex_lock(work->mutex);
+ darray_u32_clear(&work->tiles);
+ work->itile = 0;
+ mutex_unlock(work->mutex);
+}
+
+static INLINE void
+proc_work_add_tile(struct proc_work* work, const uint32_t mcode)
+{
+ mutex_lock(work->mutex);
+ CHK(darray_u32_push_back(&work->tiles, &mcode) == RES_OK);
+ mutex_unlock(work->mutex);
+}
+
+static INLINE uint32_t
+proc_work_get_tile(struct proc_work* work)
+{
+ uint32_t mcode;
+ ASSERT(work);
+ mutex_lock(work->mutex);
+ if(work->itile >= darray_u32_size_get(&work->tiles)) {
+ mcode = TILE_MCODE_NULL;
+ } else {
+ mcode = darray_u32_cdata_get(&work->tiles)[work->itile];
+ ++work->itile;
+ }
+ mutex_unlock(work->mutex);
+ return mcode;
+}
+
+static INLINE size_t
+proc_work_get_ntiles(struct proc_work* work)
+{
+ size_t sz = 0;
+ ASSERT(work);
+ mutex_lock(work->mutex);
+ sz = darray_u32_size_get(&work->tiles);
+ mutex_unlock(work->mutex);
+ return sz;
+}
+
+static void
+mpi_wait_for_request(struct htrdr* htrdr, MPI_Request* req)
+{
+ ASSERT(htrdr && req);
+
+ /* Wait for process synchronisation */
+ for(;;) {
+ struct timespec t;
+ int complete;
+ t.tv_sec = 0;
+ t.tv_nsec = 10000000; /* 10ms */
+
+ mutex_lock(htrdr->mpi_mutex);
+ MPI(Test(req, &complete, MPI_STATUS_IGNORE));
+ mutex_unlock(htrdr->mpi_mutex);
+ if(complete) break;
+
+ nanosleep(&t, NULL);
+ }
+}
+
+static void
+mpi_probe_thieves
+ (struct htrdr* htrdr,
+ struct proc_work* work,
+ ATOMIC* probe_thieves)
+{
+ uint32_t tiles[UINT8_MAX];
+ struct timespec t;
+ ASSERT(htrdr && work && probe_thieves);
+
+ if(htrdr->mpi_nprocs == 1) /* The process is alone. No thief is possible */
+ return;
+
+ t.tv_sec = 0;
+
+ /* Protect MPI calls of multiple invocations from concurrent threads */
+ #define P_MPI(Func) { \
+ mutex_lock(htrdr->mpi_mutex); \
+ MPI(Func); \
+ mutex_unlock(htrdr->mpi_mutex); \
+ } (void)0
+
+ while(ATOMIC_GET(probe_thieves)) {
+ MPI_Status status;
+ size_t itile;
+ int msg;
+
+ /* Probe if a steal request was submitted by any processes */
+ P_MPI(Iprobe(MPI_ANY_SOURCE, HTRDR_MPI_STEAL_REQUEST, MPI_COMM_WORLD, &msg,
+ &status));
+
+ if(msg) { /* A steal request was posted */
+ MPI_Request req;
+ uint8_t ntiles_to_steal;
+
+ /* Asynchronously receive the steal request */
+ P_MPI(Irecv(&ntiles_to_steal, 1, MPI_UINT8_T, status.MPI_SOURCE,
+ HTRDR_MPI_STEAL_REQUEST, MPI_COMM_WORLD, &req));
+
+ /* Wait for the completion of the steal request */
+ mpi_wait_for_request(htrdr, &req);
+
+ /* Thief some tiles */
+ FOR_EACH(itile, 0, ntiles_to_steal) {
+ tiles[itile] = proc_work_get_tile(work);
+ }
+ P_MPI(Send(&tiles, ntiles_to_steal, MPI_UINT32_T, status.MPI_SOURCE,
+ HTRDR_MPI_WORK_STEALING, MPI_COMM_WORLD));
+ }
+ t.tv_nsec = 500000000; /* 500ms */
+ nanosleep(&t, NULL);
+ }
+ #undef P_MPI
+}
+
+static int
+mpi_sample_working_process(struct htrdr* htrdr, struct ssp_rng* rng)
+{
+ int iproc, i;
+ int dst_rank;
+ ASSERT(htrdr && rng && htrdr->mpi_nworking_procs);
+
+ /* Sample the index of the 1st active process */
+ iproc = (int)(ssp_rng_canonical(rng) * (double)htrdr->mpi_nworking_procs);
+
+ /* Find the rank of the sampled active process. Use a simple linear search
+ * since the overall number of processes should be quite low; at most few
+ * dozens. */
+ i = 0;
+ FOR_EACH(dst_rank, 0, htrdr->mpi_nprocs) {
+ if(htrdr->mpi_working_procs[dst_rank] == 0) continue; /* Inactive process */
+ if(i == iproc) break; /* The rank of the sampled process is found */
+ ++i;
+ }
+ ASSERT(dst_rank < htrdr->mpi_nprocs);
+ return dst_rank;
+}
+
+/* Return the number of stolen tiles */
+static size_t
+mpi_steal_work
+ (struct htrdr* htrdr,
+ struct ssp_rng* rng,
+ struct proc_work* work)
+{
+ MPI_Request req;
+ size_t itile;
+ size_t nthieves = 0;
+ uint32_t tiles[UINT8_MAX]; /* Morton code of the stolen tile */
+ int proc_to_steal; /* Process to steal */
+ uint8_t ntiles_to_steal = MMIN((uint8_t)(htrdr->nthreads*2), 16);
+ ASSERT(htrdr && rng && work && htrdr->nthreads < UINT8_MAX);
+
+ /* Protect MPI calls of multiple invocations from concurrent threads */
+ #define P_MPI(Func) { \
+ mutex_lock(htrdr->mpi_mutex); \
+ MPI(Func); \
+ mutex_unlock(htrdr->mpi_mutex); \
+ } (void)0
+
+ /* No more working process => nohting to steal */
+ if(!htrdr->mpi_nworking_procs) return 0;
+
+ /* Sample a process to steal */
+ proc_to_steal = mpi_sample_working_process(htrdr, rng);
+
+ /* Send a steal request to the sampled process and wait for a response */
+ P_MPI(Send(&ntiles_to_steal, 1, MPI_UINT8_T, proc_to_steal,
+ HTRDR_MPI_STEAL_REQUEST, MPI_COMM_WORLD));
+
+ /* Receive the stolen tile from the sampled process */
+ P_MPI(Irecv(tiles, ntiles_to_steal, MPI_UINT32_T, proc_to_steal,
+ HTRDR_MPI_WORK_STEALING, MPI_COMM_WORLD, &req));
+
+ mpi_wait_for_request(htrdr, &req);
+
+ FOR_EACH(itile, 0, ntiles_to_steal) {
+ if(tiles[itile] == TILE_MCODE_NULL) {
+ ASSERT(htrdr->mpi_working_procs[proc_to_steal] != 0);
+ htrdr->mpi_working_procs[proc_to_steal] = 0;
+ htrdr->mpi_nworking_procs--;
+ break;
+ }
+ proc_work_add_tile(work, tiles[itile]);
+ ++nthieves;
+ }
+ #undef P_MPI
+ return nthieves;
+}
+
+static res_T
+mpi_gather_tiles
+ (struct htrdr* htrdr,
+ const struct htrdr_sensor* sensor,
+ struct htrdr_buffer* buf,
+ const size_t ntiles,
+ struct list_node* tiles)
+{
+ /* Compute the size of the tile_data */
+ const size_t msg_sz =
+ sizeof(struct tile_data) - sizeof(union pixel)/*dummy*/
+ + TILE_SIZE*TILE_SIZE*sizeof(union pixel);
+
+ struct list_node* node = NULL;
+ struct tile* tile = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && tiles);
+ ASSERT(htrdr->mpi_rank != 0 || buf);
+ (void)ntiles;
+
+ if(htrdr->mpi_rank != 0) { /* Non master process */
+ /* Send the computed tile to the master process */
+ LIST_FOR_EACH(node, tiles) {
+ struct tile* t = CONTAINER_OF(node, struct tile, node);
+ MPI(Send(&t->data, (int)msg_sz, MPI_CHAR, 0,
+ HTRDR_MPI_TILE_DATA, MPI_COMM_WORLD));
+ }
+ } else { /* Master process */
+ size_t itile = 0;
+
+ LIST_FOR_EACH(node, tiles) {
+ struct tile* t = CONTAINER_OF(node, struct tile, node);
+ write_tile_data(htrdr, sensor, 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 */
+ pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
+ tile = tile_create(htrdr->allocator, pixfmt);
+ if(!tile) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "could not allocate the temporary tile used to gather the process "
+ "output data -- %s.\n", res_to_cstr(res));
+ goto error;
+ }
+
+ /* Receive the tile data of the concurrent MPI processes */
+ 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(htrdr, sensor, buf, &tile->data);
+ }
+ }
+ }
+
+exit:
+ if(tile) tile_ref_put(tile);
+ return res;
+error:
+ goto exit;
+}
+
+static void
+draw_pixel_image
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ const size_t ipix[2],
+ const double pix_sz[2], /* Size of a pixel in the normalized image plane */
+ const struct htrdr_camera* cam,
+ const size_t spp,
+ struct ssp_rng* rng,
+ 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(htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ);
+
+ /* Reset accumulators */
+ XYZ[0] = HTRDR_ACCUM_NULL;
+ XYZ[1] = HTRDR_ACCUM_NULL;
+ XYZ[2] = HTRDR_ACCUM_NULL;
+ time = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(ichannel, 0, 3) {
+ size_t isamp;
+
+ FOR_EACH(isamp, 0, spp) {
+ struct time t0, t1;
+ double pix_samp[2];
+ double ray_org[3];
+ double ray_dir[3];
+ 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;
+
+ /* Begin the registration of the time spent to in the realisation */
+ time_current(&t0);
+
+ /* Sample a position into the pixel, in the normalized image plane */
+ pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
+ pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
+
+ /* Generate a ray starting from the pinhole camera and passing through the
+ * pixel sample */
+ htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir);
+
+ r0 = ssp_rng_canonical(rng);
+ r1 = ssp_rng_canonical(rng);
+ r2 = ssp_rng_canonical(rng);
+
+ /* Sample a spectral band and a quadrature point */
+ switch(ichannel) {
+ case 0: wlen = htrdr_cie_xyz_sample_X(htrdr->cie, r0, 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 radiance in W/m^2/sr/m */
+ weight = htrdr_compute_radiance_sw(htrdr, ithread, rng,
+ HTRDR_RADIANCE_ALL, 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;
+
+ /* Update the pixel accumulator of the current channel */
+ XYZ[ichannel].sum_weights += weight;
+ XYZ[ichannel].sum_weights_sqr += weight*weight;
+ XYZ[ichannel].nweights += 1;
+
+ /* Update the pixel accumulator of per realisation time */
+ time.sum_weights += usec;
+ time.sum_weights_sqr += usec*usec;
+ time.nweights += 1;
+ }
+ }
+
+ /* Flush pixel data */
+ htrdr_accum_get_estimation(XYZ+0, &pixel->X);
+ htrdr_accum_get_estimation(XYZ+1, &pixel->Y);
+ htrdr_accum_get_estimation(XYZ+2, &pixel->Z);
+ pixel->time = time;
+}
+
+static void
+draw_pixel_flux
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ const size_t ipix[2],
+ const double pix_sz[2], /* Size of a pixel in the normalized image plane */
+ const struct htrdr_sensor* sensor,
+ const size_t spp,
+ struct ssp_rng* rng,
+ struct htrdr_pixel_flux* pixel)
+{
+ struct htrdr_accum flux;
+ struct htrdr_accum time;
+ size_t isamp;
+ ASSERT(ipix && ipix && pix_sz && sensor && rng && pixel);
+ ASSERT(sensor->type == HTRDR_SENSOR_RECTANGLE);
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_LW
+ || htrdr->spectral_type == HTRDR_SPECTRAL_SW);
+
+ /* Reset the pixel accumulators */
+ flux = HTRDR_ACCUM_NULL;
+ time = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(isamp, 0, spp) {
+ struct time t0, t1;
+ double ray_org[3];
+ double ray_dir[3];
+ double weight;
+ double r0, r1, r2;
+ double wlen;
+ size_t iband;
+ size_t iquad;
+ double usec;
+ double band_pdf;
+ res_T res = RES_OK;
+
+ /* Begin the registration of the time spent in the realisation */
+ time_current(&t0);
+
+ res = htrdr_sensor_sample_primary_ray(&htrdr->sensor, htrdr, ipix,
+ pix_sz, rng, ray_org, ray_dir);
+ if(res != RES_OK) continue; /* Reject the current sample */
+
+ r0 = ssp_rng_canonical(rng);
+ r1 = ssp_rng_canonical(rng);
+ r2 = ssp_rng_canonical(rng);
+
+ /* Sample a wavelength */
+ 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, r2, iband);
+
+ if(htrdr->spectral_type == HTRDR_SPECTRAL_LW) {
+ weight = htrdr_compute_radiance_lw(htrdr, ithread, rng, ray_org,
+ ray_dir, wlen, iband, iquad);
+ weight *= PI / band_pdf; /* Transform weight from W/m^2/sr/m to W/m^2 */
+ } else {
+ double sun_dir[3];
+ double N[3];
+ double L_direct;
+ double L_diffuse;
+ double cos_N_sun_dir;
+ double sun_solid_angle;
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_SW);
+
+ /* Compute direct contribution if necessary */
+ htrdr_sun_sample_direction(htrdr->sun, rng, sun_dir);
+ htrdr_rectangle_get_normal(sensor->rectangle, N);
+ cos_N_sun_dir = d3_dot(N, sun_dir);
+
+ if(cos_N_sun_dir <= 0) {
+ L_direct = 0;
+ } else {
+ L_direct = htrdr_compute_radiance_sw(htrdr, ithread, rng,
+ HTRDR_RADIANCE_DIRECT, ray_org, sun_dir, wlen, iband, iquad);
+ }
+
+ /* Compute diffuse contribution */
+ L_diffuse = htrdr_compute_radiance_sw(htrdr, ithread, rng,
+ HTRDR_RADIANCE_DIFFUSE, ray_org, ray_dir, wlen, iband, iquad);
+
+ sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
+
+ /* Compute the weight in W/m^2/m */
+ weight = cos_N_sun_dir * sun_solid_angle * L_direct + PI * L_diffuse;
+
+ /* Importance sampling: correct weight with pdf */
+ weight /= band_pdf; /* In W/m^2 */
+ }
+
+ /* End the registration of the per realisation time */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update the pixel accumulator of the flux */
+ flux.sum_weights += weight;
+ flux.sum_weights_sqr += weight*weight;
+ flux.nweights += 1;
+
+ /* Update the pixel accumulator of per realisation time */
+ time.sum_weights += usec;
+ time.sum_weights_sqr += usec*usec;
+ time.nweights += 1;
+ }
+
+ /* Save the per realisation integration time */
+ pixel->flux = flux;
+ 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_xwave
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ const size_t ipix[2],
+ const double pix_sz[2], /* Size of a pixel in the normalized image plane */
+ const struct htrdr_sensor* sensor,
+ const size_t spp,
+ struct ssp_rng* rng,
+ 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 && sensor && rng && pixel);
+ ASSERT(sensor->type == HTRDR_SENSOR_CAMERA);
+ ASSERT(htrdr->spectral_type == HTRDR_SPECTRAL_LW
+ || htrdr->spectral_type == HTRDR_SPECTRAL_SW);
+
+ /* Reset the pixel accumulators */
+ radiance = HTRDR_ACCUM_NULL;
+ time = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(isamp, 0, spp) {
+ struct time t0, t1;
+ double ray_org[3];
+ double ray_dir[3];
+ double weight;
+ double r0, r1, r2;
+ double wlen;
+ size_t iband;
+ size_t iquad;
+ double usec;
+ double band_pdf;
+ res_T res = RES_OK;
+
+ /* Begin the registration of the time spent in the realisation */
+ time_current(&t0);
+
+ res = htrdr_sensor_sample_primary_ray(sensor, htrdr, ipix,
+ pix_sz, rng, ray_org, ray_dir);
+ if(res != RES_OK) continue; /* Reject the current sample */
+
+ r0 = ssp_rng_canonical(rng);
+ r1 = ssp_rng_canonical(rng);
+ r2 = ssp_rng_canonical(rng);
+
+ /* Sample a wavelength */
+ 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, 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,
+ HTRDR_RADIANCE_ALL, 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;
+
+ /* Update the pixel accumulator of the current channel */
+ radiance.sum_weights += weight;
+ radiance.sum_weights_sqr += weight*weight;
+ radiance.nweights += 1;
+
+ /* Update the pixel accumulator of per realisation time */
+ time.sum_weights += usec;
+ time.sum_weights_sqr += usec*usec;
+ time.nweights += 1;
+ }
+
+ /* Compute the estimation of the pixel radiance */
+ htrdr_accum_get_estimation(&radiance, &pixel->radiance);
+
+ /* Save the per realisation integration time */
+ pixel->time = time;
+
+ /* Compute the brightness_temperature of the pixel and estimate its standard
+ * 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
+draw_tile
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ const int64_t tile_mcode, /* For debug only */
+ const size_t tile_org[2], /* Origin of the tile in pixel space */
+ const size_t tile_sz[2], /* Definition of the tile */
+ const double pix_sz[2], /* Size of a pixel in the normalized image plane */
+ const struct htrdr_sensor* sensor,
+ const size_t spp, /* #samples per pixel */
+ struct ssp_rng* rng,
+ struct tile* tile)
+{
+ size_t npixels;
+ size_t mcode; /* Morton code of tile pixel */
+ ASSERT(htrdr && tile_org && tile_sz && pix_sz && sensor && spp && tile);
+ (void)tile_mcode;
+ /* Adjust the #pixels to process them wrt a morton order */
+ npixels = round_up_pow2(MMAX(tile_sz[0], tile_sz[1]));
+ npixels *= npixels;
+
+ FOR_EACH(mcode, 0, npixels) {
+ union pixel* pixel;
+ size_t ipix_tile[2]; /* Pixel coord in the tile */
+ size_t ipix[2]; /* Pixel coord in the buffer */
+
+ ipix_tile[0] = morton2D_decode((uint32_t)(mcode>>0));
+ if(ipix_tile[0] >= tile_sz[0]) continue; /* Pixel is out of tile */
+ ipix_tile[1] = morton2D_decode((uint32_t)(mcode>>1));
+ if(ipix_tile[1] >= tile_sz[1]) continue; /* Pixel is out of tile */
+
+ /* Fetch and reset the pixel accumulator */
+ pixel = tile_at(tile, ipix_tile[0], ipix_tile[1]);
+
+ /* Compute the pixel coordinate */
+ ipix[0] = tile_org[0] + ipix_tile[0];
+ ipix[1] = tile_org[1] + ipix_tile[1];
+
+ /* Draw the pixel */
+ switch(sensor->type) {
+ case HTRDR_SENSOR_RECTANGLE:
+ draw_pixel_flux
+ (htrdr, ithread, ipix, pix_sz, sensor, spp, rng, &pixel->flux);
+ break;
+ case HTRDR_SENSOR_CAMERA:
+ switch(htrdr->spectral_type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ draw_pixel_xwave
+ (htrdr, ithread, ipix, pix_sz, sensor, spp, rng, &pixel->xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ draw_pixel_image
+ (htrdr, ithread, ipix, pix_sz, sensor->camera, spp, rng, &pixel->image);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+ return RES_OK;
+}
+
+static res_T
+draw_image
+ (struct htrdr* htrdr,
+ const struct htrdr_sensor* sensor,
+ const size_t width, /* Image width */
+ const size_t height, /* Image height */
+ const size_t spp,
+ const size_t ntiles_x,
+ const size_t ntiles_y,
+ const size_t ntiles_adjusted,
+ const double pix_sz[2], /* Pixel size in the normalized image plane */
+ struct proc_work* work,
+ struct list_node* tiles)
+{
+ struct ssp_rng* rng_proc = NULL;
+ size_t nthreads = 0;
+ size_t nthieves = 0;
+ size_t proc_ntiles = 0;
+ enum pixel_format pixfmt;
+ ATOMIC nsolved_tiles = 0;
+ ATOMIC res = RES_OK;
+ ASSERT(htrdr && sensor && spp && ntiles_adjusted && work && tiles);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+ ASSERT(width && height);
+ (void)ntiles_x, (void)ntiles_y;
+
+ res = ssp_rng_create(htrdr->allocator, &ssp_rng_mt19937_64, &rng_proc);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr, "could not create the RNG used to sample a process "
+ "to steal -- %s.\n", res_to_cstr((res_T)res));
+ goto error;
+ }
+
+ proc_ntiles = proc_work_get_ntiles(work);
+ nthreads = MMIN(htrdr->nthreads, proc_ntiles);
+
+ /* The process is not considered as a working process for himself */
+ htrdr->mpi_working_procs[htrdr->mpi_rank] = 0;
+ --htrdr->mpi_nworking_procs;
+
+ pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
+
+ omp_set_num_threads((int)nthreads);
+ #pragma omp parallel
+ for(;;) {
+ const int ithread = omp_get_thread_num();
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng* rng;
+ struct tile* tile;
+ uint32_t mcode = TILE_MCODE_NULL;
+ size_t tile_org[2];
+ size_t tile_sz[2];
+ size_t n;
+ res_T res_local = RES_OK;
+ int32_t pcent;
+
+ /* Get a tile to draw */
+ #pragma omp critical
+ {
+ mcode = proc_work_get_tile(work);
+ if(mcode == TILE_MCODE_NULL) { /* No more work on this process */
+ /* Try to steal works to concurrent processes */
+ proc_work_reset(work);
+ nthieves = mpi_steal_work(htrdr, rng_proc, work);
+ if(nthieves != 0) {
+ mcode = proc_work_get_tile(work);
+ }
+ }
+ }
+ if(mcode == TILE_MCODE_NULL) break; /* No more work */
+
+ /* Decode the morton code to retrieve the tile index */
+ tile_org[0] = morton2D_decode((uint32_t)(mcode>>0));
+ tile_org[1] = morton2D_decode((uint32_t)(mcode>>1));
+ ASSERT(tile_org[0] < ntiles_x && tile_org[1] < ntiles_y);
+
+ /* Create the tile */
+ tile = tile_create(htrdr->allocator, pixfmt);
+ if(!tile) {
+ ATOMIC_SET(&res, RES_MEM_ERR);
+ htrdr_log_err(htrdr,
+ "could not allocate the memory space of the tile (%lu, %lu) -- %s.\n",
+ (unsigned long)tile_org[0], (unsigned long)tile_org[1],
+ res_to_cstr((res_T)ATOMIC_GET(&res)));
+ break;
+ }
+
+ /* Register the tile */
+ #pragma omp critical
+ list_add_tail(tiles, &tile->node);
+
+ tile->data.x = (uint16_t)tile_org[0];
+ tile->data.y = (uint16_t)tile_org[1];
+
+ /* Define the tile origin in pixel space */
+ tile_org[0] *= TILE_SIZE;
+ tile_org[1] *= TILE_SIZE;
+
+ /* Compute the size of the tile clamped by the borders of the buffer */
+ tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]);
+ tile_sz[1] = MMIN(TILE_SIZE, height - tile_org[1]);
+
+ /* Create a proxy RNG for the current tile. This proxy is used for the
+ * current thread only and thus it has to manage only one RNG. This proxy
+ * is initialised in order to ensure that an unique and predictable set of
+ * random numbers is used for the current tile. */
+ SSP(rng_proxy_create2
+ (&htrdr->lifo_allocators[ithread],
+ &ssp_rng_threefry,
+ RNG_SEQUENCE_SIZE * (size_t)mcode, /* Offset */
+ RNG_SEQUENCE_SIZE, /* Size */
+ RNG_SEQUENCE_SIZE * (size_t)ntiles_adjusted, /* Pitch */
+ 1, &rng_proxy));
+ SSP(rng_proxy_create_rng(rng_proxy, 0, &rng));
+
+ /* Launch the tile rendering */
+ res_local = draw_tile(htrdr, (size_t)ithread, mcode, tile_org, tile_sz,
+ pix_sz, sensor, spp, rng, tile);
+
+ SSP(rng_proxy_ref_put(rng_proxy));
+ SSP(rng_ref_put(rng));
+
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ break;
+ }
+
+ /* Update the progress status */
+ n = (size_t)ATOMIC_INCR(&nsolved_tiles);
+ pcent = (int32_t)((double)n * 100.0 / (double)proc_ntiles + 0.5/*round*/);
+
+ #pragma omp critical
+ if(pcent > htrdr->mpi_progress_render[0]) {
+ htrdr->mpi_progress_render[0] = pcent;
+ if(htrdr->mpi_rank == 0) {
+ update_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING);
+ } else { /* Send the progress percentage to the master process */
+ send_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING, pcent);
+ }
+ }
+ }
+
+ if(ATOMIC_GET(&res) != RES_OK) goto error;
+
+ /* Asynchronously wait for processes completion. Use an asynchronous barrier to
+ * avoid a dead lock with the `mpi_probe_thieves' thread that requires also
+ * the `mpi_mutex'. */
+ {
+ MPI_Request req;
+
+ mutex_lock(htrdr->mpi_mutex);
+ MPI(Ibarrier(MPI_COMM_WORLD, &req));
+ mutex_unlock(htrdr->mpi_mutex);
+
+ mpi_wait_for_request(htrdr, &req);
+ }
+
+exit:
+ if(rng_proc) SSP(rng_ref_put(rng_proc));
+ return (res_T)res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_draw_map
+ (struct htrdr* htrdr,
+ const struct htrdr_sensor* sensor,
+ const size_t width,
+ const size_t height,
+ const size_t spp,
+ struct htrdr_buffer* buf)
+{
+ char strbuf[128];
+ struct time t0, t1;
+ struct list_node tiles;
+ size_t ntiles_x, ntiles_y, ntiles, ntiles_adjusted;
+ size_t itile;
+ struct proc_work work;
+ struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL;
+ size_t proc_ntiles_adjusted;
+ double pix_sz[2];
+ ATOMIC probe_thieves = 1;
+ ATOMIC res = RES_OK;
+ ASSERT(htrdr && sensor && width && height);
+ ASSERT(htrdr->mpi_rank != 0 || buf);
+
+ list_init(&tiles);
+ proc_work_init(htrdr->allocator, &work);
+
+ if(htrdr->mpi_rank == 0) {
+ const size_t pixsz = htrdr_spectral_type_get_pixsz
+ (htrdr->spectral_type, sensor->type);
+ const size_t pixal = htrdr_spectral_type_get_pixal
+ (htrdr->spectral_type, sensor->type);
+
+ htrdr_buffer_get_layout(buf, &layout);
+ ASSERT(layout.width || layout.height || layout.elmt_size);
+ ASSERT(layout.width == width && layout.height == height);
+
+ if(layout.elmt_size != pixsz || layout.alignment < pixal) {
+ htrdr_log_err(htrdr, "%s: invalid buffer layout.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ /* Compute the overall number of tiles */
+ ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
+ ntiles_y = (height+ (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
+ ntiles = ntiles_x * ntiles_y;
+
+ /* Compute the pixel size in the normalized image plane */
+ pix_sz[0] = 1.0 / (double)width;
+ pix_sz[1] = 1.0 / (double)height;
+
+ /* Adjust the #tiles for the morton-encoding procedure */
+ ntiles_adjusted = round_up_pow2(MMAX(ntiles_x, ntiles_y));
+ ntiles_adjusted *= ntiles_adjusted;
+
+ /* Define the initial number of tiles of the current process */
+ proc_ntiles_adjusted = ntiles_adjusted / (size_t)htrdr->mpi_nprocs;
+ if(htrdr->mpi_rank == 0) { /* Affect the remaining tiles to the master proc */
+ proc_ntiles_adjusted +=
+ ntiles_adjusted - proc_ntiles_adjusted*(size_t)htrdr->mpi_nprocs;
+ }
+
+ /* Define the initial list of tiles of the process */
+ FOR_EACH(itile, 0, proc_ntiles_adjusted) {
+ uint32_t mcode;
+ uint16_t tile_org[2];
+
+ mcode = (uint32_t)itile*(uint32_t)htrdr->mpi_nprocs
+ + (uint32_t)htrdr->mpi_rank;
+
+ tile_org[0] = morton2D_decode(mcode>>0);
+ if(tile_org[0] >= ntiles_x) continue;
+ tile_org[1] = morton2D_decode(mcode>>1);
+ if(tile_org[1] >= ntiles_y) continue;
+ proc_work_add_tile(&work, mcode);
+ }
+
+ if(htrdr->mpi_rank == 0) {
+ fetch_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING);
+ print_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING);
+ }
+
+ time_current(&t0);
+
+ omp_set_nested(1); /* Enable nested threads for draw_image */
+ #pragma omp parallel sections num_threads(2)
+ {
+ #pragma omp section
+ mpi_probe_thieves(htrdr, &work, &probe_thieves);
+
+ #pragma omp section
+ {
+ draw_image(htrdr, sensor, width, height, spp, ntiles_x, ntiles_y,
+ ntiles_adjusted, pix_sz, &work, &tiles);
+ /* The processes have no more work to do. Stop probing for thieves */
+ ATOMIC_SET(&probe_thieves, 0);
+ }
+ }
+
+ if(htrdr->mpi_rank == 0) {
+ update_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING);
+ fprintf(stderr, "\n"); /* Add a new line after the progress statuses */
+ }
+
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, strbuf, sizeof(strbuf));
+ htrdr_log(htrdr, "Rendering time: %s\n", strbuf);
+
+ /* Gather accum buffers from the group of processes */
+ time_current(&t0);
+ res = mpi_gather_tiles(htrdr, sensor, buf, ntiles, &tiles);
+ if(res != RES_OK) goto error;
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, strbuf, sizeof(strbuf));
+ htrdr_log(htrdr, "Image gathering time: %s\n", strbuf);
+
+exit:
+ { /* Free allocated tiles */
+ struct list_node* node;
+ struct list_node* tmp;
+ LIST_FOR_EACH_SAFE(node, tmp, &tiles) {
+ struct tile* tile = CONTAINER_OF(node, struct tile, node);
+ list_del(node);
+ tile_ref_put(tile);
+ }
+ }
+ proc_work_release(&work);
+ return (res_T)res;
+error:
+ goto exit;
+}
+
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -1,1081 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#define _POSIX_C_SOURCE 200112L /* nanosleep && nextafter */
-
-#include "htrdr.h"
-#include "htrdr_c.h"
-#include "htrdr_buffer.h"
-#include "htrdr_camera.h"
-#include "htrdr_cie_xyz.h"
-#include "htrdr_ran_wlen.h"
-#include "htrdr_solve.h"
-
-#include <high_tune/htsky.h>
-
-#include <rsys/algorithm.h>
-#include <rsys/clock_time.h>
-#include <rsys/cstr.h>
-#include <rsys/dynamic_array_u32.h>
-#include <rsys/math.h>
-#include <rsys/mutex.h>
-#include <star/ssp.h>
-
-#include <omp.h>
-#include <mpi.h>
-#include <time.h>
-#include <unistd.h>
-
-#define RNG_SEQUENCE_SIZE 10000
-
-#define TILE_MCODE_NULL UINT32_MAX
-#define TILE_SIZE 32 /* Definition in X & Y of a tile */
-STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
-
-enum pixel_format {
- PIXEL_XWAVE,
- PIXEL_IMAGE
-};
-
-union pixel {
- struct htrdr_pixel_xwave xwave;
- struct htrdr_pixel_image image;
-};
-
-/* Tile of row ordered image pixels */
-struct tile {
- struct list_node node;
- struct mem_allocator* allocator;
- ref_T ref;
-
- struct tile_data {
- uint16_t x, y; /* 2D coordinates of the tile in tile space */
- enum pixel_format format;
- /* Simulate the flexible array member of the C99 standard. */
- union pixel pixels[1/*Dummy element*/];
- } data;
-};
-
-/* List of tile to compute onto the MPI process. */
-struct proc_work {
- struct mutex* mutex;
- struct darray_u32 tiles; /* #tiles to render */
- size_t itile; /* Next tile to render in the above list of tiles */
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static FINLINE uint16_t
-morton2D_decode(const uint32_t u32)
-{
- uint32_t x = u32 & 0x55555555;
- x = (x | (x >> 1)) & 0x33333333;
- x = (x | (x >> 2)) & 0x0F0F0F0F;
- x = (x | (x >> 4)) & 0x00FF00FF;
- x = (x | (x >> 8)) & 0x0000FFFF;
- return (uint16_t)x;
-}
-
-static FINLINE uint32_t
-morton2D_encode(const uint16_t u16)
-{
- uint32_t u32 = u16;
- u32 = (u32 | (u32 << 8)) & 0x00FF00FF;
- u32 = (u32 | (u32 << 4)) & 0X0F0F0F0F;
- u32 = (u32 | (u32 << 2)) & 0x33333333;
- u32 = (u32 | (u32 << 1)) & 0x55555555;
- 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)
-{
- struct tile* tile;
- const size_t tile_sz =
- sizeof(struct tile) - sizeof(union pixel)/*rm dummy pixel*/;
- const size_t buf_sz = /* Flexiblbe array element */
- TILE_SIZE*TILE_SIZE*sizeof(union pixel);
- ASSERT(allocator);
-
- tile = MEM_ALLOC(allocator, tile_sz+buf_sz);
- if(!tile) return NULL;
-
- tile->data.format = fmt;
-
- ref_init(&tile->ref);
- list_init(&tile->node);
- tile->allocator = allocator;
- ASSERT(IS_ALIGNED(&tile->data.pixels, ALIGNOF(union pixel)));
-
- return tile;
-}
-
-static INLINE void
-tile_ref_get(struct tile* tile)
-{
- ASSERT(tile);
- tile_ref_get(tile);
-}
-
-static INLINE void
-release_tile(ref_T* ref)
-{
- struct tile* tile = CONTAINER_OF(ref, struct tile, ref);
- ASSERT(ref);
- MEM_RM(tile->allocator, tile);
-}
-
-static INLINE void
-tile_ref_put(struct tile* tile)
-{
- ASSERT(tile);
- ref_put(&tile->ref, release_tile);
-}
-
-static FINLINE union pixel*
-tile_at
- (struct tile* tile,
- const size_t x, /* In tile space */
- const size_t y) /* In tile space */
-{
- ASSERT(tile && x < TILE_SIZE && y < TILE_SIZE);
- return tile->data.pixels + (y*TILE_SIZE + x);
-}
-
-static void
-write_tile_data
- (struct htrdr* htrdr,
- struct htrdr_buffer* buf,
- const struct tile_data* tile_data)
-{
- struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL;
- size_t icol, irow;
- size_t irow_tile;
- size_t ncols_tile, nrows_tile;
- char* buf_mem;
- ASSERT(htrdr && buf && tile_data);
- (void)htrdr;
-
- htrdr_buffer_get_layout(buf, &layout);
- buf_mem = htrdr_buffer_get_data(buf);
- 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;
- irow = tile_data->y * (size_t)TILE_SIZE;
-
- /* Define the number of tile row/columns to write into the buffer */
- ncols_tile = MMIN(icol + TILE_SIZE, layout.width) - icol;
- nrows_tile = MMIN(irow + TILE_SIZE, layout.height) - irow;
-
- /* Copy the row ordered tile data */
- FOR_EACH(irow_tile, 0, nrows_tile) {
- char* buf_row = buf_mem + (irow + irow_tile) * layout.pitch;
- char* buf_col = buf_row + icol * layout.elmt_size;
- const union pixel* tile_row = tile_data->pixels + irow_tile*TILE_SIZE;
- size_t x;
-
- FOR_EACH(x, 0, ncols_tile) {
- switch(tile_data->format) {
- case PIXEL_XWAVE:
- ((struct htrdr_pixel_xwave*)buf_col)[x] = tile_row[x].xwave;
- break;
- case PIXEL_IMAGE:
- ((struct htrdr_pixel_image*)buf_col)[x] = tile_row[x].image;
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- }
- }
-}
-
-static INLINE void
-proc_work_init(struct mem_allocator* allocator, struct proc_work* work)
-{
- ASSERT(work);
- darray_u32_init(allocator, &work->tiles);
- work->itile = 0;
- CHK(work->mutex = mutex_create());
-}
-
-static INLINE void
-proc_work_release(struct proc_work* work)
-{
- darray_u32_release(&work->tiles);
- mutex_destroy(work->mutex);
-}
-
-static INLINE void
-proc_work_reset(struct proc_work* work)
-{
- ASSERT(work);
- mutex_lock(work->mutex);
- darray_u32_clear(&work->tiles);
- work->itile = 0;
- mutex_unlock(work->mutex);
-}
-
-static INLINE void
-proc_work_add_tile(struct proc_work* work, const uint32_t mcode)
-{
- mutex_lock(work->mutex);
- CHK(darray_u32_push_back(&work->tiles, &mcode) == RES_OK);
- mutex_unlock(work->mutex);
-}
-
-static INLINE uint32_t
-proc_work_get_tile(struct proc_work* work)
-{
- uint32_t mcode;
- ASSERT(work);
- mutex_lock(work->mutex);
- if(work->itile >= darray_u32_size_get(&work->tiles)) {
- mcode = TILE_MCODE_NULL;
- } else {
- mcode = darray_u32_cdata_get(&work->tiles)[work->itile];
- ++work->itile;
- }
- mutex_unlock(work->mutex);
- return mcode;
-}
-
-static INLINE size_t
-proc_work_get_ntiles(struct proc_work* work)
-{
- size_t sz = 0;
- ASSERT(work);
- mutex_lock(work->mutex);
- sz = darray_u32_size_get(&work->tiles);
- mutex_unlock(work->mutex);
- return sz;
-}
-
-static void
-mpi_wait_for_request(struct htrdr* htrdr, MPI_Request* req)
-{
- ASSERT(htrdr && req);
-
- /* Wait for process synchronisation */
- for(;;) {
- struct timespec t;
- int complete;
- t.tv_sec = 0;
- t.tv_nsec = 10000000; /* 10ms */
-
- mutex_lock(htrdr->mpi_mutex);
- MPI(Test(req, &complete, MPI_STATUS_IGNORE));
- mutex_unlock(htrdr->mpi_mutex);
- if(complete) break;
-
- nanosleep(&t, NULL);
- }
-}
-
-static void
-mpi_probe_thieves
- (struct htrdr* htrdr,
- struct proc_work* work,
- ATOMIC* probe_thieves)
-{
- uint32_t tiles[UINT8_MAX];
- struct timespec t;
- ASSERT(htrdr && work && probe_thieves);
-
- if(htrdr->mpi_nprocs == 1) /* The process is alone. No thief is possible */
- return;
-
- t.tv_sec = 0;
-
- /* Protect MPI calls of multiple invocations from concurrent threads */
- #define P_MPI(Func) { \
- mutex_lock(htrdr->mpi_mutex); \
- MPI(Func); \
- mutex_unlock(htrdr->mpi_mutex); \
- } (void)0
-
- while(ATOMIC_GET(probe_thieves)) {
- MPI_Status status;
- size_t itile;
- int msg;
-
- /* Probe if a steal request was submitted by any processes */
- P_MPI(Iprobe(MPI_ANY_SOURCE, HTRDR_MPI_STEAL_REQUEST, MPI_COMM_WORLD, &msg,
- &status));
-
- if(msg) { /* A steal request was posted */
- MPI_Request req;
- uint8_t ntiles_to_steal;
-
- /* Asynchronously receive the steal request */
- P_MPI(Irecv(&ntiles_to_steal, 1, MPI_UINT8_T, status.MPI_SOURCE,
- HTRDR_MPI_STEAL_REQUEST, MPI_COMM_WORLD, &req));
-
- /* Wait for the completion of the steal request */
- mpi_wait_for_request(htrdr, &req);
-
- /* Thief some tiles */
- FOR_EACH(itile, 0, ntiles_to_steal) {
- tiles[itile] = proc_work_get_tile(work);
- }
- P_MPI(Send(&tiles, ntiles_to_steal, MPI_UINT32_T, status.MPI_SOURCE,
- HTRDR_MPI_WORK_STEALING, MPI_COMM_WORLD));
- }
- t.tv_nsec = 500000000; /* 500ms */
- nanosleep(&t, NULL);
- }
- #undef P_MPI
-}
-
-static int
-mpi_sample_working_process(struct htrdr* htrdr, struct ssp_rng* rng)
-{
- int iproc, i;
- int dst_rank;
- ASSERT(htrdr && rng && htrdr->mpi_nworking_procs);
-
- /* Sample the index of the 1st active process */
- iproc = (int)(ssp_rng_canonical(rng) * (double)htrdr->mpi_nworking_procs);
-
- /* Find the rank of the sampled active process. Use a simple linear search
- * since the overall number of processes should be quite low; at most few
- * dozens. */
- i = 0;
- FOR_EACH(dst_rank, 0, htrdr->mpi_nprocs) {
- if(htrdr->mpi_working_procs[dst_rank] == 0) continue; /* Inactive process */
- if(i == iproc) break; /* The rank of the sampled process is found */
- ++i;
- }
- ASSERT(dst_rank < htrdr->mpi_nprocs);
- return dst_rank;
-}
-
-/* Return the number of stolen tiles */
-static size_t
-mpi_steal_work
- (struct htrdr* htrdr,
- struct ssp_rng* rng,
- struct proc_work* work)
-{
- MPI_Request req;
- size_t itile;
- size_t nthieves = 0;
- uint32_t tiles[UINT8_MAX]; /* Morton code of the stolen tile */
- int proc_to_steal; /* Process to steal */
- uint8_t ntiles_to_steal = MMIN((uint8_t)(htrdr->nthreads*2), 16);
- ASSERT(htrdr && rng && work && htrdr->nthreads < UINT8_MAX);
-
- /* Protect MPI calls of multiple invocations from concurrent threads */
- #define P_MPI(Func) { \
- mutex_lock(htrdr->mpi_mutex); \
- MPI(Func); \
- mutex_unlock(htrdr->mpi_mutex); \
- } (void)0
-
- /* No more working process => nohting to steal */
- if(!htrdr->mpi_nworking_procs) return 0;
-
- /* Sample a process to steal */
- proc_to_steal = mpi_sample_working_process(htrdr, rng);
-
- /* Send a steal request to the sampled process and wait for a response */
- P_MPI(Send(&ntiles_to_steal, 1, MPI_UINT8_T, proc_to_steal,
- HTRDR_MPI_STEAL_REQUEST, MPI_COMM_WORLD));
-
- /* Receive the stolen tile from the sampled process */
- P_MPI(Irecv(tiles, ntiles_to_steal, MPI_UINT32_T, proc_to_steal,
- HTRDR_MPI_WORK_STEALING, MPI_COMM_WORLD, &req));
-
- mpi_wait_for_request(htrdr, &req);
-
- FOR_EACH(itile, 0, ntiles_to_steal) {
- if(tiles[itile] == TILE_MCODE_NULL) {
- ASSERT(htrdr->mpi_working_procs[proc_to_steal] != 0);
- htrdr->mpi_working_procs[proc_to_steal] = 0;
- htrdr->mpi_nworking_procs--;
- break;
- }
- proc_work_add_tile(work, tiles[itile]);
- ++nthieves;
- }
- #undef P_MPI
- return nthieves;
-}
-
-static res_T
-mpi_gather_tiles
- (struct htrdr* htrdr,
- struct htrdr_buffer* buf,
- const size_t ntiles,
- struct list_node* tiles)
-{
- /* Compute the size of the tile_data */
- const size_t msg_sz =
- sizeof(struct tile_data) - sizeof(union pixel)/*dummy*/
- + TILE_SIZE*TILE_SIZE*sizeof(union pixel);
-
- struct list_node* node = NULL;
- struct tile* tile = NULL;
- res_T res = RES_OK;
- ASSERT(htrdr && tiles);
- ASSERT(htrdr->mpi_rank != 0 || buf);
- (void)ntiles;
-
- if(htrdr->mpi_rank != 0) { /* Non master process */
- /* Send the computed tile to the master process */
- LIST_FOR_EACH(node, tiles) {
- struct tile* t = CONTAINER_OF(node, struct tile, node);
- MPI(Send(&t->data, (int)msg_sz, MPI_CHAR, 0,
- HTRDR_MPI_TILE_DATA, MPI_COMM_WORLD));
- }
- } else { /* Master process */
- size_t itile = 0;
-
- LIST_FOR_EACH(node, tiles) {
- struct tile* t = CONTAINER_OF(node, struct tile, node);
- 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 */
- pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
- tile = tile_create(htrdr->allocator, pixfmt);
- if(!tile) {
- res = RES_MEM_ERR;
- htrdr_log_err(htrdr,
- "could not allocate the temporary tile used to gather the process "
- "output data -- %s.\n", res_to_cstr(res));
- goto error;
- }
-
- /* Receive the tile data of the concurrent MPI processes */
- 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(htrdr, buf, &tile->data);
- }
- }
- }
-
-exit:
- if(tile) tile_ref_put(tile);
- return res;
-error:
- goto exit;
-}
-
-static void
-draw_pixel_image
- (struct htrdr* htrdr,
- const size_t ithread,
- const size_t ipix[2],
- const double pix_sz[2], /* Size of a pixel in the normalized image plane */
- const struct htrdr_camera* cam,
- const size_t spp,
- struct ssp_rng* rng,
- 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(htrdr->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ);
-
- /* Reset accumulators */
- XYZ[0] = HTRDR_ACCUM_NULL;
- XYZ[1] = HTRDR_ACCUM_NULL;
- XYZ[2] = HTRDR_ACCUM_NULL;
- time = HTRDR_ACCUM_NULL;
-
- FOR_EACH(ichannel, 0, 3) {
- size_t isamp;
-
- FOR_EACH(isamp, 0, spp) {
- struct time t0, t1;
- double pix_samp[2];
- double ray_org[3];
- double ray_dir[3];
- 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;
-
- /* Begin the registration of the time spent to in the realisation */
- time_current(&t0);
-
- /* Sample a position into the pixel, in the normalized image plane */
- pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
- pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
-
- /* Generate a ray starting from the pinhole camera and passing through the
- * pixel sample */
- htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir);
-
- r0 = ssp_rng_canonical(rng);
- r1 = ssp_rng_canonical(rng);
- r2 = ssp_rng_canonical(rng);
-
- /* Sample a spectral band and a quadrature point */
- switch(ichannel) {
- case 0: wlen = htrdr_cie_xyz_sample_X(htrdr->cie, r0, 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 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;
-
- /* Update the pixel accumulator of the current channel */
- XYZ[ichannel].sum_weights += weight;
- XYZ[ichannel].sum_weights_sqr += weight*weight;
- XYZ[ichannel].nweights += 1;
-
- /* Update the pixel accumulator of per realisation time */
- time.sum_weights += usec;
- time.sum_weights_sqr += usec*usec;
- time.nweights += 1;
- }
- }
-
- /* Flush pixel data */
- htrdr_accum_get_estimation(XYZ+0, &pixel->X);
- htrdr_accum_get_estimation(XYZ+1, &pixel->Y);
- htrdr_accum_get_estimation(XYZ+2, &pixel->Z);
- 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_xwave
- (struct htrdr* htrdr,
- const size_t ithread,
- const size_t ipix[2],
- const double pix_sz[2], /* Size of a pixel in the normalized image plane */
- const struct htrdr_camera* cam,
- const size_t spp,
- struct ssp_rng* rng,
- 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(htrdr->spectral_type == HTRDR_SPECTRAL_LW
- || htrdr->spectral_type == HTRDR_SPECTRAL_SW);
-
- /* Reset the pixel accumulators */
- radiance = HTRDR_ACCUM_NULL;
- time = HTRDR_ACCUM_NULL;
-
- FOR_EACH(isamp, 0, spp) {
- struct time t0, t1;
- double pix_samp[2];
- double ray_org[3];
- double ray_dir[3];
- double weight;
- double r0, r1, r2;
- double wlen;
- size_t iband;
- size_t iquad;
- double usec;
- double band_pdf;
-
- /* Begin the registration of the time spent in the realisation */
- time_current(&t0);
-
- /* Sample a position into the pixel, in the normalized image plane */
- pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
- pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
-
- /* Generate a ray starting from the pinhole camera and passing through the
- * pixel sample */
- htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir);
-
- r0 = ssp_rng_canonical(rng);
- r1 = ssp_rng_canonical(rng);
- r2 = ssp_rng_canonical(rng);
-
- /* Sample a wavelength */
- 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, 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;
-
- /* Update the pixel accumulator of the current channel */
- radiance.sum_weights += weight;
- radiance.sum_weights_sqr += weight*weight;
- radiance.nweights += 1;
-
- /* Update the pixel accumulator of per realisation time */
- time.sum_weights += usec;
- time.sum_weights_sqr += usec*usec;
- time.nweights += 1;
- }
-
- /* Compute the estimation of the pixel radiance */
- htrdr_accum_get_estimation(&radiance, &pixel->radiance);
-
- /* Save the per realisation integration time */
- pixel->time = time;
-
- /* Compute the brightness_temperature of the pixel and estimate its standard
- * 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
-draw_tile
- (struct htrdr* htrdr,
- const size_t ithread,
- const int64_t tile_mcode, /* For debug only */
- const size_t tile_org[2], /* Origin of the tile in pixel space */
- const size_t tile_sz[2], /* Definition of the tile */
- const double pix_sz[2], /* Size of a pixel in the normalized image plane */
- const struct htrdr_camera* cam,
- const size_t spp, /* #samples per pixel */
- struct ssp_rng* rng,
- struct tile* tile)
-{
- size_t npixels;
- size_t mcode; /* Morton code of tile pixel */
- ASSERT(htrdr && tile_org && tile_sz && pix_sz && cam && spp && tile);
- (void)tile_mcode;
- /* Adjust the #pixels to process them wrt a morton order */
- npixels = round_up_pow2(MMAX(tile_sz[0], tile_sz[1]));
- npixels *= npixels;
-
- FOR_EACH(mcode, 0, npixels) {
- union pixel* pixel;
- size_t ipix_tile[2]; /* Pixel coord in the tile */
- size_t ipix[2]; /* Pixel coord in the buffer */
-
- ipix_tile[0] = morton2D_decode((uint32_t)(mcode>>0));
- if(ipix_tile[0] >= tile_sz[0]) continue; /* Pixel is out of tile */
- ipix_tile[1] = morton2D_decode((uint32_t)(mcode>>1));
- if(ipix_tile[1] >= tile_sz[1]) continue; /* Pixel is out of tile */
-
- /* Fetch and reset the pixel accumulator */
- pixel = tile_at(tile, ipix_tile[0], ipix_tile[1]);
-
- /* Compute the pixel coordinate */
- ipix[0] = tile_org[0] + ipix_tile[0];
- ipix[1] = tile_org[1] + ipix_tile[1];
-
- /* Draw the pixel */
- 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;
-}
-
-static res_T
-draw_image
- (struct htrdr* htrdr,
- const struct htrdr_camera* cam,
- const size_t width, /* Image width */
- const size_t height, /* Image height */
- const size_t spp,
- const size_t ntiles_x,
- const size_t ntiles_y,
- const size_t ntiles_adjusted,
- const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct proc_work* work,
- struct list_node* tiles)
-{
- struct ssp_rng* rng_proc = NULL;
- size_t nthreads = 0;
- size_t nthieves = 0;
- size_t proc_ntiles = 0;
- enum pixel_format pixfmt;
- ATOMIC nsolved_tiles = 0;
- ATOMIC res = RES_OK;
- ASSERT(htrdr && cam && spp && ntiles_adjusted && work && tiles);
- ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
- ASSERT(width && height);
- (void)ntiles_x, (void)ntiles_y;
-
- res = ssp_rng_create(htrdr->allocator, &ssp_rng_mt19937_64, &rng_proc);
- if(res != RES_OK) {
- htrdr_log_err(htrdr, "could not create the RNG used to sample a process "
- "to steal -- %s.\n", res_to_cstr((res_T)res));
- goto error;
- }
-
- proc_ntiles = proc_work_get_ntiles(work);
- nthreads = MMIN(htrdr->nthreads, proc_ntiles);
-
- /* The process is not considered as a working process for himself */
- htrdr->mpi_working_procs[htrdr->mpi_rank] = 0;
- --htrdr->mpi_nworking_procs;
-
- pixfmt = spectral_type_to_pixfmt(htrdr->spectral_type);
-
- omp_set_num_threads((int)nthreads);
- #pragma omp parallel
- for(;;) {
- const int ithread = omp_get_thread_num();
- struct ssp_rng_proxy* rng_proxy = NULL;
- struct ssp_rng* rng;
- struct tile* tile;
- uint32_t mcode = TILE_MCODE_NULL;
- size_t tile_org[2];
- size_t tile_sz[2];
- size_t n;
- res_T res_local = RES_OK;
- int32_t pcent;
-
- /* Get a tile to draw */
- #pragma omp critical
- {
- mcode = proc_work_get_tile(work);
- if(mcode == TILE_MCODE_NULL) { /* No more work on this process */
- /* Try to steal works to concurrent processes */
- proc_work_reset(work);
- nthieves = mpi_steal_work(htrdr, rng_proc, work);
- if(nthieves != 0) {
- mcode = proc_work_get_tile(work);
- }
- }
- }
- if(mcode == TILE_MCODE_NULL) break; /* No more work */
-
- /* Decode the morton code to retrieve the tile index */
- tile_org[0] = morton2D_decode((uint32_t)(mcode>>0));
- tile_org[1] = morton2D_decode((uint32_t)(mcode>>1));
- ASSERT(tile_org[0] < ntiles_x && tile_org[1] < ntiles_y);
-
- /* Create the tile */
- tile = tile_create(htrdr->allocator, pixfmt);
- if(!tile) {
- ATOMIC_SET(&res, RES_MEM_ERR);
- htrdr_log_err(htrdr,
- "could not allocate the memory space of the tile (%lu, %lu) -- %s.\n",
- (unsigned long)tile_org[0], (unsigned long)tile_org[1],
- res_to_cstr((res_T)ATOMIC_GET(&res)));
- break;
- }
-
- /* Register the tile */
- #pragma omp critical
- list_add_tail(tiles, &tile->node);
-
- tile->data.x = (uint16_t)tile_org[0];
- tile->data.y = (uint16_t)tile_org[1];
-
- /* Define the tile origin in pixel space */
- tile_org[0] *= TILE_SIZE;
- tile_org[1] *= TILE_SIZE;
-
- /* Compute the size of the tile clamped by the borders of the buffer */
- tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]);
- tile_sz[1] = MMIN(TILE_SIZE, height - tile_org[1]);
-
- /* Create a proxy RNG for the current tile. This proxy is used for the
- * current thread only and thus it has to manage only one RNG. This proxy
- * is initialised in order to ensure that an unique and predictable set of
- * random numbers is used for the current tile. */
- SSP(rng_proxy_create2
- (&htrdr->lifo_allocators[ithread],
- &ssp_rng_threefry,
- RNG_SEQUENCE_SIZE * (size_t)mcode, /* Offset */
- RNG_SEQUENCE_SIZE, /* Size */
- RNG_SEQUENCE_SIZE * (size_t)ntiles_adjusted, /* Pitch */
- 1, &rng_proxy));
- SSP(rng_proxy_create_rng(rng_proxy, 0, &rng));
-
- /* Launch the tile rendering */
- res_local = draw_tile(htrdr, (size_t)ithread, mcode, tile_org, tile_sz,
- pix_sz, cam, spp, rng, tile);
-
- SSP(rng_proxy_ref_put(rng_proxy));
- SSP(rng_ref_put(rng));
-
- if(res_local != RES_OK) {
- ATOMIC_SET(&res, res_local);
- break;
- }
-
- /* Update the progress status */
- n = (size_t)ATOMIC_INCR(&nsolved_tiles);
- pcent = (int32_t)((double)n * 100.0 / (double)proc_ntiles + 0.5/*round*/);
-
- #pragma omp critical
- if(pcent > htrdr->mpi_progress_render[0]) {
- htrdr->mpi_progress_render[0] = pcent;
- if(htrdr->mpi_rank == 0) {
- update_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING);
- } else { /* Send the progress percentage to the master process */
- send_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING, pcent);
- }
- }
- }
-
- if(ATOMIC_GET(&res) != RES_OK) goto error;
-
- /* Asynchronously wait for processes completion. Use an asynchronous barrier to
- * avoid a dead lock with the `mpi_probe_thieves' thread that requires also
- * the `mpi_mutex'. */
- {
- MPI_Request req;
-
- mutex_lock(htrdr->mpi_mutex);
- MPI(Ibarrier(MPI_COMM_WORLD, &req));
- mutex_unlock(htrdr->mpi_mutex);
-
- mpi_wait_for_request(htrdr, &req);
- }
-
-exit:
- if(rng_proc) SSP(rng_ref_put(rng_proc));
- return (res_T)res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_draw_radiance
- (struct htrdr* htrdr,
- const struct htrdr_camera* cam,
- const size_t width,
- const size_t height,
- const size_t spp,
- struct htrdr_buffer* buf)
-{
- char strbuf[128];
- struct time t0, t1;
- struct list_node tiles;
- size_t ntiles_x, ntiles_y, ntiles, ntiles_adjusted;
- size_t itile;
- struct proc_work work;
- struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL;
- size_t proc_ntiles_adjusted;
- double pix_sz[2];
- ATOMIC probe_thieves = 1;
- ATOMIC res = RES_OK;
- ASSERT(htrdr && cam && width && height);
- ASSERT(htrdr->mpi_rank != 0 || buf);
-
- list_init(&tiles);
- proc_work_init(htrdr->allocator, &work);
-
- if(htrdr->mpi_rank == 0) {
- 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(layout.elmt_size != pixsz || layout.alignment < pixal) {
- htrdr_log_err(htrdr, "%s: invalid buffer layout.\n", FUNC_NAME);
- res = RES_BAD_ARG;
- goto error;
- }
- }
-
- /* Compute the overall number of tiles */
- ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
- ntiles_y = (height+ (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
- ntiles = ntiles_x * ntiles_y;
-
- /* Compute the pixel size in the normalized image plane */
- pix_sz[0] = 1.0 / (double)width;
- pix_sz[1] = 1.0 / (double)height;
-
- /* Adjust the #tiles for the morton-encoding procedure */
- ntiles_adjusted = round_up_pow2(MMAX(ntiles_x, ntiles_y));
- ntiles_adjusted *= ntiles_adjusted;
-
- /* Define the initial number of tiles of the current process */
- proc_ntiles_adjusted = ntiles_adjusted / (size_t)htrdr->mpi_nprocs;
- if(htrdr->mpi_rank == 0) { /* Affect the remaining tiles to the master proc */
- proc_ntiles_adjusted +=
- ntiles_adjusted - proc_ntiles_adjusted*(size_t)htrdr->mpi_nprocs;
- }
-
- /* Define the initial list of tiles of the process */
- FOR_EACH(itile, 0, proc_ntiles_adjusted) {
- uint32_t mcode;
- uint16_t tile_org[2];
-
- mcode = (uint32_t)itile*(uint32_t)htrdr->mpi_nprocs
- + (uint32_t)htrdr->mpi_rank;
-
- tile_org[0] = morton2D_decode(mcode>>0);
- if(tile_org[0] >= ntiles_x) continue;
- tile_org[1] = morton2D_decode(mcode>>1);
- if(tile_org[1] >= ntiles_y) continue;
- proc_work_add_tile(&work, mcode);
- }
-
- if(htrdr->mpi_rank == 0) {
- fetch_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING);
- print_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING);
- }
-
- time_current(&t0);
-
- omp_set_nested(1); /* Enable nested threads for draw_image */
- #pragma omp parallel sections num_threads(2)
- {
- #pragma omp section
- mpi_probe_thieves(htrdr, &work, &probe_thieves);
-
- #pragma omp section
- {
- draw_image(htrdr, cam, width, height, spp, ntiles_x, ntiles_y,
- ntiles_adjusted, pix_sz, &work, &tiles);
- /* The processes have no more work to do. Stop probing for thieves */
- ATOMIC_SET(&probe_thieves, 0);
- }
- }
-
- if(htrdr->mpi_rank == 0) {
- update_mpi_progress(htrdr, HTRDR_MPI_PROGRESS_RENDERING);
- fprintf(stderr, "\n"); /* Add a new line after the progress statuses */
- }
-
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, strbuf, sizeof(strbuf));
- htrdr_log(htrdr, "Rendering time: %s\n", strbuf);
-
- /* Gather accum buffers from the group of processes */
- time_current(&t0);
- res = mpi_gather_tiles(htrdr, buf, ntiles, &tiles);
- if(res != RES_OK) goto error;
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, strbuf, sizeof(strbuf));
- htrdr_log(htrdr, "Image gathering time: %s\n", strbuf);
-
-exit:
- { /* Free allocated tiles */
- struct list_node* node;
- struct list_node* tmp;
- LIST_FOR_EACH_SAFE(node, tmp, &tiles) {
- struct tile* tile = CONTAINER_OF(node, struct tile, node);
- list_del(node);
- tile_ref_put(tile);
- }
- }
- proc_work_release(&work);
- return (res_T)res;
-error:
- goto exit;
-}
-
diff --git a/src/htrdr_ground.c b/src/htrdr_ground.c
@@ -19,7 +19,7 @@
#include "htrdr.h"
#include "htrdr_interface.h"
#include "htrdr_ground.h"
-#include "htrdr_mtl.h"
+#include "htrdr_materials.h"
#include "htrdr_slab.h"
#include <aw.h>
@@ -165,9 +165,16 @@ parse_shape_interface
struct htrdr_interface* interf)
{
struct str str;
+ char* mtl_name0 = NULL;
+ char* mtl_name1 = NULL;
+ char* mtl_name2 = NULL;
char* mtl_name_front = NULL;
+ char* mtl_name_thin = NULL;
char* mtl_name_back = NULL;
char* tk_ctx = NULL;
+ int has_front = 0;
+ int has_thin = 0;
+ int has_back = 0;
res_T res = RES_OK;
ASSERT(htrdr && name && interf);
@@ -182,33 +189,70 @@ parse_shape_interface
goto error;
}
- /* Parse the name of the front/back faces */
- mtl_name_front = strtok_r(str_get(&str), ":", &tk_ctx);
- ASSERT(mtl_name_front); /* This can't be NULL */
- mtl_name_back = strtok_r(NULL, ":", &tk_ctx);
+ /* Reset the interface */
+ memset(interf, 0, sizeof(*interf));
+
+ /* Parse the name of the front/back/thin materials */
+ mtl_name0 = strtok_r(str_get(&str), ":", &tk_ctx);
+ mtl_name1 = strtok_r(NULL, ":", &tk_ctx);
+ mtl_name2 = strtok_r(NULL, ":", &tk_ctx);
+ ASSERT(mtl_name0); /* This can't be NULL */
+ mtl_name_front = mtl_name0;
+ if(mtl_name2) {
+ mtl_name_thin = mtl_name1;
+ mtl_name_back = mtl_name2;
+ } else {
+ mtl_name_thin = NULL;
+ mtl_name_back = mtl_name1;
+ }
+
if(!mtl_name_back) {
htrdr_log_err(htrdr,
- "The material name of the shape back faces are missing `%s'.\n", name);
+ "The back material name is missing `%s'.\n", name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Fetch the interface material if any */
+ if(mtl_name_thin) {
+ has_thin = htrdr_materials_find_mtl
+ (htrdr->mats, mtl_name_thin, &interf->mtl_thin);
+ if(!has_thin) {
+ htrdr_log_err(htrdr,
+ "Invalid interface `%s'. The interface material `%s' is unknown.\n",
+ name, mtl_name_thin);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ /* Fetch the front material */
+ has_front = htrdr_materials_find_mtl
+ (htrdr->mats, mtl_name_front, &interf->mtl_front);
+ if(!has_front) {
+ htrdr_log_err(htrdr,
+ "Invalid interface `%s'. The front material `%s' is unknown.\n",
+ name, mtl_name_front);
res = RES_BAD_ARG;
goto error;
}
- /* Fetch the front/back materials */
- interf->mtl_front = htrdr_mtl_get(htrdr->mtl, mtl_name_front);
- interf->mtl_back = htrdr_mtl_get(htrdr->mtl, mtl_name_back);
- if(!interf->mtl_front && !interf->mtl_back) {
+ /* Fetch the back material */
+ has_back = htrdr_materials_find_mtl
+ (htrdr->mats, mtl_name_back, &interf->mtl_back);
+ if(!has_back) {
htrdr_log_err(htrdr,
- "Invalid interface `%s:%s'. "
- "The front and the back materials are both uknown.\n",
- mtl_name_front, mtl_name_back);
+ "Invalid interface `%s'. The back material `%s' is unknown.\n",
+ name, mtl_name_back);
res = RES_BAD_ARG;
goto error;
}
+
exit:
str_release(&str);
return res;
error:
- interf->mtl_front = interf->mtl_back = NULL;
+ *interf = HTRDR_INTERFACE_NULL;
goto exit;
}
@@ -649,3 +693,61 @@ error:
goto exit;
}
+res_T
+htrdr_ground_find_closest_point
+ (struct htrdr_ground* ground,
+ const double pos[3],
+ const double radius,
+ struct s3d_hit* hit)
+{
+ float query_pos[3];
+ float query_radius;
+ float ground_sz[3];
+ res_T res = RES_OK;
+ ASSERT(ground && pos && hit);
+
+ if(!ground->view) { /* No ground geometry */
+ *hit = S3D_HIT_NULL;
+ goto exit;
+ }
+
+ query_radius = (float)radius;
+ f3_set_d3(query_pos, pos);
+
+ if(ground->repeat) {
+ float translation[2] = {0, 0};
+ int64_t xy[2];
+ ground_sz[0] = ground->upper[0] - ground->lower[0];
+ ground_sz[1] = ground->upper[1] - ground->lower[1];
+ ground_sz[2] = ground->upper[2] - ground->lower[2];
+
+ /* Define the 2D index of the current ground instance. (0,0) is the index
+ * of the non instantiated ground */
+ xy[0] = (int64_t)floor((query_pos[0] - ground->lower[0]) / ground_sz[0]);
+ xy[1] = (int64_t)floor((query_pos[1] - ground->lower[1]) / ground_sz[1]);
+
+ /* Define the translation along the X and Y axis from world space to local
+ * ground geometry space */
+ translation[0] = -(float)xy[0] * ground_sz[0];
+ translation[1] = -(float)xy[1] * ground_sz[1];
+
+ /* Transform the query pos in local ground geometry space */
+ query_pos[0] += translation[0];
+ query_pos[1] += translation[1];
+ }
+
+ /* Closest point query */
+ res = s3d_scene_view_closest_point
+ (ground->view, query_pos, query_radius, NULL, hit);
+ if(res != RES_OK) {
+ htrdr_log_err(ground->htrdr,
+ "%s: could not query the closest point to the ground geometry -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/htrdr_ground.h b/src/htrdr_ground.h
@@ -67,5 +67,12 @@ htrdr_ground_trace_ray
const struct s3d_hit* prev_hit,/* Previous hit. Avoid self hit. May be NULL*/
struct s3d_hit* hit);
+extern LOCAL_SYM res_T
+htrdr_ground_find_closest_point
+ (struct htrdr_ground* ground,
+ const double position[3],
+ const double radius,
+ struct s3d_hit* hit);
+
#endif /* HTRDR_GROUND_H */
diff --git a/src/htrdr_interface.c b/src/htrdr_interface.c
@@ -1,187 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "htrdr.h"
-#include "htrdr_interface.h"
-
-#include <modradurb/mrumtl.h>
-
-#include <star/s3d.h>
-#include <star/ssf.h>
-#include <star/ssp.h>
-
-#include <rsys/cstr.h>
-#include <rsys/double3.h>
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static res_T
-create_bsdf_diffuse
- (struct htrdr* htrdr,
- const struct mrumtl_brdf* brdf,
- const size_t ithread,
- struct ssf_bsdf** out_bsdf)
-{
- struct ssf_bsdf* bsdf = NULL;
- double reflectivity = 0;
- res_T res = RES_OK;
- ASSERT(htrdr && brdf && out_bsdf);
- ASSERT(mrumtl_brdf_get_type(brdf) == MRUMTL_BRDF_LAMBERTIAN);
- ASSERT(ithread < htrdr->nthreads);
-
- res = ssf_bsdf_create
- (&htrdr->lifo_allocators[ithread], &ssf_lambertian_reflection, &bsdf);
- if(res != RES_OK) goto error;
-
- reflectivity = mrumtl_brdf_lambertian_get_reflectivity(brdf);
- res = ssf_lambertian_reflection_setup(bsdf, reflectivity);
- if(res != RES_OK) goto error;
-
-exit:
- *out_bsdf = bsdf;
- return res;
-error:
- if(bsdf) { SSF(bsdf_ref_put(bsdf)); bsdf = NULL; }
- goto exit;
-}
-
-static res_T
-create_bsdf_specular
- (struct htrdr* htrdr,
- const struct mrumtl_brdf* brdf,
- const size_t ithread,
- struct ssf_bsdf** out_bsdf)
-{
- struct ssf_bsdf* bsdf = NULL;
- struct ssf_fresnel* fresnel = NULL;
- double reflectivity = 0;
- res_T res = RES_OK;
- ASSERT(htrdr && brdf && out_bsdf);
- ASSERT(mrumtl_brdf_get_type(brdf) == MRUMTL_BRDF_SPECULAR);
- ASSERT(ithread < htrdr->nthreads);
-
- res = ssf_bsdf_create
- (&htrdr->lifo_allocators[ithread], &ssf_specular_reflection, &bsdf);
- if(res != RES_OK) goto error;
-
- res = ssf_fresnel_create
- (&htrdr->lifo_allocators[ithread], &ssf_fresnel_constant, &fresnel);
- if(res != RES_OK) goto error;
-
- reflectivity = mrumtl_brdf_specular_get_reflectivity(brdf);
- res = ssf_fresnel_constant_setup(fresnel, reflectivity);
- if(res != RES_OK) goto error;
-
- res = ssf_specular_reflection_setup(bsdf, fresnel);
- if(res != RES_OK) goto error;
-
-exit:
- if(fresnel) SSF(fresnel_ref_put(fresnel));
- *out_bsdf = bsdf;
- return res;
-error:
- if(bsdf) { SSF(bsdf_ref_put(bsdf)); bsdf = NULL; }
- goto exit;
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_interface_create_bsdf
- (struct htrdr* htrdr,
- const struct htrdr_interface* interf,
- const size_t ithread,
- const double wavelength,
- const double pos[3],
- const double dir[3],
- struct ssp_rng* rng,
- struct s3d_hit* hit,
- struct ssf_bsdf** out_bsdf)
-{
- enum { FRONT, BACK };
- struct ssf_bsdf* bsdf = NULL;
- const struct mrumtl_brdf* brdf = NULL;
- const struct mrumtl* mat = NULL;
- double N[3];
- double r;
- int hit_side;
- res_T res = RES_OK;
- (void)pos;
- ASSERT(htrdr && pos && hit && out_bsdf);
- ASSERT(interf && (interf->mtl_front || interf->mtl_back));
-
- ASSERT(htrdr && interf && pos && dir && hit && out_bsdf);
- ASSERT(d3_is_normalized(dir));
-
- d3_normalize(N, d3_set_f3(N, hit->normal));
-
- hit_side = d3_dot(N, dir) < 0 ? FRONT : BACK;
-
- /* Retrieve the brdf of the material on the other side of the hit side */
- switch(hit_side) {
- case BACK: mat = interf->mtl_front; break;
- case FRONT: mat = interf->mtl_back; break;
- default: FATAL("Unreachable code.\n"); break;
- }
-
- /* Due to numerical issue the hit side might be wrong and thus the fetched
- * material might be undefined (e.g. semi-transparent materials). Handle this
- * issue by fetching the other material. */
- if(!mat) {
- switch(hit_side) {
- case BACK: mat = interf->mtl_back; break;
- case FRONT: mat = interf->mtl_front; break;
- default: FATAL("Unreachable code.\n"); break;
- }
- }
- ASSERT(mat);
-
- r = ssp_rng_canonical(rng);
-
- res = mrumtl_fetch_brdf2(mat, wavelength, r, &brdf);
- if(res != RES_OK) {
- htrdr_log_err(htrdr,
- "%s: error retreiving the MruMtl BRDF for the wavelength %g.\n",
- FUNC_NAME, wavelength);
- res = RES_BAD_ARG;
- goto error;
- }
-
- switch(mrumtl_brdf_get_type(brdf)) {
- case MRUMTL_BRDF_LAMBERTIAN:
- res = create_bsdf_diffuse(htrdr, brdf, ithread, &bsdf);
- break;
- case MRUMTL_BRDF_SPECULAR:
- res = create_bsdf_specular(htrdr, brdf, ithread, &bsdf);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- if(res != RES_OK) {
- htrdr_log_err(htrdr, "%s: could not create the BSDF -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
- }
-
-exit:
- *out_bsdf = bsdf;
- return res;
-error:
- if(bsdf) { SSF(bsdf_ref_put(bsdf)); bsdf = NULL; }
- goto exit;
-}
-
diff --git a/src/htrdr_interface.h b/src/htrdr_interface.h
@@ -17,7 +17,9 @@
#ifndef HTRDR_INTERFACE_H
#define HTRDR_INTERFACE_H
-#include <star/ssf.h>
+#include "htrdr_materials.h"
+#include <star/s3d.h>
+#include <rsys/double3.h>
/* Forward declaration of external data type */
struct mrumtl;
@@ -26,22 +28,55 @@ struct ssf_bsdf;
struct ssp_rng;
struct htrdr_interface {
- const struct mrumtl* mtl_front;
- const struct mrumtl* mtl_back;
+ struct htrdr_mtl mtl_front;
+ struct htrdr_mtl mtl_back;
+ struct htrdr_mtl mtl_thin; /* != NULL <=> thin material */
};
static const struct htrdr_interface HTRDR_INTERFACE_NULL;
-extern LOCAL_SYM res_T
-htrdr_interface_create_bsdf
- (struct htrdr* htrdr,
- const struct htrdr_interface* interf,
- const size_t ithread,
- const double wavelength,
- const double pos[3],
- const double dir[3], /* Normalized incoming direction */
- struct ssp_rng* rng,
- struct s3d_hit* hit,
- struct ssf_bsdf** bsdf);
+static INLINE const struct htrdr_mtl*
+htrdr_interface_fetch_hit_mtl
+ (const struct htrdr_interface* interf,
+ const double dir[3], /* Incoming ray */
+ struct s3d_hit* hit)
+{
+ const struct htrdr_mtl* mtl = NULL;
+ enum { FRONT, BACK };
+ ASSERT(interf && dir && d3_is_normalized(dir) && hit && !S3D_HIT_NONE(hit));
+ ASSERT(interf->mtl_front.mrumtl
+ || interf->mtl_back.mrumtl
+ || interf->mtl_thin.mrumtl);
+
+ if(interf->mtl_thin.mrumtl) {
+ mtl = &interf->mtl_thin;
+ } else {
+ double N[3];
+ int hit_side;
+ d3_normalize(N, d3_set_f3(N, hit->normal));
+ hit_side = d3_dot(N, dir) < 0 ? FRONT : BACK;
+
+ /* Retrieve the brdf of the material on the *other side* of the hit side */
+ switch(hit_side) {
+ case BACK: mtl = &interf->mtl_front; break;
+ case FRONT: mtl = &interf->mtl_back; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ /* Due to numerical issue the hit side might be wrong and thus the fetched
+ * material might be undefined (e.g. semi-transparent materials). Handle this
+ * issue by fetching the other material. */
+ if(!mtl->mrumtl) {
+ switch(hit_side) {
+ case BACK: mtl = &interf->mtl_back; break;
+ case FRONT: mtl = &interf->mtl_front; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+ ASSERT(mtl->mrumtl);
+ }
+
+ return mtl;
+}
#endif /* HTRDR_INTERFACE_H */
diff --git a/src/htrdr_materials.c b/src/htrdr_materials.c
@@ -0,0 +1,449 @@
+/* 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 /* strtok_r and wordexp support */
+
+#include "htrdr.h"
+#include "htrdr_materials.h"
+
+#include <modradurb/mrumtl.h>
+#include <star/ssf.h>
+#include <star/ssp.h>
+
+#include <rsys/cstr.h>
+#include <rsys/double3.h>
+#include <rsys/hash_table.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+#include <rsys/str.h>
+#include <rsys/text_reader.h>
+
+#include <string.h>
+#include <wordexp.h>
+
+struct mtl {
+ struct mrumtl* mrumtl;
+ double temperature; /* In Kelvin */
+};
+static const struct mtl MTL_NULL = {NULL, -1};
+
+/* Generate the hash table that maps a material name to its data */
+#define HTABLE_NAME name2mtl
+#define HTABLE_DATA struct mtl
+#define HTABLE_KEY struct str
+#define HTABLE_KEY_FUNCTOR_INIT str_init
+#define HTABLE_KEY_FUNCTOR_RELEASE str_release
+#define HTABLE_KEY_FUNCTOR_COPY str_copy
+#define HTABLE_KEY_FUNCTOR_COPY_AND_RELEASE str_copy_and_release
+#define HTABLE_KEY_FUNCTOR_HASH str_hash
+#define HTABLE_KEY_FUNCTOR_EQ str_eq
+#include <rsys/hash_table.h>
+
+struct htrdr_materials {
+ struct htable_name2mtl name2mtl;
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+static res_T
+parse_material
+ (struct htrdr_materials* mats,
+ struct txtrdr* txtrdr,
+ struct str* str) /* Scratch string */
+{
+ wordexp_t wexp;
+ char* tk = NULL;
+ char* tk_ctx = NULL;
+ struct mtl mtl = MTL_NULL;
+ int err = 0;
+ int wexp_is_allocated = 0;
+ res_T res = RES_OK;
+ ASSERT(mats && txtrdr);
+
+ tk = strtok_r(txtrdr_get_line(txtrdr), " \t", &tk_ctx);
+ ASSERT(tk);
+
+ res = str_set(str, tk);
+ if(res != RES_OK) {
+ htrdr_log_err(mats->htrdr,
+ "%s:%lu: could not copy the material name `%s' -- %s.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr), tk,
+ res_to_cstr(res));
+ goto error;
+ }
+
+ tk = strtok_r(NULL, "", &tk_ctx);
+ if(!tk) {
+ htrdr_log_err(mats->htrdr,
+ "%s:%lu: missing the MruMtl file for the material `%s'.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
+ str_cget(str));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ err = wordexp(tk, &wexp, 0);
+ if(err) {
+ htrdr_log_err(mats->htrdr,
+ "%s:%lu: error in word expension of the mrumtl path.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ wexp_is_allocated = 1;
+
+ if(wexp.we_wordc < 1) {
+ htrdr_log_err(mats->htrdr,
+ "%s:%lu: missing the MruMtl file for the material `%s'.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
+ str_cget(str));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Parse the mrumtl file if any */
+ if(strcmp(wexp.we_wordv[0], "none")) {
+ res = mrumtl_create(&mats->htrdr->logger, mats->htrdr->allocator,
+ mats->htrdr->verbose, &mtl.mrumtl);
+ if(res != RES_OK) {
+ htrdr_log_err(mats->htrdr,
+ "%s:%lu: error creating the MruMtl loader for the material `%s'-- %s.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
+ str_cget(str), res_to_cstr(res));
+ goto error;
+ }
+
+ res = mrumtl_load(mtl.mrumtl, wexp.we_wordv[0]);
+ if(res != RES_OK) goto error;
+ }
+
+ if(wexp.we_wordc < 2) {
+ if(mtl.mrumtl) {
+ htrdr_log_err(mats->htrdr,
+ "%s:%lu: missing temperature for the material `%s'.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
+ str_cget(str));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ } else {
+ /* Parse the temperature */
+ res = cstr_to_double(wexp.we_wordv[1], &mtl.temperature);
+ if(res != RES_OK) {
+ htrdr_log_err(mats->htrdr,
+ "%s:%lu: error parsing the temperature `%s' for the material `%s' "
+ "-- %s.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
+ wexp.we_wordv[1], str_cget(str), res_to_cstr(res));
+ goto error;
+ }
+ }
+
+ /* Register the material */
+ res = htable_name2mtl_set(&mats->name2mtl, str, &mtl);
+ if(res != RES_OK) {
+ htrdr_log_err(mats->htrdr,
+ "%s:%lu: could not register the material `%s' -- %s.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
+ str_cget(str), res_to_cstr(res));
+ goto error;
+ }
+
+ if(wexp.we_wordc > 2) {
+ htrdr_log_warn(mats->htrdr, "%s:%lu: unexpected text `%s'.\n",
+ txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
+ wexp.we_wordv[2]);
+ }
+
+exit:
+ if(wexp_is_allocated) wordfree(&wexp);
+ return res;
+error:
+ if(mtl.mrumtl) MRUMTL(ref_put(mtl.mrumtl));
+ goto exit;
+}
+
+static res_T
+parse_materials_list
+ (struct htrdr_materials* mats,
+ const char* filename,
+ const char* func_name)
+{
+ struct txtrdr* txtrdr = NULL;
+ struct str str;
+ res_T res = RES_OK;
+ ASSERT(mats && filename && func_name);
+
+ str_init(mats->htrdr->allocator, &str);
+
+ res = txtrdr_file(mats->htrdr->allocator, filename, '#', &txtrdr);
+ if(res != RES_OK) {
+ htrdr_log_err(mats->htrdr,
+ "%s: could not create the text reader for the material file `%s' -- %s.\n",
+ func_name, filename, res_to_cstr(res));
+ goto error;
+ }
+
+ for(;;) {
+ res = txtrdr_read_line(txtrdr);
+ if(res != RES_OK) {
+ htrdr_log_err(mats->htrdr,
+ "%s: error reading a line in the material file `%s' -- %s.\n",
+ func_name, filename, res_to_cstr(res));
+ goto error;
+ }
+
+ if(!txtrdr_get_cline(txtrdr)) break;
+
+ res = parse_material(mats, txtrdr, &str);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ str_release(&str);
+ if(txtrdr) txtrdr_ref_put(txtrdr);
+ return res;
+error:
+ goto exit;
+}
+
+static void
+materials_release(ref_T* ref)
+{
+ struct htable_name2mtl_iterator it, it_end;
+ struct htrdr_materials* mats;
+ ASSERT(ref);
+ mats = CONTAINER_OF(ref, struct htrdr_materials, ref);
+
+ htable_name2mtl_begin(&mats->name2mtl, &it);
+ htable_name2mtl_end(&mats->name2mtl, &it_end);
+ while(!htable_name2mtl_iterator_eq(&it, &it_end)) {
+ struct mtl* mtl = htable_name2mtl_iterator_data_get(&it);
+ /* The mrumtl can be NULL for semi transparent materials */
+ if(mtl->mrumtl) MRUMTL(ref_put(mtl->mrumtl));
+ htable_name2mtl_iterator_next(&it);
+ }
+ htable_name2mtl_release(&mats->name2mtl);
+ MEM_RM(mats->htrdr->allocator, mats);
+}
+
+static res_T
+create_bsdf_diffuse
+ (struct htrdr* htrdr,
+ const struct mrumtl_brdf* brdf,
+ const size_t ithread,
+ struct ssf_bsdf** out_bsdf)
+{
+ struct ssf_bsdf* bsdf = NULL;
+ double reflectivity = 0;
+ res_T res = RES_OK;
+ ASSERT(htrdr && brdf && out_bsdf);
+ ASSERT(mrumtl_brdf_get_type(brdf) == MRUMTL_BRDF_LAMBERTIAN);
+ ASSERT(ithread < htrdr->nthreads);
+
+ res = ssf_bsdf_create
+ (&htrdr->lifo_allocators[ithread], &ssf_lambertian_reflection, &bsdf);
+ if(res != RES_OK) goto error;
+
+ reflectivity = mrumtl_brdf_lambertian_get_reflectivity(brdf);
+ res = ssf_lambertian_reflection_setup(bsdf, reflectivity);
+ if(res != RES_OK) goto error;
+
+exit:
+ *out_bsdf = bsdf;
+ return res;
+error:
+ if(bsdf) { SSF(bsdf_ref_put(bsdf)); bsdf = NULL; }
+ goto exit;
+}
+
+static res_T
+create_bsdf_specular
+ (struct htrdr* htrdr,
+ const struct mrumtl_brdf* brdf,
+ const size_t ithread,
+ struct ssf_bsdf** out_bsdf)
+{
+ struct ssf_bsdf* bsdf = NULL;
+ struct ssf_fresnel* fresnel = NULL;
+ double reflectivity = 0;
+ res_T res = RES_OK;
+ ASSERT(htrdr && brdf && out_bsdf);
+ ASSERT(mrumtl_brdf_get_type(brdf) == MRUMTL_BRDF_SPECULAR);
+ ASSERT(ithread < htrdr->nthreads);
+
+ res = ssf_bsdf_create
+ (&htrdr->lifo_allocators[ithread], &ssf_specular_reflection, &bsdf);
+ if(res != RES_OK) goto error;
+
+ res = ssf_fresnel_create
+ (&htrdr->lifo_allocators[ithread], &ssf_fresnel_constant, &fresnel);
+ if(res != RES_OK) goto error;
+
+ reflectivity = mrumtl_brdf_specular_get_reflectivity(brdf);
+ res = ssf_fresnel_constant_setup(fresnel, reflectivity);
+ if(res != RES_OK) goto error;
+
+ res = ssf_specular_reflection_setup(bsdf, fresnel);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(fresnel) SSF(fresnel_ref_put(fresnel));
+ *out_bsdf = bsdf;
+ return res;
+error:
+ if(bsdf) { SSF(bsdf_ref_put(bsdf)); bsdf = NULL; }
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local symbol
+ ******************************************************************************/
+res_T
+htrdr_materials_create
+ (struct htrdr* htrdr,
+ const char* filename,
+ struct htrdr_materials** out_mtl)
+{
+ struct htrdr_materials* mats = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && filename && out_mtl);
+
+ mats = MEM_CALLOC(htrdr->allocator, 1, sizeof(*mats));
+ if(!mats) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: could not allocate the mats data structure -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ ref_init(&mats->ref);
+ mats->htrdr = htrdr;
+ htable_name2mtl_init(htrdr->allocator, &mats->name2mtl);
+
+ res = parse_materials_list(mats, filename, FUNC_NAME);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(out_mtl) *out_mtl = mats;
+ return res;
+error:
+ if(mats) {
+ htrdr_materials_ref_put(mats);
+ mats = NULL;
+ }
+ goto exit;
+}
+
+void
+htrdr_materials_ref_get(struct htrdr_materials* mats)
+{
+ ASSERT(mats);
+ ref_get(&mats->ref);
+}
+
+void
+htrdr_materials_ref_put(struct htrdr_materials* mats)
+{
+ ASSERT(mats);
+ ref_put(&mats->ref, materials_release);
+}
+
+int
+htrdr_materials_find_mtl
+ (struct htrdr_materials* mats,
+ const char* name,
+ struct htrdr_mtl* htrdr_mtl)
+{
+ struct str str;
+ struct htable_name2mtl_iterator it, it_end;
+ int found = 0;
+ ASSERT(mats && name && htrdr_mtl);
+
+ str_init(mats->htrdr->allocator, &str);
+ CHK(str_set(&str, name) == RES_OK);
+
+ htable_name2mtl_find_iterator(&mats->name2mtl, &str, &it);
+ htable_name2mtl_end(&mats->name2mtl, &it_end);
+ if(htable_name2mtl_iterator_eq(&it, &it_end)) { /* No material found */
+ *htrdr_mtl = HTRDR_MTL_NULL;
+ found = 0;
+ } else {
+ struct mtl* mtl = htable_name2mtl_iterator_data_get(&it);
+ ASSERT(mtl != NULL);
+ htrdr_mtl->name = str_cget(htable_name2mtl_iterator_key_get(&it));
+ htrdr_mtl->mrumtl = mtl->mrumtl;
+ htrdr_mtl->temperature = mtl->temperature;
+ found = 1;
+ }
+ str_release(&str);
+
+ return found;
+}
+
+res_T
+htrdr_mtl_create_bsdf
+ (struct htrdr* htrdr,
+ const struct htrdr_mtl* mtl,
+ const size_t ithread,
+ const double wavelength,
+ struct ssp_rng* rng,
+ struct ssf_bsdf** out_bsdf)
+{
+ struct ssf_bsdf* bsdf = NULL;
+ const struct mrumtl_brdf* brdf = NULL;
+ double r;
+ res_T res = RES_OK;
+ ASSERT(htrdr && mtl && wavelength && rng && out_bsdf);
+ ASSERT(ithread < htrdr->nthreads);
+
+ r = ssp_rng_canonical(rng);
+
+ res = mrumtl_fetch_brdf2(mtl->mrumtl, wavelength, r, &brdf);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr,
+ "%s: error retrieving the MruMtl BRDF for the wavelength %g.\n",
+ FUNC_NAME, wavelength);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ switch(mrumtl_brdf_get_type(brdf)) {
+ case MRUMTL_BRDF_LAMBERTIAN:
+ res = create_bsdf_diffuse(htrdr, brdf, ithread, &bsdf);
+ break;
+ case MRUMTL_BRDF_SPECULAR:
+ res = create_bsdf_specular(htrdr, brdf, ithread, &bsdf);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr, "%s: could not create the BSDF -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ *out_bsdf = bsdf;
+ return res;
+error:
+ if(bsdf) { SSF(bsdf_ref_put(bsdf)); bsdf = NULL; }
+ goto exit;
+}
+
diff --git a/src/htrdr_materials.h b/src/htrdr_materials.h
@@ -0,0 +1,67 @@
+/* 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_MATERIALS_H
+#define HTRDR_MATERIALS_H
+
+#include <rsys/rsys.h>
+
+/* Forward declarations */
+struct htrdr_materials;
+struct mrumtl;
+struct s3d_hit;
+struct ssf_bsdf;
+struct ssp_rng;
+
+struct htrdr_mtl {
+ const char* name;
+ const struct mrumtl* mrumtl;
+ double temperature;
+};
+static const struct htrdr_mtl HTRDR_MTL_NULL;
+
+extern LOCAL_SYM res_T
+htrdr_materials_create
+ (struct htrdr* htrdr,
+ const char* filename,
+ struct htrdr_materials** mats);
+
+extern LOCAL_SYM void
+htrdr_materials_ref_get
+ (struct htrdr_materials* mats);
+
+extern LOCAL_SYM void
+htrdr_materials_ref_put
+ (struct htrdr_materials* mats);
+
+/* Return 1 if the material exist and 0 otherwise */
+extern LOCAL_SYM int
+htrdr_materials_find_mtl
+ (struct htrdr_materials* mats,
+ const char* mtl_name,
+ struct htrdr_mtl* mtl);
+
+extern LOCAL_SYM res_T
+htrdr_mtl_create_bsdf
+ (struct htrdr* htrdr,
+ const struct htrdr_mtl* mtl,
+ const size_t ithread,
+ const double wavelength,
+ struct ssp_rng* rng,
+ struct ssf_bsdf** bsdf);
+
+#endif /* HTRDR_MATERIALS_H */
+
diff --git a/src/htrdr_mtl.c b/src/htrdr_mtl.c
@@ -1,281 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#define _POSIX_C_SOURCE 200112L /* strtok_r and wordexp support */
-
-#include "htrdr.h"
-#include "htrdr_mtl.h"
-
-#include <modradurb/mrumtl.h>
-
-#include <rsys/cstr.h>
-#include <rsys/hash_table.h>
-#include <rsys/mem_allocator.h>
-#include <rsys/ref_count.h>
-#include <rsys/str.h>
-#include <rsys/text_reader.h>
-
-#include <string.h>
-#include <wordexp.h>
-
-/* Generate the hash table that maps a material name to its data */
-#define HTABLE_NAME name2mtl
-#define HTABLE_DATA struct mrumtl*
-#define HTABLE_KEY struct str
-#define HTABLE_KEY_FUNCTOR_INIT str_init
-#define HTABLE_KEY_FUNCTOR_RELEASE str_release
-#define HTABLE_KEY_FUNCTOR_COPY str_copy
-#define HTABLE_KEY_FUNCTOR_COPY_AND_RELEASE str_copy_and_release
-#define HTABLE_KEY_FUNCTOR_HASH str_hash
-#define HTABLE_KEY_FUNCTOR_EQ str_eq
-#include <rsys/hash_table.h>
-
-struct htrdr_mtl {
- struct htable_name2mtl name2mtl;
- struct htrdr* htrdr;
- ref_T ref;
-};
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-static res_T
-parse_material
- (struct htrdr_mtl* mtl,
- struct txtrdr* txtrdr,
- struct str* str) /* Scratch string */
-{
- wordexp_t wexp;
- char* tk = NULL;
- char* tk_ctx = NULL;
- struct mrumtl* mrumtl = NULL;
- int err = 0;
- int wexp_is_allocated = 0;
- res_T res = RES_OK;
- ASSERT(mtl && txtrdr);
-
- tk = strtok_r(txtrdr_get_line(txtrdr), " \t", &tk_ctx);
- ASSERT(tk);
-
- res = str_set(str, tk);
- if(res != RES_OK) {
- htrdr_log_err(mtl->htrdr,
- "%s:%lu: could not copy the material name `%s' -- %s.\n",
- txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr), tk,
- res_to_cstr(res));
- goto error;
- }
-
- tk = strtok_r(NULL, "", &tk_ctx);
- if(!tk) {
- htrdr_log_err(mtl->htrdr,
- "%s:%lu: missing the MruMtl file for the material `%s'.\n",
- txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
- str_cget(str));
- res = RES_BAD_ARG;
- goto error;
- }
-
- err = wordexp(tk, &wexp, 0);
- if(err) {
- htrdr_log_err(mtl->htrdr,
- "%s:%lu: error in word expension of the mrumtl path.\n",
- txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr));
- res = RES_BAD_ARG;
- goto error;
- }
- wexp_is_allocated = 1;
-
- if(wexp.we_wordc < 1) {
- htrdr_log_err(mtl->htrdr,
- "%s:%lu: missing the MruMtl file for the material `%s'.\n",
- txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
- str_cget(str));
- res = RES_BAD_ARG;
- goto error;
- }
-
- res = mrumtl_create
- (&mtl->htrdr->logger, mtl->htrdr->allocator, mtl->htrdr->verbose, &mrumtl);
- if(res != RES_OK) {
- htrdr_log_err(mtl->htrdr,
- "%s:%lu: error creating the MruMtl loader for the material `%s'-- %s.\n",
- txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
- str_cget(str), res_to_cstr(res));
- goto error;
- }
-
- res = mrumtl_load(mrumtl, wexp.we_wordv[0]);
- if(res != RES_OK) goto error;
-
- /* Register the material */
- res = htable_name2mtl_set(&mtl->name2mtl, str, &mrumtl);
- if(res != RES_OK) {
- htrdr_log_err(mtl->htrdr,
- "%s:%lu: could not register the material `%s' -- %s.\n",
- txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
- str_cget(str), res_to_cstr(res));
- goto error;
- }
-
- if(wexp.we_wordc > 1) {
- htrdr_log_warn(mtl->htrdr, "%s:%lu: unexpected text `%s'.\n",
- txtrdr_get_name(txtrdr), (unsigned long)txtrdr_get_line_num(txtrdr),
- wexp.we_wordv[1]);
- }
-
-exit:
- if(wexp_is_allocated) wordfree(&wexp);
- return res;
-error:
- if(mrumtl) MRUMTL(ref_put(mrumtl));
- goto exit;
-}
-
-static res_T
-parse_materials_list
- (struct htrdr_mtl* mtl,
- const char* filename,
- const char* func_name)
-{
- struct txtrdr* txtrdr = NULL;
- struct str str;
- res_T res = RES_OK;
- ASSERT(mtl && filename && func_name);
-
- str_init(mtl->htrdr->allocator, &str);
-
- res = txtrdr_file(mtl->htrdr->allocator, filename, '#', &txtrdr);
- if(res != RES_OK) {
- htrdr_log_err(mtl->htrdr,
- "%s: could not create the text reader for the material file `%s' -- %s.\n",
- func_name, filename, res_to_cstr(res));
- goto error;
- }
-
- for(;;) {
- res = txtrdr_read_line(txtrdr);
- if(res != RES_OK) {
- htrdr_log_err(mtl->htrdr,
- "%s: error reading a line in the material file `%s' -- %s.\n",
- func_name, filename, res_to_cstr(res));
- goto error;
- }
-
- if(!txtrdr_get_cline(txtrdr)) break;
-
- res = parse_material(mtl, txtrdr, &str);
- if(res != RES_OK) goto error;
- }
-
-exit:
- str_release(&str);
- if(txtrdr) txtrdr_ref_put(txtrdr);
- return res;
-error:
- goto exit;
-}
-
-static void
-mtl_release(ref_T* ref)
-{
- struct htable_name2mtl_iterator it, it_end;
- struct htrdr_mtl* mtl;
- ASSERT(ref);
- mtl = CONTAINER_OF(ref, struct htrdr_mtl, ref);
-
- htable_name2mtl_begin(&mtl->name2mtl, &it);
- htable_name2mtl_end(&mtl->name2mtl, &it_end);
- while(!htable_name2mtl_iterator_eq(&it, &it_end)) {
- struct mrumtl* mrumtl = *htable_name2mtl_iterator_data_get(&it);
- MRUMTL(ref_put(mrumtl));
- htable_name2mtl_iterator_next(&it);
- }
- htable_name2mtl_release(&mtl->name2mtl);
- MEM_RM(mtl->htrdr->allocator, mtl);
-}
-
-/*******************************************************************************
- * Local symbol
- ******************************************************************************/
-res_T
-htrdr_mtl_create
- (struct htrdr* htrdr,
- const char* filename,
- struct htrdr_mtl** out_mtl)
-{
- struct htrdr_mtl* mtl = NULL;
- res_T res = RES_OK;
- ASSERT(htrdr && filename && out_mtl);
-
- mtl = MEM_CALLOC(htrdr->allocator, 1, sizeof(*mtl));
- if(!mtl) {
- res = RES_MEM_ERR;
- htrdr_log_err(htrdr,
- "%s: could not allocate the mtl data structure -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
- }
- ref_init(&mtl->ref);
- mtl->htrdr = htrdr;
- htable_name2mtl_init(htrdr->allocator, &mtl->name2mtl);
-
- res = parse_materials_list(mtl, filename, FUNC_NAME);
- if(res != RES_OK) goto error;
-
-exit:
- if(out_mtl) *out_mtl = mtl;
- return res;
-error:
- if(mtl) {
- htrdr_mtl_ref_put(mtl);
- mtl = NULL;
- }
- goto exit;
-}
-
-void
-htrdr_mtl_ref_get(struct htrdr_mtl* mtl)
-{
- ASSERT(mtl);
- ref_get(&mtl->ref);
-}
-
-void
-htrdr_mtl_ref_put(struct htrdr_mtl* mtl)
-{
- ASSERT(mtl);
- ref_put(&mtl->ref, mtl_release);
-}
-
-const struct mrumtl*
-htrdr_mtl_get(struct htrdr_mtl* mtl, const char* name)
-{
- struct str str;
- struct mrumtl** pmrumtl = NULL;
- struct mrumtl* mrumtl = NULL;
- ASSERT(mtl && name);
-
- str_init(mtl->htrdr->allocator, &str);
- CHK(str_set(&str, name) == RES_OK);
-
- pmrumtl = htable_name2mtl_find(&mtl->name2mtl, &str);
- if(pmrumtl) mrumtl = *pmrumtl;
-
- str_release(&str);
- return mrumtl;
-}
-
diff --git a/src/htrdr_mtl.h b/src/htrdr_mtl.h
@@ -1,44 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTRDR_MTL_H
-#define HTRDR_MTL_H
-
-struct htrdr_mtl;
-struct mrumtl;
-
-extern LOCAL_SYM res_T
-htrdr_mtl_create
- (struct htrdr* htrdr,
- const char* filename,
- struct htrdr_mtl** mtl);
-
-extern LOCAL_SYM void
-htrdr_mtl_ref_get
- (struct htrdr_mtl* mtl);
-
-extern LOCAL_SYM void
-htrdr_mtl_ref_put
- (struct htrdr_mtl* mtl);
-
-/* Return NULL if the material name does not exist */
-extern const struct mrumtl*
-htrdr_mtl_get
- (struct htrdr_mtl* mtl,
- const char* mtl_name);
-
-#endif /* HTRDR_MTL_H */
-
diff --git a/src/htrdr_rectangle.c b/src/htrdr_rectangle.c
@@ -26,6 +26,7 @@ struct htrdr_rectangle {
double axis_x[3];
double axis_y[3];
+ double normal[3];
double position[3]; /* Center of the rectangle */
struct htrdr* htrdr;
ref_T ref;
@@ -49,10 +50,10 @@ rectangle_release(ref_T* ref)
res_T
htrdr_rectangle_create
(struct htrdr* htrdr,
+ const double sz[2],
const double pos[3],
const double tgt[3],
- const double up[2],
- const double sz[2],
+ const double up[3],
struct htrdr_rectangle** out_rect)
{
struct htrdr_rectangle* rect = NULL;
@@ -91,6 +92,7 @@ htrdr_rectangle_create
d3_muld(rect->axis_x, x, sz[0]*0.5);
d3_muld(rect->axis_y, y, sz[1]*0.5);
+ d3_set(rect->normal, z);
d3_set(rect->position, pos);
exit:
@@ -131,3 +133,10 @@ htrdr_rectangle_ref_put(struct htrdr_rectangle* rect)
ref_put(&rect->ref, rectangle_release);
}
+void
+htrdr_rectangle_get_normal(const struct htrdr_rectangle* rect, double normal[3])
+{
+ ASSERT(rect && normal);
+ d3_set(normal, rect->normal);
+}
+
diff --git a/src/htrdr_rectangle.h b/src/htrdr_rectangle.h
@@ -28,7 +28,7 @@ htrdr_rectangle_create
const double sz[2], /* Size of the rectangle along its local X and Y axis */
const double pos[3], /* World space position of the rectangle center */
const double tgt[3], /* Vector orthognal to the rectangle plane */
- const double up[2], /* vector orthogonal to the rectangle X axis */
+ const double up[3], /* vector orthogonal to the rectangle X axis */
struct htrdr_rectangle** rect);
extern LOCAL_SYM void
@@ -45,5 +45,10 @@ htrdr_rectangle_sample_pos
const double sample[2], /* In [0, 1[ */
double pos[3]);
+extern LOCAL_SYM void
+htrdr_rectangle_get_normal
+ (const struct htrdr_rectangle* rect,
+ double normal[3]);
+
#endif /* HTRDR_RECTANGLE_H */
diff --git a/src/htrdr_sensor.c b/src/htrdr_sensor.c
@@ -0,0 +1,136 @@
+/* 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_camera.h"
+#include "htrdr_ground.h"
+#include "htrdr_interface.h"
+#include "htrdr_rectangle.h"
+#include "htrdr_sensor.h"
+
+#include <star/s3d.h>
+#include <star/ssp.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+sample_camera_ray
+ (struct htrdr_camera* cam,
+ const size_t ipix[2],
+ const double pix_sz[2],
+ struct ssp_rng* rng,
+ double ray_org[3],
+ double ray_dir[3])
+{
+ double pix_samp[2];
+ ASSERT(cam && ipix && pix_sz && rng && ray_org && ray_dir);
+
+ /* Sample a position into the pixel, in the normalized image plane */
+ pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
+ pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
+
+ /* Generate a ray starting from the pinhole camera and passing through the
+ * pixel sample */
+ htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir);
+
+ return RES_OK;
+}
+
+static res_T
+sample_rectangle_ray
+ (struct htrdr_rectangle* rect,
+ struct htrdr* htrdr,
+ const size_t ipix[2],
+ const double pix_sz[2],
+ struct ssp_rng* rng,
+ double ray_org[3],
+ double ray_dir[3])
+{
+ struct s3d_hit hit = S3D_HIT_NULL;
+ double pix_samp[2];
+ const double up_dir[3] = {0,0,1};
+ const double range[2] = {0, INF};
+ double normal[3];
+ ASSERT(rect && htrdr && ipix && pix_sz && rng && ray_org && ray_dir);
+
+ /* Sample a position into the pixel, in the normalized image plane */
+ pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
+ pix_samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
+
+ /* Retrieve the world space position of pix_samp */
+ htrdr_rectangle_sample_pos(rect, pix_samp, ray_org);
+
+ /* Check that `ray_org' is not included into a geometry */
+ HTRDR(ground_trace_ray(htrdr->ground, ray_org, up_dir, range, NULL, &hit));
+
+ /* Up direction is occluded. Check if the sample must be rejected, i.e. does it
+ * lies inside a geometry? */
+ if(!S3D_HIT_NONE(&hit)) {
+ struct htrdr_interface interf = HTRDR_INTERFACE_NULL;
+ const struct htrdr_mtl* mtl = NULL;
+ float N[3]; /* Normalized normal of the hit */
+ float wi[3];
+ float cos_wi_N;
+
+ /* Compute the cosine between the up direction and the hit normal */
+ f3_set_d3(wi, up_dir);
+ f3_normalize(N, hit.normal);
+ cos_wi_N = f3_dot(wi, N);
+
+ /* Fetch the hit interface and retrieve the material into which the ray was
+ * traced */
+ htrdr_ground_get_interface(htrdr->ground, &hit, &interf);
+ mtl = cos_wi_N < 0 ? &interf.mtl_front : &interf.mtl_back;
+
+ /* Reject the sample if the incident direction do not travel into the sky */
+ if(strcmp(mtl->name, htrdr->sky_mtl_name) != 0) return RES_BAD_OP;
+ }
+
+ /* Sample a ray direction */
+ htrdr_rectangle_get_normal(rect, normal);
+ ssp_ran_hemisphere_cos(rng, normal, ray_dir, NULL);
+
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_sensor_sample_primary_ray
+ (const struct htrdr_sensor* sensor,
+ struct htrdr* htrdr,
+ const size_t ipix[2],
+ const double pix_sz[2],
+ struct ssp_rng* rng,
+ double ray_org[3],
+ double ray_dir[3])
+{
+ res_T res = RES_OK;
+ switch(sensor->type) {
+ case HTRDR_SENSOR_CAMERA:
+ res = sample_camera_ray(sensor->camera, ipix, pix_sz, rng, ray_org, ray_dir);
+ break;
+ case HTRDR_SENSOR_RECTANGLE:
+ res = sample_rectangle_ray(sensor->rectangle, htrdr, ipix,
+ pix_sz, rng, ray_org, ray_dir);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return res;
+}
+
diff --git a/src/htrdr_sensor.h b/src/htrdr_sensor.h
@@ -0,0 +1,48 @@
+/* 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_SENSOR_H
+#define HTRDR_SENSOR_H
+
+#include <rsys/rsys.h>
+
+/* Forward declarations */
+struct htrdr;
+struct ssp_rng;
+
+enum htrdr_sensor_type {
+ HTRDR_SENSOR_CAMERA,
+ HTRDR_SENSOR_RECTANGLE
+};
+
+struct htrdr_sensor {
+ struct htrdr_camera* camera;
+ struct htrdr_rectangle* rectangle;
+ enum htrdr_sensor_type type;
+};
+
+extern LOCAL_SYM res_T
+htrdr_sensor_sample_primary_ray
+ (const struct htrdr_sensor* sensor,
+ struct htrdr* htrdr,
+ const size_t ipix[2],
+ const double pix_sz[2],
+ struct ssp_rng* rng,
+ double ray_org[3],
+ double ray_dir[3]);
+
+#endif /* HTRDR_SENSOR_H */
+
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -19,6 +19,13 @@
#include <rsys/rsys.h>
+/* Define the radiance component */
+enum htrdr_radiance_cpnt_flag {
+ HTRDR_RADIANCE_DIRECT = BIT(0),
+ HTRDR_RADIANCE_DIFFUSE = BIT(1),
+ HTRDR_RADIANCE_ALL = HTRDR_RADIANCE_DIRECT | HTRDR_RADIANCE_DIFFUSE
+};
+
/* Monte carlo accumulator */
struct htrdr_accum {
double sum_weights; /* Sum of Monte-Carlo weights */
@@ -42,6 +49,12 @@ struct htrdr_pixel_xwave {
};
static const struct htrdr_pixel_xwave HTRDR_PIXEL_XWAVE_NULL;
+struct htrdr_pixel_flux {
+ struct htrdr_accum flux;
+ struct htrdr_accum time;
+};
+static const struct htrdr_pixel_flux HTRDR_PIXEL_FLUX;
+
struct htrdr_pixel_image {
struct htrdr_estimate X; /* In W/m^2/sr */
struct htrdr_estimate Y; /* In W/m^2/sr */
@@ -62,6 +75,7 @@ htrdr_compute_radiance_sw
(struct htrdr* htrdr,
const size_t ithread,
struct ssp_rng* rng,
+ const int cpnt_mask, /* Combination of enum htrdr_radiance_cpnt_flag */
const double pos[3],
const double dir[3],
const double wlen, /* In nanometer */
@@ -81,9 +95,9 @@ htrdr_compute_radiance_lw
const size_t iquad); /* Index of the quadrature point into the band */
extern LOCAL_SYM res_T
-htrdr_draw_radiance
+htrdr_draw_map
(struct htrdr* htrdr,
- const struct htrdr_camera* cam,
+ const struct htrdr_sensor* sensor,
const size_t width, /* Image width */
const size_t height, /* Image height */
const size_t spp, /* #samples per pixel, i.e. #realisations */
@@ -121,35 +135,49 @@ htrdr_accum_get_estimation
}
static INLINE size_t
-htrdr_spectral_type_get_pixsz(const enum htrdr_spectral_type type)
+htrdr_spectral_type_get_pixsz
+ (const enum htrdr_spectral_type spectral_type,
+ const enum htrdr_sensor_type sensor_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;
+ if(sensor_type == HTRDR_SENSOR_RECTANGLE) {
+ sz = sizeof(struct htrdr_pixel_flux);
+ } else {
+ ASSERT(sensor_type == HTRDR_SENSOR_CAMERA);
+ switch(spectral_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)
+htrdr_spectral_type_get_pixal
+ (const enum htrdr_spectral_type spectral_type,
+ const enum htrdr_sensor_type sensor_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;
+ if(sensor_type == HTRDR_SENSOR_RECTANGLE) {
+ al = ALIGNOF(struct htrdr_pixel_flux);
+ } else {
+ ASSERT(sensor_type == HTRDR_SENSOR_CAMERA);
+ switch(spectral_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;
}
diff --git a/src/htrdr_sun.c b/src/htrdr_sun.c
@@ -105,7 +105,7 @@ htrdr_sun_set_direction(struct htrdr_sun* sun, const double dir[3])
d33_basis(sun->frame, dir);
}
-double*
+double
htrdr_sun_sample_direction
(struct htrdr_sun* sun,
struct ssp_rng* rng,
@@ -113,7 +113,8 @@ htrdr_sun_sample_direction
{
ASSERT(sun && rng && dir);
ssp_ran_sphere_cap_uniform_local(rng, sun->cos_half_angle, dir, NULL);
- return d33_muld3(dir, sun->frame, dir);
+ d33_muld3(dir, sun->frame, dir);
+ return 1.0 / htrdr_sun_get_solid_angle(sun);
}
double
diff --git a/src/htrdr_sun.h b/src/htrdr_sun.h
@@ -43,8 +43,8 @@ htrdr_sun_set_direction
(struct htrdr_sun* sun,
const double direction[3]); /* Must be normalized */
-/* Return a direction that points *toward* the sun */
-extern LOCAL_SYM double*
+/* Return a pdf of the sampled dir */
+extern LOCAL_SYM double
htrdr_sun_sample_direction
(struct htrdr_sun* sun,
struct ssp_rng* rng,
