commit f5071d2e814514fdf8a12d8f8ad3002eac9ced6e
parent d8f9063b539c78539677999b4f6ef1d373716665
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 18 May 2018 09:21:19 +0200
Add the -r <rectangle> and -I <image> options
Diffstat:
11 files changed, 609 insertions(+), 25 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -49,11 +49,15 @@ set(HTRDR_FILES_SRC
htrdr.c
htrdr_main.c
htrdr_args.c
- htrdr_clouds.c)
+ htrdr_clouds.c
+ htrdr_rectangle.c
+ htrdr_transmission.c)
set(HTRDR_FILES_INC
htrdr.h
htrdr_args.h
- htrdr_clouds.h)
+ htrdr_clouds.h
+ htrdr_rectangle.h
+ htrdr_solve.h)
set(HTDRD_FILES_DOC COPYING README.md)
# Prepend each file in the `HTRDR_FILES_<SRC|INC>' list by `HTRDR_SOURCE_DIR'
@@ -63,6 +67,9 @@ rcmake_prepend_path(HTRDR_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_executable(htrdr ${HTRDR_FILES_SRC} ${HTRDR_FILES_INC})
target_link_libraries(htrdr HTCP RSys SVX)
+if(CMAKE_COMPILER_IS_GNUCC)
+ target_link_libraries(htrdr m)
+endif()
set_target_properties(htrdr PROPERTIES
DEFINE_SYMBOL HTRDR_SHARED_BUILD
diff --git a/doc/cli.txt b/doc/cli.txt
@@ -0,0 +1,20 @@
+<image> ::= <imgopt>[:<imgopt> ... ]
+<rectangle> ::= <rectopt>[:<rectopt> ... ]
+
+<imgopt> ::= <definition>
+ | <#samples>
+<definition> ::= def=INTEGERxINTEGER # default 32x32
+<#samples> ::= ssp=INTEGER # default 1
+
+<rectopt> ::= <target>
+ | <position>
+ | <up>
+ | <size>
+<target> ::= tgt=<float3>
+<position> ::= pos=<float3>
+<up> ::= up=<float3>
+<size> ::= sz=<float2>
+
+<float3> ::= FLOAT,FLOAT,FLOAT
+<float2> ::= FLOAT,FLOAT
+<int2> ::= INTEGER
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -18,6 +18,8 @@
#include "htrdr.h"
#include "htrdr_args.h"
#include "htrdr_clouds.h"
+#include "htrdr_rectangle.h"
+#include "htrdr_solve.h"
#include <rsys/mem_allocator.h>
@@ -119,21 +121,40 @@ error:
goto exit;
}
+static void
+release_htrdr(ref_T* ref)
+{
+ struct htrdr* htrdr = CONTAINER_OF(ref, struct htrdr, ref);
+ ASSERT(htrdr);
+ if(htrdr->svx) SVX(device_ref_put(htrdr->svx));
+ if(htrdr->clouds) SVX(tree_ref_put(htrdr->clouds));
+ if(htrdr->rect) htrdr_rectangle_ref_put(htrdr->rect);
+ logger_release(&htrdr->logger);
+ MEM_RM(htrdr->allocator, htrdr);
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
res_T
-htrdr_init
- (struct mem_allocator* allocator,
+htrdr_create
+ (struct mem_allocator* mem_allocator,
const struct htrdr_args* args,
- struct htrdr* htrdr)
+ struct htrdr** out_htrdr)
{
+ struct htrdr* htrdr = NULL;
+ struct mem_allocator* allocator = NULL;
res_T res = RES_OK;
- ASSERT(args && htrdr);
-
- memset(htrdr, 0, sizeof(*htrdr));
+ ASSERT(args && out_htrdr);
- htrdr->allocator = allocator ? allocator : &mem_default_allocator;
+ allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
+ htrdr = MEM_CALLOC(allocator, 1, sizeof(*htrdr));
+ if(!htrdr) {
+ fprintf(stderr, "Could not allocat the htrdr main structure.\n");
+ goto error;
+ }
+ ref_init(&htrdr->ref);
+ htrdr->allocator = allocator;
logger_init(htrdr->allocator, &htrdr->logger);
logger_set_stream(&htrdr->logger, LOG_ERROR, print_err, NULL);
@@ -142,6 +163,11 @@ htrdr_init
htrdr->dump_vtk = args->dump_vtk;
htrdr->verbose = args->verbose;
+ /* Integration plane */
+ res = htrdr_rectangle_create(htrdr, args->rectangle.pos, args->rectangle.tgt,
+ args->rectangle.up, args->rectangle.sz, &htrdr->rect);
+ if(res != RES_OK) goto error;
+
res = svx_device_create(&htrdr->logger, allocator, args->verbose, &htrdr->svx);
if(res != RES_OK) {
htrdr_log_err(htrdr, "could not create the Star-VoXel device.\n");
@@ -158,19 +184,28 @@ htrdr_init
if(res != RES_OK) goto error;
exit:
+ *out_htrdr = htrdr;
return res;
error:
- htrdr_release(htrdr);
+ if(htrdr) {
+ htrdr_ref_put(htrdr);
+ htrdr = NULL;
+ }
goto exit;
}
void
-htrdr_release(struct htrdr* htrdr)
+htrdr_ref_get(struct htrdr* htrdr)
{
ASSERT(htrdr);
- if(htrdr->svx) SVX(device_ref_put(htrdr->svx));
- if(htrdr->clouds) SVX(tree_ref_put(htrdr->clouds));
- logger_release(&htrdr->logger);
+ ref_get(&htrdr->ref);
+}
+
+void
+htrdr_ref_put(struct htrdr* htrdr)
+{
+ ASSERT(htrdr);
+ ref_put(&htrdr->ref, release_htrdr);
}
res_T
@@ -180,6 +215,13 @@ htrdr_run(struct htrdr* htrdr)
if(htrdr->dump_vtk) {
res = clouds_dump_vtk(htrdr, htrdr->output);
if(res != RES_OK) goto error;
+ } else {
+ double pos[3] = {3.4,2.2,0};
+ double dir[3] = {0,0,1};
+ double val = 0;
+ res = htrdr_solve_transmission(htrdr, pos, dir, &val);
+ if(res != RES_OK) goto error;
+ printf(">>>> T = %g\n", val);
}
exit:
return res;
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -17,6 +17,7 @@
#define HTRDR_H
#include <rsys/logger.h>
+#include <rsys/ref_count.h>
/* Forward declaration */
struct svx_device;
@@ -28,22 +29,30 @@ struct htrdr {
struct svx_device* svx;
struct svx_tree* clouds;
+ struct htrdr_rectangle* rect;
+
FILE* output;
int dump_vtk;
int verbose;
+
struct logger logger;
struct mem_allocator* allocator;
+ ref_T ref;
};
extern LOCAL_SYM res_T
-htrdr_init
+htrdr_create
(struct mem_allocator* allocator,
const struct htrdr_args* args,
- struct htrdr* htrdr);
+ struct htrdr** htrdr);
+
+extern LOCAL_SYM void
+htrdr_ref_get
+ (struct htrdr* htrdr);
extern LOCAL_SYM void
-htrdr_release
+htrdr_ref_put
(struct htrdr* htrdr);
extern LOCAL_SYM res_T
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -13,8 +13,14 @@
* 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 2 /* strtok_r support */
+
#include "htrdr_args.h"
+
+#include <rsys/cstr.h>
+
#include <getopt.h>
+#include <string.h>
/*******************************************************************************
* Helper functions
@@ -33,9 +39,13 @@ print_help(const char* cmd)
printf(
" -i INPUT path of the input HTCP file.\n");
printf(
+" -I <image> define the image to compute.\n");
+ printf(
" -o OUTPUT file where data are written. If not defined, data are\n"
" written to standard output.\n");
printf(
+" -r <rectangle> define the integration plane.\n");
+ printf(
" -v make the program more verbose.\n");
printf("\n");
printf(
@@ -45,6 +55,166 @@ print_help(const char* cmd)
cmd);
}
+static INLINE res_T
+parse_doubleX(const char* str, double* val, const size_t sz)
+{
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(str && val);
+ res = cstr_to_list_double(str, ',', val, &len, sz);
+ if(res == RES_OK && len != sz) res = RES_BAD_ARG;
+ return res;
+}
+
+static INLINE res_T
+parse_definition(const char* str, unsigned val[2])
+{
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(str && val);
+ res = cstr_to_list_uint(str, 'x', val, &len, 2);
+ if(res != RES_OK) return res;
+ if(len != 2) return RES_BAD_ARG;
+ if(val[0] > 16384 || val[1] > 16384) return RES_BAD_ARG;
+ return RES_OK;
+}
+
+static res_T
+parse_rectangle_parameter(struct htrdr_args* args, const char* str)
+{
+ char buf[128];
+ char* key;
+ char* val;
+ char* ctx;
+ res_T res;
+ ASSERT(str && args);
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr,
+ "Could not duplicate the rectangle option string `%s'.\n", str);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ strncpy(buf, str, sizeof(buf));
+
+ key = strtok_r(buf, "=", &ctx);
+ val = strtok_r(NULL, "", &ctx);
+
+ if(!val) {
+ fprintf(stderr, "Missing a value to the rectangle option `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ #define PARSE(Name, Func) \
+ res = Func; \
+ if(res != RES_OK) { \
+ fprintf(stderr, "Invalid rectangle "Name" `%s'.\n", val); \
+ goto error; \
+ } (void)0
+ if(!strcmp(key, "pos")) {
+ PARSE("position", parse_doubleX(val, args->rectangle.pos, 3));
+ } else if(!strcmp(key, "tgt")) {
+ PARSE("target", parse_doubleX(val, args->rectangle.tgt, 3));
+ } else if(!strcmp(key, "up")) {
+ PARSE("up vector", parse_doubleX(val, args->rectangle.up, 3));
+ } else if(!strcmp(key, "sz")) {
+ PARSE("size", parse_doubleX(val, args->rectangle.sz, 2));
+ } else {
+ fprintf(stderr, "Invalid rectangle parameter `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ #undef PARSE
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_image_parameter(struct htrdr_args* args, const char* str)
+{
+ char buf[128];
+ char* key;
+ char* val;
+ char* ctx;
+ res_T res;
+ ASSERT(str && args);
+
+ if(strlen(str) >= sizeof(buf) -1/*NULL char*/) {
+ fprintf(stderr,
+ "Could not duplicate the image 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(!val) {
+ fprintf(stderr, "Missing a value to the image option `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ #define PARSE(Name, Func) \
+ res = Func; \
+ if(res != RES_OK) { \
+ fprintf(stderr, "Invalid image "Name" `%s'.\n", val); \
+ goto error; \
+ } (void)0
+ if(!strcmp(key, "def")) {
+ PARSE("definition", parse_definition(val, args->image.definition));
+ } else if(!strcmp(key, "ssp")) {
+ PARSE("#samples per pixel", cstr_to_uint(val, &args->image.spp));
+ } else {
+ fprintf(stderr, "Invalid image parameter `%s'.\n", key);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ #undef PARSE
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+parse_multiple_parameters
+ (struct htrdr_args* args,
+ const char* str,
+ res_T (*parse_parameter)(struct htrdr_args* args, const char* str))
+{
+ char buf[512];
+ char* tk;
+ char* ctx;
+ res_T res = RES_OK;
+ ASSERT(args && str);
+
+ if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) {
+ fprintf(stderr, "Could not duplicate the option string `%s'.\n", str);
+ res = RES_MEM_ERR;
+ 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);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -57,7 +227,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
*args = HTRDR_ARGS_DEFAULT;
- while((opt = getopt(argc, argv, "dfhi:o:v")) != -1) {
+ while((opt = getopt(argc, argv, "dfhI:i:o:r:v")) != -1) {
switch(opt) {
case 'd': args->dump_vtk = 1; break;
case 'f': args->force_overwriting = 1; break;
@@ -66,8 +236,16 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
htrdr_args_release(args);
args->quit = 1;
goto exit;
+ case 'I':
+ res = parse_multiple_parameters
+ (args, optarg, parse_image_parameter);
+ break;
case 'i': args->input = optarg; break;
case 'o': args->output = optarg; break;
+ case 'r':
+ res = parse_multiple_parameters
+ (args, optarg, parse_rectangle_parameter);
+ break;
case 'v': args->verbose = 1; break;
default: res = RES_BAD_ARG; break;
}
diff --git a/src/htrdr_args.h b/src/htrdr_args.h
@@ -22,12 +22,41 @@ struct htrdr_args {
const char* input;
const char* output;
+ struct {
+ double pos[3]; /* Position */
+ double tgt[3]; /* Target */
+ double up[3]; /* Up vector */
+ double sz[2]; /* Plane size in world space */
+ } rectangle;
+
+ struct {
+ unsigned definition[2]; /* #pixels in X and Y */
+ unsigned spp; /* #samples per pixel */
+ } image;
+
int force_overwriting;
int dump_vtk; /* Dump the loaded cloud properties in a VTK file */
int verbose; /* Verbosity level */
int quit; /* Qui the application */
};
-#define HTRDR_ARGS_DEFAULT__ {NULL, NULL, 0, 0, 0, 0}
+
+#define HTRDR_ARGS_DEFAULT__ { \
+ NULL, /* Input filename */ \
+ NULL, /* Output filename */ \
+ { \
+ {0, 0, 0}, /* plane position */ \
+ {0, 0, 1}, /* plane target */ \
+ {0, 1, 0}, /* plane up */ \
+ {1, 1}, /* plane size */ \
+ }, { \
+ {32, 32}, /* image definition */ \
+ 1 /* #samples per pixel */ \
+ }, \
+ 0, /* Force overwriting */ \
+ 0, /* dump VTK */ \
+ 0, /* Verbose flag */ \
+ 0 /* Quit the application */ \
+}
static const struct htrdr_args HTRDR_ARGS_DEFAULT = HTRDR_ARGS_DEFAULT__;
extern LOCAL_SYM res_T
diff --git a/src/htrdr_main.c b/src/htrdr_main.c
@@ -24,10 +24,9 @@
int
main(int argc, char** argv)
{
- struct htrdr htrdr;
+ struct htrdr* htrdr = NULL;
struct htrdr_args args;
size_t memsz = 0;
- int htrdr_is_init = 0;
int err = 0;
res_T res = RES_OK;
@@ -35,15 +34,15 @@ main(int argc, char** argv)
if(res != RES_OK) goto error;
if(args.quit) goto exit;
- res = htrdr_init(NULL, &args, &htrdr);
+ res = htrdr_create(NULL, &args, &htrdr);
if(res != RES_OK) goto error;
- htrdr_is_init = 1;
- res = htrdr_run(&htrdr);
+ res = htrdr_run(htrdr);
if(res != RES_OK) goto error;
exit:
- if(htrdr_is_init) htrdr_release(&htrdr);
+ if(htrdr) htrdr_ref_put(htrdr);
+
htrdr_args_release(&args);
if((memsz = mem_allocated_size()) != 0) {
fprintf(stderr, "Memory leaks: %lu Bytes\n", (unsigned long)memsz);
diff --git a/src/htrdr_rectangle.c b/src/htrdr_rectangle.c
@@ -0,0 +1,136 @@
+/* Copyright (C) 2018 Université Paul Sabatier, |Meso|Star>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "htrdr.h"
+#include "htrdr_rectangle.h"
+
+#include <rsys/double3.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+struct htrdr_rectangle {
+ /* Frame of the rectangle in world space */
+ double axis_x[3];
+ double axis_y[3];
+
+ double position[3]; /* Center of the rectangle */
+ struct htrdr* htrdr;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+rectangle_release(ref_T* ref)
+{
+ struct htrdr_rectangle* rect;
+ struct htrdr* htrdr;
+ ASSERT(ref);
+ rect = CONTAINER_OF(ref, struct htrdr_rectangle, ref);
+ htrdr = rect->htrdr;
+ MEM_RM(htrdr->allocator, rect);
+ htrdr_ref_put(htrdr);
+}
+
+/*******************************************************************************
+ * Local fuuction
+ ******************************************************************************/
+res_T
+htrdr_rectangle_create
+ (struct htrdr* htrdr,
+ const double pos[3],
+ const double tgt[3],
+ const double up[2],
+ const double sz[2],
+ struct htrdr_rectangle** out_rect)
+{
+ struct htrdr_rectangle* rect = NULL;
+ double x[3], y[3], z[3];
+ res_T res = RES_OK;
+ ASSERT(htrdr && pos && tgt && up && sz && out_rect);
+
+ rect = MEM_CALLOC(htrdr->allocator, 1, sizeof(*rect));
+ if(!rect) {
+ htrdr_log_err(htrdr, "could not allocate the rectangle data structure.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&rect->ref);
+ htrdr_ref_get(htrdr);
+ rect->htrdr = htrdr;
+
+ if(sz[0] <= 0 || sz[1] <= 1) {
+ htrdr_log_err(htrdr,
+ "invalid rectangle size `%g %g'. It must be strictly positive.\n",
+ SPLIT2(sz));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ 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"
+ "\tposition = %g %g %g\n"
+ "\ttarget = %g %g %g\n"
+ "\tup = %g %g %g\n",
+ SPLIT3(pos), SPLIT3(tgt), SPLIT3(up));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ d3_muld(rect->axis_x, x, sz[0]*0.5);
+ d3_muld(rect->axis_y, y, sz[1]*0.5);
+ d3_set(rect->position, pos);
+
+exit:
+ *out_rect = rect;
+ return res;
+error:
+ if(rect) {
+ htrdr_rectangle_ref_put(rect);
+ rect = NULL;
+ }
+ goto exit;
+}
+
+void
+htrdr_rectangle_sample_pos
+ (const struct htrdr_rectangle* rect,
+ const double sample[2], /* In [0, 1[ */
+ double pos[3])
+{
+ double x[3], y[3];
+ ASSERT(rect && sample && pos);
+ d3_muld(x, rect->axis_x, sample[0]*2.0 - 1.0);
+ d3_muld(y, rect->axis_y, sample[1]*2.0 - 1.0);
+ d3_add(pos, d3_add(pos, rect->position, x), y);
+}
+
+void
+htrdr_rectangle_ref_get(struct htrdr_rectangle* rect)
+{
+ ASSERT(rect);
+ ref_get(&rect->ref);
+}
+
+void
+htrdr_rectangle_ref_put(struct htrdr_rectangle* rect)
+{
+ ASSERT(rect);
+ ref_put(&rect->ref, rectangle_release);
+}
+
diff --git a/src/htrdr_rectangle.h b/src/htrdr_rectangle.h
@@ -0,0 +1,48 @@
+/* Copyright (C) 2018 Université Paul Sabatier, |Meso|Star>
+ *
+ * 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_RECTANGLE_H
+#define HTRDR_RECTANGLE_H
+
+#include <rsys/rsys.h>
+
+struct htrdr;
+struct htrdr_rectangle; /* 2D rectangle transformed in 3D */
+
+extern LOCAL_SYM res_T
+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 up[2], /* vector orthogonal to the rectangle X axis */
+ struct htrdr_rectangle** rect);
+
+extern LOCAL_SYM void
+htrdr_rectangle_ref_get
+ (struct htrdr_rectangle* rect);
+
+extern LOCAL_SYM void
+htrdr_rectangle_ref_put
+ (struct htrdr_rectangle* rect);
+
+extern LOCAL_SYM void
+htrdr_rectangle_sample_pos
+ (const struct htrdr_rectangle* rect,
+ const double sample[2], /* In [0, 1[ */
+ double pos[3]);
+
+#endif /* HTRDR_RECTANGLE_H */
+
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -0,0 +1,32 @@
+/* Copyright (C) 2018 Université Paul Sabatier, |Meso|Star>
+ *
+ * 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_SOLVE_H
+#define HTRDR_SOLVE_H
+
+#include <rsys/rsys.h>
+
+/* Forward declaration */
+struct htrdr;
+
+extern LOCAL_SYM res_T
+htrdr_solve_transmission
+ (struct htrdr* htrdr,
+ const double pos[3],
+ const double dir[3],
+ double *val);
+
+#endif /* HTRDR_SOLVE_H */
+
diff --git a/src/htrdr_transmission.c b/src/htrdr_transmission.c
@@ -0,0 +1,84 @@
+/* Copyright (C) 2018 Université Paul Sabatier, |Meso|Star>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "htrdr.h"
+#include "htrdr_solve.h"
+
+#include <star/svx.h>
+#include <rsys/math.h>
+
+struct transmit_context {
+ double tau;
+};
+static const struct transmit_context TRANSMIT_CONTEXT_NULL = {0};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static int
+discard_hit
+ (const struct svx_hit* hit,
+ const double org[3],
+ const double dir[3],
+ const double range[2],
+ void* context)
+{
+ struct transmit_context* ctx = context;
+ double dst;
+ double k;
+ ASSERT(hit && ctx && !SVX_HIT_NONE(hit));
+ (void)org, (void)dir, (void)range;
+
+ k = *((double*)hit->voxel.data);
+ dst = hit->distance[1] - hit->distance[0];
+ ASSERT(dst >= 0);
+ ctx->tau += k*dst;
+ return 1;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+htrdr_solve_transmission
+ (struct htrdr* htrdr,
+ const double pos[3],
+ const double dir[3],
+ double *val)
+{
+ struct svx_hit hit = SVX_HIT_NULL;
+ struct transmit_context ctx = TRANSMIT_CONTEXT_NULL;
+ const double range[2] = {0, INF};
+ res_T res = RES_OK;
+ ASSERT(htrdr && pos && dir && val);
+
+ res = svx_octree_trace_ray(htrdr->clouds, pos, dir, range, NULL,
+ discard_hit, &ctx, &hit);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr,
+ "error computing the transmission for the position `%g %g %g' "
+ "along the direction `%g %g %g'.\n",
+ SPLIT3(pos), SPLIT3(dir));
+ goto error;
+ }
+ ASSERT(SVX_HIT_NONE(&hit));
+ *val = ctx.tau ? exp(-ctx.tau) : 1.0;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
