htrdr

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

commit dddf1a2fc28a59662aebb33e05ed27d3cfe4f40d
parent 69297c12dc4074c2d6544df1a1ca4076cce1014c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 24 Jul 2018 11:54:25 +0200

Setup and use the camera

Replace the rectangle CLI option by a camera. Remove the
htrdr_solve_transmission function. Invoke the htrdr_draw_radiance_sw
function in legacy mode.

Diffstat:
Mcmake/CMakeLists.txt | 3+--
Mdoc/cli.txt | 7++++---
Msrc/htrdr.c | 59+++++++++++++++++++++++++++++------------------------------
Msrc/htrdr.h | 3+--
Msrc/htrdr_args.c | 124+++++++++++++++++++++++++++++++++++++++++++++----------------------------------
Msrc/htrdr_args.h | 12++++++++++++
Msrc/htrdr_camera.h | 2+-
Msrc/htrdr_solve.h | 4----
Dsrc/htrdr_transmission.c | 203-------------------------------------------------------------------------------
9 files changed, 119 insertions(+), 298 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -66,8 +66,7 @@ set(HTRDR_FILES_SRC htrdr_main.c htrdr_rectangle.c htrdr_sky.c - htrdr_sun.c - htrdr_transmission.c) + htrdr_sun.c) set(HTRDR_FILES_INC htrdr.h htrdr_c.h diff --git a/doc/cli.txt b/doc/cli.txt @@ -1,19 +1,20 @@ <image> ::= <imgopt>[:<imgopt> ... ] -<rectangle> ::= <rectopt>[:<rectopt> ... ] +<camera> ::= <camopt>[:<camopt> ... ] <imgopt> ::= <definition> | <#samples> <definition> ::= def=INTEGERxINTEGER # default 32x32 <#samples> ::= ssp=INTEGER # default 1 -<rectopt> ::= <target> +<camopt> ::= <target> | <position> | <up> - | <size> + | <fov> <target> ::= tgt=<float3> <position> ::= pos=<float3> <up> ::= up=<float3> <size> ::= sz=<float2> +<fov> ::= fov=FLOAT # In degrees. In [30,120] <float3> ::= FLOAT,FLOAT,FLOAT <float2> ::= FLOAT,FLOAT diff --git a/src/htrdr.c b/src/htrdr.c @@ -18,7 +18,7 @@ #include "htrdr.h" #include "htrdr_args.h" #include "htrdr_buffer.h" -#include "htrdr_rectangle.h" +#include "htrdr_camera.h" #include "htrdr_sky.h" #include "htrdr_sun.h" #include "htrdr_solve.h" @@ -218,6 +218,7 @@ htrdr_init const struct htrdr_args* args, struct htrdr* htrdr) { + double proj_ratio; res_T res = RES_OK; ASSERT(args && htrdr); @@ -233,6 +234,7 @@ htrdr_init htrdr->dump_vtk = args->dump_vtk; htrdr->verbose = args->verbose; htrdr->nthreads = MMIN(args->nthreads, (unsigned)omp_get_num_procs()); + htrdr->spp = args->image.spp; res = svx_device_create (&htrdr->logger, htrdr->allocator, args->verbose, &htrdr->svx); @@ -248,36 +250,26 @@ htrdr_init 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; - } + proj_ratio = + (double)args->image.definition[0] + / (double)args->image.definition[1]; + res = htrdr_camera_create(htrdr, args->camera.pos, args->camera.tgt, + args->camera.up, proj_ratio, args->camera.fov_x, &htrdr->cam); + if(res != RES_OK) goto error; - if(!args->output) { - htrdr->output = stdout; - } else { - res = open_output_stream(htrdr, args); - if(res != RES_OK) goto error; - } + res = htrdr_buffer_create(htrdr, + args->image.definition[0], /* Width */ + args->image.definition[1], /* Height */ + args->image.definition[0]*sizeof(struct htrdr_accum), /* Pitch */ + sizeof(struct htrdr_accum), /* Element size */ + 16, /* Alignement */ + &htrdr->buf); + if(res != RES_OK) goto error; res = htrdr_sun_create(htrdr, &htrdr->sun); if(res != RES_OK) goto error; + htrdr_sun_set_direction(htrdr->sun, args->main_dir); + res = htrdr_sky_create(htrdr, htrdr->sun, args->filename_les, args->filename_mie, &htrdr->sky); if(res != RES_OK) goto error; @@ -308,6 +300,13 @@ htrdr_init res = setup_geometry(htrdr, args->filename_obj); if(res != RES_OK) goto error; + if(!args->output) { + htrdr->output = stdout; + } else { + res = open_output_stream(htrdr, args); + if(res != RES_OK) goto error; + } + exit: return res; error: @@ -327,8 +326,8 @@ htrdr_release(struct htrdr* htrdr) if(htrdr->svx) SVX(device_ref_put(htrdr->svx)); if(htrdr->sky) htrdr_sky_ref_put(htrdr->sky); if(htrdr->sun) htrdr_sun_ref_put(htrdr->sun); + if(htrdr->cam) htrdr_camera_ref_put(htrdr->cam); if(htrdr->buf) htrdr_buffer_ref_put(htrdr->buf); - if(htrdr->rect) htrdr_rectangle_ref_put(htrdr->rect); logger_release(&htrdr->logger); } @@ -344,13 +343,13 @@ htrdr_run(struct htrdr* htrdr) char buf[128]; time_current(&t0); - res = htrdr_solve_transmission_buffer(htrdr); + res = htrdr_draw_radiance_sw(htrdr, htrdr->cam, htrdr->spp, htrdr->buf); if(res != RES_OK) goto error; time_sub(&t0, time_current(&t1), &t0); time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); htrdr_log(htrdr, "Elapsed time: %s\n", buf); - dump_buffer(htrdr); + /*dump_buffer(htrdr);*/ } exit: return res; diff --git a/src/htrdr.h b/src/htrdr.h @@ -45,9 +45,8 @@ struct htrdr { struct ssf_phase* phase_hg; /* Henyey & Greenstein phase function */ struct ssf_phase* phase_rayleigh; /* Rayleigh phase function */ + struct htrdr_camera* cam; 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.c b/src/htrdr_args.c @@ -34,9 +34,11 @@ print_help(const char* cmd) printf( " -c FILENAME path of the HTCP cloud properties file.\n"); printf( +" -C <camera> define the rendering point of view.\n"); + printf( " -d dump octree data to OUTPUT wrt the VTK ASCII file format.\n"); printf( -" -D X,Y,Z sun direction.\n"); +" -D X,Y,Z direction toward the sun.\n"); printf( " -f overwrite the OUTPUT file if it already exists.\n"); printf( @@ -44,15 +46,13 @@ print_help(const char* cmd) printf( " -h display this help and exit.\n"); printf( -" -I <image> define the image to compute.\n"); +" -i <image> define the image to compute.\n"); printf( " -m FILENAME path of the Mie data file.\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( " -t THREADS hint on the number of threads to use. By default use as\n" " many threads as CPU cores.\n"); printf( @@ -90,18 +90,38 @@ parse_definition(const char* str, unsigned val[2]) } static res_T -parse_rectangle_parameter(struct htrdr_args* args, const char* str) +parse_fov(const char* str, double* out_fov) +{ + double fov; + res_T res = RES_OK; + ASSERT(str && out_fov); + + res = cstr_to_double(str, &fov); + if(res != RES_OK) { + fprintf(stderr, "Invalid field of view `%s'.\n", str); + return RES_BAD_ARG; + } + if(fov < 30 || fov > 120) { + fprintf(stderr, "The field of view %g is not in [30, 120].\n", fov); + return RES_BAD_ARG; + } + *out_fov = fov; + return RES_OK; +} + +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; + res_T res = RES_OK; ASSERT(str && args); if(strlen(str) >= sizeof(buf) -1/*NULL char*/) { fprintf(stderr, - "Could not duplicate the rectangle option string `%s'.\n", str); + "Could not duplicate the image option string `%s'.\n", str); res = RES_MEM_ERR; goto error; } @@ -111,7 +131,7 @@ parse_rectangle_parameter(struct htrdr_args* args, const char* str) val = strtok_r(NULL, "", &ctx); if(!val) { - fprintf(stderr, "Missing a value to the rectangle option `%s'.\n", key); + fprintf(stderr, "Missing a value to the image option `%s'.\n", key); res = RES_BAD_ARG; goto error; } @@ -119,24 +139,31 @@ parse_rectangle_parameter(struct htrdr_args* args, const char* str) #define PARSE(Name, Func) \ res = Func; \ if(res != RES_OK) { \ - fprintf(stderr, "Invalid rectangle "Name" `%s'.\n", val); \ + fprintf(stderr, "Invalid image "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)); + if(!strcmp(key, "def")) { + PARSE("definition", parse_definition(val, args->image.definition)); + } else if(!strcmp(key, "spp")) { + PARSE("#samples per pixel", cstr_to_uint(val, &args->image.spp)); } else { - fprintf(stderr, "Invalid rectangle parameter `%s'.\n", key); + fprintf(stderr, "Invalid image parameter `%s'.\n", key); res = RES_BAD_ARG; goto error; } #undef PARSE + if(!args->image.definition[0] || !args->image.definition[1]) { + fprintf(stderr, "The image definition cannot be null.n"); + res = RES_BAD_ARG; + goto error; + } + if(!args->image.spp) { + fprintf(stderr, "The number of samples per pixel cannot be null.\n"); + res = RES_BAD_ARG; + goto error; + } + exit: return res; error: @@ -144,18 +171,17 @@ error: } static res_T -parse_image_parameter(struct htrdr_args* args, const char* str) +parse_camera_parameter(struct htrdr_args* args, const char* str) { char buf[128]; char* key; char* val; char* ctx; - res_T res; - ASSERT(str && args); + res_T res = RES_OK; if(strlen(str) >= sizeof(buf) -1/*NULL char*/) { fprintf(stderr, - "Could not duplicate the image option string `%s'.\n", str); + "Could not duplicate the rectangle option string `%s'.\n", str); res = RES_MEM_ERR; goto error; } @@ -165,39 +191,31 @@ parse_image_parameter(struct htrdr_args* args, const char* str) val = strtok_r(NULL, "", &ctx); if(!val) { - fprintf(stderr, "Missing a value to the image option `%s'.\n", key); + fprintf(stderr, "Missing value to the camera 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); \ + #define PARSE(Name, Func) { \ + if(RES_OK != (res = Func)) { \ + fprintf(stderr, "Invalid camera "Name" `%s'.\n", val); \ goto error; \ - } (void)0 - if(!strcmp(key, "def")) { - PARSE("definition", parse_definition(val, args->image.definition)); - } else if(!strcmp(key, "spp")) { - PARSE("#samples per pixel", cstr_to_uint(val, &args->image.spp)); + } \ + } (void)0 + if(!strcmp(key, "pos")) { + PARSE("position", parse_doubleX(val, args->camera.pos, 3)); + } else if(!strcmp(key, "tgt")) { + PARSE("target", parse_doubleX(val, args->camera.tgt, 3)); + } else if(!strcmp(key, "up")) { + PARSE("up vector", parse_doubleX(val, args->camera.up, 3)); + } else if(!strcmp(key, "fov")) { + PARSE("field-of-view", parse_fov(val, &args->camera.fov_x)); } else { - fprintf(stderr, "Invalid image parameter `%s'.\n", key); + fprintf(stderr, "Invalid camera parameter `%s'.\n", key); res = RES_BAD_ARG; goto error; } #undef PARSE - - if(!args->image.definition[0] || !args->image.definition[1]) { - fprintf(stderr, "The image definition cannot be null.n"); - res = RES_BAD_ARG; - goto error; - } - if(!args->image.spp) { - fprintf(stderr, "The number of samples per pixel cannot be null.\n"); - res = RES_BAD_ARG; - goto error; - } - exit: return res; error: @@ -274,8 +292,12 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv) *args = HTRDR_ARGS_DEFAULT; - while((opt = getopt(argc, argv, "c:D:dfhI:m:o:r:t:v")) != -1) { + while((opt = getopt(argc, argv, "C:c:D:dfhi:m:o:t:v")) != -1) { switch(opt) { + case 'C': + res = parse_multiple_parameters + (args, optarg, parse_camera_parameter); + break; case 'c': args->filename_les = optarg; break; case 'D': res = parse_sun_dir(args, optarg); break; case 'd': args->dump_vtk = 1; break; @@ -286,16 +308,12 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv) htrdr_args_release(args); args->quit = 1; goto exit; - case 'I': + case 'i': res = parse_multiple_parameters (args, optarg, parse_image_parameter); break; case 'm': args->filename_mie = optarg; break; case 'o': args->output = optarg; break; - case 'r': - res = parse_multiple_parameters - (args, optarg, parse_rectangle_parameter); - 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; @@ -312,13 +330,13 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv) } } if(!args->filename_les) { - fprintf(stderr, + fprintf(stderr, "Missing the path of the cloud properties file -- option '-c'\n"); res = RES_BAD_ARG; goto error; } if(!args->filename_mie) { - fprintf(stderr, + fprintf(stderr, "Missing the path toward the file of the Mie's data -- option '-m'\n"); res = RES_BAD_ARG; goto error; diff --git a/src/htrdr_args.h b/src/htrdr_args.h @@ -32,6 +32,13 @@ struct htrdr_args { } rectangle; struct { + double pos[3]; + double tgt[3]; + double up[3]; + double fov_x; /* In degrees */ + } camera; + + struct { unsigned definition[2]; /* #pixels in X and Y */ unsigned spp; /* #samples per pixel */ } image; @@ -56,6 +63,11 @@ struct htrdr_args { {0, 1, 0}, /* plane up */ \ {1, 1}, /* plane size */ \ }, { \ + {0, 0, 0}, /* Camera position */ \ + {0, 0, 1}, /* Camera target */ \ + {0, 1, 0}, /* Camera up */ \ + 70, /* Horizontal field of view */ \ + }, { \ {32, 32}, /* image definition */ \ 1 /* #samples per pixel */ \ }, \ diff --git a/src/htrdr_camera.h b/src/htrdr_camera.h @@ -28,7 +28,7 @@ htrdr_camera_create const double position[3], const double target[3], const double up[3], - const double proj_ratio, + const double proj_ratio, /* Width / Height */ const double fov, /* In radian */ struct htrdr_camera** cam); diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h @@ -33,10 +33,6 @@ struct htrdr; struct htrdr_camera; struct ssp_rng; -extern LOCAL_SYM res_T -htrdr_solve_transmission_buffer - (struct htrdr* htrdr); - extern LOCAL_SYM double htrdr_compute_radiance_sw (struct htrdr* htrdr, diff --git a/src/htrdr_transmission.c b/src/htrdr_transmission.c @@ -1,203 +0,0 @@ -/* 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_buffer.h" -#include "htrdr_rectangle.h" -#include "htrdr_sky.h" -#include "htrdr_solve.h" - -#include <star/ssp.h> -#include <star/svx.h> -#include <rsys/math.h> - -#include <omp.h> - -struct transmit_context { - const struct htrdr_sky* sky; - double tau_sampled; - double tau_max_min; - double tau_min; -}; -static const struct transmit_context TRANSMIT_CONTEXT_NULL = {NULL,0,0,0}; - -/******************************************************************************* - * Helper functions - ******************************************************************************/ -static FINLINE uint16_t -morton2D_decode(const uint32_t u32) -{ - uint32_t x = u32 & 0x55555555; - x = (x | (x >> 1)) & 0x33333333; - x = (x | (x >> 2)) & 0x0F0F0F0F; - x = (x | (x >> 4)) & 0x00FF00FF; - x = (x | (x >> 8)) & 0x0000FFFF; - return (uint16_t)x; -} - -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_ext_min; - double k_ext_max; - ASSERT(hit && ctx && !SVX_HIT_NONE(hit)); - (void)org, (void)dir, (void)range; - - k_ext_min = htrdr_sky_fetch_svx_voxel_property(ctx->sky, HTRDR_Kext, - HTRDR_SVX_MIN, HTRDR_ALL_COMPONENTS, -1/*FIXME*/, &hit->voxel); - k_ext_max = htrdr_sky_fetch_svx_voxel_property(ctx->sky, HTRDR_Kext, - HTRDR_SVX_MAX, HTRDR_ALL_COMPONENTS, -1/*FIXME*/, &hit->voxel); - - dst = hit->distance[1] - hit->distance[0]; - ASSERT(dst >= 0); - ctx->tau_min += k_ext_min*dst; - ctx->tau_max_min += (k_ext_max - k_ext_min)*dst; - return ctx->tau_max_min < ctx->tau_sampled; -} - -static res_T -transmission_realisation - (struct htrdr* htrdr, - struct ssp_rng* rng, - const double pos[3], - const double dir[3], - double *val) -{ - struct svx_hit hit = SVX_HIT_NULL; - struct svx_tree* svx_tree = NULL; - struct transmit_context ctx = TRANSMIT_CONTEXT_NULL; - const double range[2] = {0, INF}; - res_T res = RES_OK; - ASSERT(htrdr && pos && dir && val); - - ctx.tau_sampled = ssp_ran_exp(rng, 1.0); - ctx.sky = htrdr->sky; - svx_tree = htrdr_sky_get_svx_tree(htrdr->sky); - res = svx_octree_trace_ray(svx_tree, 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; - } - if(!SVX_HIT_NONE(&hit)) { - *val = 0; - } else { - *val = ctx.tau_min ? exp(-ctx.tau_min) : 1.0; - } - -exit: - return res; -error: - goto exit; -} - -/******************************************************************************* - * Local functions - ******************************************************************************/ -res_T -htrdr_solve_transmission_buffer(struct htrdr* htrdr) -{ - struct htrdr_buffer_layout buf_layout; - struct ssp_rng* rng = NULL; - double pixsz[2]; - int64_t mcode_max; - int64_t mcode; - ATOMIC 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; - - /* Compute the maximum morton code */ - mcode_max = (int64_t)round_up_pow2(MMAX(buf_layout.height, buf_layout.width)); - mcode_max = mcode_max*mcode_max; - - /* Setup the number of threads */ - omp_set_num_threads((int)htrdr->nthreads); - - #pragma omp parallel for schedule(dynamic, 32/*chunck size*/) - for(mcode=0; mcode<mcode_max; ++mcode) { - size_t ispp; - double* val; - double x, y; - double accum; - uint16_t ix, iy; - res_T res_local = RES_OK; - - if(ATOMIC_GET(&res) != RES_OK) continue; - - ix = morton2D_decode((uint32_t)(mcode>>0)); - if(ix > buf_layout.width) continue; - iy = morton2D_decode((uint32_t)(mcode>>1)); - if(iy > buf_layout.height) continue; - - val = htrdr_buffer_at(htrdr->buf, ix, iy); - - /* Define lower left the pixel coordinate in the normalized image space */ - x = (double)ix*pixsz[0]; - y = (double)iy*pixsz[1]; - - 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_local = transmission_realisation(htrdr, rng, pos, htrdr->main_dir, &T); - if(res_local != RES_OK) { - ATOMIC_SET(&res, res_local); - continue; - } - - accum += T; - } - *val = accum / (double)htrdr->spp; - } - -exit: - if(rng) SSP(rng_ref_put(rng)); - return (res_T)res; -error: - goto exit; -} -