htrdr

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

commit f753b2366a3fd74f233a693e309ea64aecc6f269
parent a1aef46023def6863cbdcda68cec023c96c1a463
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 16 Oct 2018 16:22:11 +0200

Add the htrdr_slab API

Diffstat:
Mcmake/CMakeLists.txt | 2++
Msrc/htrdr_c.h | 4++++
Msrc/htrdr_ground.c | 13++++++++++++-
Msrc/htrdr_sky.c | 122+++++++++++++++++++++++++++----------------------------------------------------
Asrc/htrdr_slab.c | 123+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/htrdr_slab.h | 45+++++++++++++++++++++++++++++++++++++++++++++
6 files changed, 228 insertions(+), 81 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -69,6 +69,7 @@ set(HTRDR_FILES_SRC htrdr_ground.c htrdr_main.c htrdr_rectangle.c + htrdr_slab.c htrdr_sky.c htrdr_sun.c) set(HTRDR_FILES_INC @@ -80,6 +81,7 @@ set(HTRDR_FILES_INC 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.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_ground.c b/src/htrdr_ground.c @@ -33,8 +33,11 @@ static const struct ray_context RAY_CONTEXT_NULL = { struct htrdr_ground { struct s3d_scene_view* view; - struct htrdr* htrdr; + 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; }; @@ -176,6 +179,14 @@ htrdr_ground_create goto error; } + res = s3d_scene_view_get_aabb(ground->view, ground->lower, ground->upper); + if(res != RES_OK) { + htrdr_log_err(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)); 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,23 @@ 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 && ctx->hit); + res = svx_tree_trace_ray(ctx->clouds, org, dir, range, ctx->challenge, + ctx->filter, ctx->context, ctx->hit); + *hit = !SVX_HIT_NONE(ctx->hit); + return res; +} + static FINLINE float voxel_get (const float* data, @@ -2152,86 +2181,19 @@ 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, &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_slab.c b/src/htrdr_slab.c @@ -0,0 +1,123 @@ +/* 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, + 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]; + 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]); + + /* 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] + range[0]*dir[0]; + pos[1] = org[1] + range[0]*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(;;) { + int iaxis; + int hit; + + /* Transform the ray origin in the local cell space */ + org_cs[0] = cell_low[0] + (org[0] - (double)xy[0]*cell_sz[0]); + org_cs[1] = cell_low[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,45 @@ +/* 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, + void* trace_cell_context); + +#endif /* HTRDR_SLAB_H */ +