commit 25782a5af46e920d45cdf4d2b7904e967b3053e4
parent 0064745d8e23e3a31bacde3aaac258adcf09397c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 4 May 2020 10:19:29 +0200
Merge branch 'feature_mrumtl' into develop
Diffstat:
27 files changed, 3861 insertions(+), 1230 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -24,14 +24,15 @@ option(NO_TEST "Do not build the tests" OFF)
################################################################################
# Check dependencies
################################################################################
+find_package(AW 2.0 REQUIRED)
+find_package(HTSky 0.1 REQUIRED)
+find_package(MruMtl REQUIRED)
find_package(RCMake 0.3 REQUIRED)
find_package(RSys 0.7 REQUIRED)
find_package(Star3D 0.6 REQUIRED)
-find_package(Star3DAW 0.3.1 REQUIRED)
find_package(StarSF 0.6 REQUIRED)
find_package(StarSP 0.8 REQUIRED)
find_package(StarVX REQUIRED)
-find_package(HTSky REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
find_package(MPI 1 REQUIRED)
@@ -42,9 +43,10 @@ include(rcmake_runtime)
set(CMAKE_C_COMPILER ${MPI_C_COMPILER})
include_directories(
+ ${AW_INCLUDE_DIR}
+ ${MruMtl_INCLUDE_DIR}
${RSys_INCLUDE_DIR}
${Star3D_INCLUDE_DIR}
- ${Star3DAW_INCLUDE_DIR}
${StarSF_INCLUDE_DIR}
${StarSP_INCLUDE_DIR}
${StarVX_INCLUDE_DIR}
@@ -65,8 +67,6 @@ 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_GROUND_BSDF "HTRDR_BSDF_DIFFUSE")
-set(HTRDR_ARGS_DEFAULT_GROUND_REFLECTIVITY "0.5")
set(HTRDR_ARGS_DEFAULT_IMG_WIDTH "320")
set(HTRDR_ARGS_DEFAULT_IMG_HEIGHT "240")
set(HTRDR_ARGS_DEFAULT_IMG_SPP "1")
@@ -91,11 +91,16 @@ set(HTRDR_FILES_SRC
htrdr_args.c
htrdr_buffer.c
htrdr_camera.c
+ htrdr_cie_xyz.c
htrdr_compute_radiance_sw.c
- htrdr_draw_radiance_sw.c
+ htrdr_compute_radiance_lw.c
+ htrdr_draw_radiance.c
htrdr_grid.c
htrdr_ground.c
+ htrdr_interface.c
htrdr_main.c
+ htrdr_mtl.c
+ htrdr_ran_lw.c
htrdr_rectangle.c
htrdr_slab.c
htrdr_sun.c)
@@ -104,8 +109,12 @@ set(HTRDR_FILES_INC
htrdr_c.h
htrdr_buffer.h
htrdr_camera.h
+ htrdr_cie_xyz.h
htrdr_grid.h
htrdr_ground.h
+ htrdr_interface.h
+ htrdr_mtl.h
+ htrdr_ran_lw.h
htrdr_rectangle.h
htrdr_slab.h
htrdr_sun.h
@@ -122,7 +131,7 @@ rcmake_prepend_path(HTRDR_FILES_INC2 ${CMAKE_CURRENT_BINARY_DIR})
rcmake_prepend_path(HTRDR_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_executable(htrdr ${HTRDR_FILES_SRC} ${HTRDR_FILES_INC} ${HTRDR_FILES_INC2})
-target_link_libraries(htrdr HTSky RSys Star3D Star3DAW StarSF StarSP)
+target_link_libraries(htrdr AW HTSky MruMtl RSys Star3D StarSF StarSP)
if(CMAKE_COMPILER_IS_GNUCC)
target_link_libraries(htrdr m)
diff --git a/cmake/doc/CMakeLists.txt b/cmake/doc/CMakeLists.txt
@@ -29,7 +29,7 @@ endif()
################################################################################
# Copy doc files
################################################################################
-set(MAN_NAMES htrdr-image.5)
+set(MAN_NAMES htrdr-image.5 htrdr-materials.5 htrdr-obj.5)
set(MAN_FILES)
foreach(_name IN LISTS MAN_NAMES)
diff --git a/doc/htrdr-image.5.txt b/doc/htrdr-image.5.txt
@@ -28,14 +28,26 @@ The *htrdr-image* is a raw image file format where data are stored in plain
text. Characters after the '#' character are considered as comments and are
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 4 pairs of floating points data
-representing the pixel color encoded in the CIE 1931 XYZ color space and the
-time spent per realisation to estimate the pixel. The first, second and third
-pair encodes the estimated radiance in W.sr^-1.m^-2 of the X, Y and Z pixel
-components, respectively. The first value of each pair is the expected value of
-the estimated radiance while the second one is its associated standard
-deviation. The fourth pair saves the estimate in microseconds of the per
-radiative path computation time and its standard error.
+line and per column. Then each line stores 8 pixel components.
+
+If the image is a regular rendering in the visible part of the spectrum, the
+pixel components are actually 4 pairs of floating points data representing the
+pixel color encoded in the CIE 1931 XYZ color space and the per radiative path
+computation time. The first, second and third pair encode the estimated
+integrated radiance in W.sr^-1.m^-2 of the X, Y and Z pixel components,
+respectively. The first value of each pair is the expected value of the
+estimated radiance while the second one is its associated standard deviation.
+The fourth pair saves the estimate in microseconds of the per radiative path
+computation time and its standard error.
+
+If the image is an infrared rendering, the first and second pixel components
+store the expected value and the standard error, respectively, of the
+estimated brightness temperature in Kelvin. The third and fourth component
+save the expected value and standard deviation of the pixel radiance in
+W.sr^-1.m^-2.nm^-1. The fifth and sixth pixel components are unused. Finally
+the last 2 pixel components save, as for a regular rendering, the estimate in
+microseconds of the per radiative path computation time and its standard
+error.
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
@@ -43,7 +55,7 @@ pixels per line and M the overall number of lines in the image, the first N
pixels correspond to the pixels of the top line of the image, the following N
pixels are the pixels of the second line and so on.
-Note that the *htpp*(1) program can be used to convert an *htrdr-image* in a
+Note that the *htpp*(1) program can be used to convert an *htrdr-image* into a
regular PPM image [1].
GRAMMAR
@@ -59,11 +71,18 @@ GRAMMAR
<width> ::= INTEGER
<height> ::= INTEGER
-<pixel> ::= <X> <Y> <Z> <time>
+<pixel> ::= <pixel-sw>
+ | <pixel-lw>
+
+<pixel-sw> ::= <X> <Y> <Z> <time>
+<pixel-lw> ::= <temperature> <radiance> 0 0 <time>
+
<X> ::= <estimate>
<Y> ::= <estimate>
<Z> ::= <estimate>
<time> ::= <estimate>
+<temperature> ::= <estimate>
+<radiance> ::= <estimate>
<estimate> ::= <expected-value> <standard-error>
<expected-value> ::= REAL
@@ -72,12 +91,13 @@ GRAMMAR
EXAMPLE
-------
-The following output is emitted by *htrdr*(1) invoked to render an image of
-*800* by *600* pixels. Note that actually the comments or the blank lines are
-not necessarily written as it by *htrdr*(1); they are used here only to help
-in understanding the data layout. The comment after each pixel gives the two
-dimensional index of the pixel in the image: the first and second integer is
-the index of the line and the column of the pixel in the image, respectively.
+The following output was produced by *htrdr*(1) invoked to render an image of
+*800* by *600* pixels. Note that actually the comments and blank lines were
+not necessarily written by *htrdr*(1); they are used here only to help the
+reader understand the data layout. The comment after each pixel gives the
+two-dimensional index of the pixel in the image: the first and second integer
+is the index of the line and the column of the pixel in the image,
+respectively.
[verse]
------
diff --git a/doc/htrdr-materials.5.txt b/doc/htrdr-materials.5.txt
@@ -0,0 +1,66 @@
+// 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/>.
+:toc:
+
+htrdr-materials(5)
+=================
+
+NAME
+----
+htrdr-materials - list of materials used by the ground geometry in htrdr(1)
+
+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.
+
+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
+and tabulations, are simply skipped.
+
+GRAMMAR
+-------
+
+[verse]
+-------
+<htrdr-materials> ::= <material>
+ [ <material> ... ]
+<material> ::= <name> <mrumtl-path>
+<name> ::= STRING
+<mrumtl-path> ::= PATH
+-------
+
+EXAMPLE
+-------
+
+The following file lists 2 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.
+
+[verse]
+-------
+grass /opt/materials/A001.mrumtl
+sand /opt/materials/B002.mrumtl
+-------
+
+SEE ALSO
+--------
+*htrdr*(1), *htrdr-obj*(5), *mrumtl*(5)
diff --git a/doc/htrdr-obj.5.txt b/doc/htrdr-obj.5.txt
@@ -0,0 +1,103 @@
+// 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/>.
+:toc:
+
+htrdr-obj(5)
+============
+
+NAME
+----
+htrdr-obj - file format of the ground geometry in htrdr(1)
+
+DESCRIPTION
+-----------
+A *htrdr-obj* file is a regular OBJ [1] composed only of triangular meshes.
+Each triangle must be included in a material group as defined by the 'usemtl'
+directive. The name of the material group must be of the form
+"<__front-mtl-name__>:<__back-mtl-name__>", where <__front-mtl-name__> and
+<__back-mtl-name__> are strings separated by a colon (:) defining the name of
+the front and back facing materials, respectively. These names can be composed
+of any characters expected for spaces and tabulations. Note that regarding the
+*htrdr*(1) convention, a triangle side is said "front facing" when its
+vertices are clock-wise ordered.
+
+Note that to be a valid *htrdr-obj*(5) file for *htrdr*(1), the front and the
+back facing names must reference a material listed in *htrdr-materials*(5)
+file submitted to the *htrdr*(1) command line.
+
+The grammar of a *htrdr-obj* file is thus a subset of the OBJ file
+format [1] with only a specific convention regarding the material name.
+As a consequence, any software supporting the OBJ file format can be
+used to create or visualise an *htrdr-obj* file.
+
+EXAMPLES
+--------
+Define a quad at the interface between the air medium and the concrete
+material.
+
+[verse]
+-------
+v -5.0 -5.0 0
+v -5.0 5.0 0
+v 5.0 -5.0 0
+v 5.0 5.0 0
+
+usemtl air:concrete
+f 1 2 3
+f 3 2 4
+-------
+
+Define a wooden cube whose faces Z-aligned faces are against a brick material.
+The remaining faces are in contact with the air.
+
+[verse]
+-------
+v 0 0 0
+v 1 0 0
+v 0 1 0
+v 1 1 0
+v 0 0 1
+v 1 0 1
+v 0 1 1
+v 1 1 1
+
+usemtl wood:air
+f 1 3 2
+f 2 3 4
+f 1 5 3
+f 3 5 7
+f 5 6 7
+f 7 6 8
+f 4 8 2
+f 2 8 6
+
+usemtl wood:brick
+f 3 7 4
+f 4 7 8
+f 1 2 5
+f 5 2 6
+-------
+
+NOTES
+-----
+
+1. OBJ file format -
+ <http://www.martinreddy.net/gfx/3d/OBJ.spec>
+
+SEE ALSO
+--------
+
+*htrdr*(1)
diff --git a/doc/htrdr.1.txt.in b/doc/htrdr.1.txt.in
@@ -29,43 +29,49 @@ SYNOPSIS
DESCRIPTION
-----------
-*htrdr* is an image renderer in the visible part of the spectrum, for scenes
-composed of an atmospheric gas mixture, clouds, and a ground. 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 OBJ [2]
-representing the ground geometry with a diffuse reflectivity (*-g* _ground_).
-Both, the clouds and the ground, can be infinitely repeated along the X and Y
-axis by setting the *-r* and the *-R* options, respectively.
-
-*htrdr* evaluates the intensity incoming on each pixel of the sensor
-array. The underlying algorithm is based on a Monte-Carlo method: it consists
-in simulating a given number of optical paths originating from the camera,
+*htrdr* is an image renderer of scenes composed of an atmospheric gas mixture,
+clouds, and a ground. Images can be rendered in the visible or the infrared
+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.
+
+*htrdr* evaluates the intensity incoming on each pixel of the sensor array. The
+underlying algorithm is based on a Monte-Carlo method: it consists in
+simulating a given number of optical paths originating from the camera,
directed into the atmosphere, taking into account light absorption and
-scattering phenomena. The computation is performed over the whole visible part
-of the spectrum, 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 required
-observation position.
-
-In *htrdr* the spatial unit 1.0 corresponds to one meter. The estimated
-radiance of each pixel component is given in W.sr^-1.m^-2. The pixels are
-written into the _output_ file or to the standard output whether the *-o*
-option is defined or not, respectively. 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, the estimated radiance of a pixel
-component is provided with its numerical accuracy.
-
-During simulation, *htrdr* dynamically loads/unloads cloud properties to handle
-clouds whose data do not feat in main memory. *htrdr* also supports shared
-memory parallelism and relies on the Message Passing Interface specification
-[4] to parallelise its computations in a distribute memory environment; it can
-thus be run either directly or through a MPI process launcher like *mpirun*(1).
+scattering phenomena. When rendering in the visible part of the spectrum, the
+computation is performed over the visible part of the spectrum in [380, 780]
+nanometers, for the three components of the CIE 1931 XYZ colorimetric space
+that are subsequently recombined in order to obtain the final color for each
+pixel, and finally the whole image of the scene as seen from the set
+observation position. In long waves, the rendering is performed for the range
+of wavelengths defined by the *-l* option. The estimated radiance per pixel is
+then converted to its brightness temperature and both are saved in the output
+image.
+
+In *htrdr* the spatial unit 1.0 corresponds to one meter, the estimated
+radiance is given in W.sr^-1.m^-2 and the temperatures are expressed in Kelvin.
+The results are written to the output file if the *-o* option is defined and
+the standard output otherwise. The output image is a list of raw ASCII data
+formatted with respect to the *htrdr-image*(5) file format. Since *htrdr*
+relies on the Monte-Carlo method, each estimation is given with its numerical
+accuracy.
+
+During the simulation, *htrdr* dynamically loads/unloads cloud properties to
+handle clouds whose data that do not feat in main memory. *htrdr* also supports
+shared memory parallelism and relies on the Message Passing Interface
+specification [3] to parallelise its computations in a distributed memory
+environment; it can thus be run either directly or through a MPI process
+launcher like *mpirun*(1).
OPTIONS
-------
@@ -73,13 +79,6 @@ OPTIONS
Path toward the file containing the gas optical properties of the atmosphere.
Data must be formatted according to the fileformat described in [1].
-*-b* <diffuse|specular>::
- Define how light is reflected at the ground surface. When set to diffuse, the
- reflections are lambertian while a specular surface simulates a perfect
- mirror. By default, the ground surface is diffuse. Note that this option
- controls only the directional part of the surface reflections; the amount of
- reflected energy is defined by the *-e* option.
-
*-c* _clouds_::
Submit a *htcp*(5) file describing the properties of the clouds. If not
defined, only the atmosphere properties submitted through the *-a* option
@@ -117,22 +116,16 @@ OPTIONS
*-d*::
Write in _output_ the space partitioning data structures used to speed up
the radiative transfer computations in the clouds. The written data are
- octrees saved in the VTK file format [3]. Each octree node stores the minimum
+ octrees saved in the VTK file format [2]. Each octree node stores the minimum
and the maximum of the extinction coefficients of the cloud cells overlapped
- by the octree node. In the _output_ file, each octree is separated of the
+ by the octree node. In the _output_ file, each octree is separated from the
previous one by a line with three minus characters, i.e. '---'.
-*-e* _reflectivity_::
- Reflectivity of the ground geometry in [0, 1]. By default it is set to
- @HTRDR_ARGS_DEFAULT_GROUND_REFLECTIVITY@. This parameter is fixed for the
- whole visible range and defines the amount of reflected energy. Refer to the
- *-b* option to control how the light is reflected on the ground surface.
-
*-f*::
Force overwrite of the _output_ file.
*-g* _ground_::
- Path toward an OBJ file [2] representing the ground geometry.
+ Path toward a *htrdr-obj*(5) representing the ground geometry.
*-h*::
List short help and exit.
@@ -145,10 +138,17 @@ OPTIONS
@HTRDR_ARGS_DEFAULT_IMG_WIDTH@x@HTRDR_ARGS_DEFAULT_IMG_HEIGHT@.
**spp**=**_samples-count_**;;
- Number of samples per pixel and per component. i.e. the estimation of a
- pixel will use "3 * _samples-count_" Monte-Carlo realisations, one set of
+ 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. By default, *spp* is set to @HTRDR_ARGS_DEFAULT_IMG_SPP@.
+ XYZ color space. In long wave rendering (*-l* option) only one set of
+ _samples-count_ is used to estimate the pixel radiance and the resulting
+ brightness temperature for the submitted range of wavelengths. By default,
+ *spp* is set to @HTRDR_ARGS_DEFAULT_IMG_SPP@.
+
+*-l* __wlen-min__,__wlen-max__::
+ Switch in infrared rendering for the long wave interval in [__wlen-min__,
+ __wlen-max__] nanometers.
*-R*::
Infinitely repeat the _ground_ along the X and Y axis.
@@ -156,6 +156,9 @@ OPTIONS
*-r*::
Infinitely repeat the _clouds_ along the X and Y axis.
+*-M* _materials_::
+ Path toward a *htrdr-materials*(5) file listing the ground materials.
+
*-m* _mie_::
Path toward a *htmie*(5) file defining the optical properties of water
droplets.
@@ -165,7 +168,7 @@ OPTIONS
created and filled with the sky data built from the _clouds_, the
_atmosphere_ and the _mie_ input files. This cached data can then be reused
in the next runs as long as the input files provided on the command are the
- same of the ones used to setup the cache; leading to a significant speed-up
+ same as the ones used to setup the cache; leading to a significant speed-up
of the pre-processing step. If _cache_ contains data generated from input
files that are not the ones submitted on the command line, an error is
notified and the execution is stopped, avoiding the use of wrong cached data.
@@ -180,6 +183,7 @@ OPTIONS
*-T* _threshold_::
Optical thickness used as threshold criteria to partition the properties of
the clouds. By default its value is @HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@.
+ This option is ignored if a cache file is used (option *-O*).
*-t* _threads-count_::
Hint on the number of threads to use. By default use as many threads as CPU
@@ -188,7 +192,8 @@ OPTIONS
*-V* **_x_**,**_y_**,**_z_**::
Define the maximum definition of the acceleration data structure that
partitions the cloud properties. By default the finest definition is the
- definition of the submitted *htcp*(5) file.
+ definition of the submitted *htcp*(5) file. This option is ignored if a cache
+ file is used (option *-O*).
*-v*::
Make *htrdr* verbose.
@@ -201,33 +206,46 @@ EXAMPLES
Render a clear sky scene, i.e. a scene without any cloud, whose sun is at
zenith. The vertical atmospheric gaz mixture along the Z axis is described in
-the *gas.txt* file. The Mie data are provided through the *Mie.nc* file and
-the ground geometry is a quad repeated to the infinity. The camera is
-positioned at *400* meters high and looks toward the positive Y axis. The
-definition of the rendered image is *800* by *600* pixels and the radiance of
-each pixel component is estimated with *64* Monte-Carlo realisations. The
-resulting image is written to *output* excepted if the file already exists; in
-this case an error is notified, the program stops and the *output* file
-remains unchanged:
+the *gas.txt* file. the ground geometry is a quad repeated to the infinity
+whose materials are listed in the *material.mtl* file. The camera is positioned
+at *400* meters height and looks toward the positive Y axis. The definition of
+the rendered image is *800* by *600* pixels and the radiance of each pixel
+component is estimated with *64* Monte-Carlo realisations. The resulting image
+is written to *output* excepted if the file already exists; in this case an
+error is notified, the program stops and the *output* file remains unchanged:
$ htrdr -D0,90 -a gas.txt -m Mie.nc -g quad.obj -R \
+ -M materials.mtl \
-C pos=0,0,400:tgt=0,1,0:up=0,0,1 \
-i def=800x600:spp=64 \
-o output
Add clouds to the previous scene and use a more complex geometry to represent
-the ground. The ground geometry was carefully designed to be cyclic and can be
-thus repeated to the infinity without visual glitches. Use the *-f* option to
-write the rendered image to *output* even though the file already exists.
-Fianlly, use the *htpp*(1) command to convert the *htrdr-image*(5) saved in
-output in a regular PPM image [5]:
-
- $ htrdr -D0,90 -a gas.txt -m Mie.nc -g mountains.obj -R -c clouds.htcp \
+the ground. The Mie data are provided through the *Mie.nc* file. The ground
+geometry was carefully designed to be cyclic and can be thus infinitely
+repeated without visual glitches. Use the *-f* option to write the rendered
+image to *output* even though the file already exists. Finally, use the
+*htpp*(1) command to convert the *htrdr-image*(5) saved in output in a regular
+PPM image [4]:
+
+ $ htrdr -D0,90 -a gas.txt -m Mie.nc -g mountains.obj -R \
+ -M materials.mtl \
+ -c clouds.htcp \
-C pos=0,0,400:tgt=0,1,0:up=0,0,1 \
-i def=800x600:spp=64 \
-f -o output
$ htpp -o image.ppm output
+Render the previous scene in infrared for the wavelengths in [9200, 10000]
+nanometers:
+
+ $ htrdr -l 9200,10000 -a gas.txt -m Mie.nc -g mountains.obj -R \
+ -M materials.mtl \
+ -c clouds.htcp \
+ -C pos=0,0,400:tgt=0,1,0:up=0,0,1 \
+ -i def=800x600:spp=64 \
+ -f -o output
+
Move the sun by setting its azimuthal and elevation angles to *120* and *40*
degrees respectively. Use the *-O* option to enable the cache mechanism on
sky data. Increase the image definition to *1280* by *720* and set the
@@ -235,7 +253,9 @@ number of samples per pixel component to *1024*. Write results on standard
output and convert the resulting image in PPM before visualising it through the
*feh*(1) image viewer:
- $ htrdr -D120,40 -a gas.txt -m Mie.nc -g mountains.obj -R -c clouds.htcp \
+ $ htrdr -D120,40 -a gas.txt -m Mie.nc -g mountains.obj -R \
+ -M materials.mtl \
+ -c clouds.htcp \
-O my_cache \
-C pos=0,0,400:tgt=0,1,0:up=0,0,1 \
-i def=1280x720:spp=1024 | htpp | feh -
@@ -253,7 +273,9 @@ Use *mpirun*(1) to launch *htrdr* on several hosts defined in the *my_hosts*
file. Make the clouds infinite along the X and Y axis:
$ mpirun --hostfile my_hosts htrdr \
- -D120,40 -a gas.txt -m Mie.nc -g mountains.obj -R -c clouds.htcp -r \
+ -D120,40 -a gas.txt -m Mie.nc -g mountains.obj -R \
+ -M materials.mtl \
+ -c clouds.htcp -r \
-C pos=0,0,400:tgt=0,1,0:up=0,0,1 \
-i def=1280x720:spp=1024 \
-f -o output
@@ -262,12 +284,10 @@ NOTES
-----
1. High-Tune: Gas Optical Properties file format -
<https://www.meso-star.com/projects/high-tune/downloads/gas_opt_prop_en.pdf>
-2. OBJ file format -
- <http://www.martinreddy.net/gfx/3d/OBJ.spec>
-3. VTK file format -
+2. VTK file format -
<http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf>
-4. MPI specifications - <https://www.mpi-forum.org/docs/>
-5. Portable PixMap - <http://netpbm.sourceforge.net/doc/ppm.html>
+3. MPI specifications - <https://www.mpi-forum.org/docs/>
+4. Portable PixMap - <http://netpbm.sourceforge.net/doc/ppm.html>
COPYRIGHT
---------
@@ -286,4 +306,6 @@ SEE ALSO
*htcp*(5),
*htmie*(5),
*htpp*(1),
-*htrdr-image*(5)
+*htrdr-image*(5),
+*htrdr-materials*(5)
+*htrdr-obj*(5)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -20,8 +20,11 @@
#include "htrdr_c.h"
#include "htrdr_args.h"
#include "htrdr_buffer.h"
+#include "htrdr_cie_xyz.h"
#include "htrdr_camera.h"
#include "htrdr_ground.h"
+#include "htrdr_mtl.h"
+#include "htrdr_ran_lw.h"
#include "htrdr_sun.h"
#include "htrdr_solve.h"
@@ -100,7 +103,7 @@ log_msg
}
static res_T
-dump_accum_buffer
+dump_buffer
(struct htrdr* htrdr,
struct htrdr_buffer* buf,
struct htrdr_accum* time_acc, /* May be NULL */
@@ -108,18 +111,23 @@ dump_accum_buffer
FILE* stream)
{
struct htrdr_buffer_layout layout;
+ size_t pixsz, pixal;
size_t x, y;
res_T res = RES_OK;
ASSERT(htrdr && buf && stream_name && stream);
(void)stream_name;
+ if(htsky_is_long_wave(htrdr->sky)) {
+ pixsz = sizeof(struct htrdr_pixel_lw);
+ pixal = ALIGNOF(struct htrdr_pixel_lw);
+ } else {
+ pixsz = sizeof(struct htrdr_pixel_sw);
+ pixal = ALIGNOF(struct htrdr_pixel_sw);
+ }
+
htrdr_buffer_get_layout(buf, &layout);
- if(layout.elmt_size != sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__])
- || layout.alignment < ALIGNOF(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__])) {
- htrdr_log_err(htrdr,
- "%s: invalid buffer layout. "
- "The pixel size must be the size of %lu accumulators.\n",
- FUNC_NAME, (unsigned long)HTRDR_ESTIMATES_COUNT__);
+ if(layout.elmt_size != pixsz || layout.alignment != pixal) {
+ htrdr_log_err(htrdr, "%s: invalid buffer layout. ", FUNC_NAME);
res = RES_BAD_ARG;
goto error;
}
@@ -129,20 +137,33 @@ dump_accum_buffer
if(time_acc) *time_acc = HTRDR_ACCUM_NULL;
FOR_EACH(y, 0, layout.height) {
FOR_EACH(x, 0, layout.width) {
- const struct htrdr_accum* accums = htrdr_buffer_at(buf, x, y);
- int i;
- FOR_EACH(i, 0, HTRDR_ESTIMATES_COUNT__) {
- double E, SE;
-
- htrdr_accum_get_estimation(&accums[i], &E, &SE);
- fprintf(stream, "%g %g ", E, SE);
+ struct htrdr_estimate pix_time = HTRDR_ESTIMATE_NULL;
+ const struct htrdr_accum* pix_time_acc = NULL;
+
+ if(htsky_is_long_wave(htrdr->sky)) {
+ const struct htrdr_pixel_lw* pix = htrdr_buffer_at(buf, x, y);
+ 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_sw* 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);
+ fprintf(stream, "%g %g\n", pix_time.E, pix_time.SE);
+
if(time_acc) {
- time_acc->sum_weights += accums[HTRDR_ESTIMATE_TIME].sum_weights;
- time_acc->sum_weights_sqr += accums[HTRDR_ESTIMATE_TIME].sum_weights_sqr;
- time_acc->nweights += accums[HTRDR_ESTIMATE_TIME].nweights;
+ time_acc->sum_weights += pix_time_acc->sum_weights;
+ time_acc->sum_weights_sqr += pix_time_acc->sum_weights_sqr;
+ time_acc->nweights += pix_time_acc->nweights;
}
- fprintf(stream, "\n");
}
fprintf(stream, "\n");
}
@@ -382,9 +403,7 @@ htrdr_init
logger_set_stream(&htrdr->logger, LOG_OUTPUT, print_out, NULL);
logger_set_stream(&htrdr->logger, LOG_ERROR, print_err, NULL);
logger_set_stream(&htrdr->logger, LOG_WARNING, print_warn, NULL);
-
str_init(htrdr->allocator, &htrdr->output_name);
-
nthreads_max = MMAX(omp_get_max_threads(), omp_get_num_procs());
htrdr->dump_vtk = args->dump_vtk;
htrdr->verbose = args->verbose;
@@ -402,6 +421,9 @@ htrdr_init
if(!args->output) {
htrdr->output = stdout;
output_name = "<stdout>";
+ } else if(htrdr->mpi_rank != 0) {
+ htrdr->output = NULL;
+ output_name = "<null>";
} else {
res = open_output_stream
(htrdr, args->output, 0/*read*/, args->force_overwriting, &htrdr->output);
@@ -428,8 +450,14 @@ htrdr_init
goto error;
}
- res = htrdr_ground_create(htrdr, args->filename_obj, args->ground_bsdf_type,
- args->ground_reflectivity, args->repeat_ground, &htrdr->ground);
+ /* Materials are necessary only if a ground geometry is defined */
+ if(args->filename_obj) {
+ res = htrdr_mtl_create(htrdr, args->filename_mtl, &htrdr->mtl);
+ if(res != RES_OK) goto error;
+ }
+
+ res = htrdr_ground_create(htrdr, args->filename_obj, args->repeat_ground,
+ &htrdr->ground);
if(res != RES_OK) goto error;
proj_ratio =
@@ -439,22 +467,6 @@ htrdr_init
args->camera.up, proj_ratio, MDEG2RAD(args->camera.fov_y), &htrdr->cam);
if(res != RES_OK) goto error;
- /* 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) {
- /* 4 accums: X, Y, Z components and one more for the per realisation time */
- const size_t pixsz = sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]);
- const size_t pixal = ALIGNOF(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]);
- res = htrdr_buffer_create(htrdr,
- args->image.definition[0], /* Width */
- args->image.definition[1], /* Height */
- args->image.definition[0] * pixsz, /* Pitch */
- pixsz, /* Size of a pixel */
- pixal, /* Alignment of a pixel */
- &htrdr->buf);
- if(res != RES_OK) goto error;
- }
-
res = htrdr_sun_create(htrdr, &htrdr->sun);
if(res != RES_OK) goto error;
spherical_to_cartesian_dir
@@ -472,9 +484,33 @@ htrdr_init
htsky_args.nthreads = htrdr->nthreads;
htsky_args.repeat_clouds = args->repeat_clouds;
htsky_args.verbose = htrdr->mpi_rank == 0 ? args->verbose : 0;
+ htsky_args.wlen_lw_range[0] = args->wlen_lw_range[0];
+ htsky_args.wlen_lw_range[1] = args->wlen_lw_range[1];
res = htsky_create(&htrdr->logger, htrdr->allocator, &htsky_args, &htrdr->sky);
if(res != RES_OK) goto error;
+ if(!htsky_is_long_wave(htrdr->sky)) { /* Short wave random variate */
+ const double* range = HTRDR_CIE_XYZ_RANGE_DEFAULT;
+ size_t n;
+
+ n = (size_t)(range[1] - range[0]);
+ res = htrdr_cie_xyz_create(htrdr, range, n, &htrdr->cie);
+ if(res != RES_OK) goto error;
+
+ } else { /* Long Wave random variate */
+ const double Tref = 290; /* In Kelvin */
+ size_t n;
+
+ htrdr->wlen_range_m[0] = args->wlen_lw_range[0]*1e-9; /* Convert in meters */
+ htrdr->wlen_range_m[1] = args->wlen_lw_range[1]*1e-9; /* Convert in meters */
+ ASSERT(htrdr->wlen_range_m[0] <= htrdr->wlen_range_m[1]);
+
+ n = (size_t)(args->wlen_lw_range[1] - args->wlen_lw_range[0]);
+ res = htrdr_ran_lw_create
+ (htrdr, args->wlen_lw_range, n, Tref, &htrdr->ran_lw);
+ if(res != RES_OK) goto error;
+ }
+
htrdr->lifo_allocators = MEM_CALLOC
(htrdr->allocator, htrdr->nthreads, sizeof(*htrdr->lifo_allocators));
if(!htrdr->lifo_allocators) {
@@ -496,6 +532,29 @@ htrdr_init
}
}
+ /* Create the image buffer only on the master process; the image parts
+ * rendered by the processes are gathered onto the master process. */
+ if(!htrdr->dump_vtk && htrdr->mpi_rank == 0) {
+ size_t pixsz = 0; /* sizeof(pixel) */
+ size_t pixal = 0; /* alignof(pixel) */
+
+ if(htsky_is_long_wave(htrdr->sky)) {
+ pixsz = sizeof(struct htrdr_pixel_lw);
+ pixal = ALIGNOF(struct htrdr_pixel_lw);
+ } else {
+ pixsz = sizeof(struct htrdr_pixel_sw);
+ pixal = ALIGNOF(struct htrdr_pixel_sw);
+ }
+ res = htrdr_buffer_create(htrdr,
+ args->image.definition[0], /* Width */
+ args->image.definition[1], /* Height */
+ args->image.definition[0] * pixsz, /* Pitch */
+ pixsz, /* Size of a pixel */
+ pixal, /* Alignment of a pixel */
+ &htrdr->buf);
+ if(res != RES_OK) goto error;
+ }
+
exit:
return res;
error:
@@ -514,6 +573,10 @@ htrdr_release(struct htrdr* htrdr)
if(htrdr->sun) htrdr_sun_ref_put(htrdr->sun);
if(htrdr->cam) htrdr_camera_ref_put(htrdr->cam);
if(htrdr->buf) htrdr_buffer_ref_put(htrdr->buf);
+ if(htrdr->mtl) htrdr_mtl_ref_put(htrdr->mtl);
+ if(htrdr->cie) htrdr_cie_xyz_ref_put(htrdr->cie);
+ if(htrdr->ran_lw) htrdr_ran_lw_ref_put(htrdr->ran_lw);
+ if(htrdr->output && htrdr->output != stdout) fclose(htrdr->output);
if(htrdr->lifo_allocators) {
size_t i;
FOR_EACH(i, 0, htrdr->nthreads) {
@@ -530,15 +593,15 @@ htrdr_run(struct htrdr* htrdr)
{
res_T res = RES_OK;
if(htrdr->dump_vtk) {
- const size_t nbands = htsky_get_sw_spectral_bands_count(htrdr->sky);
+ const size_t nbands = htsky_get_spectral_bands_count(htrdr->sky);
size_t i;
/* Nothing to do */
if(htrdr->mpi_rank != 0) goto exit;
FOR_EACH(i, 0, nbands) {
- const size_t iband = htsky_get_sw_spectral_band_id(htrdr->sky, i);
- const size_t nquads = htsky_get_sw_spectral_band_quadrature_length
+ const size_t iband = htsky_get_spectral_band_id(htrdr->sky, i);
+ const size_t nquads = htsky_get_spectral_band_quadrature_length
(htrdr->sky, iband);
size_t iquad;
FOR_EACH(iquad, 0, nquads) {
@@ -548,20 +611,22 @@ htrdr_run(struct htrdr* htrdr)
}
}
} else {
- res = htrdr_draw_radiance_sw(htrdr, htrdr->cam, htrdr->width,
+ res = htrdr_draw_radiance(htrdr, htrdr->cam, 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;
- double E, SE;
+ struct htrdr_estimate path_time;
- res = dump_accum_buffer(htrdr, htrdr->buf, &path_time_acc,
+ res = dump_buffer(htrdr, htrdr->buf, &path_time_acc,
str_cget(&htrdr->output_name), htrdr->output);
if(res != RES_OK) goto error;
- htrdr_accum_get_estimation(&path_time_acc, &E, &SE);
+ htrdr_accum_get_estimation(&path_time_acc, &path_time);
htrdr_log(htrdr,
- "Time per radiative path (in micro seconds): %g +/- %g\n", E, SE);
+ "Time per radiative path (in micro seconds): %g +/- %g\n",
+ path_time.E,
+ path_time.SE);
}
}
exit:
@@ -644,7 +709,67 @@ htrdr_fflush(struct htrdr* htrdr, FILE* stream)
/*******************************************************************************
* Local functions
******************************************************************************/
-extern LOCAL_SYM res_T
+res_T
+brightness_temperature
+ (struct htrdr* htrdr,
+ const double lambda_min,
+ const double lambda_max,
+ const double radiance,
+ double* temperature)
+{
+ const size_t MAX_ITER = 100;
+ const double epsilon_T = 1e-4; /* In K */
+ const double epsilon_B = radiance * 1e-8;
+ double T, T0, T1, T2;
+ double B, B0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(temperature && lambda_min <= lambda_max);
+
+ /* Search for a brightness temperature whose radiance is greater than or
+ * equal to the estimated radiance */
+ T2 = 200;
+ FOR_EACH(i, 0, MAX_ITER) {
+ const double B2 = planck(lambda_min, lambda_max, T2);
+ if(B2 >= radiance) break;
+ T2 *= 2;
+ }
+ if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+
+ B0 = T0 = T1 = 0;
+ FOR_EACH(i, 0, MAX_ITER) {
+ T = (T1+T2)*0.5;
+ B = planck(lambda_min, lambda_max, T);
+
+ if(B < radiance) {
+ T1 = T;
+ } else {
+ T2 = T;
+ }
+
+ if(fabs(T-T0) < epsilon_T || fabs(B-B0) < epsilon_B)
+ break;
+
+ T0 = T;
+ B0 = B;
+ }
+ if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; }
+
+ *temperature = T;
+
+exit:
+ return res;
+error:
+ htrdr_log_err(htrdr,
+ "Could not compute the brightness temperature for the radiance %g "
+ "estimated over [%g, %g] nanometers.\n",
+ radiance,
+ lambda_min*1e9,
+ lambda_max*1e9);
+ goto exit;
+}
+
+res_T
open_output_stream
(struct htrdr* htrdr,
const char* filename,
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -34,6 +34,7 @@
struct htsky;
struct htrdr_args;
struct htrdr_buffer;
+struct htrdr_cie_xyz;
struct htrdr_rectangle;
struct mem_allocator;
struct mutext;
@@ -46,12 +47,16 @@ struct htrdr {
struct s3d_device* s3d;
struct htrdr_ground* ground;
+ struct htrdr_mtl* mtl;
struct htrdr_sun* sun;
+ struct htrdr_cie_xyz* cie;
+ struct htrdr_ran_lw* ran_lw;
struct htrdr_camera* cam;
struct htrdr_buffer* buf;
struct htsky* sky;
+ double wlen_range_m[2]; /* Integration range in *meters* */
size_t spp; /* #samples per pixel */
size_t width; /* Image width */
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -28,36 +28,21 @@
/*******************************************************************************
* Helper functions
******************************************************************************/
-static const char*
-bsdf_type_to_string(const enum htrdr_bsdf_type type)
-{
- const char* str = "<none>";
- switch(type) {
- case HTRDR_BSDF_DIFFUSE: str = "diffuse"; break;
- case HTRDR_BSDF_SPECULAR: str = "specular"; break;
- default: FATAL("Unreachable code.\n"); break;
- }
- return str;
-}
-
static void
print_help(const char* cmd)
{
ASSERT(cmd);
printf("Usage: %s [OPION]... -a ATMOSPHERE\n", cmd);
printf(
-"Render an image in the visible part of the spectrum, for scenes composed of an\n"
-"atmospheric gaz mixture, clouds and a ground.\n\n");
+"Render an image for scenes composed of an atmospheric gas mixture,\n"
+"clouds and a ground.\n\n");
printf(
" -a ATMOSPHERE gas optical properties of the atmosphere.\n");
printf(
-" -b <diffuse|specular>\n"
-" BSDF of the ground. Default value is %s.\n",
- bsdf_type_to_string(HTRDR_ARGS_DEFAULT.ground_bsdf_type));
- printf(
" -c CLOUDS properties of the clouds.\n");
printf(
-" -C <camera> define the rendering point of view.\n");
+" -C <camera> define the rendering point of view. Refer to the\n"
+" htrdr man page for the list of camera options.\n");
printf(
" -D AZIMUTH,ELEVATION\n"
" direction in degrees toward the sun center. By default\n"
@@ -67,24 +52,31 @@ print_help(const char* cmd)
printf(
" -d dump octrees data to OUTPUT and exit.\n");
printf(
-" -e REFLECT ground reflectivity in [0, 1]. Default value is %g.\n",
- HTRDR_ARGS_DEFAULT.ground_reflectivity);
- printf(
" -f overwrite the OUTPUT file if it already exists.\n");
printf(
" -g GROUND ground geometry.\n");
printf(
" -h display this help and exit.\n");
printf(
-" -i <image> define the image to compute.\n");
+" -i <image> define the image to compute. Refer to the htrdr man\n"
+" page for the list of image options\n");
+ printf(
+" -l WLEN_MIN,WLEN_MAX\n"
+" switch in infrared rendering for the long waves in\n"
+" [WLEN_MIN, WLEN_MAX], in nanometers. By default, the\n"
+" rendering is performed for the visible part of the\n"
+" spectrum in [380, 780] nanometers.\n");
printf(
" -R infinitely repeat the ground along the X and Y axis.\n");
printf(
" -r infinitely repeat the clouds along the X and Y axis.\n");
printf(
+" -M MATERIALS file listing the ground materials.\n");
+ printf(
" -m MIE file of Mie's data.\n");
printf(
-" -O CACHE name of the cache file used to store/restore the sky data.\n");
+" -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");
@@ -93,21 +85,23 @@ print_help(const char* cmd)
" building. By default its value is `%g'.\n",
HTRDR_ARGS_DEFAULT.optical_thickness);
printf(
-" -t THREADS hint on the number of threads to use. By default use as\n"
-" many threads as CPU cores.\n");
+" -t THREADS hint on the number of threads to use. By default use\n"
+" as many threads as CPU cores.\n");
printf(
-" -V X,Y,Z maximum definition of the cloud acceleration grids along\n"
-" the 3 axis. By default use the definition of the clouds\n");
+" -V X,Y,Z maximum definition of the cloud acceleration grids\n"
+" along the 3 axis. By default use the definition of\n"
+" the clouds\n");
printf(
" -v make the program verbose.\n");
printf(
" --version display version information and exit.\n");
printf("\n");
printf(
-"Copyright (C) 2018, 2019, 2020 |Meso|Star> <contact@meso-star.com>. Copyright\n"
-"(C) 2018, 2019 CNRS, Université Paul Sabatier. htrdr is free software released\n"
-"under the GNU GPL license, version 3 or later. You are free to change or\n"
-"redistribute it under certain conditions <http://gnu.org/licenses/gpl.html>\n");
+"Copyright (C) 2018, 2019, 2020 |Meso|Star> <contact@meso-star.com>.\n"
+"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");
}
static INLINE res_T
@@ -372,19 +366,24 @@ error:
}
static res_T
-parse_bsdf_type(struct htrdr_args* args, const char* str)
+parse_lw_range(struct htrdr_args* args, const char* str)
{
+ double range[2];
+ size_t len;
res_T res = RES_OK;
- if(!strcmp(str, "diffuse")) {
- args->ground_bsdf_type = HTRDR_BSDF_DIFFUSE;
- } else if(!strcmp(str, "specular")) {
- args->ground_bsdf_type = HTRDR_BSDF_SPECULAR;
- } else {
- fprintf(stderr, "Invalid BRDF type `%s'.\n", str);
- res = RES_BAD_ARG;
+ ASSERT(args && str);
+
+ res = cstr_to_list_double(str, ',', range, &len, 2);
+ if(res == RES_OK && len != 2) res = RES_BAD_ARG;
+ if(res == RES_OK && range[0] > range[1]) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ fprintf(stderr, "Invalid long wave range `%s'.\n", str);
goto error;
}
+ args->wlen_lw_range[0] = range[0];
+ args->wlen_lw_range[1] = range[1];
+
exit:
return res;
error:
@@ -415,25 +414,16 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
}
}
- while((opt = getopt(argc, argv, "a:b:C:c:D:de:fg:hi:m:O:o:RrT:t:V:v")) != -1) {
+ while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:l:M:m:O:o:RrT:t:V:v")) != -1) {
switch(opt) {
case 'a': args->filename_gas = optarg; break;
- case 'b':
- res = parse_bsdf_type(args, optarg);
- break;
- case 'C':
+ case 'C':
res = parse_multiple_parameters
(args, optarg, parse_camera_parameter);
break;
case 'c': args->filename_les = optarg; break;
case 'D': res = parse_sun_dir(args, optarg); break;
case 'd': args->dump_vtk = 1; break;
- case 'e':
- res = cstr_to_double(optarg, &args->ground_reflectivity);
- if(args->ground_reflectivity < 0 || args->ground_reflectivity > 1) {
- res = RES_BAD_ARG;
- }
- break;
case 'f': args->force_overwriting = 1; break;
case 'g': args->filename_obj = optarg; break;
case 'h':
@@ -445,6 +435,10 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
res = parse_multiple_parameters
(args, optarg, parse_image_parameter);
break;
+ case 'l':
+ res = parse_lw_range(args, optarg);
+ break;
+ case 'M': args->filename_mtl = optarg; break;
case 'm': args->filename_mie = optarg; break;
case 'O': args->cache = optarg; break;
case 'o': args->output = optarg; break;
@@ -476,6 +470,12 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
res = RES_BAD_ARG;
goto error;
}
+ if(args->filename_obj && !args->filename_mtl) {
+ fprintf(stderr,
+ "Missing the path of the file listing the ground materials -- option '-M'\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
if(args->filename_les && !args->filename_mie) {
fprintf(stderr,
"Missing the path toward the file of the Mie's data -- option '-m'\n");
diff --git a/src/htrdr_args.h.in b/src/htrdr_args.h.in
@@ -18,6 +18,7 @@
#include "htrdr_ground.h"
+#include <float.h>
#include <limits.h>
#include <rsys/rsys.h>
@@ -26,6 +27,7 @@ struct htrdr_args {
const char* filename_les; /* Path of the HTCP file */
const char* filename_mie; /* Path of the Mie properties */
const char* filename_obj; /* Path of the 3D geometry */
+ const char* filename_mtl; /* Path of the materials */
const char* cache;
const char* output;
@@ -51,10 +53,10 @@ struct htrdr_args {
double sun_azimuth; /* In degrees */
double sun_elevation; /* In degrees */
double optical_thickness; /* Threshold used during octree building */
- enum htrdr_bsdf_type ground_bsdf_type;
- double ground_reflectivity; /* Reflectivity of the ground */
unsigned grid_max_definition[3]; /* Maximum definition of the grid */
+ double wlen_lw_range[2]; /* Long Wave range to handle in nm */
+
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 */
@@ -69,6 +71,7 @@ struct htrdr_args {
NULL, /* LES filename */ \
NULL, /* Mie filename */ \
NULL, /* Obj filename */ \
+ NULL, /* Mtl filename */ \
NULL, /* Cache filename */ \
NULL, /* Output filename */ \
{ \
@@ -88,9 +91,8 @@ struct htrdr_args {
0, /* Sun azimuth */ \
90, /* Sun elevation */ \
@HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@, /* Optical thickness */ \
- @HTRDR_ARGS_DEFAULT_GROUND_BSDF@, \
- @HTRDR_ARGS_DEFAULT_GROUND_REFLECTIVITY@, /* Ground reflectivity */ \
{UINT_MAX, UINT_MAX, UINT_MAX}, /* Maximum definition of the grid */ \
+ {DBL_MAX, -DBL_MAX}, /* Long wave range. Degenerated <=> short wave */ \
(unsigned)~0, /* #threads */ \
0, /* Force overwriting */ \
0, /* dump VTK */ \
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -18,6 +18,7 @@
#define HTRDR_C_H
#include <rsys/rsys.h>
+#include <rsys/math.h> /* PI support */
#ifndef NDEBUG
#define MPI(Func) ASSERT(MPI_##Func == MPI_SUCCESS)
@@ -33,14 +34,6 @@ enum htrdr_mpi_message {
HTRDR_MPI_TILE_DATA
};
-enum htrdr_estimate {
- HTRDR_ESTIMATE_X,
- HTRDR_ESTIMATE_Y,
- HTRDR_ESTIMATE_Z,
- HTRDR_ESTIMATE_TIME, /* Time per realisation */
- HTRDR_ESTIMATES_COUNT__
-};
-
struct htrdr;
#define SW_WAVELENGTH_MIN 380 /* In nanometer */
@@ -103,6 +96,111 @@ morton_xyz_decode_u21(const uint64_t code, uint32_t xyz[3])
xyz[2] = (uint32_t)morton3D_decode_u21(code >> 0);
}
+static INLINE double
+wiebelt(const double v)
+{
+ int m;
+ double w, v2, v4;
+ /*.153989717364e+00;*/
+ const double fifteen_over_pi_power_4 = 15.0/(PI*PI*PI*PI);
+ const double z0 = 1.0/3.0;
+ const double z1 = 1.0/8.0;
+ const double z2 = 1.0/60.0;
+ const double z4 = 1.0/5040.0;
+ const double z6 = 1.0/272160.0;
+ const double z8 = 1.0/13305600.0;
+
+ if(v >= 2.) {
+ w = 0;
+ for(m=1; m<6 ;m++)
+ w+=exp(-m*v)/(m*m*m*m) * (((m*v+3)*m*v+6)*m*v+6);
+ w = w * fifteen_over_pi_power_4;
+ } else {
+ v2 = v*v;
+ v4 = v2*v2;
+ w = z0 - z1*v + z2*v2 - z4*v2*v2 + z6*v4*v2 - z8*v4*v4;
+ w = 1. - fifteen_over_pi_power_4*v2*v*w;
+ }
+ ASSERT(w >= 0.0 && w <= 1.0);
+ return w;
+}
+
+static INLINE double
+blackbody_fraction
+ (const double lambda0, /* In meters */
+ const double lambda1, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ const double C2 = 1.43877735e-2; /* m.K */
+ double x0 = C2 / lambda0;
+ double x1 = C2 / lambda1;
+ double v0 = x0 / temperature;
+ double v1 = x1 / temperature;
+ double w0 = wiebelt(v0);
+ double w1 = wiebelt(v1);
+ return w1 - w0;
+}
+
+
+
+/* Return the Planck value in W/m^2/sr/m at a given wavelength */
+static INLINE double
+planck_monochromatic
+ (const double lambda, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ const double c = 299792458; /* m/s */
+ const double h = 6.62607015e-34; /* J.s */
+ const double k = 1.380649e-23; /* J/K */
+ const double lambda2 = lambda*lambda;
+ const double lambda5 = lambda2*lambda2*lambda;
+ const double B = ((2.0 * h * c*c) / lambda5) /* W/m^2/sr/m */
+ / (exp(h*c/(lambda*k*temperature))-1.0);
+ ASSERT(temperature > 0);
+ return B;
+}
+
+/* Return the average Planck value in W/m^2/sr/m over the
+ * [lambda_min, lambda_max] interval. */
+static INLINE double
+planck_interval
+ (const double lambda_min, /* In meters */
+ const double lambda_max, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ const double T2 = temperature*temperature;
+ const double T4 = T2*T2;
+ const double BOLTZMANN_CONSTANT = 5.6696e-8; /* W/m^2/K^4 */
+ ASSERT(lambda_min < lambda_max && temperature > 0);
+ return blackbody_fraction(lambda_min, lambda_max, temperature)
+ * BOLTZMANN_CONSTANT * T4 / (PI * (lambda_max-lambda_min)); /* In W/m^2/sr/m */
+}
+
+/* Invoke planck_monochromatic or planck_interval whether the submitted
+ * interval is null or not, respectively. The returned value is in W/m^2/sr/m */
+static FINLINE double
+planck
+ (const double lambda_min, /* In meters */
+ const double lambda_max, /* In meters */
+ const double temperature) /* In Kelvin */
+{
+ ASSERT(lambda_min <= lambda_max && temperature > 0);
+ if(lambda_min == lambda_max) {
+ return planck_monochromatic(lambda_min, temperature);
+ } else {
+ return planck_interval(lambda_min, lambda_max, temperature);
+ }
+}
+
+extern LOCAL_SYM res_T
+brightness_temperature
+ (struct htrdr* htrdr,
+ const double lambda_min, /* In meters */
+ const double lambda_max, /* In meters */
+ /* Integrated over [lambda_min, lambda_max]. In W/m^2/sr/m */
+ const double radiance,
+ double* temperature);
+
extern LOCAL_SYM res_T
open_output_stream
(struct htrdr* htrdr,
@@ -146,5 +244,13 @@ update_mpi_progress(struct htrdr* htrdr, const enum htrdr_mpi_message progress)
print_mpi_progress(htrdr, progress);
}
+static FINLINE int
+cmp_dbl(const void* a, const void* b)
+{
+ const double d0 = *((const double*)a);
+ const double d1 = *((const double*)b);
+ return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0);
+}
+
#endif /* HTRDR_C_H */
diff --git a/src/htrdr_cie_xyz.c b/src/htrdr_cie_xyz.c
@@ -0,0 +1,355 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "htrdr.h"
+#include "htrdr_c.h"
+#include "htrdr_cie_xyz.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <math.h> /* nextafter */
+
+struct htrdr_cie_xyz {
+ struct darray_double cdf_X;
+ struct darray_double cdf_Y;
+ struct darray_double cdf_Z;
+ double range[2]; /* Boundaries of the handled CIE XYZ color space */
+ double band_len; /* Length in nanometers of a band */
+
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE double
+trapezoidal_integration
+ (const double lambda_lo, /* Integral lower bound. In nanometer */
+ const double lambda_hi, /* Integral upper bound. In nanometer */
+ double (*f_bar)(const double lambda)) /* Function to integrate */
+{
+ double dlambda;
+ size_t i, n;
+ double integral = 0;
+ ASSERT(lambda_lo <= lambda_hi);
+ ASSERT(lambda_lo > 0);
+
+ n = (size_t)(lambda_hi - lambda_lo) + 1;
+ dlambda = (lambda_hi - lambda_lo) / (double)n;
+
+ FOR_EACH(i, 0, n) {
+ const double lambda1 = lambda_lo + dlambda*(double)(i+0);
+ const double lambda2 = lambda_hi + dlambda*(double)(i+1);
+ const double f1 = f_bar(lambda1);
+ const double f2 = f_bar(lambda2);
+ integral += (f1 + f2)*dlambda*0.5;
+ }
+ return integral;
+}
+
+/* The following 3 functions are used to fit the CIE Xbar, Ybar and Zbar curved
+ * has defined by the 1931 standard. These analytical fits are propsed by C.
+ * Wyman, P. P. Sloan & P. Shirley in "Simple Analytic Approximations to the
+ * CIE XYZ Color Matching Functions" - JCGT 2013. */
+static INLINE double
+fit_x_bar_1931(const double lambda)
+{
+ const double a = (lambda - 442.0) * (lambda < 442.0 ? 0.0624 : 0.0374);
+ const double b = (lambda - 599.8) * (lambda < 599.8 ? 0.0264 : 0.0323);
+ const double c = (lambda - 501.1) * (lambda < 501.1 ? 0.0490 : 0.0382);
+ return 0.362*exp(-0.5*a*a) + 1.056*exp(-0.5f*b*b) - 0.065*exp(-0.5*c*c);
+}
+
+static FINLINE double
+fit_y_bar_1931(const double lambda)
+{
+ const double a = (lambda - 568.8) * (lambda < 568.8 ? 0.0213 : 0.0247);
+ const double b = (lambda - 530.9) * (lambda < 530.9 ? 0.0613 : 0.0322);
+ return 0.821*exp(-0.5*a*a) + 0.286*exp(-0.5*b*b);
+}
+
+static FINLINE double
+fit_z_bar_1931(const double lambda)
+{
+ const double a = (lambda - 437.0) * (lambda < 437.0 ? 0.0845 : 0.0278);
+ const double b = (lambda - 459.0) * (lambda < 459.0 ? 0.0385 : 0.0725);
+ return 1.217*exp(-0.5*a*a) + 0.681*exp(-0.5*b*b);
+}
+
+static INLINE double
+sample_cie_xyz
+ (const struct htrdr_cie_xyz* cie,
+ const double* cdf,
+ const size_t cdf_length,
+ double (*f_bar)(const double lambda), /* Function to integrate */
+ const double r0, /* Canonical number in [0, 1[ */
+ const double r1) /* Canonical number in [0, 1[ */
+{
+ double r0_next = nextafter(r0, DBL_MAX);
+ double* find;
+ double f_min, f_max; /* CIE 1931 value for the band boundaries */
+ double lambda; /* Sampled wavelength */
+ double lambda_min, lambda_max; /* Boundaries of the sampled band */
+ double lambda_1, lambda_2; /* Solutions if the equation to solve */
+ double a, b, c, d; /* Equation parameters */
+ double delta, sqrt_delta;
+ size_t iband; /* Index of the sampled band */
+ ASSERT(cie && cdf && cdf_length);
+ ASSERT(0 <= r0 && r0 < 1);
+ ASSERT(0 <= r1 && r1 < 1);
+
+ /* Use r_next rather than r in order to find the first entry that is not less
+ * than *or equal* to r */
+ find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ ASSERT(find);
+
+ /* Define and check the sampled band */
+ iband = (size_t)(find - cdf);
+ ASSERT(iband < cdf_length);
+ ASSERT(cdf[iband] > r0 && (!iband || cdf[iband-1] <= r0));
+
+ /* Define the boundaries of the sampled band */
+ lambda_min = cie->range[0] + cie->band_len * (double)iband;
+ lambda_max = lambda_min + cie->band_len;
+
+ /* Define the value of the CIE 1931 function for the boudaries of the sampled
+ * band */
+ f_min = f_bar(lambda_min);
+ f_max = f_bar(lambda_max);
+
+ /* Compute the equation constants */
+ a = 0.5 * (f_max - f_min) / cie->band_len;
+ b = (lambda_max * f_min - lambda_min * f_max) / cie->band_len;
+ c = -lambda_min * f_min + lambda_min*lambda_min * a;
+ d = 0.5 * (f_max + f_min) * cie->band_len;
+
+ delta = b*b - 4*a*(c-d*r1);
+ if(delta < 0 && eq_eps(delta, 0, 1.e-6)) {
+ delta = 0;
+ }
+ ASSERT(delta > 0);
+ sqrt_delta = sqrt(delta);
+
+ /* Compute the roots that solve the equation */
+ lambda_1 = (-b - sqrt_delta) / (2*a);
+ lambda_2 = (-b + sqrt_delta) / (2*a);
+
+ /* Select the solution */
+ if(lambda_min <= lambda_1 && lambda_1 < lambda_max) {
+ lambda = lambda_1;
+ } else if(lambda_min <= lambda_2 && lambda_2 < lambda_max) {
+ lambda = lambda_2;
+ } else {
+ FATAL("Unexpected error.\n");
+ }
+
+ return lambda;
+}
+
+static res_T
+setup_cie_xyz
+ (struct htrdr_cie_xyz* cie,
+ const char* func_name,
+ const double range[2],
+ const size_t nbands)
+{
+ enum { X, Y, Z }; /* Helper constant */
+ double* pdf[3] = {NULL, NULL, NULL};
+ double* cdf[3] = {NULL, NULL, NULL};
+ double sum[3] = {0,0,0};
+ size_t i;
+ res_T res = RES_OK;
+
+ ASSERT(cie && func_name && range && nbands);
+ ASSERT(range[0] >= HTRDR_CIE_XYZ_RANGE_DEFAULT[0]);
+ ASSERT(range[1] <= HTRDR_CIE_XYZ_RANGE_DEFAULT[1]);
+ ASSERT(range[0] < range[1]);
+
+ /* Allocate and reset the memory space for the tristimulus CDF */
+ #define SETUP_STIMULUS(Stimulus) { \
+ res = darray_double_resize(&cie->cdf_ ## Stimulus, nbands); \
+ if(res != RES_OK) { \
+ htrdr_log_err(cie->htrdr, \
+ "%s: Could not reserve the memory space for the CDF " \
+ "of the "STR(X)" stimulus -- %s.\n", func_name, res_to_cstr(res)); \
+ goto error; \
+ } \
+ cdf[Stimulus] = darray_double_data_get(&cie->cdf_ ## Stimulus); \
+ pdf[Stimulus] = cdf[Stimulus]; \
+ memset(cdf[Stimulus], 0, nbands*sizeof(double)); \
+ } (void)0
+ SETUP_STIMULUS(X);
+ SETUP_STIMULUS(Y);
+ SETUP_STIMULUS(Z);
+ #undef SETUP_STIMULUS
+
+ /* Compute the *unormalized* pdf of the tristimulus */
+ cie->range[0] = range[0];
+ cie->range[1] = range[1];
+ cie->band_len = (range[1] - range[0]) / (double)nbands;
+ FOR_EACH(i, 0, nbands) {
+ const double lambda_lo = range[0] + (double)i * cie->band_len;
+ const double lambda_hi = lambda_lo + cie->band_len;
+ ASSERT(lambda_lo <= lambda_hi);
+ pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
+ pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931);
+ pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931);
+ sum[X] += pdf[X][i];
+ sum[Y] += pdf[Y][i];
+ sum[Z] += pdf[Z][i];
+ }
+
+ FOR_EACH(i, 0, nbands) {
+ /* Normalize the pdf */
+ pdf[X][i] /= sum[X];
+ pdf[Y][i] /= sum[Y];
+ pdf[Z][i] /= sum[Z];
+ /* Setup the cumulative */
+ if(i == 0) {
+ cdf[X][i] = pdf[X][i];
+ cdf[Y][i] = pdf[Y][i];
+ cdf[Z][i] = pdf[Z][i];
+ } else {
+ cdf[X][i] = pdf[X][i] + cdf[X][i-1];
+ cdf[Y][i] = pdf[Y][i] + cdf[Y][i-1];
+ cdf[Z][i] = pdf[Z][i] + cdf[Z][i-1];
+ ASSERT(cdf[X][i] >= cdf[X][i-1]);
+ ASSERT(cdf[Y][i] >= cdf[Y][i-1]);
+ ASSERT(cdf[Z][i] >= cdf[Z][i-1]);
+ }
+ }
+ ASSERT(eq_eps(cdf[X][nbands-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Y][nbands-1], 1, 1.e-6));
+ ASSERT(eq_eps(cdf[Z][nbands-1], 1, 1.e-6));
+
+ /* Handle numerical issue */
+ cdf[X][nbands-1] = 1.0;
+ cdf[Y][nbands-1] = 1.0;
+ cdf[Z][nbands-1] = 1.0;
+
+exit:
+ return res;
+error:
+ darray_double_clear(&cie->cdf_X);
+ darray_double_clear(&cie->cdf_Y);
+ darray_double_clear(&cie->cdf_Z);
+ goto exit;
+}
+
+static void
+release_cie_xyz(ref_T* ref)
+{
+ struct htrdr_cie_xyz* cie = NULL;
+ ASSERT(ref);
+ cie = CONTAINER_OF(ref, struct htrdr_cie_xyz, ref);
+ darray_double_release(&cie->cdf_X);
+ darray_double_release(&cie->cdf_Y);
+ darray_double_release(&cie->cdf_Z);
+ MEM_RM(cie->htrdr->allocator, cie);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_cie_xyz_create
+ (struct htrdr* htrdr,
+ const double range[2], /* Must be included in [380, 780] nanometers */
+ const size_t nbands, /* # bands used to discretisze the CIE tristimulus */
+ struct htrdr_cie_xyz** out_cie)
+{
+ struct htrdr_cie_xyz* cie = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && range && nbands && out_cie);
+
+ cie = MEM_CALLOC(htrdr->allocator, 1, sizeof(*cie));
+ if(!cie) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: could not allocate the CIE XYZ data structure -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ ref_init(&cie->ref);
+ cie->htrdr = htrdr;
+ darray_double_init(htrdr->allocator, &cie->cdf_X);
+ darray_double_init(htrdr->allocator, &cie->cdf_Y);
+ darray_double_init(htrdr->allocator, &cie->cdf_Z);
+
+ res = setup_cie_xyz(cie, FUNC_NAME, range, nbands);
+ if(res != RES_OK) goto error;
+
+ htrdr_log(htrdr, "Short wave interval defined on [%g, %g] nanometers.\n",
+ range[0], range[1]);
+
+exit:
+ *out_cie = cie;
+ return res;
+error:
+ if(cie) htrdr_cie_xyz_ref_put(cie);
+ goto exit;
+}
+
+void
+htrdr_cie_xyz_ref_get(struct htrdr_cie_xyz* cie)
+{
+ ASSERT(cie);
+ ref_get(&cie->ref);
+}
+
+void
+htrdr_cie_xyz_ref_put(struct htrdr_cie_xyz* cie)
+{
+ ASSERT(cie);
+ ref_put(&cie->ref, release_cie_xyz);
+}
+
+double
+htrdr_cie_xyz_sample_X
+ (struct htrdr_cie_xyz* cie,
+ const double r0,
+ const double r1)
+{
+ return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
+ darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
+}
+
+double
+htrdr_cie_xyz_sample_Y
+ (struct htrdr_cie_xyz* cie,
+ const double r0,
+ const double r1)
+{
+ return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
+ darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
+}
+
+double
+htrdr_cie_xyz_sample_Z
+ (struct htrdr_cie_xyz* cie,
+ const double r0,
+ const double r1)
+{
+ return sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
+ darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
+}
+
diff --git a/src/htrdr_cie_xyz.h b/src/htrdr_cie_xyz.h
@@ -0,0 +1,61 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_CIE_XYZ_H
+#define HTRDR_CIE_XYZ_H
+
+#include <rsys/rsys.h>
+
+struct htrdr;
+struct htrdr_cie_xyz;
+
+/* Wavelength boundaries of the CIE XYZ color space in nanometers */
+static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] = {380, 780};
+
+extern LOCAL_SYM res_T
+htrdr_cie_xyz_create
+ (struct htrdr* htrdr,
+ const double range[2], /* Must be included in [380, 780] nanometers */
+ const size_t nbands, /* # bands used to discretisze the CIE tristimulus s*/
+ struct htrdr_cie_xyz** cie);
+
+extern LOCAL_SYM void
+htrdr_cie_xyz_ref_get
+ (struct htrdr_cie_xyz* cie);
+
+extern LOCAL_SYM void
+htrdr_cie_xyz_ref_put
+ (struct htrdr_cie_xyz* cie);
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_cie_xyz_sample_X
+ (struct htrdr_cie_xyz* cie,
+ const double r0, const double r1); /* Canonical numbers in [0, 1[ */
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_cie_xyz_sample_Y
+ (struct htrdr_cie_xyz* cie,
+ const double r0, const double r1); /* Canonical number in [0, 1[ */
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_cie_xyz_sample_Z
+ (struct htrdr_cie_xyz* cie,
+ const double r0, const double r1); /* Canonical number in [0, 1[ */
+
+#endif /* HTRDR_cie_xyz_H */
+
diff --git a/src/htrdr_compute_radiance_lw.c b/src/htrdr_compute_radiance_lw.c
@@ -0,0 +1,285 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "htrdr.h"
+#include "htrdr_c.h"
+#include "htrdr_interface.h"
+#include "htrdr_ground.h"
+#include "htrdr_solve.h"
+
+#include <high_tune/htsky.h>
+
+#include <star/s3d.h>
+#include <star/ssf.h>
+#include <star/ssp.h>
+#include <star/svx.h>
+
+#include <rsys/double2.h>
+#include <rsys/double3.h>
+
+enum event {
+ EVENT_ABSORPTION,
+ EVENT_SCATTERING,
+ EVENT_NONE
+};
+
+struct filter_context {
+ struct ssp_rng* rng;
+ const struct htsky* htsky;
+ size_t iband; /* Index of the spectral band */
+ size_t iquad; /* Index of the quadrature point into the band */
+
+ double Ts; /* Sampled optical thickness */
+ double traversal_dst; /* Distance traversed along the ray */
+
+ enum event event_type;
+};
+static const struct filter_context FILTER_CONTEXT_NULL = {
+ NULL, NULL, 0, 0, 0.0, 0.0, EVENT_NONE
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static int
+hit_filter
+ (const struct svx_hit* hit,
+ const double org[3],
+ const double dir[3],
+ const double range[2],
+ void* context)
+{
+ struct filter_context* ctx = context;
+ double kext_max;
+ int pursue_traversal = 1;
+ ASSERT(hit && ctx && !SVX_HIT_NONE(hit) && org && dir && range);
+ (void)range;
+
+ kext_max = htsky_fetch_svx_voxel_property(ctx->htsky, HTSKY_Kext,
+ HTSKY_SVX_MAX, HTSKY_CPNT_MASK_ALL, ctx->iband, ctx->iquad, &hit->voxel);
+
+ ctx->traversal_dst = hit->distance[0];
+
+ for(;;) {
+ double vox_dst = hit->distance[1] - ctx->traversal_dst;
+ const double T = vox_dst * kext_max;
+
+ /* A collision occurs behind `vox_dst' */
+ if(ctx->Ts > T) {
+ ctx->Ts -= T;
+ ctx->traversal_dst = hit->distance[1];
+ pursue_traversal = 1;
+ break;
+
+ /* A real/null collision occurs before `vox_dst' */
+ } else {
+ const double collision_dst = ctx->Ts / kext_max;
+ double pos[3];
+ double ks, ka;
+ double r;
+ double proba_abs;
+ double proba_sca;
+
+ /* Compute the traversed distance up to the challenged collision */
+ ctx->traversal_dst += collision_dst;
+ ASSERT(ctx->traversal_dst >= hit->distance[0]);
+ ASSERT(ctx->traversal_dst <= hit->distance[1]);
+
+ /* Compute the world space position where a collision may occur */
+ pos[0] = org[0] + ctx->traversal_dst * dir[0];
+ pos[1] = org[1] + ctx->traversal_dst * dir[1];
+ pos[2] = org[2] + ctx->traversal_dst * dir[2];
+
+ ka = htsky_fetch_raw_property(ctx->htsky, HTSKY_Ka, HTSKY_CPNT_MASK_ALL,
+ ctx->iband, ctx->iquad, pos, -DBL_MAX, DBL_MAX);
+ ks = htsky_fetch_raw_property(ctx->htsky, HTSKY_Ks, HTSKY_CPNT_MASK_ALL,
+ ctx->iband, ctx->iquad, pos, -DBL_MAX, DBL_MAX);
+
+ r = ssp_rng_canonical(ctx->rng);
+ proba_abs = ka / kext_max;
+ proba_sca = ks / kext_max;
+ if(r < proba_abs) { /* Absorption event */
+ pursue_traversal = 0;
+ ctx->event_type = EVENT_ABSORPTION;
+ break;
+ } else if(r < proba_abs + proba_sca) { /* Scattering event */
+ pursue_traversal = 0;
+ ctx->event_type = EVENT_SCATTERING;
+ break;
+ } else { /* Null collision */
+ ctx->Ts = ssp_ran_exp(ctx->rng, 1); /* Sample a new optical thickness */
+ }
+ }
+ }
+ return pursue_traversal;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+double
+htrdr_compute_radiance_lw
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ struct ssp_rng* rng,
+ const double pos_in[3],
+ const double dir_in[3],
+ const double wlen, /* In nanometer */
+ const size_t iband,
+ const size_t iquad)
+{
+ struct s3d_hit s3d_hit = S3D_HIT_NULL;
+ struct s3d_hit s3d_hit_prev = S3D_HIT_NULL;
+ struct svx_hit svx_hit = SVX_HIT_NULL;
+ struct ssf_phase* phase_hg = NULL;
+
+ double wo[3];
+ double pos[3];
+ double dir[3];
+ double range[2];
+ double pos_next[3];
+ double dir_next[3];
+ double band_bounds[2]; /* In nanometers */
+ double band_bounds_m[2]; /* In meters */
+ double temperature;
+ double g;
+ double w = 0; /* Weight */
+
+ ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads);
+
+ /* Retrieve the band boundaries */
+ htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
+
+ /* Transform the band boundaries in meters and clamp them to the integration
+ * domain */
+ band_bounds_m[0] = MMAX(band_bounds[0] * 1e-9, htrdr->wlen_range_m[0]);
+ band_bounds_m[1] = MMIN(band_bounds[1] * 1e-9, htrdr->wlen_range_m[1]);
+
+ /* Setup the phase function for this spectral band & quadrature point */
+ CHK(RES_OK == ssf_phase_create
+ (&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
+ g = htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter
+ (htrdr->sky, wlen);
+ SSF(phase_hg_setup(phase_hg, g));
+
+ /* Initialise the random walk */
+ d3_set(pos, pos_in);
+ d3_set(dir, dir_in);
+ d2(range, 0, INF);
+
+ for(;;) {
+ struct filter_context ctx = FILTER_CONTEXT_NULL;
+
+ /* Sample an optical thickness */
+ ctx.Ts = ssp_ran_exp(rng, 1);
+
+ /* Setup the remaining fields of the hit filter context */
+ ctx.rng = rng;
+ ctx.htsky = htrdr->sky;
+ ctx.iband = iband;
+ ctx.iquad = iquad;
+
+ /* Found the first intersection with the surface geometry */
+ HTRDR(ground_trace_ray
+ (htrdr->ground, pos, dir, range, &s3d_hit_prev, &s3d_hit));
+
+ /* Fit the ray range to the surface distance along the ray */
+ range[0] = 0;
+ range[1] = s3d_hit.distance;
+
+ /* Trace a ray into the participating media */
+ HTSKY(trace_ray(htrdr->sky, pos, dir, range, NULL,
+ hit_filter, &ctx, iband, iquad, &svx_hit));
+
+ /* No scattering and no surface reflection.
+ * Congratulation !! You are in space. */
+ if(S3D_HIT_NONE(&s3d_hit) && SVX_HIT_NONE(&svx_hit)) {
+ w = 0;
+ break;
+ }
+
+ /* Compute the next position */
+ pos_next[0] = pos[0] + dir[0]*ctx.traversal_dst;
+ pos_next[1] = pos[1] + dir[1]*ctx.traversal_dst;
+ pos_next[2] = pos[2] + dir[2]*ctx.traversal_dst;
+
+ /* Absorption event. Stop the realisation */
+ if(ctx.event_type == EVENT_ABSORPTION) {
+ ASSERT(!SVX_HIT_NONE(&svx_hit));
+ temperature = htsky_fetch_temperature(htrdr->sky, pos_next);
+ w = planck(band_bounds_m[0], band_bounds_m[1], temperature);
+ break;
+ }
+
+ /* Negate the incoming dir to match the convention of the SSF library */
+ d3_minus(wo, dir);
+
+ /* Scattering in the volume */
+ if(ctx.event_type == EVENT_SCATTERING) {
+ ASSERT(!SVX_HIT_NONE(&svx_hit));
+ ssf_phase_sample(phase_hg, rng, wo, dir_next, NULL);
+ s3d_hit_prev = S3D_HIT_NULL;
+
+ /* Scattering at a surface */
+ } else {
+ struct htrdr_interface interf = HTRDR_INTERFACE_NULL;
+ struct ssf_bsdf* bsdf = NULL;
+ double bounce_reflectivity = 0;
+ double N[3];
+ int type;
+ ASSERT(ctx.event_type == EVENT_NONE);
+ ASSERT(!S3D_HIT_NONE(&s3d_hit));
+
+ /* Fetch the hit interface and build its BSDF */
+ htrdr_ground_get_interface(htrdr->ground, &s3d_hit, &interf);
+ HTRDR(interface_create_bsdf
+ (htrdr, &interf, ithread, wlen, pos_next, dir, rng, &s3d_hit, &bsdf));
+
+ d3_normalize(N, d3_set_f3(N, s3d_hit.normal));
+ if(d3_dot(N, wo) < 0) d3_minus(N, N);
+
+ bounce_reflectivity = ssf_bsdf_sample
+ (bsdf, rng, wo, N, dir_next, &type, NULL);
+ if(!(type & SSF_REFLECTION)) { /* Handle only reflections */
+ bounce_reflectivity = 0;
+ }
+
+ /* Release the BSDF */
+ 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 {
+ w = planck(band_bounds_m[0], band_bounds_m[1], temperature);
+ }
+ break;
+ }
+
+ s3d_hit_prev = s3d_hit;
+ }
+ d3_set(pos, pos_next);
+ d3_set(dir, dir_next);
+ }
+ SSF(phase_ref_put(phase_hg));
+
+ return w;
+}
+
diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c
@@ -16,6 +16,7 @@
#include "htrdr.h"
#include "htrdr_c.h"
+#include "htrdr_interface.h"
#include "htrdr_ground.h"
#include "htrdr_solve.h"
#include "htrdr_sun.h"
@@ -34,7 +35,7 @@
struct scattering_context {
struct ssp_rng* rng;
const struct htsky* sky;
- size_t iband; /* Index of the spectrald */
+ size_t iband; /* Index of the spectral band */
size_t iquad; /* Index of the quadrature point into the band */
double Ts; /* Sampled optical thickness */
@@ -246,6 +247,7 @@ htrdr_compute_radiance_sw
struct ssp_rng* rng,
const double pos_in[3],
const double dir_in[3],
+ const double wlen, /* In nanometer */
const size_t iband,
const size_t iquad)
{
@@ -255,7 +257,6 @@ htrdr_compute_radiance_sw
struct svx_hit svx_hit = SVX_HIT_NULL;
struct ssf_phase* phase_hg = NULL;
struct ssf_phase* phase_rayleigh = NULL;
- struct ssf_bsdf* bsdf = NULL;
double pos[3];
double dir[3];
@@ -264,7 +265,6 @@ htrdr_compute_radiance_sw
double dir_next[3];
double band_bounds[2]; /* In nanometers */
- double wlen;
double R;
double r; /* Random number */
double wo[3]; /* -dir */
@@ -279,27 +279,23 @@ htrdr_compute_radiance_sw
double g = 0; /* Asymmetry parameter of the HG phase function */
ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads);
+ ASSERT(!htsky_is_long_wave(htrdr->sky));
CHK(RES_OK == ssf_phase_create
(&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
CHK(RES_OK == ssf_phase_create
(&htrdr->lifo_allocators[ithread], &ssf_phase_rayleigh, &phase_rayleigh));
- /* Setup the phase function for this spectral band & quadrature point */
- g = htsky_fetch_particle_phase_function_asymmetry_parameter
- (htrdr->sky, iband, iquad);
+ /* Setup the phase function for this wavelength */
+ g = htsky_fetch_per_wavelength_particle_phase_function_asymmetry_parameter
+ (htrdr->sky, wlen);
SSF(phase_hg_setup(phase_hg, g));
- /* Fetch the ground BSDF */
- bsdf = htrdr_ground_get_bsdf(htrdr->ground);
-
- /* Fetch sun properties. Arbitrarily use the wavelength at the center of the
- * band to retrieve the sun radiance of the current band. Note that the sun
- * spectral data are defined by bands that, actually are the same of the SW
- * spectral bands defined in the default "ecrad_opt_prot.txt" file provided
- * by the HTGOP project. */
- htsky_get_sw_spectral_band_bounds(htrdr->sky, iband, band_bounds);
- wlen = (band_bounds[0] + band_bounds[1]) * 0.5;
+ /* Fetch sun properties. Note that the sun spectral data are defined by bands
+ * that, actually are the same of the SW spectral bands defined in the
+ * default "ecrad_opt_prot.txt" file provided by the HTGOP project. */
+ htsky_get_spectral_band_bounds(htrdr->sky, iband, band_bounds);
+ ASSERT(band_bounds[0] <= wlen && wlen <= band_bounds[1]);
sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
@@ -389,9 +385,16 @@ htrdr_compute_radiance_sw
/* Scattering at a surface */
if(SVX_HIT_NONE(&svx_hit)) {
+ struct htrdr_interface interf = HTRDR_INTERFACE_NULL;
+ struct ssf_bsdf* bsdf = NULL;
double N[3];
int type;
+ /* Fetch the hit interface and build its BSDF */
+ htrdr_ground_get_interface(htrdr->ground, &s3d_hit, &interf);
+ HTRDR(interface_create_bsdf
+ (htrdr, &interf, ithread, wlen, pos_next, dir, rng, &s3d_hit, &bsdf));
+
d3_normalize(N, d3_set_f3(N, s3d_hit.normal));
if(d3_dot(N, wo) < 0) d3_minus(N, N);
@@ -407,6 +410,9 @@ htrdr_compute_radiance_sw
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;
@@ -414,7 +420,7 @@ htrdr_compute_radiance_sw
double ks_gas; /* Scattering coefficient of the gaz */
double ks; /* Overall scattering coefficient */
- ks_gas = htsky_fetch_raw_property(htrdr->sky, HTSKY_Ks,
+ ks_gas = htsky_fetch_raw_property(htrdr->sky, HTSKY_Ks,
HTSKY_CPNT_FLAG_GAS, iband, iquad, pos_next, -DBL_MAX, DBL_MAX);
ks_particle = htsky_fetch_raw_property(htrdr->sky, HTSKY_Ks,
HTSKY_CPNT_FLAG_PARTICLES, iband, iquad, pos_next, -DBL_MAX, DBL_MAX);
diff --git a/src/htrdr_draw_radiance.c b/src/htrdr_draw_radiance.c
@@ -0,0 +1,1039 @@
+/* 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_lw.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_LW,
+ PIXEL_SW
+};
+
+union pixel {
+ struct htrdr_pixel_lw lw;
+ struct htrdr_pixel_sw sw;
+};
+
+/* 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 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_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(buf && tile_data);
+
+ htrdr_buffer_get_layout(buf, &layout);
+ buf_mem = htrdr_buffer_get_data(buf);
+ ASSERT(layout.elmt_size == (tile_data->format == PIXEL_LW
+ ? sizeof(struct htrdr_pixel_lw)
+ : sizeof(struct htrdr_pixel_sw)));
+
+ /* 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_LW:
+ ((struct htrdr_pixel_lw*)buf_col)[x] = tile_row[x].lw;
+ break;
+ case PIXEL_SW:
+ ((struct htrdr_pixel_sw*)buf_col)[x] = tile_row[x].sw;
+ 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(buf, &t->data);
+ ++itile;
+ }
+
+ if(itile != ntiles) {
+ ASSERT(htrdr->mpi_nprocs > 1);
+
+ /* Create a temporary tile to receive the tile data computed by the
+ * concurrent MPI processes */
+ tile = tile_create
+ (htrdr->allocator, htsky_is_long_wave(htrdr->sky) ? PIXEL_LW : PIXEL_SW);
+ 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(buf, &tile->data);
+ }
+ }
+ }
+
+exit:
+ if(tile) tile_ref_put(tile);
+ return res;
+error:
+ goto exit;
+}
+
+static void
+draw_pixel_sw
+ (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_sw* pixel)
+{
+ struct htrdr_accum XYZ[3]; /* X, Y, and Z */
+ struct htrdr_accum time;
+ size_t ichannel;
+ ASSERT(ipix && ipix && pix_sz && cam && rng && pixel);
+ ASSERT(!htsky_is_long_wave(htrdr->sky));
+
+ /* 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 */
+ 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); break;
+ case 1: wlen = htrdr_cie_xyz_sample_Y(htrdr->cie, r0, r1); break;
+ case 2: wlen = htrdr_cie_xyz_sample_Z(htrdr->cie, r0, r1); break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ iband = htsky_find_spectral_band(htrdr->sky, wlen);
+ iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r2, iband);
+
+ /* Compute the luminance */
+ weight = htrdr_compute_radiance_sw
+ (htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
+ ASSERT(weight >= 0);
+
+ /* 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_lw
+ (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_lw* pixel)
+{
+ struct htrdr_accum radiance;
+ struct htrdr_accum time;
+ size_t isamp;
+ double temp_min, temp_max;
+ ASSERT(ipix && ipix && pix_sz && cam && rng && pixel);
+ ASSERT(htsky_is_long_wave(htrdr->sky));
+
+ /* 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;
+ 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);
+
+ /* Sample a wavelength */
+ wlen = htrdr_ran_lw_sample(htrdr->ran_lw, r0);
+
+ /* Select the associated band and sample a quadrature point */
+ iband = htsky_find_spectral_band(htrdr->sky, wlen);
+ iquad = htsky_spectral_band_sample_quadrature(htrdr->sky, r1, iband);
+
+ /* Retrieve the PDF to sample this sky band */
+ band_pdf = htrdr_ran_lw_get_sky_band_pdf(htrdr->ran_lw, iband);
+
+ /* Compute the luminance in W/m^2/sr/m */
+ weight = htrdr_compute_radiance_lw
+ (htrdr, ithread, rng, ray_org, ray_dir, wlen, iband, iquad);
+ weight /= band_pdf;
+ ASSERT(weight >= 0);
+
+ /* 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 */
+ #define BRIGHTNESS_TEMPERATURE(Radiance, Temperature) { \
+ res_T res = brightness_temperature \
+ (htrdr, \
+ htrdr->wlen_range_m[0], \
+ htrdr->wlen_range_m[1], \
+ (Radiance), \
+ &(Temperature)); \
+ if(res != RES_OK) { \
+ htrdr_log_warn(htrdr, \
+ "Could not compute the brightness temperature for the radiance %g.\n", \
+ (Radiance)); \
+ (Temperature) = 0; \
+ } \
+ } (void)0
+ BRIGHTNESS_TEMPERATURE(pixel->radiance.E, pixel->radiance_temperature.E);
+ BRIGHTNESS_TEMPERATURE(pixel->radiance.E - pixel->radiance.SE, temp_min);
+ BRIGHTNESS_TEMPERATURE(pixel->radiance.E + pixel->radiance.SE, temp_max);
+ pixel->radiance_temperature.SE = temp_max - temp_min;
+ #undef BRIGHTNESS_TEMPERATURE
+
+ /* Transform the pixel radiance from W/m^2/sr/m to W/m^/sr/nm */
+ pixel->radiance.E *= 1.0e-9;
+ pixel->radiance.SE *= 1.0e-9;
+}
+
+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 */
+ if(htsky_is_long_wave(htrdr->sky)) {
+ draw_pixel_lw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->lw);
+ } else {
+ draw_pixel_sw(htrdr, ithread, ipix, pix_sz, cam, spp, rng, &pixel->sw);
+ }
+ }
+ 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 = htsky_is_long_wave(htrdr->sky) ? PIXEL_LW : PIXEL_SW;
+
+ 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) {
+ size_t pixsz, pixal;
+ htrdr_buffer_get_layout(buf, &layout);
+ ASSERT(layout.width || layout.height || layout.elmt_size);
+ ASSERT(layout.width == width && layout.height == height);
+
+ if(htsky_is_long_wave(htrdr->sky)) {
+ pixsz = sizeof(struct htrdr_pixel_lw);
+ pixal = ALIGNOF(struct htrdr_pixel_lw);
+ } else {
+ pixsz = sizeof(struct htrdr_pixel_sw);
+ pixal = ALIGNOF(struct htrdr_pixel_sw);
+ }
+
+ if(layout.elmt_size != pixsz || layout.alignment < pixal) {
+ htrdr_log_err(htrdr, "%s: invalid buffer layout.\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ 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_draw_radiance_sw.c b/src/htrdr_draw_radiance_sw.c
@@ -1,851 +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 199309L /* nanosleep */
-
-#include "htrdr.h"
-#include "htrdr_c.h"
-#include "htrdr_buffer.h"
-#include "htrdr_camera.h"
-#include "htrdr_solve.h"
-
-#include <high_tune/htsky.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);
-
-/* 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 */
- /* Simulate the flexible array member of the C99 standard. */
- struct htrdr_accum accums[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 FINLINE struct tile*
-tile_create(struct mem_allocator* allocator)
-{
- struct tile* tile;
- const size_t tile_sz =
- sizeof(struct tile) - sizeof(struct htrdr_accum)/*rm dummy accum*/;
- const size_t buf_sz = /* Flexiblbe array element */
- TILE_SIZE*TILE_SIZE*sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]);
- ASSERT(allocator);
-
- tile = MEM_ALLOC(allocator, tile_sz+buf_sz);
- if(!tile) return NULL;
-
- ref_init(&tile->ref);
- list_init(&tile->node);
- tile->allocator = allocator;
- ASSERT(IS_ALIGNED(&tile->data.accums, ALIGNOF(struct htrdr_accum)));
-
- 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 struct htrdr_accum*
-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.accums + (y*TILE_SIZE + x) * HTRDR_ESTIMATES_COUNT__;
-}
-
-static void
-write_tile_data(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(buf && tile_data);
-
- htrdr_buffer_get_layout(buf, &layout);
- buf_mem = htrdr_buffer_get_data(buf);
- ASSERT(layout.elmt_size == sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]));
-
- /* 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 tile data, row by row */
- FOR_EACH(irow_tile, 0, nrows_tile) {
- char* buf_row = buf_mem + (irow + irow_tile) * layout.pitch;
- const struct htrdr_accum* tile_row =
- tile_data->accums + irow_tile*TILE_SIZE*HTRDR_ESTIMATES_COUNT__;
- memcpy(buf_row + icol*sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]),
- tile_row, ncols_tile*sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]));
- }
-}
-
-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(struct htrdr_accum)/*dummy*/
- + TILE_SIZE*TILE_SIZE*sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]);
-
- 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(buf, &t->data);
- ++itile;
- }
-
- if(itile != ntiles) {
- ASSERT(htrdr->mpi_nprocs > 1);
-
- /* Create a temporary tile to receive the tile data computed by the
- * concurrent MPI processes */
- tile = tile_create(htrdr->allocator);
- 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 concurret 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(buf, &tile->data);
- }
- }
- }
-
-exit:
- if(tile) tile_ref_put(tile);
- return res;
-error:
- goto exit;
-}
-
-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) {
- struct htrdr_accum* pix_accums;
- size_t ipix_tile[2]; /* Pixel coord in the tile */
- size_t ipix[2]; /* Pixel coord in the buffer */
- size_t ichannel;
-
- 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 */
- pix_accums = tile_at(tile, ipix_tile[0], ipix_tile[1]);
- pix_accums[HTRDR_ESTIMATE_TIME] = HTRDR_ACCUM_NULL; /* Reset time per radiative path */
-
- /* Compute the pixel coordinate */
- ipix[0] = tile_org[0] + ipix_tile[0];
- ipix[1] = tile_org[1] + ipix_tile[1];
-
- FOR_EACH(ichannel, 0, 3) {
- /* Check that the X, Y and Z estimates are stored in accumulators 0, 1 et
- * 2, respectively */
- STATIC_ASSERT
- ( HTRDR_ESTIMATE_X == 0
- && HTRDR_ESTIMATE_Y == 1
- && HTRDR_ESTIMATE_Z == 2,
- Unexpected_htrdr_estimate_enumerate);
- size_t isamp;
-
- pix_accums[ichannel] = 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;
- size_t iband;
- size_t iquad;
- double usec;
-
- /* 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);
-
- /* Sample a spectral band and a quadrature point */
- r0 = ssp_rng_canonical(rng);
- r1 = ssp_rng_canonical(rng);
- switch(ichannel) {
- case 0:
- htsky_sample_sw_spectral_data_CIE_1931_X
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
- case 1:
- htsky_sample_sw_spectral_data_CIE_1931_Y
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
- case 2:
- htsky_sample_sw_spectral_data_CIE_1931_Z
- (htrdr->sky, r0, r1, &iband, &iquad);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
-
- /* Compute the radiance that reach the pixel through the ray */
- time_current(&t0);
- weight = htrdr_compute_radiance_sw
- (htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
- ASSERT(weight >= 0);
- time_sub(&t0, time_current(&t1), &t0);
- usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
-
- /* Update the pixel accumulator of the current channel */
- pix_accums[ichannel].sum_weights += weight;
- pix_accums[ichannel].sum_weights_sqr += weight*weight;
- pix_accums[ichannel].nweights += 1;
-
- /* Update the pixel accumulator of per realisation time */
- pix_accums[HTRDR_ESTIMATE_TIME].sum_weights += usec;
- pix_accums[HTRDR_ESTIMATE_TIME].sum_weights_sqr += usec*usec;
- pix_accums[HTRDR_ESTIMATE_TIME].nweights += 1;
- }
- }
- }
- 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;
- 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;
-
- 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);
- 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 a 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;
-
- /* Synchronize the process */
- mutex_lock(htrdr->mpi_mutex);
- MPI(Barrier(MPI_COMM_WORLD));
- mutex_unlock(htrdr->mpi_mutex);
-
-exit:
- if(rng_proc) SSP(rng_ref_put(rng_proc));
- return (res_T)res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_draw_radiance_sw
- (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) {
- 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 != sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__])
- || layout.alignment < ALIGNOF(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__])) {
- htrdr_log_err(htrdr,
- "%s: invalid buffer layout. "
- "The pixel size must be the size of %lu accumulators.\n",
- FUNC_NAME, (unsigned long)HTRDR_ESTIMATES_COUNT__);
- 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
@@ -1,5 +1,5 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+/* Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ * Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -15,19 +15,42 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "htrdr.h"
+#include "htrdr_interface.h"
#include "htrdr_ground.h"
+#include "htrdr_mtl.h"
#include "htrdr_slab.h"
+#include <aw.h>
#include <rsys/clock_time.h>
#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/dynamic_array_size_t.h>
#include <rsys/double2.h>
#include <rsys/double3.h>
#include <rsys/float2.h>
#include <rsys/float3.h>
+#include <rsys/hash_table.h>
#include <star/s3d.h>
-#include <star/s3daw.h>
-#include <star/ssf.h>
+
+/* Define the hash table that maps an Obj vertex id to its position into the
+ * vertex buffer */
+#define HTABLE_NAME vertex
+#define HTABLE_KEY size_t /* Obj vertex id */
+#define HTABLE_DATA size_t
+#include <rsys/hash_table.h>
+
+/* Define the hash table that maps the Star-3D shape id to its interface */
+#define HTABLE_NAME interface
+#define HTABLE_KEY unsigned /* Star-3D shape id */
+#define HTABLE_DATA struct htrdr_interface
+#include <rsys/hash_table.h>
+
+struct mesh {
+ const struct darray_double* positions;
+ const struct darray_size_t* indices;
+};
+static const struct mesh MESH_NULL;
struct ray_context {
float range[2];
@@ -50,9 +73,10 @@ struct htrdr_ground {
struct s3d_scene_view* view;
float lower[3]; /* Ground lower bound */
float upper[3]; /* Ground upper bound */
- struct ssf_bsdf* bsdf;
int repeat; /* Make the ground infinite in X and Y */
+ struct htable_interface interfaces; /* Map a Star3D shape to its interface */
+
struct htrdr* htrdr;
ref_T ref;
};
@@ -131,146 +155,341 @@ trace_ground
}
static res_T
-setup_ground(struct htrdr_ground* ground, const char* obj_filename)
+parse_shape_interface
+ (struct htrdr* htrdr,
+ const char* name,
+ struct htrdr_interface* interf)
{
- struct s3d_scene* scn = NULL;
- struct s3daw* s3daw = NULL;
- size_t nshapes;
- size_t ishape;
+ struct str str;
+ char* mtl_name_front = NULL;
+ char* mtl_name_back = NULL;
+ char* tk_ctx = NULL;
res_T res = RES_OK;
- ASSERT(ground);
+ ASSERT(htrdr && name && interf);
- if(!obj_filename) {
- /* No geometry! Discard the creation of the scene view */
- htrdr_log_warn(ground->htrdr,
- "%s: the scene does not have ground geometry.\n",
- FUNC_NAME);
- goto exit;
- }
+ str_init(htrdr->allocator, &str);
- res = s3daw_create(&ground->htrdr->logger, ground->htrdr->allocator, NULL,
- NULL, ground->htrdr->s3d, ground->htrdr->verbose, &s3daw);
+ /* Locally copy the string to parse */
+ res = str_set(&str, name);
if(res != RES_OK) {
- htrdr_log_err(ground->htrdr,
- "%s: could not create the Star-3DAW device -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
+ htrdr_log_err(htrdr,
+ "Could not locally copy the shape material string `%s' -- %s.\n",
+ name, res_to_cstr(res));
goto error;
}
- res = s3daw_load(s3daw, obj_filename);
- if(res != RES_OK) {
- htrdr_log_err(ground->htrdr,
- "%s: could not load the obj file `%s' -- %s.\n",
- FUNC_NAME, obj_filename, res_to_cstr(res));
+ /* 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);
+ if(!mtl_name_back) {
+ htrdr_log_err(htrdr,
+ "The material name of the shape back faces are missing `%s'.\n", name);
+ res = RES_BAD_ARG;
goto error;
}
- res = s3d_scene_create(ground->htrdr->s3d, &scn);
- if(res != RES_OK) {
- htrdr_log_err(ground->htrdr,
- "%s: could not create the Star-3D scene of the ground -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
+ /* 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) {
+ htrdr_log_err(htrdr,
+ "Invalid interface `%s:%s'. "
+ "The front and the back materials are both uknown.\n",
+ mtl_name_front, mtl_name_back);
+ res = RES_BAD_ARG;
goto error;
}
+exit:
+ str_release(&str);
+ return res;
+error:
+ interf->mtl_front = interf->mtl_back = NULL;
+ goto exit;
+}
- S3DAW(get_shapes_count(s3daw, &nshapes));
- FOR_EACH(ishape, 0, nshapes) {
- struct s3d_shape* shape;
- S3DAW(get_shape(s3daw, ishape, &shape));
- res = s3d_mesh_set_hit_filter_function(shape, ground_filter, NULL);
- if(res != RES_OK) {
- htrdr_log_err(ground->htrdr,
- "%s: could not setup the hit filter function of the ground geometry "
- "-- %s.\n", FUNC_NAME, res_to_cstr(res));
+static res_T
+setup_mesh
+ (struct htrdr* htrdr,
+ const char* filename,
+ struct aw_obj* obj,
+ struct aw_obj_named_group* mtl,
+ struct darray_double* positions,
+ struct darray_size_t* indices,
+ struct htable_vertex* vertices) /* Scratch data structure */
+{
+ size_t iface;
+ res_T res = RES_OK;
+ ASSERT(htrdr && filename && obj && mtl && positions && indices && vertices);
+
+ darray_double_clear(positions);
+ darray_size_t_clear(indices);
+ htable_vertex_clear(vertices);
+
+ FOR_EACH(iface, mtl->face_id, mtl->face_id+mtl->faces_count) {
+ struct aw_obj_face face;
+ size_t ivertex;
+
+ AW(obj_get_face(obj, iface, &face));
+ if(face.vertices_count != 3) {
+ htrdr_log_err(htrdr,
+ "The obj `%s' has non-triangulated polygons "
+ "while currently only triangular meshes are supported.\n",
+ filename);
+ res = RES_BAD_ARG;
goto error;
}
- res = s3d_scene_attach_shape(scn, shape);
- if(res != RES_OK) {
- htrdr_log_err(ground->htrdr,
- "%s: could not attach the ground geometry to its Star-3D scene -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
+
+ FOR_EACH(ivertex, 0, face.vertices_count) {
+ struct aw_obj_vertex vertex;
+ size_t id;
+ size_t* pid;
+
+ AW(obj_get_vertex(obj, face.vertex_id + ivertex, &vertex));
+ pid = htable_vertex_find(vertices, &vertex.position_id);
+ if(pid) {
+ id = *pid;
+ } else {
+ struct aw_obj_vertex_data vdata;
+
+ id = darray_double_size_get(positions) / 3;
+
+ res = darray_double_resize(positions, id*3 + 3);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr,
+ "Could not locally copy the vertex position -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ AW(obj_get_vertex_data(obj, &vertex, &vdata));
+ darray_double_data_get(positions)[id*3+0] = vdata.position[0];
+ darray_double_data_get(positions)[id*3+1] = vdata.position[1];
+ darray_double_data_get(positions)[id*3+2] = vdata.position[2];
+
+ res = htable_vertex_set(vertices, &vertex.position_id, &id);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr,
+ "Could not register the vertex position -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+ }
+
+ res = darray_size_t_push_back(indices, &id);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr,
+ "Could not locally copy the face index -- %s\n",
+ res_to_cstr(res));
+ goto error;
+ }
}
}
+exit:
+ return res;
+error:
+ darray_double_clear(positions);
+ darray_size_t_clear(indices);
+ htable_vertex_clear(vertices);
+ goto exit;
+}
+
+static void
+get_position(const unsigned ivert, float position[3], void* ctx)
+{
+ const struct mesh* mesh = ctx;
+ const double* pos = NULL;
+ ASSERT(mesh);
+ ASSERT(ivert < darray_double_size_get(mesh->positions) / 3);
+ pos = darray_double_cdata_get(mesh->positions) + ivert*3;
+ position[0] = (float)pos[0];
+ position[1] = (float)pos[1];
+ position[2] = (float)pos[2];
+}
+
+static void
+get_indices(const unsigned itri, unsigned indices[3], void* ctx)
+{
+ const struct mesh* mesh = ctx;
+ const size_t* ids = NULL;
+ ASSERT(mesh);
+ ASSERT(itri < darray_size_t_size_get(mesh->indices) / 3);
+ ids = darray_size_t_cdata_get(mesh->indices) + itri*3;
+ indices[0] = (unsigned)ids[0];
+ indices[1] = (unsigned)ids[1];
+ indices[2] = (unsigned)ids[2];
+}
+
+static res_T
+create_s3d_shape
+ (struct htrdr* htrdr,
+ const struct darray_double* positions,
+ const struct darray_size_t* indices,
+ struct s3d_shape** out_shape)
+{
+ struct s3d_shape* shape = NULL;
+ struct s3d_vertex_data vdata = S3D_VERTEX_DATA_NULL;
+ struct mesh mesh = MESH_NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && positions && indices && out_shape);
+ ASSERT(darray_double_size_get(positions) != 0);
+ ASSERT(darray_size_t_size_get(indices) != 0);
+ ASSERT(darray_double_size_get(positions)%3 == 0);
+ ASSERT(darray_size_t_size_get(indices)%3 == 0);
- res = s3d_scene_view_create(scn, S3D_TRACE, &ground->view);
+ mesh.positions = positions;
+ mesh.indices = indices;
+
+ res = s3d_shape_create_mesh(htrdr->s3d, &shape);
if(res != RES_OK) {
- htrdr_log_err(ground->htrdr,
- "%s: could not create the Star-3D scene view of the ground geometry "
- "-- %s.\n", FUNC_NAME, res_to_cstr(res));
+ htrdr_log_err(htrdr, "Error creating a Star-3D shape -- %s.\n",
+ res_to_cstr(res));
goto error;
}
- res = s3d_scene_view_get_aabb(ground->view, ground->lower, ground->upper);
+ vdata.usage = S3D_POSITION;
+ vdata.type = S3D_FLOAT3;
+ vdata.get = get_position;
+
+ res = s3d_mesh_setup_indexed_vertices
+ (shape, (unsigned int)(darray_size_t_size_get(indices)/3), get_indices,
+ (unsigned int)(darray_double_size_get(positions)/3), &vdata, 1, &mesh);
+ if(res != RES_OK){
+ htrdr_log_err(htrdr, "Could not setup the Star-3D shape -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ res = s3d_mesh_set_hit_filter_function(shape, ground_filter, NULL);
if(res != RES_OK) {
- htrdr_log_err(ground->htrdr,
- "%s: could not get the ground bounding box -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
+ htrdr_log_err(htrdr,
+ "Could not setup the Star-3D hit filter function of the ground geometry "
+ "-- %s.\n", res_to_cstr(res));
goto error;
}
exit:
- if(s3daw) S3DAW(ref_put(s3daw));
- if(scn) S3D(scene_ref_put(scn));
+ *out_shape = shape;
return res;
error:
+ if(shape) {
+ S3D(shape_ref_put(shape));
+ shape = NULL;
+ }
goto exit;
}
static res_T
-setup_bsdf_diffuse(struct htrdr_ground* ground, const double reflectivity)
+setup_ground(struct htrdr_ground* ground, const char* obj_filename)
{
+ struct aw_obj_desc desc;
+ struct htable_vertex vertices;
+ struct darray_double positions;
+ struct darray_size_t indices;
+ struct aw_obj* obj = NULL;
+ struct s3d_shape* shape = NULL;
+ struct s3d_scene* scene = NULL;
+ size_t iusemtl;
+
res_T res = RES_OK;
- ASSERT(ground);
+ ASSERT(obj_filename);
- res = ssf_bsdf_create
- (ground->htrdr->allocator, &ssf_lambertian_reflection, &ground->bsdf);
- if(res != RES_OK) goto error;
+ htable_vertex_init(ground->htrdr->allocator, &vertices);
+ darray_double_init(ground->htrdr->allocator, &positions);
+ darray_size_t_init(ground->htrdr->allocator, &indices);
+
+ res = aw_obj_create(&ground->htrdr->logger, ground->htrdr->allocator,
+ ground->htrdr->verbose, &obj);
+ if(res != RES_OK) {
+ htrdr_log_err(ground->htrdr, "Could not create the obj loader -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ res = s3d_scene_create(ground->htrdr->s3d, &scene);
+ if(res != RES_OK) {
+ htrdr_log_err(ground->htrdr, "Could not create the Star-3D scene -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
- res = ssf_lambertian_reflection_setup(ground->bsdf, reflectivity);
+ /* Load the geometry data */
+ res = aw_obj_load(obj, obj_filename);
if(res != RES_OK) goto error;
-exit:
- return res;
-error:
- htrdr_log_err(ground->htrdr,
- "Could not setup the ground diffuse BSDF with a reflectivity of %g -- %s.\n",
- reflectivity, res_to_cstr(res));
- if(ground->bsdf) {
- SSF(bsdf_ref_put(ground->bsdf));
- ground->bsdf = NULL;
+ /* Fetch the descriptor of the loaded geometry */
+ AW(obj_get_desc(obj, &desc));
+
+ if(desc.usemtls_count == 0) {
+ htrdr_log_err(ground->htrdr, "The obj `%s' has no material.\n", obj_filename);
+ res = RES_BAD_ARG;
+ goto error;
}
- goto exit;
-}
-static res_T
-setup_bsdf_specular(struct htrdr_ground* ground, const double reflectivity)
-{
- struct ssf_fresnel* fresnel = NULL;
- res_T res = RES_OK;
- ASSERT(ground);
+ /* Setup the geometry */
+ FOR_EACH(iusemtl, 0, desc.usemtls_count) {
+ struct aw_obj_named_group mtl;
+ struct htrdr_interface interf;
+ unsigned ishape;
- res = ssf_bsdf_create
- (ground->htrdr->allocator, &ssf_specular_reflection, &ground->bsdf);
- if(res != RES_OK) goto error;
- res = ssf_fresnel_create
- (ground->htrdr->allocator, &ssf_fresnel_constant, &fresnel);
- if(res != RES_OK) goto error;
- res = ssf_fresnel_constant_setup(fresnel, reflectivity);
- if(res != RES_OK) goto error;
- res = ssf_specular_reflection_setup(ground->bsdf, fresnel);
- if(res != RES_OK) goto error;
+ AW(obj_get_mtl(obj, iusemtl , &mtl));
+
+ res = parse_shape_interface(ground->htrdr, mtl.name, &interf);
+ if(res != RES_OK) goto error;
+
+ res = setup_mesh
+ (ground->htrdr, obj_filename, obj, &mtl, &positions, &indices, &vertices);
+ if(res != RES_OK) goto error;
+
+ res = create_s3d_shape(ground->htrdr, &positions, &indices, &shape);
+ if(res != RES_OK) goto error;
+
+ S3D(shape_get_id(shape, &ishape));
+ res = htable_interface_set(&ground->interfaces, &ishape, &interf);
+ if(res != RES_OK) {
+ htrdr_log_err(ground->htrdr,
+ "Could not map the Star-3D shape to its Star-Materials -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ res = s3d_scene_attach_shape(scene, shape);
+ if(res != RES_OK) {
+ htrdr_log_err(ground->htrdr,
+ "Could not attach a Star-3D shape to the Star-3D scene -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ S3D(shape_ref_put(shape));
+ shape = NULL;
+ }
+
+ res = s3d_scene_view_create(scene, S3D_GET_PRIMITIVE|S3D_TRACE, &ground->view);
+ if(res != RES_OK) {
+ htrdr_log_err(ground->htrdr,
+ "Could not create the Star-3D scene view -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ res = s3d_scene_view_get_aabb(ground->view, ground->lower, ground->upper);
+ if(res != RES_OK) {
+ htrdr_log_err(ground->htrdr,
+ "Could not get the bounding box of the geometry -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
exit:
+ if(obj) AW(obj_ref_put(obj));
+ if(shape) S3D(shape_ref_put(shape));
+ if(scene) S3D(scene_ref_put(scene));
+ htable_vertex_release(&vertices);
+ darray_double_release(&positions);
+ darray_size_t_release(&indices);
return res;
error:
- htrdr_log_err(ground->htrdr,
- "Could not setup the ground specular BSDF with a reflectivity of %g -- %s.\n",
- reflectivity, res_to_cstr(res));
- if(ground->bsdf) {
- SSF(bsdf_ref_put(ground->bsdf));
- ground->bsdf = NULL;
- }
goto exit;
}
@@ -281,7 +500,7 @@ release_ground(ref_T* ref)
ASSERT(ref);
ground = CONTAINER_OF(ref, struct htrdr_ground, ref);
if(ground->view) S3D(scene_view_ref_put(ground->view));
- if(ground->bsdf) SSF(bsdf_ref_put(ground->bsdf));
+ htable_interface_release(&ground->interfaces);
MEM_RM(ground->htrdr->allocator, ground);
}
@@ -292,8 +511,6 @@ res_T
htrdr_ground_create
(struct htrdr* htrdr,
const char* obj_filename, /* May be NULL */
- const enum htrdr_bsdf_type bsdf_type,
- const double reflectivity,
const int repeat_ground, /* Infinitely repeat the ground in X and Y */
struct htrdr_ground** out_ground)
{
@@ -302,7 +519,6 @@ htrdr_ground_create
struct time t0, t1;
res_T res = RES_OK;
ASSERT(htrdr && out_ground);
- ASSERT(reflectivity >= 0 || reflectivity <= 1);
ground = MEM_CALLOC(htrdr->allocator, 1, sizeof(*ground));
if(!ground) {
@@ -317,18 +533,9 @@ htrdr_ground_create
ground->repeat = repeat_ground;
f3_splat(ground->lower, (float)INF);
f3_splat(ground->upper,-(float)INF);
+ htable_interface_init(ground->htrdr->allocator, &ground->interfaces);
- switch(bsdf_type) {
- case HTRDR_BSDF_DIFFUSE:
- res = setup_bsdf_diffuse(ground, reflectivity);
- break;
- case HTRDR_BSDF_SPECULAR:
- res = setup_bsdf_specular(ground, reflectivity);
- break;
- default: FATAL("Unreachable code\n");
-
- }
- if(res != RES_OK) goto error;
+ if(!obj_filename) goto exit;
time_current(&t0);
res = setup_ground(ground, obj_filename);
@@ -362,11 +569,19 @@ htrdr_ground_ref_put(struct htrdr_ground* ground)
ref_put(&ground->ref, release_ground);
}
-struct ssf_bsdf*
-htrdr_ground_get_bsdf(const struct htrdr_ground* ground)
+void
+htrdr_ground_get_interface
+ (struct htrdr_ground* ground,
+ const struct s3d_hit* hit,
+ struct htrdr_interface* out_interface)
{
- ASSERT(ground);
- return ground->bsdf;
+ struct htrdr_interface* interf = NULL;
+ ASSERT(ground && hit && out_interface);
+
+ interf = htable_interface_find(&ground->interfaces, &hit->prim.geom_id);
+ ASSERT(interf);
+
+ *out_interface = *interf;
}
res_T
diff --git a/src/htrdr_ground.h b/src/htrdr_ground.h
@@ -1,5 +1,5 @@
-/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+/* Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ * Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -19,14 +19,10 @@
#include <rsys/rsys.h>
-enum htrdr_bsdf_type {
- HTRDR_BSDF_DIFFUSE,
- HTRDR_BSDF_SPECULAR
-};
-
/* Forward declarations */
struct htrdr;
struct htrdr_ground;
+struct htrdr_interface;
struct s3d_hit;
struct ssf_bsdf;
@@ -34,8 +30,6 @@ extern LOCAL_SYM res_T
htrdr_ground_create
(struct htrdr* htrdr,
const char* obj_filename, /* May be NULL <=> No ground geometry */
- const enum htrdr_bsdf_type bsdf_type,
- const double reflectivity, /* In [0, 1] */
const int repeat_ground, /* Infinitely repeat the ground in X and Y */
struct htrdr_ground** ground);
@@ -47,9 +41,22 @@ extern LOCAL_SYM void
htrdr_ground_ref_put
(struct htrdr_ground* ground);
-extern LOCAL_SYM struct ssf_bsdf*
-htrdr_ground_get_bsdf
- (const struct htrdr_ground* ground);
+extern LOCAL_SYM void
+htrdr_ground_get_interface
+ (struct htrdr_ground* ground,
+ const struct s3d_hit* hit,
+ struct htrdr_interface* interface);
+
+extern LOCAL_SYM res_T
+htrdr_ground_create_bsdf
+ (struct htrdr_ground* ground,
+ const size_t ithread,
+ const double wavelength,
+ const double pos[3],
+ const double dir[3], /* Incoming ray */
+ const struct s3d_hit* hit,
+ struct htrdr_interface* interf, /* NULL <=> do not return the interface */
+ struct ssf_bsdf** bsdf);
extern LOCAL_SYM res_T
htrdr_ground_trace_ray
diff --git a/src/htrdr_interface.c b/src/htrdr_interface.c
@@ -0,0 +1,187 @@
+/* Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ * Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#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
@@ -0,0 +1,47 @@
+/* Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ * Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_INTERFACE_H
+#define HTRDR_INTERFACE_H
+
+#include <star/ssf.h>
+
+/* Forward declaration of external data type */
+struct mrumtl;
+struct s3d_hit;
+struct ssf_bsdf;
+struct ssp_rng;
+
+struct htrdr_interface {
+ const struct mrumtl* mtl_front;
+ const struct mrumtl* mtl_back;
+};
+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);
+
+#endif /* HTRDR_INTERFACE_H */
+
diff --git a/src/htrdr_mtl.c b/src/htrdr_mtl.c
@@ -0,0 +1,281 @@
+/* Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ * Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* 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
@@ -0,0 +1,44 @@
+/* Copyright (C) 2018, 2019 CNRS, Université Paul Sabatier
+ * Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_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_ran_lw.c b/src/htrdr_ran_lw.c
@@ -0,0 +1,417 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "htrdr.h"
+#include "htrdr_c.h"
+#include "htrdr_ran_lw.h"
+
+#include <high_tune/htsky.h>
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <math.h> /* nextafter */
+
+struct htrdr_ran_lw {
+ struct darray_double cdf;
+ /* Probas to sample a sky band overlapped by the range */
+ struct darray_double sky_bands_pdf;
+ double range[2]; /* Boundaries of the handled CIE XYZ color space */
+ double band_len; /* Length in nanometers of a band */
+ double ref_temperature; /* In Kelvin */
+ size_t nbands; /* # bands */
+
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+setup_ran_lw_cdf
+ (struct htrdr_ran_lw* ran_lw,
+ const char* func_name)
+{
+ double* pdf = NULL;
+ double* cdf = NULL;
+ double sum = 0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(ran_lw && func_name && ran_lw->nbands != HTRDR_RAN_LW_CONTINUE);
+
+ res = darray_double_resize(&ran_lw->cdf, ran_lw->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(ran_lw->htrdr,
+ "%s: Error allocating the CDF of the long wave spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+
+ cdf = darray_double_data_get(&ran_lw->cdf);
+ pdf = cdf; /* Alias the CDF by the PDF since only one array is necessary */
+
+ htrdr_log(ran_lw->htrdr,
+ "Number of bands of the long wave cumulative: %lu\n",
+ (unsigned long)ran_lw->nbands);
+
+ /* Compute the *unormalized* probability to sample a long wave band */
+ FOR_EACH(i, 0, ran_lw->nbands) {
+ double lambda_lo = ran_lw->range[0] + (double)i * ran_lw->band_len;
+ double lambda_hi = MMIN(lambda_lo + ran_lw->band_len, ran_lw->range[1]);
+ ASSERT(lambda_lo<= lambda_hi);
+ ASSERT(lambda_lo > ran_lw->range[0] || eq_eps(lambda_lo, ran_lw->range[0], 1.e-6));
+ ASSERT(lambda_lo < ran_lw->range[1] || eq_eps(lambda_lo, ran_lw->range[1], 1.e-6));
+
+ /* Convert from nanometer to meter */
+ lambda_lo *= 1.e-9;
+ lambda_hi *= 1.e-9;
+
+ /* Compute the probability of the current band */
+ pdf[i] = planck(lambda_lo, lambda_hi, ran_lw->ref_temperature);
+
+ /* Update the norm */
+ sum += pdf[i];
+ }
+
+ /* Compute the cumulative of the previously computed probabilities */
+ FOR_EACH(i, 0, ran_lw->nbands) {
+ /* Normalize the probability */
+ pdf[i] /= sum;
+
+ /* Setup the cumulative */
+ if(i == 0) {
+ cdf[i] = pdf[i];
+ } else {
+ cdf[i] = pdf[i] + cdf[i-1];
+ ASSERT(cdf[i] >= cdf[i-1]);
+ }
+ }
+ ASSERT(eq_eps(cdf[ran_lw->nbands-1], 1, 1.e-6));
+ cdf[ran_lw->nbands - 1] = 1.0; /* Handle numerical issue */
+
+exit:
+ return res;
+error:
+ darray_double_clear(&ran_lw->cdf);
+ goto exit;
+}
+
+static res_T
+compute_sky_bands_pdf
+ (struct htrdr_ran_lw* ran_lw,
+ const char* func_name,
+ double* out_min_band_len)
+{
+ double* pdf = NULL;
+ double min_band_len = DBL_MAX;
+ size_t nbands;
+ res_T res = RES_OK;
+ ASSERT(ran_lw && htsky_is_long_wave(ran_lw->htrdr->sky) && func_name);
+
+ nbands = htsky_get_spectral_bands_count(ran_lw->htrdr->sky);
+
+ res = darray_double_resize(&ran_lw->sky_bands_pdf, nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(ran_lw->htrdr,
+ "%s: error allocating the PDF of the long wave spectral bands "
+ "of the sky -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+
+ pdf = darray_double_data_get(&ran_lw->sky_bands_pdf);
+
+ if(eq_eps(ran_lw->range[0], ran_lw->range[1], 1.e-6)) {
+ ASSERT(nbands == 1);
+ pdf[0] = 1;
+ min_band_len = 0;
+ } else {
+ double sum = 0;
+ size_t i = 0;
+
+ /* Compute the *unormalized* probability to sample a long wave band */
+ FOR_EACH(i, 0, nbands) {
+ const size_t iband = htsky_get_spectral_band_id(ran_lw->htrdr->sky, i);
+ double wlens[2];
+ HTSKY(get_spectral_band_bounds(ran_lw->htrdr->sky, iband, wlens));
+
+ /* Adjust band boundaries to the submitted range */
+ wlens[0] = MMAX(wlens[0], ran_lw->range[0]);
+ wlens[1] = MMIN(wlens[1], ran_lw->range[1]);
+
+ min_band_len = MMIN(wlens[1] - wlens[0], min_band_len);
+
+ /* Convert from nanometer to meter */
+ wlens[0] *= 1.e-9;
+ wlens[1] *= 1.e-9;
+
+ /* Compute the probability of the current band */
+ pdf[i] = blackbody_fraction(wlens[0], wlens[1], ran_lw->ref_temperature);
+
+ /* Update the norm */
+ sum += pdf[i];
+ }
+
+ /* Normalize the probabilities */
+ FOR_EACH(i, 0, nbands) pdf[i] /= sum;
+ }
+
+exit:
+ *out_min_band_len = min_band_len;
+ return res;
+error:
+ darray_double_clear(&ran_lw->sky_bands_pdf);
+ goto exit;
+}
+
+static double
+ran_lw_sample_discreet
+ (const struct htrdr_ran_lw* ran_lw,
+ const double r,
+ const char* func_name)
+{
+ const double* cdf = NULL;
+ const double* find = NULL;
+ double r_next = nextafter(r, DBL_MAX);
+ double wlen = 0;
+ size_t cdf_length = 0;
+ size_t i;
+ ASSERT(ran_lw && ran_lw->nbands != HTRDR_RAN_LW_CONTINUE);
+ ASSERT(0 <= r && r < 1);
+ (void)func_name;
+
+ cdf = darray_double_cdata_get(&ran_lw->cdf);
+ cdf_length = darray_double_size_get(&ran_lw->cdf);
+ ASSERT(cdf_length > 0);
+
+ /* Use r_next rather than r in order to find the first entry that is not less
+ * than *or equal* to r */
+ find = search_lower_bound(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ ASSERT(find);
+
+ i = (size_t)(find - cdf);
+ ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
+
+ /* Return the wavelength at the center of the sampled band */
+ wlen = ran_lw->range[0] + ran_lw->band_len * ((double)i + 0.5);
+ ASSERT(ran_lw->range[0] < wlen && wlen < ran_lw->range[1]);
+
+ return wlen;
+}
+
+static double
+ran_lw_sample_continue
+ (const struct htrdr_ran_lw* ran_lw,
+ const double r,
+ const char* func_name)
+{
+ /* Numerical parameters of the dichotomy search */
+ const size_t MAX_ITER = 100;
+ const double EPSILON_LAMBDA_M = 1e-15;
+ const double EPSILON_BF = 1e-6;
+
+ /* Local variables */
+ double bf = 0;
+ double bf_prev = 0;
+ double bf_min_max = 0;
+ double lambda_m = 0;
+ double lambda_m_prev = 0;
+ double lambda_m_min = 0;
+ double lambda_m_max = 0;
+ double range_m[2] = {0, 0};
+ size_t i;
+
+ /* Check precondition */
+ ASSERT(ran_lw && ran_lw->nbands == HTRDR_RAN_LW_CONTINUE && func_name);
+ ASSERT(0 <= r && r < 1);
+
+ /* Convert the wavelength range in meters */
+ range_m[0] = ran_lw->range[0] * 1.0e-9;
+ range_m[1] = ran_lw->range[1] * 1.0e-9;
+
+ /* Setup the dichotomy search */
+ lambda_m_min = range_m[0];
+ lambda_m_max = range_m[1];
+ bf_min_max = blackbody_fraction
+ (range_m[0], range_m[1], ran_lw->ref_temperature);
+
+ /* Numerically search the lambda corresponding to the submitted canonical
+ * number */
+ FOR_EACH(i, 0, MAX_ITER) {
+ double r_test;
+ lambda_m = (lambda_m_min + lambda_m_max) * 0.5;
+ bf = blackbody_fraction(range_m[0], lambda_m, ran_lw->ref_temperature);
+
+ r_test = bf / bf_min_max;
+ if(r_test < r) {
+ lambda_m_min = lambda_m;
+ } else {
+ lambda_m_max = lambda_m;
+ }
+
+ if(fabs(lambda_m_prev - lambda_m) < EPSILON_LAMBDA_M
+ || fabs(bf_prev - bf) < EPSILON_BF)
+ break;
+
+ lambda_m_prev = lambda_m;
+ bf_prev = bf;
+ }
+ if(i >= MAX_ITER) {
+ htrdr_log_warn(ran_lw->htrdr,
+ "%s: could not sample a wavelength in the range [%g, %g] nanometers "
+ "for the reference temperature %g Kelvin.\n",
+ func_name, SPLIT2(ran_lw->range), ran_lw->ref_temperature);
+ }
+
+ return lambda_m * 1.0e9; /* Convert in nanometers */
+}
+
+static void
+release_ran_lw(ref_T* ref)
+{
+ struct htrdr_ran_lw* ran_lw = NULL;
+ ASSERT(ref);
+ ran_lw = CONTAINER_OF(ref, struct htrdr_ran_lw, ref);
+ darray_double_release(&ran_lw->cdf);
+ darray_double_release(&ran_lw->sky_bands_pdf);
+ MEM_RM(ran_lw->htrdr->allocator, ran_lw);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_ran_lw_create
+ (struct htrdr* htrdr,
+ const double range[2], /* Must be included in [1000, 100000] nanometers */
+ const size_t nbands, /* # bands used to discretized CDF */
+ const double ref_temperature,
+ struct htrdr_ran_lw** out_ran_lw)
+{
+ struct htrdr_ran_lw* ran_lw = NULL;
+ double min_band_len = 0;
+ res_T res = RES_OK;
+ ASSERT(htrdr && range && out_ran_lw && ref_temperature > 0);
+ ASSERT(ref_temperature > 0);
+ ASSERT(range[0] >= 1000);
+ ASSERT(range[1] <= 100000);
+ ASSERT(range[0] <= range[1]);
+
+ ran_lw = MEM_CALLOC(htrdr->allocator, 1, sizeof(*ran_lw));
+ if(!ran_lw) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: could not allocate long wave random variate data structure -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ ref_init(&ran_lw->ref);
+ ran_lw->htrdr = htrdr;
+ darray_double_init(htrdr->allocator, &ran_lw->cdf);
+ darray_double_init(htrdr->allocator, &ran_lw->sky_bands_pdf);
+
+ ran_lw->range[0] = range[0];
+ ran_lw->range[1] = range[1];
+ ran_lw->ref_temperature = ref_temperature;
+ ran_lw->nbands = nbands;
+
+ res = compute_sky_bands_pdf(ran_lw, FUNC_NAME, &min_band_len);
+ if(res != RES_OK) goto error;
+
+ if(nbands == HTRDR_RAN_LW_CONTINUE) {
+ ran_lw->band_len = 0;
+ } else {
+ ran_lw->band_len = (range[1] - range[0]) / (double)ran_lw->nbands;
+
+ /* Adjust the band length to ensure that each sky spectral interval is
+ * overlapped by at least one band */
+ if(ran_lw->band_len > min_band_len) {
+ ran_lw->band_len = min_band_len;
+ ran_lw->nbands = (size_t)ceil((range[1] - range[0]) / ran_lw->band_len);
+ }
+
+ res = setup_ran_lw_cdf(ran_lw, FUNC_NAME);
+ if(res != RES_OK) goto error;
+ }
+
+ htrdr_log(htrdr, "Long wave interval defined on [%g, %g] nanometers.\n",
+ range[0], range[1]);
+
+exit:
+ *out_ran_lw = ran_lw;
+ return res;
+error:
+ if(ran_lw) htrdr_ran_lw_ref_put(ran_lw);
+ goto exit;
+}
+
+void
+htrdr_ran_lw_ref_get(struct htrdr_ran_lw* ran_lw)
+{
+ ASSERT(ran_lw);
+ ref_get(&ran_lw->ref);
+}
+
+void
+htrdr_ran_lw_ref_put(struct htrdr_ran_lw* ran_lw)
+{
+ ASSERT(ran_lw);
+ ref_put(&ran_lw->ref, release_ran_lw);
+}
+
+double
+htrdr_ran_lw_sample
+ (const struct htrdr_ran_lw* ran_lw,
+ const double r)
+{
+ ASSERT(ran_lw);
+ if(ran_lw->nbands != HTRDR_RAN_LW_CONTINUE) {
+ return ran_lw_sample_discreet(ran_lw, r, FUNC_NAME);
+ } else if(eq_eps(ran_lw->range[0], ran_lw->range[1], 1.e-6)) {
+ return ran_lw->range[0];
+ } else {
+ return ran_lw_sample_continue(ran_lw, r, FUNC_NAME);
+ }
+}
+
+double
+htrdr_ran_lw_get_sky_band_pdf
+ (const struct htrdr_ran_lw* ran_lw,
+ const size_t iband)
+{
+ size_t i;
+ size_t nbands;
+ (void)nbands;
+ ASSERT(ran_lw);
+
+ nbands = htsky_get_spectral_bands_count(ran_lw->htrdr->sky);
+ ASSERT(nbands);
+ ASSERT(iband >= htsky_get_spectral_band_id(ran_lw->htrdr->sky, 0));
+ ASSERT(iband <= htsky_get_spectral_band_id(ran_lw->htrdr->sky, nbands-1));
+
+ /* Compute the index of the spectral band */
+ i = iband - htsky_get_spectral_band_id(ran_lw->htrdr->sky, 0);
+
+ /* Fetch its PDF */
+ ASSERT(i < darray_double_size_get(&ran_lw->sky_bands_pdf));
+ return darray_double_cdata_get(&ran_lw->sky_bands_pdf)[i];
+}
+
diff --git a/src/htrdr_ran_lw.h b/src/htrdr_ran_lw.h
@@ -0,0 +1,56 @@
+/* Copyright (C) 2018, 2019, 2020 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_RAN_LW_H
+#define HTRDR_RAN_LW_H
+
+#include <rsys/rsys.h>
+
+#define HTRDR_RAN_LW_CONTINUE 0
+
+struct htrdr;
+struct htrdr_ran_lw;
+
+extern LOCAL_SYM res_T
+htrdr_ran_lw_create
+ (struct htrdr* htrdr,
+ const double range[2], /* Must be included in [1000, 100000] nanometers */
+ /* # bands used to discretisze the LW domain. HTRDR_RAN_LW_CONTINUE <=> no
+ * discretisation */
+ const size_t nbands, /* Hint on #bands used to discretised th CDF */
+ const double ref_temperature, /* Reference temperature */
+ struct htrdr_ran_lw** ran_lw);
+
+extern LOCAL_SYM void
+htrdr_ran_lw_ref_get
+ (struct htrdr_ran_lw* ran_lw);
+
+extern LOCAL_SYM void
+htrdr_ran_lw_ref_put
+ (struct htrdr_ran_lw* ran_lw);
+
+/* Return a wavelength in nanometer */
+extern LOCAL_SYM double
+htrdr_ran_lw_sample
+ (const struct htrdr_ran_lw* ran_lw,
+ const double r); /* Canonical number in [0, 1[ */
+
+extern LOCAL_SYM double
+htrdr_ran_lw_get_sky_band_pdf
+ (const struct htrdr_ran_lw* ran_lw,
+ const size_t iband);
+
+#endif /* HTRDR_RAN_LW_H */
+
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -26,8 +26,29 @@ struct htrdr_accum {
size_t nweights; /* #accumlated weights */
size_t nfailures; /* #failures */
};
-#define HTRDR_ACCUM_NULL__ {0,0,0,0}
-static const struct htrdr_accum HTRDR_ACCUM_NULL = HTRDR_ACCUM_NULL__;
+static const struct htrdr_accum HTRDR_ACCUM_NULL;
+
+/* Monte carlo estimate */
+struct htrdr_estimate {
+ double E; /* Expected value */
+ double SE; /* Standard error */
+};
+static const struct htrdr_estimate HTRDR_ESTIMATE_NULL;
+
+struct htrdr_pixel_sw {
+ struct htrdr_estimate X; /* In W/m^2/sr */
+ struct htrdr_estimate Y; /* In W/m^2/sr */
+ struct htrdr_estimate Z; /* In W/m^2/sr */
+ struct htrdr_accum time; /* In microseconds */
+};
+static const struct htrdr_pixel_sw HTRDR_PIXEL_SW_NULL;
+
+struct htrdr_pixel_lw {
+ struct htrdr_estimate radiance; /* In W/m^2/sr */
+ struct htrdr_estimate radiance_temperature; /* In K */
+ struct htrdr_accum time; /* In microseconds */
+};
+static const struct htrdr_pixel_lw HTRDR_PIXEL_LW_NULL;
/* Forward declarations */
struct htrdr;
@@ -42,18 +63,30 @@ htrdr_compute_radiance_sw
struct ssp_rng* rng,
const double pos[3],
const double dir[3],
+ const double wlen, /* In nanometer */
+ const size_t iband, /* Index of the spectral band */
+ const size_t iquad); /* Index of the quadrature point into the band */
+
+extern LOCAL_SYM double
+htrdr_compute_radiance_lw
+ (struct htrdr* htrdr,
+ const size_t ithread,
+ struct ssp_rng* rng,
+ const double pos[3],
+ const double dir[3],
+ const double wlen, /* In nanometer */
const size_t iband, /* Index of the spectral band */
const size_t iquad); /* Index of the quadrature point into the band */
extern LOCAL_SYM res_T
-htrdr_draw_radiance_sw
+htrdr_draw_radiance
(struct htrdr* htrdr,
const struct htrdr_camera* cam,
const size_t width, /* Image width */
const size_t height, /* Image height */
const size_t spp, /* #samples per pixel, i.e. #realisations */
/* Buffer of struct htrdr_accum[4]. May be NULL on non master processes */
- struct htrdr_buffer* buf);
+ struct htrdr_buffer* buf);
extern LOCAL_SYM int
htrdr_ground_filter
@@ -66,14 +99,13 @@ htrdr_ground_filter
static FINLINE void
htrdr_accum_get_estimation
(const struct htrdr_accum* acc,
- double* expected_value,
- double* std_err)
+ struct htrdr_estimate* estimate)
{
- ASSERT(acc && expected_value && std_err);
+ ASSERT(acc && estimate);
if(!acc->nweights) {
- *expected_value = 0;
- *std_err = 0;
+ estimate->E = 0;
+ estimate->SE = 0;
} else {
const double N = (double)acc->nweights;
double E, V, SE;
@@ -81,8 +113,8 @@ htrdr_accum_get_estimation
V = MMAX(acc->sum_weights_sqr / N - E*E, 0);
SE = sqrt(V/N);
- *expected_value = E;
- *std_err = SE;
+ estimate->E = E;
+ estimate->SE = SE;
}
}
diff --git a/src/htrdr_sun.c b/src/htrdr_sun.c
@@ -45,14 +45,6 @@ struct htrdr_sun {
/*******************************************************************************
* Helper functions
******************************************************************************/
-static INLINE int
-cmp_dbl(const void* a, const void* b)
-{
- const double d0 = *((const double*)a);
- const double d1 = *((const double*)b);
- return d0 < d1 ? -1 : (d0 > d1 ? 1 : 0);
-}
-
static void
release_sun(ref_T* ref)
{
