commit 15fc38e1a53f458034471fb9e0e68caa82a80237
parent dfbf43fc51669a2ca706b3acd5ac0322dace7053
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 15 Apr 2021 15:33:50 +0200
Merge branch 'feature_combustion' into develop
Diffstat:
28 files changed, 2884 insertions(+), 62 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -17,6 +17,7 @@
cmake_minimum_required(VERSION 3.1)
project(htrdr C)
+enable_testing()
set(HTRDR_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src)
set(HTRDR_BUILD_DIR ${CMAKE_CURRENT_BINARY_DIR})
@@ -32,6 +33,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
add_subdirectory(core)
add_subdirectory(atmosphere)
add_subdirectory(commands)
+add_subdirectory(combustion)
add_subdirectory(doc)
################################################################################
diff --git a/cmake/combustion/CMakeLists.txt b/cmake/combustion/CMakeLists.txt
@@ -0,0 +1,103 @@
+# 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-combustion)
+
+################################################################################
+# Check dependencies
+################################################################################
+find_package(AtrSTM REQUIRED)
+find_package(RCMake 0.3 REQUIRED)
+find_package(RSys 0.11 REQUIRED)
+find_package(Star3D 0.7.1 REQUIRED)
+find_package(StarSF 0.6 REQUIRED)
+find_package(StarSP 0.8 REQUIRED)
+find_package(StarVX 0.1 REQUIRED)
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
+include(rcmake)
+include(rcmake_runtime)
+
+include_directories(
+ ${AtrSTM_INCLUDE_DIR}
+ ${RSys_INCLUDE_DIR}
+ ${Star3D_INCLUDE_DIR}
+ ${StarSF_INCLUDE_DIR}
+ ${StarSP_INCLUDE_DIR}
+ ${StarVX_INCLUDE_DIR}
+ ${HTRDR_BUILD_DIR}
+ ${HTRDR_SOURCE_DIR})
+
+################################################################################
+# Configure and define targets
+################################################################################
+set(HTRDR_COMBUSTION_FILES_SRC
+ htrdr_combustion.c
+ htrdr_combustion_args.c
+ htrdr_combustion_draw_map.c
+ htrdr_combustion_compute_radiance_sw.c
+ htrdr_combustion_laser.c
+ htrdr_combustion_main.c)
+
+set(HTRDR_COMBUSTION_FILES_INC
+ htrdr_combustion.h
+ htrdr_combustion_c.h
+ htrdr_combustion_args.h
+ htrdr_combustion_laser.h)
+
+# Prepend each file in the `HTRDR_FILES_<SRC|INC>' list by `HTRDR_SOURCE_DIR'
+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
+ ${HTRDR_COMBUSTION_FILES_SRC}
+ ${HTRDR_COMBUSTION_FILES_INC})
+target_link_libraries(htrdr-combustion htrdr-core AtrSTM RSys Star3D StarSF StarSP)
+
+set_target_properties(htrdr-combustion PROPERTIES
+ DEFINE_SYMBOL HTRDR_SHARED_BUILD
+ VERSION ${VERSION}
+ SOVERSION ${VERSION_MAJOR})
+
+################################################################################
+# Add tests
+################################################################################
+if(NOT NO_TEST)
+ function(build_test _name)
+ add_executable(${_name}
+ ${HTRDR_SOURCE_DIR}/combustion/${_name}.c)
+ target_link_libraries(${_name} htrdr-combustion ${ARGN})
+ endfunction()
+
+ function(new_test _name)
+ build_test(${_name} ${ARGN})
+ add_test(${_name} ${_name})
+ endfunction()
+
+ new_test(test_htrdr_combustion_laser)
+endif()
+
+################################################################################
+# Define output & install directories
+################################################################################
+install(TARGETS htrdr-core
+ ARCHIVE DESTINATION bin
+ LIBRARY DESTINATION lib
+ RUNTIME DESTINATION bin)
+
diff --git a/cmake/commands/CMakeLists.txt b/cmake/commands/CMakeLists.txt
@@ -30,7 +30,7 @@ include_directories(${RSys_INCLUDE_DIR})
# Configure and define targets
################################################################################
add_executable(htrdr_cmd ${HTRDR_SOURCE_DIR}/commands/htrdr_cmd.c)
-target_link_libraries(htrdr_cmd htrdr-atmosphere)
+target_link_libraries(htrdr_cmd htrdr-atmosphere htrdr-combustion)
set_target_properties(htrdr_cmd PROPERTIES
OUTPUT_NAME htrdr)
@@ -40,6 +40,12 @@ target_link_libraries(htrdr_atmosphere_cmd htrdr-atmosphere)
set_target_properties(htrdr_atmosphere_cmd PROPERTIES
OUTPUT_NAME htrdr-atmosphere)
+add_executable(htrdr_combustion_cmd
+ ${HTRDR_SOURCE_DIR}/commands/htrdr_combustion_cmd.c)
+target_link_libraries(htrdr_combustion_cmd htrdr-combustion)
+set_target_properties(htrdr_combustion_cmd PROPERTIES
+ OUTPUT_NAME htrdr-combustion)
+
################################################################################
# Define output & install directories
################################################################################
diff --git a/src/atmosphere/htrdr_atmosphere.c b/src/atmosphere/htrdr_atmosphere.c
@@ -391,7 +391,7 @@ htrdr_atmosphere_create
if(!cmd->dump_volumetric_acceleration_structure) {
- struct atmosphere_pixel_format pixfmt = ATMOSPHERE_PIXEL_FORMAT_NULL;
+ struct htrdr_pixel_format pixfmt = HTRDR_PIXEL_FORMAT_NULL;
atmosphere_get_pixel_format(cmd, &pixfmt);
/* Setup the buffer layout */
@@ -459,8 +459,9 @@ error:
void
atmosphere_get_pixel_format
(const struct htrdr_atmosphere* cmd,
- struct atmosphere_pixel_format* fmt)
+ struct htrdr_pixel_format* fmt)
{
+ ASSERT(cmd && fmt);
switch(cmd->sensor.type) {
case HTRDR_SENSOR_RECTANGLE:
fmt->size = sizeof(struct atmosphere_pixel_flux);
diff --git a/src/atmosphere/htrdr_atmosphere_c.h b/src/atmosphere/htrdr_atmosphere_c.h
@@ -36,14 +36,6 @@ enum atmosphere_radiance_cpnt_flag {
| ATMOSPHERE_RADIANCE_DIFFUSE
};
-struct atmosphere_pixel_format {
- size_t size; /* In bytes */
- size_t alignment; /* Power of two, in Bytes */
-};
-#define ATMOSPHERE_PIXEL_FORMAT_NULL__ {0, 0}
-static const struct atmosphere_pixel_format ATMOSPHERE_PIXEL_FORMAT_NULL =
- ATMOSPHERE_PIXEL_FORMAT_NULL__;
-
struct atmosphere_pixel_xwave {
struct htrdr_estimate radiance; /* In W/m^2/sr */
struct htrdr_estimate radiance_temperature; /* In K */
@@ -52,7 +44,7 @@ struct atmosphere_pixel_xwave {
#define ATMOSPHERE_PIXEL_XWAVE_NULL__ { \
HTRDR_ESTIMATE_NULL__, /* Radiance */ \
HTRDR_ESTIMATE_NULL__, /* Radiance temperature */ \
- HTRDR_ACCUM_NULL /* Time */ \
+ HTRDR_ACCUM_NULL__ /* Time */ \
}
static const struct atmosphere_pixel_xwave ATMOSPHERE_PIXEL_XWAVE_NULL =
ATMOSPHERE_PIXEL_XWAVE_NULL__;
@@ -62,8 +54,8 @@ struct atmosphere_pixel_flux {
struct htrdr_accum time;
};
#define ATMOSPHERE_PIXEL_FLUX_NULL__ { \
- HTRDR_ACCUM_NULL, \
- HTRDR_ACCUM_NULL \
+ HTRDR_ACCUM_NULL__, \
+ HTRDR_ACCUM_NULL__ \
}
static const struct atmosphere_pixel_flux ATMOSPHERE_PIXEL_FLUX_NULL =
ATMOSPHERE_PIXEL_FLUX_NULL__;
@@ -78,7 +70,7 @@ struct atmosphere_pixel_image {
HTRDR_ESTIMATE_NULL__, /* X */ \
HTRDR_ESTIMATE_NULL__, /* Y */ \
HTRDR_ESTIMATE_NULL__, /* Z */ \
- HTRDR_ACCUM_NULL /* Time */ \
+ HTRDR_ACCUM_NULL__ /* Time */ \
}
static const struct atmosphere_pixel_image ATMOSPHERE_PIXEL_IMAGE_NULL =
ATMOSPHERE_PIXEL_IMAGE_NULL__;
@@ -129,7 +121,7 @@ struct htrdr_atmosphere {
extern LOCAL_SYM void
atmosphere_get_pixel_format
(const struct htrdr_atmosphere* cmd,
- struct atmosphere_pixel_format* fmt);
+ struct htrdr_pixel_format* fmt);
extern LOCAL_SYM res_T
atmosphere_draw_map
diff --git a/src/atmosphere/htrdr_atmosphere_draw_map.c b/src/atmosphere/htrdr_atmosphere_draw_map.c
@@ -165,6 +165,7 @@ draw_pixel_image
case 2: wlen = htrdr_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 */
iband = htsky_find_spectral_band(cmd->sky, wlen);
iquad = htsky_spectral_band_sample_quadrature(cmd->sky, r2, iband);
@@ -174,7 +175,6 @@ draw_pixel_image
ATMOSPHERE_RADIANCE_ALL, ray_org, ray_dir, wlen, iband, iquad);
ASSERT(weight >= 0);
- pdf *= 1.e9; /* Transform the pdf from nm^-1 to m^-1 */
weight /= pdf; /* In W/m^2/sr */
/* End the registration of the per realisation time */
@@ -250,6 +250,7 @@ draw_pixel_flux
/* Sample a wavelength */
wlen = htrdr_ran_wlen_sample(cmd->ran_wlen, 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 */
iband = htsky_find_spectral_band(cmd->sky, wlen);
@@ -370,6 +371,7 @@ draw_pixel_xwave
/* Sample a wavelength */
wlen = htrdr_ran_wlen_sample(cmd->ran_wlen, 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 */
iband = htsky_find_spectral_band(cmd->sky, wlen);
@@ -395,7 +397,7 @@ draw_pixel_xwave
time_sub(&t0, time_current(&t1), &t0);
usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
- /* Update the pixel accumulator of the current channel */
+ /* Update the pixel accumulator */
radiance.sum_weights += weight;
radiance.sum_weights_sqr += weight*weight;
radiance.nweights += 1;
@@ -540,7 +542,7 @@ dump_buffer
struct htrdr_accum* flux_acc, /* May be NULL */
FILE* stream)
{
- struct atmosphere_pixel_format pixfmt;
+ struct htrdr_pixel_format pixfmt;
struct htrdr_buffer_layout layout;
size_t x, y;
ASSERT(cmd && buf && stream);
diff --git a/src/combustion/htrdr_combustion.c b/src/combustion/htrdr_combustion.c
@@ -0,0 +1,448 @@
+/* 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 "combustion/htrdr_combustion.h"
+#include "combustion/htrdr_combustion_args.h"
+#include "combustion/htrdr_combustion_c.h"
+#include "combustion/htrdr_combustion_laser.h"
+
+#include "core/htrdr.h"
+#include "core/htrdr_camera.h"
+#include "core/htrdr_log.h"
+#include "core/htrdr_geometry.h"
+#include "core/htrdr_materials.h"
+#include "core/htrdr_rectangle.h"
+
+#include <astoria/atrstm.h>
+
+#include <star/ssf.h>
+
+#include <rsys/cstr.h>
+#include <rsys/mem_allocator.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+release_phase_functions
+ (struct htrdr_combustion* cmd,
+ struct ssf_phase* phases[])
+{
+ size_t i;
+ ASSERT(cmd);
+
+ if(!phases) return; /* Nothing to release */
+
+ FOR_EACH(i, 0, htrdr_get_threads_count(cmd->htrdr)) {
+ if(phases[i]) {
+ SSF(phase_ref_put(phases[i]));
+ }
+ }
+ MEM_RM(htrdr_get_allocator(cmd->htrdr), phases);
+}
+
+static res_T
+setup_output
+ (struct htrdr_combustion* cmd,
+ const struct htrdr_combustion_args* args)
+{
+ const char* output_name = NULL;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ if(htrdr_get_mpi_rank(cmd->htrdr) != 0) {
+ /* No output stream on non master processes */
+ cmd->output = NULL;
+ output_name = "<null>";
+
+ } else if(!args->path_output) {
+ /* Write results to standard output when no destination file is defined */
+ cmd->output = stdout;
+ output_name = "<stdout>";
+
+ } else {
+ /* Open the output stream */
+ res = htrdr_open_output_stream(cmd->htrdr, args->path_output,
+ 0/*read*/, args->force_overwriting, &cmd->output);
+ if(res != RES_OK) goto error;
+ output_name = args->path_output;
+ }
+
+ /* Setup the output name */
+ 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",
+ output_name, res_to_cstr(res));
+ goto error;
+ }
+
+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 res_T
+setup_geometry
+ (struct htrdr_combustion* cmd,
+ const struct htrdr_combustion_args* args)
+{
+ res_T res = RES_OK;
+
+ if(!args->geom.path_obj) goto exit;
+ ASSERT(args->geom.path_mats);
+
+ res = htrdr_materials_create(cmd->htrdr, args->geom.path_mats, &cmd->mats);
+ if(res != RES_OK) goto error;
+
+ res = htrdr_geometry_create
+ (cmd->htrdr, args->geom.path_obj, cmd->mats, &cmd->geom);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ if(cmd->mats) { htrdr_materials_ref_put(cmd->mats); cmd->mats = NULL; }
+ if(cmd->geom) { htrdr_geometry_ref_put(cmd->geom); cmd->geom = NULL; }
+ goto exit;
+}
+
+static res_T
+setup_camera
+ (struct htrdr_combustion* cmd,
+ const struct htrdr_combustion_args* args)
+{
+ double proj_ratio = 0;
+ ASSERT(cmd && args && args->image.definition[0] && args->image.definition[1]);
+
+ proj_ratio =
+ (double)args->image.definition[0]
+ / (double)args->image.definition[1];
+
+ return htrdr_camera_create
+ (cmd->htrdr,
+ args->camera.position,
+ args->camera.target,
+ args->camera.up,
+ proj_ratio,
+ MDEG2RAD(args->camera.fov_y),
+ &cmd->camera);
+}
+
+static res_T
+setup_laser
+ (struct htrdr_combustion* cmd,
+ const struct htrdr_combustion_args* args)
+{
+ struct htrdr_combustion_laser_create_args laser_args =
+ HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT;
+ ASSERT(cmd && args);
+ cmd->wavelength = args->wavelength;
+ laser_args.surface = args->laser;
+ laser_args.wavelength = args->wavelength;
+ laser_args.flux_density = args->laser_flux_density;
+ return htrdr_combustion_laser_create(cmd->htrdr, &laser_args, &cmd->laser);
+}
+
+static res_T
+setup_phase_functions
+ (struct htrdr_combustion* cmd,
+ const struct ssf_phase_type* phase_type,
+ struct ssf_phase** out_phases[])
+{
+ struct mem_allocator* allocator = NULL;
+ struct ssf_phase** phases = NULL;
+ size_t nthreads;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(cmd && phase_type && out_phases);
+
+ nthreads = htrdr_get_threads_count(cmd->htrdr);
+ allocator = htrdr_get_allocator(cmd->htrdr);
+
+ /* Allocate the list of per thread phase function */
+ phases = MEM_CALLOC(allocator, nthreads, sizeof(*phases));
+ if(!phases) {
+ htrdr_log_err(cmd->htrdr,
+ "Could not allocate the per thread RDG-FA phase function.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Create the per thread phase function */
+ FOR_EACH(i, 0, nthreads) {
+ res = ssf_phase_create(allocator, phase_type, phases+i);
+ if(res != RES_OK) {
+ htrdr_log_err(cmd->htrdr,
+ "Could not create the phase function for the thread %lu -- %s.\n",
+ (unsigned long)i, res_to_cstr(res));
+ goto error;
+ }
+ }
+
+exit:
+ *out_phases = phases;
+ return res;
+error:
+ release_phase_functions(cmd, phases);
+ phases = NULL;
+ goto exit;
+}
+
+static res_T
+setup_buffer
+ (struct htrdr_combustion* cmd,
+ const struct htrdr_combustion_args* args)
+{
+ struct htrdr_pixel_format pixfmt = HTRDR_PIXEL_FORMAT_NULL;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ if(cmd->dump_volumetric_acceleration_structure) goto exit;
+
+ combustion_get_pixel_format(cmd, &pixfmt);
+
+ /* Setup the 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;
+
+ /* Create the image buffer only on the master process; the image parts
+ * rendered by the others processes are gathered onto the master process */
+ 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 res_T
+setup_medium
+ (struct htrdr_combustion* cmd,
+ const struct htrdr_combustion_args* args)
+{
+ struct atrstm_args atrstm_args = ATRSTM_ARGS_DEFAULT;
+ res_T res = RES_OK;
+ ASSERT(cmd && args);
+
+ /* Setup the semi-transaprent medium arguments */
+ atrstm_args.sth_filename = args->path_tetra;
+ atrstm_args.atrtp_filename = args->path_therm_props;
+ atrstm_args.atrri_filename = args->path_refract_ids;
+ atrstm_args.cache_filename = args->path_cache;
+ atrstm_args.spectral_type = ATRSTM_SPECTRAL_SW;
+ atrstm_args.wlen_range[0] = args->wavelength;
+ atrstm_args.wlen_range[1] = args->wavelength;
+ atrstm_args.fractal_prefactor = args->fractal_prefactor;
+ atrstm_args.fractal_dimension = args->fractal_dimension;
+ atrstm_args.optical_thickness = args->optical_thickness;
+ atrstm_args.precompute_normals = args->precompute_normals;
+ atrstm_args.nthreads = args->nthreads;
+ atrstm_args.verbose = args->verbose;
+
+ switch(args->grid.type) {
+ case HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_AUTO:
+ atrstm_args.auto_grid_definition = 1;
+ atrstm_args.auto_grid_definition_hint = args->grid.definition.hint;
+ break;
+ case HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_FIXED:
+ atrstm_args.auto_grid_definition = 0;
+ atrstm_args.grid_max_definition[0] = args->grid.definition.fixed[0];
+ atrstm_args.grid_max_definition[1] = args->grid.definition.fixed[1];
+ atrstm_args.grid_max_definition[2] = args->grid.definition.fixed[2];
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+
+ /* Here we go! Create the semi-transparent medium */
+ res = atrstm_create
+ (htrdr_get_logger(cmd->htrdr),
+ htrdr_get_allocator(cmd->htrdr),
+ &atrstm_args,
+ &cmd->medium);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ if(cmd->medium) { ATRSTM(ref_put(cmd->medium)); cmd->medium = NULL; }
+ goto exit;
+}
+
+static res_T
+dump_volumetric_acceleration_structure(struct htrdr_combustion* cmd)
+{
+ struct atrstm_dump_svx_octree_args args = ATRSTM_DUMP_SVX_OCTREE_ARGS_DEFAULT;
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ /* Nothing to do on non master process */
+ if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
+
+ htrdr_log(cmd->htrdr, "Write volumetric acceleration structure to '%s'.\n",
+ str_cget(&cmd->output_name));
+
+ res = atrstm_dump_svx_octree(cmd->medium, &args, cmd->output);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static void
+combustion_release(ref_T* ref)
+{
+ struct htrdr_combustion* cmd = CONTAINER_OF(ref, struct htrdr_combustion, ref);
+ struct htrdr* htrdr = NULL;
+ ASSERT(ref);
+
+ if(cmd->geom) htrdr_geometry_ref_put(cmd->geom);
+ if(cmd->mats) htrdr_materials_ref_put(cmd->mats);
+ if(cmd->medium) ATRSTM(ref_put(cmd->medium));
+ if(cmd->camera) htrdr_camera_ref_put(cmd->camera);
+ if(cmd->laser) htrdr_combustion_laser_ref_put(cmd->laser);
+ if(cmd->buf) htrdr_buffer_ref_put(cmd->buf);
+ if(cmd->output && cmd->output != stdout) CHK(fclose(cmd->output) == 0);
+ release_phase_functions(cmd, cmd->rdgfa_phase_functions);
+ release_phase_functions(cmd, cmd->hg_phase_functions);
+ str_release(&cmd->output_name);
+
+ htrdr = cmd->htrdr;
+ MEM_RM(htrdr_get_allocator(htrdr), cmd);
+ htrdr_ref_put(htrdr);
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+htrdr_combustion_create
+ (struct htrdr* htrdr,
+ const struct htrdr_combustion_args* args,
+ struct htrdr_combustion** out_cmd)
+{
+ struct htrdr_combustion* cmd = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && args && out_cmd);
+
+ 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;
+ goto error;
+ }
+ ref_init(&cmd->ref);
+ str_init(htrdr_get_allocator(htrdr), &cmd->output_name);
+
+ /* Get the ownership on the htrdr structure */
+ htrdr_ref_get(htrdr);
+ cmd->htrdr = htrdr;
+
+ cmd->spp = args->image.spp;
+ cmd->dump_volumetric_acceleration_structure =
+ args->dump_volumetric_acceleration_structure;
+
+ res = setup_output(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_geometry(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_camera(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_laser(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_phase_functions(cmd, &ssf_phase_rdgfa, &cmd->rdgfa_phase_functions);
+ if(res != RES_OK) goto error;
+ res = setup_phase_functions(cmd, &ssf_phase_hg, &cmd->hg_phase_functions);
+ if(res != RES_OK) goto error;
+ res = setup_buffer(cmd, args);
+ if(res != RES_OK) goto error;
+ res = setup_medium(cmd, args);
+ if(res != RES_OK) goto error;
+
+exit:
+ *out_cmd = cmd;
+ return res;
+error:
+ if(cmd) {
+ htrdr_combustion_ref_put(cmd);
+ cmd = NULL;
+ }
+ goto exit;
+}
+
+void
+htrdr_combustion_ref_get(struct htrdr_combustion* cmd)
+{
+ ASSERT(cmd);
+ ref_get(&cmd->ref);
+}
+
+void
+htrdr_combustion_ref_put(struct htrdr_combustion* cmd)
+{
+ ASSERT(cmd);
+ ref_put(&cmd->ref, combustion_release);
+}
+
+res_T
+htrdr_combustion_run(struct htrdr_combustion* cmd)
+{
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ if(cmd->dump_volumetric_acceleration_structure) {
+ res = dump_volumetric_acceleration_structure(cmd);
+ if(res != RES_OK) goto error;
+ } else {
+ res = combustion_draw_map(cmd);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+void
+combustion_get_pixel_format
+ (const struct htrdr_combustion* cmd,
+ struct htrdr_pixel_format* fmt)
+{
+ ASSERT(cmd && fmt);
+ (void)cmd;
+ fmt->size = sizeof(struct combustion_pixel);
+ fmt->alignment = ALIGNOF(struct combustion_pixel);
+}
diff --git a/src/combustion/htrdr_combustion.h b/src/combustion/htrdr_combustion.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_COMBUSTION_H
+#define HTRDR_COMBUSTION_H
+
+#include "core/htrdr.h"
+#include <rsys/rsys.h>
+
+struct htrdr;
+struct htrdr_combustion;
+struct htrdr_combustion_args;
+
+BEGIN_DECLS
+
+HTRDR_API res_T
+htrdr_combustion_create
+ (struct htrdr* htrdr,
+ const struct htrdr_combustion_args* args,
+ struct htrdr_combustion** cmd);
+
+HTRDR_API void
+htrdr_combustion_ref_get
+ (struct htrdr_combustion* cmd);
+
+HTRDR_API void
+htrdr_combustion_ref_put
+ (struct htrdr_combustion* cmd);
+
+HTRDR_API res_T
+htrdr_combustion_run
+ (struct htrdr_combustion* cmd);
+
+HTRDR_API int
+htrdr_combustion_main
+ (int argc,
+ char** argv);
+
+END_DECLS
+
+#endif /* HTRDR_COMBUSTION_H */
diff --git a/src/combustion/htrdr_combustion_args.c b/src/combustion/htrdr_combustion_args.c
@@ -0,0 +1,307 @@
+/* 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 "combustion/htrdr_combustion_args.h"
+
+#include <rsys/cstr.h>
+
+#include <getopt.h>
+#include <string.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+print_help(const char* cmd)
+{
+ ASSERT(cmd);
+ printf(
+"Usage: %s [<options>] -m TETRAHEDRA -p THERMOPROPS -r REFRACT_IDS\n",
+ cmd);
+ printf(
+"Render a monochromatic image within a sooting flame described according\n"
+"to the RDG-FA theory and lightened by a laser source.\n\n");
+
+ printf(
+" -C <camera> define the rendering point of view. Refer to the\n"
+" man page for the list of camera options.\n");
+ printf(
+" -D FLUX_DENSITY\n"
+" flux density of the laser in W/m^2. By default the\n"
+" flux density is %g W/m^2.\n",
+ HTRDR_COMBUSTION_ARGS_DEFAULT.laser_flux_density);
+ printf(
+" -d dump volumetric acceleration structures to OUTPUT\n"
+" and exit.\n");
+ printf(
+" -F <fractal-coefs>\n"
+" value of the fractal prefactor and fractal dimension\n"
+" to use in the RDG-FA model. Refer to the man page\n"
+" for the syntax of the <fractal-coefs> option. Default\n"
+" fractal prefactor is %g and the default fractal\n"
+" dimension is %g.\n",
+ HTRDR_COMBUSTION_ARGS_DEFAULT.fractal_prefactor,
+ HTRDR_COMBUSTION_ARGS_DEFAULT.fractal_dimension);
+ printf(
+" -f overwrite the OUTPUT file if it already exists.\n");
+ printf(
+" -g <geometry> define the combustion chamber geometry. Refer to the\n"
+" man page for the list of geometry options.\n");
+ printf(
+" -h display this help and exit.\n");
+ printf(
+" -i <image> define the image to compute. Refer to the man\n"
+" page for the list of image options.\n");
+ printf(
+" -l <laser> define the geometry of the laser sheet. Refer to the\n"
+" man page for the list of laser options.\n");
+ printf(
+" -m TETRAHEDRA path toward the volumetric mesh.\n");
+ printf(
+" -N precompute the tetrahedra normals.\n");
+ printf(
+" -O CACHE path of the cache file used to store/restore the\n"
+" volumetric data. By default do not use any cache.\n");
+ printf(
+" -o OUTPUT file where data are written. If not defined, data are\n"
+" written to standard output.\n");
+ printf(
+" -p THERMOPROPS path toward the thermodynamic properties.\n");
+ printf(
+" -r REFRACT_ID path toward the per wavelength refractive\n"
+" indices.\n");
+ printf(
+" -T THRESHOLD optical thickness used as threshold during the octree\n"
+" building. By default its value is %g.\n",
+ HTRDR_COMBUSTION_ARGS_DEFAULT.optical_thickness);
+ printf(
+" -t NTHREADS hint on the number of threads to use. By default use\n"
+" as many threads as CPU cores.\n");
+ printf(
+" -V <HINT|X,Y,Z>\n"
+" definition of the volumetric acceleration grids along\n"
+" the 3 axis. By default it is computed automatically\n"
+" with a hint on the expected definition set to %u.\n",
+ HTRDR_COMBUSTION_ARGS_DEFAULT.grid.definition.hint);
+ printf(
+" -v make the command verbose.\n");
+ printf(
+" -w WAVELENGTH wavelength definition of the laser in nanometer.\n"
+" By default its value is %g.\n",
+ HTRDR_COMBUSTION_ARGS_DEFAULT.wavelength);
+
+ printf("\n");
+
+ htrdr_fprint_copyright(cmd, stdout);
+ htrdr_fprint_license(cmd, stdout);
+}
+
+static res_T
+parse_grid_definition
+ (struct htrdr_combustion_args_grid_definition* grid_def,
+ const char* str)
+{
+ unsigned def[3];
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(grid_def && str);
+
+ res = cstr_to_list_uint(str, ',', def, &len, 3);
+ if(res == RES_OK && len == 2) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ fprintf(stderr, "Invalid grid definition `%s'.\n", str);
+ goto error;
+ }
+
+ if(len == 1) {
+ if(!def[0]) {
+ fprintf(stderr, "Invalid null grid definition %u.\n", def[0]);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ grid_def->type = HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_AUTO;
+ grid_def->definition.hint = def[0];
+
+ } else {
+ if(!def[0] || !def[1] || !def[2]) {
+ fprintf(stderr,
+ "Invalid null grid definition [%u, %u, %u].\n", SPLIT3(def));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ grid_def->type = HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_FIXED;
+ grid_def->definition.fixed[0] = def[0];
+ grid_def->definition.fixed[1] = def[1];
+ grid_def->definition.fixed[2] = def[2];
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_fractal_parameters(const char* str, void* ptr)
+{
+ char buf[128];
+ struct htrdr_combustion_args* args = ptr;
+ char* key;
+ char* val;
+ char* ctx;
+ res_T res = RES_OK;
+ ASSERT(ptr && str);
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr,
+ "Could not duplicate the fractal option string `%s'.\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &ctx);
+ val = strtok_r(NULL, "", &ctx);
+
+ if(!strcmp(key, "prefactor")) {
+ res = cstr_to_double(val, &args->fractal_prefactor);
+ if(res != RES_OK) goto error;
+ } else if(!strcmp(key, "dimension")) {
+ res = cstr_to_double(val, &args->fractal_dimension);
+ if(res != RES_OK) goto error;
+ } else {
+ fprintf(stderr, "Invalid fractal parameter `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_combustion_args_init
+ (struct htrdr_combustion_args* args,
+ int argc,
+ char** argv)
+{
+ int opt;
+ res_T res = RES_OK;
+ ASSERT(args && argc && argv);
+
+ *args = HTRDR_COMBUSTION_ARGS_DEFAULT;
+
+ while((opt = getopt(argc, argv, "C:D:dF:fg:hi:l:m:NO:o:p:r:T:t:V:vw:")) != -1) {
+ switch(opt) {
+ case 'C':
+ res = htrdr_args_camera_parse(&args->camera, optarg);
+ break;
+ case 'D':
+ res = cstr_to_double(optarg, &args->laser_flux_density);
+ if(res == RES_OK && args->laser_flux_density <= 0) res = RES_BAD_ARG;
+ break;
+ case 'd':
+ args->dump_volumetric_acceleration_structure = 1;
+ break;
+ case 'F':
+ res = cstr_parse_list(optarg, ':', parse_fractal_parameters, args);
+ break;
+ case 'f': args->force_overwriting = 1; break;
+ case 'g':
+ res = htrdr_args_geometry_parse(&args->geom, optarg);
+ break;
+ case 'h':
+ print_help(argv[0]);
+ htrdr_combustion_args_release(args);
+ args->quit = 1;
+ goto exit;
+ case 'i':
+ res = htrdr_args_image_parse(&args->image, optarg);
+ break;
+ case 'l':
+ res = htrdr_args_rectangle_parse(&args->laser, optarg);
+ break;
+ case 'm': args->path_tetra = optarg; break;
+ case 'N': args->precompute_normals = 1; break;
+ case 'O': args->path_cache = optarg; break;
+ case 'o': args->path_output = optarg; break;
+ case 'p': args->path_therm_props = optarg; break;
+ case 'r': args->path_refract_ids = optarg; 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': /* Submit an hint on the number of threads to use */
+ res = cstr_to_uint(optarg, &args->nthreads);
+ if(res == RES_OK && !args->nthreads) res = RES_BAD_ARG;
+ break;
+ case 'V':
+ res = parse_grid_definition(&args->grid, optarg);
+ break;
+ case 'v': args->verbose = 1; break;
+ case 'w':
+ res = cstr_to_double(optarg, &args->wavelength);
+ break;
+ default: res = RES_BAD_ARG; break;
+ }
+ if(res != RES_OK) {
+ if(optarg) {
+ fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n",
+ argv[0], optarg, opt);
+ }
+ goto error;
+ }
+ }
+
+ if(!args->path_tetra) {
+ fprintf(stderr, "Missing the volumetric mesh -- option '-m'\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(!args->path_therm_props) {
+ fprintf(stderr, "Missing the thermodynamic properties -- option '-p'\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ if(!args->path_refract_ids) {
+ fprintf(stderr, "Missing the refractive indices -- option '-r'\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ htrdr_combustion_args_release(args);
+ goto exit;
+}
+
+void
+htrdr_combustion_args_release(struct htrdr_combustion_args* args)
+{
+ ASSERT(args);
+ htrdr_args_geometry_free(&args->geom);
+ *args = HTRDR_COMBUSTION_ARGS_DEFAULT;
+}
+
diff --git a/src/combustion/htrdr_combustion_args.h b/src/combustion/htrdr_combustion_args.h
@@ -0,0 +1,127 @@
+/* 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_COMBUSTION_ARGS_H
+#define HTRDR_COMBUSTION_ARGS_H
+
+#include "core/htrdr_args.h"
+#include <rsys/rsys.h>
+
+#include <limits.h> /* UINT_MAX support */
+
+enum htrdr_combustion_args_grid_definition_type {
+ HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_FIXED,
+ HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_AUTO,
+ HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_TYPES_COUNT__
+};
+
+struct htrdr_combustion_args_grid_definition {
+ union {
+ unsigned hint; /* Hint on the grid definition to eval */
+ unsigned fixed[3]; /* Fixed grid definition along the 3 axis */
+ } definition;
+ enum htrdr_combustion_args_grid_definition_type type;
+};
+#define HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_DEFAULT__ { \
+ {256}, \
+ HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_AUTO \
+}
+static const struct htrdr_combustion_args_grid_definition
+HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_DEFAULT =
+ HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_DEFAULT__;
+
+struct htrdr_combustion_args {
+ struct htrdr_args_geometry geom; /* Combustion chamber geometry */
+
+ const char* path_tetra; /* Volumetric mesh of the medium */
+ const char* path_therm_props; /* Termodynamic properties of the medium */
+ const char* path_refract_ids; /* Refractive indices in the medium */
+
+ const char* path_cache; /* Path of the file to store/restore cached data */
+ const char* path_output; /* Name of the output file */
+
+ struct htrdr_args_camera camera; /* Pinhole Camera */
+
+ struct htrdr_args_rectangle laser; /* Laser surface emission */
+ double wavelength; /* Wavelength of the laser in nanometer */
+ double laser_flux_density; /* In W/m^2 */
+
+ struct htrdr_args_image image; /* Output Image */
+
+ /* RDG-FA parameters */
+ double fractal_prefactor;
+ double fractal_dimension;
+
+ struct htrdr_combustion_args_grid_definition grid;
+
+ double optical_thickness; /* Threshold used during octree building */
+
+ /* Miscellaneous parameters */
+ unsigned nthreads; /* Hint on the number of threads to use */
+ int precompute_normals; /* Pre-compute the tetrahedra normals */
+ int force_overwriting;
+ int dump_volumetric_acceleration_structure;
+ int verbose; /* Verbosity level */
+ int quit; /* Stop the command */
+};
+
+#define HTRDR_COMBUSTION_ARGS_DEFAULT__ { \
+ HTRDR_ARGS_GEOMETRY_NULL__, \
+ \
+ NULL, /* Tetra path */ \
+ NULL, /* Therm props path */ \
+ NULL, /* Refractive ids path */ \
+ \
+ NULL, /* Cache path */ \
+ NULL, /* Output path */ \
+ \
+ HTRDR_ARGS_CAMERA_DEFAULT__, /* Pinhole camera */ \
+ \
+ HTRDR_ARGS_RECTANGLE_DEFAULT__, /* Laser surface emission */ \
+ 532, /* Wavelength in nanometer */ \
+ 1, /* Flux density */ \
+ \
+ HTRDR_ARGS_IMAGE_DEFAULT__, /* Image */ \
+ \
+ 1.30, /* Gyration radius prefactor */ \
+ 1.80, /* Fractal dimension */ \
+ \
+ HTRDR_COMBUSTION_ARGS_GRID_DEFINITION_DEFAULT__, \
+ \
+ 1, /* Optical thickness */ \
+ \
+ UINT_MAX, /* #threads */ \
+ 0, /* Precompute normals */ \
+ 0, /* Force overwriting */ \
+ 0, /* dump volumetric acceleration structure */ \
+ 0, /* Verbose flag */ \
+ 0 /* Stop the command */ \
+}
+static const struct htrdr_combustion_args HTRDR_COMBUSTION_ARGS_DEFAULT =
+ HTRDR_COMBUSTION_ARGS_DEFAULT__;
+
+extern LOCAL_SYM res_T
+htrdr_combustion_args_init
+ (struct htrdr_combustion_args* args,
+ int argc,
+ char** argv);
+
+extern LOCAL_SYM void
+htrdr_combustion_args_release
+ (struct htrdr_combustion_args* args);
+
+#endif /* HTRDR_COMBUSTION_ARGS_H */
diff --git a/src/combustion/htrdr_combustion_c.h b/src/combustion/htrdr_combustion_c.h
@@ -0,0 +1,93 @@
+/* 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_COMBUSTION_C_H
+#define HTRDR_COMBUSTION_C_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 atrstm;
+struct htrdr;
+struct htrdr_camera;
+struct htrdr_combustion_laser;
+struct htrdr_geometry;
+struct htrdr_materials;
+struct htrdr_rectangle;
+struct ssf_phase;
+struct ssp_rng;
+
+struct combustion_pixel {
+ struct htrdr_estimate radiance; /* In W/m^2/sr */
+ struct htrdr_accum time; /* In microseconds */
+};
+#define COMBUSTION_PIXEL_NULL__ { \
+ HTRDR_ESTIMATE_NULL__, /* Radiance */ \
+ HTRDR_ACCUM_NULL__, /* Time */ \
+}
+static const struct combustion_pixel COMBUSTION_PIXEL_NULL =
+ COMBUSTION_PIXEL_NULL__;
+
+struct htrdr_combustion {
+ struct htrdr_geometry* geom; /* Combustion chamber geometry */
+ struct htrdr_materials* mats; /* Materials of the combustion chamber */
+ struct atrstm* medium; /* Semi transparent medium */
+
+ struct htrdr_camera* camera; /* Pinhole camera */
+ struct htrdr_combustion_laser* laser; /* Laser sheet */
+ double wavelength; /* Wavelength of the laser in nanometer */
+
+ struct ssf_phase** rdgfa_phase_functions; /* Per thread RDG-FA phase func */
+ struct ssf_phase** hg_phase_functions; /* Per thread Henyey-Greenstein func */
+
+ struct htrdr_buffer_layout buf_layout;
+ struct htrdr_buffer* buf; /* NULL on non master processes */
+ size_t spp; /* #samples per pixel */
+
+ FILE* output; /* Output stream */
+ struct str output_name; /* Name of the output stream */
+ int dump_volumetric_acceleration_structure;
+
+
+ ref_T ref;
+ struct htrdr* htrdr;
+};
+
+extern LOCAL_SYM void
+combustion_get_pixel_format
+ (const struct htrdr_combustion* cmd,
+ struct htrdr_pixel_format* fmt);
+
+extern LOCAL_SYM res_T
+combustion_draw_map
+ (struct htrdr_combustion* cmd);
+
+/* Return the shortwave radiance in W/m^2/sr */
+extern LOCAL_SYM double
+combustion_compute_radiance_sw
+ (struct htrdr_combustion* cmd,
+ const size_t ithread,
+ struct ssp_rng* rng,
+ const double pos_in[3],
+ const double dir_in[3]);
+
+#endif /* HTRDR_COMBUSTION_C_H */
diff --git a/src/combustion/htrdr_combustion_compute_radiance_sw.c b/src/combustion/htrdr_combustion_compute_radiance_sw.c
@@ -0,0 +1,712 @@
+/* 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 "combustion/htrdr_combustion_c.h"
+#include "combustion/htrdr_combustion_laser.h"
+
+#include "core/htrdr.h"
+
+#include <astoria/atrstm.h>
+
+#include <star/ssf.h>
+#include <star/ssp.h>
+
+#include <rsys/double2.h>
+#include <rsys/double3.h>
+#include <rsys/double4.h>
+#include <rsys/dynamic_array.h>
+
+/* Define a position along the ray into the semi-transparent medium */
+struct position {
+ double distance; /* Distance up to the position */
+ struct suvm_primitive prim; /* Volumetric primitive of the position */
+ double bcoords[4]; /* Local coordinate of the position into `prim' */
+};
+#define POSITION_NULL__ { \
+ DBL_MAX, /* Distance */ \
+ SUVM_PRIMITIVE_NULL__, /* Primitive */ \
+ {0, 0, 0, 0} /* Barycentric coordinates */ \
+}
+static const struct position POSITION_NULL = POSITION_NULL__;
+
+/* Syntactic sugar to check if the position is valid */
+#define POSITION_NONE(Pos) ((Pos)->distance >= DBL_MAX)
+
+/* Common position but preferentially sampled within a limited range. Its
+ * associated ksi variable defines the correction of the weight due to the
+ * normalization of the sampling pdf, and the recursivity associated with the
+ * null-collision technique. */
+struct position_limited {
+ struct position position;
+ double ksi;
+};
+#define POSITION_LIMITED_NULL__ {POSITION_NULL__, 1}
+static const struct position_limited POSITION_LIMITED_NULL =
+ POSITION_LIMITED_NULL__;
+
+struct sample_position_context {
+ /* Input data */
+ struct ssp_rng* rng; /* Random Number Generator */
+ struct atrstm* medium; /* Semi-transparent medium */
+ double wavelength; /* Wavelength to handel in nanometer */
+ enum atrstm_radcoef radcoef; /* Radiative coef used to sample a position */
+ double Ts; /* Sampled optical thickness */
+
+ /* Output data */
+ struct position position;
+};
+#define SAMPLE_POSITION_CONTEXT_NULL__ { \
+ NULL, /* RNG */ \
+ NULL, /* Medium */ \
+ 0, /* Wavelength */ \
+ ATRSTM_RADCOEFS_COUNT__, /* Radiative coefficient */ \
+ 0, /* Optical thickness */ \
+ POSITION_NULL__ \
+}
+static const struct sample_position_context SAMPLE_POSITION_CONTEXT_NULL =
+ SAMPLE_POSITION_CONTEXT_NULL__;
+
+struct sample_scattering_limited_context {
+ /* Input data */
+ struct ssp_rng* rng; /* Random number generator to use */
+ struct atrstm* medium; /* Semi transparent medium */
+ double wavelength; /* Wavelength to handle in nanometer */
+
+ /* Local parameters update during ray traversal */
+ double ks_2hat; /* Smallest scattering upper-field over the ray range */
+ double Tmax; /* Scattering optical thickness computed using ks_2hat */
+ double Ume; /* Normalization of the pdf for the sampled optical thickness */
+ double sampled_vox_collision_dst; /* Scattering path length */
+
+ /* Output data */
+ struct position_limited position_limited;
+};
+#define SAMPLE_SCATTERING_LIMITED_CONTEXT_NULL__ { \
+ NULL, /* RNG */ \
+ NULL, /* Medium */ \
+ -1, /* Wavelength */ \
+ -1, /* ks_2hat */ \
+ -1, /* Tau max */ \
+ -1, /* Ume */ \
+ -1, /* Sampled collision dst */ \
+ POSITION_LIMITED_NULL__ \
+}
+static const struct sample_scattering_limited_context
+SAMPLE_SCATTERING_LIMITED_CONTEXT_NULL =
+ SAMPLE_SCATTERING_LIMITED_CONTEXT_NULL__;
+
+/*******************************************************************************
+ * Sample a position along a ray into the inhomogeneous medium for a given
+ * radiative coefficient
+ ******************************************************************************/
+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)
+{
+ atrstm_radcoefs_svx_T radcoefs_svx;
+ struct atrstm_fetch_radcoefs_args fetch_raw_args;
+ struct atrstm_fetch_radcoefs_svx_voxel_args fetch_svx_args;
+ struct sample_position_context* ctx = context;
+ double k_min = 0;
+ double k_max = 0;
+ double traversal_dst = 0;
+ int pursue_traversal = 1;
+ ASSERT(hit && org && dir && range && ctx);
+ (void)range;
+
+ /* Fetch the K<min|max> of the current traversed voxel */
+ fetch_svx_args.voxel = hit->voxel;
+ fetch_svx_args.radcoefs_mask = BIT(ctx->radcoef);
+ fetch_svx_args.components_mask = ATRSTM_CPNTS_MASK_ALL;
+ fetch_svx_args.operations_mask = ATRSTM_SVX_OP_FLAG_MIN | ATRSTM_SVX_OP_FLAG_MAX;
+ ATRSTM(fetch_radcoefs_svx_voxel(ctx->medium, &fetch_svx_args, radcoefs_svx));
+
+ k_min = radcoefs_svx[ctx->radcoef][ATRSTM_SVX_OP_MIN];
+ k_max = radcoefs_svx[ctx->radcoef][ATRSTM_SVX_OP_MAX];
+
+ /* Setup the constants of the 'fetch' function for the current voxel */
+ fetch_raw_args.wavelength = ctx->wavelength;
+ fetch_raw_args.radcoefs_mask = BIT(ctx->radcoef);
+ fetch_raw_args.components_mask = ATRSTM_CPNTS_MASK_ALL;
+ fetch_raw_args.k_min[ctx->radcoef] = k_min;
+ fetch_raw_args.k_max[ctx->radcoef] = k_max;
+
+ /* Initialised the already traversed distance to the distance from which the
+ * current ray enters into the current voxel */
+ traversal_dst = hit->distance[0];
+
+ for(;;) {
+ /* Compute the optical thickness of the current leaf */
+ const double vox_dst = hit->distance[1] - traversal_dst;
+ 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_collision_dst = ctx->Ts / k_max;
+ atrstm_radcoefs_T radcoefs;
+ double x[3];
+ double k;
+ double r;
+
+ /* Compute the traversed distance up to the challenged collision */
+ traversal_dst += vox_collision_dst;
+
+ /* Compute the world space position where a collision may occur */
+ x[0] = org[0] + traversal_dst * dir[0];
+ x[1] = org[1] + traversal_dst * dir[1];
+ x[2] = org[2] + traversal_dst * dir[2];
+
+ /* Fetch the radiative coefficient at the collision position */
+ ATRSTM(at(ctx->medium, x, &fetch_raw_args.prim, fetch_raw_args.bcoords));
+ if(SUVM_PRIMITIVE_NONE(&fetch_raw_args.prim)) {
+ k = 0;
+ } else {
+ ATRSTM(fetch_radcoefs(ctx->medium, &fetch_raw_args, radcoefs));
+ k = radcoefs[ctx->radcoef];
+ }
+
+ r = ssp_rng_canonical(ctx->rng);
+
+ if(r > k/k_max) { /* Null collision */
+ ctx->Ts = ssp_ran_exp(ctx->rng, 1);
+ } else { /* Real collision */
+ /* Setup output data */
+ ctx->position.distance = traversal_dst;
+ ctx->position.prim = fetch_raw_args.prim;
+ d4_set(ctx->position.bcoords, fetch_raw_args.bcoords);
+
+ /* Stop ray traversal */
+ pursue_traversal = 0;
+ break;
+ }
+ }
+ }
+ return pursue_traversal;
+}
+
+static void
+sample_position
+ (struct htrdr_combustion* cmd,
+ struct ssp_rng* rng,
+ const enum atrstm_radcoef radcoef,
+ const double pos[3],
+ const double dir[3],
+ const double range[2],
+ struct position* position)
+{
+ struct atrstm_trace_ray_args args = ATRSTM_TRACE_RAY_ARGS_DEFAULT;
+ struct sample_position_context sample_pos_ctx = SAMPLE_POSITION_CONTEXT_NULL;
+ struct svx_hit svx_hit = SVX_HIT_NULL;
+ double wlen = 0;
+ ASSERT(cmd && rng && pos && dir && range);
+
+ wlen = htrdr_combustion_laser_get_wavelength(cmd->laser);
+
+ /* Sample an optical thickness */
+ sample_pos_ctx.Ts = ssp_ran_exp(rng, 1);
+
+ /* Setup the arguments of the function invoked per voxel (i.e. the filter
+ * function) */
+ sample_pos_ctx.rng = rng;
+ sample_pos_ctx.medium = cmd->medium;
+ sample_pos_ctx.wavelength = wlen;
+ sample_pos_ctx.radcoef = radcoef;
+
+ /* Setup ray tracing arguments */
+ d3_set(args.ray_org, pos);
+ d3_set(args.ray_dir, dir);
+ d2_set(args.ray_range, range);
+ args.filter = sample_position_hit_filter;
+ args.context = &sample_pos_ctx;
+
+ /* Trace the ray into the heterogeneous medium */
+ ATRSTM(trace_ray(cmd->medium, &args, &svx_hit));
+
+ if(SVX_HIT_NONE(&svx_hit)) {
+ /* No collision with the medium was found */
+ *position = POSITION_NULL;
+ } else {
+ /* A collision occurs into the medium */
+ *position = sample_pos_ctx.position;
+ }
+}
+
+/*******************************************************************************
+ * Preferentially sample a scattering position into an inhomogeneous medium
+ * according to a limited range
+ ******************************************************************************/
+/* Find the constant Ks max (named ks_2hat) along the traced ray */
+static int
+sample_scattering_limited_find_ks_2hat
+ (const struct svx_hit* hit,
+ const double org[3],
+ const double dir[3],
+ const double range[2],
+ void* context)
+{
+ struct sample_scattering_limited_context* ctx = context;
+ struct atrstm_fetch_radcoefs_svx_voxel_args fetch_svx_args;
+ atrstm_radcoefs_svx_T radcoefs_svx;
+ ASSERT(hit && org && dir && range && context);
+ (void)org, (void)dir;
+
+ /* In all situations, initialise ks_2hat with the ks_max of the root node */
+ if(hit->voxel.depth == 0) {
+ fetch_svx_args.voxel = hit->voxel;
+ fetch_svx_args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks;
+ fetch_svx_args.components_mask = ATRSTM_CPNTS_MASK_ALL;
+ fetch_svx_args.operations_mask = ATRSTM_SVX_OP_FLAG_MAX;
+ ATRSTM(fetch_radcoefs_svx_voxel(ctx->medium, &fetch_svx_args, radcoefs_svx));
+ ctx->ks_2hat = radcoefs_svx[ATRSTM_RADCOEF_Ks][ATRSTM_SVX_OP_MAX];
+
+ /* Update ks_2hat with the ks_max of the current descending node until the ray
+ * range was no more included into this node */
+ } else if(hit->distance[0] <= range[0] && range[1] <= hit->distance[1]) {
+ fetch_svx_args.voxel = hit->voxel;
+ fetch_svx_args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks;
+ fetch_svx_args.components_mask = ATRSTM_CPNTS_MASK_ALL;
+ fetch_svx_args.operations_mask = ATRSTM_SVX_OP_FLAG_MAX;
+ ATRSTM(fetch_radcoefs_svx_voxel(ctx->medium, &fetch_svx_args, radcoefs_svx));
+ ctx->ks_2hat = radcoefs_svx[ATRSTM_RADCOEF_Ks][ATRSTM_SVX_OP_MAX];
+ }
+
+ /* Do not stop here: keep diving up to the leaves */
+ return 0;
+}
+
+static int
+sample_scattering_limited_hit_filter
+ (const struct svx_hit* hit,
+ const double org[3],
+ const double dir[3],
+ const double range[2],
+ void* context)
+{
+ atrstm_radcoefs_svx_T radcoefs_svx;
+ struct atrstm_fetch_radcoefs_args fetch_raw_args;
+ struct atrstm_fetch_radcoefs_svx_voxel_args fetch_svx_args;
+ struct sample_scattering_limited_context* ctx = context;
+ double ks_min = 0;
+ double ks_max = 0;
+ double traversal_dst = 0;
+ int pursue_traversal = 1;
+ ASSERT(hit && org && dir && range && ctx);
+ (void)range;
+
+ /* Fetch the Ks<min|max> of the current traversed voxel */
+ fetch_svx_args.voxel = hit->voxel;
+ fetch_svx_args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks;
+ fetch_svx_args.components_mask = ATRSTM_CPNTS_MASK_ALL;
+ fetch_svx_args.operations_mask = ATRSTM_SVX_OP_FLAG_MIN | ATRSTM_SVX_OP_FLAG_MAX;
+ ATRSTM(fetch_radcoefs_svx_voxel(ctx->medium, &fetch_svx_args, radcoefs_svx));
+
+ ks_min = radcoefs_svx[ATRSTM_RADCOEF_Ks][ATRSTM_SVX_OP_MIN];
+ ks_max = radcoefs_svx[ATRSTM_RADCOEF_Ks][ATRSTM_SVX_OP_MAX];
+ ASSERT(ks_max <= ctx->ks_2hat);
+
+ /* Setup the constants of the 'fetch' function for the current voxel */
+ fetch_raw_args.wavelength = ctx->wavelength;
+ fetch_raw_args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Ks;
+ fetch_raw_args.components_mask = ATRSTM_CPNTS_MASK_ALL;
+ fetch_raw_args.k_min[ATRSTM_RADCOEF_Ks] = ks_min;
+ fetch_raw_args.k_max[ATRSTM_RADCOEF_Ks] = ks_max;
+
+ /* Initialised the already traversed distance to the distance from which the
+ * current ray enters into the current voxel */
+ traversal_dst = hit->distance[0];
+
+ /* First traversed leaf */
+ if(ctx->sampled_vox_collision_dst < 0) {
+ ctx->Tmax = (range[1] - range[0]) * ctx->ks_2hat;
+ ctx->Ume = 1 - exp(-ctx->Tmax);
+ }
+
+ for(;;) {
+ atrstm_radcoefs_T radcoefs;
+ double vox_dst;
+ double tau;
+ double ks;
+ double r;
+
+ /* Compute the remaining distance to traverse in the current voxel */
+ vox_dst = hit->distance[1] - traversal_dst;
+
+ /* A collision distance was not already sampled */
+ if(ctx->sampled_vox_collision_dst < 0) {
+ r = ssp_rng_canonical(ctx->rng);
+ tau = -log(1.0-r*(1.0-exp(-ctx->Tmax)));
+ ctx->sampled_vox_collision_dst = tau / ctx->ks_2hat;
+
+ /* Update the ksi output data */
+ ctx->position_limited.ksi *= ctx->Ume;
+ }
+
+ /* Compute the traversed distance up to the challenged collision */
+ traversal_dst = traversal_dst + ctx->sampled_vox_collision_dst;
+
+ /* The collision to challenge lies behind the current voxel */
+ if(traversal_dst > hit->distance[1]) {
+ ctx->sampled_vox_collision_dst -= vox_dst;
+ pursue_traversal = 1;
+ break;
+ }
+ ASSERT(traversal_dst >= hit->distance[0]);
+
+ r = ssp_rng_canonical(ctx->rng);
+ if(r >= ks_max / ctx->ks_2hat) { /* Null collision */
+ ctx->sampled_vox_collision_dst = -1;
+ } else { /* Collision with ks_max */
+ double x[3];
+
+ /* Compute the world space position where a collision may occur */
+ x[0] = org[0] + traversal_dst * dir[0];
+ x[1] = org[1] + traversal_dst * dir[1];
+ x[2] = org[2] + traversal_dst * dir[2];
+
+ /* Fetch the radiative coefficient at the collision position */
+ ATRSTM(at(ctx->medium, x, &fetch_raw_args.prim, fetch_raw_args.bcoords));
+ if(SUVM_PRIMITIVE_NONE(&fetch_raw_args.prim)) {
+ ks = 0;
+ } else {
+ ATRSTM(fetch_radcoefs(ctx->medium, &fetch_raw_args, radcoefs));
+ ks = radcoefs[ATRSTM_RADCOEF_Ks];
+ }
+
+ r = ssp_rng_canonical(ctx->rng);
+
+ if(r >= ks/ks_max) { /* Null collision */
+ ctx->sampled_vox_collision_dst = -1;
+ } else { /* Real collision */
+ /* Setup output data */
+ ctx->position_limited.position.distance = traversal_dst;
+ ctx->position_limited.position.prim = fetch_raw_args.prim;
+ d4_set(ctx->position_limited.position.bcoords, fetch_raw_args.bcoords);
+
+ /* Stop ray traversal */
+ pursue_traversal = 0;
+ break;
+ }
+ }
+ }
+ return pursue_traversal;
+}
+
+static void
+sample_scattering_limited
+ (struct htrdr_combustion* cmd,
+ struct ssp_rng* rng,
+ const double pos[3],
+ const double dir[3],
+ const double range[2],
+ struct position_limited* position)
+{
+ struct atrstm_trace_ray_args args = ATRSTM_TRACE_RAY_ARGS_DEFAULT;
+ struct sample_scattering_limited_context sample_scattering_limited_ctx =
+ SAMPLE_SCATTERING_LIMITED_CONTEXT_NULL;
+ struct svx_hit svx_hit = SVX_HIT_NULL;
+ double wlen = 0;
+ ASSERT(cmd && rng && pos && dir && position);
+
+ wlen = htrdr_combustion_laser_get_wavelength(cmd->laser);
+
+ /* Setup the trace ray context */
+ sample_scattering_limited_ctx.rng = rng;
+ sample_scattering_limited_ctx.medium = cmd->medium;
+ sample_scattering_limited_ctx.wavelength = wlen;
+ sample_scattering_limited_ctx.ks_2hat = 0;
+ sample_scattering_limited_ctx.Tmax = -1;
+ sample_scattering_limited_ctx.Ume = -1;
+ sample_scattering_limited_ctx.sampled_vox_collision_dst = -1;
+ sample_scattering_limited_ctx.position_limited = POSITION_LIMITED_NULL;
+
+ /* Setup the input arguments for the ray tracing into the medium */
+ d3_set(args.ray_org, pos);
+ d3_set(args.ray_dir, dir);
+ d2_set(args.ray_range, range);
+ args.challenge = sample_scattering_limited_find_ks_2hat;
+ args.filter = sample_scattering_limited_hit_filter;
+ args.context = &sample_scattering_limited_ctx;
+
+ /* Trace the ray into the medium to compute the transmissivity */
+ ATRSTM(trace_ray(cmd->medium, &args, &svx_hit));
+
+ if(SVX_HIT_NONE(&svx_hit)) {
+ /* No scattering event was found */
+ *position = POSITION_LIMITED_NULL;
+ } else {
+ /* A scattering event occurs into the medium */
+ *position = sample_scattering_limited_ctx.position_limited;
+ }
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static double
+transmissivity
+ (struct htrdr_combustion* cmd,
+ struct ssp_rng* rng,
+ const enum atrstm_radcoef radcoef,
+ const double pos[3],
+ const double dir[3],
+ const double range[2])
+{
+ struct position position = POSITION_NULL;
+
+ /* Degenerated range => no attenuation along dir */
+ if(range[1] <= range[0]) return 1;
+
+ sample_position(cmd, rng, radcoef, pos, dir, range, &position);
+
+ if(POSITION_NONE(&position)) {
+ return 1; /* No collision with the medium */
+ } else {
+ return 0; /* A collision occurs */
+ }
+}
+
+/* This function computes the contribution of the in-laser scattering for the
+ * current scattering position of the reverse optical path 'pos' and current
+ * direction 'dir'. One contribution has to be computed for each scattering
+ * position in order to compute the value of the Monte-Carlo weight for the
+ * optical path (weight of the statistical realization). */
+static double
+laser_once_scattered
+ (struct htrdr_combustion* cmd,
+ const size_t ithread,
+ struct ssp_rng* rng,
+ const double pos[3],
+ const double dir[3])
+{
+ /* RDG-FA phase function */
+ struct atrstm_rdgfa rdgfa_param = ATRSTM_RDGFA_NULL;
+ struct atrstm_fetch_rdgfa_args fetch_rdgfa_args =
+ ATRSTM_FETCH_RDGFA_ARGS_DEFAULT;
+ struct ssf_phase_rdgfa_setup_args setup_rdgfa_args =
+ SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT;
+ struct ssf_phase* phase = NULL;
+
+ /* Scattering position into the laser sheet */
+ struct position_limited sc_sample = POSITION_LIMITED_NULL;
+ double xsc[3] = {0, 0, 0}; /* World space position */
+
+ /* The transmissivities to compute */
+ double Tr_ext_pos_xin = 0;
+ double Tr_ext_xsc_lse = 0;
+ double Tr_abs_xin_xsc = 0;
+
+ /* Laser data */
+ double laser_hit_dst[2] = {0, 0};
+ double laser_flux_density = 0; /* In W/m^2 */
+ double wlen = 0; /* In nm */
+
+ /* Miscellaneous varbale */
+ double wi[3] = {0, 0, 0};
+ double wo[3] = {0, 0, 0};
+ double range[2] = {0, DBL_MAX};
+ double R = 0;
+
+ /* Radiance to compute in W/m^2/sr */
+ double L = 0;
+
+ ASSERT(cmd && rng && pos && dir);
+
+ /* Fetch the laser wavelength */
+ wlen = htrdr_combustion_laser_get_wavelength(cmd->laser);
+
+ /* Find the intersections along dir with the laser sheet */
+ htrdr_combustion_laser_trace_ray(cmd->laser, pos, dir, range, laser_hit_dst);
+
+ /* No intersection with the laser sheet => no laser contribution */
+ if(HTRDR_COMBUSTION_LASER_HIT_NONE(laser_hit_dst)) return 0;
+
+ /* Compute the transmissivity from 'pos' to 'xin' */
+ range[0] = 0;
+ range[1] = laser_hit_dst[0];
+ Tr_ext_pos_xin = transmissivity(cmd, rng, ATRSTM_RADCOEF_Kext, pos, dir, range);
+ if(Tr_ext_pos_xin == 0) return 0; /* no laser contribution */
+
+ /* Sample the scattering position into the laser sheet */
+ range[0] = laser_hit_dst[0];
+ range[1] = laser_hit_dst[1];
+ sample_scattering_limited(cmd, rng, pos, dir, range, &sc_sample);
+ if(POSITION_NONE(&sc_sample.position)) return 0; /* No laser contribution */
+ ASSERT(laser_hit_dst[0] <= sc_sample.position.distance);
+ ASSERT(laser_hit_dst[1] >= sc_sample.position.distance);
+
+ /* Compute the transmissivity from 'xin' to 'xsc' */
+ range[0] = laser_hit_dst[0]; /* <=> xin */
+ range[1] = sc_sample.position.distance; /* <=> xsc */
+ Tr_abs_xin_xsc = transmissivity(cmd, rng, ATRSTM_RADCOEF_Ka, pos, dir, range);
+ if(Tr_abs_xin_xsc == 0) return 0; /* No laser contribution */
+
+ /* Compute the scattering position */
+ xsc[0] = pos[0] + dir[0] * sc_sample.position.distance;
+ xsc[1] = pos[1] + dir[1] * sc_sample.position.distance;
+ xsc[2] = pos[2] + dir[2] * sc_sample.position.distance;
+
+ /* Retrieve the direction toward the laser surface emission */
+ htrdr_combustion_laser_get_direction(cmd->laser, wi);
+ d3_minus(wi, wi);
+
+ /* Compute the transmissivity from xsc to the laser surface emission */
+ range[0] = 0;
+ range[1] = htrdr_combustion_laser_compute_surface_plane_distance
+ (cmd->laser, xsc);
+ ASSERT(range[1] >= 0);
+ Tr_ext_xsc_lse = transmissivity(cmd, rng, ATRSTM_RADCOEF_Kext, xsc, wi, range);
+ if(Tr_ext_xsc_lse == 0) return 0; /* No laser contribution */
+
+ /* Retrieve the RDG-FA phase function parameters from the semi transparent
+ * medium */
+ fetch_rdgfa_args.wavelength = wlen;
+ fetch_rdgfa_args.prim = sc_sample.position.prim;
+ fetch_rdgfa_args.bcoords[0] = sc_sample.position.bcoords[0];
+ fetch_rdgfa_args.bcoords[1] = sc_sample.position.bcoords[1];
+ fetch_rdgfa_args.bcoords[2] = sc_sample.position.bcoords[2];
+ fetch_rdgfa_args.bcoords[3] = sc_sample.position.bcoords[3];
+ ATRSTM(fetch_rdgfa(cmd->medium, &fetch_rdgfa_args, &rdgfa_param));
+
+ /* Setup the RDG-FA phase function */
+ phase = cmd->rdgfa_phase_functions[ithread];
+ setup_rdgfa_args.wavelength = rdgfa_param.wavelength;
+ setup_rdgfa_args.fractal_dimension = rdgfa_param.fractal_dimension;
+ setup_rdgfa_args.gyration_radius = rdgfa_param.gyration_radius;
+ SSF(phase_rdgfa_setup(phase, &setup_rdgfa_args));
+
+ /* Evaluate the phase function at the scattering position */
+ d3_minus(wo, dir); /* Ensure SSF convention */
+ R = ssf_phase_eval(phase, wo, wi); /* In sr^-1 */
+
+ laser_flux_density = htrdr_combustion_laser_get_flux_density(cmd->laser);
+
+ L = Tr_ext_pos_xin
+ * Tr_abs_xin_xsc
+ * sc_sample.ksi * R
+ * Tr_ext_xsc_lse
+ * laser_flux_density;
+ return L;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+extern LOCAL_SYM double
+combustion_compute_radiance_sw
+ (struct htrdr_combustion* cmd,
+ const size_t ithread,
+ struct ssp_rng* rng,
+ const double pos_in[3],
+ const double dir_in[3])
+{
+ /* RDG-FA phase function */
+ struct atrstm_rdgfa rdgfa_param = ATRSTM_RDGFA_NULL;
+ struct atrstm_fetch_rdgfa_args fetch_rdgfa_args =
+ ATRSTM_FETCH_RDGFA_ARGS_DEFAULT;
+ struct ssf_phase_rdgfa_setup_args setup_rdgfa_args =
+ SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT;
+ struct ssf_phase* phase = NULL;
+
+ /* Transmissivity between the probe position (i.e. 'pos_in') and the current
+ * scattering position over the reverse scattering path */
+ double Tr_abs = 1;
+
+ /* Monte carlo weight of the simulated optical path */
+ double weight = 0;
+
+ /* Miscellaneous variables */
+ double pos[3] = {0, 0, 0};
+ double dir[3] = {0, 0, 0};
+ double range[2] = {0, DBL_MAX};
+ double wlen = 0;
+ ASSERT(cmd && rng && pos_in && dir_in);
+
+ d3_set(pos, pos_in);
+ d3_set(dir, dir_in);
+ d2(range, 0, FLT_MAX);
+
+ wlen = htrdr_combustion_laser_get_wavelength(cmd->laser);
+
+ Tr_abs = 1;
+ weight = 0;
+
+ for(;;) {
+ struct position scattering = POSITION_NULL;
+ double laser_sheet_emissivity = 0; /* In W/m^2/sr */
+ double Tr_abs_pos_xsc = 0;
+ double wo[3];
+ double wi[3];
+
+ /* Handle the laser contribution */
+ laser_sheet_emissivity = laser_once_scattered(cmd, ithread, rng, pos, dir);
+ weight += Tr_abs * laser_sheet_emissivity;
+
+ /* Sample a scattering position */
+ sample_position(cmd, rng, ATRSTM_RADCOEF_Ks, pos, dir, range, &scattering);
+ if(POSITION_NONE(&scattering)) break;
+
+ /* Compute the absorption transmissivity */
+ range[0] = 0;
+ range[1] = scattering.distance;
+ Tr_abs_pos_xsc = transmissivity(cmd, rng, ATRSTM_RADCOEF_Ka, pos, dir, range);
+ if(Tr_abs_pos_xsc == 0) break;
+
+ /* Update the overall absorption transmissivity of the optical path */
+ Tr_abs *= Tr_abs_pos_xsc;
+
+ /* Update the position of the optical path */
+ pos[0] = pos[0] + dir[0] * scattering.distance;
+ pos[1] = pos[1] + dir[1] * scattering.distance;
+ pos[2] = pos[2] + dir[2] * scattering.distance;
+
+ /* Retrieve the RDG-FA phase function parameters */
+ fetch_rdgfa_args.wavelength = wlen;
+ fetch_rdgfa_args.prim = scattering.prim;
+ fetch_rdgfa_args.bcoords[0] = scattering.bcoords[0];
+ fetch_rdgfa_args.bcoords[1] = scattering.bcoords[1];
+ fetch_rdgfa_args.bcoords[2] = scattering.bcoords[2];
+ fetch_rdgfa_args.bcoords[3] = scattering.bcoords[3];
+ ATRSTM(fetch_rdgfa(cmd->medium, &fetch_rdgfa_args, &rdgfa_param));
+
+ /* Setup the RDG-FA phase function corresponding to the scattering event */
+ phase = cmd->rdgfa_phase_functions[ithread];
+ setup_rdgfa_args.wavelength = rdgfa_param.wavelength;
+ setup_rdgfa_args.fractal_dimension = rdgfa_param.fractal_dimension;
+ setup_rdgfa_args.gyration_radius = rdgfa_param.gyration_radius;
+ SSF(phase_rdgfa_setup(phase, &setup_rdgfa_args));
+
+ /* Sample a new optical path direction from the phase function */
+ d3_minus(wo, dir); /* Ensure SSF convention */
+ ssf_phase_sample(phase, rng, wo, wi, NULL);
+
+ /* Update the optical path direction */
+ dir[0] = wi[0];
+ dir[1] = wi[1];
+ dir[2] = wi[2];
+ }
+ return weight;
+}
diff --git a/src/combustion/htrdr_combustion_draw_map.c b/src/combustion/htrdr_combustion_draw_map.c
@@ -0,0 +1,186 @@
+/* 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 "combustion/htrdr_combustion_c.h"
+
+#include "core/htrdr_accum.h"
+#include "core/htrdr_camera.h"
+#include "core/htrdr_draw_map.h"
+#include "core/htrdr_log.h"
+
+#include <star/ssp.h>
+
+#include <rsys/clock_time.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+draw_pixel
+ (struct htrdr* htrdr,
+ const struct htrdr_draw_pixel_args* args,
+ void* data)
+{
+ struct htrdr_accum radiance = HTRDR_ACCUM_NULL;
+ struct htrdr_accum time = HTRDR_ACCUM_NULL;
+ struct htrdr_combustion* cmd = NULL;
+ struct combustion_pixel* pixel = NULL;
+ size_t isamp;
+ ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data);
+ (void)htrdr;
+
+ cmd = args->context;
+ pixel = data;
+
+ FOR_EACH(isamp, 0, args->spp) {
+ struct time t0, t1;
+ double pix_samp[2];
+ double ray_org[3];
+ double ray_dir[3];
+ double weight;
+ double usec;
+
+ /* Begin the registration of the time spent in the realisation */
+ time_current(&t0);
+
+ /* Sample a position into the pixel, in the normalized image plane */
+ pix_samp[0] = (double)args->pixel_coord[0] + ssp_rng_canonical(args->rng);
+ pix_samp[1] = (double)args->pixel_coord[1] + ssp_rng_canonical(args->rng);
+ pix_samp[0] *= args->pixel_normalized_size[0];
+ pix_samp[1] *= args->pixel_normalized_size[1];
+
+ /* Generate a ray starting from the pinhole camera and passing through the
+ * pixel sample */
+ htrdr_camera_ray(cmd->camera, pix_samp, ray_org, ray_dir);
+
+ /* Backward trace the path */
+ weight = combustion_compute_radiance_sw(cmd, args->ithread, args->rng,
+ ray_org, ray_dir);
+
+ /* End the registration of the per realisation time */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update the pixel accumulator */
+ radiance.sum_weights += weight;
+ radiance.sum_weights_sqr += weight*weight;
+ radiance.nweights += 1;
+
+ /* Update the pixel accumulator of per realisation time */
+ time.sum_weights += usec;
+ time.sum_weights_sqr += usec*usec;
+ time.nweights += 1;
+ }
+
+ /* Compute the estimation of the pixel radiance */
+ htrdr_accum_get_estimation(&radiance, &pixel->radiance);
+
+ /* Save the per realisation time */
+ pixel->time = time;
+}
+
+static void
+dump_pixel
+ (const struct combustion_pixel* pix,
+ struct htrdr_accum* time_acc, /* May be NULL */
+ FILE* stream)
+{
+ struct htrdr_estimate time = HTRDR_ESTIMATE_NULL;
+ ASSERT(pix && stream && pix->time.nweights);
+
+ htrdr_accum_get_estimation(&pix->time, &time);
+
+ fprintf(stream, "%g %g 0 0 0 0 %g %g\n",
+ pix->radiance.E, pix->radiance.SE, time.E,time.SE);
+
+ if(time_acc) {
+ time_acc->sum_weights += pix->time.sum_weights;
+ time_acc->sum_weights_sqr += pix->time.sum_weights_sqr;
+ time_acc->nweights += pix->time.nweights;
+ }
+}
+
+static res_T
+dump_buffer
+ (struct htrdr_combustion* cmd,
+ struct htrdr_buffer* buf,
+ struct htrdr_accum* time_acc, /* May be NULL */
+ FILE* stream)
+{
+ struct htrdr_pixel_format pixfmt;
+ struct htrdr_buffer_layout layout;
+ size_t x, y;
+ ASSERT(cmd && buf && stream);
+
+ combustion_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);
+
+ if(time_acc) *time_acc = HTRDR_ACCUM_NULL;
+
+ FOR_EACH(y, 0, layout.height) {
+ FOR_EACH(x, 0, layout.width) {
+ const struct combustion_pixel* pix = htrdr_buffer_at(buf, x, y);
+ ASSERT(IS_ALIGNED(pix, pixfmt.alignment));
+ dump_pixel(pix, time_acc, stream);
+ }
+ fprintf(stream, "\n");
+ }
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+combustion_draw_map(struct htrdr_combustion* cmd)
+{
+ struct htrdr_draw_map_args args = HTRDR_DRAW_MAP_ARGS_NULL;
+ struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL;
+ struct htrdr_estimate path_time = HTRDR_ESTIMATE_NULL;
+ res_T res = RES_OK;
+ ASSERT(cmd);
+
+ /* Setup the draw map input arguments */
+ args.draw_pixel = draw_pixel;
+ args.buffer_layout = cmd->buf_layout;
+ args.spp = cmd->spp;
+ args.context = cmd;
+
+ res = htrdr_draw_map(cmd->htrdr, &args, cmd->buf);
+ if(res != RES_OK) goto error;
+
+ /* No more to do on non master processes */
+ if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
+
+ /* Write buffer to output */
+ res = dump_buffer(cmd, cmd->buf, &path_time_acc, cmd->output);
+ if(res != RES_OK) goto error;
+
+ htrdr_accum_get_estimation(&path_time_acc, &path_time);
+ htrdr_log(cmd->htrdr,
+ "Time per radiative path (in micro seconds): %g +/- %g\n",
+ path_time.E, path_time.SE);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
diff --git a/src/combustion/htrdr_combustion_laser.c b/src/combustion/htrdr_combustion_laser.c
@@ -0,0 +1,266 @@
+/* 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 "combustion/htrdr_combustion_laser.h"
+
+#include "core/htrdr.h"
+#include "core/htrdr_log.h"
+#include "core/htrdr_rectangle.h"
+
+#include <rsys/cstr.h>
+#include <rsys/double33.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+struct htrdr_combustion_laser {
+ struct htrdr_rectangle* surface; /* Surface emission */
+ double world2local[12];
+ double low_ls[3], upp_ls[3]; /* Local space AABB */
+ double flux_density; /* In W/m^2 */
+ double wavelength; /* In nm */
+ ref_T ref;
+ struct htrdr* htrdr;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+laser_release(ref_T* ref)
+{
+ struct htrdr_combustion_laser* laser;
+ struct htrdr* htrdr;
+ ASSERT(ref);
+ laser = CONTAINER_OF(ref, struct htrdr_combustion_laser, ref);
+ if(laser->surface) htrdr_rectangle_ref_put(laser->surface);
+ htrdr = laser->htrdr;
+ MEM_RM(htrdr_get_allocator(htrdr), laser);
+ htrdr_ref_put(htrdr);
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_combustion_laser_create
+ (struct htrdr* htrdr,
+ const struct htrdr_combustion_laser_create_args* args,
+ struct htrdr_combustion_laser** out_laser)
+{
+ struct htrdr_combustion_laser* laser = NULL;
+ res_T res = RES_OK;
+ ASSERT(htrdr && args && out_laser);
+
+ if(args->flux_density <= 0) {
+ htrdr_log_err(htrdr, "Invalid flux density %g W.m^2\n", args->flux_density);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(args->wavelength <= 0) {
+ htrdr_log_err(htrdr, "Invalid wavelength %g nm\n", args->wavelength);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ laser = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*laser));
+ if(!laser) {
+ htrdr_log_err(htrdr, "Could not allocate the laser data structure.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&laser->ref);
+ htrdr_ref_get(htrdr);
+ laser->htrdr = htrdr;
+ laser->flux_density = args->flux_density;
+ laser->wavelength = args->wavelength;
+
+ res = htrdr_rectangle_create
+ (laser->htrdr,
+ args->surface.size,
+ args->surface.position,
+ args->surface.target,
+ args->surface.up,
+ &laser->surface);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr, "Could not create the laser surface emission -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ htrdr_rectangle_get_transform_inverse(laser->surface, laser->world2local);
+
+ /* Define the laser sheet AABB in local space */
+ laser->upp_ls[0] = args->surface.size[0] * 0.5;
+ laser->upp_ls[1] = args->surface.size[1] * 0.5;
+ laser->upp_ls[2] = 0;
+ laser->low_ls[0] = -laser->upp_ls[0];
+ laser->low_ls[1] = -laser->upp_ls[1];
+ laser->upp_ls[2] = INF;
+
+exit:
+ *out_laser = laser;
+ return res;
+error:
+ if(laser) { htrdr_combustion_laser_ref_put(laser); laser = NULL; }
+ goto exit;
+}
+
+void
+htrdr_combustion_laser_ref_get(struct htrdr_combustion_laser* laser)
+{
+ ASSERT(laser);
+ ref_get(&laser->ref);
+}
+
+void
+htrdr_combustion_laser_ref_put(struct htrdr_combustion_laser* laser)
+{
+ ASSERT(laser);
+ ref_put(&laser->ref, laser_release);
+}
+
+void
+htrdr_combustion_laser_trace_ray
+ (struct htrdr_combustion_laser* laser,
+ const double pos[3],
+ const double dir[3],
+ const double range[2],
+ double distance[2])
+{
+ double pos_ls[3]; /* Position in local space */
+ double dir_ls[3]; /* Direction in local space */
+ double t_min;
+ double t_max;
+ int i;
+ ASSERT(laser && pos && dir && range && distance);
+
+ /* Transform the ray in local space */
+ d33_muld3(pos_ls, laser->world2local, pos);
+ d33_muld3(dir_ls, laser->world2local, dir);
+ d3_add(pos_ls, pos_ls, laser->world2local+9);
+
+ /* Reset the distance */
+ distance[0] = INF;
+ distance[1] = INF;
+
+ /* Initialise the range */
+ t_min = range[0];
+ t_max = range[1];
+
+ /* TODO one can optimise the Ray/AABB intersection test along the Z axis
+ * whose AABB interval is in [0, inf) */
+ FOR_EACH(i, 0, 3) {
+ const double rcp_dir_ls = 1.0/dir_ls[i];
+ double t0 = (laser->low_ls[i] - pos_ls[i]) * rcp_dir_ls;
+ double t1 = (laser->upp_ls[i] - pos_ls[i]) * rcp_dir_ls;
+ if(t0 > t1) SWAP(double, t0, t1);
+ t_min = MMAX(t_min, t0);
+ t_max = MMIN(t_max, t1);
+ if(t_min > t_max) return; /* No intersection */
+ }
+
+ /* Save the hit distance */
+ distance[0] = t_min;
+ distance[1] = t_max;
+}
+
+void
+htrdr_combustion_laser_get_mesh
+ (const struct htrdr_combustion_laser* laser,
+ const double extend,
+ struct htrdr_combustion_laser_mesh* mesh)
+{
+ double transform[12];
+ double low[3];
+ double upp[3];
+ ASSERT(laser && extend > 0 && mesh);
+
+ htrdr_rectangle_get_transform(laser->surface, transform);
+
+ /* Define the laser sheet vertices in local space */
+ low[0] = laser->low_ls[0];
+ low[1] = laser->low_ls[1];
+ low[2] = laser->low_ls[2];
+ upp[0] = laser->upp_ls[0];
+ upp[1] = laser->upp_ls[1];
+ upp[2] = laser->low_ls[2] + extend;
+
+ /* Transform the laser sheet vertices in world space */
+ d3(mesh->vertices + 0*3, low[0], low[1], low[2]);
+ d3(mesh->vertices + 1*3, upp[0], low[1], low[2]);
+ d3(mesh->vertices + 2*3, low[0], upp[1], low[2]);
+ d3(mesh->vertices + 3*3, upp[0], upp[1], low[2]);
+ d3(mesh->vertices + 4*3, low[0], low[1], upp[2]);
+ d3(mesh->vertices + 5*3, upp[0], low[1], upp[2]);
+ d3(mesh->vertices + 6*3, low[0], upp[1], upp[2]);
+ d3(mesh->vertices + 7*3, upp[0], upp[1], upp[2]);
+ mesh->nvertices = 8;
+
+ /* Define the laser sheet triangles */
+ mesh->triangles[0] = 0; mesh->triangles[1] = 2; mesh->triangles[2] = 1;
+ mesh->triangles[3] = 1; mesh->triangles[4] = 2; mesh->triangles[5] = 3;
+ mesh->triangles[6] = 0; mesh->triangles[7] = 4; mesh->triangles[8] = 2;
+ mesh->triangles[9] = 2; mesh->triangles[10] = 4; mesh->triangles[11] = 6;
+ mesh->triangles[12] = 1; mesh->triangles[13] = 3; mesh->triangles[14] = 5;
+ mesh->triangles[15] = 5; mesh->triangles[16] = 3; mesh->triangles[17] = 7;
+ mesh->triangles[18] = 1; mesh->triangles[19] = 5; mesh->triangles[20] = 0;
+ mesh->triangles[21] = 0; mesh->triangles[22] = 5; mesh->triangles[23] = 4;
+ mesh->triangles[24] = 3; mesh->triangles[25] = 2; mesh->triangles[26] = 7;
+ mesh->triangles[27] = 7; mesh->triangles[28] = 2; mesh->triangles[29] = 6;
+ mesh->ntriangles = 10;
+}
+
+void
+htrdr_combustion_laser_get_direction
+ (const struct htrdr_combustion_laser* laser,
+ double dir[3])
+{
+ ASSERT(laser && dir);
+ htrdr_rectangle_get_normal(laser->surface, dir);
+}
+
+double
+htrdr_combustion_laser_get_flux_density
+ (const struct htrdr_combustion_laser* laser)
+{
+ ASSERT(laser);
+ return laser->flux_density;
+}
+
+double
+htrdr_combustion_laser_get_wavelength
+ (const struct htrdr_combustion_laser* laser)
+{
+ ASSERT(laser);
+ return laser->wavelength;
+}
+
+double
+htrdr_combustion_laser_compute_surface_plane_distance
+ (const struct htrdr_combustion_laser* laser,
+ const double pos[3])
+{
+ double center[3];
+ double normal[3];
+ double vec[3];
+ ASSERT(laser && pos);
+ htrdr_rectangle_get_normal(laser->surface, normal);
+ htrdr_rectangle_get_center(laser->surface, center);
+ d3_sub(vec, pos, center);
+ return d3_dot(normal, vec);
+}
diff --git a/src/combustion/htrdr_combustion_laser.h b/src/combustion/htrdr_combustion_laser.h
@@ -0,0 +1,106 @@
+/* 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_COMBUSTION_LASER_H
+#define HTRDR_COMBUSTION_LASER_H
+
+#include "core/htrdr_args.h"
+
+#include <rsys/rsys.h>
+
+/* Monochromatic laser */
+struct htrdr_combustion_laser_create_args {
+ struct htrdr_args_rectangle surface; /* Surface emission */
+ double wavelength; /* In nanometers */
+ double flux_density; /* In W/m^2 */
+};
+#define HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT__ { \
+ HTRDR_ARGS_RECTANGLE_DEFAULT__, \
+ -1, /* Wavelength */ \
+ -1, /* Flux density */ \
+}
+static const struct htrdr_combustion_laser_create_args
+HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT =
+ HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT__;
+
+struct htrdr_combustion_laser_mesh {
+ double vertices[8/*#vertices*/*3/*#coords per vertex*/];
+ /* Triangle are clock wise ordered from the inside of the laser sheet */
+ unsigned triangles[10/*#triangles*/*3/*#vertices per triangle*/];
+ unsigned nvertices;
+ unsigned ntriangles;
+};
+
+/* Syntactic sugar to check if a laser sheet hit is valid */
+#define HTRDR_COMBUSTION_LASER_HIT_NONE(Hit) ((Hit)[0] >= DBL_MAX)
+
+/* Forward declaration */
+struct htrdr;
+struct htrdr_combustion_laser;
+
+BEGIN_DECLS
+
+HTRDR_API 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
+htrdr_combustion_laser_ref_get
+ (struct htrdr_combustion_laser* laser);
+
+HTRDR_API void
+htrdr_combustion_laser_ref_put
+ (struct htrdr_combustion_laser* laser);
+
+HTRDR_API void
+htrdr_combustion_laser_trace_ray
+ (struct htrdr_combustion_laser* laser,
+ const double pos[3],
+ const double dir[3],
+ const double range[2],
+ double distance[2]);
+
+HTRDR_API void
+htrdr_combustion_laser_get_mesh
+ (const struct htrdr_combustion_laser* laser,
+ /* Max distance of the laser mesh along its infinite dimension */
+ const double extend,
+ struct htrdr_combustion_laser_mesh* mesh);
+
+HTRDR_API void
+htrdr_combustion_laser_get_direction
+ (const struct htrdr_combustion_laser* laser,
+ double dir[3]); /* Normalized */
+
+HTRDR_API double /* In W.m^2 */
+htrdr_combustion_laser_get_flux_density
+ (const struct htrdr_combustion_laser* laser);
+
+HTRDR_API double /* In nm */
+htrdr_combustion_laser_get_wavelength
+ (const struct htrdr_combustion_laser* laser);
+
+HTRDR_API double
+htrdr_combustion_laser_compute_surface_plane_distance
+ (const struct htrdr_combustion_laser* laser,
+ const double pos[3]);
+
+END_DECLS
+
+#endif /* HTRDR_COMBUSTION_LASER_H */
diff --git a/src/combustion/htrdr_combustion_main.c b/src/combustion/htrdr_combustion_main.c
@@ -0,0 +1,87 @@
+/* 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 "combustion/htrdr_combustion.h"
+#include "combustion/htrdr_combustion_args.h"
+
+#include "core/htrdr.h"
+#include "core/htrdr_args.h"
+#include "core/htrdr_log.h"
+
+#include <rsys/mem_allocator.h>
+
+int
+htrdr_combustion_main(int argc, char** argv)
+{
+ char cmd_name[] = "htrdr-combustion";
+ struct htrdr_args htrdr_args = HTRDR_ARGS_DEFAULT;
+ 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);
+ 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_combustion_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.dump_volumetric_acceleration_structure
+ && htrdr_get_mpi_rank(htrdr) != 0) {
+ goto exit; /* Nothing to do except for the master process */
+ }
+
+ res = htrdr_combustion_create(htrdr, &cmd_args, &cmd);
+ if(res != RES_OK) goto error;
+
+ res = htrdr_combustion_run(cmd);
+ if(res != RES_OK) goto error;
+
+exit:
+ htrdr_combustion_args_release(&cmd_args);
+ if(is_mpi_init) htrdr_mpi_finalize();
+ if(htrdr) htrdr_ref_put(htrdr);
+ if(cmd) htrdr_combustion_ref_put(cmd);
+
+ /* Check memory leaks */
+ memsz_end = MEM_ALLOCATED_SIZE(&mem_default_allocator);
+ 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/combustion/test_htrdr_combustion_laser.c b/src/combustion/test_htrdr_combustion_laser.c
@@ -0,0 +1,140 @@
+/* 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 "combustion/htrdr_combustion_laser.h"
+
+#include <rsys/double2.h>
+#include <rsys/double3.h>
+
+static void
+dump_obj(const struct htrdr_combustion_laser_mesh* mesh, FILE* stream)
+{
+ unsigned i;
+ ASSERT(mesh && stream);
+
+ FOR_EACH(i, 0, mesh->nvertices) {
+ fprintf(stream, "v %g %g %g\n",
+ mesh->vertices[i*3+0],
+ mesh->vertices[i*3+1],
+ mesh->vertices[i*3+2]);
+ }
+ FOR_EACH(i, 0, mesh->ntriangles) {
+ fprintf(stream, "f %u %u %u\n",
+ mesh->triangles[i*3+0]+1,
+ mesh->triangles[i*3+1]+1,
+ mesh->triangles[i*3+2]+1);
+ }
+}
+
+int
+main(int argc, char** argv)
+{
+ struct htrdr_args args = HTRDR_ARGS_DEFAULT;
+ struct htrdr_combustion_laser_mesh mesh;
+ struct htrdr_combustion_laser_create_args laser_args =
+ HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT;
+ struct htrdr* htrdr = NULL;
+ struct htrdr_combustion_laser* laser = NULL;
+ double org[3];
+ double dir[3];
+ double range[2];
+ double hit_range[2];
+ double t[2];
+ double pt[3];
+ double x[3], y[3], z[3];
+ double plane0[4];
+ double plane1[4];
+ FILE* fp = NULL;
+
+ args.verbose = 1;
+ htrdr_mpi_init(argc, argv);
+ CHK(htrdr_create(NULL, &args, &htrdr) == RES_OK);
+
+ /* Setup the laser sheet */
+ d3(laser_args.surface.position, 0, 0, 0);
+ d3(laser_args.surface.target, 10, 10, 0);
+ d3(laser_args.surface.up, 0, 1, 0);
+ d2(laser_args.surface.size, 100, 50);
+ laser_args.wavelength = 300;
+ laser_args.flux_density = 1;
+ CHK(htrdr_combustion_laser_create(htrdr, &laser_args, &laser) == RES_OK);
+ htrdr_combustion_laser_get_mesh(laser, 100/*arbitrary extend*/, &mesh);
+
+ /* Write the laser geometry */
+ CHK(fp = fopen("laser.obj", "w"));
+ dump_obj(&mesh, fp);
+ fclose(fp);
+
+ /* Compute the frame of the surface emission */
+ d3_sub(z, laser_args.surface.target, laser_args.surface.position);
+ d3_cross(x, z, laser_args.surface.up);
+ d3_cross(y, x, z);
+ CHK(d3_normalize(y, y) != 0);
+
+ /* Compute the bottom plane equation of the laser sheet */
+ pt[0] = laser_args.surface.position[0] - y[0]*laser_args.surface.size[1]*0.5;
+ pt[1] = laser_args.surface.position[1] - y[1]*laser_args.surface.size[1]*0.5;
+ pt[2] = laser_args.surface.position[2] - y[2]*laser_args.surface.size[1]*0.5;
+ plane0[0] = y[0];
+ plane0[1] = y[1];
+ plane0[2] = y[2];
+ plane0[3] = -d3_dot(y, pt);
+
+ /* Compute the top plane equation of the laser sheet */
+ pt[0] = laser_args.surface.position[0] + y[0]*laser_args.surface.size[1]*0.5;
+ pt[1] = laser_args.surface.position[1] + y[1]*laser_args.surface.size[1]*0.5;
+ pt[2] = laser_args.surface.position[2] + y[2]*laser_args.surface.size[1]*0.5;
+ plane1[0] = y[0];
+ plane1[1] = y[1];
+ plane1[2] = y[2];
+ plane1[3] = -d3_dot(y, pt);
+
+ /* Trace a ray that misses the laser sheet */
+ d3(org, 50, 0, 0);
+ d3(dir, 0, -1, 0);
+ d2(range, 0, INF);
+ htrdr_combustion_laser_trace_ray(laser, org, dir, range, hit_range);
+ CHK(hit_range[0] > DBL_MAX);
+ CHK(hit_range[1] > DBL_MAX);
+
+ /* Trace a ray that intersects both bottom and top laser planes */
+ d3(dir, 0, 1, 0);
+ htrdr_combustion_laser_trace_ray(laser, org, dir, range, hit_range);
+ CHK(hit_range[0] < hit_range[1]);
+
+ /* Compute the intersection of the ray with the bottom/top laser planes */
+ t[0] = (-d3_dot(plane0, org) - plane0[3]) / d3_dot(plane0, dir);
+ t[1] = (-d3_dot(plane1, org) - plane1[3]) / d3_dot(plane1, dir);
+
+ /* Check the returned distances against the computed ones */
+ CHK(eq_eps(hit_range[0], t[0], 1.e-6*hit_range[0]));
+ CHK(eq_eps(hit_range[1], t[1], 1.e-6*hit_range[1]));
+
+ /* Trace a ray that starts into the laser sheet */
+ range[0] = 0.5*(hit_range[0] + hit_range[1]);
+ htrdr_combustion_laser_trace_ray(laser, org, dir, range, hit_range);
+ CHK(hit_range[0] < hit_range[1]);
+ CHK(hit_range[0] == range[0]);
+ CHK(eq_eps(hit_range[1], t[1], 1.e-6*hit_range[1]));
+
+ htrdr_ref_put(htrdr);
+ htrdr_combustion_laser_ref_put(laser);
+ htrdr_mpi_finalize();
+
+ return 0;
+}
+
diff --git a/src/commands/htrdr_cmd.c b/src/commands/htrdr_cmd.c
@@ -16,6 +16,7 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "atmosphere/htrdr_atmosphere.h"
+#include "combustion/htrdr_combustion.h"
#include "core/htrdr_version.h"
#include <string.h>
@@ -76,7 +77,8 @@ main(int argc, char** argv)
/* Combustion mode */
} else if(!strcmp(argv[1], "combustion")) {
- /* TODO */
+ err = htrdr_combustion_main(argc-1, argv+1);
+ if(err) goto error;
/* Version */
} else if(!strcmp(argv[1], "--version")) {
diff --git a/src/commands/htrdr_combustion_cmd.c b/src/commands/htrdr_combustion_cmd.c
@@ -0,0 +1,24 @@
+/* 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 "combustion/htrdr_combustion.h"
+
+int
+main(int argc, char** argv)
+{
+ return htrdr_combustion_main(argc, argv);
+}
diff --git a/src/core/htrdr_args.c b/src/core/htrdr_args.c
@@ -23,6 +23,7 @@
#include <rsys/cstr.h>
#include <rsys/double3.h>
+#include <rsys/str.h>
#include <string.h>
@@ -74,7 +75,7 @@ parse_fov(const char* str, double* out_fov)
}
static res_T
-parse_image_parameter(void* args, const char* str)
+parse_image_parameter(const char* str, void* args)
{
char buf[128];
struct htrdr_args_image* img = args;
@@ -136,7 +137,7 @@ error:
}
static res_T
-parse_camera_parameter(void* args, const char* str)
+parse_camera_parameter(const char* str, void* args)
{
char buf[128];
struct htrdr_args_camera* cam = args;
@@ -190,7 +191,7 @@ error:
}
static res_T
-parse_rectangle_parameter(void* args, const char* str)
+parse_rectangle_parameter(const char* str, void* args)
{
char buf[128];
struct htrdr_args_rectangle* rect = args;
@@ -208,7 +209,6 @@ parse_rectangle_parameter(void* args, const char* str)
}
strncpy(buf, str, sizeof(buf));
- /* pos=0,0,10.1; key <- pos, val <- 0,0,10 */
key = strtok_r(buf, "=", &ctx);
val = strtok_r(NULL, "", &ctx);
@@ -269,7 +269,7 @@ error:
}
static res_T
-parse_spectral_parameter(void* ptr, const char* str)
+parse_spectral_parameter(const char* str, void* ptr)
{
char buf[128];
struct htrdr_args_spectral* args = ptr;
@@ -329,34 +329,64 @@ error:
}
static res_T
-parse_multiple_parameters
- (void* args,
- const char* str,
- res_T (*parse_parameter)(void* args, const char* str))
+parse_geometry_parameter(const char* str, void* ptr)
{
- char buf[512];
- char* tk;
+ struct str buf;
+ struct htrdr_args_geometry* geom = ptr;
+ char* key;
+ char* val;
char* ctx;
res_T res = RES_OK;
- ASSERT(args && str);
+ ASSERT(geom && str);
- if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
- fprintf(stderr, "Could not duplicate the option string `%s'.\n", str);
- res = RES_MEM_ERR;
+ str_init(NULL, &buf);
+ res = str_set(&buf, str);
+ if(res != RES_OK) {
+ fprintf(stderr,
+ "Could not duplicate the geometry option string `%s' -- %s.\n",
+ str, res_to_cstr(res));
goto error;
}
- strncpy(buf, str, sizeof(buf));
- tk = strtok_r(buf, ":", &ctx);
- do {
- res = parse_parameter(args, tk);
- if(res != RES_OK) goto error;
- tk = strtok_r(NULL, ":", &ctx);
- } while(tk);
+ key = strtok_r(str_get(&buf), "=", &ctx);
+ val = strtok_r(NULL, "", &ctx);
+
+ if(!val) {
+ fprintf(stderr, "Missing value to the geometry parameter `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ htrdr_args_geometry_free(geom);
+
+ #define SET_VALUE(Key, Val, Str) { \
+ const size_t len = strlen(Val) + 1; \
+ Str = mem_alloc(len); \
+ if(!Str) { \
+ fprintf(stderr, \
+ "Could not duplicate the value `%s' of the geometry parameter `%s'.\n",\
+ (Val), (Key)); \
+ res = RES_MEM_ERR; \
+ goto error; \
+ } \
+ strncpy(Str, Val, len); \
+ } (void)0
+ if(!strcmp(key, "obj")) {
+ SET_VALUE(key, val, geom->path_obj);
+ } else if(!strcmp(key, "mats")) {
+ SET_VALUE(key, val, geom->path_mats);
+ } else {
+ fprintf(stderr, "Invalid geometry parameter `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ #undef SET_VALUE
exit:
+ str_release(&buf);
return res;
error:
+ htrdr_args_geometry_free(geom);
goto exit;
}
@@ -367,27 +397,68 @@ res_T
htrdr_args_camera_parse(struct htrdr_args_camera* cam, const char* str)
{
if(!cam || !str) return RES_BAD_ARG;
- return parse_multiple_parameters(cam, str, parse_camera_parameter);
+ return cstr_parse_list(str, ':', parse_camera_parameter, cam);
}
res_T
htrdr_args_rectangle_parse(struct htrdr_args_rectangle* rect, const char* str)
{
if(!rect || !str) return RES_BAD_ARG;
- return parse_multiple_parameters(rect, str, parse_rectangle_parameter);
+ return cstr_parse_list(str, ':', parse_rectangle_parameter, rect);
}
-
+
res_T
htrdr_args_image_parse(struct htrdr_args_image* img, const char* str)
{
if(!img || !str) return RES_BAD_ARG;
- return parse_multiple_parameters(img, str, parse_image_parameter);
+ return cstr_parse_list(str, ':', parse_image_parameter, img);
}
res_T
htrdr_args_spectral_parse(struct htrdr_args_spectral* spectral, const char* str)
{
if(!spectral || !str) return RES_BAD_ARG;
- return parse_multiple_parameters(spectral, str, parse_spectral_parameter);
+ return cstr_parse_list(str, ':', parse_spectral_parameter, spectral);
+}
+
+void
+htrdr_args_geometry_free(struct htrdr_args_geometry* geom)
+{
+ ASSERT(geom);
+ if(geom->path_obj) mem_rm(geom->path_obj);
+ if(geom->path_mats) mem_rm(geom->path_mats);
+ *geom = HTRDR_ARGS_GEOMETRY_NULL;
+}
+
+res_T
+htrdr_args_geometry_parse(struct htrdr_args_geometry* geom, const char* str)
+{
+ res_T res = RES_OK;
+
+ if(!geom || !str) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = cstr_parse_list(str, ':', parse_geometry_parameter, geom);
+ if(res != RES_OK) goto error;
+
+ if(!geom->path_obj) {
+ fprintf(stderr, "Missing the `obj' geometry parameter.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(!geom->path_mats) {
+ fprintf(stderr, "Missing the `mats' geometry parameter.\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ if(geom) htrdr_args_geometry_free(geom);
+ goto exit;
}
diff --git a/src/core/htrdr_args.h.in b/src/core/htrdr_args.h.in
@@ -67,6 +67,15 @@ struct htrdr_args_image {
static const struct htrdr_args_image HTRDR_ARGS_IMAGE_DEFAULT =
HTRDR_ARGS_IMAGE_DEFAULT__;
+/* Arguments of a geometry */
+struct htrdr_args_geometry {
+ char* path_obj; /* Path toward a htrdr-obj */
+ char* path_mats; /* Path toward a htrdr-materials */
+};
+#define HTRDR_ARGS_GEOMETRY_NULL__ {NULL, NULL}
+static const struct htrdr_args_geometry HTRDR_ARGS_GEOMETRY_NULL =
+ HTRDR_ARGS_GEOMETRY_NULL__;
+
/* Arguments of the spectral domain */
struct htrdr_args_spectral {
double wlen_range[2]; /* Spectral range of integration in nm */
@@ -84,18 +93,18 @@ static const struct htrdr_args_spectral HTRDR_ARGS_SPECTRAL_DEFAULT =
/*******************************************************************************
* Exported functions
******************************************************************************/
-BEGIN_DECLS
+BEGIN_DECLS
HTRDR_CORE_API res_T
htrdr_args_camera_parse
- (struct htrdr_args_camera* cam,
+ (struct htrdr_args_camera* cam,
const char* str);
HTRDR_CORE_API res_T
htrdr_args_rectangle_parse
(struct htrdr_args_rectangle* rect,
const char* str);
-
+
HTRDR_CORE_API res_T
htrdr_args_image_parse
(struct htrdr_args_image* img,
@@ -106,6 +115,15 @@ htrdr_args_spectral_parse
(struct htrdr_args_spectral* spectral,
const char* str);
+HTRDR_CORE_API void
+htrdr_args_geometry_free
+ (struct htrdr_args_geometry* geom);
+
+HTRDR_CORE_API res_T
+htrdr_args_geometry_parse
+ (struct htrdr_args_geometry* geom,
+ const char* str);
+
END_DECLS
#endif /* HTRDR_ARGS_H */
diff --git a/src/core/htrdr_buffer.h b/src/core/htrdr_buffer.h
@@ -27,6 +27,14 @@
* Row major ordered 2D buffer
*/
+struct htrdr_pixel_format {
+ size_t size; /* In bytes */
+ size_t alignment; /* Power of two, in Bytes */
+};
+#define HTRDR_PIXEL_FORMAT_NULL__ {0, 0}
+static const struct htrdr_pixel_format HTRDR_PIXEL_FORMAT_NULL =
+ HTRDR_PIXEL_FORMAT_NULL__;
+
struct htrdr_buffer_layout {
size_t width; /* #elements in X */
size_t height; /* #elements in Y */
diff --git a/src/core/htrdr_draw_map.c b/src/core/htrdr_draw_map.c
@@ -635,7 +635,7 @@ draw_map
* is initialised in order to ensure that an unique and predictable set of
* random numbers is used for the current tile. */
SSP(rng_proxy_create2
- (&htrdr->lifo_allocators[ithread],
+ (htrdr_get_thread_allocator(htrdr, (size_t)ithread),
&ssp_rng_threefry,
RNG_SEQUENCE_SIZE * (size_t)mcode, /* Offset */
RNG_SEQUENCE_SIZE, /* Size */
diff --git a/src/core/htrdr_draw_map.h b/src/core/htrdr_draw_map.h
@@ -52,7 +52,6 @@ typedef void
struct htrdr_draw_map_args {
htrdr_draw_pixel_T draw_pixel;
struct htrdr_buffer_layout buffer_layout;
- struct ssp_rng* rng;
size_t spp; /* Samples per pixel */
void* context; /* User defined data */
};
@@ -60,7 +59,6 @@ struct htrdr_draw_map_args {
#define HTRDR_DRAW_MAP_ARGS_NULL__ { \
NULL, /* Draw pixel functor */ \
HTRDR_BUFFER_LAYOUT_NULL__, /* Layout of the destination buffer */ \
- NULL, /* Random Number Generator */ \
0, /* #Samples per pixel */ \
NULL /* User defined data */ \
}
diff --git a/src/core/htrdr_ran_wlen.c b/src/core/htrdr_ran_wlen.c
@@ -201,9 +201,10 @@ wlen_ran_sample_continue
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.0e9; /* Convert in nanometers */
+ return lambda_m * 1.e9; /* Convert in nanometers */
}
static double
@@ -249,7 +250,7 @@ wlen_ran_sample_discrete
/* Uniformly sample a wavelength in the sampled band */
lambda = band_range[0] + (band_range[1] - band_range[0]) * r1;
- pdf_continue = 1.0 / ((band_range[1] - band_range[0])*1.e-9);
+ pdf_continue = 1.0 / (band_range[1] - band_range[0]);
if(pdf) {
*pdf = pdf_band * pdf_continue;
diff --git a/src/core/htrdr_ran_wlen.h b/src/core/htrdr_ran_wlen.h
@@ -53,7 +53,7 @@ htrdr_ran_wlen_sample
(const struct htrdr_ran_wlen* wlen_ran,
const double r0, /* Canonical number in [0, 1[ */
const double r1, /* Canonical number in [0, 1[ */
- double* pdf); /* May be NULL */
+ double* pdf); /* In nm^-1. May be NULL */
END_DECLS
diff --git a/src/core/htrdr_rectangle.c b/src/core/htrdr_rectangle.c
@@ -20,6 +20,7 @@
#include "core/htrdr_rectangle.h"
#include <rsys/double3.h>
+#include <rsys/double33.h>
#include <rsys/mem_allocator.h>
#include <rsys/ref_count.h>
@@ -27,8 +28,11 @@ struct htrdr_rectangle {
/* Frame of the rectangle in world space */
double axis_x[3];
double axis_y[3];
-
double normal[3];
+
+ double local2world[12]; /* Rectangle to world transformation matrix */
+ double world2local[12]; /* World to rectangle transformation matrix */
+
double position[3]; /* Center of the rectangle */
struct htrdr* htrdr;
ref_T ref;
@@ -50,7 +54,7 @@ rectangle_release(ref_T* ref)
}
/*******************************************************************************
- * Local fuuction
+ * Exported functions
******************************************************************************/
res_T
htrdr_rectangle_create
@@ -63,12 +67,13 @@ htrdr_rectangle_create
{
struct htrdr_rectangle* rect = NULL;
double x[3], y[3], z[3];
+ double trans[3];
res_T res = RES_OK;
ASSERT(htrdr && pos && tgt && up && sz && out_rect);
rect = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*rect));
if(!rect) {
- htrdr_log_err(htrdr, "could not allocate the rectangle data structure.\n");
+ htrdr_log_err(htrdr, "Could not allocate the rectangle data structure.\n");
res = RES_MEM_ERR;
goto error;
}
@@ -78,7 +83,7 @@ htrdr_rectangle_create
if(sz[0] <= 0 || sz[1] <= 0) {
htrdr_log_err(htrdr,
- "invalid rectangle size `%g %g'. It must be strictly positive.\n",
+ "Invalid rectangle size `%g %g'. It must be strictly positive.\n",
SPLIT2(sz));
res = RES_BAD_ARG;
goto error;
@@ -87,7 +92,7 @@ htrdr_rectangle_create
if(d3_normalize(z, d3_sub(z, tgt, pos)) <= 0
|| d3_normalize(x, d3_cross(x, z, up)) <= 0
|| d3_normalize(y, d3_cross(y, z, x)) <= 0) {
- htrdr_log_err(htrdr, "invalid rectangle frame:\n"
+ htrdr_log_err(htrdr, "Invalid rectangle frame:\n"
"\tposition = %g %g %g\n"
"\ttarget = %g %g %g\n"
"\tup = %g %g %g\n",
@@ -96,6 +101,21 @@ htrdr_rectangle_create
goto error;
}
+ /* Setup the local to world transformation matrix */
+ d3_set(rect->local2world+0, x);
+ d3_set(rect->local2world+3, y);
+ d3_set(rect->local2world+6, z);
+ d3_set(rect->local2world+9, pos);
+
+ /* Inverse the local to world transformation matrix. Note that since the
+ * represented frame is orthonormal one can transpose the original rotation
+ * matrix to cheaply compute its inverse */
+ d33_transpose(rect->world2local, rect->local2world);
+
+ /* Compute the affine inverse transform */
+ d3_minus(trans, pos);
+ d33_muld3(rect->world2local+9, rect->world2local, trans);
+
d3_muld(rect->axis_x, x, sz[0]*0.5);
d3_muld(rect->axis_y, y, sz[1]*0.5);
d3_set(rect->normal, z);
@@ -146,3 +166,35 @@ htrdr_rectangle_get_normal(const struct htrdr_rectangle* rect, double normal[3])
d3_set(normal, rect->normal);
}
+void
+htrdr_rectangle_get_center(const struct htrdr_rectangle* rect, double pos[3])
+{
+ ASSERT(rect && pos);
+ d3_set(pos, rect->position);
+}
+
+double*
+htrdr_rectangle_get_transform
+ (const struct htrdr_rectangle* rect,
+ double transform[12])
+{
+ ASSERT(rect && transform);
+ d3_set(transform+0, rect->local2world+0);
+ d3_set(transform+3, rect->local2world+3);
+ d3_set(transform+6, rect->local2world+6);
+ d3_set(transform+9, rect->local2world+9);
+ return transform;
+}
+
+double*
+htrdr_rectangle_get_transform_inverse
+ (const struct htrdr_rectangle* rect,
+ double transform_inverse[12])
+{
+ ASSERT(rect && transform_inverse);
+ d3_set(transform_inverse+0, rect->world2local+0);
+ d3_set(transform_inverse+3, rect->world2local+3);
+ d3_set(transform_inverse+6, rect->world2local+6);
+ d3_set(transform_inverse+9, rect->world2local+9);
+ return transform_inverse;
+}
diff --git a/src/core/htrdr_rectangle.h b/src/core/htrdr_rectangle.h
@@ -32,7 +32,7 @@ htrdr_rectangle_create
(struct htrdr* htrdr,
const double sz[2], /* Size of the rectangle along its local X and Y axis */
const double pos[3], /* World space position of the rectangle center */
- const double tgt[3], /* Vector orthognal to the rectangle plane */
+ const double tgt[3], /* Targeted point */
const double up[3], /* vector orthogonal to the rectangle X axis */
struct htrdr_rectangle** rect);
@@ -55,6 +55,21 @@ htrdr_rectangle_get_normal
(const struct htrdr_rectangle* rect,
double normal[3]);
+HTRDR_CORE_API void
+htrdr_rectangle_get_center
+ (const struct htrdr_rectangle* rect,
+ double center[3]);
+
+HTRDR_CORE_API double*
+htrdr_rectangle_get_transform
+ (const struct htrdr_rectangle* rect,
+ double transform[12]);
+
+HTRDR_CORE_API double*
+htrdr_rectangle_get_transform_inverse
+ (const struct htrdr_rectangle* rect,
+ double transform_inverse[12]);
+
END_DECLS
#endif /* HTRDR_RECTANGLE_H */
