commit 0eadf0f2c4a207253b895b74463057945d78043b parent 15fc38e1a53f458034471fb9e0e68caa82a80237 Author: Vincent Forest <vincent.forest@meso-star.com> Date: Fri, 23 Apr 2021 12:00:55 +0200 Merge branch 'feature_combustion_chamber' into develop Diffstat:
16 files changed, 639 insertions(+), 106 deletions(-)
diff --git a/cmake/atmosphere/CMakeLists.txt b/cmake/atmosphere/CMakeLists.txt @@ -46,8 +46,8 @@ include_directories( ################################################################################ # Generate files ################################################################################ -set(HTRDR_ATMOSPHERE_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD "1") -set(HTRDR_ATMOSPHERE_ARGS_DEFAULT_SKY_MTL_NAME "\"air\"") +set(HTRDR_ATMOSPHERE_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD "1" CACHE INTERNAL "") +set(HTRDR_ATMOSPHERE_ARGS_DEFAULT_SKY_MTL_NAME "\"air\"" CACHE INTERNAL "") configure_file(${HTRDR_SOURCE_DIR}/atmosphere/htrdr_atmosphere_args.h.in ${HTRDR_BUILD_DIR}/atmosphere/htrdr_atmosphere_args.h @ONLY) diff --git a/cmake/combustion/CMakeLists.txt b/cmake/combustion/CMakeLists.txt @@ -51,6 +51,7 @@ set(HTRDR_COMBUSTION_FILES_SRC htrdr_combustion_args.c htrdr_combustion_draw_map.c htrdr_combustion_compute_radiance_sw.c + htrdr_combustion_geometry_ray_filter.c htrdr_combustion_laser.c htrdr_combustion_main.c) @@ -58,6 +59,7 @@ set(HTRDR_COMBUSTION_FILES_INC htrdr_combustion.h htrdr_combustion_c.h htrdr_combustion_args.h + htrdr_combustion_geometry_ray_filter.h htrdr_combustion_laser.h) # Prepend each file in the `HTRDR_FILES_<SRC|INC>' list by `HTRDR_SOURCE_DIR' diff --git a/cmake/core/CMakeLists.txt b/cmake/core/CMakeLists.txt @@ -53,17 +53,17 @@ include_directories( ################################################################################ # Generate files ################################################################################ -set(HTRDR_ARGS_DEFAULT_CAMERA_POS "0,0,0") -set(HTRDR_ARGS_DEFAULT_CAMERA_TGT "0,1,0") -set(HTRDR_ARGS_DEFAULT_CAMERA_UP "0,0,1") -set(HTRDR_ARGS_DEFAULT_CAMERA_FOV "70") -set(HTRDR_ARGS_DEFAULT_RECTANGLE_POS "0,0,0") -set(HTRDR_ARGS_DEFAULT_RECTANGLE_TGT "0,0,1") -set(HTRDR_ARGS_DEFAULT_RECTANGLE_UP "0,1,0") -set(HTRDR_ARGS_DEFAULT_RECTANGLE_SZ "1,1") -set(HTRDR_ARGS_DEFAULT_IMG_WIDTH "320") -set(HTRDR_ARGS_DEFAULT_IMG_HEIGHT "240") -set(HTRDR_ARGS_DEFAULT_IMG_SPP "1") +set(HTRDR_ARGS_DEFAULT_CAMERA_POS "0,0,0" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_CAMERA_TGT "0,1,0" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_CAMERA_UP "0,0,1" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_CAMERA_FOV "70" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_RECTANGLE_POS "0,0,0" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_RECTANGLE_TGT "0,0,1" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_RECTANGLE_UP "0,1,0" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_RECTANGLE_SZ "1,1" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_IMG_WIDTH "320" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_IMG_HEIGHT "240" CACHE INTERNAL "") +set(HTRDR_ARGS_DEFAULT_IMG_SPP "1" CACHE INTERNAL "") configure_file(${HTRDR_SOURCE_DIR}/core/htrdr_args.h.in ${HTRDR_BUILD_DIR}/core/htrdr_args.h @ONLY) diff --git a/cmake/doc/CMakeLists.txt b/cmake/doc/CMakeLists.txt @@ -1,4 +1,6 @@ -# Copyright (C) 2018-2019 CNRS, |Meso|Star>, Université Paul Sabatier +# 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 diff --git a/doc/htrdr-atmosphere.1.txt.in b/doc/htrdr-atmosphere.1.txt.in @@ -176,7 +176,7 @@ OPTIONS *-n* _sky-mtl_:: Name of the material representing the sky in the *htrdr-materials*(5) file. - By default, _sky-mtl_ is @HTRDR_ARGS_DEFAULT_SKY_MTL_NAME@. + By default, _sky-mtl_ is @HTRDR_ATMOSPHERE_ARGS_DEFAULT_SKY_MTL_NAME@. *-O* _cache_:: File used to cache the sky data. If the _cache_ file does not exists, it is @@ -256,7 +256,8 @@ OPTIONS *-T* _threshold_:: Optical thickness used as threshold criteria to partition the properties of - the clouds. By default its value is @HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@. + the clouds. By default its value is + @HTRDR_ATMOSPHERE_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@. This option is ignored if a cache file is used (option *-O*). *-t* _threads-count_:: diff --git a/src/atmosphere/htrdr_atmosphere_ground.c b/src/atmosphere/htrdr_atmosphere_ground.c @@ -57,11 +57,16 @@ trace_slab int* hit) { struct trace_slab_context* ctx = context; + struct htrdr_geometry_trace_ray_args rt_args = + HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL; res_T res = RES_OK; ASSERT(org && dir && range && context && hit); - res = htrdr_geometry_trace_ray - (ctx->geom, org, dir, range, ctx->hit_prev, ctx->hit); + d3_set(rt_args.ray_org, org); + d3_set(rt_args.ray_dir, dir); + d3_set(rt_args.ray_range, range); + rt_args.hit_from = ctx->hit_prev ? *ctx->hit_prev : S3D_HIT_NULL; + res = htrdr_geometry_trace_ray(ctx->geom, &rt_args, ctx->hit); if(res != RES_OK) return res; *hit = !S3D_HIT_NONE(ctx->hit); @@ -172,9 +177,16 @@ htrdr_atmosphere_ground_trace_ray } if(!ground->repeat) { - res = htrdr_geometry_trace_ray - (ground->geom, org, dir, range, prev_hit, hit); + struct htrdr_geometry_trace_ray_args rt_args = + HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL; + + d3_set(rt_args.ray_org, org); + d3_set(rt_args.ray_dir, dir); + d3_set(rt_args.ray_range, range); + if(prev_hit) rt_args.hit_from = *prev_hit; + res = htrdr_geometry_trace_ray(ground->geom, &rt_args, hit); if(res != RES_OK) goto error; + } else { struct trace_slab_context slab_ctx = TRACE_SLAB_CONTEXT_NULL; double low[3], upp[3]; diff --git a/src/combustion/htrdr_combustion_c.h b/src/combustion/htrdr_combustion_c.h @@ -67,7 +67,6 @@ struct htrdr_combustion { struct str output_name; /* Name of the output stream */ int dump_volumetric_acceleration_structure; - ref_T ref; struct htrdr* htrdr; }; @@ -81,13 +80,13 @@ extern LOCAL_SYM res_T combustion_draw_map (struct htrdr_combustion* cmd); -/* Return the shortwave radiance in W/m^2/sr */ -extern LOCAL_SYM double +extern LOCAL_SYM res_T combustion_compute_radiance_sw (struct htrdr_combustion* cmd, const size_t ithread, struct ssp_rng* rng, const double pos_in[3], - const double dir_in[3]); + const double dir_in[3], + double* out_weigh); /* Shortwave radiance in W/m^2/sr */ #endif /* HTRDR_COMBUSTION_C_H */ diff --git a/src/combustion/htrdr_combustion_compute_radiance_sw.c b/src/combustion/htrdr_combustion_compute_radiance_sw.c @@ -15,10 +15,17 @@ * 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 200112L /* nextafterf */ + #include "combustion/htrdr_combustion_c.h" #include "combustion/htrdr_combustion_laser.h" +#include "combustion/htrdr_combustion_geometry_ray_filter.h" #include "core/htrdr.h" +#include "core/htrdr_log.h" +#include "core/htrdr_geometry.h" +#include "core/htrdr_interface.h" +#include "core/htrdr_materials.h" #include <astoria/atrstm.h> @@ -30,6 +37,14 @@ #include <rsys/double4.h> #include <rsys/dynamic_array.h> +#include <math.h> /* nextafterf */ + +enum scattering_type { + SCATTERING_IN_VOLUME, + SCATTERING_AT_SURFACE, + SCATTERING_NONE +}; + /* Define a position along the ray into the semi-transparent medium */ struct position { double distance; /* Distance up to the position */ @@ -44,7 +59,7 @@ struct position { static const struct position POSITION_NULL = POSITION_NULL__; /* Syntactic sugar to check if the position is valid */ -#define POSITION_NONE(Pos) ((Pos)->distance >= DBL_MAX) +#define POSITION_NONE(Pos) ((Pos)->distance >= FLT_MAX) /* Common position but preferentially sampled within a limited range. Its * associated ksi variable defines the correction of the weight due to the @@ -492,15 +507,27 @@ transmissivity /* This function computes the contribution of the in-laser scattering for the * current scattering position of the reverse optical path 'pos' and current * direction 'dir'. One contribution has to be computed for each scattering - * position in order to compute the value of the Monte-Carlo weight for the - * optical path (weight of the statistical realization). */ + * position 'xsc' in order to compute the value of the Monte-Carlo weight for + * the optical path (weight of the statistical realization). + * __ + * /\ dir + * xout / + * +---------------------x--------------- - - - + * lse | ---> wi / + * (laser surface | ---> <----x xsc + * emission) | ---> / + * +-----------------x------------------- - - - + * / xin + * / + * o pos */ static double laser_once_scattered (struct htrdr_combustion* cmd, const size_t ithread, struct ssp_rng* rng, const double pos[3], - const double dir[3]) + const double dir[3], + const double range_in[2]) { /* RDG-FA phase function */ struct atrstm_rdgfa rdgfa_param = ATRSTM_RDGFA_NULL; @@ -510,6 +537,12 @@ laser_once_scattered SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT; struct ssf_phase* phase = NULL; + /* Surface ray tracing */ + struct htrdr_geometry_trace_ray_args rt_args = + HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL; + struct geometry_ray_filter_context rt_ctx = GEOMETRY_RAY_FILTER_CONTEXT_NULL; + struct s3d_hit hit = S3D_HIT_NULL; + /* Scattering position into the laser sheet */ struct position_limited sc_sample = POSITION_LIMITED_NULL; double xsc[3] = {0, 0, 0}; /* World space position */ @@ -524,7 +557,7 @@ laser_once_scattered double laser_flux_density = 0; /* In W/m^2 */ double wlen = 0; /* In nm */ - /* Miscellaneous varbale */ + /* Miscellaneous variables */ double wi[3] = {0, 0, 0}; double wo[3] = {0, 0, 0}; double range[2] = {0, DBL_MAX}; @@ -533,16 +566,21 @@ laser_once_scattered /* Radiance to compute in W/m^2/sr */ double L = 0; - ASSERT(cmd && rng && pos && dir); + ASSERT(cmd && rng && pos && dir && range_in); + ASSERT(range_in[0] <= range_in[1]); /* Fetch the laser wavelength */ wlen = htrdr_combustion_laser_get_wavelength(cmd->laser); /* Find the intersections along dir with the laser sheet */ + range[0] = range_in[0]; + range[1] = range_in[1]; htrdr_combustion_laser_trace_ray(cmd->laser, pos, dir, range, laser_hit_dst); /* No intersection with the laser sheet => no laser contribution */ if(HTRDR_COMBUSTION_LASER_HIT_NONE(laser_hit_dst)) return 0; + ASSERT(laser_hit_dst[0] <= laser_hit_dst[1]); + ASSERT(laser_hit_dst[1] <= range_in[1]); /* Compute the transmissivity from 'pos' to 'xin' */ range[0] = 0; @@ -573,6 +611,29 @@ laser_once_scattered htrdr_combustion_laser_get_direction(cmd->laser, wi); d3_minus(wi, wi); + /* Find the intersection with the combustion chamber */ + if(cmd->geom) { /* Is there a combustion chamber? */ + /* Setup the ray to trace */ + d3_set(rt_args.ray_org, xsc); + d3_set(rt_args.ray_dir, wi); + d2_set(rt_args.ray_range, range); + rt_args.hit_from = S3D_HIT_NULL; + + /* Configure the "X-ray" surface ray trace filtering. This filtering function + * helps to avoid intersections with the side of the surfaces facing inside + * the combustion chamber to allow the rays to exit out. */ + rt_args.filter = geometry_ray_filter_discard_medium_interface; + rt_args.filter_context = &rt_ctx; + rt_ctx.geom = cmd->geom; + rt_ctx.medium_name = "chamber"; + + HTRDR(geometry_trace_ray(cmd->geom, &rt_args, &hit)); + + /* If a surface was intersected then the laser surface emissions is + * occluded along `wi' and thus there is no laser contribution */ + if(!S3D_HIT_NONE(&hit)) return 0; + } + /* Compute the transmissivity from xsc to the laser surface emission */ range[0] = 0; range[1] = htrdr_combustion_laser_compute_surface_plane_distance @@ -612,16 +673,15 @@ laser_once_scattered return L; } -/******************************************************************************* - * Local functions - ******************************************************************************/ -extern LOCAL_SYM double -combustion_compute_radiance_sw +static void +sample_scattering_direction (struct htrdr_combustion* cmd, const size_t ithread, struct ssp_rng* rng, - const double pos_in[3], - const double dir_in[3]) + const struct position* scattering, + const double wlen, + const double incoming_dir[3], + double scattering_dir[3]) { /* RDG-FA phase function */ struct atrstm_rdgfa rdgfa_param = ATRSTM_RDGFA_NULL; @@ -631,6 +691,206 @@ combustion_compute_radiance_sw SSF_PHASE_RDGFA_SETUP_ARGS_DEFAULT; struct ssf_phase* phase = NULL; + /* Miscellaneous variable */ + double wo[3]; + + ASSERT(cmd && rng && scattering && incoming_dir && scattering_dir); + ASSERT(!POSITION_NONE(scattering)); + + /* Retrieve the RDG-FA phase function parameters */ + fetch_rdgfa_args.wavelength = wlen; + fetch_rdgfa_args.prim = scattering->prim; + fetch_rdgfa_args.bcoords[0] = scattering->bcoords[0]; + fetch_rdgfa_args.bcoords[1] = scattering->bcoords[1]; + fetch_rdgfa_args.bcoords[2] = scattering->bcoords[2]; + fetch_rdgfa_args.bcoords[3] = scattering->bcoords[3]; + ATRSTM(fetch_rdgfa(cmd->medium, &fetch_rdgfa_args, &rdgfa_param)); + + /* Setup the RDG-FA phase function corresponding to the scattering event */ + phase = cmd->rdgfa_phase_functions[ithread]; + setup_rdgfa_args.wavelength = rdgfa_param.wavelength; + setup_rdgfa_args.fractal_dimension = rdgfa_param.fractal_dimension; + setup_rdgfa_args.gyration_radius = rdgfa_param.gyration_radius; + SSF(phase_rdgfa_setup(phase, &setup_rdgfa_args)); + + /* Sample a new optical path direction from the phase function */ + d3_minus(wo, incoming_dir); /* Ensure SSF convention */ + ssf_phase_sample(phase, rng, wo, scattering_dir, NULL); +} + +/* Return the bounce reflectivity */ +static double +sample_bounce_direction + (struct htrdr_combustion* cmd, + const size_t ithread, + struct ssp_rng* rng, + const struct s3d_hit* hit, + const double wlen, + const double incoming_dir[3], + double bounce_dir[3]) +{ + struct htrdr_interface interf = HTRDR_INTERFACE_NULL; + const struct htrdr_mtl* mtl = NULL; + struct ssf_bsdf* bsdf = NULL; + double wo[3]; + double N[3]; + double bounce_reflectivity; + int bsdf_type = 0; + + ASSERT(cmd && rng && hit && incoming_dir && bounce_dir); + ASSERT(!S3D_HIT_NONE(hit)); + + /* Recover the bsdf of the intersected interface */ + htrdr_geometry_get_interface(cmd->geom, hit, &interf); + mtl = htrdr_interface_fetch_hit_mtl(&interf, incoming_dir, hit); + HTRDR(mtl_create_bsdf(cmd->htrdr, mtl, ithread, wlen, rng, &bsdf)); + + /* Surface normal */ + d3_set_f3(N, hit->normal); + d3_normalize(N, N); + if(d3_dot(N, incoming_dir) > 0) d3_minus(N, N); /* Ensure SSF convention */ + + /* Sample a new optical path direction from the brdf function */ + d3_minus(wo, incoming_dir); /* Ensure SSF convention */ + bounce_reflectivity = ssf_bsdf_sample + (bsdf, rng, wo, N, bounce_dir, &bsdf_type, NULL); + + if(!(bsdf_type & SSF_REFLECTION)) { /* Handle only reflections */ + bounce_reflectivity = 0; + } + + SSF(bsdf_ref_put(bsdf)); + + return bounce_reflectivity; +} + +static res_T +move_to_scattering_position + (struct htrdr_combustion* cmd, + const double pos[3], + const double dir[3], + const double sc_distance, + const enum scattering_type sc_type, + double out_pos[3]) +{ + res_T res = RES_OK; + ASSERT(cmd && pos && dir && sc_distance >= 0 && out_pos); + ASSERT(sc_type == SCATTERING_IN_VOLUME || sc_type == SCATTERING_AT_SURFACE); + + /* The scattering position is at a surface or in an open air volume */ + if(sc_type == SCATTERING_AT_SURFACE || cmd->geom == NULL) { + out_pos[0] = pos[0] + dir[0] * sc_distance; + out_pos[1] = pos[1] + dir[1] * sc_distance; + out_pos[2] = pos[2] + dir[2] * sc_distance; + + /* The scattering position is in a volume surrounded by the geometry of the + * combustion chamber. Be careful when moving along 'dir'; due to numerical + * uncertainty, the diffusion position could be outside the combustion + * chamber */ + } else { + const int MAX_ATTEMPTS = 10; /* Max #attempts before reporting an error */ + int iattempt = 0; /* Index of the current attempt */ + + double tmp[3]; /* Temporary vector */ + double pos_to_challenge[3]; /* Scattering position to challenge */ + double dst; /* Distance up to the scattering position */ + + /* Search distance of a near position onto the geometry of the combustion + * chamber */ + double search_radius; + + search_radius = htrdr_geometry_get_epsilon(cmd->geom) * 10.0; + dst = sc_distance; + + while(iattempt < MAX_ATTEMPTS) { + struct htrdr_interface interf = HTRDR_INTERFACE_NULL; + struct s3d_hit hit = S3D_HIT_NULL; + double hit_pos[3]; + double N[3]; + double sign; + const struct htrdr_mtl* mtl = NULL; + + /* Move to the scattering position */ + pos_to_challenge[0] = pos[0] + dir[0] * dst; + pos_to_challenge[1] = pos[1] + dir[1] * dst; + pos_to_challenge[2] = pos[2] + dir[2] * dst; + + /* Find the geometry near the scattering position */ + HTRDR(geometry_find_closest_point + (cmd->geom, pos_to_challenge, search_radius, &hit)); + + /* No geometry => the scattering position to challenge is necessarily + * inside the combustion chamber. */ + if(S3D_HIT_NONE(&hit)) break; + + /* Retrieve the property of the geometry near the scattering position */ + htrdr_geometry_get_hit_position(cmd->geom, &hit, hit_pos); + htrdr_geometry_get_interface(cmd->geom, &hit, &interf); + + /* Fetch the material looking toward the scattering position */ + d3_normalize(N, d3_set_f3(N, hit.normal)); + sign = d3_dot(N, d3_sub(tmp, pos_to_challenge, hit_pos)); + mtl = sign < 0 ? &interf.mtl_back : &interf.mtl_front; + + /* The scattering position is inside the combustion chamber. Everything + * is fine. */ + if(!strcmp(mtl->name, "chamber")) break; + + /* The scattering position is outside the combustion chamber! Due to + * numerical uncertainty, while moving along 'dir' we accidentally + * crossed the combustion chamber geometry. To deal with this problem, we + * try to move slightly less than expected. By "slightly less" we mean + * the next representable single precision floating point value, directly + * less than the previous challenge distance. */ + dst = (double)nextafterf((float)dst, 0); + + /* Next trial */ + iattempt += 1; + } + + /* Max attempts is reached. Report an error */ + if(iattempt == MAX_ATTEMPTS) { + htrdr_log_warn(cmd->htrdr, + "The scattering position {%g, %g, %g} " + "is outside the combustion chamber.\n", + pos_to_challenge[0], + pos_to_challenge[1], + pos_to_challenge[2]); + res = RES_BAD_OP; + goto error; + } + + /* A valid scattering position was found */ + out_pos[0] = pos_to_challenge[0]; + out_pos[1] = pos_to_challenge[1]; + out_pos[2] = pos_to_challenge[2]; + } + +exit: + return res; +error: + goto exit; +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +extern LOCAL_SYM res_T +combustion_compute_radiance_sw + (struct htrdr_combustion* cmd, + const size_t ithread, + struct ssp_rng* rng, + const double pos_in[3], + const double dir_in[3], + double* out_weight) +{ + /* Surface ray tracing */ + struct htrdr_geometry_trace_ray_args rt_args = + HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL; + struct geometry_ray_filter_context rt_ctx = GEOMETRY_RAY_FILTER_CONTEXT_NULL; + struct s3d_hit hit_curr = S3D_HIT_NULL; /* Current surface hit */ + struct s3d_hit hit_prev = S3D_HIT_NULL; /* Previous surface hit */ + /* Transmissivity between the probe position (i.e. 'pos_in') and the current * scattering position over the reverse scattering path */ double Tr_abs = 1; @@ -643,7 +903,9 @@ combustion_compute_radiance_sw double dir[3] = {0, 0, 0}; double range[2] = {0, DBL_MAX}; double wlen = 0; - ASSERT(cmd && rng && pos_in && dir_in); + enum scattering_type sc_type = SCATTERING_NONE; + res_T res = RES_OK; + ASSERT(cmd && rng && pos_in && dir_in && out_weight); d3_set(pos, pos_in); d3_set(dir, dir_in); @@ -654,24 +916,61 @@ combustion_compute_radiance_sw Tr_abs = 1; weight = 0; + /* Configure the "X-ray" surface ray trace filtering. This filtering function + * helps to avoid intersections with the side of the surfaces facing outside + * the combustion chamber to allow primary rays to enter. */ + rt_args.filter = geometry_ray_filter_discard_medium_interface; + rt_args.filter_context = &rt_ctx; + rt_ctx.geom = cmd->geom; + rt_ctx.medium_name = "air"; + for(;;) { struct position scattering = POSITION_NULL; double laser_sheet_emissivity = 0; /* In W/m^2/sr */ double Tr_abs_pos_xsc = 0; - double wo[3]; double wi[3]; + double sc_distance; + + if(cmd->geom) { /* Is there a combustion chamber? */ + /* Find the intersection with the combustion chamber geometry */ + d3_set(rt_args.ray_org, pos); + d3_set(rt_args.ray_dir, dir); + d2_set(rt_args.ray_range, range); + rt_args.hit_from = hit_prev; /* Avoid self intersection */ + HTRDR(geometry_trace_ray(cmd->geom, &rt_args, &hit_curr)); + + /* Deactivate the "X-ray" filtering function. It is only used during the + * first step of the random walk to allow primary rays (i.e. camera rays) + * to enter the combustion chamber. */ + rt_args.filter = NULL; + } /* Handle the laser contribution */ - laser_sheet_emissivity = laser_once_scattered(cmd, ithread, rng, pos, dir); + range[0] = 0; + range[1] = hit_curr.distance; + + laser_sheet_emissivity = laser_once_scattered(cmd, ithread, rng, pos, dir, range); weight += Tr_abs * laser_sheet_emissivity; /* Sample a scattering position */ + range[0] = 0; + range[1] = hit_curr.distance; sample_position(cmd, rng, ATRSTM_RADCOEF_Ks, pos, dir, range, &scattering); - if(POSITION_NONE(&scattering)) break; + + if(!POSITION_NONE(&scattering)) { /* Scattering event in volume */ + sc_type = SCATTERING_IN_VOLUME; + sc_distance = scattering.distance; + } else if(!S3D_HIT_NONE(&hit_curr)) { /* Scattering event at surface */ + sc_type = SCATTERING_AT_SURFACE; + sc_distance = hit_curr.distance; + } else { /* No scattering event. Stop the random walk */ + sc_type = SCATTERING_NONE; + break; + } /* Compute the absorption transmissivity */ range[0] = 0; - range[1] = scattering.distance; + range[1] = sc_distance; Tr_abs_pos_xsc = transmissivity(cmd, rng, ATRSTM_RADCOEF_Ka, pos, dir, range); if(Tr_abs_pos_xsc == 0) break; @@ -679,34 +978,41 @@ combustion_compute_radiance_sw Tr_abs *= Tr_abs_pos_xsc; /* Update the position of the optical path */ - pos[0] = pos[0] + dir[0] * scattering.distance; - pos[1] = pos[1] + dir[1] * scattering.distance; - pos[2] = pos[2] + dir[2] * scattering.distance; - - /* Retrieve the RDG-FA phase function parameters */ - fetch_rdgfa_args.wavelength = wlen; - fetch_rdgfa_args.prim = scattering.prim; - fetch_rdgfa_args.bcoords[0] = scattering.bcoords[0]; - fetch_rdgfa_args.bcoords[1] = scattering.bcoords[1]; - fetch_rdgfa_args.bcoords[2] = scattering.bcoords[2]; - fetch_rdgfa_args.bcoords[3] = scattering.bcoords[3]; - ATRSTM(fetch_rdgfa(cmd->medium, &fetch_rdgfa_args, &rdgfa_param)); - - /* Setup the RDG-FA phase function corresponding to the scattering event */ - phase = cmd->rdgfa_phase_functions[ithread]; - setup_rdgfa_args.wavelength = rdgfa_param.wavelength; - setup_rdgfa_args.fractal_dimension = rdgfa_param.fractal_dimension; - setup_rdgfa_args.gyration_radius = rdgfa_param.gyration_radius; - SSF(phase_rdgfa_setup(phase, &setup_rdgfa_args)); - - /* Sample a new optical path direction from the phase function */ - d3_minus(wo, dir); /* Ensure SSF convention */ - ssf_phase_sample(phase, rng, wo, wi, NULL); + res = move_to_scattering_position(cmd, pos, dir, sc_distance, sc_type, pos); + if(res != RES_OK) goto error; + + if(sc_type == SCATTERING_IN_VOLUME) { + /* Sample a new optical path direction from the medium phase function */ + sample_scattering_direction(cmd, ithread, rng, &scattering, wlen, dir, wi); + + } else { + double bounce_reflectivity; + double r; + ASSERT(sc_type == SCATTERING_AT_SURFACE); + + /* Sample a new optical path direction from the surface BSDF */ + bounce_reflectivity = sample_bounce_direction + (cmd, ithread, rng, &hit_curr, wlen, dir, wi); + + /* Russian roulette wrt surface scattering */ + r = ssp_rng_canonical(rng); + if(r > bounce_reflectivity) break; + + /* Register the current hit to avoid a self intersection for the next + * optical path direction */ + hit_prev = hit_curr; + } /* Update the optical path direction */ dir[0] = wi[0]; dir[1] = wi[1]; dir[2] = wi[2]; } - return weight; + +exit: + *out_weight = weight; + return res; +error: + weight = NaN; + goto exit; } diff --git a/src/combustion/htrdr_combustion_draw_map.c b/src/combustion/htrdr_combustion_draw_map.c @@ -40,6 +40,7 @@ draw_pixel struct htrdr_combustion* cmd = NULL; struct combustion_pixel* pixel = NULL; size_t isamp; + res_T res = RES_OK; ASSERT(htrdr && htrdr_draw_pixel_args_check(args) && data); (void)htrdr; @@ -68,8 +69,9 @@ draw_pixel htrdr_camera_ray(cmd->camera, pix_samp, ray_org, ray_dir); /* Backward trace the path */ - weight = combustion_compute_radiance_sw(cmd, args->ithread, args->rng, - ray_org, ray_dir); + res = combustion_compute_radiance_sw(cmd, args->ithread, args->rng, + ray_org, ray_dir, &weight); + if(res != RES_OK) continue; /* Reject the path */ /* End the registration of the per realisation time */ time_sub(&t0, time_current(&t1), &t0); diff --git a/src/combustion/htrdr_combustion_geometry_ray_filter.c b/src/combustion/htrdr_combustion_geometry_ray_filter.c @@ -0,0 +1,68 @@ +/* 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_geometry_ray_filter.h" + +#include "core/htrdr_interface.h" +#include "core/htrdr_geometry.h" + +#include <rsys/float3.h> + +#include <string.h> + +/******************************************************************************* + * Local functions + ******************************************************************************/ +int +geometry_ray_filter_discard_medium_interface + (const struct s3d_hit* hit, + const float ray_org[3], + const float ray_dir[3], + void* ray_data, + void* filter_data) +{ + struct geometry_ray_filter_context* ctx = ray_data; + struct htrdr_interface interf = HTRDR_INTERFACE_NULL; + const struct htrdr_mtl* mtl = NULL; + float N[3]; + int hit_front; + int discard = 0; + ASSERT(hit && ray_org && ray_dir && ray_data && !S3D_HIT_NONE(hit)); + (void)ray_org, (void)ray_dir, (void)filter_data; + + /* Recover the interface of the intersected surface */ + htrdr_geometry_get_interface(ctx->geom, hit, &interf); + + /* Define if the ray intersects the front face of the surface */ + f3_normalize(N, hit->normal); /* Limit the numerical instabilities */ + hit_front = f3_dot(ray_dir, N) < 0; + + /* Recover the material in which the ray is traced */ + if(hit_front) { + mtl = &interf.mtl_front; + } else { + mtl = &interf.mtl_back; + } + + /* The material should be semi-transparent and thus could not have a BRDF */ + ASSERT(mtl->mrumtl == NULL); + + /* Discard the intersection if the ray comes from the material to filter */ + discard = !strcmp(mtl->name, ctx->medium_name); + return discard; +} + diff --git a/src/combustion/htrdr_combustion_geometry_ray_filter.h b/src/combustion/htrdr_combustion_geometry_ray_filter.h @@ -0,0 +1,47 @@ +/* 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_GEOMETRY_RAY_FILTER_H +#define HTRDR_COMBUSTION_GEOMETRY_RAY_FILTER_H + +#include <star/s3d.h> +#include <rsys/rsys.h> + +struct geometry_ray_filter_context { + struct htrdr_geometry* geom; + + /* Name of the medium attached to the interface to filter. An intersection is + * ignorer if the side of the intersected surface points toward the medium + * whose name is `medium_name`. */ + const char* medium_name; +}; + +#define GEOMETRY_RAY_FILTER_CONTEXT_NULL__ {NULL, NULL} +static const struct geometry_ray_filter_context +GEOMETRY_RAY_FILTER_CONTEXT_NULL = GEOMETRY_RAY_FILTER_CONTEXT_NULL__; + +/* Filter function used to ignore the intersections with surfaces pointing + * toward a user defined medium */ +extern LOCAL_SYM int +geometry_ray_filter_discard_medium_interface + (const struct s3d_hit* hit, + const float ray_org[3], + const float ray_dir[3], + void* ray_data, /* Must point to a struct geometry_ray_filter_context */ + void* filter_data); + +#endif /* HTRDR_COMBUSTION_GEOMETRY_RAY_FILTER_H */ diff --git a/src/combustion/htrdr_combustion_laser.h b/src/combustion/htrdr_combustion_laser.h @@ -46,7 +46,7 @@ struct htrdr_combustion_laser_mesh { }; /* Syntactic sugar to check if a laser sheet hit is valid */ -#define HTRDR_COMBUSTION_LASER_HIT_NONE(Hit) ((Hit)[0] >= DBL_MAX) +#define HTRDR_COMBUSTION_LASER_HIT_NONE(Hit) ((Hit)[0] >= FLT_MAX) /* Forward declaration */ struct htrdr; diff --git a/src/core/htrdr_args.c b/src/core/htrdr_args.c @@ -357,8 +357,6 @@ parse_geometry_parameter(const char* str, void* ptr) goto error; } - htrdr_args_geometry_free(geom); - #define SET_VALUE(Key, Val, Str) { \ const size_t len = strlen(Val) + 1; \ Str = mem_alloc(len); \ @@ -440,6 +438,8 @@ htrdr_args_geometry_parse(struct htrdr_args_geometry* geom, const char* str) goto error; } + htrdr_args_geometry_free(geom); + res = cstr_parse_list(str, ':', parse_geometry_parameter, geom); if(res != RES_OK) goto error; diff --git a/src/core/htrdr_geometry.c b/src/core/htrdr_geometry.c @@ -61,9 +61,12 @@ static const struct mesh MESH_NULL; struct ray_context { float range[2]; - struct s3d_hit hit_prev; + struct s3d_hit hit_from; + + s3d_hit_filter_function_T user_filter; /* May be NULL */ + void* user_filter_data; /* May be NULL */ }; -#define RAY_CONTEXT_NULL__ {{0,INF}, S3D_HIT_NULL__} +#define RAY_CONTEXT_NULL__ {{0,INF}, S3D_HIT_NULL__, NULL, NULL} static const struct ray_context RAY_CONTEXT_NULL = RAY_CONTEXT_NULL__; struct htrdr_geometry { @@ -72,6 +75,10 @@ struct htrdr_geometry { float upper[3]; /* Ground upper bound */ int repeat; /* Make the geom infinite in X and Y */ + /* A empirical value relative to the extent of the geometry that represents + * the threshold below which a numerical problem could occur. */ + float epsilon; + struct htable_interface interfaces; /* Map a Star3D shape to its interface */ struct htrdr* htrdr; @@ -99,6 +106,44 @@ hit_on_edge(const struct s3d_hit* hit) || eq_epsf(w, 1.f, on_edge_eps); } +/* Return 1 if the submitted hit is actually a self hit, i.e. if the ray starts + * the primitive from which it starts */ +static INLINE int +self_hit + (const struct s3d_hit* hit, + const float ray_org[3], + const float ray_dir[3], + const float ray_range[2], + const struct s3d_hit* hit_from) +{ + ASSERT(hit && hit_from); + (void)ray_org, (void)ray_dir; + + /* Internally, Star-3D relies on Embree that, due to numerically inaccuracy, + * can find hits whose distances are not strictly included in the submitted + * ray range. Discard these hits. */ + if(hit->distance <= ray_range[0] || hit->distance >= ray_range[1]) + return 1; + + /* The hit primitive is the one from which the the ray starts. Discard these + * hits */ + if(S3D_PRIMITIVE_EQ(&hit->prim, &hit_from->prim)) + return 1; + + /* If the targeted point is near of the origin, check that it lies on an + * edge/vertex shared by the 2 primitives. In such situation, we assume that + * it is a self intersection and discard this hit */ + if(!S3D_HIT_NONE(hit_from) + && eq_epsf(hit->distance, 0, 1.e-1f) + && hit_on_edge(hit_from) + && hit_on_edge(hit)) { + return 1; + } + + /* No self hit */ + return 0; +} + static int geometry_filter (const struct s3d_hit* hit, @@ -108,20 +153,19 @@ geometry_filter void* filter_data) { const struct ray_context* ray_ctx = ray_data; - (void)ray_org, (void)ray_dir, (void)filter_data; + (void)filter_data; - if(!ray_ctx) return 0; + if(!ray_ctx) /* Nothing to do */ + return 0; - if(S3D_PRIMITIVE_EQ(&hit->prim, &ray_ctx->hit_prev.prim)) return 1; + if(self_hit(hit, ray_org, ray_dir, ray_ctx->range, &ray_ctx->hit_from)) + return 1; /* Discard this hit */ - 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); - } + if(!ray_ctx->user_filter) /* That's all */ + return 0; - return hit->distance <= ray_ctx->range[0] - || hit->distance >= ray_ctx->range[1]; + return ray_ctx->user_filter /* Invoke user filtering */ + (hit, ray_org, ray_dir, ray_ctx->user_filter_data, filter_data); } static res_T @@ -541,6 +585,10 @@ htrdr_geometry_create { char buf[128]; struct htrdr_geometry* geom = NULL; + double low[3]; + double upp[3]; + double tmp[3]; + double extent; struct time t0, t1; res_T res = RES_OK; ASSERT(htrdr && obj_filename && mats && out_ground); @@ -568,6 +616,10 @@ htrdr_geometry_create time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); htrdr_log(geom->htrdr, "Setup geom in %s\n", buf); + htrdr_geometry_get_aabb(geom, low, upp); + extent = d3_len(d3_sub(tmp, upp, low)); + geom->epsilon = MMAX((float)(extent * 1e-6), FLT_EPSILON); + exit: *out_ground = geom; return res; @@ -608,25 +660,40 @@ htrdr_geometry_get_interface *out_interface = *interf; } +void +htrdr_geometry_get_hit_position + (const struct htrdr_geometry* geom, + const struct s3d_hit* hit, + double position[3]) +{ + struct s3d_attrib attr; + ASSERT(geom && hit && position && !S3D_HIT_NONE(hit)); + (void)geom; + + S3D(primitive_get_attrib(&hit->prim, S3D_POSITION, hit->uv, &attr)); + position[0] = attr.value[0]; + position[1] = attr.value[1]; + position[2] = attr.value[2]; +} + res_T htrdr_geometry_trace_ray (struct htrdr_geometry* geom, - const double org[3], - const double dir[3], /* Must be normalized */ - const double range[2], - const struct s3d_hit* prev_hit, + const struct htrdr_geometry_trace_ray_args* args, struct s3d_hit* hit) { struct ray_context ray_ctx = RAY_CONTEXT_NULL; float ray_org[3]; float ray_dir[3]; res_T res = RES_OK; - ASSERT(geom && org && dir && range && hit); + ASSERT(geom && args && hit); - 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; + f3_set_d3(ray_org, args->ray_org); + f3_set_d3(ray_dir, args->ray_dir); + f2_set_d2(ray_ctx.range, args->ray_range); + ray_ctx.hit_from = args->hit_from; + ray_ctx.user_filter = args->filter; + ray_ctx.user_filter_data = args->filter_context; res = s3d_scene_view_trace_ray (geom->view, ray_org, ray_dir, ray_ctx.range, &ray_ctx, hit); @@ -684,3 +751,10 @@ htrdr_geometry_get_aabb d3_set_f3(lower, geom->lower); d3_set_f3(upper, geom->upper); } + +double +htrdr_geometry_get_epsilon(const struct htrdr_geometry* geom) +{ + ASSERT(geom); + return geom->epsilon; +} diff --git a/src/core/htrdr_geometry.h b/src/core/htrdr_geometry.h @@ -20,6 +20,8 @@ #include "core/htrdr.h" +#include <star/s3d.h> + /* Forware declarations */ struct htrdr; struct htrdr_geometry; @@ -28,6 +30,26 @@ struct htrdr_materials; struct s3d_hit; struct ssf_bsdf; +struct htrdr_geometry_trace_ray_args { + double ray_org[3]; + double ray_dir[3]; + double ray_range[2]; + struct s3d_hit hit_from; /* Hit from which the ray starts */ + s3d_hit_filter_function_T filter; /* NULL <=> no user defined filter */ + void* filter_context; +}; + +#define HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL__ { \ + {0,0,0}, /* Ray origin */ \ + {0,0,1}, /* Ray direction */ \ + {0,DBL_MAX}, /* Ray range */ \ + S3D_HIT_NULL, /* Hit from */ \ + NULL, /* User defined filter function */ \ + NULL /* User filter function */ \ +} +static const struct htrdr_geometry_trace_ray_args +HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL = HTRDR_GEOMETRY_TRACE_RAY_ARGS_NULL__; + BEGIN_DECLS HTRDR_CORE_API res_T @@ -51,24 +73,16 @@ htrdr_geometry_get_interface const struct s3d_hit* hit, struct htrdr_interface* interface); -HTRDR_CORE_API res_T -htrdr_geometry_create_bsdf - (struct htrdr_geometry* geom, - const size_t ithread, - const double wavelength, - const double pos[3], - const double dir[3], /* Incoming ray */ +HTRDR_CORE_API void +htrdr_geometry_get_hit_position + (const struct htrdr_geometry* geom, const struct s3d_hit* hit, - struct htrdr_interface* interf, /* NULL <=> do not return the interface */ - struct ssf_bsdf** bsdf); + double position[3]); HTRDR_CORE_API res_T htrdr_geometry_trace_ray (struct htrdr_geometry* geom, - 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*/ + const struct htrdr_geometry_trace_ray_args* args, struct s3d_hit* hit); HTRDR_CORE_API res_T @@ -84,6 +98,12 @@ htrdr_geometry_get_aabb double lower[3], double upper[3]); +/* Empirical value relative to the extent of the geometry that represents the + * threshold below which a numerical problem could occur. */ +HTRDR_CORE_API double +htrdr_geometry_get_epsilon + (const struct htrdr_geometry* geom); + END_DECLS #endif /* HTRDR_GEOMETRY_H */ diff --git a/src/core/htrdr_interface.h b/src/core/htrdr_interface.h @@ -39,7 +39,7 @@ static INLINE const struct htrdr_mtl* htrdr_interface_fetch_hit_mtl (const struct htrdr_interface* interf, const double dir[3], /* Incoming ray */ - struct s3d_hit* hit) + const struct s3d_hit* hit) { const struct htrdr_mtl* mtl = NULL; enum { FRONT, BACK };