commit dd428386847e549d91aa7aa38b7c865e5805903e
parent 4785a3d5f5faa607bf552d5c3f8990568b65509b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 10 Jan 2023 12:28:08 +0100
Merge branch 'feature_planeto' into develop
Diffstat:
45 files changed, 5899 insertions(+), 957 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -29,6 +29,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
option(HTRDR_BUILD_ATMOSPHERE "Build the htrdr-atmosphere program" ON)
option(HTRDR_BUILD_COMBUSTION "Build the htrdr-combustion program" ON)
+option(HTRDR_BUILD_PLANETO "Build the htrdr-planeto program" ON)
set(HTRDR_BUILD_COMMANDS "")
if(HTRDR_BUILD_ATMOSPHERE)
@@ -37,6 +38,9 @@ endif()
if(HTRDR_BUILD_COMBUSTION)
set(HTRDR_BUILD_COMMANDS "${HTRDR_BUILD_COMMANDS};HTRDR_BUILD_COMBUSTION")
endif()
+if(HTRDR_BUILD_COMBUSTION)
+ set(HTRDR_BUILD_COMMANDS "${HTRDR_BUILD_COMMANDS};HTRDR_BUILD_PLANETO")
+endif()
################################################################################
# Add sub projects
@@ -49,6 +53,9 @@ endif()
if(HTRDR_BUILD_COMBUSTION)
add_subdirectory(combustion)
endif()
+if(HTRDR_BUILD_PLANETO)
+ add_subdirectory(planeto)
+endif()
add_subdirectory(commands)
add_subdirectory(doc)
diff --git a/cmake/atmosphere/CMakeLists.txt b/cmake/atmosphere/CMakeLists.txt
@@ -82,7 +82,7 @@ rcmake_prepend_path(HTRDR_ATMOSPHERE_FILES_INC ${HTRDR_SOURCE_DIR}/atmosphere)
rcmake_prepend_path(HTRDR_ATMOSPHERE_FILES_INC2 ${HTRDR_BUILD_DIR}/atmosphere)
# Atmosphere library
-add_library(htrdr-atmosphere SHARED
+add_library(htrdr-atmosphere STATIC
${HTRDR_ATMOSPHERE_FILES_SRC}
${HTRDR_ATMOSPHERE_FILES_INC}
${HTRDR_ATMOSPHERE_FILES_INC2})
@@ -94,7 +94,7 @@ if(CMAKE_COMPILER_IS_GNUCC)
endif()
set_target_properties(htrdr-atmosphere PROPERTIES
- DEFINE_SYMBOL HTRDR_SHARED_BUILD
+ DEFINE_SYMBOL HTRDR_STATIC
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
diff --git a/cmake/combustion/CMakeLists.txt b/cmake/combustion/CMakeLists.txt
@@ -81,15 +81,15 @@ set(HTRDR_COMBUSTION_FILES_INC
rcmake_prepend_path(HTRDR_COMBUSTION_FILES_SRC ${HTRDR_SOURCE_DIR}/combustion)
rcmake_prepend_path(HTRDR_COMBUSTION_FILES_INC ${HTRDR_SOURCE_DIR}/combustion)
-# Atmosphere library
-add_library(htrdr-combustion SHARED
+# Combustion library
+add_library(htrdr-combustion STATIC
${HTRDR_COMBUSTION_FILES_SRC}
${HTRDR_COMBUSTION_FILES_INC})
target_link_libraries(htrdr-combustion
htrdr-core AtrSTM RSys Star3D StarCamera StarSF StarSP)
set_target_properties(htrdr-combustion PROPERTIES
- DEFINE_SYMBOL HTRDR_SHARED_BUILD
+ DEFINE_SYMBOL HTRDR_STATIC
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
diff --git a/cmake/commands/CMakeLists.txt b/cmake/commands/CMakeLists.txt
@@ -28,6 +28,9 @@ endif()
if(HTRDR_BUILD_COMBUSTION)
list(APPEND _link_libraries htrdr-combustion)
endif()
+if(HTRDR_BUILD_PLANETO)
+ list(APPEND _lin_libraries htrdr-planeto)
+endif()
################################################################################
# Check dependencies
@@ -62,10 +65,19 @@ if(HTRDR_BUILD_COMBUSTION)
target_link_libraries(htrdr_combustion_cmd htrdr-combustion)
endif()
+add_executable(htrdr_planeto_cmd
+ ${HTRDR_SOURCE_DIR}/commands/htrdr_planeto_cmd.c)
+set_target_properties(htrdr_planeto_cmd PROPERTIES
+ COMPILE_DEFINITIONS "${HTRDR_BUILD_COMMANDS}"
+ OUTPUT_NAME htrdr-planeto)
+if(HTRDR_BUILD_PLANETO)
+ target_link_libraries(htrdr_planeto_cmd htrdr-planeto)
+endif()
+
################################################################################
# Define output & install directories
################################################################################
-install(TARGETS htrdr_cmd htrdr_atmosphere_cmd htrdr_combustion_cmd
+install(TARGETS htrdr_cmd htrdr_atmosphere_cmd htrdr_combustion_cmd htrdr_planeto_cmd
ARCHIVE DESTINATION bin
LIBRARY DESTINATION lib
RUNTIME DESTINATION bin)
diff --git a/cmake/core/CMakeLists.txt b/cmake/core/CMakeLists.txt
@@ -84,12 +84,13 @@ set(HTRDR_CORE_FILES_SRC
htrdr.c
htrdr_args.c
htrdr_buffer.c
- htrdr_cie_xyz.c
htrdr_draw_map.c
htrdr_geometry.c
htrdr_log.c
htrdr_materials.c
- htrdr_ran_wlen.c
+ htrdr_ran_wlen_cie_xyz.c
+ htrdr_ran_wlen_discrete.c
+ htrdr_ran_wlen_planck.c
htrdr_rectangle.c
htrdr_slab.c
htrdr_spectral.c)
@@ -98,13 +99,14 @@ set(HTRDR_CORE_FILES_INC
htrdr_accum.h
htrdr_buffer.h
htrdr_c.h
- htrdr_cie_xyz.h
htrdr_draw_map.h
htrdr_geometry.c
htrdr_interface.h
htrdr_log.h
htrdr_materials.h
- htrdr_ran_wlen.h
+ htrdr_ran_wlen_cie_xyz.h
+ htrdr_ran_wlen_discrete.h
+ htrdr_ran_wlen_planck.h
htrdr_rectangle.h
htrdr_sensor.h
htrdr_slab.h
diff --git a/cmake/doc/CMakeLists.txt b/cmake/doc/CMakeLists.txt
@@ -49,6 +49,12 @@ if(HTRDR_BUILD_COMBUSTION)
list(APPEND MAN_SOURCES ${CMAKE_CURRENT_BINARY_DIR}/htrdr-combustion.1.scd)
endif()
+if(HTRDR_BUILD_PLANETO)
+ configure_file(${HTRDR_DOC_DIR}/htrdr-planeto.1.scd.in
+ ${CMAKE_CURRENT_BINARY_DIR}/htrdr-planeto.1.scd @ONLY)
+ list(APPEND MAN_SOURCES ${CMAKE_CURRENT_BINARY_DIR}/htrdr-planeto.1.scd)
+endif()
+
set(MAN_FILES)
set(MAN5_FILES)
set(MAN1_FILES)
diff --git a/cmake/planeto/CMakeLists.txt b/cmake/planeto/CMakeLists.txt
@@ -0,0 +1,118 @@
+# Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+# Copyright (C) 2018, 2019, 2021 CNRS
+# Copyright (C) 2018, 2019 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/>.
+
+cmake_minimum_required(VERSION 3.1)
+project(htrdr-planeto)
+
+################################################################################
+# Check dependencies
+################################################################################
+find_package(RCMake 0.4 REQUIRED)
+find_package(RNATM 0.0 REQUIRED)
+find_package(RNGRD 0.0 REQUIRED)
+find_package(RSys 0.11 REQUIRED)
+find_package(Star3D 0.8 REQUIRED)
+find_package(StarBuffer 0.0 REQUIRED)
+find_package(StarCamera 0.0 REQUIRED)
+find_package(StarSF 0.6 REQUIRED)
+find_package(StarSP 0.12 REQUIRED)
+find_package(StarVX 0.1 REQUIRED)
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
+include(rcmake)
+include(rcmake_runtime)
+
+include_directories(
+ ${RNATM_INCLUDE_DIR}
+ ${RNGRD}
+ ${RSys_INCLUDE_DIR}
+ ${Star3D_INCLUDE_DIR}
+ ${StarBuffer_INCLUDE_DIR}
+ ${StarCamera_INCLUDE_DIR}
+ ${StarSF_INCLUDE_DIR}
+ ${StarSP_INCLUDE_DIR}
+ ${StarVX_INCLUDE_DIR}
+ ${HTRDR_BUILD_DIR}
+ ${HTRDR_SOURCE_DIR})
+
+################################################################################
+# Generate files
+################################################################################
+set(HTRDR_PLANETO_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD "1" CACHE INTERNAL "")
+set(HTRDR_PLANETO_ARGS_DEFAULT_GRID_DEFINITION_HINT "512" CACHE INTERNAL "")
+
+configure_file(${HTRDR_SOURCE_DIR}/planeto/htrdr_planeto_args.h.in
+ ${HTRDR_BUILD_DIR}/planeto/htrdr_planeto_args.h @ONLY)
+
+################################################################################
+# Configure and define targets
+################################################################################
+set(HTRDR_PLANETO_FILES_SRC
+ htrdr_planeto.c
+ htrdr_planeto_args.c
+ htrdr_planeto_compute_radiance.c
+ htrdr_planeto_draw_map.c
+ htrdr_planeto_main.c
+ htrdr_planeto_source.c)
+
+set(HTRDR_PLANETO_FILES_INC
+ htrdr_planeto.h
+ htrdr_planeto_c.h
+ htrdr_planeto_source.h)
+
+# Prepend each file in the `HTRDR_FILES_<SRC|INC>' list by `HTRDR_SOURCE_DIR'
+rcmake_prepend_path(HTRDR_PLANETO_FILES_SRC ${HTRDR_SOURCE_DIR}/planeto)
+rcmake_prepend_path(HTRDR_PLANETO_FILES_INC ${HTRDR_SOURCE_DIR}/planeto)
+
+# Planeto library
+add_library(htrdr-planeto STATIC
+ ${HTRDR_PLANETO_FILES_SRC}
+ ${HTRDR_PLANETO_FILES_INC})
+target_link_libraries(htrdr-planeto
+ htrdr-core RNATM RNGRD RSys Star3D StarBuffer StarCamera StarSF StarSP)
+
+set_target_properties(htrdr-planeto PROPERTIES
+ DEFINE_SYMBOL HTRDR_STATIC
+ VERSION ${VERSION}
+ SOVERSION ${VERSION_MAJOR})
+
+################################################################################
+# Add tests
+################################################################################
+if(NOT NO_TEST)
+ function(build_test _name)
+ add_executable(${_name}
+ ${HTRDR_SOURCE_DIR}/planeto/${_name}.c)
+ target_link_libraries(${_name} htrdr-core htrdr-planeto ${ARGN})
+ endfunction()
+
+ function(new_test _name)
+ build_test(${_name} ${ARGN})
+ add_test(${_name} ${_name})
+ endfunction()
+
+ new_test(test_htrdr_planeto_source StarBuffer)
+endif()
+
+################################################################################
+# Define output & install directories
+################################################################################
+install(TARGETS htrdr-planeto
+ ARCHIVE DESTINATION bin
+ LIBRARY DESTINATION lib
+ RUNTIME DESTINATION bin)
+
diff --git a/doc/htrdr-planeto.1.scd.in b/doc/htrdr-planeto.1.scd.in
@@ -0,0 +1,407 @@
+htrdr-planeto(1)
+
+; Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+; Copyright (C) 2018, 2019, 2021 CNRS
+; Copyright (C) 2018, 2019 Université Paul Sabatier
+; (contact-edstar@laplace.univ-tlse.fr)
+;
+; 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/>.
+
+# NAME
+
+htrdr-planeto - simulate radiative transfer in 3D planetary atmosphere
+
+# SYNOPSIS
+
+htrdr-planeto [_option_] ... -G _ground_ -g _gas_
+
+# DESCRIPTION
+
+*htrdr-planeto* simulates the radiative transfer of a terrestrial planet in the
+visible or the infrared part of the spectrum. The planet's soil (option *-G*)
+can be any set of triangles with BRDFs and temperatures defined per triangle.
+The atmosphere is composed of a gas mixture (option *-g*) and a potentially
+empty set of aerosols (option *-a*). Both can have arbitrary tetrahedral meshes
+with per-node radiative properties. Rayleigh is used as a gas phase function and
+the temperature of the gas is defined by the node of the mesh. Aerosol phase
+functions (Henyey and Greenstein or measured) are also defined per node.
+
+*htrdr-planeto* is mainly a renderer that calculates an image (option *-i*)
+for a given observation position (option *-C*). Its internal rendering algorithm
+is based on Monte-Carlo integration, which consists for each pixel of simulating
+a given number of optical paths from the sensor, taking into account the
+phenomena of light absorption and scattering.
+
+*htrdr-planeto* offers three ways to perform spectral integration (*-s* option).
+By default, it calculates an image for the visible part of the spectrum between
+380 and 780 nanometers, for the three components of the CIE 1931 XYZ color space
+which are then recombined to obtain the final color for each pixel, and finally
+the entire image of the scene as seen from the observation position. The other
+two methods are to explicitly define the longwave or shortwave spectral range to
+be integrated and continuously sample a wavelength in this range. In fact,
+longwave and shortwave are keywords that mean that the source of radiation is
+either external or internal to the medium, respectively. In shortwave, only
+radiance is evaluated and stored in the output image. For longwave rendering,
+this estimated radiance is then converted to brightness temperature and both are
+recorded in the image.
+
+In *htrdr-planeto*, the spatial unit 1.0 corresponds to one meter and
+temperatures are expressed in Kelvin. The estimated radiances are given in
+W/sr/m², except for monochromatic calculations where the calculated spectral
+radiance is defined in W/sr/m²/nm.
+
+*htrdr-planeto* implements mixed parallelism. On a single computer (that is, a
+node), it uses shared memory parallelism while it relies on Message Passing
+Interface (MPI) to parallelize calculations between multiple nodes.
+*htrdr-planeto* can therefore be launched either directly or via a process
+launcher such as *mpirun*(1) to distribute the rendering on several computers.
+
+# OPTIONS
+
+*-a* <_aerosol-parameter_:...>
+ Define an aerosol. Use this option as many times as there are aerosols to be
+ defined. Available aerosol parameters are:
+
+ *mesh*=_path_
+ Path to the *smsh*(5) file that stores the aerosol tetrahedral mesh.
+
+ *name*=_string_
+ Name assigned to the aerosol.
+
+ *radprop*=_path_
+ Path to the *sars*(5) file that stores the radiative properties of the
+ aerosol. Radiative properties are defined per volumetric mesh node and,
+ therefore, this file and the volumetric mesh file (see *mesh* parameter)
+ must be consistent with each other.
+
+ *phasefn*=_path_
+ Path to the *rnsl*(5) file that lists the *rnsf*(5) files to load; each of
+ theses files stores an aerosol phase function. The phase function to be used
+ per volumetric mesh node is defined in another file (see *phaseids*
+ parameter).
+
+ *phaseids*=_path_
+ Path to the *rnpfi*(5) file that stores the index of the phase function to
+ be used per volumetric mesh node; the list of phase function is defined in
+ another file (see *phasefn* parameter). Note that this file must be
+ consistent with the volumetric mesh defined in the *mesh* parameter.
+
+*-C* <_camera-parameter_:...>
+ Configure a perspective camera. Available parameters are:
+
+ *focal-dst*=_dst_
+ Distance to focus on with a thin lens camera, that is, a camera whose
+ *lens-radius* is not zero. The default focal distance is
+ @HTRDR_ARGS_DEFAULT_CAMERA_PERSPECTIVE_FOCAL_DST@ meter.
+
+ *focal-length*=_length_
+ Focal length of a camera lens. It is another way to control the field of
+ view of a thin lens camera. By default, the field of view is directly set
+ through the **fov** parameter.
+
+ *fov*=_angle_
+ Vertical field of view of the camera in
+ \]@HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MIN@,
+ @HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MAX@[ degrees. By
+ default _angle_ is set to
+ @HTRDR_ARGS_DEFAULT_CAMERA_PERSPECTIVE_FOV@ degrees.
+
+ *lens-radius*=_radius_
+ Radius of the camera lens. A non-zero radius means that the camera is a thin
+ lens camera while a zero radius defines a pinhole camera. By default the
+ lens radius is @HTRDR_ARGS_DEFAULT_CAMERA_PERSPECTIVE_LENS_RADIUS@.
+
+ *pos*=_x_,_y_,_z_
+ Camera lens position. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_CAMERA_POS@}.
+
+ *tgt*=_x_,_y_,_z_
+ Position targeted by the camera. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_CAMERA_TGT@}.
+
+ *up*=_x_,_y_,_z_
+ Up vector of the camera. By default it is set to
+ {@HTRDR_ARGS_DEFAULT_CAMERA_UP@}.
+
+*-d*
+ Write atmospheric acceleration structures to _output_. Each structure is saved
+ in VTK ASCII file format [1]. To divide the resulting output into _n_ files
+ (_n_ > 1), each storing an acceleration structure, one can use the *csplit*(1)
+ command as below:
+
+ ```
+ csplit -f octree -k output %^#\\ vtk% /^#\\ vtk/ {$(($(grep -ce "^#\\ vtk" output)-2))}
+ ```
+
+*-f*
+ Force overwrite the _output_ file.
+
+*-G* <_ground-parameter_:...>
+ Define the ground of the planet. Available ground parameters are:
+
+ *brdf*=_path_
+ Path to the *rnsl*(5) file that lists the *mrumtl*(5) files to load; each of
+ these files stores a ground BRDF. The BRDF to be used per ground node is
+ defined in another file (see *prop* parameter).
+
+ *mesh*=_path_
+ Path to the *smsh*(5) file which stores the triangular mesh of the ground.
+
+ *name*=_string_
+ Name assigned to the ground.
+
+ *prop*=_path_
+ Path to the *rnsp*(5) file that stores ground surface properties. The
+ properties (BRDF index and temperature) are defined per triangle and,
+ therefore, this file and the mesh file (see *mesh* parameter) must be
+ consistent with each other.
+
+*-g* <_gas-parameter_:...>
+ Define the gas mixture. Available gas parameters are:
+
+ *mesh*=_path_
+ Path to the *smsh*(5) file that stores the gas tetrahedral mesh.
+
+ *ck*=_path_
+ Path to the *sck*(5) file that stores the correlated K of the gas. The CKs
+ are defined per volumetric mesh node and, therefore, this file and the
+ volumetric mesh file (see *mesh* parameter) must be consistent with each
+ other.
+
+ *temp*=_path_
+ Path to the *rngt*(5) file that stores the temperature of the gas. The
+ temperature is defined per volumetric mesh node and, therefore, this file
+ and the volumetric mesh file (see *mesh* parameter) must be consistent with
+ each other.
+
+*-h*
+ Display short help and exit.
+
+*-i* <_image-parameter_:...>
+ Define the image to compute. Available image parameters are:
+
+ *def*=<_width_>x<_height_>
+ Image definition. By default the image definition is
+ @HTRDR_ARGS_DEFAULT_IMG_WIDTH@x@HTRDR_ARGS_DEFAULT_IMG_HEIGHT@.
+
+ *spp*=_samples-count_
+ Number of samples to estimate one pixel, i.e. number of radiative paths
+ sampled per pixel. By default, *spp* is set to
+ @HTRDR_ARGS_DEFAULT_IMG_SPP@.
+
+*-N*
+ Precalculate tetrahedron normals. This speeds up runtime performance by
+ calculating normals once and for all rather than re-evaluating them every time
+ a tetrahedron is queried at a given position. In return, the memory space used
+ to store normals increases the memory footprint.
+
+*-O* _storage_
+ File where atmospheric acceleration structures are stored/loaded. If _storage_
+ does not exist, it is created and is used to store the built acceleration
+ structures.
+
+ If _storage_ exists, acceleration structures are loaded from it rather than
+ built from scratch, resulting in significant acceleration of the preprocessing
+ step. Note that if the data structures stored in _storage_ are not as expected
+ (that is, the input atmospheric data or construction parameters are
+ different), an error is notified and execution is stopped, thus avoiding the
+ use of incorrect acceleration structures.
+
+*-o* _output_
+ File to write the output data. The output data is either an image or atmospheric
+ acceleration structures if the *-d* option is set. If it is not defined, the
+ data is written to the standard output.
+
+*-S* <_source-parameter_:...>
+ Define the external source. Available source parameters are:
+
+ *lat*=_real_
+ The latitude of the source, i.e. its angle between [-90, 90] degrees about
+ the x-axis. The default latitude is set to 0.
+
+ *lon*=_real_
+ The longitude of the source, i.e. its angle between [-180, 180] degrees
+ about the z-axis. The default longitude is set to 0.
+
+ *dst*=_real_
+ Distance in meters from source to origin. The default distance is 0.
+
+ *radius*=_real_
+ Source radius in meters.
+
+ *temp*=_real_
+ Source temperature in Kelvin.
+
+*-s* <_spectral-parameter_:...>
+ Configure spectral integration. Available spectral parameters are:
+
+ *cie_xyz*
+ The radiance is calculated for the visible part of the spectrum between
+ 380 nm and 780 nm by sampling the wavelength relative to the XYZ CIE 1931
+ color space. This is the default spectral integration.
+
+ *lw*=_wlen-min_,_wlen-max_
+ Perform continuous spectral sampling in the wavelength range [_wlen-min_, _wlen-max_]
+ (wavelengths must be provided in nanometers) according to the Planck
+ function for a reference temperature which is the maximum ground
+ temperature, which is assumed to be the maximum scene temperature. If
+ _wlen-min_ and _wlen-max_ are equal, the calculation is monochromatic. *lw*
+ means longwaves but is here a code word that actually means "calculation of
+ radiance using the internal source of radiation": in other words, radiation
+ is emitted by the medium and its limits (ground and space).
+
+ *sw*=_wlen-min_,_wlen-max_
+ Perform continuous spectral sampling in the wavelength range [_wlen-min_, _wlen-max_]
+ (wavelengths must be provided in nanometers) according to the Planck
+ function for a reference temperatura which is the source temperature. If
+ _wlen-min_ and _wlen-max_ are equal, the calculation is monochromatic. Here,
+ *sw* means shortwaves, i.e. the radiation source is external to the medium
+ (option *-S*).
+
+*-T* _optical-thickness_
+ Optical thickness used as a criterion to construct atmospheric acceleration
+ structures. Its default value is
+ @HTRDR_PLANETO_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@.
+
+*-t* _threads-count_
+ Hint on the number of threads to use. Default assumes as many threads as CPU
+ cores.
+
+*-V* _definition_
+ Advice on the definition of the atmospheric acceleration structures. Its
+ default value is
+ @HTRDR_PLANETO_ARGS_DEFAULT_GRID_DEFINITION_HINT@.
+
+*-v*
+ Make the command verbose.
+
+# OUTPUT IMAGE
+
+Images calculated by *htrdr-planeto* are saved in *htrdr-image*(5) file format.
+This section describes the nature and arrangement of the output data depending
+on the type of calculation.
+
+## XYZ image
+
+For image rendering in the visible part of the spectrum (default behavior or
+when the *-s cie_xyz* option is set), each pixel stores 4 estimates, or 8
+floating-point values. The first, second and third pairs of numbers store the
+estimated radiation in W/sr/m² for the X, Y, and Z components of the CIE 1931
+XYZ color space. For each pair, the first corresponds to the expected value and
+the second its standard error. Finally, the fourth and last pair records the
+estimate of the calculation time in µs of a radiative path (expected value and
+standard error).
+
+## Longwave image
+
+If the image is an infrared rendering (option *-s lw*=_wlen-min_,_wlen-max_),
+the first and second pixel values store the expected value and standard error of
+the estimated brightness temperature in Kelvin. The third and fourth values
+record the expected value and standard error of the estimated radiance, which is
+either integrated radiance in W/sr/m² or spectral radiance in W/sr/m²/nm
+depending on whether this radiance was calculated for a spectral range or at a
+single wavelength. The fifth and sixth values are not used and are therefore set
+to 0. Finally, the last 2 components of the pixel record the expected value and
+the standard error in µs of the calculation time per radiative path.
+
+## Shortwave image
+
+For shortwave renderings (option *-s sw*=_wlen-min_,_wlen-max_), the image
+layout is the same as for infrared rendering, except for the first and second
+pixel values that are not used. That is, the third and fourth values record the
+estimated radiance per pixel and the seventh and eighth values store the
+estimate of the calculation time by radiative path. The other values are set to
+0.
+
+# EXAMPLES
+
+The following command line runs *htrdr-planeto* in a verbose way (option *-v*)
+to calculate an _800_ by _600_ pixel image by sampling _64_ radiative paths per
+pixel for the 3 components of the CIE XYZ 1931 color space. The external source
+is positioned at _-45_ degrees longitude and _50_ degrees latitude relative to
+the absolute referential. The camera looks at the origin (*tgt=*_0_,_0_,_0_)and
+is positioned at _1.5e7_ meters along the Y axis with an image plane aligned
+along the Z axis (*up=*_0_,_0_,_1_). Its vertical field of view is _70_ degrees.
+The gas of the planetary atmosphere is described by the tetrahedral mesh
+recorded in the _gas.smsh_ file, while its spectral data and temperature are
+given by the files _gas.sck_ and _gas.rngt_, respectively. Two aerosols complete
+the planetary atmosphere: one for _clouds_ and one for _haze_. Their respective
+meshes are stored in the _a<1|2>.smsh_ files while their radiative properties
+are given by the _a<1|2>.sars_ files. Finally, their phase functions are
+described by a set of 2 files: the _a<1|2>.rnsf_ file which lists the aerosol
+phase functions and the _a<1|2>.rnpfi_ file which references them by volumetric
+mesh node. To speed up rendering time, the normals of the tetrahedral meshes of
+the gas and aerosols are precalculated once and for all (option *-N*). The
+ground geometry is stored in the _ground.smsh_ file with its triangle properties
+(temperature and BRDF) defined in the _ground.rnsp_ file. The referenced BRDFs
+are listed in the _ground.rnsl_ file. The definition of acceleration structures
+cannot exceed _512³_ and its voxels can be merged until their optical thickness
+is greater than _10_. These structures are either reloaded from
+_storage_cie.bin_ or built from scratch and stored in _storage_cie.bin_
+depending on whether that file exists or not. Finally, the calculated images are
+stored in the _image_CIE_XYZ.ht_ file even if the file already exists (option
+*-f*).
+
+```
+htrdr-planeto -v -N \
+ -i def=800x600:spp=64 \
+ -s cie_xyz \
+ -S lon=-45:lat=50:dst=1.5e8:radius=6.9e5:temp=5778 \
+ -C pos=0,1.5e7,0:tgt=0,0,0:up=0,0,1:fov=70 \
+ -g mesh=gas.smsh:ck=gas.sck:temp=gas.rngt \
+ -a mesh=a1.smsh:radprop=a1.sars:phasefn=a1.rnsf:phaseids=a1.rnpfi:name=clouds \
+ -a mesh=a2.smsh:radprop=a2.sars:phasefn=a2.rnsf:phaseids=a2.rnpfi:name=haze \
+ -G mesh=ground.smsh:prop=ground.rnsp:brdf=ground.rnsl:name=namek \
+ -V 512 -T 10 -O storage_cie.bin \
+ -f -o image_CIE_XYZ.ht
+```
+
+The next command line is the same as the previous one, except that it calculates
+an infrared image between _10,000_ nm and _20,000_ nm (option *-s*). Note that
+the acceleration structure storage file is no longer the same (_storage_lw.bin_
+rather than _storage_cie.bin_). Indeed, the previous one records the
+acceleration structures for the spectral range of the CIE XYZ color space, while
+one wants to store/reload the acceleration structures for a spectral range
+between 10 and 20 µm (see *-O* option). In any case, if the previous storage,
+incompatible with the current spectral range, had been submitted, the command
+would have stopped with an error message, thus avoiding the use of the wrong
+accelerartion structures.
+
+```
+htrdr-planeto -v -N \
+ -i def=800x600:spp=64 \
+ -s lw=10000,20000 \
+ -C pos=0,1.5e7,0:tgt=0,0,0:up=0,0,1:fov=70 \
+ -g mesh=gas.smsh:ck=gas.sck:temp=gas.rngt \
+ -a mesh=a1.smsh:radprop=a1.sars:phasefn=a1.rnsf:phaseids=a1.rnpfi:name=clouds \
+ -a mesh=a2.smsh:radprop=a2.sars:phasefn=a2.rnsf:phaseids=a2.rnpfi:name=haze \
+ -G mesh=ground.smsh:prop=ground.rnsp:brdf=ground.rnsl:name=namek \
+ -V 512 -T 10 -O storage_lw.bin \
+ -f -o image_infrared.ht
+```
+
+# SEE ALSO
+
+. VTK file format -
+ <http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf>
+
+*htpp*(1),
+*htrdr-image*(5),
+*mpirun*(1),
+*mrumtl*(5),
+*rnpfi*(5),
+*rnsf*(5),
+*rnsl*(5),
+*sars*(5),
+*smsh*(5)
diff --git a/doc/htrdr.1.scd b/doc/htrdr.1.scd
@@ -55,6 +55,9 @@ The available _htrdr-<mode>_ commands are:
*htrdr-combustion*(1)
Radiative transfer computations in a combustion medium.
+*htrdr-planeto*(1)
+ Radiative transfer computations in 3D planetary atmosphere.
+
# COPYRIGHT
Copyright © 2018, 2019, 2020, 2021 |Méso|Star> <contact@meso-star.com>++
@@ -70,4 +73,5 @@ redistribute it. There is NO WARRANTY, to the extent permitted by law.
# SEE ALSO
*htrdr-atmosphere*(1),
-*htrdr-combustion*(1)
+*htrdr-combustion*(1),
+*htrdr-planeto*(1)
diff --git a/doc/rnrl.scd b/doc/rnrl.scd
@@ -0,0 +1,85 @@
+rnrl(5)
+
+; Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+; Copyright (C) 2018, 2019, 2021 CNRS
+; Copyright (C) 2018, 2019 Université Paul Sabatier
+; (contact-edstar@laplace.univ-tlse.fr)
+; 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/>.
+
+# NAME
+
+rnsr - Rad-Net Radiance List file format
+
+# DESCRIPTION
+
+*rnrl* is a binary file format for storing a list of spectrally varying
+radiances. The radiances (in W/m²/sr/m) are listed by wavelength (in nm)
+sorted in ascending order.
+
+A *rnrl* file is actually a Star-Buffer file (see *sbuf*(5)). It starts with a
+header of 4 64-bit integers describing the layout of the data. The first integer
+is a power of two (usually 4096) that defines the size of the memory page in
+bytes to which the list of per wavelength radiances aligns (_pagesize_). The
+second integer is the _size_ of the array, i.e. the number of wavelength for
+which a radiance is defined. Finally, the 2 remaining integers store the memory
+size (16 bytes, i.e. 2 double precision values) and the memory alignment (16
+bytes) of a per wavelength radiance.
+
+An *rnrl* file is actually a Star-Buffer file (see *sbuf*(5)). It starts with a
+header of 4 64-bit integers describing the layout of the data. The first integer
+is a power of two (usually 4096) which defines the size of the memory page in
+bytes on which the list of radiances by wavelength aligns (_pagesize_). The
+second integer is the _size_ of the array, that is, the number of wavelengths
+for which a radiance is defined. Finally, the remaining 2 integers store the
+memory size (16 bytes) and memory alignment (16 bytes) of a radiance per
+wavelength.
+
+The fill bytes follow the file header to align the radiances by wavelength to
+_pagesize_.
+
+Each item in the list is composed of 2 double-precision floating-point values:
+the wavelength in nanometers and its corresponding radiance in W/m²/sr/m.
+
+The end of the file is eventually padded with dummy bytes to ensure that the
+overall file size is a multiple of _pagesize_.
+
+# BINARY FILE FORMAT
+
+Data are encoded with respect to the little endian bytes ordering, i.e. least
+significant bytes are stored first.
+
+```
+<rnsr> ::= <pagesize> <size> 16 16
+ <padding>
+ <radiances>
+ <padding>
+
+<pagesize> ::= UINT64
+<size> ::= UINT64 # Number of items stored
+
+---
+
+<radiances> ::= <property> ... ...
+<property> ::= <wavelength> <radiance>
+<wavelength> ::= DOUBLE # in nanometer
+<radiance> ::= DOUBLE # in W/m²/sr/m
+
+---
+
+<padding> ::= [ BYTE ... ] # Ensure alignement
+```
+
+# SEE ALSO
+
+*sbuf*(5)
diff --git a/src/atmosphere/htrdr_atmosphere.c b/src/atmosphere/htrdr_atmosphere.c
@@ -24,10 +24,10 @@
#include "atmosphere/htrdr_atmosphere_sun.h"
#include "core/htrdr_buffer.h"
-#include "core/htrdr_cie_xyz.h"
#include "core/htrdr_log.h"
#include "core/htrdr_materials.h"
-#include "core/htrdr_ran_wlen.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_ran_wlen_planck.h"
#include "core/htrdr_rectangle.h"
#include <high_tune/htsky.h>
@@ -306,8 +306,8 @@ atmosphere_release(ref_T* ref)
if(cmd->ground) htrdr_atmosphere_ground_ref_put(cmd->ground);
if(cmd->mats) htrdr_materials_ref_put(cmd->mats);
if(cmd->sun) htrdr_atmosphere_sun_ref_put(cmd->sun);
- if(cmd->cie) htrdr_cie_xyz_ref_put(cmd->cie);
- if(cmd->ran_wlen) htrdr_ran_wlen_ref_put(cmd->ran_wlen);
+ if(cmd->cie) htrdr_ran_wlen_cie_xyz_ref_put(cmd->cie);
+ if(cmd->planck) htrdr_ran_wlen_planck_ref_put(cmd->planck);
if(cmd->camera) SCAM(ref_put(cmd->camera));
if(cmd->flux_map) htrdr_rectangle_ref_put(cmd->flux_map);
if(cmd->buf) htrdr_buffer_ref_put(cmd->buf);
@@ -439,7 +439,8 @@ htrdr_atmosphere_create
nintervals = compute_spectral_bands_count(cmd);
if(cmd->spectral_type == HTRDR_SPECTRAL_SW_CIE_XYZ) {
- res = htrdr_cie_xyz_create(htrdr, spectral_range, nintervals, &cmd->cie);
+ res = htrdr_ran_wlen_cie_xyz_create
+ (htrdr, spectral_range, nintervals, &cmd->cie);
if(res != RES_OK) goto error;
} else {
if(cmd->ref_temperature <= 0) {
@@ -448,8 +449,8 @@ htrdr_atmosphere_create
res = RES_BAD_ARG;
goto error;
}
- res = htrdr_ran_wlen_create
- (htrdr, spectral_range, nintervals, cmd->ref_temperature, &cmd->ran_wlen);
+ res = htrdr_ran_wlen_planck_create
+ (htrdr, spectral_range, nintervals, cmd->ref_temperature, &cmd->planck);
if(res != RES_OK) goto error;
}
diff --git a/src/atmosphere/htrdr_atmosphere_c.h b/src/atmosphere/htrdr_atmosphere_c.h
@@ -81,17 +81,17 @@ struct htsky;
struct htrdr;
struct htrdr_atmosphere_args;
struct htrdr_buffer;
-struct htrdr_cie_xyz;
struct htrdr_materials;
-struct htrdr_ran_wlen;
+struct htrdr_ran_wlen_cie_xyz;
+struct htrdr_ran_wlen_planck;
struct ssp_rng;
struct htrdr_atmosphere {
struct htrdr_atmosphere_ground* ground;
struct htrdr_atmosphere_sun* sun;
struct htrdr_materials* mats;
- struct htrdr_cie_xyz* cie;
- struct htrdr_ran_wlen* ran_wlen;
+ struct htrdr_ran_wlen_cie_xyz* cie;
+ struct htrdr_ran_wlen_planck* planck;
struct scam* camera; /* Camera */
struct htrdr_rectangle* flux_map; /* Flux map */
diff --git a/src/atmosphere/htrdr_atmosphere_draw_map.c b/src/atmosphere/htrdr_atmosphere_draw_map.c
@@ -21,11 +21,11 @@
#include "core/htrdr.h"
#include "core/htrdr_buffer.h"
-#include "core/htrdr_cie_xyz.h"
#include "core/htrdr_draw_map.h"
#include "core/htrdr_interface.h"
#include "core/htrdr_log.h"
-#include "core/htrdr_ran_wlen.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_ran_wlen_planck.h"
#include "core/htrdr_rectangle.h"
#include <high_tune/htsky.h>
@@ -160,9 +160,9 @@ draw_pixel_image
/* Sample a spectral band and a quadrature point */
switch(ichannel) {
- case 0: wlen = htrdr_cie_xyz_sample_X(cmd->cie, r0, r1, &pdf); break;
- case 1: wlen = htrdr_cie_xyz_sample_Y(cmd->cie, r0, r1, &pdf); break;
- case 2: wlen = htrdr_cie_xyz_sample_Z(cmd->cie, r0, r1, &pdf); break;
+ case 0: wlen = htrdr_ran_wlen_cie_xyz_sample_X(cmd->cie, r0, r1, &pdf); break;
+ case 1: wlen = htrdr_ran_wlen_cie_xyz_sample_Y(cmd->cie, r0, r1, &pdf); break;
+ case 2: wlen = htrdr_ran_wlen_cie_xyz_sample_Z(cmd->cie, r0, r1, &pdf); break;
default: FATAL("Unreachable code.\n"); break;
}
pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
@@ -249,7 +249,7 @@ draw_pixel_flux
r2 = ssp_rng_canonical(args->rng);
/* Sample a wavelength */
- wlen = htrdr_ran_wlen_sample(cmd->ran_wlen, r0, r1, &band_pdf);
+ wlen = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &band_pdf);
band_pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
/* Select the associated band and sample a quadrature point */
@@ -370,7 +370,7 @@ draw_pixel_xwave
r2 = ssp_rng_canonical(args->rng);
/* Sample a wavelength */
- wlen = htrdr_ran_wlen_sample(cmd->ran_wlen, r0, r1, &band_pdf);
+ wlen = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &band_pdf);
band_pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
/* Select the associated band and sample a quadrature point */
diff --git a/src/atmosphere/htrdr_atmosphere_main.c b/src/atmosphere/htrdr_atmosphere_main.c
@@ -31,7 +31,7 @@ htrdr_atmosphere_main(int argc, char** argv)
struct htrdr_atmosphere_args cmd_args = HTRDR_ATMOSPHERE_ARGS_DEFAULT;
struct htrdr* htrdr = NULL;
struct htrdr_atmosphere* cmd = NULL;
- const size_t memsz_begin = MEM_ALLOCATED_SIZE(&mem_default_allocator);
+ const size_t memsz_begin = mem_allocated_size();
size_t memsz_end;
int is_mpi_init = 0;
int err = 0;
@@ -71,7 +71,7 @@ exit:
if(cmd) htrdr_atmosphere_ref_put(cmd);
/* Check memory leaks */
- memsz_end = MEM_ALLOCATED_SIZE(&mem_default_allocator);
+ memsz_end = mem_allocated_size();
if(memsz_begin != memsz_end) {
ASSERT(memsz_end >= memsz_begin);
fprintf(stderr, HTRDR_LOG_WARNING_PREFIX"Memory leaks: %lu Bytes\n",
diff --git a/src/atmosphere/htrdr_atmosphere_sun.c b/src/atmosphere/htrdr_atmosphere_sun.c
@@ -74,6 +74,7 @@ htrdr_atmosphere_sun_create
goto error;
}
ref_init(&sun->ref);
+
sun->half_angle = 4.6524e-3;
sun->temperature = 5778;
sun->cos_half_angle = cos(sun->half_angle);
diff --git a/src/combustion/htrdr_combustion.c b/src/combustion/htrdr_combustion.c
@@ -83,7 +83,7 @@ setup_output
}
/* Setup the output name */
- str_set(&cmd->output_name, output_name);
+ res = str_set(&cmd->output_name, output_name);
if(res != RES_OK) {
htrdr_log_err(cmd->htrdr,
"Could not store the name of the output stream `%s' -- %s.\n",
@@ -589,7 +589,7 @@ htrdr_combustion_create
cmd = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cmd));
if(!cmd) {
htrdr_log_err(htrdr, "Could not allocate the htrdr_combustion data.\n");
- res = RES_BAD_ARG;
+ res = RES_MEM_ERR;
goto error;
}
ref_init(&cmd->ref);
diff --git a/src/combustion/htrdr_combustion_draw_map.c b/src/combustion/htrdr_combustion_draw_map.c
@@ -287,7 +287,7 @@ combustion_draw_map(struct htrdr_combustion* cmd)
res = htrdr_draw_map(cmd->htrdr, &args, cmd->buf);
if(res != RES_OK) goto error;
- /* No more to do on non master processes */
+ /* Nothing more to do on non master processes */
if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
/* Write buffer to output */
diff --git a/src/combustion/htrdr_combustion_laser.h b/src/combustion/htrdr_combustion_laser.h
@@ -54,21 +54,21 @@ struct htrdr_combustion_laser;
BEGIN_DECLS
-HTRDR_API res_T
+extern LOCAL_SYM res_T
htrdr_combustion_laser_create
(struct htrdr* htrdr,
const struct htrdr_combustion_laser_create_args* args,
struct htrdr_combustion_laser** laser);
-HTRDR_API void
+extern LOCAL_SYM void
htrdr_combustion_laser_ref_get
(struct htrdr_combustion_laser* laser);
-HTRDR_API void
+extern LOCAL_SYM void
htrdr_combustion_laser_ref_put
(struct htrdr_combustion_laser* laser);
-HTRDR_API void
+extern LOCAL_SYM void
htrdr_combustion_laser_trace_ray
(struct htrdr_combustion_laser* laser,
const double pos[3],
@@ -76,32 +76,32 @@ htrdr_combustion_laser_trace_ray
const double range[2],
double distance[2]);
-HTRDR_API void
+extern LOCAL_SYM void
htrdr_combustion_laser_get_mesh
(const struct htrdr_combustion_laser* laser,
/* Max distance of the laser mesh along its infinite dimension */
const double extent,
struct htrdr_combustion_laser_mesh* mesh);
-HTRDR_API void
+extern LOCAL_SYM void
htrdr_combustion_laser_get_position
(const struct htrdr_combustion_laser* laser,
double pos[3]);
-HTRDR_API void
+extern LOCAL_SYM void
htrdr_combustion_laser_get_direction
(const struct htrdr_combustion_laser* laser,
double dir[3]); /* Normalized */
-HTRDR_API double /* In W.m^2 */
+extern LOCAL_SYM double /* In W.m^2 */
htrdr_combustion_laser_get_flux_density
(const struct htrdr_combustion_laser* laser);
-HTRDR_API double /* In nm */
+extern LOCAL_SYM double /* In nm */
htrdr_combustion_laser_get_wavelength
(const struct htrdr_combustion_laser* laser);
-HTRDR_API double
+extern LOCAL_SYM double
htrdr_combustion_laser_compute_surface_plane_distance
(const struct htrdr_combustion_laser* laser,
const double pos[3]);
diff --git a/src/combustion/htrdr_combustion_main.c b/src/combustion/htrdr_combustion_main.c
@@ -32,7 +32,7 @@ htrdr_combustion_main(int argc, char** argv)
struct htrdr_combustion_args cmd_args = HTRDR_COMBUSTION_ARGS_DEFAULT;
struct htrdr* htrdr = NULL;
struct htrdr_combustion* cmd = NULL;
- const size_t memsz_begin = MEM_ALLOCATED_SIZE(&mem_default_allocator);
+ const size_t memsz_begin = mem_allocated_size();
size_t memsz_end;
int is_mpi_init = 0;
res_T res = RES_OK;
@@ -72,7 +72,7 @@ exit:
if(cmd) htrdr_combustion_ref_put(cmd);
/* Check memory leaks */
- memsz_end = MEM_ALLOCATED_SIZE(&mem_default_allocator);
+ memsz_end = mem_allocated_size();
if(memsz_begin != memsz_end) {
ASSERT(memsz_end >= memsz_begin);
fprintf(stderr, HTRDR_LOG_WARNING_PREFIX"Memory leaks: %lu Bytes\n",
diff --git a/src/commands/htrdr_planeto_cmd.c b/src/commands/htrdr_planeto_cmd.c
@@ -0,0 +1,35 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019 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/>. */
+
+#ifdef HTRDR_BUILD_PLANETO
+ #include "planeto/htrdr_planeto.h"
+#else
+ #include <stdio.h>
+#endif
+
+int
+main(int argc, char** argv)
+{
+#ifdef HTRDR_BUILD_PLANETO
+ return htrdr_planeto_main(argc, argv);
+#else
+ (void)argc, (void)argv;
+ fprintf(stderr,
+ "The htrdr-planeto command is not available in this htrdr build.\n");
+ return 1;
+#endif
+}
diff --git a/src/core/htrdr.h b/src/core/htrdr.h
@@ -67,9 +67,9 @@ htrdr_fprint_copyright(const char* cmd, FILE* stream)
{
(void)cmd;
fprintf(stream,
-"Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> <contact@meso-star.com>.\n"
-"Copyright (C) 2018, 2019, 2021 CNRS.\n"
-"Copyright (C) 2018, 2019 Université Paul Sabatier.\n");
+"Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> <contact@meso-star.com>\n"
+"Copyright (C) 2018, 2019, 2021 CNRS\n"
+"Copyright (C) 2018, 2019 Université Paul Sabatier\n");
}
static INLINE void
@@ -79,7 +79,7 @@ htrdr_fprint_license(const char* cmd, FILE* stream)
fprintf(stream,
"%s is free software released under the GNU GPL license,\n"
"version 3 or later. You are free to change or redistribute it\n"
-"under certain conditions <http://gnu.org/licenses/gpl.html>.\n",
+"under certain conditions <http://gnu.org/licenses/gpl.html>\n",
cmd);
}
diff --git a/src/core/htrdr_args.c b/src/core/htrdr_args.c
@@ -68,10 +68,10 @@ parse_fov(const char* str, double* out_fov)
fprintf(stderr, "Invalid field of view `%s'.\n", str);
return RES_BAD_ARG;
}
- if(fov <= HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MIN
+ if(fov <= HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MIN
|| fov >= HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MAX) {
fprintf(stderr, "The field of view %g is not in ]%g, %g[.\n", fov,
- HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MIN,
+ HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MIN,
HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MAX);
return RES_BAD_ARG;
}
@@ -148,13 +148,13 @@ parse_image_plane_height(const char* str, double* out_height)
res = cstr_to_double(str, &height);
if(res != RES_OK) {
- fprintf(stderr,
+ fprintf(stderr,
"Invalid height `%s' of the image plane of the orthographic camera.\n",
str);
return RES_BAD_ARG;
}
if(height <= 0) {
- fprintf(stderr,
+ fprintf(stderr,
"Invalid negative or null height of the image plane "
"of the orthographic camera.\n");
return RES_BAD_ARG;
@@ -444,8 +444,8 @@ parse_spectral_parameter(const char* str, void* ptr)
if(!strcmp(key, "cie_xyz")) {
args->spectral_type = HTRDR_SPECTRAL_SW_CIE_XYZ;
- args->wlen_range[0] = HTRDR_CIE_XYZ_RANGE_DEFAULT[0];
- args->wlen_range[1] = HTRDR_CIE_XYZ_RANGE_DEFAULT[1];
+ args->wlen_range[0] = HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[0];
+ args->wlen_range[1] = HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[1];
} else {
if(!val) {
fprintf(stderr, "Missing value to the spectral option `%s'.\n", key);
@@ -557,7 +557,7 @@ htrdr_args_camera_perspective_parse
res = cstr_parse_list(str, ':', parse_camera_perspective_parameter, cam);
if(res != RES_OK) goto error;
-
+
if(cam->focal_length < 0) {
ASSERT(cam->fov_y > 0);
res = scam_field_of_view_to_focal_length
@@ -565,7 +565,7 @@ htrdr_args_camera_perspective_parse
if(res != RES_OK) {
fprintf(stderr,
"Cannot compute the focal length from the lens radius %g "
- "and the field of view %g -- %s.\n",
+ "and the field of view %g -- %s.\n",
cam->lens_radius,
cam->fov_y,
res_to_cstr(res));
@@ -579,7 +579,7 @@ htrdr_args_camera_perspective_parse
if(res != RES_OK) {
fprintf(stderr,
"Cannot compute the field of view from the lens radius %g "
- "and the focal length %g -- %s.\n",
+ "and the focal length %g -- %s.\n",
cam->lens_radius,
cam->focal_length,
res_to_cstr(res));
@@ -608,6 +608,36 @@ error:
}
res_T
+htrdr_args_camera_perspective_check
+ (const struct htrdr_args_camera_perspective* cam)
+{
+ if(!cam) return RES_BAD_ARG;
+
+ /* Invalid Fov */
+ if(cam->fov_y <= HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MIN
+ || cam->fov_y >= HTRDR_ARGS_CAMERA_PERSPECTIVE_FOV_EXCLUSIVE_MAX) {
+ /* Is the fov defined by the focal length? */
+ if(cam->focal_length < 0) return RES_BAD_ARG;
+ }
+
+ /* Invalid focal length */
+ if(cam->focal_length <= 0) {
+ /* Is the focal length defined by the fov? */
+ if(cam->fov_y < 0) return RES_BAD_ARG;
+ }
+
+ /* Invalid lens radius */
+ if(cam->lens_radius < 0)
+ return RES_BAD_ARG;
+
+ /* Invalid focal distance */
+ if(cam->lens_radius > 0/*!pinhole*/ && cam->focal_dst < 0)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+res_T
htrdr_args_camera_orthographic_parse
(struct htrdr_args_camera_orthographic* cam,
const char* str)
@@ -631,6 +661,22 @@ htrdr_args_image_parse(struct htrdr_args_image* img, const char* str)
}
res_T
+htrdr_args_image_check(const struct htrdr_args_image* img)
+{
+ if(!img) return RES_BAD_ARG;
+
+ /* Invalid definition */
+ if(!img->definition[0] || !img->definition[1])
+ return RES_BAD_ARG;
+
+ /* Invalid number of samples per pixel */
+ if(!img->spp)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+res_T
htrdr_args_spectral_parse(struct htrdr_args_spectral* spectral, const char* str)
{
if(!spectral || !str) return RES_BAD_ARG;
diff --git a/src/core/htrdr_args.h.in b/src/core/htrdr_args.h.in
@@ -18,7 +18,7 @@
#ifndef HTRDR_ARGS_H
#define HTRDR_ARGS_H
-#include "core/htrdr_cie_xyz.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
#include "core/htrdr_spectral.h"
#include <rsys/rsys.h>
@@ -115,7 +115,7 @@ struct htrdr_args_spectral {
enum htrdr_spectral_type spectral_type;
};
#define HTRDR_ARGS_SPECTRAL_DEFAULT__ { \
- HTRDR_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
+ HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
-1, /* Reference temperature */ \
HTRDR_SPECTRAL_SW_CIE_XYZ, /* Spectral type */ \
}
@@ -133,6 +133,10 @@ htrdr_args_camera_perspective_parse
const char* str);
HTRDR_CORE_API res_T
+htrdr_args_camera_perspective_check
+ (const struct htrdr_args_camera_perspective* cam);
+
+HTRDR_CORE_API res_T
htrdr_args_camera_orthographic_parse
(struct htrdr_args_camera_orthographic* cam,
const char* str);
@@ -148,6 +152,10 @@ htrdr_args_image_parse
const char* str);
HTRDR_CORE_API res_T
+htrdr_args_image_check
+ (const struct htrdr_args_image* img);
+
+HTRDR_CORE_API res_T
htrdr_args_spectral_parse
(struct htrdr_args_spectral* spectral,
const char* str);
diff --git a/src/core/htrdr_cie_xyz.c b/src/core/htrdr_cie_xyz.c
@@ -1,389 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020, 2021 |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 "core/htrdr.h"
-#include "core/htrdr_c.h"
-#include "core/htrdr_cie_xyz.h"
-#include "core/htrdr_log.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 rcp_integral_X;
- double rcp_integral_Y;
- double rcp_integral_Z;
- double range[2]; /* Boundaries of the handled CIE XYZ color space */
- double band_len; /* Length in nanometers of a band */
-
- 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_lo + 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 {
- htrdr_log_warn(cie->htrdr,
- "%s: cannot sample a wavelength in [%g, %g[. The possible wavelengths"
- "were %g and %g.\n",
- FUNC_NAME, lambda_min, lambda_max, lambda_1, lambda_2);
- /* Arbitrarly choose the wavelength at the center of the sampled band */
- lambda = (lambda_min + lambda_max)*0.5;
- }
-
- return lambda;
-}
-
-static res_T
-setup_cie_xyz
- (struct htrdr_cie_xyz* cie,
- const char* func_name,
- 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 && nbands);
- ASSERT(cie->range[0] >= HTRDR_CIE_XYZ_RANGE_DEFAULT[0]);
- ASSERT(cie->range[1] <= HTRDR_CIE_XYZ_RANGE_DEFAULT[1]);
- ASSERT(cie->range[0] < cie->range[1]);
-
- /* Allocate and reset the memory space for the tristimulus CDF */
- #define SETUP_STIMULUS(Stimulus) { \
- 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 */
- FOR_EACH(i, 0, nbands) {
- const double lambda_lo = cie->range[0] + (double)i * cie->band_len;
- const double lambda_hi = MMIN(lambda_lo + cie->band_len, cie->range[1]);
- ASSERT(lambda_lo <= lambda_hi);
- ASSERT(lambda_lo >= cie->range[0]);
- ASSERT(lambda_hi <= cie->range[1]);
- pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
- pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931);
- pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931);
- sum[X] += pdf[X][i];
- sum[Y] += pdf[Y][i];
- sum[Z] += pdf[Z][i];
- }
- #define CHK_SUM(Sum, Range, Fit) \
- ASSERT(eq_eps(Sum, trapezoidal_integration(Range[0], Range[1], Fit), 1.e-3))
- CHK_SUM(sum[X], cie->range, fit_x_bar_1931);
- CHK_SUM(sum[Y], cie->range, fit_y_bar_1931);
- CHK_SUM(sum[Z], cie->range, fit_z_bar_1931);
- #undef CHK_SUM
- cie->rcp_integral_X = 1.0 / sum[X];
- cie->rcp_integral_Y = 1.0 / sum[Y];
- cie->rcp_integral_Z = 1.0 / sum[Z];
-
- FOR_EACH(i, 0, nbands) {
- /* Normalize the pdf */
- 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;
- struct htrdr* htrdr = 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);
- htrdr = cie->htrdr;
- MEM_RM(htrdr_get_allocator(cie->htrdr), cie);
- htrdr_ref_put(htrdr);
-}
-
-/*******************************************************************************
- * 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 bands_count, /* # bands used to discretize the CIE tristimulus */
- struct htrdr_cie_xyz** out_cie)
-{
- struct htrdr_cie_xyz* cie = NULL;
- size_t nbands = bands_count;
- res_T res = RES_OK;
- ASSERT(htrdr && range && nbands && out_cie);
-
- cie = MEM_CALLOC(htrdr_get_allocator(htrdr), 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);
- darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_X);
- darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Y);
- darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Z);
- cie->range[0] = range[0];
- cie->range[1] = range[1];
- htrdr_ref_get(htrdr);
- cie->htrdr = htrdr;
-
- cie->band_len = (range[1] - range[0]) / (double)nbands;
-
- res = setup_cie_xyz(cie, FUNC_NAME, nbands);
- if(res != RES_OK) goto error;
-
- htrdr_log(htrdr, "CIE XYZ spectral 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,
- double* pdf)
-{
- const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
- darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
- if(pdf) *pdf = cie->rcp_integral_X;
- return wlen;
-}
-
-double
-htrdr_cie_xyz_sample_Y
- (struct htrdr_cie_xyz* cie,
- const double r0,
- const double r1,
- double* pdf)
-{
- const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
- darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
- if(pdf) *pdf = cie->rcp_integral_Y;
- return wlen;
-}
-
-double
-htrdr_cie_xyz_sample_Z
- (struct htrdr_cie_xyz* cie,
- const double r0,
- const double r1,
- double* pdf)
-{
- const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
- darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
- if(pdf) *pdf = cie->rcp_integral_Z;
- return wlen;
-}
-
diff --git a/src/core/htrdr_cie_xyz.h b/src/core/htrdr_cie_xyz.h
@@ -1,72 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020, 2021 |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 "core/htrdr.h"
-#include <rsys/rsys.h>
-
-/* Wavelength boundaries of the CIE XYZ color space in nanometers */
-#define HTRDR_CIE_XYZ_RANGE_DEFAULT__ {380, 780}
-static const double HTRDR_CIE_XYZ_RANGE_DEFAULT[2] =
- HTRDR_CIE_XYZ_RANGE_DEFAULT__;
-
-/* Forward declarations */
-struct htrdr;
-struct htrdr_cie_xyz;
-
-BEGIN_DECLS
-
-HTRDR_CORE_API 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);
-
-HTRDR_CORE_API void
-htrdr_cie_xyz_ref_get
- (struct htrdr_cie_xyz* cie);
-
-HTRDR_CORE_API void
-htrdr_cie_xyz_ref_put
- (struct htrdr_cie_xyz* cie);
-
-/* Return a wavelength in nanometer */
-HTRDR_CORE_API double
-htrdr_cie_xyz_sample_X
- (struct htrdr_cie_xyz* cie,
- const double r0, const double r1, /* Canonical numbers in [0, 1[ */
- double* pdf); /* In nm^-1. May be NULL */
-
-/* Return a wavelength in nanometer */
-HTRDR_CORE_API double
-htrdr_cie_xyz_sample_Y
- (struct htrdr_cie_xyz* cie,
- const double r0, const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* In nm^-1. May be NULL */
-
-/* Return a wavelength in nanometer */
-HTRDR_CORE_API double
-htrdr_cie_xyz_sample_Z
- (struct htrdr_cie_xyz* cie,
- const double r0, const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* In nm^-1. May be NULL */
-
-END_DECLS
-
-#endif /* HTRDR_cie_xyz_H */
-
diff --git a/src/core/htrdr_ran_wlen.c b/src/core/htrdr_ran_wlen.c
@@ -1,364 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019, 2021 CNRS
- * Copyright (C) 2018, 2019 Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#define _POSIX_C_SOURCE 200112L /* nextafter */
-
-#include "core/htrdr.h"
-#include "core/htrdr_c.h"
-#include "core/htrdr_log.h"
-#include "core/htrdr_ran_wlen.h"
-#include "core/htrdr_spectral.h"
-
-#include <rsys/algorithm.h>
-#include <rsys/cstr.h>
-#include <rsys/dynamic_array_double.h>
-#include <rsys/mem_allocator.h>
-#include <rsys/ref_count.h>
-
-#include <math.h> /* nextafter */
-
-struct htrdr_ran_wlen {
- struct darray_double pdf;
- struct darray_double cdf;
- double range[2]; /* Boundaries of the spectral integration interval */
- double band_len; /* Length in nanometers of a band */
- double ref_temperature; /* In Kelvin */
- size_t nbands; /* # bands */
-
- struct htrdr* htrdr;
- ref_T ref;
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static res_T
-setup_wlen_ran_cdf
- (struct htrdr_ran_wlen* wlen_ran,
- const char* func_name)
-{
- double* pdf = NULL;
- double* cdf = NULL;
- double sum = 0;
- size_t i;
- res_T res = RES_OK;
- ASSERT(wlen_ran && func_name && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
-
- res = darray_double_resize(&wlen_ran->cdf, wlen_ran->nbands);
- if(res != RES_OK) {
- htrdr_log_err(wlen_ran->htrdr,
- "%s: Error allocating the CDF of the spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
- res = darray_double_resize(&wlen_ran->pdf, wlen_ran->nbands);
- if(res != RES_OK) {
- htrdr_log_err(wlen_ran->htrdr,
- "%s: Error allocating the pdf of the spectral bands -- %s.\n",
- func_name, res_to_cstr(res));
- goto error;
- }
-
- cdf = darray_double_data_get(&wlen_ran->cdf);
- pdf = darray_double_data_get(&wlen_ran->pdf); /* Now save pdf to correct weight */
-
- htrdr_log(wlen_ran->htrdr,
- "Number of bands of the spectrum cumulative: %lu\n",
- (unsigned long)wlen_ran->nbands);
-
- /* Compute the *unnormalized* probability to sample a small band */
- FOR_EACH(i, 0, wlen_ran->nbands) {
- double lambda_lo = wlen_ran->range[0] + (double)i * wlen_ran->band_len;
- double lambda_hi = MMIN(lambda_lo + wlen_ran->band_len, wlen_ran->range[1]);
- ASSERT(lambda_lo<= lambda_hi);
- ASSERT(lambda_lo > wlen_ran->range[0]
- || eq_eps(lambda_lo, wlen_ran->range[0], 1.e-6));
- ASSERT(lambda_lo < wlen_ran->range[1]
- || eq_eps(lambda_lo, wlen_ran->range[1], 1.e-6));
-
- /* Convert from nanometer to meter */
- lambda_lo *= 1.e-9;
- lambda_hi *= 1.e-9;
-
- /* Compute the probability of the current band */
- pdf[i] = htrdr_blackbody_fraction
- (lambda_lo, lambda_hi, wlen_ran->ref_temperature);
-
- /* Update the norm */
- sum += pdf[i];
- }
-
- /* Compute the cumulative of the previously computed probabilities */
- FOR_EACH(i, 0, wlen_ran->nbands) {
- /* Normalize the probability */
- pdf[i] /= sum;
-
- /* Setup the cumulative */
- if(i == 0) {
- cdf[i] = pdf[i];
- } else {
- cdf[i] = pdf[i] + cdf[i-1];
- ASSERT(cdf[i] >= cdf[i-1]);
- }
- }
- ASSERT(eq_eps(cdf[wlen_ran->nbands-1], 1, 1.e-6));
- cdf[wlen_ran->nbands - 1] = 1.0; /* Handle numerical issue */
-
-exit:
- return res;
-error:
- darray_double_clear(&wlen_ran->cdf);
- darray_double_clear(&wlen_ran->pdf);
- goto exit;
-}
-
-static double
-wlen_ran_sample_continue
- (const struct htrdr_ran_wlen* wlen_ran,
- const double r,
- const double range[2], /* In nanometer */
- const char* func_name,
- double* pdf)
-{
- /* Numerical parameters of the dichotomy search */
- const size_t MAX_ITER = 100;
- const double EPSILON_LAMBDA_M = 1e-15;
- const double EPSILON_BF = 1e-6;
-
- /* Local variables */
- double bf = 0;
- double bf_prev = 0;
- double bf_min_max = 0;
- double lambda_m = 0;
- double lambda_m_prev = 0;
- double lambda_m_min = 0;
- double lambda_m_max = 0;
- double range_m[2] = {0, 0};
- size_t i;
-
- /* Check precondition */
- ASSERT(wlen_ran && func_name);
- ASSERT(range && range[0] < range[1]);
- ASSERT(0 <= r && r < 1);
-
- /* Convert the wavelength range in meters */
- range_m[0] = range[0] * 1.0e-9;
- range_m[1] = range[1] * 1.0e-9;
-
- /* Setup the dichotomy search */
- lambda_m_min = range_m[0];
- lambda_m_max = range_m[1];
- bf_min_max = htrdr_blackbody_fraction
- (range_m[0], range_m[1], wlen_ran->ref_temperature);
-
- /* Numerically search the lambda corresponding to the submitted canonical
- * number */
- FOR_EACH(i, 0, MAX_ITER) {
- double r_test;
- lambda_m = (lambda_m_min + lambda_m_max) * 0.5;
- bf = htrdr_blackbody_fraction
- (range_m[0], lambda_m, wlen_ran->ref_temperature);
-
- r_test = bf / bf_min_max;
- if(r_test < r) {
- lambda_m_min = lambda_m;
- } else {
- lambda_m_max = lambda_m;
- }
-
- if(fabs(lambda_m_prev - lambda_m) < EPSILON_LAMBDA_M
- || fabs(bf_prev - bf) < EPSILON_BF)
- break;
-
- lambda_m_prev = lambda_m;
- bf_prev = bf;
- }
- if(i >= MAX_ITER) {
- htrdr_log_warn(wlen_ran->htrdr,
- "%s: could not sample a wavelength in the range [%g, %g] nanometers "
- "for the reference temperature %g Kelvin.\n",
- func_name, SPLIT2(range), wlen_ran->ref_temperature);
- }
-
- if(pdf) {
- const double Tref = wlen_ran->ref_temperature;
- const double B_lambda = htrdr_planck(lambda_m, lambda_m, Tref);
- const double B_mean = htrdr_planck(range_m[0], range_m[1], Tref);
- *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0]));
- *pdf *= 1.e-9; /* Transform from m^-1 to nm^-1 */
- }
-
- return lambda_m * 1.e9; /* Convert in nanometers */
-}
-
-static double
-wlen_ran_sample_discrete
- (const struct htrdr_ran_wlen* wlen_ran,
- const double r0,
- const double r1,
- const char* func_name,
- double* pdf)
-{
- const double* cdf = NULL;
- const double* find = NULL;
- double r0_next = nextafter(r0, DBL_MAX);
- double band_range[2];
- double lambda = 0;
- double pdf_continue = 0;
- double pdf_band = 0;
- size_t cdf_length = 0;
- size_t i;
- ASSERT(wlen_ran && wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE);
- ASSERT(0 <= r0 && r0 < 1);
- ASSERT(0 <= r1 && r1 < 1);
- (void)func_name;
- (void)pdf_band;
-
- cdf = darray_double_cdata_get(&wlen_ran->cdf);
- cdf_length = darray_double_size_get(&wlen_ran->cdf);
- ASSERT(cdf_length > 0);
-
- /* Use r_next rather than r0 in order to find the first entry that is not less
- * than *or equal* to r0 */
- find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
- ASSERT(find);
-
- i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r0 && (!i || cdf[i-1] <= r0));
-
- band_range[0] = wlen_ran->range[0] + (double)i*wlen_ran->band_len;
- band_range[1] = band_range[0] + wlen_ran->band_len;
-
- /* Fetch the pdf of the sampled band */
- pdf_band = darray_double_cdata_get(&wlen_ran->pdf)[i];
-
- /* Uniformly sample a wavelength in the sampled band */
- lambda = band_range[0] + (band_range[1] - band_range[0]) * r1;
- pdf_continue = 1.0 / (band_range[1] - band_range[0]);
-
- if(pdf) {
- *pdf = pdf_band * pdf_continue;
- }
-
- return lambda;
-}
-
-static void
-release_wlen_ran(ref_T* ref)
-{
- struct htrdr_ran_wlen* wlen_ran = NULL;
- struct htrdr* htrdr = NULL;
- ASSERT(ref);
- wlen_ran = CONTAINER_OF(ref, struct htrdr_ran_wlen, ref);
- darray_double_release(&wlen_ran->cdf);
- darray_double_release(&wlen_ran->pdf);
- htrdr = wlen_ran->htrdr;
- MEM_RM(htrdr_get_allocator(htrdr), wlen_ran);
- htrdr_ref_put(wlen_ran->htrdr);
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-res_T
-htrdr_ran_wlen_create
- (struct htrdr* htrdr,
- /* range must be included in [200,1000] nm for shortwave or in [1000,100000]
- * nanometers for longwave (thermal) */
- const double range[2],
- const size_t nbands, /* # bands used to discretized CDF */
- const double ref_temperature,
- struct htrdr_ran_wlen** out_wlen_ran)
-{
- struct htrdr_ran_wlen* wlen_ran = NULL;
- res_T res = RES_OK;
- ASSERT(htrdr && range && out_wlen_ran && ref_temperature > 0);
- ASSERT(ref_temperature > 0);
- ASSERT(range[0] <= range[1]);
-
- wlen_ran = MEM_CALLOC(htrdr->allocator, 1, sizeof(*wlen_ran));
- if(!wlen_ran) {
- res = RES_MEM_ERR;
- htrdr_log_err(htrdr,
- "%s: could not allocate longwave random variate data structure -- %s.\n",
- FUNC_NAME, res_to_cstr(res));
- goto error;
- }
- ref_init(&wlen_ran->ref);
- darray_double_init(htrdr->allocator, &wlen_ran->cdf);
- darray_double_init(htrdr->allocator, &wlen_ran->pdf);
- htrdr_ref_get(htrdr);
- wlen_ran->htrdr = htrdr;
-
- wlen_ran->range[0] = range[0];
- wlen_ran->range[1] = range[1];
- wlen_ran->ref_temperature = ref_temperature;
- wlen_ran->nbands = nbands;
-
- if(nbands == HTRDR_WLEN_RAN_CONTINUE) {
- wlen_ran->band_len = 0;
- } else {
- wlen_ran->band_len = (range[1] - range[0]) / (double)wlen_ran->nbands;
-
- res = setup_wlen_ran_cdf(wlen_ran, FUNC_NAME);
- if(res != RES_OK) goto error;
- }
-
- htrdr_log(htrdr, "Spectral interval defined on [%g, %g] nanometers.\n",
- range[0], range[1]);
-
-exit:
- *out_wlen_ran = wlen_ran;
- return res;
-error:
- if(wlen_ran) htrdr_ran_wlen_ref_put(wlen_ran);
- goto exit;
-}
-
-void
-htrdr_ran_wlen_ref_get(struct htrdr_ran_wlen* wlen_ran)
-{
- ASSERT(wlen_ran);
- ref_get(&wlen_ran->ref);
-}
-
-void
-htrdr_ran_wlen_ref_put(struct htrdr_ran_wlen* wlen_ran)
-{
- ASSERT(wlen_ran);
- ref_put(&wlen_ran->ref, release_wlen_ran);
-}
-
-double
-htrdr_ran_wlen_sample
- (const struct htrdr_ran_wlen* wlen_ran,
- const double r0,
- const double r1,
- double* pdf)
-{
- ASSERT(wlen_ran);
- if(wlen_ran->nbands != HTRDR_WLEN_RAN_CONTINUE) { /* Discrete */
- return wlen_ran_sample_discrete(wlen_ran, r0, r1, FUNC_NAME, pdf);
- } else if(eq_eps(wlen_ran->range[0], wlen_ran->range[1], 1.e-6)) {
- if(pdf) *pdf = 1;
- return wlen_ran->range[0];
- } else { /* Continue */
- return wlen_ran_sample_continue
- (wlen_ran, r0, wlen_ran->range, FUNC_NAME, pdf);
- }
-}
-
diff --git a/src/core/htrdr_ran_wlen.h b/src/core/htrdr_ran_wlen.h
@@ -1,61 +0,0 @@
-/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
- * Copyright (C) 2018, 2019, 2021 CNRS
- * Copyright (C) 2018, 2019 Université Paul Sabatier
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTRDR_RAN_WLEN_H
-#define HTRDR_RAN_WLEN_H
-
-#include "core/htrdr.h"
-#include <rsys/rsys.h>
-
-#define HTRDR_WLEN_RAN_CONTINUE 0
-
-/* Forward declarations */
-struct htrdr;
-struct htrdr_ran_wlen;
-
-BEGIN_DECLS
-
-HTRDR_CORE_API res_T
-htrdr_ran_wlen_create
- (struct htrdr* htrdr,
- const double range[2],
- /* # bands used to discretisze the spectral domain. HTRDR_WLEN_RAN_CONTINUE
- * <=> no discretisation */
- const size_t nbands, /* Hint on #bands used to discretised th CDF */
- const double ref_temperature, /* Reference temperature */
- struct htrdr_ran_wlen** wlen_ran);
-
-HTRDR_CORE_API void
-htrdr_ran_wlen_ref_get
- (struct htrdr_ran_wlen* wlen_ran);
-
-HTRDR_CORE_API void
-htrdr_ran_wlen_ref_put
- (struct htrdr_ran_wlen* wlen_ran);
-
-/* Return a wavelength in nanometer */
-HTRDR_CORE_API double
-htrdr_ran_wlen_sample
- (const struct htrdr_ran_wlen* wlen_ran,
- const double r0, /* Canonical number in [0, 1[ */
- const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* In nm^-1. May be NULL */
-
-END_DECLS
-
-#endif /* HTRDR_RAN_WLEN_H */
-
diff --git a/src/core/htrdr_ran_wlen_cie_xyz.c b/src/core/htrdr_ran_wlen_cie_xyz.c
@@ -0,0 +1,391 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "core/htrdr.h"
+#include "core/htrdr_c.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_log.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <math.h> /* nextafter */
+
+struct htrdr_ran_wlen_cie_xyz {
+ struct darray_double cdf_X;
+ struct darray_double cdf_Y;
+ struct darray_double cdf_Z;
+ double rcp_integral_X;
+ double rcp_integral_Y;
+ double rcp_integral_Z;
+ double range[2]; /* Boundaries of the handled CIE XYZ color space */
+ double band_len; /* Length in nanometers of a band */
+
+ 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_lo + 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_ran_wlen_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 {
+ htrdr_log_warn(cie->htrdr,
+ "%s: cannot sample a wavelength in [%g, %g[. The possible wavelengths"
+ "were %g and %g.\n",
+ FUNC_NAME, lambda_min, lambda_max, lambda_1, lambda_2);
+ /* Arbitrarly choose the wavelength at the center of the sampled band */
+ lambda = (lambda_min + lambda_max)*0.5;
+ }
+
+ return lambda;
+}
+
+static res_T
+setup_cie_xyz
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const char* func_name,
+ 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 && nbands);
+ ASSERT(cie->range[0] >= HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[0]);
+ ASSERT(cie->range[1] <= HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[1]);
+ ASSERT(cie->range[0] < cie->range[1]);
+
+ /* Allocate and reset the memory space for the tristimulus CDF */
+ #define SETUP_STIMULUS(Stimulus) { \
+ 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 */
+ FOR_EACH(i, 0, nbands) {
+ const double lambda_lo = cie->range[0] + (double)i * cie->band_len;
+ const double lambda_hi = MMIN(lambda_lo + cie->band_len, cie->range[1]);
+ ASSERT(lambda_lo <= lambda_hi);
+ ASSERT(lambda_lo >= cie->range[0]);
+ ASSERT(lambda_hi <= cie->range[1]);
+ pdf[X][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_x_bar_1931);
+ pdf[Y][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_y_bar_1931);
+ pdf[Z][i] = trapezoidal_integration(lambda_lo, lambda_hi, fit_z_bar_1931);
+ sum[X] += pdf[X][i];
+ sum[Y] += pdf[Y][i];
+ sum[Z] += pdf[Z][i];
+ }
+ #define CHK_SUM(Sum, Range, Fit) \
+ ASSERT(eq_eps(Sum, trapezoidal_integration(Range[0], Range[1], Fit), 1.e-3))
+ CHK_SUM(sum[X], cie->range, fit_x_bar_1931);
+ CHK_SUM(sum[Y], cie->range, fit_y_bar_1931);
+ CHK_SUM(sum[Z], cie->range, fit_z_bar_1931);
+ #undef CHK_SUM
+ cie->rcp_integral_X = 1.0 / sum[X];
+ cie->rcp_integral_Y = 1.0 / sum[Y];
+ cie->rcp_integral_Z = 1.0 / sum[Z];
+
+ FOR_EACH(i, 0, nbands) {
+ /* Normalize the pdf */
+ 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_ran_wlen_cie_xyz* cie = NULL;
+ struct htrdr* htrdr = NULL;
+ ASSERT(ref);
+ cie = CONTAINER_OF(ref, struct htrdr_ran_wlen_cie_xyz, ref);
+ darray_double_release(&cie->cdf_X);
+ darray_double_release(&cie->cdf_Y);
+ darray_double_release(&cie->cdf_Z);
+ htrdr = cie->htrdr;
+ MEM_RM(htrdr_get_allocator(cie->htrdr), cie);
+ htrdr_ref_put(htrdr);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_ran_wlen_cie_xyz_create
+ (struct htrdr* htrdr,
+ const double range[2], /* Must be included in [380, 780] nanometers */
+ const size_t bands_count, /* # bands used to discretize the CIE tristimulus */
+ struct htrdr_ran_wlen_cie_xyz** out_cie)
+{
+ struct htrdr_ran_wlen_cie_xyz* cie = NULL;
+ size_t nbands = bands_count;
+ res_T res = RES_OK;
+ ASSERT(htrdr && range && nbands && out_cie);
+
+ cie = MEM_CALLOC(htrdr_get_allocator(htrdr), 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);
+ darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_X);
+ darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Y);
+ darray_double_init(htrdr_get_allocator(htrdr), &cie->cdf_Z);
+ cie->range[0] = range[0];
+ cie->range[1] = range[1];
+ htrdr_ref_get(htrdr);
+ cie->htrdr = htrdr;
+
+ cie->band_len = (range[1] - range[0]) / (double)nbands;
+
+ res = setup_cie_xyz(cie, FUNC_NAME, nbands);
+ if(res != RES_OK) goto error;
+
+ htrdr_log(htrdr, "CIE XYZ spectral interval defined on [%g, %g] nanometers.\n",
+ range[0], range[1]);
+
+exit:
+ *out_cie = cie;
+ return res;
+error:
+ if(cie) htrdr_ran_wlen_cie_xyz_ref_put(cie);
+ goto exit;
+}
+
+void
+htrdr_ran_wlen_cie_xyz_ref_get(struct htrdr_ran_wlen_cie_xyz* cie)
+{
+ ASSERT(cie);
+ ref_get(&cie->ref);
+}
+
+void
+htrdr_ran_wlen_cie_xyz_ref_put(struct htrdr_ran_wlen_cie_xyz* cie)
+{
+ ASSERT(cie);
+ ref_put(&cie->ref, release_cie_xyz);
+}
+
+double
+htrdr_ran_wlen_cie_xyz_sample_X
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_X),
+ darray_double_size_get(&cie->cdf_X), fit_x_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_X;
+ return wlen;
+}
+
+double
+htrdr_ran_wlen_cie_xyz_sample_Y
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Y),
+ darray_double_size_get(&cie->cdf_Y), fit_y_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_Y;
+ return wlen;
+}
+
+double
+htrdr_ran_wlen_cie_xyz_sample_Z
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ const double wlen = sample_cie_xyz(cie, darray_double_cdata_get(&cie->cdf_Z),
+ darray_double_size_get(&cie->cdf_Z), fit_z_bar_1931, r0, r1);
+ if(pdf) *pdf = cie->rcp_integral_Z;
+ return wlen;
+}
+
diff --git a/src/core/htrdr_ran_wlen_cie_xyz.h b/src/core/htrdr_ran_wlen_cie_xyz.h
@@ -0,0 +1,74 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_RAN_WLEN_CIE_XYZ_H
+#define HTRDR_RAN_WLEN_CIE_XYZ_H
+
+#include "core/htrdr.h"
+#include <rsys/rsys.h>
+
+/* Wavelength boundaries of the CIE XYZ color space in nanometers */
+#define HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT__ {380, 780}
+static const double HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[2] =
+ HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT__;
+
+/* Forward declarations */
+struct htrdr;
+struct htrdr_ran_wlen_cie_xyz;
+
+BEGIN_DECLS
+
+HTRDR_CORE_API res_T
+htrdr_ran_wlen_cie_xyz_create
+ (struct htrdr* htrdr,
+ const double range[2], /* Must be included in [380, 780] nanometers */
+ const size_t nbands, /* # bands used to discretize the CIE tristimulus */
+ struct htrdr_ran_wlen_cie_xyz** cie);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_cie_xyz_ref_get
+ (struct htrdr_ran_wlen_cie_xyz* cie);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_cie_xyz_ref_put
+ (struct htrdr_ran_wlen_cie_xyz* cie);
+
+/* Return a wavelength in nanometer */
+HTRDR_CORE_API double
+htrdr_ran_wlen_cie_xyz_sample_X
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0, const double r1, /* Canonical numbers in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
+
+/* Return a wavelength in nanometer */
+HTRDR_CORE_API double
+htrdr_ran_wlen_cie_xyz_sample_Y
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
+
+/* Return a wavelength in nanometer */
+HTRDR_CORE_API double
+htrdr_ran_wlen_cie_xyz_sample_Z
+ (struct htrdr_ran_wlen_cie_xyz* cie,
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
+
+END_DECLS
+
+#endif /* HTRDR_RAN_WLEN_CIE_XYZ_H */
+
diff --git a/src/core/htrdr_ran_wlen_discrete.c b/src/core/htrdr_ran_wlen_discrete.c
@@ -0,0 +1,306 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "core/htrdr_c.h"
+#include "core/htrdr_log.h"
+#include "core/htrdr_ran_wlen_discrete.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <math.h>
+
+struct htrdr_ran_wlen_discrete {
+ struct darray_double cumul;
+ struct darray_double proba;
+ struct darray_double radia; /* In W/m²/sr/m */
+ struct darray_double wlens; /* In nm */
+ double range[2]; /* Boundaries of the spectral integration interval in nm */
+ size_t nbands; /* #bands */
+
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE res_T
+check_htrdr_ran_wlen_discrete_create_args
+ (const struct htrdr_ran_wlen_discrete_create_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Invalid number of wavelength */
+ if(!args->nwavelengths)
+ return RES_BAD_ARG;
+
+ /* Invalid functor */
+ if(!args->get)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+static res_T
+setup_per_wlen_radiance
+ (struct htrdr_ran_wlen_discrete* ran,
+ const struct htrdr_ran_wlen_discrete_create_args* args)
+{
+ double* wlens = NULL;
+ double* radia = NULL;
+ size_t iwlen = 0;
+ res_T res = RES_OK;
+ ASSERT(ran && args);
+
+ res = darray_double_resize(&ran->wlens, args->nwavelengths);
+ if(res != RES_OK) {
+ htrdr_log_err(ran->htrdr,
+ "Error allocating discrete wavelength distribution wavelength list.\n");
+ goto error;
+ }
+ res = darray_double_resize(&ran->radia, args->nwavelengths);
+ if(res != RES_OK) {
+ htrdr_log_err(ran->htrdr,
+ "Error allocating discrete wavelength distribution radiance list.\n");
+ goto error;
+ }
+
+ wlens = darray_double_data_get(&ran->wlens);
+ radia = darray_double_data_get(&ran->radia);
+
+ /* Store the discrete values */
+ FOR_EACH(iwlen, 0, args->nwavelengths) {
+ args->get(args->context, iwlen, wlens+iwlen, radia+iwlen);
+
+ if(iwlen > 0 && wlens[iwlen] <= wlens[iwlen-1]) {
+ htrdr_log_err(ran->htrdr,
+ "Failed to calculate discrete luminance distribution probabilities. "
+ "Wavelengths are not sorted in ascending order.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ /* Setup the spectral range */
+ ran->range[0] = wlens[0];
+ ran->range[1] = wlens[args->nwavelengths-1];
+
+exit:
+ return res;
+error:
+ darray_double_clear(&ran->wlens);
+ darray_double_clear(&ran->radia);
+ goto exit;
+}
+
+static res_T
+setup_distribution
+ (struct htrdr_ran_wlen_discrete* ran,
+ const struct htrdr_ran_wlen_discrete_create_args* args)
+{
+ double* cumul = NULL;
+ double* proba = NULL;
+ double sum = 0;
+ size_t iband;
+ res_T res = RES_OK;
+ ASSERT(ran && check_htrdr_ran_wlen_discrete_create_args(args) == RES_OK);
+ ASSERT(ran->nbands >= 1); /* At least one band */
+
+ res = darray_double_resize(&ran->cumul, ran->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(ran->htrdr,
+ "Error allocating the cumulative discrete wavelength distribution.\n");
+ goto error;
+ }
+ res = darray_double_resize(&ran->proba, ran->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(ran->htrdr,
+ "Error allocating the discrete wavelength distribution probabilities.\n");
+ goto error;
+ }
+
+ cumul = darray_double_data_get(&ran->cumul);
+ proba = darray_double_data_get(&ran->proba);
+
+ /* Compute the unormalized probabilities to sample a band */
+ FOR_EACH(iband, 0, ran->nbands) {
+ const size_t iw0 = iband+0;
+ const size_t iw1 = iband+1;
+ double w0, L0;
+ double w1, L1;
+ double area;
+
+ args->get(args->context, iw0, &w0, &L0);
+ args->get(args->context, iw1, &w1, &L1);
+ ASSERT(w0 < w1);
+
+ area = (L0 + L1) * (w1-w0) * 0.5;
+
+ proba[iband] = area;
+ sum += area;
+ }
+
+ htrdr_log(ran->htrdr, "Discrete radiance integral = %g W/m²/sr\n", sum);
+
+ /* Normalize the probabilities and setup the cumulative */
+ FOR_EACH(iband, 0, ran->nbands) {
+ proba[iband] /= sum;
+ if(iband == 0) {
+ cumul[iband] = proba[iband];
+ } else {
+ cumul[iband] = proba[iband] + cumul[iband-1];
+ ASSERT(cumul[iband] >= cumul[iband-1]);
+ }
+ }
+ ASSERT(eq_eps(cumul[ran->nbands-1], 1, 1e-6));
+ cumul[ran->nbands-1] = 1.0; /* Fix numerical imprecision */
+
+exit:
+ return res;
+error:
+ darray_double_clear(&ran->proba);
+ darray_double_clear(&ran->cumul);
+ goto exit;
+}
+
+static void
+release_discrete(ref_T* ref)
+{
+ struct htrdr_ran_wlen_discrete* ran = NULL;
+ struct htrdr* htrdr = NULL;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, struct htrdr_ran_wlen_discrete, ref);
+ darray_double_release(&ran->cumul);
+ darray_double_release(&ran->proba);
+ darray_double_release(&ran->radia);
+ darray_double_release(&ran->wlens);
+ htrdr = ran->htrdr;
+ MEM_RM(htrdr_get_allocator(htrdr), ran);
+ htrdr_ref_put(htrdr);
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+htrdr_ran_wlen_discrete_create
+ (struct htrdr* htrdr,
+ const struct htrdr_ran_wlen_discrete_create_args* args,
+ struct htrdr_ran_wlen_discrete** out_ran)
+{
+ struct htrdr_ran_wlen_discrete* ran = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && args && out_ran);
+ ASSERT(check_htrdr_ran_wlen_discrete_create_args(args) == RES_OK);
+
+ ran = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*ran));
+ if(!ran) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: error allocating discrete wavelength distribution\n",
+ FUNC_NAME);
+ goto error;
+ }
+
+ ref_init(&ran->ref);
+ darray_double_init(htrdr_get_allocator(htrdr), &ran->cumul);
+ darray_double_init(htrdr_get_allocator(htrdr), &ran->proba);
+ darray_double_init(htrdr_get_allocator(htrdr), &ran->radia);
+ darray_double_init(htrdr_get_allocator(htrdr), &ran->wlens);
+ htrdr_ref_get(htrdr);
+ ran->htrdr = htrdr;
+
+ ran->nbands = args->nwavelengths - 1;
+
+ res = setup_per_wlen_radiance(ran, args);
+ if(res != RES_OK) goto error;
+ res = setup_distribution(ran, args);
+ if(res != RES_OK) goto error;
+
+exit:
+ *out_ran = ran;
+ return res;
+error:
+ if(ran) { htrdr_ran_wlen_discrete_ref_put(ran); ran = NULL; }
+ goto exit;
+}
+
+void
+htrdr_ran_wlen_discrete_ref_get(struct htrdr_ran_wlen_discrete* ran)
+{
+ ASSERT(ran);
+ ref_get(&ran->ref);
+}
+
+void
+htrdr_ran_wlen_discrete_ref_put(struct htrdr_ran_wlen_discrete* ran)
+{
+ ASSERT(ran);
+ ref_put(&ran->ref, release_discrete);
+}
+
+double
+htrdr_ran_wlen_discrete_sample
+ (struct htrdr_ran_wlen_discrete* ran,
+ const double r0,
+ const double r1,
+ double* pdf) /* In nm⁻¹ */
+{
+ double* find = NULL;
+ const double* proba = NULL;
+ const double* cumul = NULL;
+ const double* wlens = NULL;
+ double w0, w1; /* Wavelengths of the sampled band in nm */
+ double lambda; /* Sampled wavelength in nm */
+ size_t iband = 0;
+ double r0_next = nextafter(r0, DBL_MAX);
+ ASSERT(0 <= r0 && r0 < 1);
+ ASSERT(0 <= r1 && r1 < 1);
+
+ cumul = darray_double_cdata_get(&ran->cumul);
+ proba = darray_double_cdata_get(&ran->proba);
+ wlens = darray_double_cdata_get(&ran->wlens);
+
+ /* Sample a band. Use r0_next rather than r0 to find the first entry that is
+ * not less than *or equal* to r0 */
+ find = search_lower_bound
+ (&r0_next, cumul, ran->nbands, sizeof(double), cmp_dbl);
+ ASSERT(find);
+
+ iband = (size_t)(find - cumul);
+ ASSERT(iband < ran->nbands);
+ ASSERT(cumul[iband] > r0 && (!iband || cumul[iband-1] <= r0));
+
+ /* Retrieve the boundaries of the sampled band */
+ w0 = wlens[iband+0];
+ w1 = wlens[iband+1];
+
+ /* Uniformely sample the wavelength in [w0, w1[ */
+ lambda = w0 + r1 * (w1 - w0);
+
+ if(pdf) {
+ const double pdf_wlen = 1.f / (w1-w0);
+ *pdf = proba[iband] * pdf_wlen;
+ }
+
+ return lambda;
+}
diff --git a/src/core/htrdr_ran_wlen_discrete.h b/src/core/htrdr_ran_wlen_discrete.h
@@ -0,0 +1,66 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_RAN_WLEN_DISCRETE_H
+#define HTRDR_RAN_WLEN_DISCRETE_H
+
+#include "core/htrdr.h"
+#include <rsys/rsys.h>
+
+struct htrdr_ran_wlen_discrete_create_args {
+ void (*get)
+ (void* ctx,
+ const size_t i,
+ double* wlen, /* In nanometer */
+ double* radiance); /* In W/m²/sr/m */
+ size_t nwavelengths;
+ void* context; /* User defined data */
+};
+#define HTRDR_RAN_WLEN_DISCRETE_CREATE_ARGS_NULL__ {NULL, 0, NULL}
+static const struct htrdr_ran_wlen_discrete_create_args
+HTRDR_RAN_WLEN_DISCRETE_CREATE_ARGS_NULL =
+ HTRDR_RAN_WLEN_DISCRETE_CREATE_ARGS_NULL__;
+
+/* Forward declarations */
+struct htrdr;
+struct htrdr_ran_wlen_discrete;
+
+BEGIN_DECLS
+
+HTRDR_CORE_API res_T
+htrdr_ran_wlen_discrete_create
+ (struct htrdr* htrdr,
+ const struct htrdr_ran_wlen_discrete_create_args* args,
+ struct htrdr_ran_wlen_discrete** ran);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_discrete_ref_get
+ (struct htrdr_ran_wlen_discrete* ran);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_discrete_ref_put
+ (struct htrdr_ran_wlen_discrete* ran);
+
+HTRDR_CORE_API double /* wavelength in nanometer */
+htrdr_ran_wlen_discrete_sample
+ (struct htrdr_ran_wlen_discrete* ran,
+ const double r0, const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm⁻¹ */
+
+END_DECLS
+
+#endif /* HTRDR_RAN_WLEN_DISCRETE_H */
diff --git a/src/core/htrdr_ran_wlen_planck.c b/src/core/htrdr_ran_wlen_planck.c
@@ -0,0 +1,368 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019 Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* nextafter */
+
+#include "core/htrdr.h"
+#include "core/htrdr_c.h"
+#include "core/htrdr_log.h"
+#include "core/htrdr_ran_wlen_planck.h"
+#include "core/htrdr_spectral.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <math.h> /* nextafter */
+
+struct htrdr_ran_wlen_planck {
+ struct darray_double pdf;
+ struct darray_double cdf;
+ double range[2]; /* Boundaries of the spectral integration interval */
+ double band_len; /* Length in nanometers of a band */
+ double ref_temperature; /* In Kelvin */
+ size_t nbands; /* # bands */
+
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+setup_wlen_planck_ran_cdf
+ (struct htrdr_ran_wlen_planck* planck,
+ const char* func_name)
+{
+ double* pdf = NULL;
+ double* cdf = NULL;
+ double sum = 0;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(planck && func_name);
+ ASSERT(planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE);
+
+ res = darray_double_resize(&planck->cdf, planck->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(planck->htrdr,
+ "%s: Error allocating the CDF of the spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+ res = darray_double_resize(&planck->pdf, planck->nbands);
+ if(res != RES_OK) {
+ htrdr_log_err(planck->htrdr,
+ "%s: Error allocating the pdf of the spectral bands -- %s.\n",
+ func_name, res_to_cstr(res));
+ goto error;
+ }
+
+ cdf = darray_double_data_get(&planck->cdf);
+ pdf = darray_double_data_get(&planck->pdf); /* Now save pdf to correct weight */
+
+ htrdr_log(planck->htrdr,
+ "Number of bands used to speed up Planck distribution: %lu\n",
+ (unsigned long)planck->nbands);
+
+ /* Compute the *unnormalized* probability to sample a small band */
+ FOR_EACH(i, 0, planck->nbands) {
+ double lambda_lo = planck->range[0] + (double)i * planck->band_len;
+ double lambda_hi = MMIN(lambda_lo + planck->band_len, planck->range[1]);
+ ASSERT(lambda_lo<= lambda_hi);
+ ASSERT(lambda_lo > planck->range[0]
+ || eq_eps(lambda_lo, planck->range[0], 1.e-6));
+ ASSERT(lambda_lo < planck->range[1]
+ || eq_eps(lambda_lo, planck->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] = htrdr_blackbody_fraction
+ (lambda_lo, lambda_hi, planck->ref_temperature);
+
+ /* Update the norm */
+ sum += pdf[i];
+ }
+
+ /* Compute the cumulative of the previously computed probabilities */
+ FOR_EACH(i, 0, planck->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[planck->nbands-1], 1, 1.e-6));
+ cdf[planck->nbands - 1] = 1.0; /* Handle numerical issue */
+
+exit:
+ return res;
+error:
+ darray_double_clear(&planck->cdf);
+ darray_double_clear(&planck->pdf);
+ goto exit;
+}
+
+static double
+wlen_ran_sample_continue
+ (const struct htrdr_ran_wlen_planck* planck,
+ const double r,
+ const double range[2], /* In nanometer */
+ const char* func_name,
+ double* pdf)
+{
+ /* Numerical parameters of the dichotomy search */
+ const size_t MAX_ITER = 100;
+ const double EPSILON_LAMBDA_M = 1e-15;
+ const double EPSILON_BF = 1e-6;
+
+ /* Local variables */
+ double bf = 0;
+ double bf_prev = 0;
+ double bf_min_max = 0;
+ double lambda_m = 0;
+ double lambda_m_prev = 0;
+ double lambda_m_min = 0;
+ double lambda_m_max = 0;
+ double range_m[2] = {0, 0};
+ size_t i;
+
+ /* Check precondition */
+ ASSERT(planck && func_name);
+ ASSERT(range && range[0] < range[1]);
+ ASSERT(0 <= r && r < 1);
+
+ /* Convert the wavelength range in meters */
+ range_m[0] = range[0] * 1.0e-9;
+ range_m[1] = range[1] * 1.0e-9;
+
+ /* Setup the dichotomy search */
+ lambda_m_min = range_m[0];
+ lambda_m_max = range_m[1];
+ bf_min_max = htrdr_blackbody_fraction
+ (range_m[0], range_m[1], planck->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 = htrdr_blackbody_fraction
+ (range_m[0], lambda_m, planck->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(planck->htrdr,
+ "%s: could not sample a wavelength in the range [%g, %g] nanometers "
+ "for the reference temperature %g Kelvin.\n",
+ func_name, SPLIT2(range), planck->ref_temperature);
+ }
+
+ if(pdf) {
+ const double Tref = planck->ref_temperature; /* K */
+
+ /* W/m²/sr/m */
+ const double B_lambda = htrdr_planck(lambda_m, lambda_m, Tref);
+ const double B_mean = htrdr_planck(range_m[0], range_m[1], Tref);
+
+ *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0]));
+ *pdf *= 1.e-9; /* Transform from m⁻¹ to nm⁻¹ */
+ }
+
+ return lambda_m * 1.e9; /* Convert in nanometers */
+}
+
+static double
+wlen_ran_sample_discrete
+ (const struct htrdr_ran_wlen_planck* planck,
+ const double r0,
+ const double r1,
+ const char* func_name,
+ double* pdf)
+{
+ const double* cdf = NULL;
+ const double* find = NULL;
+ double r0_next = nextafter(r0, DBL_MAX);
+ double band_range[2];
+ double lambda = 0;
+ double pdf_continue = 0;
+ double pdf_band = 0;
+ size_t cdf_length = 0;
+ size_t i;
+ ASSERT(planck && planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE);
+ ASSERT(0 <= r0 && r0 < 1);
+ ASSERT(0 <= r1 && r1 < 1);
+ (void)func_name;
+ (void)pdf_band;
+
+ cdf = darray_double_cdata_get(&planck->cdf);
+ cdf_length = darray_double_size_get(&planck->cdf);
+ ASSERT(cdf_length > 0);
+
+ /* Use r_next rather than r0 in order to find the first entry that is not less
+ * than *or equal* to r0 */
+ find = search_lower_bound(&r0_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ ASSERT(find);
+
+ i = (size_t)(find - cdf);
+ ASSERT(i < cdf_length && cdf[i] > r0 && (!i || cdf[i-1] <= r0));
+
+ band_range[0] = planck->range[0] + (double)i*planck->band_len;
+ band_range[1] = band_range[0] + planck->band_len;
+
+ /* Fetch the pdf of the sampled band */
+ pdf_band = darray_double_cdata_get(&planck->pdf)[i];
+
+ /* Uniformly sample a wavelength in the sampled band */
+ lambda = band_range[0] + (band_range[1] - band_range[0]) * r1;
+ pdf_continue = 1.0 / (band_range[1] - band_range[0]);
+
+ if(pdf) {
+ *pdf = pdf_band * pdf_continue;
+ }
+
+ return lambda;
+}
+
+static void
+release_wlen_planck_ran(ref_T* ref)
+{
+ struct htrdr_ran_wlen_planck* planck = NULL;
+ struct htrdr* htrdr = NULL;
+ ASSERT(ref);
+ planck = CONTAINER_OF(ref, struct htrdr_ran_wlen_planck, ref);
+ darray_double_release(&planck->cdf);
+ darray_double_release(&planck->pdf);
+ htrdr = planck->htrdr;
+ MEM_RM(htrdr_get_allocator(htrdr), planck);
+ htrdr_ref_put(planck->htrdr);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_ran_wlen_planck_create
+ (struct htrdr* htrdr,
+ /* range must be included in [200,1000] nm for shortwave or in [1000,100000]
+ * nanometers for longwave (thermal) */
+ const double range[2],
+ const size_t nbands, /* # bands used to discretized CDF */
+ const double ref_temperature,
+ struct htrdr_ran_wlen_planck** out_wlen_planck_ran)
+{
+ struct htrdr_ran_wlen_planck* planck = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && range && out_wlen_planck_ran && ref_temperature > 0);
+ ASSERT(ref_temperature > 0);
+ ASSERT(range[0] <= range[1]);
+
+ planck = MEM_CALLOC(htrdr->allocator, 1, sizeof(*planck));
+ if(!planck) {
+ res = RES_MEM_ERR;
+ htrdr_log_err(htrdr,
+ "%s: could not allocate Planck distribution -- %s.\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ ref_init(&planck->ref);
+ darray_double_init(htrdr->allocator, &planck->cdf);
+ darray_double_init(htrdr->allocator, &planck->pdf);
+ htrdr_ref_get(htrdr);
+ planck->htrdr = htrdr;
+
+ planck->range[0] = range[0];
+ planck->range[1] = range[1];
+ planck->ref_temperature = ref_temperature;
+ planck->nbands = nbands;
+
+ if(nbands == HTRDR_WLEN_RAN_PLANCK_CONTINUE) {
+ planck->band_len = 0;
+ } else {
+ planck->band_len = (range[1] - range[0]) / (double)planck->nbands;
+
+ res = setup_wlen_planck_ran_cdf(planck, FUNC_NAME);
+ if(res != RES_OK) goto error;
+ }
+
+ htrdr_log(htrdr, "Spectral interval defined on [%g, %g] nanometers.\n",
+ range[0], range[1]);
+
+exit:
+ *out_wlen_planck_ran = planck;
+ return res;
+error:
+ if(planck) htrdr_ran_wlen_planck_ref_put(planck);
+ goto exit;
+}
+
+void
+htrdr_ran_wlen_planck_ref_get(struct htrdr_ran_wlen_planck* planck)
+{
+ ASSERT(planck);
+ ref_get(&planck->ref);
+}
+
+void
+htrdr_ran_wlen_planck_ref_put(struct htrdr_ran_wlen_planck* planck)
+{
+ ASSERT(planck);
+ ref_put(&planck->ref, release_wlen_planck_ran);
+}
+
+double
+htrdr_ran_wlen_planck_sample
+ (const struct htrdr_ran_wlen_planck* planck,
+ const double r0,
+ const double r1,
+ double* pdf)
+{
+ ASSERT(planck);
+ if(planck->nbands != HTRDR_WLEN_RAN_PLANCK_CONTINUE) { /* Discrete */
+ return wlen_ran_sample_discrete(planck, r0, r1, FUNC_NAME, pdf);
+ } else if(eq_eps(planck->range[0], planck->range[1], 1.e-6)) {
+ if(pdf) *pdf = 1;
+ return planck->range[0];
+ } else { /* Continue */
+ return wlen_ran_sample_continue
+ (planck, r0, planck->range, FUNC_NAME, pdf);
+ }
+}
+
diff --git a/src/core/htrdr_ran_wlen_planck.h b/src/core/htrdr_ran_wlen_planck.h
@@ -0,0 +1,61 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019 Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_RAN_WLEN_PLANCK_H
+#define HTRDR_RAN_WLEN_PLANCK_H
+
+#include "core/htrdr.h"
+#include <rsys/rsys.h>
+
+#define HTRDR_WLEN_RAN_PLANCK_CONTINUE 0
+
+/* Forward declarations */
+struct htrdr;
+struct htrdr_ran_wlen_planck;
+
+BEGIN_DECLS
+
+HTRDR_CORE_API res_T
+htrdr_ran_wlen_planck_create
+ (struct htrdr* htrdr,
+ const double range[2],
+ /* # bands used to discretisze the spectral domain. HTRDR_WLEN_RAN_CONTINUE
+ * <=> no discretisation */
+ const size_t nbands, /* Hint on #bands used to discretised th CDF */
+ const double ref_temperature, /* Reference temperature */
+ struct htrdr_ran_wlen_planck** planck);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_planck_ref_get
+ (struct htrdr_ran_wlen_planck* planck);
+
+HTRDR_CORE_API void
+htrdr_ran_wlen_planck_ref_put
+ (struct htrdr_ran_wlen_planck* planck);
+
+/* Return a wavelength in nanometer */
+HTRDR_CORE_API double
+htrdr_ran_wlen_planck_sample
+ (const struct htrdr_ran_wlen_planck* planck,
+ const double r0, /* Canonical number in [0, 1[ */
+ const double r1, /* Canonical number in [0, 1[ */
+ double* pdf); /* In nm^-1. May be NULL */
+
+END_DECLS
+
+#endif /* HTRDR_RAN_WLEN_PLANCK_H */
+
diff --git a/src/core/htrdr_spectral.h b/src/core/htrdr_spectral.h
@@ -36,12 +36,12 @@ enum htrdr_spectral_type {
struct htrdr;
static FINLINE double /* In nanometer */
-htrdr_wavenumber_to_wavelength(const double nu/*In cm^-1*/)
+htrdr_wavenumber_to_wavelength(const double nu/*In cm⁻¹*/)
{
return 1.e7 / nu;
}
-static FINLINE double /* In cm^-1 */
+static FINLINE double /* In cm⁻¹ */
wavelength_to_wavenumber(const double lambda/*In nanometer*/)
{
return htrdr_wavenumber_to_wavelength(lambda);
@@ -92,7 +92,7 @@ htrdr_blackbody_fraction
return w1 - w0;
}
-/* Return the Planck value in W/m^2/sr/m at a given wavelength */
+/* Return the Planck value in W/m²/sr/m at a given wavelength */
static INLINE double
htrdr_planck_monochromatic
(const double lambda, /* In meters */
@@ -103,13 +103,13 @@ htrdr_planck_monochromatic
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 */
+ const double B = ((2.0 * h * c*c) / lambda5) /* W/m²/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
+/* Return the average Planck value in W/m²/sr/m over the
* [lambda_min, lambda_max] interval. */
static INLINE double
htrdr_planck_interval
@@ -119,14 +119,14 @@ htrdr_planck_interval
{
const double T2 = temperature*temperature;
const double T4 = T2*T2;
- const double BOLTZMANN_CONSTANT = 5.6696e-8; /* W/m^2/K^4 */
+ const double BOLTZMANN_CONSTANT = 5.6696e-8; /* W/m²/K⁴ */
ASSERT(lambda_min < lambda_max && temperature > 0);
return htrdr_blackbody_fraction(lambda_min, lambda_max, temperature)
- * BOLTZMANN_CONSTANT * T4 / (PI * (lambda_max-lambda_min)); /* In W/m^2/sr/m */
+ * BOLTZMANN_CONSTANT * T4 / (PI * (lambda_max-lambda_min)); /* In W/m²/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 */
+ * interval is null or not, respectively. The returned value is in W/m²/sr/m */
static FINLINE double
htrdr_planck
(const double lambda_min, /* In meters */
@@ -148,7 +148,7 @@ htrdr_brightness_temperature
(struct htrdr* htrdr,
const double lambda_min, /* In meters */
const double lambda_max, /* In meters */
- /* Averaged over [lambda_min, lambda_max]. In W/m^2/sr/m */
+ /* Averaged over [lambda_min, lambda_max]. In W/m²/sr/m */
const double radiance,
double* temperature);
@@ -157,7 +157,7 @@ htrdr_radiance_temperature
(struct htrdr* htrdr,
const double lambda_min, /* In meters */
const double lambda_max, /* In meters */
- const double radiance); /* In W/m^2/sr */
+ const double radiance); /* In W/m²/sr */
END_DECLS
diff --git a/src/planeto/htrdr_planeto.c b/src/planeto/htrdr_planeto.c
@@ -0,0 +1,637 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, 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 /* fdopen, nextafter, rint */
+
+#include "core/htrdr.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_ran_wlen_discrete.h"
+#include "core/htrdr_ran_wlen_planck.h"
+#include "core/htrdr_log.h"
+
+#include "planeto/htrdr_planeto.h"
+#include "planeto/htrdr_planeto_args.h"
+#include "planeto/htrdr_planeto_c.h"
+#include "planeto/htrdr_planeto_source.h"
+
+#include <rad-net/rnatm.h>
+#include <rad-net/rngrd.h>
+
+#include <star/scam.h>
+
+#include <rsys/cstr.h>
+#include <rsys/double3.h>
+#include <rsys/mem_allocator.h>
+
+#include <fcntl.h> /* open */
+#include <math.h> /* nextafter, rint */
+#include <unistd.h> /* close */
+#include <sys/stat.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+/* Calculate the minimum length of the atmospheric spectral bands for the
+ * spectral domain considered */
+static double
+compute_min_band_len(const struct htrdr_planeto* cmd)
+{
+ const double* range = NULL; /* In nm */
+ double len = DBL_MAX;
+ size_t ibands[2];
+ size_t i;
+ ASSERT(cmd);
+
+ range = cmd->spectral_domain.wlen_range;
+
+ /* The spectral range is degenerate to a wavelength */
+ if(eq_eps(range[0], range[1], 1.e-6)) {
+ return 0;
+ }
+
+ RNATM(find_bands(cmd->atmosphere, cmd->spectral_domain.wlen_range, ibands));
+
+ /* At least one band must be overlaped by the spectral domain */
+ ASSERT(ibands[0]<=ibands[1]);
+ FOR_EACH(i, ibands[0], ibands[1]+1) {
+ struct rnatm_band_desc band;
+ double band_range[2];
+ RNATM(band_get_desc(cmd->atmosphere, i, &band));
+
+ /* Make the upper bound inclusive */
+ band_range[0] = band.lower;
+ band_range[1] = nextafter(band.upper, 0);
+
+ /* Clamp the band range to the spectral domain */
+ band_range[0] = MMAX(band_range[0], range[0]);
+ band_range[1] = MMIN(band_range[1], range[1]);
+ len = MMIN(band_range[1] - band_range[0], len);
+ }
+ return len;
+}
+
+/* Calculate the number of fixed size spectral intervals to use for the
+ * cumulative */
+static size_t
+compute_nintervals_for_spectral_cdf(const struct htrdr_planeto* cmd)
+{
+ double range_size;
+ double interval_len;
+ size_t nintervals;
+ ASSERT(cmd);
+
+ range_size =
+ cmd->spectral_domain.wlen_range[1]
+ - cmd->spectral_domain.wlen_range[0];
+
+ /* Initially assume ~one interval per nanometer */
+ nintervals = (size_t)rint(range_size);
+
+ /* Calculate the minimum length of the atmospheric spectral bands fixed to
+ * the spectral integration domain. We ensure that an interval of the
+ * spectral cdf cannot be greater than this length */
+ interval_len = compute_min_band_len(cmd);
+ if(interval_len < (range_size / (double)nintervals)) {
+ nintervals = (size_t)ceil(range_size / interval_len);
+ }
+
+ return nintervals;
+}
+
+static res_T
+setup_octree_storage
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args,
+ struct rnatm_create_args* rnatm_args)
+{
+ struct stat file_stat;
+ int fd = -1;
+ int err = 0;
+ res_T res = RES_OK;
+ ASSERT(cmd && args && rnatm_args);
+
+ rnatm_args->octrees_storage = NULL;
+ rnatm_args->load_octrees_from_storage = 0;
+
+ if(!args->octrees_storage) goto exit;
+
+ fd = open(args->octrees_storage, O_CREAT|O_RDWR, S_IRUSR|S_IWUSR);
+ if(fd < 0) { res = RES_IO_ERR; goto error; }
+
+ rnatm_args->octrees_storage = fdopen(fd, "w+");
+ if(!rnatm_args->octrees_storage) { res = RES_IO_ERR; goto error; }
+
+ /* From now on, manage the opened file from its pointer and not from its
+ * descriptor */
+ fd = -1;
+
+ err = stat(args->octrees_storage, &file_stat);
+ if(err < 0) { res = RES_IO_ERR; goto error; }
+
+ if(file_stat.st_size != 0) {
+ /* The file is not empty and therefore must contain valid octrees */
+ rnatm_args->load_octrees_from_storage = 1;
+ }
+
+exit:
+ cmd->octrees_storage = rnatm_args->octrees_storage;
+ return res;
+error:
+ htrdr_log_err(cmd->htrdr, "error opening the octree storage `%s' -- %s\n",
+ args->octrees_storage, res_to_cstr(res));
+
+ if(fd >= 0) CHK(close(fd) == 0);
+ if(rnatm_args->octrees_storage) CHK(fclose(rnatm_args->octrees_storage) == 0);
+ rnatm_args->octrees_storage = NULL;
+ rnatm_args->load_octrees_from_storage = 1;
+ goto exit;
+}
+
+static res_T
+setup_atmosphere
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args)
+{
+ struct rnatm_create_args rnatm_args = RNATM_CREATE_ARGS_DEFAULT;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ rnatm_args.gas = args->gas;
+ rnatm_args.aerosols = args->aerosols;
+ rnatm_args.naerosols = args->naerosols;
+ rnatm_args.name = "atmosphere";
+ rnatm_args.spectral_range[0] = args->spectral_domain.wlen_range[0];
+ rnatm_args.spectral_range[1] = args->spectral_domain.wlen_range[1];
+ rnatm_args.optical_thickness = args->optical_thickness;
+ rnatm_args.grid_definition_hint = args->octree_definition_hint;
+ rnatm_args.precompute_normals = args->precompute_normals;
+ rnatm_args.logger = htrdr_get_logger(cmd->htrdr);
+ rnatm_args.allocator = htrdr_get_allocator(cmd->htrdr);
+ rnatm_args.nthreads = args->nthreads;
+ rnatm_args.verbose = args->verbose;
+
+ res = setup_octree_storage(cmd, args, &rnatm_args);
+ if(res != RES_OK) goto error;
+
+ res = rnatm_create(&rnatm_args, &cmd->atmosphere);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ if(cmd->atmosphere) {
+ RNATM(ref_put(cmd->atmosphere));
+ cmd->atmosphere = NULL;
+ }
+ goto exit;
+}
+
+static res_T
+setup_ground
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args)
+{
+ struct rngrd_create_args rngrd_args = RNGRD_CREATE_ARGS_DEFAULT;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ if(cmd->output_type == HTRDR_PLANETO_ARGS_OUTPUT_OCTREES)
+ goto exit;
+
+ rngrd_args.smsh_filename = args->ground.smsh_filename;
+ rngrd_args.props_filename = args->ground.props_filename;
+ rngrd_args.mtllst_filename = args->ground.mtllst_filename;
+ rngrd_args.name = args->ground.name;
+ rngrd_args.logger = htrdr_get_logger(cmd->htrdr);
+ rngrd_args.allocator = htrdr_get_allocator(cmd->htrdr);
+ rngrd_args.verbose = args->verbose;
+
+ res = rngrd_create(&rngrd_args, &cmd->ground);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ if(cmd->ground) {
+ RNGRD(ref_put(cmd->ground));
+ cmd->ground = NULL;
+ }
+ goto exit;
+}
+
+static res_T
+setup_spectral_domain_sw
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args)
+{
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+ ASSERT(cmd->spectral_domain.type == HTRDR_SPECTRAL_SW);
+
+ /* Discrete distribution */
+ if(args->source.rnrl_filename) {
+ struct htrdr_planeto_source_spectrum spectrum;
+ struct htrdr_ran_wlen_discrete_create_args discrete_args;
+
+ res = htrdr_planeto_source_get_spectrum
+ (cmd->source, cmd->spectral_domain.wlen_range, &spectrum);
+ if(res != RES_OK) goto error;
+
+ discrete_args.get = htrdr_planeto_source_spectrum_at;
+ discrete_args.nwavelengths = spectrum.size;
+ discrete_args.context = &spectrum;
+ res = htrdr_ran_wlen_discrete_create
+ (cmd->htrdr, &discrete_args, &cmd->discrete);
+ if(res != RES_OK) goto error;
+
+ /* Planck distribution */
+ } else {
+ const size_t nintervals = compute_nintervals_for_spectral_cdf(cmd);
+
+ /* Use the source temperature as the reference temperature of the Planck
+ * distribution */
+ res = htrdr_ran_wlen_planck_create(cmd->htrdr,
+ cmd->spectral_domain.wlen_range, nintervals, args->source.temperature,
+ &cmd->planck);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
+setup_spectral_domain
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args)
+{
+ double ground_T_range[2];
+ size_t nintervals;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ cmd->spectral_domain = args->spectral_domain;
+
+ /* Configure the spectral distribution */
+ switch(cmd->spectral_domain.type) {
+
+ case HTRDR_SPECTRAL_LW:
+ res = rngrd_get_temperature_range(cmd->ground, ground_T_range);
+ if(res != RES_OK) goto error;
+
+ /* Use as the reference temperature of the Planck distribution the
+ * maximum scene temperature which, in fact, should be the maximum ground
+ * temperature */
+ nintervals = compute_nintervals_for_spectral_cdf(cmd);
+ res = htrdr_ran_wlen_planck_create(cmd->htrdr,
+ cmd->spectral_domain.wlen_range, nintervals, ground_T_range[1],
+ &cmd->planck);
+ if(res != RES_OK) goto error;
+ break;
+
+ case HTRDR_SPECTRAL_SW:
+ res = setup_spectral_domain_sw(cmd, args);
+ if(res != RES_OK) goto error;
+ break;
+
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ /* CIE XYZ distribution */
+ nintervals = compute_nintervals_for_spectral_cdf(cmd);
+ res = htrdr_ran_wlen_cie_xyz_create(cmd->htrdr,
+ cmd->spectral_domain.wlen_range, nintervals, &cmd->cie);
+ if(res != RES_OK) goto error;
+ break;
+
+ default: FATAL("Unreachable code\n"); break;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_output
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args)
+{
+ const char* output_name = NULL;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ /* No output stream on non master processes */
+ if(htrdr_get_mpi_rank(cmd->htrdr) != 0) {
+ cmd->output = NULL;
+ output_name = "<null>";
+
+ /* Write results on stdout */
+ } else if(!args->output) {
+ cmd->output = stdout;
+ output_name = "<stdout>";
+
+ /* Open the output stream */
+ } else {
+ res = htrdr_open_output_stream(cmd->htrdr, args->output, 0/*read*/,
+ args->force_output_overwrite, &cmd->output);
+ if(res != RES_OK) goto error;
+ output_name = args->output;
+ }
+
+ res = str_set(&cmd->output_name, output_name);
+ if(res != RES_OK) {
+ htrdr_log_err(cmd->htrdr, "error storing output stream name `%s' -- %s\n",
+ output_name, res_to_cstr(res));
+ goto error;
+ }
+
+ cmd->output_type = args->output_type;
+
+exit:
+ return res;
+error:
+ str_clear(&cmd->output_name);
+ if(cmd->output && cmd->output != stdout) {
+ CHK(fclose(cmd->output) == 0);
+ cmd->output = NULL;
+ }
+ goto exit;
+}
+
+static INLINE res_T
+setup_source
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args)
+{
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ if(cmd->output_type == HTRDR_PLANETO_ARGS_OUTPUT_OCTREES)
+ goto exit;
+
+ res = htrdr_planeto_source_create(cmd->htrdr, &args->source, &cmd->source);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_camera
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args)
+{
+ struct scam_perspective_args cam_args = SCAM_PERSPECTIVE_ARGS_DEFAULT;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ if(cmd->output_type != HTRDR_PLANETO_ARGS_OUTPUT_IMAGE)
+ goto exit;
+
+ ASSERT(htrdr_args_camera_perspective_check(&args->cam_persp) == RES_OK);
+ ASSERT(htrdr_args_image_check(&args->image) == RES_OK);
+
+ d3_set(cam_args.position, args->cam_persp.position);
+ d3_set(cam_args.target, args->cam_persp.target);
+ d3_set(cam_args.up, args->cam_persp.up);
+ cam_args.field_of_view = MDEG2RAD(args->cam_persp.fov_y);
+ cam_args.lens_radius = args->cam_persp.lens_radius;
+ cam_args.focal_distance = args->cam_persp.focal_dst;
+ cam_args.aspect_ratio =
+ (double)args->image.definition[0]
+ / (double)args->image.definition[1];
+
+ res = scam_create_perspective
+ (htrdr_get_logger(cmd->htrdr),
+ htrdr_get_allocator(cmd->htrdr),
+ htrdr_get_verbosity_level(cmd->htrdr),
+ &cam_args,
+ &cmd->camera);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_buffer
+ (struct htrdr_planeto* cmd,
+ const struct htrdr_planeto_args* args)
+{
+ struct htrdr_pixel_format pixfmt = HTRDR_PIXEL_FORMAT_NULL;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ if(cmd->output_type != HTRDR_PLANETO_ARGS_OUTPUT_IMAGE)
+ goto exit;
+
+ planeto_get_pixel_format(cmd, &pixfmt);
+
+ /* Setup buffer layout */
+ cmd->buf_layout.width = args->image.definition[0];
+ cmd->buf_layout.height = args->image.definition[1];
+ cmd->buf_layout.pitch = args->image.definition[0] * pixfmt.size;
+ cmd->buf_layout.elmt_size = pixfmt.size;
+ cmd->buf_layout.alignment = pixfmt.alignment;
+
+ /* Save the number of samples per pixel */
+ cmd->spp = args->image.spp;
+
+ /* Create the image buffer only on the master process; Image parts rendered
+ * by other processes are collected there */
+ if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
+
+ res = htrdr_buffer_create(cmd->htrdr, &cmd->buf_layout, &cmd->buf);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ if(cmd->buf) { htrdr_buffer_ref_put(cmd->buf); cmd->buf = NULL; }
+ goto exit;
+}
+
+static INLINE res_T
+write_vtk_octrees(const struct htrdr_planeto* cmd)
+{
+ size_t octrees_range[2];
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ /* Nothing to do on non master process */
+ if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
+
+ octrees_range[0] = 0;
+ octrees_range[1] = rnatm_get_spectral_items_count(cmd->atmosphere) - 1;
+
+ res = rnatm_write_vtk_octrees(cmd->atmosphere, octrees_range, cmd->output);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static void
+planeto_release(ref_T* ref)
+{
+ struct htrdr_planeto* cmd = CONTAINER_OF(ref, struct htrdr_planeto, ref);
+ struct htrdr* htrdr = NULL;
+ ASSERT(ref);
+
+ if(cmd->atmosphere) RNATM(ref_put(cmd->atmosphere));
+ if(cmd->ground) RNGRD(ref_put(cmd->ground));
+ if(cmd->source) htrdr_planeto_source_ref_put(cmd->source);
+ if(cmd->cie) htrdr_ran_wlen_cie_xyz_ref_put(cmd->cie);
+ if(cmd->discrete) htrdr_ran_wlen_discrete_ref_put(cmd->discrete);
+ if(cmd->planck) htrdr_ran_wlen_planck_ref_put(cmd->planck);
+ if(cmd->octrees_storage) CHK(fclose(cmd->octrees_storage) == 0);
+ if(cmd->output && cmd->output != stdout) CHK(fclose(cmd->output) == 0);
+ if(cmd->buf) htrdr_buffer_ref_put(cmd->buf);
+ if(cmd->camera) SCAM(ref_put(cmd->camera));
+ str_release(&cmd->output_name);
+
+ htrdr = cmd->htrdr;
+ MEM_RM(htrdr_get_allocator(htrdr), cmd);
+ htrdr_ref_put(htrdr);
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+htrdr_planeto_create
+ (struct htrdr* htrdr,
+ const struct htrdr_planeto_args* args,
+ struct htrdr_planeto** out_cmd)
+{
+ struct htrdr_planeto* cmd = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && out_cmd);
+
+ res = htrdr_planeto_args_check(args);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr, "Invalid htrdr_planeto arguments -- %s\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ cmd = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cmd));
+ if(!cmd) {
+ htrdr_log_err(htrdr, "Error allocating htrdr_planeto command\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&cmd->ref);
+ htrdr_ref_get(htrdr);
+ cmd->htrdr = htrdr;
+ str_init(htrdr_get_allocator(htrdr), &cmd->output_name);
+
+ res = setup_output(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_source(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_camera(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_ground(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_atmosphere(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_spectral_domain(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_buffer(cmd, args);
+ if(res != RES_OK) goto error;
+
+exit:
+ *out_cmd = cmd;
+ return res;
+error:
+ if(cmd) {
+ htrdr_planeto_ref_put(cmd);
+ cmd = NULL;
+ }
+ goto exit;
+}
+
+void
+htrdr_planeto_ref_get(struct htrdr_planeto* cmd)
+{
+ ASSERT(cmd);
+ ref_get(&cmd->ref);
+}
+
+void
+htrdr_planeto_ref_put(struct htrdr_planeto* cmd)
+{
+ ASSERT(cmd);
+ ref_put(&cmd->ref, planeto_release);
+}
+
+res_T
+htrdr_planeto_run(struct htrdr_planeto* cmd)
+{
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ switch(cmd->output_type) {
+ case HTRDR_PLANETO_ARGS_OUTPUT_IMAGE:
+ res = planeto_draw_map(cmd);
+ break;
+ case HTRDR_PLANETO_ARGS_OUTPUT_OCTREES:
+ res = write_vtk_octrees(cmd);
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+void
+planeto_get_pixel_format
+ (const struct htrdr_planeto* cmd,
+ struct htrdr_pixel_format* fmt)
+{
+ ASSERT(cmd && fmt && cmd->output_type == HTRDR_PLANETO_ARGS_OUTPUT_IMAGE);
+ (void)cmd;
+
+ switch(cmd->spectral_domain.type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ fmt->size = sizeof(struct planeto_pixel_xwave);
+ fmt->alignment = ALIGNOF(struct planeto_pixel_xwave);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ fmt->size = sizeof(struct planeto_pixel_image);
+ fmt->alignment = ALIGNOF(struct planeto_pixel_image);
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+}
diff --git a/src/planeto/htrdr_planeto.h b/src/planeto/htrdr_planeto.h
@@ -0,0 +1,55 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_PLANETO_H
+#define HTRDR_PLANETO_H
+
+#include "core/htrdr.h"
+#include <rsys/rsys.h>
+
+struct htrdr;
+struct htrdr_planeto;
+struct htrdr_planeto_args;
+
+BEGIN_DECLS
+
+HTRDR_API res_T
+htrdr_planeto_create
+ (struct htrdr* htrdr,
+ const struct htrdr_planeto_args* args,
+ struct htrdr_planeto** cmd);
+
+HTRDR_API void
+htrdr_planeto_ref_get
+ (struct htrdr_planeto* cmd);
+
+HTRDR_API void
+htrdr_planeto_ref_put
+ (struct htrdr_planeto* cmd);
+
+HTRDR_API res_T
+htrdr_planeto_run
+ (struct htrdr_planeto* cmd);
+
+HTRDR_API int
+htrdr_planeto_main
+ (int argc,
+ char** argv);
+
+END_DECLS
+
+#endif /* HTRDR_PLANETO_H */
diff --git a/src/planeto/htrdr_planeto_args.c b/src/planeto/htrdr_planeto_args.c
@@ -0,0 +1,781 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#define _POSIX_C_SOURCE 200112L /* strtok_r support */
+
+#include "planeto/htrdr_planeto_args.h"
+
+#include <rsys/cstr.h>
+#include <rsys/stretchy_array.h>
+#include <rsys/mem_allocator.h>
+
+#include <getopt.h>
+#include <string.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE res_T
+check_gas_args(const struct rnatm_gas_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Filenames cannot be NULL */
+ if(!args->smsh_filename
+ || !args->sck_filename
+ || !args->temperatures_filename)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+static INLINE res_T
+check_aerosol_args(const struct rnatm_aerosol_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Filenames cannot be NULL */
+ if(!args->smsh_filename
+ || !args->sars_filename
+ || !args->phase_fn_ids_filename
+ || !args->phase_fn_lst_filename)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+static INLINE res_T
+check_ground_args(const struct htrdr_planeto_ground_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Filenames cannot be NULL */
+ if(!args->smsh_filename
+ || !args->props_filename
+ || !args->mtllst_filename)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+static INLINE res_T
+check_spectral_args(const struct htrdr_planeto_spectral_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Invalid type */
+ switch(args->type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ /* Nothing to be done */
+ break;
+ default:
+ return RES_BAD_ARG;
+ }
+
+ /* Invalid spectral range */
+ if(args->wlen_range[0] < 0
+ || args->wlen_range[1] < 0
+ || args->wlen_range[0] > args->wlen_range[1])
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+static void
+print_help(const char* cmd)
+{
+ ASSERT(cmd);
+ printf(
+"Usage: %s [-dfhv] [-s spectral_domain] [-t threads]\n"
+" [-T optical_thickness] [-V octree_definition]\n"
+" [-O octrees_storage] [-o output] [-C camera]\n"
+" [-S source] [-a aerosol]... -G ground -g gas\n", cmd);
+ printf(
+"Simulate radiative transfer in heterogeneous 3D planetary atmosphere.\n"
+"See htrdr-planeto(1) man page for details\n\n");
+ printf(
+" -a aerosol define an aerosol\n");
+ printf(
+" -C camera configure a perspective camera\n");
+ printf(
+" -d write the atmospheric acceleration structures\n");
+ printf(
+" -f force overwrite the output file\n");
+ printf(
+" -G ground define the ground of the planet\n");
+ printf(
+" -g gas define the gas mixture\n");
+ printf(
+" -h display this help and exit\n");
+ printf(
+" -i image image to compute\n");
+ printf(
+" -N precompute tetrahedron normals\n");
+ printf(
+" -O octrees_storage\n"
+" file where atmospheric acceleration structures are\n"
+" stored/loaded\n");
+ printf(
+" -o output file where the result is written. If not defined,\n"
+" the result is written to standard output\n");
+ printf(
+" -S source define the source\n");
+ printf(
+" -s spectral_domain\n"
+" define the spectral domain of integration\n");
+ printf(
+" -T optical_thickness\n"
+" optical thickness criteria for octree building.\n"
+" Default is %g\n",
+ HTRDR_PLANETO_ARGS_DEFAULT.optical_thickness);
+ printf(
+" -t threads hint on the number of threads to use.\n"
+" Default assumes as many threads as CPU cores\n");
+ printf(
+" -V octree_definition\n"
+" advice on the definition of the atmospheric\n"
+" acceleration structures. Default is %u\n",
+ HTRDR_PLANETO_ARGS_DEFAULT.octree_definition_hint);
+ printf(
+" -v make the command verbose\n");
+ printf("\n");
+ htrdr_fprint_license(cmd, stdout);
+}
+
+static INLINE char*
+str_dup(const char* str)
+{
+ size_t len = 0;
+ char* dup = NULL;
+ ASSERT(str);
+ len = strlen(str) + 1/*NULL char*/;
+ dup = mem_alloc(len);
+ if(!dup) {
+ return NULL;
+ } else {
+ return memcpy(dup, str, len);
+ }
+}
+
+static res_T
+parse_aerosol_parameters(const char* str, void* ptr)
+{
+ enum { MESH, NAME, RADPROP, PHASEFN, PHASEIDS } iparam;
+ struct rnatm_aerosol_args* aerosol = NULL;
+ char buf[BUFSIZ];
+ struct htrdr_planeto_args* args = ptr;
+ char* key;
+ char* val;
+ char* tk_ctx;
+ res_T res = RES_OK;
+ ASSERT(args && str);
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr, "Could not duplicate the aerosol parameter `%s'\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &tk_ctx);
+ val = strtok_r(NULL, "", &tk_ctx);
+
+ if(!strcmp(key, "mesh")) iparam = MESH;
+ else if(!strcmp(key, "name")) iparam = NAME;
+ else if(!strcmp(key, "radprop")) iparam = RADPROP;
+ else if(!strcmp(key, "phasefn")) iparam = PHASEFN;
+ else if(!strcmp(key, "phaseids")) iparam = PHASEIDS;
+ else {
+ fprintf(stderr, "Invalid aerosol parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(!val) {
+ fprintf(stderr, "Invalid null value for aerosol parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ ASSERT(args->naerosols);
+ aerosol = args->aerosols + (args->naerosols - 1);
+
+ #define SET_STR(Dst) { \
+ if(Dst) mem_rm(Dst); \
+ if(!((Dst) = str_dup(val))) res = RES_MEM_ERR; \
+ } (void)0
+ switch(iparam) {
+ case MESH: SET_STR(aerosol->smsh_filename); break;
+ case NAME: SET_STR(aerosol->name); break;
+ case RADPROP: SET_STR(aerosol->sars_filename); break;
+ case PHASEFN: SET_STR(aerosol->phase_fn_lst_filename); break;
+ case PHASEIDS: SET_STR(aerosol->phase_fn_ids_filename); break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ #undef SET_STR
+ if(res != RES_OK) {
+ fprintf(stderr, "Unable to parse the aerosol parameter `%s' -- %s\n",
+ str, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_ground_parameters(const char* str, void* ptr)
+{
+ enum { BRDF, MESH, NAME, PROP } iparam;
+ char buf[BUFSIZ];
+ struct htrdr_planeto_args* args = ptr;
+ char* key;
+ char* val;
+ char* tk_ctx;
+ res_T res = RES_OK;
+ ASSERT(args && str);
+
+ if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
+ fprintf(stderr, "Could not duplicate the ground parameter `%s'\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &tk_ctx);
+ val = strtok_r(NULL, "", &tk_ctx);
+
+ if(!strcmp(key, "brdf")) iparam = BRDF;
+ else if(!strcmp(key, "mesh")) iparam = MESH;
+ else if(!strcmp(key, "name")) iparam = NAME;
+ else if(!strcmp(key, "prop")) iparam = PROP;
+ else {
+ fprintf(stderr, "Invalid ground parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(!val) {
+ fprintf(stderr, "Invalid null value for ground parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ #define SET_STR(Dst) { \
+ if(Dst) mem_rm(Dst); \
+ if(!((Dst) = str_dup(val))) res = RES_MEM_ERR; \
+ } (void)0
+ switch(iparam) {
+ case BRDF: SET_STR(args->ground.mtllst_filename); break;
+ case MESH: SET_STR(args->ground.smsh_filename); break;
+ case NAME: SET_STR(args->ground.name); break;
+ case PROP: SET_STR(args->ground.props_filename); break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ #undef SET_STR
+ if(res != RES_OK) {
+ fprintf(stderr, "Unable to parse the ground parameter `%s' -- %s\n",
+ str, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_gas_parameters(const char* str, void* ptr)
+{
+ enum { MESH, CK, TEMP } iparam;
+ char buf[BUFSIZ];
+ struct htrdr_planeto_args* args = ptr;
+ char* key;
+ char* val;
+ char* tk_ctx;
+ res_T res = RES_OK;
+ ASSERT(args && str);
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr, "Could not duplicate the gas parameter `%s'\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &tk_ctx);
+ val = strtok_r(NULL, "", &tk_ctx);
+
+ if(!strcmp(key, "mesh")) iparam = MESH;
+ else if(!strcmp(key, "ck")) iparam = CK;
+ else if(!strcmp(key, "temp")) iparam = TEMP;
+ else {
+ fprintf(stderr, "Invalid gas parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(!val) {
+ fprintf(stderr, "Invalid null value for gas parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ #define SET_STR(Dst) { \
+ if(Dst) mem_rm(Dst); \
+ if(!((Dst) = str_dup(val))) res = RES_MEM_ERR; \
+ } (void)0
+ switch(iparam) {
+ case MESH: SET_STR(args->gas.smsh_filename); break;
+ case CK: SET_STR(args->gas.sck_filename); break;
+ case TEMP: SET_STR(args->gas.temperatures_filename); break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ #undef SET_STR
+ if(res != RES_OK) {
+ fprintf(stderr, "Unable to parse the gas parameter `%s' -- %s\n",
+ str, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_source_parameters(const char* str, void* ptr)
+{
+ enum {LAT, LON, DST, RADIUS, TEMP, RAD} iparam;
+ char buf[BUFSIZ];
+ struct htrdr_planeto_args* args = ptr;
+ struct htrdr_planeto_source_args* src = NULL;
+ char* key;
+ char* val;
+ char* tk_ctx;
+ res_T res = RES_OK;
+ ASSERT(str && ptr);
+
+ src = &args->source;
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr, "Could not duplicate the source parameter `%s'\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &tk_ctx);
+ val = strtok_r(NULL, "", &tk_ctx);
+
+ if(!strcmp(key, "lat")) iparam = LAT;
+ else if(!strcmp(key, "lon")) iparam = LON;
+ else if(!strcmp(key, "dst")) iparam = DST;
+ else if(!strcmp(key, "rad")) iparam = RAD;
+ else if(!strcmp(key, "radius")) iparam = RADIUS;
+ else if(!strcmp(key, "temp")) iparam = TEMP;
+ else {
+ fprintf(stderr, "Invalid source parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(!val) {
+ fprintf(stderr, "Invalid null value for the source parameter`%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ switch(iparam) {
+ case LAT:
+ res = cstr_to_double(val, &src->latitude);
+ if(res == RES_OK && (src->latitude < -90 || src->latitude > 90)) {
+ res = RES_BAD_ARG;
+ }
+ break;
+ case LON:
+ res = cstr_to_double(val, &src->longitude);
+ if(res == RES_OK && (src->longitude < -180 || src->longitude > 180)) {
+ res = RES_BAD_ARG;
+ }
+ break;
+ case DST:
+ res = cstr_to_double(val, &src->distance);
+ if(res == RES_OK && src->distance < 0) res = RES_BAD_ARG;
+ break;
+ case RAD:
+ /* Use a per wavelength radiance rather than a constant temperature */
+ src->temperature = -1;
+ if(src->rnrl_filename) mem_rm(src->rnrl_filename);
+ src->rnrl_filename = str_dup(val);
+ if(!src->rnrl_filename) res = RES_MEM_ERR;
+ break;
+ case RADIUS:
+ res = cstr_to_double(val, &src->radius);
+ if(res == RES_OK && src->radius < 0) res = RES_BAD_ARG;
+ break;
+ case TEMP:
+ /* Use a constant temperature rather than a per wavelength radiance */
+ if(src->rnrl_filename) {
+ mem_rm(src->rnrl_filename);
+ src->rnrl_filename = NULL;
+ }
+ res = cstr_to_double(val, &src->temperature);
+ if(res == RES_OK && src->temperature < 0) res = RES_BAD_ARG;
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ if(res != RES_OK) {
+ fprintf(stderr, "Unable to parse the source parameter `%s' -- %s\n",
+ str, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
+parse_spectral_range(const char* str, double wlen_range[2])
+{
+ double range[2];
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(wlen_range && str);
+
+ res = cstr_to_list_double(str, ',', range, &len, 2);
+ if(res == RES_OK && len != 2) res = RES_BAD_ARG;
+ if(res == RES_OK && range[0] > range[1]) res = RES_BAD_ARG;
+ if(res == RES_OK && (range[0] < 0 || range[1] < 0)) res = RES_BAD_ARG;
+ if(res != RES_OK) goto error;
+
+ wlen_range[0] = range[0];
+ wlen_range[1] = range[1];
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_spectral_parameters(const char* str, void* ptr)
+{
+ enum {CIE_XYZ, LW, SW} iparam;
+ char buf[BUFSIZ];
+ struct htrdr_planeto_args* args = ptr;
+ struct htrdr_planeto_spectral_args* spectral = NULL;
+ char* key;
+ char* val;
+ char* tk_ctx;
+ res_T res = RES_OK;
+ ASSERT(str && ptr);
+
+ spectral = &args->spectral_domain;
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr, "Could not duplicate the spectral parameter `%s'\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &tk_ctx);
+ val = strtok_r(NULL, "", &tk_ctx);
+
+ if(!strcmp(key, "cie_xyz")) iparam = CIE_XYZ;
+ else if(!strcmp(key, "lw")) iparam = LW;
+ else if(!strcmp(key, "sw")) iparam = SW;
+ else {
+ fprintf(stderr, "Invalid spectral parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if((iparam == LW || iparam == SW) && !val) {
+ fprintf(stderr,
+ "Invalid null value for the spectral parameter `%s'\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ switch(iparam) {
+ case CIE_XYZ:
+ spectral->type = HTRDR_SPECTRAL_SW_CIE_XYZ;
+ spectral->wlen_range[0] = HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[0];
+ spectral->wlen_range[1] = HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT[1];
+ break;
+ case LW:
+ spectral->type = HTRDR_SPECTRAL_LW;
+ res = parse_spectral_range(val, spectral->wlen_range);
+ break;
+ case SW:
+ spectral->type = HTRDR_SPECTRAL_SW;
+ res = parse_spectral_range(val, spectral->wlen_range);
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ if(res != RES_OK) {
+ fprintf(stderr, "Unable to parse the spectral parameter `%s' -- %s\n",
+ str, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_planeto_args_init(struct htrdr_planeto_args* args, int argc, char** argv)
+{
+ int opt;
+ res_T res = RES_OK;
+ ASSERT(args && argc && argv);
+
+ *args = HTRDR_PLANETO_ARGS_DEFAULT;
+
+ while((opt = getopt(argc, argv, "a:C:dfG:g:hi:NO:o:S:s:T:t:V:v")) != -1) {
+ switch(opt) {
+ case 'a':
+ (void)sa_add(args->aerosols, 1);
+ args->aerosols[args->naerosols] = RNATM_AEROSOL_ARGS_NULL;
+ args->naerosols += 1;
+ res = cstr_parse_list(optarg, ':', parse_aerosol_parameters, args);
+ if(res == RES_OK) {
+ res = check_aerosol_args(args->aerosols+args->naerosols-1);
+ }
+ break;
+ case 'C':
+ res = htrdr_args_camera_perspective_parse(&args->cam_persp, optarg);
+ args->output_type = HTRDR_PLANETO_ARGS_OUTPUT_IMAGE;
+ break;
+ case 'd':
+ args->output_type = HTRDR_PLANETO_ARGS_OUTPUT_OCTREES;
+ break;
+ case 'f':
+ args->force_output_overwrite = 1;
+ break;
+ case 'G':
+ res = cstr_parse_list(optarg, ':', parse_ground_parameters, args);
+ if(res == RES_OK) {
+ res = check_ground_args(&args->ground);
+ }
+ break;
+ case 'g':
+ res = cstr_parse_list(optarg, ':', parse_gas_parameters, args);
+ if(res == RES_OK) {
+ res = check_gas_args(&args->gas);
+ }
+ break;
+ case 'h':
+ print_help(argv[0]);
+ htrdr_planeto_args_release(args);
+ args->quit = 1;
+ goto exit;
+ case 'i':
+ res = htrdr_args_image_parse(&args->image, optarg);
+ break;
+ case 'N': args->precompute_normals = 1; break;
+ case 'O': args->octrees_storage = optarg; break;
+ case 'o': args->output = optarg; break;
+ case 'S':
+ res = cstr_parse_list(optarg, ':', parse_source_parameters, args);
+ break;
+ case 's':
+ res = cstr_parse_list(optarg, ':', parse_spectral_parameters, args);
+ break;
+ case 'T':
+ res = cstr_to_double(optarg, &args->optical_thickness);
+ if(res != RES_OK && args->optical_thickness < 0) res = RES_BAD_ARG;
+ break;
+ case 't':
+ res = cstr_to_uint(optarg, &args->nthreads);
+ if(res != RES_OK && !args->nthreads) res = RES_BAD_ARG;
+ break;
+ case 'V':
+ res = cstr_to_uint(optarg, &args->octree_definition_hint);
+ if(res != RES_OK && !args->octree_definition_hint) res = RES_BAD_ARG;
+ break;
+ case 'v': args->verbose = 1; break;
+ default: res = RES_BAD_ARG; goto error;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+
+ res = check_gas_args(&args->gas);
+ if(res != RES_OK) {
+ fprintf(stderr, "Missing gas definition -- option '-g'\n");
+ goto error;
+ }
+
+ res = check_ground_args(&args->ground);
+ if(res != RES_OK) {
+ fprintf(stderr, "Missing ground definition -- option '-G'\n");
+ goto error;
+ }
+
+ if(args->output_type != HTRDR_PLANETO_ARGS_OUTPUT_OCTREES) {
+ res = check_ground_args(&args->ground);
+ if(res != RES_OK) {
+ fprintf(stderr, "Missing ground definition -- option '-G'\n");
+ goto error;
+ }
+
+ /* Check the source */
+ if(args->spectral_domain.type == HTRDR_SPECTRAL_SW
+ || args->spectral_domain.type == HTRDR_SPECTRAL_SW_CIE_XYZ) {
+ res = htrdr_planeto_source_args_check(&args->source);
+ if(res != RES_OK) {
+ fprintf(stderr, "Missing source definition -- option '-S'\n");
+ goto error;
+ }
+ }
+ }
+
+exit:
+ return res;
+error:
+ htrdr_planeto_args_release(args);
+ goto exit;
+}
+
+void
+htrdr_planeto_args_release(struct htrdr_planeto_args* args)
+{
+ size_t i;
+ ASSERT(args);
+
+ if(args->gas.smsh_filename) mem_rm(args->gas.smsh_filename);
+ if(args->gas.sck_filename) mem_rm(args->gas.sck_filename);
+ if(args->gas.temperatures_filename) mem_rm(args->gas.temperatures_filename);
+ if(args->ground.smsh_filename) mem_rm(args->ground.smsh_filename);
+ if(args->ground.props_filename) mem_rm(args->ground.props_filename);
+ if(args->ground.mtllst_filename) mem_rm(args->ground.mtllst_filename);
+ if(args->ground.name) mem_rm(args->ground.name);
+ if(args->source.rnrl_filename) mem_rm(args->source.rnrl_filename);
+
+ FOR_EACH(i, 0, args->naerosols) {
+ struct rnatm_aerosol_args* aerosol = args->aerosols + i;
+ if(aerosol->name) mem_rm(aerosol->name);
+ if(aerosol->smsh_filename) mem_rm(aerosol->smsh_filename);
+ if(aerosol->sars_filename) mem_rm(aerosol->sars_filename);
+ if(aerosol->phase_fn_ids_filename) mem_rm(aerosol->phase_fn_ids_filename);
+ if(aerosol->phase_fn_lst_filename) mem_rm(aerosol->phase_fn_lst_filename);
+ }
+ sa_release(args->aerosols);
+
+ *args = HTRDR_PLANETO_ARGS_DEFAULT;
+}
+
+res_T
+htrdr_planeto_args_check(const struct htrdr_planeto_args* args)
+{
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!args) return RES_BAD_ARG;
+
+ /* Check the gas */
+ res = check_gas_args(&args->gas);
+ if(res != RES_OK) return res;
+
+ /* Check the aerosols */
+ FOR_EACH(i, 0, args->naerosols) {
+ res = check_aerosol_args(args->aerosols+i);
+ if(res != RES_OK) return res;
+ }
+
+ /* Check the octree parameters */
+ if(args->octree_definition_hint == 0
+ || args->optical_thickness < 0)
+ return RES_BAD_ARG;
+
+ /* Check the spectral domain */
+ res = check_spectral_args(&args->spectral_domain);
+ if(res != RES_OK) return res;
+
+ if(args->output_type != HTRDR_PLANETO_ARGS_OUTPUT_OCTREES) {
+ /* Check the ground */
+ res = check_ground_args(&args->ground);
+ if(res != RES_OK) return res;
+
+ /* Check the source */
+ if(args->spectral_domain.type == HTRDR_SPECTRAL_SW
+ || args->spectral_domain.type == HTRDR_SPECTRAL_SW_CIE_XYZ) {
+ res = htrdr_planeto_source_args_check(&args->source);
+ if(res != RES_OK) return res;
+ }
+ }
+
+ if(args->output_type != HTRDR_PLANETO_ARGS_OUTPUT_IMAGE) {
+ res = htrdr_args_camera_perspective_check(&args->cam_persp);
+ if(res != RES_OK) return res;
+
+ res = htrdr_args_image_check(&args->image);
+ if(res != RES_OK) return res;
+ }
+
+ /* Check miscalleneous parameters */
+ if(args->nthreads == 0
+ || (unsigned)args->output_type >= HTRDR_PLANETO_ARGS_OUTPUT_TYPES_COUNT__)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+res_T
+htrdr_planeto_source_args_check(const struct htrdr_planeto_source_args* args)
+{
+ if(!args) return RES_BAD_ARG;
+
+ /* Invalid position */
+ if(args->latitude <-90
+ || args->latitude > 90
+ || args->longitude <-180
+ || args->longitude > 180
+ || args->distance < 0)
+ return RES_BAD_ARG;
+
+ /* Invalid radius */
+ if(args->radius < 0)
+ return RES_BAD_ARG;
+
+ /* Invalid radiance */
+ if((args->temperature < 0 && !args->rnrl_filename) /* Both are invalids */
+ || (args->temperature >=0 && args->rnrl_filename)) /* Both are valids */
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
diff --git a/src/planeto/htrdr_planeto_args.h.in b/src/planeto/htrdr_planeto_args.h.in
@@ -0,0 +1,144 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_PLANETO_ARGS_H
+#define HTRDR_PLANETO_ARGS_H
+
+#include "core/htrdr_args.h"
+
+#include <rad-net/rnatm.h>
+#include <rsys/rsys.h>
+
+#include <limits.h> /* UINT_MAX */
+
+enum htrdr_planeto_args_output_type {
+ HTRDR_PLANETO_ARGS_OUTPUT_IMAGE,
+ HTRDR_PLANETO_ARGS_OUTPUT_OCTREES,
+ HTRDR_PLANETO_ARGS_OUTPUT_TYPES_COUNT__
+};
+
+struct htrdr_planeto_spectral_args {
+ double wlen_range[2]; /* Spectral range in nm */
+ enum htrdr_spectral_type type;
+};
+#define HTRDR_PLANETO_SPECTRAL_ARGS_DEFAULT__ { \
+ HTRDR_RAN_WLEN_CIE_XYZ_RANGE_DEFAULT__, /* Spectral range */ \
+ HTRDR_SPECTRAL_SW_CIE_XYZ, /* Spectral type */ \
+}
+static const struct htrdr_planeto_spectral_args
+HTRDR_PLANETO_SPECTRAL_ARGS_DEFAULT = HTRDR_PLANETO_SPECTRAL_ARGS_DEFAULT__;
+
+struct htrdr_planeto_source_args {
+ /* Radiance of the source per wavelength. May be NULL */
+ char* rnrl_filename;
+
+ double longitude; /* In [-180, 180] degrees */
+ double latitude; /* In [-90, 90] degrees */
+ double distance; /* In km */
+ double radius; /* In km */
+ double temperature; /* In Kelvin */
+};
+#define HTRDR_PLANETO_SOURCE_ARGS_NULL__ {NULL,0,0,0,-1,-1}
+static const struct htrdr_planeto_source_args HTRDR_PLANETO_SOURCE_ARGS_NULL =
+ HTRDR_PLANETO_SOURCE_ARGS_NULL__;
+
+struct htrdr_planeto_ground_args {
+ char* smsh_filename; /* The Star-Mesh geometry */
+ char* props_filename; /* Per triangle physical properties */
+ char* mtllst_filename; /* List of used materials */
+ char* name;
+};
+#define HTRDR_PLANETO_GROUND_ARGS_NULL__ {NULL,NULL,NULL,NULL}
+static const struct htrdr_planeto_ground_args HTRDR_PLANETO_GROUND_ARGS_NULL =
+ HTRDR_PLANETO_GROUND_ARGS_NULL__;
+
+struct htrdr_planeto_args {
+ /* System data */
+ struct rnatm_gas_args gas;
+ struct rnatm_aerosol_args* aerosols;
+ size_t naerosols;
+ struct htrdr_planeto_ground_args ground;
+
+ /* Read/Write file where octrees are stored. May be NULL => octres are built
+ * at runtime and kept in memory */
+ char* octrees_storage;
+
+ unsigned octree_definition_hint; /* Hint on octree definition */
+ double optical_thickness; /* Threshold used during octree building */
+
+ char* output; /* File where the result is written */
+ struct htrdr_planeto_spectral_args spectral_domain; /* Integration spectral domain */
+ struct htrdr_planeto_source_args source;
+ struct htrdr_args_image image;
+
+ struct htrdr_args_camera_perspective cam_persp; /* Perspective camera */
+
+ /* Miscellaneous arguments */
+ unsigned nthreads; /* Hint on the nimber of threads to use */
+ enum htrdr_planeto_args_output_type output_type;
+ int precompute_normals; /* Pre-compute tetrahedron normals */
+ int force_output_overwrite; /* Replace output if it exists */
+ int verbose; /* Verbose level */
+ int quit; /* Stop the command */
+};
+#define HTRDR_PLANETO_ARGS_DEFAULT__ { \
+ RNATM_GAS_ARGS_NULL__, /* Gas */ \
+ NULL, /* List of aerosols */ \
+ 0, /* Number of aerosols */ \
+ HTRDR_PLANETO_GROUND_ARGS_NULL__, /* Ground */ \
+ \
+ NULL, /* File where to dump octrees */ \
+ \
+ @HTRDR_PLANETO_ARGS_DEFAULT_GRID_DEFINITION_HINT@, /* octree definition */ \
+ @HTRDR_PLANETO_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@, \
+ \
+ NULL, /* Ouput file */ \
+ HTRDR_PLANETO_SPECTRAL_ARGS_DEFAULT__, /* Spectral domain */ \
+ HTRDR_PLANETO_SOURCE_ARGS_NULL__, /* Source */ \
+ HTRDR_ARGS_IMAGE_DEFAULT__, /* Image */ \
+ \
+ HTRDR_ARGS_CAMERA_PERSPECTIVE_DEFAULT__, /* Perspective camera */ \
+ \
+ UINT_MAX, /* Number of threads */ \
+ HTRDR_PLANETO_ARGS_OUTPUT_IMAGE, \
+ 0, /* Force output overwrite */ \
+ 0, /* Precompute normals */ \
+ 0, /* Verbosity level */ \
+ 0 /* Stop the command */ \
+}
+static const struct htrdr_planeto_args HTRDR_PLANETO_ARGS_DEFAULT =
+ HTRDR_PLANETO_ARGS_DEFAULT__;
+
+extern LOCAL_SYM res_T
+htrdr_planeto_args_init
+ (struct htrdr_planeto_args* args,
+ int argc,
+ char** argv);
+
+extern LOCAL_SYM void
+htrdr_planeto_args_release
+ (struct htrdr_planeto_args* args);
+
+extern LOCAL_SYM res_T
+htrdr_planeto_args_check
+ (const struct htrdr_planeto_args* args);
+
+extern LOCAL_SYM res_T
+htrdr_planeto_source_args_check
+ (const struct htrdr_planeto_source_args* args);
+
+#endif /* HTRDR_PLANETO_ARGS_H */
diff --git a/src/planeto/htrdr_planeto_c.h b/src/planeto/htrdr_planeto_c.h
@@ -0,0 +1,121 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_PLANETO_C_H
+#define HTRDR_PLANETO_C_H
+
+#include "planeto/htrdr_planeto_args.h"
+
+#include "core/htrdr_accum.h"
+#include "core/htrdr_args.h"
+#include "core/htrdr_buffer.h"
+
+#include <rsys/ref_count.h>
+#include <rsys/str.h>
+
+/* Forward declarations */
+struct htrdr;
+struct htrdr_pixel_format;
+struct htrdr_ran_wlen_cie_xyz;
+struct htrdr_ran_wlen_planck;
+struct rnatm;
+struct rngrd;
+struct scam;
+
+struct planeto_pixel_xwave {
+ struct htrdr_accum radiance; /* In W/m²/sr */
+ struct htrdr_accum time; /* In µs */
+ struct htrdr_estimate radiance_temperature; /* In W/m²/sr */
+};
+#define PLANETO_PIXEL_XWAVE_NULL__ { \
+ HTRDR_ACCUM_NULL__, \
+ HTRDR_ACCUM_NULL__, \
+ HTRDR_ESTIMATE_NULL__ \
+}
+static const struct planeto_pixel_xwave PLANETO_PIXEL_XWAVE_NULL =
+ PLANETO_PIXEL_XWAVE_NULL__;
+
+struct planeto_pixel_image {
+ struct htrdr_estimate X; /* In W/m²/sr */
+ struct htrdr_estimate Y; /* In W/m²/sr */
+ struct htrdr_estimate Z; /* In W/m²/sr */
+ struct htrdr_accum time; /* In µs */
+};
+#define PLANETO_PIXEL_IMAGE_NULL__ { \
+ HTRDR_ESTIMATE_NULL__, \
+ HTRDR_ESTIMATE_NULL__, \
+ HTRDR_ESTIMATE_NULL__, \
+ HTRDR_ACCUM_NULL__ \
+}
+
+struct planeto_compute_radiance_args {
+ struct ssp_rng* rng;
+ size_t ithread; /* Index of the thread executing the function */
+
+ double path_org[3]; /* Origin of the path to trace */
+ double path_dir[3]; /* Initial direction of the path to trace */
+
+ double wlen; /* In nm */
+ size_t iband; /* Spectral band index */
+ size_t iquad; /* Quadrature point */
+};
+#define PLANETO_COMPUTE_RADIANCE_ARGS_NULL__ {NULL, 0, {0,0,0}, {0,0,0}, 0, 0, 0}
+static const struct planeto_compute_radiance_args
+PLANETO_COMPUTE_RADIANCE_ARGS_NULL = PLANETO_COMPUTE_RADIANCE_ARGS_NULL__;
+
+struct htrdr_planeto {
+ struct rnatm* atmosphere;
+ struct rngrd* ground;
+ struct htrdr_planeto_source* source;
+
+ struct htrdr_planeto_spectral_args spectral_domain;
+ struct htrdr_ran_wlen_cie_xyz* cie; /* HTRDR_SPECTRAL_SW_CIE_XYZ */
+ struct htrdr_ran_wlen_planck* planck; /* HTRDR_SPECTRAL_<LW|SW> */
+ struct htrdr_ran_wlen_discrete* discrete; /* HTRDR_SPECTRAL_SW */
+
+ FILE* octrees_storage;
+
+ FILE* output;
+ struct str output_name;
+ enum htrdr_planeto_args_output_type output_type;
+
+ struct scam* camera;
+
+ struct htrdr_buffer_layout buf_layout;
+ struct htrdr_buffer* buf; /* NULL on non master processes */
+ size_t spp; /* Samples per pixel */
+
+ ref_T ref;
+ struct htrdr* htrdr;
+};
+
+extern LOCAL_SYM res_T
+planeto_draw_map
+ (struct htrdr_planeto* cmd);
+
+extern LOCAL_SYM void
+planeto_get_pixel_format
+ (const struct htrdr_planeto* cmd,
+ struct htrdr_pixel_format* fmt);
+
+/* Return the radiance in W/m²/sr/m */
+extern LOCAL_SYM double
+planeto_compute_radiance
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args);
+
+#endif /* HTRDR_PLANETO_C_H */
diff --git a/src/planeto/htrdr_planeto_compute_radiance.c b/src/planeto/htrdr_planeto_compute_radiance.c
@@ -0,0 +1,665 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, 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 "planeto/htrdr_planeto_c.h"
+#include "planeto/htrdr_planeto_source.h"
+
+#include <rad-net/rnatm.h>
+#include <rad-net/rngrd.h>
+
+#include <star/s3d.h>
+#include <star/ssf.h>
+#include <star/ssp.h>
+#include <star/suvm.h>
+#include <star/svx.h>
+
+#include <rsys/double2.h>
+#include <rsys/double3.h>
+
+/* Syntactic sugar */
+#define DISTANCE_NONE(Dst) ((Dst) >= FLT_MAX)
+#define SURFACE_EVENT(Event) (!S3D_HIT_NONE(&(Event)->hit))
+
+struct event {
+ /* Set to S3D_HIT_NULL if the event occurs in volume.*/
+ struct s3d_hit hit;
+
+ /* The surface normal is defined only if event is on the surface. It is
+ * normalized and looks towards the incoming direction */
+ double normal[3];
+
+ /* Cells in which the event position is located. It makes sense only for an
+ * event in volume */
+ struct rnatm_cell_pos cells[RNATM_MAX_COMPONENTS_COUNT];
+
+ double distance; /* Distance from ray origin to scattering position */
+};
+#define EVENT_NULL__ { \
+ S3D_HIT_NULL__, {0,0,0}, {RNATM_CELL_POS_NULL__}, DBL_MAX \
+}
+static const struct event EVENT_NULL = EVENT_NULL__;
+
+/* Arguments of the filtering function used to sample a position */
+struct sample_distance_context {
+ struct ssp_rng* rng;
+ struct rnatm* atmosphere;
+ size_t iband;
+ size_t iquad;
+ double wavelength; /* In nm */
+ enum rnatm_radcoef radcoef;
+ double Ts; /* Sample optical thickness */
+
+ /* Output data */
+ struct rnatm_cell_pos* cells;
+ double distance;
+};
+#define SAMPLE_DISTANCE_CONTEXT_NULL__ { \
+ NULL, NULL, 0, 0, 0, RNATM_RADCOEFS_COUNT__, 0, NULL, DBL_MAX \
+}
+static const struct sample_distance_context SAMPLE_DISTANCE_CONTEXT_NULL =
+ SAMPLE_DISTANCE_CONTEXT_NULL__;
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE res_T
+check_planeto_compute_radiance_args
+ (const struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args)
+{
+ struct rnatm_band_desc band = RNATM_BAND_DESC_NULL;
+ res_T res = RES_OK;
+
+ if(!args || !args->rng)
+ return RES_BAD_ARG;
+
+ /* Invalid thread index */
+ if(args->ithread >= htrdr_get_threads_count(cmd->htrdr))
+ return RES_BAD_ARG;
+
+ /* Invalid input direction */
+ if(!d3_is_normalized(args->path_dir))
+ return RES_BAD_ARG;
+
+ /* Invalid wavelength */
+ if(args->wlen < cmd->spectral_domain.wlen_range[0]
+ || args->wlen > cmd->spectral_domain.wlen_range[1])
+ return RES_BAD_ARG;
+
+ res = rnatm_band_get_desc(cmd->atmosphere, args->iband, &band);
+ if(res != RES_OK) return res;
+
+ /* Inconsistent spectral dimension */
+ if(args->wlen < band.lower
+ || args->wlen >= band.upper /* Exclusive */
+ || args->iquad >= band.quad_pts_count)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+static int
+sample_position_hit_filter
+ (const struct svx_hit* hit,
+ const double org[3],
+ const double dir[3],
+ const double range[2],
+ void* context)
+{
+ struct rnatm_get_radcoef_args get_k_args = RNATM_GET_RADCOEF_ARGS_NULL;
+ struct sample_distance_context* ctx = context;
+ double k_min = 0;
+ double k_max = 0;
+ double dst_travelled = 0;
+ int pursue_traversal = 1;
+ ASSERT(hit && org && range && context);
+ (void)range;
+
+ dst_travelled = hit->distance[0];
+
+ /* Get the k_min, k_max of the voxel */
+ k_min = rnatm_get_k_svx_voxel
+ (ctx->atmosphere, &hit->voxel, ctx->radcoef, RNATM_SVX_OP_MIN);
+ k_max = rnatm_get_k_svx_voxel
+ (ctx->atmosphere, &hit->voxel, ctx->radcoef, RNATM_SVX_OP_MAX);
+
+ /* Configure the common input arguments to retrieve the radiative coefficient
+ * of a given position */
+ get_k_args.cells = ctx->cells;
+ get_k_args.iband = ctx->iband;
+ get_k_args.iquad = ctx->iquad;
+ get_k_args.radcoef = ctx->radcoef;
+ get_k_args.k_min = k_min;
+ get_k_args.k_max = k_max;
+
+ for(;;) {
+ /* Compute the optical thickness of the voxel */
+ const double vox_dst = hit->distance[1] - dst_travelled;
+ const double T = vox_dst * k_max;
+
+ /* A collision occurs behind the voxel */
+ if(ctx->Ts > T) {
+ ctx->Ts -= T;
+ pursue_traversal = 1;
+ break;
+
+ /* A collision occurs in the voxel */
+ } else {
+ const double vox_dst_collision = ctx->Ts / k_max;
+ double pos[3];
+ double k, r;
+
+ /* Calculate the distance travelled to the collision to be queried */
+ dst_travelled += vox_dst_collision;
+
+ /* Retrieve the radiative coefficient at the collision position */
+ pos[0] = org[0] + dst_travelled * dir[0];
+ pos[1] = org[1] + dst_travelled * dir[1];
+ pos[2] = org[2] + dst_travelled * dir[2];
+ RNATM(fetch_cell_list(ctx->atmosphere, pos, get_k_args.cells, NULL));
+ RNATM(get_radcoef(ctx->atmosphere, &get_k_args, &k));
+
+ r = ssp_rng_canonical(ctx->rng);
+
+ /* Null collision */
+ if(r > k/k_max) {
+ ctx->Ts = ssp_ran_exp(ctx->rng, 1);
+
+ /* Real collision */
+ } else {
+ ctx->distance = dst_travelled;
+ pursue_traversal = 0;
+ break;
+ }
+ }
+ }
+
+ return pursue_traversal;
+}
+
+static double
+sample_distance
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ struct rnatm_cell_pos* cells,
+ const enum rnatm_radcoef radcoef,
+ const double pos[3],
+ const double dir[3],
+ const double range[2])
+{
+ struct rnatm_trace_ray_args rt = RNATM_TRACE_RAY_ARGS_NULL;
+ struct sample_distance_context ctx = SAMPLE_DISTANCE_CONTEXT_NULL;
+ struct svx_hit hit;
+ ASSERT(cmd && args && cells && pos && dir && d3_is_normalized(dir) && range);
+ ASSERT((unsigned)radcoef < RNATM_RADCOEFS_COUNT__);
+ ASSERT(range[0] < range[1]);
+
+ /* Sample an optical thickness */
+ ctx.Ts = ssp_ran_exp(args->rng, 1);
+
+ /* Setup the remaining arguments of the RT context */
+ ctx.rng = args->rng;
+ ctx.atmosphere = cmd->atmosphere;
+ ctx.iband = args->iband;
+ ctx.iquad = args->iquad;
+ ctx.wavelength = args->wlen;
+ ctx.radcoef = radcoef;
+ ctx.cells = cells;
+
+ /* Trace the ray into the atmosphere */
+ d3_set(rt.ray_org, pos);
+ d3_set(rt.ray_dir, dir);
+ rt.ray_range[0] = range[0];
+ rt.ray_range[1] = range[1];
+ rt.filter = sample_position_hit_filter;
+ rt.context = &ctx;
+ rt.iband = args->iband;
+ rt.iquad = args->iquad;
+ RNATM(trace_ray(cmd->atmosphere, &rt, &hit));
+
+ if(SVX_HIT_NONE(&hit)) { /* No collision found */
+ return INF;
+ } else { /* A (real) collision occured */
+ return ctx.distance;
+ }
+}
+
+static INLINE double
+transmissivity
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ const enum rnatm_radcoef radcoef,
+ const double pos[3],
+ const double dir[3],
+ const double range_max)
+{
+ struct rnatm_cell_pos cells[RNATM_MAX_COMPONENTS_COUNT];
+ double range[2];
+ double dst = 0;
+ ASSERT(range_max >= 0);
+
+ range[0] = 0;
+ range[1] = range_max;
+ dst = sample_distance(cmd, args, cells, radcoef, pos, dir, range);
+
+ if(DISTANCE_NONE(dst)) {
+ return 1.0; /* No collision in the medium */
+ } else {
+ return 0.0; /* A (real) collision occurs */
+ }
+}
+
+static double
+direct_contribution
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ const double pos[3],
+ const double dir[3],
+ const struct s3d_hit* hit_from)
+{
+ struct rngrd_trace_ray_args rt = RNGRD_TRACE_RAY_ARGS_DEFAULT;
+ struct s3d_hit hit;
+ double Tr;
+ double Ld;
+ double src_dst;
+ ASSERT(cmd && args && pos && dir);
+
+ /* Is the source hidden? */
+ d3_set(rt.ray_org, pos);
+ d3_set(rt.ray_dir, dir);
+ if(hit_from) rt.hit_from = *hit_from;
+ RNGRD(trace_ray(cmd->ground, &rt, &hit));
+ if(!S3D_HIT_NONE(&hit)) return 0;
+
+ /* Calculate the distance between the source and `pos' */
+ src_dst = htrdr_planeto_source_distance_to(cmd->source, pos);
+ ASSERT(src_dst >= 0);
+
+ Tr = transmissivity(cmd, args, RNATM_RADCOEF_Kext, pos, dir, src_dst);
+ Ld = htrdr_planeto_source_get_radiance(cmd->source, args->wlen);
+ return Ld * Tr;
+}
+
+static void
+find_event
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ const enum rnatm_radcoef radcoef,
+ const double pos[3],
+ const double dir[3],
+ const struct s3d_hit* hit_from,
+ struct event* evt)
+{
+ struct rngrd_trace_ray_args rt = RNGRD_TRACE_RAY_ARGS_DEFAULT;
+ struct s3d_hit hit;
+ double range[2];
+ double dst;
+ ASSERT(cmd && args && pos && dir && hit_from && evt);
+
+ *evt = EVENT_NULL;
+
+ /* Look for a surface intersection */
+ d3_set(rt.ray_org, pos);
+ d3_set(rt.ray_dir, dir);
+ d2(rt.ray_range, 0, INF);
+ rt.hit_from = *hit_from;
+ RNGRD(trace_ray(cmd->ground, &rt, &hit));
+
+ /* Look for an atmospheric collision */
+ range[0] = 0;
+ range[1] = hit.distance;
+ dst = sample_distance(cmd, args, evt->cells, radcoef, pos, dir, range);
+
+ /* Event occurs in volume */
+ if(!DISTANCE_NONE(dst)) {
+ evt->distance = dst;
+ evt->hit = S3D_HIT_NULL;
+
+ /* Event is on surface */
+ } else if(!S3D_HIT_NONE(&hit)) {
+ /* Normalize the normal and ensure that it points to `dir' */
+ d3_normalize(evt->normal, d3_set_f3(evt->normal, hit.normal));
+ if(d3_dot(evt->normal, dir) > 0) d3_minus(evt->normal, evt->normal);
+
+ evt->distance = hit.distance;
+ evt->hit = hit;
+
+ /* No event */
+ } else {
+ evt->distance = INF;
+ evt->hit = S3D_HIT_NULL;
+ }
+}
+
+static INLINE struct ssf_bsdf*
+create_bsdf
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ const struct s3d_hit* hit)
+{
+ struct rngrd_create_bsdf_args bsdf_args = RNGRD_CREATE_BSDF_ARGS_NULL;
+ struct ssf_bsdf* bsdf = NULL;
+ ASSERT(!S3D_HIT_NONE(hit));
+
+ /* Retrieve the BSDF at the intersected surface position */
+ bsdf_args.prim = hit->prim;
+ bsdf_args.barycentric_coords[0] = hit->uv[0];
+ bsdf_args.barycentric_coords[1] = hit->uv[1];
+ bsdf_args.barycentric_coords[2] = 1 - hit->uv[0] - hit->uv[1];
+ bsdf_args.wavelength = args->wlen;
+ bsdf_args.r = ssp_rng_canonical(args->rng);
+ RNGRD(create_bsdf(cmd->ground, &bsdf_args, &bsdf));
+
+ return bsdf;
+}
+
+static INLINE struct ssf_phase*
+create_phase_fn
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ const struct rnatm_cell_pos* cells) /* Cells in which scattering occurs */
+{
+ struct rnatm_sample_component_args sample_args =
+ RNATM_SAMPLE_COMPONENT_ARGS_NULL;
+ struct rnatm_cell_create_phase_fn_args phase_fn_args =
+ RNATM_CELL_CREATE_PHASE_FN_ARGS_NULL;
+ struct ssf_phase* phase = NULL;
+ size_t cpnt;
+ ASSERT(cmd && args && cells);
+
+ /* Sample the atmospheric scattering component */
+ sample_args.cells = cells;
+ sample_args.iband = args->iband;
+ sample_args.iquad = args->iquad;
+ sample_args.radcoef = RNATM_RADCOEF_Ks;
+ sample_args.r = ssp_rng_canonical(args->rng);
+ RNATM(sample_component(cmd->atmosphere, &sample_args, &cpnt));
+
+ /* Retrieve the component cell in which the scattering position is located */
+ phase_fn_args.cell = RNATM_GET_COMPONENT_CELL(cells, cpnt);
+ ASSERT(!SUVM_PRIMITIVE_NONE(&phase_fn_args.cell.prim));
+
+ /* Retrieve the component phase function */
+ phase_fn_args.wavelength = args->wlen;
+ phase_fn_args.r[0] = ssp_rng_canonical(args->rng);
+ phase_fn_args.r[1] = ssp_rng_canonical(args->rng);
+ RNATM(cell_create_phase_fn(cmd->atmosphere, &phase_fn_args, &phase));
+
+ return phase;
+}
+
+/* In shortwave, return the contribution of the external source at the bounce
+ * position. In longwave, simply return 0 */
+static double
+surface_bounce
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ const struct event* sc,
+ const double sc_pos[3], /* Scattering position */
+ const double in_dir[3], /* Incident direction */
+ double sc_dir[3], /* Sampled scattering direction */
+ double *out_refl) /* Surface reflectivity */
+{
+ struct ssf_bsdf* bsdf = NULL;
+ const double* N = NULL;
+ double wo[3] = {0,0,0};
+ double reflectivity = 0; /* Surface reflectivity */
+ double L = 0;
+ int mask = 0;
+ ASSERT(cmd && args && sc && sc_pos && in_dir && sc_dir && out_refl);
+ ASSERT(d3_dot(sc->normal, in_dir) < 0 && d3_is_normalized(sc->normal));
+
+ bsdf = create_bsdf(cmd, args, &sc->hit);
+ N = sc->normal;
+ d3_minus(wo, in_dir); /* Match StarSF convention */
+ ASSERT(d3_dot(wo, N) > 0);
+
+ /* Sample the scattering direction */
+ reflectivity = ssf_bsdf_sample(bsdf, args->rng, wo, N, sc_dir, &mask, NULL);
+
+ /* Fully absorbs transmissions */
+ if(mask & SSF_TRANSMISSION) reflectivity = 0;
+
+ /* No external source in longwave */
+ if(cmd->spectral_domain.type == HTRDR_SPECTRAL_LW)
+ goto exit;
+
+ /* Calculate direct contribution for specular reflection */
+ if((mask & SSF_SPECULAR)
+ && (mask & SSF_REFLECTION)
+ && htrdr_planeto_source_is_targeted(cmd->source, sc_pos, sc_dir)) {
+ const double Ld = direct_contribution(cmd, args, sc_pos, sc_dir, &sc->hit);
+ L = Ld * reflectivity;
+
+ /* Calculate direct contribution in general case */
+ } else if(!(mask & SSF_SPECULAR)) {
+ double wi[3], pdf;
+
+ /* Sample a direction toward the source */
+ pdf = htrdr_planeto_source_sample_direction(cmd->source, args->rng, sc_pos, wi);
+
+ /* The source is behind the surface */
+ if(d3_dot(wi, N) <= 0) {
+ L = 0;
+
+ /* The source is above the surface */
+ } else {
+ const double Ld = direct_contribution(cmd, args, sc_pos, wi, &sc->hit);
+ const double rho = ssf_bsdf_eval(bsdf, wo, N, wi);
+ const double cos_N_wi = d3_dot(N, wi);
+ ASSERT(cos_N_wi > 0);
+
+ L = Ld * rho * cos_N_wi / pdf;
+ }
+ }
+
+exit:
+ SSF(bsdf_ref_put(bsdf));
+ *out_refl = reflectivity;
+ return L;
+}
+
+/* In shortwave, return the contribution at the scattering position of the
+ * external source. Returns 0 in long wave */
+static INLINE double
+volume_scattering
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ const struct event* sc,
+ const double sc_pos[3], /* Scattering position */
+ const double in_dir[3], /* Incident direction */
+ double sc_dir[3]) /* Sampled scattering direction */
+{
+ struct ssf_phase* phase = NULL;
+ double wo[3] = {0,0,0};
+ double wi[3] = {0,0,0};
+ double L = 0;
+ double pdf = 0;
+ double rho = 0;
+ double Ld = 0;
+ ASSERT(cmd && args && sc && sc_pos && in_dir && sc_dir);
+
+ phase = create_phase_fn(cmd, args, sc->cells);
+ d3_minus(wo, in_dir); /* Match StarSF convention */
+
+ ssf_phase_sample(phase, args->rng, wo, sc_dir, NULL);
+
+ /* Sample a direction toward the source */
+ pdf = htrdr_planeto_source_sample_direction(cmd->source, args->rng, sc_pos, wi);
+
+ /* In short wave, manage the contribution of the external source */
+ switch(cmd->spectral_domain.type) {
+ case HTRDR_SPECTRAL_LW:
+ /* Nothing to be done */
+ break;
+
+ case HTRDR_SPECTRAL_SW:
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ /* Calculate the direct contribution at the scattering position */
+ Ld = direct_contribution(cmd, args, sc_pos, wi, NULL);
+ rho = ssf_phase_eval(phase, wo, wi);
+ L = Ld * rho / pdf;
+ break;
+
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ SSF(phase_ref_put(phase));
+ return L;
+}
+
+static INLINE enum rnatm_radcoef
+sample_volume_event_type
+ (const struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args,
+ struct event* evt)
+{
+ struct rnatm_get_radcoef_args get_k_args = RNATM_GET_RADCOEF_ARGS_NULL;
+ double ka, kext;
+ double r;
+ ASSERT(cmd && args && evt);
+
+ get_k_args.cells = evt->cells;
+ get_k_args.iband = args->iband;
+ get_k_args.iquad = args->iquad;
+
+ /* Retrieve the absorption coefficient */
+ get_k_args.radcoef = RNATM_RADCOEF_Ka;
+ RNATM(get_radcoef(cmd->atmosphere, &get_k_args, &ka));
+
+ /* Retrieve the extinction coefficient */
+ get_k_args.radcoef = RNATM_RADCOEF_Kext;
+ RNATM(get_radcoef(cmd->atmosphere, &get_k_args, &kext));
+
+ r = ssp_rng_canonical(args->rng);
+ if(r < ka / kext) {
+ return RNATM_RADCOEF_Ka; /* Absorption */
+ } else {
+ return RNATM_RADCOEF_Ks; /* Scattering */
+ }
+}
+
+static INLINE double
+get_temperature
+ (const struct htrdr_planeto* cmd,
+ const struct event* evt)
+{
+ double T = 0;
+ ASSERT(cmd && evt);
+
+ if(!SURFACE_EVENT(evt)) {
+ const struct rnatm_cell_pos* cell = NULL;
+
+ /* Get gas temperature */
+ cell = &RNATM_GET_COMPONENT_CELL(evt->cells, RNATM_GAS);
+ RNATM(cell_get_gas_temperature(cmd->atmosphere, cell, &T));
+
+ } else {
+ struct rngrd_get_temperature_args temp_args = RNGRD_GET_TEMPERATURE_ARGS_NULL;
+
+ /* Get ground temperature */
+ temp_args.prim = evt->hit.prim;
+ temp_args.barycentric_coords[0] = evt->hit.uv[0];
+ temp_args.barycentric_coords[1] = evt->hit.uv[1];
+ temp_args.barycentric_coords[2] = 1 - evt->hit.uv[0] - evt->hit.uv[1];
+ RNGRD(get_temperature(cmd->ground, &temp_args, &T));
+
+ }
+ return T;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+double /* Radiance in W/m²/sr/m */
+planeto_compute_radiance
+ (struct htrdr_planeto* cmd,
+ const struct planeto_compute_radiance_args* args)
+{
+ struct s3d_hit hit_from = S3D_HIT_NULL;
+ struct event ev;
+ double pos[3];
+ double dir[3];
+ double L = 0; /* Radiance in W/m²/sr/m */
+ size_t nsc = 0; /* Number of surface or volume scatterings (for debug) */
+ int longwave = 0;
+ ASSERT(cmd && check_planeto_compute_radiance_args(cmd, args) == RES_OK);
+
+ d3_set(pos, args->path_org);
+ d3_set(dir, args->path_dir);
+ longwave = cmd->spectral_domain.type == HTRDR_SPECTRAL_LW;
+
+ if(!longwave && htrdr_planeto_source_is_targeted(cmd->source, pos, dir)) {
+ L = direct_contribution(cmd, args, pos, dir, NULL); /* In W/m²/sr/m */
+ }
+
+ for(;;) {
+ double ev_pos[3];
+ double sc_dir[3];
+
+ find_event(cmd, args, RNATM_RADCOEF_Kext, pos, dir, &hit_from, &ev);
+
+ /* No event on surface or in volume. Stop the path */
+ if(DISTANCE_NONE(ev.distance)) break;
+
+ /* Compute the event position */
+ ev_pos[0] = pos[0] + dir[0] * ev.distance;
+ ev_pos[1] = pos[1] + dir[1] * ev.distance;
+ ev_pos[2] = pos[2] + dir[2] * ev.distance;
+
+ /* The event occurs on the surface */
+ if(SURFACE_EVENT(&ev)) {
+ double refl; /* Surface reflectivity */
+ L += surface_bounce(cmd, args, &ev, ev_pos, dir, sc_dir, &refl);
+
+ /* Check the absorption of the surface with a Russian roulette against
+ * the reflectivity of the surface */
+ if(ssp_rng_canonical(args->rng) >= refl) break;
+
+ /* Save current intersected surface to avoid self-intersection when
+ * sampling next event */
+ hit_from = ev.hit;
+
+ /* The event occurs in the volume */
+ } else {
+ enum rnatm_radcoef ev_type = sample_volume_event_type(cmd, args, &ev);
+ ASSERT(ev_type == RNATM_RADCOEF_Ka || ev_type == RNATM_RADCOEF_Ks);
+
+ /* Absorption. Stop the path */
+ if(ev_type == RNATM_RADCOEF_Ka) break;
+
+ L += volume_scattering(cmd, args, &ev, ev_pos, dir, sc_dir);
+ hit_from = S3D_HIT_NULL; /* Reset surface intersection */
+ }
+
+ d3_set(pos, ev_pos);
+ d3_set(dir, sc_dir);
+
+ ++nsc;
+ }
+
+ /* From there, a valid event means that the path has stopped in surface or
+ * volume absorption. In longwave, add the contribution of the internal
+ * source */
+ if(longwave && !DISTANCE_NONE(ev.distance)) {
+ const double wlen_m = args->wlen * 1.e-9; /* wavelength in meters */
+ const double temperature = get_temperature(cmd, &ev); /* In K */
+ L += htrdr_planck_monochromatic(wlen_m, temperature);
+ }
+
+ return L; /* In W/m²/sr/m */
+}
diff --git a/src/planeto/htrdr_planeto_draw_map.c b/src/planeto/htrdr_planeto_draw_map.c
@@ -0,0 +1,451 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, 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 "planeto/htrdr_planeto_c.h"
+#include "planeto/htrdr_planeto_source.h"
+
+#include "core/htrdr.h"
+#include "core/htrdr_accum.h"
+#include "core/htrdr_buffer.h"
+#include "core/htrdr_draw_map.h"
+#include "core/htrdr_log.h"
+#include "core/htrdr_ran_wlen_cie_xyz.h"
+#include "core/htrdr_ran_wlen_discrete.h"
+#include "core/htrdr_ran_wlen_planck.h"
+
+#include <rad-net/rnatm.h>
+#include <star/scam.h>
+#include <star/ssp.h>
+
+#include <rsys/clock_time.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+draw_pixel_xwave
+ (struct htrdr* htrdr,
+ const struct htrdr_draw_pixel_args* args,
+ void* data)
+{
+ struct planeto_compute_radiance_args rad_args =
+ PLANETO_COMPUTE_RADIANCE_ARGS_NULL;
+
+ struct htrdr_accum radiance;
+ struct htrdr_accum time;
+ struct htrdr_planeto* cmd;
+ struct planeto_pixel_xwave* pixel = data;
+ size_t isamp = 0;
+ ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data);
+ (void)htrdr;
+
+ cmd = args->context;
+ ASSERT(cmd);
+ ASSERT(cmd->spectral_domain.type == HTRDR_SPECTRAL_SW
+ || cmd->spectral_domain.type == HTRDR_SPECTRAL_LW);
+ ASSERT(cmd->output_type == HTRDR_PLANETO_ARGS_OUTPUT_IMAGE);
+
+ /* Reset accumulators */
+ radiance = HTRDR_ACCUM_NULL;
+ time = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(isamp, 0, args->spp) {
+ struct time t0, t1;
+ struct scam_sample sample = SCAM_SAMPLE_NULL;
+ struct scam_ray ray = SCAM_RAY_NULL;
+ double weight;
+ double r0, r1, r2;
+ double wlen[2]; /* Sampled wavelength */
+ double pdf;
+ size_t iband[2];
+ size_t iquad;
+ 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 */
+ sample.film[0] = (double)args->pixel_coord[0]+ssp_rng_canonical(args->rng);
+ sample.film[1] = (double)args->pixel_coord[1]+ssp_rng_canonical(args->rng);
+ sample.film[0] *= args->pixel_normalized_size[0];
+ sample.film[1] *= args->pixel_normalized_size[1];
+ sample.lens[0] = ssp_rng_canonical(args->rng);
+ sample.lens[1] = ssp_rng_canonical(args->rng);
+
+ /* Generate a camera ray */
+ scam_generate_ray(cmd->camera, &sample, &ray);
+
+ r0 = ssp_rng_canonical(args->rng);
+ r1 = ssp_rng_canonical(args->rng);
+ r2 = ssp_rng_canonical(args->rng);
+
+ /* Sample a wavelength */
+ switch(cmd->spectral_domain.type) {
+ case HTRDR_SPECTRAL_LW:
+ wlen[0] = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &pdf);
+ break;
+ case HTRDR_SPECTRAL_SW:
+ if(htrdr_planeto_source_does_radiance_vary_spectrally(cmd->source)) {
+ wlen[0] = htrdr_ran_wlen_discrete_sample(cmd->discrete, r0, r1, &pdf);
+ } else {
+ wlen[0] = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &pdf);
+ }
+ break;
+ default: FATAL("Unreachable code\n"); break;
+
+ }
+ wlen[1] = wlen[0];
+ pdf *= 1.e9; /* Transform the pdf from nm⁻¹ to m⁻¹ */
+
+ /* Find the band the wavelength belongs to and sample a quadrature point */
+ RNATM(find_bands(cmd->atmosphere, wlen, iband));
+ RNATM(band_sample_quad_pt(cmd->atmosphere, r2, iband[0], &iquad));
+ ASSERT(iband[0] == iband[1]);
+
+ /* Compute the radiance in W/m²/sr/m */
+ d3_set(rad_args.path_org, ray.org);
+ d3_set(rad_args.path_dir, ray.dir);
+ rad_args.rng = args->rng;
+ rad_args.ithread = args->ithread;
+ rad_args.wlen = wlen[0];
+ rad_args.iband = iband[0];
+ rad_args.iquad = iquad;
+ weight = planeto_compute_radiance(cmd, &rad_args);
+ ASSERT(weight >= 0);
+
+ weight /= pdf; /* In W/m²/sr */
+
+ /* End of time recording per realisation */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update the radiance accumulator */
+ radiance.sum_weights += weight;
+ radiance.sum_weights_sqr += weight*weight;
+ radiance.nweights += 1;
+
+ /* Update the per realisation time accumulator */
+ time.sum_weights += usec;
+ time.sum_weights_sqr += usec*usec;
+ time.nweights += 1;
+ }
+
+ /* Flush pixel data */
+ pixel->radiance = radiance;
+ pixel->time = time;
+
+ if(cmd->spectral_domain.type == HTRDR_SPECTRAL_SW) {
+ pixel->radiance_temperature.E = 0;
+ pixel->radiance_temperature.SE = 0;
+ } else {
+ struct htrdr_estimate L;
+
+ /* Wavelength in meters */
+ const double wmin_m = cmd->spectral_domain.wlen_range[0] * 1.e-9;
+ const double wmax_m = cmd->spectral_domain.wlen_range[1] * 1.e-9;
+
+ /* Brightness temperatures in W/m²/sr */
+ double T, Tmin, Tmax;
+
+ htrdr_accum_get_estimation(&radiance, &L);
+
+ T = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E);
+ Tmin = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E - L.SE);
+ Tmax = htrdr_radiance_temperature(cmd->htrdr, wmin_m, wmax_m, L.E + L.SE);
+
+ pixel->radiance_temperature.E = T;
+ pixel->radiance_temperature.SE = Tmax - Tmin;
+ }
+}
+
+static void
+draw_pixel_image
+ (struct htrdr* htrdr,
+ const struct htrdr_draw_pixel_args* args,
+ void* data)
+{
+ struct planeto_compute_radiance_args rad_args =
+ PLANETO_COMPUTE_RADIANCE_ARGS_NULL;
+
+ struct htrdr_accum XYZ[3]; /* X, Y, and Z */
+ struct htrdr_accum time;
+ struct htrdr_planeto* cmd;
+ struct planeto_pixel_image* pixel = data;
+ size_t ichannel;
+ ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data);
+ (void)htrdr;
+
+ cmd = args->context;
+ ASSERT(cmd);
+ ASSERT(cmd->spectral_domain.type == HTRDR_SPECTRAL_SW_CIE_XYZ);
+ ASSERT(cmd->output_type == HTRDR_PLANETO_ARGS_OUTPUT_IMAGE);
+
+ /* 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, args->spp) {
+ struct time t0, t1;
+ struct scam_sample sample = SCAM_SAMPLE_NULL;
+ struct scam_ray ray = SCAM_RAY_NULL;
+ double weight;
+ double r0, r1, r2;
+ double wlen[2]; /* Sampled wavelength into the spectral band */
+ double pdf;
+ size_t iband[2]; /* 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 */
+ sample.film[0] = (double)args->pixel_coord[0]+ssp_rng_canonical(args->rng);
+ sample.film[1] = (double)args->pixel_coord[1]+ssp_rng_canonical(args->rng);
+ sample.film[0] *= args->pixel_normalized_size[0];
+ sample.film[1] *= args->pixel_normalized_size[1];
+ sample.lens[0] = ssp_rng_canonical(args->rng);
+ sample.lens[1] = ssp_rng_canonical(args->rng);
+
+ /* Generate a camera ray */
+ SCAM(generate_ray(cmd->camera, &sample, &ray));
+
+ r0 = ssp_rng_canonical(args->rng);
+ r1 = ssp_rng_canonical(args->rng);
+ r2 = ssp_rng_canonical(args->rng);
+
+ /* Sample a wavelength according to the CIE XYZ color space */
+ switch(ichannel) {
+ case 0:
+ wlen[0] = htrdr_ran_wlen_cie_xyz_sample_X(cmd->cie, r0, r1, &pdf);
+ break;
+ case 1:
+ wlen[0] = htrdr_ran_wlen_cie_xyz_sample_Y(cmd->cie, r0, r1, &pdf);
+ break;
+ case 2:
+ wlen[0] = htrdr_ran_wlen_cie_xyz_sample_Z(cmd->cie, r0, r1, &pdf);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ pdf *= 1.e9; /* Transform the pdf from nm⁻¹ to m⁻¹ */
+
+ /* Find the band the wavelength belongs to and sample a quadrature point */
+ wlen[1] = wlen[0];
+ RNATM(find_bands(cmd->atmosphere, wlen, iband));
+ RNATM(band_sample_quad_pt(cmd->atmosphere, r2, iband[0], &iquad));
+ ASSERT(iband[0] == iband[1]);
+
+ /* Compute the radiance in W/m²/sr/m */
+ d3_set(rad_args.path_org, ray.org);
+ d3_set(rad_args.path_dir, ray.dir);
+ rad_args.rng = args->rng;
+ rad_args.ithread = args->ithread;
+ rad_args.wlen = wlen[0];
+ rad_args.iband = iband[0];
+ rad_args.iquad = iquad;
+ weight = planeto_compute_radiance(cmd, &rad_args);
+ ASSERT(weight >= 0);
+
+ weight /= pdf; /* In W/m²/sr */
+
+ /* End of time recording per realisation */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update pixel channel accumulators */
+ XYZ[ichannel].sum_weights += weight;
+ XYZ[ichannel].sum_weights_sqr += weight*weight;
+ XYZ[ichannel].nweights += 1;
+
+ /* Update the per realisation time accumulator */
+ time.sum_weights += usec;
+ time.sum_weights_sqr += usec*usec;
+ time.nweights += 1;
+ }
+ }
+
+ /* Flush pixel data */
+ htrdr_accum_get_estimation(XYZ+0, &pixel->X);
+ htrdr_accum_get_estimation(XYZ+1, &pixel->Y);
+ htrdr_accum_get_estimation(XYZ+2, &pixel->Z);
+ pixel->time = time;
+}
+
+
+static INLINE void
+write_accum
+ (const struct htrdr_accum* acc, /* Accum to dump */
+ struct htrdr_accum* out_acc, /* May be NULL */
+ FILE* stream)
+{
+ ASSERT(acc && stream);
+
+ if(acc->nweights == 0) {
+ fprintf(stream, "0 0 ");
+ } else {
+ struct htrdr_estimate estimate = HTRDR_ESTIMATE_NULL;
+
+ htrdr_accum_get_estimation(acc, &estimate);
+ fprintf(stream, "%g %g ", estimate.E, estimate.SE);
+
+ if(out_acc) {
+ out_acc->sum_weights += acc->sum_weights;
+ out_acc->sum_weights_sqr += acc->sum_weights_sqr;
+ out_acc->nweights += acc->nweights;
+ }
+ }
+}
+
+static INLINE void
+write_pixel_image
+ (const struct planeto_pixel_image* pix,
+ struct htrdr_accum* time_acc, /* May be NULL */
+ FILE* stream)
+{
+ ASSERT(pix && stream);
+ 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);
+ write_accum(&pix->time, time_acc, stream);
+ fprintf(stream, "\n");
+}
+
+static INLINE void
+write_pixel_xwave
+ (const struct planeto_pixel_xwave* pix,
+ struct htrdr_accum* radiance_acc, /* May be NULL */
+ struct htrdr_accum* time_acc, /* May be NULL */
+ FILE* stream)
+{
+ ASSERT(pix && stream);
+ fprintf(stream, "%g %g ",
+ pix->radiance_temperature.E,
+ pix->radiance_temperature.SE);
+ write_accum(&pix->radiance, radiance_acc, stream);
+ fprintf(stream, "0 0 ");
+ write_accum(&pix->time, time_acc, stream);
+ fprintf(stream, "\n");
+}
+
+static res_T
+write_buffer
+ (struct htrdr_planeto* cmd,
+ struct htrdr_buffer* buf,
+ struct htrdr_accum* radiance_acc, /* May be NULL */
+ struct htrdr_accum* time_acc,
+ FILE* stream)
+{
+ struct htrdr_pixel_format pixfmt;
+ struct htrdr_buffer_layout layout;
+ size_t x, y;
+ ASSERT(cmd && buf && time_acc && stream);
+ ASSERT(cmd->output_type == HTRDR_PLANETO_ARGS_OUTPUT_IMAGE);
+
+ planeto_get_pixel_format(cmd, &pixfmt);
+ htrdr_buffer_get_layout(buf, &layout);
+ CHK(pixfmt.size == layout.elmt_size);
+
+ fprintf(stream, "%lu %lu\n", layout.width, layout.height);
+ *time_acc = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(y, 0, layout.height) {
+ FOR_EACH(x, 0, layout.width) {
+ void* pix_raw = htrdr_buffer_at(buf, x, y);
+ ASSERT(IS_ALIGNED(pix_raw, pixfmt.alignment));
+
+ switch(cmd->spectral_domain.type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ write_pixel_xwave(pix_raw, radiance_acc, time_acc, stream);
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ write_pixel_image(pix_raw, time_acc, stream);
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ }
+ }
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+planeto_draw_map(struct htrdr_planeto* cmd)
+{
+ struct htrdr_draw_map_args args = HTRDR_DRAW_MAP_ARGS_NULL;
+ struct htrdr_estimate path_time = HTRDR_ESTIMATE_NULL;
+ struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL;
+ struct htrdr_accum radiance_acc = HTRDR_ACCUM_NULL;
+ res_T res = RES_OK;
+ ASSERT(cmd && cmd->output_type == HTRDR_PLANETO_ARGS_OUTPUT_IMAGE);
+
+ args.buffer_layout = cmd->buf_layout;
+ args.spp = cmd->spp;
+ args.context = cmd;
+ switch(cmd->spectral_domain.type) {
+ case HTRDR_SPECTRAL_LW:
+ case HTRDR_SPECTRAL_SW:
+ args.draw_pixel = draw_pixel_xwave;
+ break;
+ case HTRDR_SPECTRAL_SW_CIE_XYZ:
+ args.draw_pixel = draw_pixel_image;
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
+ res = htrdr_draw_map(cmd->htrdr, &args, cmd->buf);
+ if(res != RES_OK) goto error;
+
+ /* Nothing more to do on non master processes */
+ if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
+
+ /* Write output image */
+ res = write_buffer(cmd, cmd->buf, &radiance_acc, &path_time_acc, cmd->output);
+ if(res != RES_OK) goto error;
+
+ CHK(fflush(cmd->output) == 0);
+
+ /* Log time per realisation */
+ htrdr_accum_get_estimation(&path_time_acc, &path_time);
+ htrdr_log(cmd->htrdr, "Time per radiative path (in µs): %g +/- %g\n",
+ path_time.E, path_time.SE);
+
+ /* Log measured radiance on the whole image */
+ if(cmd->spectral_domain.type == HTRDR_SPECTRAL_LW
+ || cmd->spectral_domain.type == HTRDR_SPECTRAL_SW) {
+ struct htrdr_estimate L;
+ double omega; /* Solid angle of the camera */
+
+ htrdr_accum_get_estimation(&radiance_acc, &L);
+ SCAM(perspective_get_solid_angle(cmd->camera, &omega));
+ htrdr_log(cmd->htrdr, "Radiance in W/m²/sr: %g +/- %g\n", L.E, L.SE);
+ htrdr_log(cmd->htrdr, "Radiance in W/m² (solid angle = %g sr): %g +/- %g\n",
+ omega, L.E*omega, L.SE*omega);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/planeto/htrdr_planeto_main.c b/src/planeto/htrdr_planeto_main.c
@@ -0,0 +1,85 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, 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 "planeto/htrdr_planeto.h"
+#include "planeto/htrdr_planeto_args.h"
+
+#include "core/htrdr_log.h"
+
+#include <rsys/mem_allocator.h>
+
+int
+htrdr_planeto_main(int argc, char** argv)
+{
+ char cmd_name[] = "htrdr-planeto";
+ struct htrdr_args htrdr_args = HTRDR_ARGS_DEFAULT;
+ struct htrdr_planeto_args cmd_args = HTRDR_PLANETO_ARGS_DEFAULT;
+ struct htrdr* htrdr = NULL;
+ struct htrdr_planeto* cmd = NULL;
+ const size_t memsz_begin = mem_allocated_size();
+ size_t memsz_end;
+ int is_mpi_init = 0;
+ res_T res = RES_OK;
+ int err = 0;
+
+ /* Overwrite command name */
+ argv[0] = cmd_name;
+
+ res = htrdr_mpi_init(argc, argv);
+ if(res != RES_OK) goto error;
+ is_mpi_init = 1;
+
+ res = htrdr_planeto_args_init(&cmd_args, argc, argv);
+ if(res != RES_OK) goto error;
+ if(cmd_args.quit) goto exit;
+
+ htrdr_args.nthreads = cmd_args.nthreads;
+ htrdr_args.verbose = cmd_args.verbose;
+ res = htrdr_create(&mem_default_allocator, &htrdr_args, &htrdr);
+ if(res != RES_OK) goto error;
+
+ if(cmd_args.output_type == HTRDR_PLANETO_ARGS_OUTPUT_OCTREES
+ && htrdr_get_mpi_rank(htrdr) != 0) {
+ goto exit; /* Nothing to do except for the master process */
+ }
+
+ res = htrdr_planeto_create(htrdr, &cmd_args, &cmd);
+ if(res != RES_OK) goto error;
+
+ res = htrdr_planeto_run(cmd);
+ if(res != RES_OK) goto error;
+
+exit:
+ htrdr_planeto_args_release(&cmd_args);
+ if(is_mpi_init) htrdr_mpi_finalize();
+ if(htrdr) htrdr_ref_put(htrdr);
+ if(cmd) htrdr_planeto_ref_put(cmd);
+
+ /* Check memory leaks */
+ memsz_end = mem_allocated_size();
+ if(memsz_begin != memsz_end) {
+ ASSERT(memsz_end >= memsz_begin);
+ fprintf(stderr, HTRDR_LOG_WARNING_PREFIX"Memory leaks: %lu Bytes\n",
+ (unsigned long)(memsz_end - memsz_begin));
+ err = -1;
+ }
+ return err;
+error:
+ err = -1;
+ goto exit;
+}
+
diff --git a/src/planeto/htrdr_planeto_source.c b/src/planeto/htrdr_planeto_source.c
@@ -0,0 +1,485 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, 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 "planeto/htrdr_planeto_c.h"
+#include "planeto/htrdr_planeto_source.h"
+
+#include "core/htrdr.h"
+#include "core/htrdr_log.h"
+
+#include <star/sbuf.h>
+#include <star/ssp.h>
+
+#include <rsys/algorithm.h>
+#include <rsys/cstr.h>
+#include <rsys/double3.h>
+#include <rsys/ref_count.h>
+
+typedef struct ALIGN(16) {
+ double wavelength; /* in nm */
+ double radiance; /* in W/m²/sr/m */
+} source_radiance_T;
+
+struct htrdr_planeto_source {
+ double position[3]; /* In m */
+
+ double radius; /* In m */
+
+ /* In Kelvin. Defined if the radiances by wavelength is no set */
+ double temperature;
+
+ struct sbuf* per_wlen_radiances; /* List of radiances by wavelength */
+
+ ref_T ref;
+ struct htrdr* htrdr;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE res_T
+check_per_wlen_radiance_sbuf_desc
+ (const struct htrdr_planeto_source* src,
+ const struct sbuf_desc* desc)
+{
+ const source_radiance_T* spectrum = NULL;
+ size_t i;
+ ASSERT(src && desc);
+
+ /* Invalid size */
+ if(desc->size == 0) {
+ htrdr_log_err(src->htrdr, "invalid empty source spectrum\n");
+ return RES_BAD_ARG;
+ }
+
+ /* Invalid memory layout */
+ if(desc->szitem != 16 || desc->alitem != 16 || desc->pitch != 16) {
+ htrdr_log_err(src->htrdr, "unexpected layout of source spectrum\n");
+ return RES_BAD_ARG;
+ }
+
+ /* Data must be sorted */
+ spectrum = desc->buffer;
+ FOR_EACH(i, 1, desc->size) {
+ if(spectrum[i-1].wavelength >= spectrum[i].wavelength) {
+ htrdr_log_err(src->htrdr,
+ "the source spectrum is not sorted in ascending order "
+ "with respect to wavelengths\n");
+ return RES_BAD_ARG;
+ }
+ }
+
+ return RES_OK;
+}
+
+static res_T
+setup_per_wavelength_radiances
+ (struct htrdr_planeto_source* src,
+ const struct htrdr_planeto_source_args* args)
+{
+ struct sbuf_create_args sbuf_args;
+ struct sbuf_desc desc;
+ res_T res = RES_OK;
+ ASSERT(src && args && args->rnrl_filename && args->temperature < 0);
+
+ sbuf_args.logger = htrdr_get_logger(src->htrdr);
+ sbuf_args.allocator = htrdr_get_allocator(src->htrdr);
+ sbuf_args.verbose = htrdr_get_verbosity_level(src->htrdr);
+ res = sbuf_create(&sbuf_args, &src->per_wlen_radiances);
+ if(res != RES_OK) goto error;
+
+ res = sbuf_load(src->per_wlen_radiances, args->rnrl_filename);
+ if(res != RES_OK) goto error;
+ res = sbuf_get_desc(src->per_wlen_radiances, &desc);
+ if(res != RES_OK) goto error;
+ res = check_per_wlen_radiance_sbuf_desc(src, &desc);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ htrdr_log_err(src->htrdr, "error loading %s -- %s\n",
+ args->rnrl_filename, res_to_cstr(res));
+ goto exit;
+}
+
+static INLINE int
+cmp_wlen(const void* a, const void* b)
+{
+ const double wlen = *((double*)a);
+ const source_radiance_T* src_rad1 = b;
+ ASSERT(a && b);
+
+ if(wlen < src_rad1->wavelength) {
+ return -1;
+ } else if(wlen > src_rad1->wavelength) {
+ return +1;
+ } else {
+ return 0;
+ }
+}
+
+static double
+get_radiance
+ (const struct htrdr_planeto_source* src,
+ const double wlen)
+{
+ struct sbuf_desc desc;
+ const source_radiance_T* spectrum;
+ const source_radiance_T* find;
+ size_t id;
+ ASSERT(src && src->per_wlen_radiances);
+
+ SBUF(get_desc(src->per_wlen_radiances, &desc));
+ spectrum = desc.buffer;
+
+ if(wlen < spectrum[0].wavelength) {
+ htrdr_log_warn(src->htrdr,
+ "The wavelength %g nm is before the spectrum of the source\n", wlen);
+ return spectrum[0].radiance;
+ }
+ if(wlen > spectrum[desc.size-1].wavelength) {
+ htrdr_log_warn(src->htrdr,
+ "The wavelength %g nm is above the spectrum of the source\n", wlen);
+ return spectrum[desc.size-1].radiance;
+ }
+
+ /* Look for the first item whose wavelength is not less than 'wlen' */
+ find = search_lower_bound(&wlen, spectrum, desc.size, desc.pitch, cmp_wlen);
+ ASSERT(find);
+ id = (size_t)(find - spectrum);
+ ASSERT(id < desc.size);
+
+ if(id == 0) {
+ ASSERT(wlen == spectrum[0].wavelength);
+ return spectrum[0].radiance;
+ } else {
+ const double w0 = spectrum[id-1].wavelength;
+ const double w1 = spectrum[id-0].wavelength;
+ const double L0 = spectrum[id-1].radiance;
+ const double L1 = spectrum[id-0].radiance;
+ const double u = (wlen - w0) / (w1 - w0);
+ const double L = L0 + u*(L1 - L0); /* Linear interpolation */
+ return L;
+ }
+}
+
+static void
+release_source(ref_T* ref)
+{
+ struct htrdr_planeto_source* source;
+ struct htrdr* htrdr;
+ ASSERT(ref);
+
+ source = CONTAINER_OF(ref, struct htrdr_planeto_source, ref);
+ htrdr = source->htrdr;
+ if(source->per_wlen_radiances) SBUF(ref_put(source->per_wlen_radiances));
+ MEM_RM(htrdr_get_allocator(htrdr), source);
+ htrdr_ref_put(htrdr);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_planeto_source_create
+ (struct htrdr* htrdr,
+ const struct htrdr_planeto_source_args* args,
+ struct htrdr_planeto_source** out_source)
+{
+ struct htrdr_planeto_source* src = NULL;
+ double dst; /* In m */
+ double lat; /* In radians */
+ double lon; /* In radians */
+ res_T res = RES_OK;
+ ASSERT(htrdr && out_source);
+ ASSERT(htrdr_planeto_source_args_check(args) == RES_OK);
+
+ src = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*src));
+ if(!src) {
+ htrdr_log_err(htrdr, "error allocating source\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&src->ref);
+ htrdr_ref_get(htrdr);
+ src->htrdr = htrdr;
+ src->radius = args->radius * 1e3/*From km to m*/;
+
+ if(!args->rnrl_filename) {
+ src->temperature = args->temperature;
+ } else {
+ res = setup_per_wavelength_radiances(src, args);
+ if(res != RES_OK) goto error;
+ src->temperature = -1; /* Not used */
+ }
+
+ /* Convert latitude and longitude to radians and distance in m */
+ lat = MDEG2RAD(args->latitude);
+ lon = MDEG2RAD(args->longitude);
+ dst = args->distance * 1e3/*From km to m*/;
+
+ /* Compute the position of the source */
+ src->position[0] = dst * cos(lat) * cos(lon);
+ src->position[1] = dst * cos(lat) * sin(lon);
+ src->position[2] = dst * sin(lat);
+
+exit:
+ *out_source = src;
+ return res;
+error:
+ if(src) { htrdr_planeto_source_ref_put(src); src = NULL; }
+ goto exit;
+}
+
+void
+htrdr_planeto_source_ref_get(struct htrdr_planeto_source* source)
+{
+ ASSERT(source);
+ ref_get(&source->ref);
+}
+
+void htrdr_planeto_source_ref_put(struct htrdr_planeto_source* source)
+{
+ ASSERT(source);
+ ref_put(&source->ref, release_source);
+}
+
+double
+htrdr_planeto_source_sample_direction
+ (const struct htrdr_planeto_source* source,
+ struct ssp_rng* rng,
+ const double pos[3],
+ double dir[3])
+{
+ double main_dir[3];
+ double half_angle; /* In radians */
+ double cos_half_angle;
+ double dst; /* In m */
+ double pdf;
+ ASSERT(source && rng && pos && dir);
+
+ /* compute the direction of `pos' toward the center of the source */
+ d3_sub(main_dir, source->position, pos);
+
+ /* Normalize the direction and keep the distance from `pos' to the center of
+ * the source */
+ dst = d3_normalize(main_dir, main_dir);
+ CHK(dst > source->radius);
+
+ /* Sample the source according to its solid angle,
+ * i.e. 2*PI*(1 - cos(half_angle)) */
+ half_angle = asin(source->radius/dst);
+ cos_half_angle = cos(half_angle);
+ ssp_ran_sphere_cap_uniform(rng, main_dir, cos_half_angle, dir, &pdf);
+
+ return pdf;
+}
+
+double /* In W/m²/sr/m */
+htrdr_planeto_source_get_radiance
+ (const struct htrdr_planeto_source* source,
+ const double wlen)
+{
+ if(source->per_wlen_radiances) {
+ return get_radiance(source, wlen);
+ } else {
+ return htrdr_planck_monochromatic
+ (wlen*1e-9/*From nm to m*/, source->temperature);
+ }
+}
+
+double
+htrdr_planeto_source_distance_to
+ (const struct htrdr_planeto_source* source,
+ const double pos[3])
+{
+ double vec[3];
+ double dst;
+ ASSERT(source && pos);
+
+ d3_sub(vec, source->position, pos);
+ dst = d3_len(vec);
+ return dst - source->radius;
+}
+
+int
+htrdr_planeto_source_is_targeted
+ (const struct htrdr_planeto_source* source,
+ const double pos[3],
+ const double dir[3])
+{
+ double main_dir[3];
+ double half_angle; /* In radians */
+ double dst; /* In m */
+ ASSERT(source && dir && d3_is_normalized(dir));
+
+ /* compute the direction of `pos' toward the center of the source */
+ d3_sub(main_dir, source->position, pos);
+
+ /* Normalize the direction and keep the distance from `pos' to the center of
+ * the source */
+ dst = d3_normalize(main_dir, main_dir);
+ CHK(dst > source->radius);
+
+ /* Compute the the half angle of the source as seen from pos */
+ half_angle = asin(source->radius/dst);
+
+ return d3_dot(dir, main_dir) >= cos(half_angle);
+}
+
+res_T
+htrdr_planeto_source_get_spectral_range
+ (const struct htrdr_planeto_source* source,
+ double range[2])
+{
+ res_T res = RES_OK;
+ ASSERT(source && range);
+
+ if(!source->per_wlen_radiances) {
+ range[0] = 0;
+ range[1] = INF;
+ } else {
+ struct sbuf_desc desc = SBUF_DESC_NULL;
+ const source_radiance_T* spectrum = NULL;
+
+ res = sbuf_get_desc(source->per_wlen_radiances, &desc);
+ if(res != RES_OK) goto error;
+
+ spectrum = desc.buffer;
+ range[0] = spectrum[0].wavelength;
+ range[1] = spectrum[desc.size-1].wavelength;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+int
+htrdr_planeto_source_does_radiance_vary_spectrally
+ (const struct htrdr_planeto_source* source)
+{
+ ASSERT(source);
+ return source->per_wlen_radiances != NULL;
+}
+
+res_T
+htrdr_planeto_source_get_spectrum
+ (const struct htrdr_planeto_source* source,
+ const double range[2], /* In nm. Limits are inclusive */
+ struct htrdr_planeto_source_spectrum* source_spectrum)
+{
+ double full_range[2];
+ res_T res = RES_OK;
+ ASSERT(source && range && source_spectrum && range[0] <= range[1]);
+
+ if(!htrdr_planeto_source_does_radiance_vary_spectrally(source)) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = htrdr_planeto_source_get_spectral_range(source, full_range);
+ if(res != RES_OK) goto error;
+
+ if(range[0] < full_range[0] || full_range[1] < range[1]) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ source_spectrum->source = source;
+ source_spectrum->range[0] = range[0];
+ source_spectrum->range[1] = range[1];
+
+ if(range[0] == range[1]) {
+ /* Degenerated spectral range */
+ source_spectrum->size = 1;
+ source_spectrum->buffer = NULL;
+
+ } else {
+ const source_radiance_T* spectrum;
+ const source_radiance_T* low;
+ const source_radiance_T* upp;
+ struct sbuf_desc desc;
+
+ res = sbuf_get_desc(source->per_wlen_radiances, &desc);
+ if(res != RES_OK) goto error;
+
+ spectrum = desc.buffer;
+ low = search_lower_bound(&range[0], spectrum, desc.size, desc.pitch, cmp_wlen);
+ upp = search_lower_bound(&range[1], spectrum, desc.size, desc.pitch, cmp_wlen);
+ ASSERT(low && upp);
+
+ if(low == upp) {
+ /* The range is fully included in a band */
+ ASSERT(low->radiance > range[0] && upp->radiance >= range[1]);
+ source_spectrum->size = 2;
+ source_spectrum->buffer = NULL;
+
+ } else {
+ source_spectrum->size =
+ 2/* Boundaries */ + (size_t)(upp - low)/*discrete items*/;
+
+ if(low->wavelength == range[0]) {
+ /* The lower limit coincide with a discrete element.
+ * Remove the discrete element */
+ source_spectrum->size -= 1;
+ source_spectrum->buffer = low + 1;
+ } else {
+ source_spectrum->buffer = low;
+ }
+
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+void
+htrdr_planeto_source_spectrum_at
+ (void* source_spectrum,
+ const size_t i, /* between [0, spectrum->size[ */
+ double* wavelength, /* In nm */
+ double* radiance) /* In W/m²/sr/m */
+{
+ struct htrdr_planeto_source_spectrum* spectrum = source_spectrum;
+ ASSERT(spectrum && i < spectrum->size && wavelength && radiance);
+
+ /* Lower limit */
+ if(i == 0) {
+ *wavelength = spectrum->range[0];
+ *radiance = htrdr_planeto_source_get_radiance
+ (spectrum->source, spectrum->range[0]);
+
+ /* Upper limit */
+ } else if(i == spectrum->size-1) {
+ *wavelength = spectrum->range[1];
+ *radiance = htrdr_planeto_source_get_radiance
+ (spectrum->source, spectrum->range[1]);
+
+ /* Discrete element */
+ } else {
+ const source_radiance_T* item =
+ (const source_radiance_T*)spectrum->buffer + (i-1);
+ *wavelength = item->wavelength;
+ *radiance = item->radiance;
+ }
+}
diff --git a/src/planeto/htrdr_planeto_source.h b/src/planeto/htrdr_planeto_source.h
@@ -0,0 +1,110 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, Université Paul Sabatier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef HTRDR_PLANETO_SOURCE_H
+#define HTRDR_PLANETO_SOURCE_H
+
+#include <rsys/rsys.h>
+
+/* Forward declarations */
+struct htrdr;
+struct htrdr_planeto_source;
+struct htrdr_planeto_source_args;
+struct ssp_rng;
+
+struct htrdr_planeto_source_spectrum {
+ const struct htrdr_planeto_source* source;
+ double range[2]; /* In nm. Limits are inclusive */
+ size_t size; /* Number of elements representing the spectrum */
+ const void* buffer; /* Pointer toward the spectrum data */
+};
+#define HTRDR_PLANETO_SOURCE_SPECTRUM_NULL__ {NULL, {0,0}, 0, NULL}
+static const struct htrdr_planeto_source_spectrum
+HTRDR_PLANETO_SOURCE_SPECTRUM_NULL = HTRDR_PLANETO_SOURCE_SPECTRUM_NULL__;
+
+extern LOCAL_SYM res_T
+htrdr_planeto_source_create
+ (struct htrdr* htrdr,
+ const struct htrdr_planeto_source_args* args,
+ struct htrdr_planeto_source** source);
+
+extern LOCAL_SYM void
+htrdr_planeto_source_ref_get
+ (struct htrdr_planeto_source* source);
+
+extern LOCAL_SYM void
+htrdr_planeto_source_ref_put
+ (struct htrdr_planeto_source* source);
+
+/* Return the pdf of the sampled direction */
+extern LOCAL_SYM double
+htrdr_planeto_source_sample_direction
+ (const struct htrdr_planeto_source* source,
+ struct ssp_rng* rng,
+ const double pos[3], /* Position from which direction is sampled */
+ double dir[3]);
+
+extern LOCAL_SYM double /* In W/m²/sr/m */
+htrdr_planeto_source_get_radiance
+ (const struct htrdr_planeto_source* source,
+ const double wlen); /* In nanometers */
+
+/* Return the distance between the source surface and the input position. Can
+ * be negative if the position is in the source */
+extern LOCAL_SYM double /* In m */
+htrdr_planeto_source_distance_to
+ (const struct htrdr_planeto_source* source,
+ const double pos[3]);
+
+/* Return 1 if the source is targeted by the submitted ray and 0 otherwise */
+extern LOCAL_SYM int
+htrdr_planeto_source_is_targeted
+ (const struct htrdr_planeto_source* source,
+ const double pos[3], /* Ray origin */
+ const double dir[3]);/* Ray direction */
+
+extern LOCAL_SYM res_T
+htrdr_planeto_source_get_spectral_range
+ (const struct htrdr_planeto_source* source,
+ double range[2]); /* In nm. Limits are inclusive */
+
+extern LOCAL_SYM int
+htrdr_planeto_source_does_radiance_vary_spectrally
+ (const struct htrdr_planeto_source* source);
+
+/* Get discrete spectrum data for a given range. If the boundaries of the
+ * spectral range do not coincide with a discrete element, their radiance is
+ * recovered from the htrdr_planeto_source_get_radiance function. Note that
+ * this function returns an error if the radiance from the source does not vary
+ * spectrally, that is, its radiance is recovered from a constant temperature */
+extern LOCAL_SYM res_T
+htrdr_planeto_source_get_spectrum
+ (const struct htrdr_planeto_source* source,
+ const double range[2], /* In nm. Limits are inclusive */
+ struct htrdr_planeto_source_spectrum* spectrum);
+
+/* Note that the following function profile corresponds to the type expected by
+ * the discrete wavelength distribution
+ * (see htrdr_ran_wlen_discrete_create_args structure) */
+extern LOCAL_SYM void
+htrdr_planeto_source_spectrum_at
+ (void* spectrum,
+ size_t i, /* between [0, spectrum->size[ */
+ double* wavelength, /* In nm */
+ double* radiance); /* In W/m²/sr/m */
+
+#endif /* HTRDR_PLANETO_SOURCE_H */
diff --git a/src/planeto/test_htrdr_planeto_source.c b/src/planeto/test_htrdr_planeto_source.c
@@ -0,0 +1,296 @@
+/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2018, 2019, 2021 CNRS
+ * Copyright (C) 2018, 2019, 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 "planeto/htrdr_planeto_args.h"
+#include "planeto/htrdr_planeto_source.h"
+
+#include "core/htrdr.h"
+#include "core/htrdr_ran_wlen_discrete.h"
+
+#include <rsys/math.h>
+#include <rsys/mem_allocator.h>
+
+#include <stdio.h>
+
+static void
+write_per_wlen_radiances
+ (FILE* fp,
+ const size_t pagesize,
+ const size_t size,
+ const size_t szelmt,
+ const size_t alelmt)
+{
+ const char byte = 0;
+ size_t i;
+
+ CHK(fp);
+
+ /* Header */
+ CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1);
+ CHK(fwrite(&size, sizeof(size), 1, fp) == 1);
+ CHK(fwrite(&szelmt, sizeof(szelmt), 1, fp) == 1);
+ CHK(fwrite(&alelmt, sizeof(alelmt), 1, fp) == 1);
+
+ /* Padding */
+ CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET) == 0);
+
+ FOR_EACH(i, 0, size) {
+ const double w = (double)i;
+ const double L = (double)(100 + i);
+
+ CHK(fwrite(&w, sizeof(w), 1, fp) == 1);
+ CHK(fwrite(&L, sizeof(L), 1, fp) == 1);
+ }
+
+ /* Padding. Write one char to position the EOF indicator */
+ CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0);
+ CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1);
+
+ CHK(fflush(fp) == 0);
+}
+
+static void
+test_spectrum(struct htrdr* htrdr)
+{
+ struct htrdr_planeto_source_args source_args = HTRDR_PLANETO_SOURCE_ARGS_NULL;
+ struct htrdr_planeto_source_spectrum spectrum;
+ struct htrdr_planeto_source* source = NULL;
+
+ FILE* fp = NULL;
+ char rnrl_filename[] = "rnrl.bin";
+ double range[2];
+ double w, L;
+
+ CHK(fp = fopen(rnrl_filename, "w"));
+ write_per_wlen_radiances(fp, 4096, 10, 16, 16);
+ CHK(fclose(fp) == 0);
+
+ source_args.rnrl_filename = rnrl_filename;
+ source_args.longitude = 0;
+ source_args.latitude = 0;
+ source_args.distance = 0;
+ source_args.radius = 1e8;
+ source_args.temperature = -1;
+ CHK(htrdr_planeto_source_create(htrdr, &source_args, &source) == RES_OK);
+ CHK(htrdr_planeto_source_does_radiance_vary_spectrally(source) == 1);
+ CHK(htrdr_planeto_source_get_spectral_range(source, range) == RES_OK);
+ CHK(range[0] == 0);
+ CHK(range[1] == 9);
+
+ range[0] = 0; range[1] = 10;
+ CHK(htrdr_planeto_source_get_spectrum(source, range, &spectrum) == RES_BAD_ARG);
+
+ range[0] = 1; range[1] = 3;
+ CHK(htrdr_planeto_source_get_spectrum(source, range, &spectrum) == RES_OK);
+ CHK(spectrum.source == source);
+ CHK(spectrum.range[0] == 1);
+ CHK(spectrum.range[1] == 3);
+ CHK(spectrum.size == 3);
+
+ htrdr_planeto_source_spectrum_at(&spectrum, 0, &w, &L);
+ CHK(w == 1 && L == 101);
+ htrdr_planeto_source_spectrum_at(&spectrum, 1, &w, &L);
+ CHK(w == 2 && L == 102);
+ htrdr_planeto_source_spectrum_at(&spectrum, 2, &w, &L);
+ CHK(w == 3 && L == 103);
+
+ range[0] = 1.7; range[1] = 1.95;
+ CHK(htrdr_planeto_source_get_spectrum(source, range, &spectrum) == RES_OK);
+ CHK(spectrum.source == source);
+ CHK(spectrum.range[0] = 1.7);
+ CHK(spectrum.range[1] = 1.95);
+ CHK(spectrum.size == 2);
+ htrdr_planeto_source_spectrum_at(&spectrum, 0, &w, &L);
+ CHK(w == 1.7 && eq_eps(L, 101.7, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 1, &w, &L);
+ CHK(w == 1.95 && eq_eps(L, 101.95, 1.e-6));
+
+ range[0] = 2; range[1] = 2.01;
+ CHK(htrdr_planeto_source_get_spectrum(source, range, &spectrum) == RES_OK);
+ CHK(spectrum.size == 2);
+ htrdr_planeto_source_spectrum_at(&spectrum, 0, &w, &L);
+ CHK(w == 2 && L == 102);
+ htrdr_planeto_source_spectrum_at(&spectrum, 1, &w, &L);
+ CHK(w == 2.01 && eq_eps(L, 102.01, 1.e-6));
+
+ range[0] = 5.1; range[1] = 6;
+ CHK(htrdr_planeto_source_get_spectrum(source, range, &spectrum) == RES_OK);
+ CHK(spectrum.size == 2);
+ htrdr_planeto_source_spectrum_at(&spectrum, 0, &w, &L);
+ CHK(w == 5.1 && eq_eps(L, 105.1, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 1, &w, &L);
+ CHK(w == 6 && L == 106);
+
+ range[0] = 7.5; range[1] = 9;
+ CHK(htrdr_planeto_source_get_spectrum(source, range, &spectrum) == RES_OK);
+ CHK(spectrum.size == 3);
+ htrdr_planeto_source_spectrum_at(&spectrum, 0, &w, &L);
+ CHK(w == 7.5 && eq_eps(L, 107.5, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 1, &w, &L);
+ CHK(w == 8 && L == 108);
+ htrdr_planeto_source_spectrum_at(&spectrum, 2, &w, &L);
+ CHK(w == 9 && L == 109);
+
+ range[0] = 0.9; range[1] = 7.456;
+ CHK(htrdr_planeto_source_get_spectrum(source, range, &spectrum) == RES_OK);
+ CHK(spectrum.size == 9);
+ htrdr_planeto_source_spectrum_at(&spectrum, 0, &w, &L);
+ CHK(w == 0.9 && eq_eps(L, 100.9, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 1, &w, &L);
+ CHK(w == 1 && eq_eps(L, 101, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 2, &w, &L);
+ CHK(w == 2 && eq_eps(L, 102, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 3, &w, &L);
+ CHK(w == 3 && eq_eps(L, 103, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 4, &w, &L);
+ CHK(w == 4 && eq_eps(L, 104, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 5, &w, &L);
+ CHK(w == 5 && eq_eps(L, 105, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 6, &w, &L);
+ CHK(w == 6 && eq_eps(L, 106, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 7, &w, &L);
+ CHK(w == 7 && eq_eps(L, 107, 1.e-6));
+ htrdr_planeto_source_spectrum_at(&spectrum, 8, &w, &L);
+ CHK(w == 7.456 && eq_eps(L, 107.456, 1.e-6));
+
+ htrdr_planeto_source_ref_put(source);
+}
+
+static void
+test_spectrum_fail(struct htrdr* htrdr)
+{
+ struct htrdr_planeto_source_args source_args = HTRDR_PLANETO_SOURCE_ARGS_NULL;
+ struct htrdr_planeto_source* source = NULL;
+ FILE* fp = NULL;
+ char rnrl_filename[] = "rnrl.bin";
+ double w, L;
+
+ source_args.rnrl_filename = rnrl_filename;
+ source_args.longitude = 0;
+ source_args.latitude = 0;
+ source_args.distance = 0;
+ source_args.radius = 1e8;
+ source_args.temperature = -1;
+
+ /* Wrong item size */
+ CHK(fp = fopen(rnrl_filename, "w"));
+ write_per_wlen_radiances(fp, 4096, 10, 8, 16);
+ CHK(fclose(fp) == 0);
+ CHK(htrdr_planeto_source_create(htrdr, &source_args, &source) == RES_BAD_ARG);
+
+ /* Wrong item alignment */
+ CHK(fp = fopen(rnrl_filename, "w"));
+ write_per_wlen_radiances(fp, 4096, 10, 16, 32);
+ CHK(fclose(fp) == 0);
+ CHK(htrdr_planeto_source_create(htrdr, &source_args, &source) == RES_BAD_ARG);
+
+ CHK(fp = fopen(rnrl_filename, "w"));
+ write_per_wlen_radiances(fp, 4096, 4, 16, 16);
+
+ /* Overwrite sorted items by unsorted items */
+ CHK(fseek(fp, 4096, SEEK_SET) == 0);
+ w = 10; L = 1;
+ CHK(fwrite(&w, sizeof(w), 1, fp) == 1);
+ CHK(fwrite(&L, sizeof(L), 1, fp) == 1);
+ w = 11; L = 2;
+ CHK(fwrite(&w, sizeof(w), 1, fp) == 1);
+ CHK(fwrite(&L, sizeof(L), 1, fp) == 1);
+ w = 9; L = 3;
+ CHK(fwrite(&w, sizeof(w), 1, fp) == 1);
+ CHK(fwrite(&L, sizeof(L), 1, fp) == 1);
+ w = 12; L = 4;
+ CHK(fwrite(&w, sizeof(w), 1, fp) == 1);
+ CHK(fwrite(&L, sizeof(L), 1, fp) == 1);
+ CHK(fclose(fp) == 0);
+
+ /* Unsorted items */
+ CHK(htrdr_planeto_source_create(htrdr, &source_args, &source) == RES_BAD_ARG);
+}
+
+static void
+test_spectrum_from_files(struct htrdr* htrdr, int argc, char** argv)
+{
+ struct htrdr_ran_wlen_discrete_create_args distrib_args =
+ HTRDR_RAN_WLEN_DISCRETE_CREATE_ARGS_NULL;
+ struct htrdr_ran_wlen_discrete* distrib = NULL;
+
+ struct htrdr_planeto_source_args source_args = HTRDR_PLANETO_SOURCE_ARGS_NULL;
+ struct htrdr_planeto_source_spectrum spectrum = HTRDR_PLANETO_SOURCE_SPECTRUM_NULL;
+ struct htrdr_planeto_source* source = NULL;
+ size_t i;
+
+ source_args.longitude = 0;
+ source_args.latitude = 0;
+ source_args.distance = 0;
+ source_args.radius = 1e8;
+ source_args.temperature = -1;
+
+ FOR_EACH(i, 1, argc) {
+ double range[2];
+ double lambda, pdf;
+ source_args.rnrl_filename = argv[i];
+
+ CHK(htrdr_planeto_source_create(htrdr, &source_args, &source) == RES_OK);
+ CHK(htrdr_planeto_source_does_radiance_vary_spectrally(source));
+ CHK(htrdr_planeto_source_get_spectral_range(source, range) == RES_OK);
+
+ range[0] = 250;
+ range[1] = 850;
+ CHK(htrdr_planeto_source_get_spectrum(source, range, &spectrum) == RES_OK);
+
+ printf("`%s' stores %lu entries between [%g, %g] nm\n",
+ argv[i], spectrum.size, SPLIT2(range));
+
+ distrib_args.get = htrdr_planeto_source_spectrum_at;
+ distrib_args.nwavelengths = spectrum.size;
+ distrib_args.context = &spectrum;
+ CHK(htrdr_ran_wlen_discrete_create(htrdr, &distrib_args, &distrib) == RES_OK);
+
+ lambda = htrdr_ran_wlen_discrete_sample(distrib, 0.3, 0.5, &pdf);
+ printf("lambda = %g nm; pdf = %f nm⁻¹\n", lambda, pdf);
+
+ htrdr_planeto_source_ref_put(source);
+ htrdr_ran_wlen_discrete_ref_put(distrib);
+ }
+}
+
+int
+main(int argc, char** argv)
+{
+ struct htrdr_args args = HTRDR_ARGS_DEFAULT;
+ struct htrdr* htrdr = NULL;
+ size_t memsz = 0;
+
+ args.verbose = 1;
+ htrdr_mpi_init(argc, argv);
+ CHK(htrdr_create(NULL, &args, &htrdr) == RES_OK);
+
+ memsz = MEM_ALLOCATED_SIZE(htrdr_get_allocator(htrdr));
+
+ if(argc > 1) {
+ test_spectrum_from_files(htrdr, argc, argv);
+ } else {
+ test_spectrum(htrdr);
+ test_spectrum_fail(htrdr);
+ }
+
+ CHK(MEM_ALLOCATED_SIZE(htrdr_get_allocator(htrdr)) == memsz);
+
+ htrdr_ref_put(htrdr);
+ htrdr_mpi_finalize();
+ return 0;
+}
