htrdr

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

commit 986e8bd7cdef8f96b2ee38bd7997a2b45478f5e1
parent 12137598531c69e2df1c6a74e40dd0317f736e9e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 24 May 2018 15:15:07 +0200

Implement the htrdr_solve_transmission_buffer function

Diffstat:
Mcmake/CMakeLists.txt | 3++-
Msrc/htrdr.c | 24++++++++++++++++++++++++
Msrc/htrdr.h | 9+++++++--
Msrc/htrdr_args.h | 3+++
Msrc/htrdr_buffer.c | 31+++++++++++++++++++++++++++++--
Msrc/htrdr_solve.h | 4++++
Msrc/htrdr_transmission.c | 63+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7 files changed, 132 insertions(+), 5 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -25,6 +25,7 @@ option(NO_TEST "Do not build the tests" OFF) ################################################################################ find_package(RCMake 0.3 REQUIRED) find_package(RSys 0.6 REQUIRED) +find_package(StarSP 0.7 REQUIRED) find_package(SVX REQUIRED) find_package(HTCP REQUIRED) @@ -68,7 +69,7 @@ rcmake_prepend_path(HTRDR_FILES_INC ${HTRDR_SOURCE_DIR}) 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) +target_link_libraries(htrdr HTCP RSys SVX StarSP) if(CMAKE_COMPILER_IS_GNUCC) target_link_libraries(htrdr m) endif() diff --git a/src/htrdr.c b/src/htrdr.c @@ -17,6 +17,7 @@ #include "htrdr.h" #include "htrdr_args.h" +#include "htrdr_buffer.h" #include "htrdr_clouds.h" #include "htrdr_rectangle.h" #include "htrdr_solve.h" @@ -128,6 +129,7 @@ release_htrdr(ref_T* ref) ASSERT(htrdr); if(htrdr->svx) SVX(device_ref_put(htrdr->svx)); if(htrdr->clouds) SVX(tree_ref_put(htrdr->clouds)); + if(htrdr->buf) htrdr_buffer_ref_put(htrdr->buf); if(htrdr->rect) htrdr_rectangle_ref_put(htrdr->rect); logger_release(&htrdr->logger); MEM_RM(htrdr->allocator, htrdr); @@ -183,6 +185,28 @@ htrdr_create res = clouds_load(htrdr, args->verbose, args->input); if(res != RES_OK) goto error; + if(!args->dump_vtk) { /* Legacy mode */ + const size_t elmtsz = sizeof(double); + const size_t pitch = elmtsz * args->image.definition[0]; + + /* Create the image buffer */ + res = htrdr_buffer_create(htrdr, args->image.definition[0], + args->image.definition[1], pitch, elmtsz, 16, &htrdr->buf); + if(res != RES_OK) goto error; + + /* TODO check the validity of the parameters */ + htrdr->main_dir[0] = args->main_dir[0]; + htrdr->main_dir[1] = args->main_dir[1]; + htrdr->main_dir[2] = args->main_dir[2]; + htrdr->spp = args->image.spp; + + /* Create the plane on which the image buffer lies */ + 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; + + } + exit: *out_htrdr = htrdr; return res; diff --git a/src/htrdr.h b/src/htrdr.h @@ -20,16 +20,21 @@ #include <rsys/ref_count.h> /* Forward declaration */ +struct htrdr_args; +struct htrdr_buffer; +struct htrdr_rectangle; +struct mem_allocator; struct svx_device; struct svx_tree; -struct mem_allocator; -struct htrdr_args; struct htrdr { struct svx_device* svx; struct svx_tree* clouds; + struct htrdr_buffer* buf; struct htrdr_rectangle* rect; + double main_dir[3]; /* Main direction */ + size_t spp; /* #samples per pixel */ FILE* output; diff --git a/src/htrdr_args.h b/src/htrdr_args.h @@ -34,6 +34,8 @@ struct htrdr_args { unsigned spp; /* #samples per pixel */ } image; + double main_dir[3]; + int force_overwriting; int dump_vtk; /* Dump the loaded cloud properties in a VTK file */ int verbose; /* Verbosity level */ @@ -52,6 +54,7 @@ struct htrdr_args { {32, 32}, /* image definition */ \ 1 /* #samples per pixel */ \ }, \ + {0, 0, -1}, /* Main direction */ \ 0, /* Force overwriting */ \ 0, /* dump VTK */ \ 0, /* Verbose flag */ \ diff --git a/src/htrdr_buffer.c b/src/htrdr_buffer.c @@ -64,8 +64,35 @@ htrdr_buffer_create struct htrdr_buffer* buf = NULL; size_t memsz = 0; res_T res = RES_OK; - ASSERT(htrdr && width && height && pitch && elmt_size && align && out_buf); - ASSERT(IS_POW2(align) && width <= pitch); + ASSERT(htrdr && out_buf); + + if(!width || !height) { + htrdr_log_err(htrdr, "invalid buffer definition %lux%lu.\n", + (unsigned long)width, (unsigned long)height); + res = RES_BAD_ARG; + goto error; + } + if(pitch < width) { + htrdr_log_err(htrdr, + "invalid buffer pitch `%lu' wrt the buffer width `%lu'. " + "The buffer pitch cannot be less than the buffer width.\n", + (unsigned long)pitch, (unsigned long)width); + res = RES_BAD_ARG; + goto error; + } + if(elmtsz) { + htrdr_log_err(htrdr, + "the size of the buffer's elements cannot be null.\n"); + res = RES_BAD_ARG; + goto error; + } + if(!IS_POW2(align)) { + htrdr_log_err(htrdr, + "invalid buffer alignment `%lu'. It must be a power of 2.\n", + (unsigned long)align); + res = RES_BAD_ARG; + goto error; + } buf = MEM_CALLOC(htrdr->allocator, 1, sizeof(*buf)); if(!buf) { diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h @@ -28,5 +28,9 @@ htrdr_solve_transmission const double dir[3], double *val); +extern LOCAL_SYM res_T +htrdr_solve_transmission_buffer + (struct htrdr* htrdr); + #endif /* HTRDR_SOLVE_H */ diff --git a/src/htrdr_transmission.c b/src/htrdr_transmission.c @@ -14,8 +14,11 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "htrdr.h" +#include "htrdr_buffer.h" +#include "htrdr_rectangle.h" #include "htrdr_solve.h" +#include <star/ssp.h> #include <star/svx.h> #include <rsys/math.h> @@ -82,3 +85,63 @@ error: goto exit; } +res_T +htrdr_solve_transmission_buffer(struct htrdr* htrdr) +{ + struct htrdr_buffer_layout buf_layout; + struct ssp_rng* rng = NULL; + double pixsz[2]; + size_t ix, iy, ispp; + res_T res = RES_OK; + ASSERT(htrdr && htrdr->rect && htrdr->buf); + + htrdr_buffer_get_layout(htrdr->buf, &buf_layout); + + res = ssp_rng_create(htrdr->allocator, &ssp_rng_mt19937_64, &rng); + if(res != RES_OK) { + htrdr_log_err(htrdr, "could not allocate a RNG.\n"); + goto error; + } + + /* Compute the pixel size in the normalized image space */ + pixsz[0] = 1.0 / (double)buf_layout.width; + pixsz[1] = 1.0 / (double)buf_layout.height; + + FOR_EACH(iy, 0, buf_layout.height) { + const double y = (double)iy * pixsz[1]; /* Y in normalized image space */ + + FOR_EACH(ix, 0, buf_layout.width) { + const double x = (double)ix * pixsz[0]; /* X in normalized image space */ + double* val = htrdr_buffer_at(htrdr->buf, ix, iy); + double accum = 0; + + FOR_EACH(ispp, 0, htrdr->spp) { + double pixsamp[2]; /* Pixel sample */ + double pos[3]; /* World space position */ + double T; + + /* Uniformaly sample the pixel in the normalized image space */ + pixsamp[0] = x + ssp_rng_canonical(rng)*pixsz[0]; + pixsamp[1] = y + ssp_rng_canonical(rng)*pixsz[1]; + + /* Compute the world space position of the sample according to the + * image plane */ + htrdr_rectangle_sample_pos(htrdr->rect, pixsamp, pos); + + /* Solve the transmission for the sample position wrt the main dir */ + res = htrdr_solve_transmission(htrdr, pos, htrdr->main_dir, &T); + if(res != RES_OK) goto error; + + accum += T; + } + *val = accum / (double)htrdr->spp; + } + } + +exit: + if(rng) SSP(rng_ref_put(rng)); + return res; +error: + goto exit; +} +