htrdr

Solving radiative transfer in heterogeneous media
git clone git://git.meso-star.fr/htrdr.git
Log | Files | Refs | README | LICENSE

commit 2a055919de25ebec3e983bc8b2d53686ee450795
parent db50d7467f11d45a20c61cc8473c4dc1ab9dc277
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 12 Mar 2021 16:45:38 +0100

Add the htrdr_combustion_laser data structure

Diffstat:
Mcmake/combustion/CMakeLists.txt | 4+++-
Msrc/combustion/htrdr_combustion.c | 14++++++--------
Msrc/combustion/htrdr_combustion_c.h | 3++-
Msrc/combustion/htrdr_combustion_compute_radiance_sw.c | 5+++--
Asrc/combustion/htrdr_combustion_laser.c | 165+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/combustion/htrdr_combustion_laser.h | 61+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/core/htrdr_rectangle.c | 53+++++++++++++++++++++++++++++++++++++++++++++++++----
Msrc/core/htrdr_rectangle.h | 10++++++++++
8 files changed, 299 insertions(+), 16 deletions(-)

diff --git a/cmake/combustion/CMakeLists.txt b/cmake/combustion/CMakeLists.txt @@ -51,12 +51,14 @@ set(HTRDR_COMBUSTION_FILES_SRC 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_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) diff --git a/src/combustion/htrdr_combustion.c b/src/combustion/htrdr_combustion.c @@ -18,6 +18,7 @@ #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" @@ -133,15 +134,12 @@ 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; - return htrdr_rectangle_create - (cmd->htrdr, - args->laser.size, - args->laser.position, - args->laser.target, - args->laser.up, - &cmd->laser); + laser_args.surface = args->laser; + return htrdr_combustion_laser_create(cmd->htrdr, &laser_args, &cmd->laser); } static res_T @@ -264,7 +262,7 @@ combustion_release(ref_T* ref) 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_rectangle_ref_put(cmd->laser); + 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); str_release(&cmd->output_name); diff --git a/src/combustion/htrdr_combustion_c.h b/src/combustion/htrdr_combustion_c.h @@ -29,6 +29,7 @@ struct atrstm; struct htrdr; struct htrdr_camera; +struct htrdr_combustion_laser; struct htrdr_geometry; struct htrdr_materials; struct htrdr_rectangle; @@ -51,7 +52,7 @@ struct htrdr_combustion { struct atrstm* medium; /* Semi transparent medium */ struct htrdr_camera* camera; /* Pinhole camera */ - struct htrdr_rectangle* laser; /* Laser surface emission */ + struct htrdr_combustion_laser* laser; /* Laser sheet */ double wavelength; /* Wavelength of the laser in nanometer */ struct htrdr_buffer_layout buf_layout; diff --git a/src/combustion/htrdr_combustion_compute_radiance_sw.c b/src/combustion/htrdr_combustion_compute_radiance_sw.c @@ -16,6 +16,7 @@ * 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 <astoria/atrstm.h> @@ -340,8 +341,8 @@ laser_once_scattered double L = 0; /* Radiance in W/m^2/sr/m */ ASSERT(cmd && rng && wlen >= 0 && pos && dir); - /* Find the intersections along dir with the laser sheet TODO */ - /* htrdr_combustion_laser_trace_ray(cmd->laser, pos, dir, range, laser_hit_dst); */ + /* 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(HIT_NONE(laser_hit_dst[0])) return 0; diff --git a/src/combustion/htrdr_combustion_laser.c b/src/combustion/htrdr_combustion_laser.c @@ -0,0 +1,165 @@ +/* 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 */ + 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); + + 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_BAD_ARG; + goto error; + } + ref_init(&laser->ref); + htrdr_ref_get(htrdr); + laser->htrdr = htrdr; + + 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[2] = -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; +} + diff --git a/src/combustion/htrdr_combustion_laser.h b/src/combustion/htrdr_combustion_laser.h @@ -0,0 +1,61 @@ +/* Copyright (C) 2018, 2019, 2020, 2021 |Meso|Star> (contact@meso-star.com) + * Copyright (C) 2018, 2019, 2021 CNRS + * Copyright (C) 2018, 2019, Université Paul Sabatier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#ifndef HTRDR_COMBUSTION_LASER_H +#define HTRDR_COMBUSTION_LASER_H + +#include "core/htrdr_args.h" + +#include <rsys/rsys.h> + +struct htrdr_combustion_laser_create_args { + struct htrdr_args_rectangle surface; /* Surface emission */ +}; +#define HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT__ { \ + HTRDR_ARGS_RECTANGLE_DEFAULT__ \ +} +static const struct htrdr_combustion_laser_create_args +HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT = + HTRDR_COMBUSTION_LASER_CREATE_ARGS_DEFAULT__; + +/* Forward declaration */ +struct htrdr; +struct htrdr_combustion_laser; + +extern LOCAL_SYM res_T +htrdr_combustion_laser_create + (struct htrdr* htrdr, + const struct htrdr_combustion_laser_create_args* args, + struct htrdr_combustion_laser** laser); + +extern LOCAL_SYM void +htrdr_combustion_laser_ref_get + (struct htrdr_combustion_laser* laser); + +extern LOCAL_SYM void +htrdr_combustion_laser_ref_put + (struct htrdr_combustion_laser* laser); + +extern LOCAL_SYM 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]); + +#endif /* HTRDR_COMBUSTION_LASER_H */ 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; @@ -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,28 @@ htrdr_rectangle_get_normal(const struct htrdr_rectangle* rect, double normal[3]) d3_set(normal, rect->normal); } +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 @@ -55,6 +55,16 @@ htrdr_rectangle_get_normal (const struct htrdr_rectangle* rect, double normal[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 */