htrdr

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

commit f98a86907d867a2f96f52b49506ec7998917c7fb
parent ee7c9e6439f70ec2a427f9f2a2321cc6a6fa3379
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 19 Oct 2018 15:53:36 +0200

Merge branch 'feature_ground'

Diffstat:
Mcmake/CMakeLists.txt | 6+++++-
Msrc/htrdr.c | 9++++++---
Msrc/htrdr.h | 3+--
Msrc/htrdr_args.c | 17++++++++++-------
Msrc/htrdr_args.h | 4+++-
Msrc/htrdr_c.h | 4++++
Msrc/htrdr_compute_radiance_sw.c | 143++++++++++++++++++++++---------------------------------------------------------
Asrc/htrdr_ground.c | 333+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/htrdr_ground.h | 51+++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/htrdr_sky.c | 126+++++++++++++++++++++++++++++--------------------------------------------------
Msrc/htrdr_sky.h | 2+-
Asrc/htrdr_slab.c | 133+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/htrdr_slab.h | 46++++++++++++++++++++++++++++++++++++++++++++++
13 files changed, 678 insertions(+), 199 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -24,7 +24,7 @@ option(NO_TEST "Do not build the tests" OFF) # Check dependencies ################################################################################ find_package(RCMake 0.3 REQUIRED) -find_package(RSys 0.6 REQUIRED) +find_package(RSys 0.7 REQUIRED) find_package(Star3D 0.5 REQUIRED) find_package(Star3DAW 0.3.1 REQUIRED) find_package(StarSF 0.6 REQUIRED) @@ -66,8 +66,10 @@ set(HTRDR_FILES_SRC htrdr_compute_radiance_sw.c htrdr_draw_radiance_sw.c htrdr_grid.c + htrdr_ground.c htrdr_main.c htrdr_rectangle.c + htrdr_slab.c htrdr_sky.c htrdr_sun.c) set(HTRDR_FILES_INC @@ -77,7 +79,9 @@ set(HTRDR_FILES_INC htrdr_buffer.h htrdr_camera.h htrdr_grid.h + htrdr_ground.h htrdr_rectangle.h + htrdr_slab.h htrdr_sky.h htrdr_sun.h htrdr_solve.h) diff --git a/src/htrdr.c b/src/htrdr.c @@ -20,6 +20,7 @@ #include "htrdr_args.h" #include "htrdr_buffer.h" #include "htrdr_camera.h" +#include "htrdr_ground.h" #include "htrdr_sky.h" #include "htrdr_sun.h" #include "htrdr_solve.h" @@ -29,7 +30,6 @@ #include <rsys/str.h> #include <star/s3d.h> -#include <star/s3daw.h> #include <star/ssf.h> #include <star/svx.h> @@ -149,6 +149,7 @@ error: goto exit; } +#if 0 static res_T setup_geometry(struct htrdr* htrdr, const char* filename) { @@ -210,6 +211,7 @@ exit: error: goto exit; } +#endif static res_T open_file_stamp @@ -370,7 +372,8 @@ htrdr_init goto error; } - res = setup_geometry(htrdr, args->filename_obj); + res = htrdr_ground_create + (htrdr, args->filename_obj, args->repeat_ground, &htrdr->ground); if(res != RES_OK) goto error; proj_ratio = @@ -431,9 +434,9 @@ void htrdr_release(struct htrdr* htrdr) { ASSERT(htrdr); - if(htrdr->s3d_scn_view) S3D(scene_view_ref_put(htrdr->s3d_scn_view)); if(htrdr->s3d) S3D(device_ref_put(htrdr->s3d)); if(htrdr->svx) SVX(device_ref_put(htrdr->svx)); + if(htrdr->ground) htrdr_ground_ref_put(htrdr->ground); 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); diff --git a/src/htrdr.h b/src/htrdr.h @@ -45,8 +45,7 @@ struct htrdr { struct svx_device* svx; struct s3d_device* s3d; - struct s3d_scene_view* s3d_scn_view; - + struct htrdr_ground* ground; struct htrdr_sky* sky; struct htrdr_sun* sun; diff --git a/src/htrdr_args.c b/src/htrdr_args.c @@ -38,22 +38,24 @@ print_help(const char* cmd) 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 AZIMUTH,ELEVATION\n" " sun direction in degrees. Following the right-handed\n" " convention, the azimuthal rotation is counter-clockwise\n" " around the Z axis, with 0 aligned on the X axis. The\n" " elevation rotation starts from 0 up to 90 at zenith.\n"); printf( +" -d dump octree data to OUTPUT wrt the VTK ASCII file format.\n"); + printf( " -f overwrite the OUTPUT file if it already exists.\n"); printf( -" -g FILENAME path of the OBJ geometry file.\n"); +" -g FILENAME path of an OBJ file representing the ground geometry.\n"); printf( " -h display this help and exit.\n"); printf( " -i <image> define the image to compute.\n"); printf( +" -R infinitely repeat the ground along the X and Y axis.\n"); + printf( " -r infinitely repeat the clouds along the X and Y axis.\n"); printf( " -m FILENAME path of the Mie data file.\n"); @@ -61,13 +63,13 @@ print_help(const char* cmd) " -o OUTPUT file where data are written. If not defined, data are\n" " written to standard output.\n"); printf( -" -t THREADS hint on the number of threads to use. By default use as\n" -" many threads as CPU cores.\n"); - printf( " -T THRESHOLD optical thickness used as threshold during the octree\n" " building. By default its value is `%g'.\n", HTRDR_ARGS_DEFAULT.optical_thickness); printf( +" -t THREADS hint on the number of threads to use. By default use as\n" +" many threads as CPU cores.\n"); + printf( " -v make the program more verbose.\n"); printf("\n"); printf( @@ -318,7 +320,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv) *args = HTRDR_ARGS_DEFAULT; - while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:m:o:rT:t:v")) != -1) { + while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:m:o:RrT:t:v")) != -1) { switch(opt) { case 'a': args->filename_gas = optarg; break; case 'C': @@ -342,6 +344,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv) case 'm': args->filename_mie = optarg; break; case 'o': args->output = optarg; break; case 'r': args->repeat_clouds = 1; break; + case 'R': args->repeat_ground = 1; break; case 'T': res = cstr_to_double(optarg, &args->optical_thickness); if(res == RES_OK && args->optical_thickness < 0) res = RES_BAD_ARG; diff --git a/src/htrdr_args.h b/src/htrdr_args.h @@ -52,7 +52,8 @@ struct htrdr_args { int force_overwriting; int dump_vtk; /* Dump the loaded cloud properties in a VTK file */ int verbose; /* Verbosity level */ - int repeat_clouds; /* Make clouds infinite in X and Y */ + int repeat_clouds; /* Make the clouds infinite in X and Y */ + int repeat_ground; /* Make the ground infinite in X and Y */ int quit; /* Quit the application */ }; @@ -84,6 +85,7 @@ struct htrdr_args { 0, /* dump VTK */ \ 0, /* Verbose flag */ \ 0, /* Repeat clouds */ \ + 0, /* Repeat ground */ \ 0 /* Quit the application */ \ } static const struct htrdr_args HTRDR_ARGS_DEFAULT = HTRDR_ARGS_DEFAULT__; diff --git a/src/htrdr_c.h b/src/htrdr_c.h @@ -16,6 +16,10 @@ #ifndef HTRDR_C_H #define HTRDR_C_H +#include <rsys/rsys.h> + +struct htrdr; + #define SW_WAVELENGTH_MIN 380 /* In nanometer */ #define SW_WAVELENGTH_MAX 780 /* In nanometer */ diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c @@ -15,6 +15,7 @@ #include "htrdr.h" #include "htrdr_c.h" +#include "htrdr_ground.h" #include "htrdr_solve.h" #include "htrdr_sky.h" #include "htrdr_sun.h" @@ -28,14 +29,6 @@ #include <rsys/float2.h> #include <rsys/float3.h> -struct ray_context { - float range[2]; - struct s3d_hit hit_prev; -}; -static const struct ray_context RAY_CONTEXT_NULL = { - {0, INF}, S3D_HIT_NULL__ -}; - struct scattering_context { struct ssp_rng* rng; const struct htrdr_sky* sky; @@ -68,24 +61,6 @@ static const struct transmissivity_context TRANSMISSION_CONTEXT_NULL = { /******************************************************************************* * Helper functions ******************************************************************************/ -/* Check that `hit' roughly lies on an edge. For triangular primitives, a - * simple but approximative way is to test that its position have at least one - * barycentric coordinate roughly equal to 0 or 1. */ -static FINLINE int -hit_on_edge(const struct s3d_hit* hit) -{ - const float on_edge_eps = 1.e-4f; - float w; - ASSERT(hit && !S3D_HIT_NONE(hit)); - w = 1.f - hit->uv[0] - hit->uv[1]; - return eq_epsf(hit->uv[0], 0.f, on_edge_eps) - || eq_epsf(hit->uv[0], 1.f, on_edge_eps) - || eq_epsf(hit->uv[1], 0.f, on_edge_eps) - || eq_epsf(hit->uv[1], 1.f, on_edge_eps) - || eq_epsf(w, 0.f, on_edge_eps) - || eq_epsf(w, 1.f, on_edge_eps); -} - static int scattering_hit_filter (const struct svx_hit* hit, @@ -95,7 +70,6 @@ scattering_hit_filter void* context) { struct scattering_context* ctx = context; - double vox_dst; /* Distance to traverse into the voxel */ double ks_max; int pursue_traversal = 1; ASSERT(hit && ctx && !SVX_HIT_NONE(hit) && org && dir && range); @@ -104,13 +78,14 @@ scattering_hit_filter ks_max = htrdr_sky_fetch_svx_voxel_property(ctx->sky, HTRDR_Ks, HTRDR_SVX_MAX, HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, &hit->voxel); - /* Compute the distance to traverse into the voxel */ - vox_dst = hit->distance[1] - hit->distance[0]; + ctx->traversal_dst = hit->distance[0]; /* Iterate until a collision occurs into the voxel or until the ray * does not collide the voxel */ for(;;) { - const double T = vox_dst * ks_max; /* Compute tau for the current leaf */ + /* Compute tau for the current leaf */ + const double vox_dst = hit->distance[1] - ctx->traversal_dst; + const double T = vox_dst * ks_max; /* A collision occurs behind `vox_dst' */ if(ctx->Ts > T) { @@ -127,7 +102,7 @@ scattering_hit_filter const double collision_dst = ctx->Ts / ks_max; /* Compute the traversed distance up to the challenged collision */ - ctx->traversal_dst = hit->distance[0] + collision_dst; + ctx->traversal_dst += collision_dst; ASSERT(ctx->traversal_dst >= hit->distance[0]); ASSERT(ctx->traversal_dst <= hit->distance[1]); @@ -147,8 +122,6 @@ scattering_hit_filter break; } else { /* Null collision */ ctx->Ts = ssp_ran_exp(ctx->rng, 1); /* Sample a new optical thickness */ - /* Compute the remaining distance to traverse into the voxel */ - vox_dst = hit->distance[1] - ctx->traversal_dst; } } } @@ -165,7 +138,6 @@ transmissivity_hit_filter { struct transmissivity_context* ctx = context; int comp_mask = HTRDR_ALL_COMPONENTS; - double vox_dst; /* Distance to traverse into the voxel */ double k_max; double k_min; int pursue_traversal = 1; @@ -178,13 +150,13 @@ transmissivity_hit_filter HTRDR_SVX_MAX, comp_mask, ctx->iband, ctx->iquad, &hit->voxel); ASSERT(k_min <= k_max); - /* Compute the distance to traverse into the voxel */ - vox_dst = hit->distance[1] - hit->distance[0]; - ctx->Tmin += vox_dst * k_min; + ctx->Tmin += (hit->distance[1] - hit->distance[0]) * k_min; + ctx->traversal_dst = hit->distance[0]; /* Iterate until a collision occurs into the voxel or until the ray * does not collide the voxel */ for(;;) { + const double vox_dst = hit->distance[1] - ctx->traversal_dst; const double Tdif = vox_dst * (k_max-k_min); /* A collision occurs behind `vox_dst' */ @@ -202,7 +174,7 @@ transmissivity_hit_filter double collision_dst = ctx->Ts / (k_max - k_min); /* Compute the traversed distance up to the challenged collision */ - ctx->traversal_dst = hit->distance[0] + collision_dst; + ctx->traversal_dst += collision_dst; ASSERT(ctx->traversal_dst >= hit->distance[0]); ASSERT(ctx->traversal_dst <= hit->distance[1]); @@ -222,8 +194,6 @@ transmissivity_hit_filter break; } else { /* Null collision */ ctx->Ts = ssp_ran_exp(ctx->rng, 1); /* Sample a new optical thickness */ - /* Compute the remaining distance to traverse into the voxel */ - vox_dst = hit->distance[1] - ctx->traversal_dst; } } } @@ -239,30 +209,13 @@ transmissivity const size_t iquad, const double pos[3], const double dir[3], - const double range[2], - const struct s3d_hit* hit_prev) /* May be NULL */ + const double range[2]) { - struct s3d_hit s3d_hit; struct svx_hit svx_hit; struct transmissivity_context transmissivity_ctx = TRANSMISSION_CONTEXT_NULL; - struct ray_context rctx = RAY_CONTEXT_NULL; - float ray_pos[3]; - float ray_dir[3]; ASSERT(htrdr && rng && pos && dir && range); - /* Find the first intersection with a surface */ - f3_set_d3(ray_pos, pos); - f3_set_d3(ray_dir, dir); - f2_set_d2(rctx.range, range); - rctx.hit_prev = hit_prev ? *hit_prev : S3D_HIT_NULL; - S3D(scene_view_trace_ray(htrdr->s3d_scn_view, ray_pos, ray_dir, rctx.range, - &rctx, &s3d_hit)); - - if(!S3D_HIT_NONE(&s3d_hit)) { - return 0; - } - transmissivity_ctx.rng = rng; transmissivity_ctx.sky = htrdr->sky; transmissivity_ctx.iband = iband; @@ -295,8 +248,9 @@ htrdr_compute_radiance_sw const size_t iquad) { struct s3d_hit s3d_hit = S3D_HIT_NULL; - struct svx_hit svx_hit = SVX_HIT_NULL; + struct s3d_hit s3d_hit_tmp = S3D_HIT_NULL; struct s3d_hit s3d_hit_prev = S3D_HIT_NULL; + struct svx_hit svx_hit = SVX_HIT_NULL; struct ssf_phase* phase_hg = NULL; struct ssf_phase* phase_rayleigh = NULL; struct ssf_bsdf* bsdf = NULL; @@ -322,9 +276,6 @@ htrdr_compute_radiance_sw double w = 0; /* MC weight */ double g = 0; /* Asymmetry parameter of the HG phase function */ - float ray_pos[3]; - float ray_dir[3]; - ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads); CHK(RES_OK == ssf_bsdf_create @@ -358,25 +309,26 @@ htrdr_compute_radiance_sw if(htrdr_sun_is_dir_in_solar_cone(htrdr->sun, dir)) { /* Add the direct contribution of the sun */ d2(range, 0, FLT_MAX); - Tr = transmissivity - (htrdr, rng, HTRDR_Kext, iband, iquad , pos, dir, range, NULL); - w = L_sun * Tr; + + /* Check that the ray is not occlude along the submitted range */ + HTRDR(ground_trace_ray(htrdr->ground, pos, dir, range, NULL, &s3d_hit_tmp)); + if(!S3D_HIT_NONE(&s3d_hit_tmp)) { + Tr = 0; + } else { + Tr = transmissivity + (htrdr, rng, HTRDR_Kext, iband, iquad , pos, dir, range); + w = L_sun * Tr; + } } /* Radiative random walk */ for(;;) { struct scattering_context scattering_ctx = SCATTERING_CONTEXT_NULL; - struct ray_context rctx = RAY_CONTEXT_NULL; - - /* Setup the ray to trace */ - f3_set_d3(ray_pos, pos); - f3_set_d3(ray_dir, dir); - f2(rctx.range, 0, FLT_MAX); - rctx.hit_prev = s3d_hit_prev; /* Find the first intersection with a surface */ - S3D(scene_view_trace_ray(htrdr->s3d_scn_view, ray_pos, ray_dir, rctx.range, - &rctx, &s3d_hit)); + d2(range, 0, DBL_MAX); + HTRDR(ground_trace_ray + (htrdr->ground, pos, dir, range, &s3d_hit_prev, &s3d_hit)); /* Sample an optical thickness */ scattering_ctx.Ts = ssp_ran_exp(rng, 1); @@ -400,7 +352,7 @@ htrdr_compute_radiance_sw || ( svx_hit.distance[0] <= scattering_ctx.traversal_dst && svx_hit.distance[1] >= scattering_ctx.traversal_dst)); - /* Negative the incoming dir to match the convention of the SSF library */ + /* Negate the incoming dir to match the convention of the SSF library */ d3_minus(wo, dir); /* Compute the new position */ @@ -412,7 +364,7 @@ htrdr_compute_radiance_sw * next position */ d2(range, 0, scattering_ctx.traversal_dst); Tr_abs = transmissivity - (htrdr, rng, HTRDR_Ka, iband, iquad, pos, dir, range, &s3d_hit_prev); + (htrdr, rng, HTRDR_Ka, iband, iquad, pos, dir, range); if(Tr_abs <= 0) break; /* Sample a sun direction */ @@ -422,8 +374,16 @@ htrdr_compute_radiance_sw s3d_hit_prev = SVX_HIT_NONE(&svx_hit) ? s3d_hit : S3D_HIT_NULL; /* Check that the sun is visible from the new position */ - Tr = transmissivity(htrdr, rng, HTRDR_Kext, iband, iquad, pos_next, - sun_dir, range, &s3d_hit_prev); + HTRDR(ground_trace_ray + (htrdr->ground, pos_next, sun_dir, range, &s3d_hit_prev, &s3d_hit_tmp)); + + /* Compute the sun transmissivity */ + if(!S3D_HIT_NONE(&s3d_hit_tmp)) { + Tr = 0; + } else { + Tr = transmissivity + (htrdr, rng, HTRDR_Kext, iband, iquad, pos_next, sun_dir, range); + } /* Scattering at a surface */ if(SVX_HIT_NONE(&svx_hit)) { @@ -440,7 +400,7 @@ htrdr_compute_radiance_sw if(d3_dot(N, sun_dir) < 0) { /* Below the ground */ R = 0; } else { - R = ssf_bsdf_eval(bsdf, wo, N, sun_dir); + R = ssf_bsdf_eval(bsdf, wo, N, sun_dir) * d3_dot(N, sun_dir); } /* Scattering in the medium */ @@ -481,28 +441,3 @@ htrdr_compute_radiance_sw return w; } -int -htrdr_ground_filter - (const struct s3d_hit* hit, - const float ray_org[3], - const float ray_dir[3], - void* ray_data, - void* filter_data) -{ - const struct ray_context* ray_ctx = ray_data; - (void)ray_org, (void)ray_dir, (void)filter_data; - - if(!ray_ctx) return 0; - - if(S3D_PRIMITIVE_EQ(&hit->prim, &ray_ctx->hit_prev.prim)) return 1; - - if(!S3D_HIT_NONE(&ray_ctx->hit_prev) && eq_epsf(hit->distance, 0, 1.e-1f)) { - /* If the targeted point is near of the origin, check that it lies on an - * edge/vertex shared by the 2 primitives. */ - return hit_on_edge(&ray_ctx->hit_prev) && hit_on_edge(hit); - } - - return hit->distance <= ray_ctx->range[0] - || hit->distance >= ray_ctx->range[1]; -} - diff --git a/src/htrdr_ground.c b/src/htrdr_ground.c @@ -0,0 +1,333 @@ +/* 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_ground.h" +#include "htrdr_slab.h" + +#include <rsys/clock_time.h> +#include <rsys/cstr.h> +#include <rsys/double2.h> +#include <rsys/double3.h> +#include <rsys/float2.h> +#include <rsys/float3.h> + +#include <star/s3d.h> +#include <star/s3daw.h> + +struct ray_context { + float range[2]; + struct s3d_hit hit_prev; + int id[2]; +}; +#define RAY_CONTEXT_NULL__ {{0,INF}, S3D_HIT_NULL__, {0,0}} +static const struct ray_context RAY_CONTEXT_NULL = RAY_CONTEXT_NULL__; + +struct trace_ground_context { + struct s3d_scene_view* view; + struct ray_context context; + struct s3d_hit* hit; +}; +static const struct trace_ground_context TRACE_GROUND_CONTEXT_NULL = { + NULL, RAY_CONTEXT_NULL__, NULL +}; + +struct htrdr_ground { + struct s3d_scene_view* view; + float lower[3]; /* Ground lower bound */ + float upper[3]; /* Ground upper bound */ + int repeat; /* Make the ground infinite in X and Y */ + + struct htrdr* htrdr; + ref_T ref; +}; + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +/* Check that `hit' roughly lies on an edge. For triangular primitives, a + * simple but approximative way is to test that its position have at least one + * barycentric coordinate roughly equal to 0 or 1. */ +static FINLINE int +hit_on_edge(const struct s3d_hit* hit) +{ + const float on_edge_eps = 1.e-4f; + float w; + ASSERT(hit && !S3D_HIT_NONE(hit)); + w = 1.f - hit->uv[0] - hit->uv[1]; + return eq_epsf(hit->uv[0], 0.f, on_edge_eps) + || eq_epsf(hit->uv[0], 1.f, on_edge_eps) + || eq_epsf(hit->uv[1], 0.f, on_edge_eps) + || eq_epsf(hit->uv[1], 1.f, on_edge_eps) + || eq_epsf(w, 0.f, on_edge_eps) + || eq_epsf(w, 1.f, on_edge_eps); +} + +static int +ground_filter + (const struct s3d_hit* hit, + const float ray_org[3], + const float ray_dir[3], + void* ray_data, + void* filter_data) +{ + const struct ray_context* ray_ctx = ray_data; + (void)ray_org, (void)ray_dir, (void)filter_data; + + if(!ray_ctx) return 0; + + if(S3D_PRIMITIVE_EQ(&hit->prim, &ray_ctx->hit_prev.prim)) return 1; + + if(!S3D_HIT_NONE(&ray_ctx->hit_prev) && eq_epsf(hit->distance, 0, 1.e-1f)) { + /* If the targeted point is near of the origin, check that it lies on an + * edge/vertex shared by the 2 primitives. */ + return hit_on_edge(&ray_ctx->hit_prev) && hit_on_edge(hit); + } + + return hit->distance <= ray_ctx->range[0] + || hit->distance >= ray_ctx->range[1]; +} + +static INLINE res_T +trace_ground + (const double org[3], + const double dir[3], + const double range[2], + void* context, + int* hit) +{ + struct trace_ground_context* ctx = context; + float ray_org[3]; + float ray_dir[3]; + float ray_range[2]; + res_T res = RES_OK; + ASSERT(org && dir && range && context && hit); + + f3_set_d3(ray_org, org); + f3_set_d3(ray_dir, dir); + f2_set_d2(ray_range, range); + + res = s3d_scene_view_trace_ray + (ctx->view, ray_org, ray_dir, ray_range, &ctx->context, ctx->hit); + if(res != RES_OK) return res; + + *hit = !S3D_HIT_NONE(ctx->hit); + return RES_OK; +} + +static res_T +setup_ground(struct htrdr_ground* ground, const char* obj_filename) +{ + struct s3d_scene* scn = NULL; + struct s3daw* s3daw = NULL; + size_t nshapes; + size_t ishape; + res_T res = RES_OK; + ASSERT(ground && obj_filename); + + res = s3daw_create(&ground->htrdr->logger, ground->htrdr->allocator, NULL, + NULL, ground->htrdr->s3d, ground->htrdr->verbose, &s3daw); + if(res != RES_OK) { + htrdr_log_err(ground->htrdr, + "%s: could not create the Star-3DAW device -- %s.\n", + FUNC_NAME, res_to_cstr(res)); + goto error; + } + + res = s3daw_load(s3daw, obj_filename); + if(res != RES_OK) { + htrdr_log_err(ground->htrdr, + "%s: could not load the obj file `%s' -- %s.\n", + FUNC_NAME, obj_filename, res_to_cstr(res)); + goto error; + } + + res = s3d_scene_create(ground->htrdr->s3d, &scn); + if(res != RES_OK) { + htrdr_log_err(ground->htrdr, + "%s: could not create the Star-3D scene of the ground -- %s.\n", + FUNC_NAME, res_to_cstr(res)); + goto error; + } + + S3DAW(get_shapes_count(s3daw, &nshapes)); + FOR_EACH(ishape, 0, nshapes) { + struct s3d_shape* shape; + S3DAW(get_shape(s3daw, ishape, &shape)); + res = s3d_mesh_set_hit_filter_function(shape, ground_filter, NULL); + if(res != RES_OK) { + htrdr_log_err(ground->htrdr, + "%s: could not setup the hit filter function of the ground geometry " + "-- %s.\n", FUNC_NAME, res_to_cstr(res)); + goto error; + } + res = s3d_scene_attach_shape(scn, shape); + if(res != RES_OK) { + htrdr_log_err(ground->htrdr, + "%s: could not attach the ground geometry to its Star-3D scene -- %s.\n", + FUNC_NAME, res_to_cstr(res)); + goto error; + } + } + + res = s3d_scene_view_create(scn, S3D_TRACE, &ground->view); + if(res != RES_OK) { + htrdr_log_err(ground->htrdr, + "%s: could not create the Star-3D scene view of the ground geometry " + "-- %s.\n", FUNC_NAME, res_to_cstr(res)); + goto error; + } + + res = s3d_scene_view_get_aabb(ground->view, ground->lower, ground->upper); + if(res != RES_OK) { + htrdr_log_err(ground->htrdr, + "%s: could not get the ground bounding box -- %s.\n", + FUNC_NAME, res_to_cstr(res)); + goto error; + } + +exit: + if(s3daw) S3DAW(ref_put(s3daw)); + if(scn) S3D(scene_ref_put(scn)); + return res; +error: + goto exit; +} + +static void +release_ground(ref_T* ref) +{ + struct htrdr_ground* ground; + ASSERT(ref); + ground = CONTAINER_OF(ref, struct htrdr_ground, ref); + if(ground->view) S3D(scene_view_ref_put(ground->view)); + MEM_RM(ground->htrdr->allocator, ground); +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +htrdr_ground_create + (struct htrdr* htrdr, + const char* obj_filename, + const int repeat_ground, /* Infinitely repeat the ground in X and Y */ + struct htrdr_ground** out_ground) +{ + char buf[128]; + struct htrdr_ground* ground = NULL; + struct time t0, t1; + res_T res = RES_OK; + ASSERT(htrdr && obj_filename && out_ground); + + ground = MEM_CALLOC(htrdr->allocator, 1, sizeof(*ground)); + if(!ground) { + res = RES_MEM_ERR; + htrdr_log_err(htrdr, + "%s: could not allocate the ground data structure -- %s.\n", + FUNC_NAME, res_to_cstr(res)); + goto error; + } + ref_init(&ground->ref); + ground->htrdr = htrdr; + ground->repeat = repeat_ground; + + time_current(&t0); + res = setup_ground(ground, obj_filename); + if(res != RES_OK) goto error; + time_sub(&t0, time_current(&t1), &t0); + time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); + htrdr_log(ground->htrdr, "Setup ground in %s\n", buf); + +exit: + *out_ground = ground; + return res; +error: + if(ground) { + htrdr_ground_ref_put(ground); + ground = NULL; + } + goto exit; +} + +void +htrdr_ground_ref_get(struct htrdr_ground* ground) +{ + ASSERT(ground); + ref_get(&ground->ref); +} + +void +htrdr_ground_ref_put(struct htrdr_ground* ground) +{ + ASSERT(ground); + ref_put(&ground->ref, release_ground); +} + +res_T +htrdr_ground_trace_ray + (struct htrdr_ground* ground, + const double org[3], + const double dir[3], /* Must be normalized */ + const double range[2], + const struct s3d_hit* prev_hit, + struct s3d_hit* hit) +{ + res_T res = RES_OK; + ASSERT(ground && org && dir && range && hit); + + if(!ground->repeat) { + struct ray_context ray_ctx = RAY_CONTEXT_NULL; + float ray_org[3]; + float ray_dir[3]; + + f3_set_d3(ray_org, org); + f3_set_d3(ray_dir, dir); + f2_set_d2(ray_ctx.range, range); + ray_ctx.hit_prev = prev_hit ? *prev_hit : S3D_HIT_NULL; + + res = s3d_scene_view_trace_ray + (ground->view, ray_org, ray_dir, ray_ctx.range, &ray_ctx, hit); + if(res != RES_OK) { + htrdr_log_err(ground->htrdr, + "%s: could not trace the ray against the ground geometry -- %s.\n", + FUNC_NAME, res_to_cstr(res)); + goto error; + } + } else { + struct trace_ground_context slab_ctx = TRACE_GROUND_CONTEXT_NULL; + double low[3], upp[3]; + + *hit = S3D_HIT_NULL; + slab_ctx.view = ground->view; + slab_ctx.context.range[0] = (float)range[0]; + slab_ctx.context.range[1] = (float)range[1]; + slab_ctx.context.hit_prev = prev_hit ? *prev_hit : S3D_HIT_NULL; + slab_ctx.hit = hit; + + d3_set_f3(low, ground->lower); + d3_set_f3(upp, ground->upper); + + res = htrdr_slab_trace_ray(ground->htrdr, org, dir, range, low, upp, + trace_ground, 32, &slab_ctx); + if(res != RES_OK) goto error; + } + +exit: + return res; +error: + goto exit; +} + diff --git a/src/htrdr_ground.h b/src/htrdr_ground.h @@ -0,0 +1,51 @@ +/* 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_GROUND_H +#define HTRDR_GROUND_H + +#include <rsys/rsys.h> + +/* Forward declarations */ +struct htrdr; +struct htrdr_ground; +struct s3d_hit; + +extern LOCAL_SYM res_T +htrdr_ground_create + (struct htrdr* htrdr, + const char* obj_filename, + const int repeat_ground, /* Infinitely repeat the ground in X and Y */ + struct htrdr_ground** ground); + +extern LOCAL_SYM void +htrdr_ground_ref_get + (struct htrdr_ground* ground); + +extern LOCAL_SYM void +htrdr_ground_ref_put + (struct htrdr_ground* ground); + +extern LOCAL_SYM res_T +htrdr_ground_trace_ray + (struct htrdr_ground* ground, + const double ray_origin[3], + const double ray_direction[3], /* Must be normalized */ + const double ray_range[2], + const struct s3d_hit* prev_hit,/* Previous hit. Avoid self hit. May be NULL*/ + struct s3d_hit* hit); + +#endif /* HTRDR_GROUND_H */ + diff --git a/src/htrdr_sky.c b/src/htrdr_sky.c @@ -18,6 +18,7 @@ #include "htrdr.h" #include "htrdr_c.h" #include "htrdr_grid.h" +#include "htrdr_slab.h" #include "htrdr_sky.h" #include "htrdr_sun.h" @@ -113,6 +114,17 @@ struct octree_data { const struct htrdr_sky* sky; }; +struct trace_cloud_context { + struct svx_tree* clouds; + struct svx_hit* hit; + svx_hit_challenge_T challenge; + svx_hit_filter_T filter; + void* context; +}; +static const struct trace_cloud_context TRACE_CLOUD_CONTEXT_NULL = { + NULL, NULL, NULL, NULL, NULL +}; + /* Properties of a short wave spectral band */ struct sw_band_prop { /* Average cross section in the band */ @@ -235,6 +247,26 @@ world_to_cloud return d3_set(out_pos_cs, pos_cs); } +static INLINE res_T +trace_cloud + (const double org[3], + const double dir[3], + const double range[2], + void* context, + int* hit) +{ + struct trace_cloud_context* ctx = context; + res_T res = RES_OK; + ASSERT(org && dir && range && context && hit); + + res = svx_tree_trace_ray(ctx->clouds, org, dir, range, ctx->challenge, + ctx->filter, ctx->context, ctx->hit); + if(res != RES_OK) return res; + + *hit = !SVX_HIT_NONE(ctx->hit); + return RES_OK; +} + static FINLINE float voxel_get (const float* data, @@ -2152,86 +2184,20 @@ htrdr_sky_trace_ray /* Clouds are infinitely repeated along the X and Y axis */ } else { - double pos[2]; - double org_cs[3]; /* Origin of the ray transformed in local cloud space */ - double cloud_low_ws[3]; /* Cloud lower bound in world space */ - double cloud_upp_ws[3]; /* Cloud upper bound in world space */ - double cloud_sz[3]; /* Size of the cloud */ - double t_max[3], t_delta[2]; - int xy[2]; /* 2D index of the repeated cloud */ - int incr[2]; /* Index increment */ - ASSERT(cloud_range[0] < cloud_range[1]); - - /* Compute the size of the cloud */ - cloud_sz[0] = sky->htcp_desc.upper[0] - sky->htcp_desc.lower[0]; - cloud_sz[1] = sky->htcp_desc.upper[1] - sky->htcp_desc.lower[1]; - cloud_sz[2] = sky->htcp_desc.upper[2] - sky->htcp_desc.lower[2]; - - /* Define the 2D index of the current cloud. (0,0) is the index of the - * non duplicated cloud */ - pos[0] = org[0] + cloud_range[0]*dir[0]; - pos[1] = org[1] + cloud_range[0]*dir[1]; - xy[0] = (int)floor((pos[0]-sky->htcp_desc.lower[0])/cloud_sz[0]); - xy[1] = (int)floor((pos[1]-sky->htcp_desc.lower[1])/cloud_sz[1]); - - /* Define the 2D index increment wrt dir sign */ - incr[0] = dir[0] < 0 ? -1 : 1; - incr[1] = dir[1] < 0 ? -1 : 1; - - /* Compute the world space AABB of the repeated cloud currently hit */ - cloud_low_ws[0] = sky->htcp_desc.lower[0] + (double)xy[0]*cloud_sz[0]; - cloud_low_ws[1] = sky->htcp_desc.lower[1] + (double)xy[1]*cloud_sz[1]; - cloud_low_ws[2] = sky->htcp_desc.lower[2]; - cloud_upp_ws[0] = cloud_low_ws[0] + cloud_sz[0]; - cloud_upp_ws[1] = cloud_low_ws[1] + cloud_sz[1]; - cloud_upp_ws[2] = sky->htcp_desc.upper[2]; - - /* Compute the max ray intersection with the current cloud AABB */ - t_max[0] = ((dir[0]<0 ? cloud_low_ws[0] : cloud_upp_ws[0]) - org[0])/dir[0]; - t_max[1] = ((dir[1]<0 ? cloud_low_ws[1] : cloud_upp_ws[1]) - org[1])/dir[1]; - t_max[2] = ((dir[2]<0 ? cloud_low_ws[2] : cloud_upp_ws[2]) - org[2])/dir[2]; - ASSERT(t_max[0] >= 0 && t_max[1] >= 0 && t_max[2] >= 0); - - /* Compute the distance along the ray to traverse in order to move of a - * distance equal to the cloud size along the X and Y axis */ - t_delta[0] = (dir[0]<0 ? -cloud_sz[0]:cloud_sz[0]) / dir[0]; - t_delta[1] = (dir[1]<0 ? -cloud_sz[1]:cloud_sz[1]) / dir[1]; - ASSERT(t_delta[0] >= 0 && t_delta[1] >= 0); - - org_cs[2] = org[2]; - - for(;;) { - int iaxis; - - /* Transform the ray origin in the local cloud space */ - org_cs[0] = sky->htcp_desc.lower[0] + (org[0]-(double)xy[0]*cloud_sz[0]); - org_cs[1] = sky->htcp_desc.lower[1] + (org[1]-(double)xy[1]*cloud_sz[1]); - - /* Trace the ray in the repeated cloud */ - res = svx_tree_trace_ray(clouds, org_cs, dir, cloud_range, challenge, - filter, context, hit); - if(res != RES_OK) { - htrdr_log_err(sky->htrdr, - "%s: could not trace the ray in the repeated clouds.\n", - FUNC_NAME); - goto error; - } - if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */ - - /* Define the next axis to traverse */ - iaxis = t_max[0] < t_max[1] - ? (t_max[0] < t_max[2] ? 0 : 2) - : (t_max[1] < t_max[2] ? 1 : 2); - - if(iaxis == 2) break; /* The ray traverse the infinite cloud slice */ - - if(t_max[iaxis] >= cloud_range[1]) break; /* Out of bound */ - - t_max[iaxis] += t_delta[iaxis]; - - /* Define the 2D index of the next traversed cloud */ - xy[iaxis] += incr[iaxis]; - } + struct trace_cloud_context slab_ctx = TRACE_CLOUD_CONTEXT_NULL; + + slab_ctx.clouds = clouds; + slab_ctx.challenge = challenge; + slab_ctx.filter = filter; + slab_ctx.context = context; + slab_ctx.hit = hit; + + res = htrdr_slab_trace_ray(sky->htrdr, org, dir, cloud_range, + sky->htcp_desc.lower, sky->htcp_desc.upper, trace_cloud, 32, + &slab_ctx); + if(res != RES_OK) goto error; + + if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */ } /* Pursue ray traversal into the atmosphere */ diff --git a/src/htrdr_sky.h b/src/htrdr_sky.h @@ -64,7 +64,7 @@ htrdr_sky_create const char* htgop_filename, const char* htmie_filename, const double optical_thickness, /* Threshold used during octree building */ - const int repeat_clouds, /* Infinitly repeat the clouds in X and Y */ + const int repeat_clouds, /* Infinitely repeat the clouds in X and Y */ struct htrdr_sky** sky); extern LOCAL_SYM void diff --git a/src/htrdr_slab.c b/src/htrdr_slab.c @@ -0,0 +1,133 @@ +/* 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_slab.h" + +#include <rsys/cstr.h> +#include <math.h> + +/******************************************************************************* + * Local function + ******************************************************************************/ +res_T +htrdr_slab_trace_ray + (struct htrdr* htrdr, + const double org[3], + const double dir[3], + const double range[2], + const double cell_low[2], + const double cell_upp[2], + htrdr_trace_cell_T trace_cell, + const size_t max_steps, + void* trace_cell_context) +{ + double pos[2]; + double org_cs[3]; /* Origin of the ray transformed in local cell space */ + double cell_low_ws[3]; /* Cell lower bound in world space */ + double cell_upp_ws[3]; /* Cell upper bound in world space */ + double cell_sz[3]; /* Size of a cell */ + double t_max[3], t_delta[2], t_min_z; + size_t istep; + int xy[2]; /* 2D index of the repeated cell */ + int incr[2]; /* Index increment */ + res_T res = RES_OK; + ASSERT(htrdr && org && dir && range && cell_low && cell_upp && trace_cell); + ASSERT(range[0] < range[1]); + + /* Check that the ray intersects the slab */ + t_min_z = (cell_low[2] - org[2]) / dir[2]; + t_max[2] = (cell_upp[2] - org[2]) / dir[2]; + if(t_min_z > t_max[2]) SWAP(double, t_min_z, t_max[2]); + t_min_z = MMAX(t_min_z, range[0]); + t_max[2] = MMIN(t_max[2], range[1]); + if(t_min_z > t_max[2]) return RES_OK; + + /* Compute the size of a cell */ + cell_sz[0] = cell_upp[0] - cell_low[0]; + cell_sz[1] = cell_upp[1] - cell_low[1]; + cell_sz[2] = cell_upp[2] - cell_low[2]; + + /* Define the 2D index of the current cell. (0,0) is the index of the + * non duplicated cell */ + pos[0] = org[0] + t_min_z*dir[0]; + pos[1] = org[1] + t_min_z*dir[1]; + xy[0] = (int)floor((pos[0] - cell_low[0]) / cell_sz[0]); + xy[1] = (int)floor((pos[1] - cell_low[1]) / cell_sz[1]); + + /* Define the 2D index increment wrt dir sign */ + incr[0] = dir[0] < 0 ? -1 : 1; + incr[1] = dir[1] < 0 ? -1 : 1; + + /* Compute the world space AABB of the repeated cell currently hit */ + cell_low_ws[0] = cell_low[0] + (double)xy[0]*cell_sz[0]; + cell_low_ws[1] = cell_low[1] + (double)xy[1]*cell_sz[1]; + cell_low_ws[2] = cell_low[2]; + cell_upp_ws[0] = cell_low_ws[0] + cell_sz[0]; + cell_upp_ws[1] = cell_low_ws[1] + cell_sz[1]; + cell_upp_ws[2] = cell_upp[2]; + + /* Compute the max ray intersection with the current cell */ + t_max[0] = ((dir[0]<0 ? cell_low_ws[0] : cell_upp_ws[0]) - org[0]) / dir[0]; + t_max[1] = ((dir[1]<0 ? cell_low_ws[1] : cell_upp_ws[1]) - org[1]) / dir[1]; + /*t_max[2] = ((dir[2]<0 ? cell_low_ws[2] : cell_upp_ws[2]) - org[2]) / dir[2];*/ + ASSERT(t_max[0] >= 0 && t_max[1] >= 0 && t_max[2] >= 0); + + /* Compute the distance along the ray to traverse in order to move of a + * distance equal to the cloud size along the X and Y axis */ + t_delta[0] = (dir[0]<0 ? -cell_sz[0] : cell_sz[0]) / dir[0]; + t_delta[1] = (dir[1]<0 ? -cell_sz[1] : cell_sz[1]) / dir[1]; + ASSERT(t_delta[0] >= 0 && t_delta[1] >= 0); + + org_cs[2] = org[2]; + FOR_EACH(istep, 0, max_steps) { + int iaxis; + int hit; + + /* Transform the ray origin in the local cell space */ + org_cs[0] = org[0] - (double)xy[0]*cell_sz[0]; + org_cs[1] = org[1] - (double)xy[1]*cell_sz[1]; + + res = trace_cell(org_cs, dir, range, trace_cell_context, &hit); + if(res != RES_OK) { + htrdr_log_err(htrdr, + "%s: could not trace the ray in the repeated cells -- %s.\n", + FUNC_NAME, res_to_cstr(res)); + goto error; + } + if(hit) goto exit; + + /* Define the next axis to traverse */ + iaxis = t_max[0] < t_max[1] + ? (t_max[0] < t_max[2] ? 0 : 2) + : (t_max[1] < t_max[2] ? 1 : 2); + + if(iaxis == 2) break; /* The ray traverse the slab */ + + if(t_max[iaxis] >= range[1]) break; /* Out of bound */ + + t_max[iaxis] += t_delta[iaxis]; + + /* Define the 2D index of the next traversed cloud */ + xy[iaxis] += incr[iaxis]; + } + +exit: + return res; +error: + goto exit; +} + + diff --git a/src/htrdr_slab.h b/src/htrdr_slab.h @@ -0,0 +1,46 @@ +/* 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_SLAB_H +#define HTRDR_SLAB_H + +#include <rsys/rsys.h> + +/* Forward declaration */ +struct htrdr; + +typedef res_T +(*htrdr_trace_cell_T) + (const double org[3], /* Ray origin */ + const double dir[3], /* Ray direction. Must be normalized */ + const double range[2], /* Ray range */ + void* ctx, /* User defined data */ + int* hit); /* Hit something ? */ + +/* Trace a ray into a slab composed of a cell infinitely repeated in X and Y */ +extern LOCAL_SYM res_T +htrdr_slab_trace_ray + (struct htrdr* htrdr, + const double org[3], + const double dir[3], + const double range[2], + const double cell_low[2], + const double cell_upp[2], + htrdr_trace_cell_T trace_cell, + const size_t max_steps, /* Max traversed cell */ + void* trace_cell_context); + +#endif /* HTRDR_SLAB_H */ +