commit 61583fa27a1d20c938460c5fa3f01039a27cde06
parent f8f35d5bf11c7697e87401c5d7a4723bca952c0b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 1 Oct 2018 09:38:53 +0200
Merge branch 'feature_atmosphere'
Diffstat:
9 files changed, 864 insertions(+), 260 deletions(-)
diff --git a/doc/cli.txt b/doc/cli.txt
@@ -9,7 +9,7 @@
<camopt> ::= <target>
| <position>
| <up>
- | <fov>
+ | <fov> # Vertical field of view, in degrees
<target> ::= tgt=<float3>
<position> ::= pos=<float3>
<up> ::= up=<float3>
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -189,6 +189,76 @@ error:
goto exit;
}
+static res_T
+open_file_stamp
+ (struct htrdr* htrdr,
+ const char* filename,
+ struct stat* out_stat, /* Stat of the submitted filename */
+ int* out_fd, /* Descriptor of the opened file. Must be closed by the caller */
+ struct str* stamp_filename)
+{
+ struct stat statbuf;
+ struct str str;
+ int err;
+ int fd = -1;
+ res_T res = RES_OK;
+ ASSERT(htrdr && filename && out_fd && out_stat && stamp_filename);
+
+ str_init(htrdr->allocator, &str);
+
+ err = stat(filename, &statbuf);
+ if(err) {
+ htrdr_log_err(htrdr, "%s: could not stat the file -- %s.\n",
+ filename, strerror(errno));
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ if(!S_ISREG(statbuf.st_mode)) {
+ htrdr_log_err(htrdr, "%s: not a regular file.\n", filename);
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ res = create_directory(htrdr, ".htrdr/");
+ if(res != RES_OK) goto error;
+
+ #define CHK_STR(Func, ErrMsg) { \
+ res = str_##Func; \
+ if(res != RES_OK) { \
+ htrdr_log_err(htrdr, "%s: "ErrMsg"\n", filename); \
+ goto error; \
+ } \
+ } (void)0
+ CHK_STR(set(&str, filename), "could not copy the filename");
+ CHK_STR(set(&str, basename(str_get(&str))), "could not setup the basename");
+ CHK_STR(insert(&str, 0, ".htrdr/"), "could not setup the stamp directory");
+ CHK_STR(append(&str, ".stamp"), "could not setup the stamp extension");
+ #undef CHK_STR
+
+ fd = open(str_cget(&str), O_CREAT|O_RDWR, S_IRUSR|S_IWUSR);
+ if(fd < 0) {
+ htrdr_log_err(htrdr, "%s: could not open/create the file -- %s.\n",
+ str_cget(&str), strerror(errno));
+ res = RES_IO_ERR;
+ goto error;
+ }
+
+ CHK(str_copy_and_clear(stamp_filename, &str) == RES_OK);
+
+exit:
+ str_release(&str);
+ *out_fd = fd;
+ *out_stat = statbuf;
+ return res;
+error:
+ if(fd >= 0) {
+ CHK(close(fd) == 0);
+ fd = -1;
+ }
+ goto exit;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -253,11 +323,14 @@ htrdr_init
goto error;
}
+ res = setup_geometry(htrdr, args->filename_obj);
+ 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, MDEG2RAD(args->camera.fov_x), &htrdr->cam);
+ args->camera.up, proj_ratio, MDEG2RAD(args->camera.fov_y), &htrdr->cam);
if(res != RES_OK) goto error;
res = htrdr_buffer_create(htrdr,
@@ -274,7 +347,8 @@ htrdr_init
htrdr_sun_set_direction(htrdr->sun, args->main_dir);
res = htrdr_sky_create(htrdr, htrdr->sun, args->filename_les,
- args->filename_gas, args->filename_mie, &htrdr->sky);
+ args->filename_gas, args->filename_mie, args->optical_thickness,
+ &htrdr->sky);
if(res != RES_OK) goto error;
htrdr->lifo_allocators = MEM_CALLOC
@@ -297,10 +371,7 @@ htrdr_init
}
}
- res = setup_geometry(htrdr, args->filename_obj);
- if(res != RES_OK) goto error;
-
-exit:
+ exit:
return res;
error:
htrdr_release(htrdr);
@@ -338,9 +409,14 @@ htrdr_run(struct htrdr* htrdr)
size_t i;
FOR_EACH(i, 0, nbands) {
const size_t iband = htrdr_sky_get_sw_spectral_band_id(htrdr->sky, i);
- res = htrdr_sky_dump_clouds_vtk(htrdr->sky, iband, 0, htrdr->output);
- if(res != RES_OK) goto error;
- fprintf(htrdr->output, "---\n");
+ const size_t nquads = htrdr_sky_get_sw_spectral_band_quadrature_length
+ (htrdr->sky, iband);
+ size_t iquad;
+ FOR_EACH(iquad, 0, nquads) {
+ res = htrdr_sky_dump_clouds_vtk(htrdr->sky, iband, iquad, htrdr->output);
+ if(res != RES_OK) goto error;
+ fprintf(htrdr->output, "---\n");
+ }
}
} else {
struct time t0, t1;
@@ -396,7 +472,7 @@ htrdr_log_warn(struct htrdr* htrdr, const char* msg, ...)
/*******************************************************************************
* Local functions
******************************************************************************/
-extern LOCAL_SYM res_T
+extern LOCAL_SYM res_T
open_output_stream
(struct htrdr* htrdr,
const char* filename,
@@ -454,90 +530,94 @@ error:
res_T
is_file_updated(struct htrdr* htrdr, const char* filename, int* out_upd)
{
- struct str str;
+ struct str stamp_filename;
struct stat statbuf;
ssize_t n;
off_t size;
struct timespec mtime;
int fd = -1;
- int err;
int upd = 1;
res_T res = RES_OK;
ASSERT(htrdr && filename && out_upd);
- str_init(htrdr->allocator, &str);
+ str_init(htrdr->allocator, &stamp_filename);
- err = stat(filename, &statbuf);
- if(err) {
- htrdr_log_err(htrdr, "%s: could not stat the file -- %s.\n",
- filename, strerror(errno));
+ res = open_file_stamp(htrdr, filename, &statbuf, &fd, &stamp_filename);
+ if(res != RES_OK) goto error;
+
+ n = read(fd, &mtime, sizeof(mtime));
+ if(n < 0) {
+ htrdr_log_err(htrdr, "%s: could not read the `mtime' data -- %s.\n",
+ str_cget(&stamp_filename), strerror(errno));
res = RES_IO_ERR;
goto error;
}
- if(!S_ISREG(statbuf.st_mode)) {
- htrdr_log_err(htrdr, "%s: not a regular file.\n", filename);
- res = RES_IO_ERR;
- goto error;
+ upd = (size_t)n != sizeof(mtime)
+ ||mtime.tv_nsec != statbuf.st_mtim.tv_nsec
+ ||mtime.tv_sec != statbuf.st_mtim.tv_sec;
+
+ if(!upd) {
+ n = read(fd, &size, sizeof(size));
+ if(n < 0) {
+ htrdr_log_err(htrdr, "%s: could not read the `size' data -- %s.\n",
+ str_cget(&stamp_filename), strerror(errno));
+ res = RES_IO_ERR;
+ goto error;
+ }
+ upd = (size_t)n != sizeof(size) || statbuf.st_size != size;
}
- res = create_directory(htrdr, ".htrdr/");
- if(res != RES_OK) goto error;
+exit:
+ *out_upd = upd;
+ str_release(&stamp_filename);
+ if(fd >= 0) CHK(close(fd) == 0);
+ return res;
+error:
+ goto exit;
+}
- #define CHK_STR(Func, ErrMsg) { \
- res = str_##Func; \
- if(res != RES_OK) { \
- htrdr_log_err(htrdr, "%s: "ErrMsg"\n", filename); \
- goto error; \
- } \
- } (void)0
- CHK_STR(set(&str, filename), "could not copy the filename");
- CHK_STR(set(&str, basename(str_get(&str))), "could not setup the basename");
- CHK_STR(insert(&str, 0, ".htrdr/"), "could not setup the stamp directory");
- CHK_STR(append(&str, ".stamp"), "could not setup the stamp extension");
- #undef CHK_STR
+
+res_T
+update_file_stamp(struct htrdr* htrdr, const char* filename)
+{
+ struct str stamp_filename;
+ struct stat statbuf;
+ int fd = -1;
+ ssize_t n;
+ res_T res = RES_OK;
+ ASSERT(htrdr && filename);
+
+ str_init(htrdr->allocator, &stamp_filename);
+
+ res = open_file_stamp(htrdr, filename, &statbuf, &fd, &stamp_filename);
+ if(res != RES_OK) goto error;
#define CHK_IO(Func, ErrMsg) { \
if((Func) < 0) { \
htrdr_log_err(htrdr, "%s: "ErrMsg" -- %s.\n", \
- str_cget(&str), strerror(errno)); \
+ str_cget(&stamp_filename), strerror(errno)); \
res = RES_IO_ERR; \
goto error; \
} \
} (void) 0
- fd = open(str_cget(&str), O_CREAT|O_RDWR, S_IRUSR|S_IWUSR);
- CHK_IO(fd, "could not open/create the file");
-
- CHK_IO(n = read(fd, &mtime, sizeof(mtime)), "could not read the `mtime' data");
-
- upd = (size_t)n != sizeof(mtime)
- ||mtime.tv_nsec != statbuf.st_mtim.tv_nsec
- ||mtime.tv_sec != statbuf.st_mtim.tv_sec;
-
- if(!upd) {
- CHK_IO(n = read(fd, &size, sizeof(size)), "could not read the `size' data");
- upd = (size_t)n != sizeof(size) || statbuf.st_size != size;
- }
-
- if(upd) {
- CHK_IO(lseek(fd, 0, SEEK_SET), "could not rewind the file descriptor");
+ CHK_IO(lseek(fd, 0, SEEK_SET), "could not rewind the file descriptor");
- /* NOTE: Ignore n >=0 but != sizeof(DATA). In such case stamp is currupted
- * and on the next invocation on the same filename, this function will
- * return 1 */
- n = write(fd, &statbuf.st_mtim, sizeof(statbuf.st_mtim));
- CHK_IO(n, "could not update the `mtime' data");
- n = write(fd, &statbuf.st_size, sizeof(statbuf.st_size));
- CHK_IO(n, "could not update the `size' data");
+ /* NOTE: Ignore n >=0 but != sizeof(DATA). In such case stamp is currupted
+ * and on the next invocation on the same filename, this function will
+ * return 1 */
+ n = write(fd, &statbuf.st_mtim, sizeof(statbuf.st_mtim));
+ CHK_IO(n, "could not update the `mtime' data");
+ n = write(fd, &statbuf.st_size, sizeof(statbuf.st_size));
+ CHK_IO(n, "could not update the `size' data");
- CHK_IO(fsync(fd), "could not sync the file with storage device");
- }
+ CHK_IO(fsync(fd), "could not sync the file with storage device");
#undef CHK_IO
+
exit:
- *out_upd = upd;
- str_release(&str);
+ str_release(&stamp_filename);
if(fd >= 0) CHK(close(fd) == 0);
return res;
error:
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -58,6 +58,10 @@ print_help(const char* cmd)
" -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(
" -v make the program more verbose.\n");
printf("\n");
printf(
@@ -211,7 +215,7 @@ parse_camera_parameter(struct htrdr_args* args, const char* str)
} 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));
+ PARSE("field-of-view", parse_fov(val, &args->camera.fov_y));
} else {
fprintf(stderr, "Invalid camera parameter `%s'.\n", key);
res = RES_BAD_ARG;
@@ -294,7 +298,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:t:v")) != -1) {
+ while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:m:o:T:t:v")) != -1) {
switch(opt) {
case 'a': args->filename_gas = optarg; break;
case 'C':
@@ -317,6 +321,10 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
break;
case 'm': args->filename_mie = optarg; break;
case 'o': args->output = optarg; break;
+ case 'T':
+ res = cstr_to_double(optarg, &args->optical_thickness);
+ if(res == RES_OK && args->optical_thickness < 0) res = RES_BAD_ARG;
+ 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;
diff --git a/src/htrdr_args.h b/src/htrdr_args.h
@@ -36,7 +36,7 @@ struct htrdr_args {
double pos[3];
double tgt[3];
double up[3];
- double fov_x; /* In degrees */
+ double fov_y; /* In degrees */
} camera;
struct {
@@ -45,6 +45,7 @@ struct htrdr_args {
} image;
double main_dir[3];
+ double optical_thickness; /* Threshold used during octree building */
unsigned nthreads; /* Hint on the number of threads to use */
int force_overwriting;
@@ -74,6 +75,7 @@ struct htrdr_args {
1 /* #samples per pixel */ \
}, \
{0, 0, 1}, /* Main direction */ \
+ 1.0, /* Optical thickness */ \
(unsigned)~0, /* #threads */ \
0, /* Force overwriting */ \
0, /* dump VTK */ \
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -91,6 +91,11 @@ is_file_updated
int* is_upd);
extern LOCAL_SYM res_T
+update_file_stamp
+ (struct htrdr* htrdr,
+ const char* filename);
+
+extern LOCAL_SYM res_T
create_directory
(struct htrdr* htrdt,
const char* path);
diff --git a/src/htrdr_camera.c b/src/htrdr_camera.c
@@ -27,8 +27,6 @@ struct htrdr_camera {
double axis_z[3];
double position[3];
- double fov_x; /* Field of view in radians */
- double rcp_proj_ratio; /* height / width */
ref_T ref;
struct htrdr* htrdr;
@@ -79,14 +77,12 @@ htrdr_camera_create
res = RES_BAD_ARG;
goto error;
}
- cam->fov_x = fov;
if(proj_ratio <= 0) {
htrdr_log_err(htrdr, "invalid projection ratio `%g'\n", proj_ratio);
res = RES_BAD_ARG;
goto error;
}
- cam->rcp_proj_ratio = 1.0 / proj_ratio;
if(d3_normalize(z, d3_sub(z, target, position)) <= 0
|| d3_normalize(x, d3_cross(x, z, up)) <= 0
@@ -101,9 +97,9 @@ htrdr_camera_create
goto error;
}
- img_plane_depth = 1.0/tan(fov);
- d3_set(cam->axis_x, x);
- d3_muld(cam->axis_y, y, cam->rcp_proj_ratio);
+ img_plane_depth = 1.0/tan(fov*0.5);
+ d3_muld(cam->axis_x, x, proj_ratio);
+ d3_set(cam->axis_y, y);
d3_muld(cam->axis_z, z, img_plane_depth);
d3_set(cam->position, position);
diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c
@@ -138,6 +138,7 @@ transmissivity_hit_filter
void* context)
{
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;
@@ -146,9 +147,9 @@ transmissivity_hit_filter
(void)range;
k_min = htrdr_sky_fetch_svx_voxel_property(ctx->sky, ctx->prop,
- HTRDR_SVX_MIN, HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, &hit->voxel);
+ HTRDR_SVX_MIN, comp_mask, ctx->iband, ctx->iquad, &hit->voxel);
k_max = htrdr_sky_fetch_svx_voxel_property(ctx->sky, ctx->prop,
- HTRDR_SVX_MAX, HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, &hit->voxel);
+ HTRDR_SVX_MAX, comp_mask, ctx->iband, ctx->iquad, &hit->voxel);
ASSERT(k_min <= k_max);
/* Compute the distance to traverse into the voxel */
@@ -185,10 +186,9 @@ transmissivity_hit_filter
x[2] = org[2] + ctx->traversal_dst * dir[2];
k = htrdr_sky_fetch_raw_property(ctx->sky, ctx->prop,
- HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, x, k_min, k_max);
+ comp_mask, ctx->iband, ctx->iquad, x, k_min, k_max);
ASSERT(k >= k_min && k <= k_max);
- /* Handle the case that k_max is not *really* the max */
proba = (k - k_min) / (k_max - k_min);
if(ssp_rng_canonical(ctx->rng) < proba) { /* Collide */
@@ -218,7 +218,6 @@ transmissivity
{
struct s3d_hit s3d_hit;
struct svx_hit svx_hit;
- struct svx_tree* svx_tree;
struct transmissivity_context transmissivity_ctx = TRANSMISSION_CONTEXT_NULL;
struct s3d_hit s3d_hit_prev = hit_prev ? *hit_prev : S3D_HIT_NULL;
@@ -245,9 +244,8 @@ transmissivity
transmissivity_ctx.prop = prop;
/* Compute the transmissivity */
- svx_tree = htrdr_sky_get_svx_tree(htrdr->sky, iband, iquad);
- SVX(tree_trace_ray(svx_tree, pos, dir, range, NULL,
- transmissivity_hit_filter, &transmissivity_ctx, &svx_hit));
+ HTRDR(sky_trace_ray(htrdr->sky, pos, dir, range, NULL,
+ transmissivity_hit_filter, &transmissivity_ctx, iband, iquad, &svx_hit));
if(SVX_HIT_NONE(&svx_hit)) {
return transmissivity_ctx.Tmin ? exp(-transmissivity_ctx.Tmin) : 1.0;
@@ -272,7 +270,6 @@ htrdr_compute_radiance_sw
struct s3d_hit s3d_hit = S3D_HIT_NULL;
struct svx_hit svx_hit = SVX_HIT_NULL;
struct s3d_hit s3d_hit_prev = S3D_HIT_NULL;
- struct svx_tree* svx_tree = NULL;
struct ssf_phase* phase_hg = NULL;
struct ssf_phase* phase_rayleigh = NULL;
struct ssf_bsdf* bsdf = NULL;
@@ -311,7 +308,7 @@ htrdr_compute_radiance_sw
CHK(RES_OK == ssf_phase_create
(&htrdr->lifo_allocators[ithread], &ssf_phase_rayleigh, &phase_rayleigh));
- SSF(lambertian_reflection_setup(bsdf, 0.02));
+ SSF(lambertian_reflection_setup(bsdf, 0.5));
/* Setup the phase function for this spectral band & quadrature point */
g = htrdr_sky_fetch_particle_phase_function_asymmetry_parameter
@@ -328,8 +325,6 @@ htrdr_compute_radiance_sw
sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
- svx_tree = htrdr_sky_get_svx_tree(htrdr->sky, iband, iquad);
-
d3_set(pos, pos_in);
d3_set(dir, dir_in);
@@ -363,7 +358,7 @@ htrdr_compute_radiance_sw
S3D(scene_view_trace_ray(htrdr->s3d_scn_view, ray_pos, ray_dir, ray_range,
&s3d_hit_prev, &s3d_hit));
- /* Sample an optical thicknes */
+ /* Sample an optical thickness */
scattering_ctx.Ts = ssp_ran_exp(rng, 1);
/* Setup the remaining scattering context fields */
@@ -374,12 +369,11 @@ htrdr_compute_radiance_sw
/* Define if a scattering event occurs */
d2(range, 0, s3d_hit.distance);
- SVX(tree_trace_ray(svx_tree, pos, dir, range, NULL,
- scattering_hit_filter, &scattering_ctx, &svx_hit));
+ HTRDR(sky_trace_ray(htrdr->sky, pos, dir, range, NULL,
+ scattering_hit_filter, &scattering_ctx, iband, iquad, &svx_hit));
/* No scattering and no surface reflection. Stop the radiative random walk */
if(S3D_HIT_NONE(&s3d_hit) && SVX_HIT_NONE(&svx_hit)) {
- w *= ksi;
break;
}
ASSERT(SVX_HIT_NONE(&svx_hit)
diff --git a/src/htrdr_sky.c b/src/htrdr_sky.c
@@ -27,6 +27,7 @@
#include <high_tune/htgop.h>
#include <high_tune/htmie.h>
+#include <rsys/clock_time.h>
#include <rsys/dynamic_array_double.h>
#include <rsys/dynamic_array_size_t.h>
#include <rsys/hash_table.h>
@@ -57,21 +58,22 @@ struct split {
#define DARRAY_DATA struct split
#include <rsys/dynamic_array.h>
-struct build_octree_context {
+struct build_tree_context {
const struct htrdr_sky* sky;
- double vxsz[3]; /* Size of a SVX voxel */
- double dst_max; /* Max traversal distance */
+ double vxsz[3];
double tau_threshold; /* Threshold criteria for the merge process */
size_t iband; /* Index of the band that overlaps the CIE XYZ color space */
+ size_t quadrature_range[2]; /* Range of quadrature point indices to handle */
/* Precomputed voxel data of the finest level. May be NULL <=> compute the
- * voxel data at runtime. */
+ * voxel data at runtime. This structure is not used during the construction
+ * of the binary tree of the atmosphere */
struct htrdr_grid* grid;
};
-#define BUILD_OCTREE_CONTEXT_NULL__ { NULL, {0,0,0}, 0, 0, 0, NULL }
-static const struct build_octree_context BUILD_OCTREE_CONTEXT_NULL =
- BUILD_OCTREE_CONTEXT_NULL__;
+#define BUILD_TREE_CONTEXT_NULL__ {NULL,{0,0,0},0,0,{SIZE_MAX,SIZE_MAX},NULL}
+static const struct build_tree_context BUILD_TREE_CONTEXT_NULL =
+ BUILD_TREE_CONTEXT_NULL__;
struct vertex {
double x;
@@ -115,8 +117,8 @@ struct sw_band_prop {
double g_avg;
};
-/* Encompass the hierarchical data structure of the cloud data and its
- * associated descriptor.
+/* Encompass the hierarchical data structure of the cloud/atmospheric data and
+ * its associated descriptor.
*
* For each SVX voxel, the data of the optical property are stored
* linearly as N single precision floating point data, with N computed as
@@ -137,9 +139,16 @@ struct cloud {
struct svx_tree* octree;
struct svx_tree_desc octree_desc;
};
+struct atmosphere {
+ struct svx_tree* bitree;
+ struct svx_tree_desc bitree_desc;
+};
struct htrdr_sky {
- struct cloud* clouds; /* Per sw_band cloud data structure */
+ struct cloud** clouds; /* Per sw_band cloud data structure */
+
+ /* Per sw_band and per quadrature point atmosphere data structure */
+ struct atmosphere** atmosphere;
struct htrdr_sun* sun; /* Sun attached to the sky */
@@ -216,6 +225,39 @@ cloud_water_vapor_molar_fraction
return rvt / (rvt + H2O_MOLAR_MASS/DRY_AIR_MOLAR_MASS);
}
+/* Smits intersection: "Efficiency issues for ray tracing" */
+static FINLINE void
+ray_intersect_aabb
+ (const double org[3],
+ const double dir[3],
+ const double range[2],
+ const double low[3],
+ const double upp[3],
+ double hit_range[2])
+{
+ double hit[2];
+ int i;
+ ASSERT(org && dir && range && low && upp && hit_range);
+ ASSERT(low[0] < upp[0]);
+ ASSERT(low[1] < upp[1]);
+ ASSERT(low[2] < upp[2]);
+
+ hit_range[0] = INF;
+ hit_range[1] =-INF;
+ hit[0] = range[0];
+ hit[1] = range[1];
+ FOR_EACH(i, 0, 3) {
+ double t_min = (low[i] - org[i]) / dir[i];
+ double t_max = (upp[i] - org[i]) / dir[i];
+ if(t_min > t_max) SWAP(double, t_min, t_max);
+ hit[0] = MMAX(t_min, hit[0]);
+ hit[1] = MMIN(t_max, hit[1]);
+ if(hit[0] > hit[1]) return;
+ }
+ hit_range[0] = hit[0];
+ hit_range[1] = hit[1];
+}
+
static INLINE void
octree_data_init
(const struct htrdr_sky* sky,
@@ -296,10 +338,72 @@ register_leaf
}
static void
-vox_get_particle
+clean_clouds(struct htrdr_sky* sky)
+{
+ size_t nbands;
+ size_t i;
+ ASSERT(sky);
+
+ if(!sky->clouds) return;
+
+ nbands = htrdr_sky_get_sw_spectral_bands_count(sky);
+ FOR_EACH(i, 0, nbands) {
+ struct htgop_spectral_interval band;
+ size_t iband;
+ size_t iquad;
+
+ iband = sky->sw_bands_range[0] + i;
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+
+ if(!sky->clouds[i]) continue;
+
+ FOR_EACH(iquad, 0, band.quadrature_length) {
+ if(sky->clouds[i][iquad].octree) {
+ SVX(tree_ref_put(sky->clouds[i][iquad].octree));
+ sky->clouds[i][iquad].octree = NULL;
+ }
+ }
+ MEM_RM(sky->htrdr->allocator, sky->clouds[i]);
+ }
+ MEM_RM(sky->htrdr->allocator, sky->clouds);
+}
+
+static void
+clean_atmosphere(struct htrdr_sky* sky)
+{
+ size_t nbands;
+ size_t i;
+ ASSERT(sky);
+
+ if(!sky->atmosphere) return;
+
+ nbands = htrdr_sky_get_sw_spectral_bands_count(sky);
+ FOR_EACH(i, 0, nbands) {
+ struct htgop_spectral_interval band;
+ size_t iband;
+ size_t iquad;
+
+ iband = sky->sw_bands_range[0] + i;
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+
+ if(!sky->atmosphere[i]) continue;
+
+ FOR_EACH(iquad, 0, band.quadrature_length) {
+ if(sky->atmosphere[i][iquad].bitree) {
+ SVX(tree_ref_put(sky->atmosphere[i][iquad].bitree));
+ sky->atmosphere[i][iquad].bitree = NULL;
+ }
+ }
+ MEM_RM(sky->htrdr->allocator, sky->atmosphere[i]);
+ }
+ MEM_RM(sky->htrdr->allocator, sky->atmosphere);
+}
+
+static void
+cloud_vox_get_particle
(const size_t xyz[3],
float vox[],
- const struct build_octree_context* ctx)
+ const struct build_tree_context* ctx)
{
double rct;
double ka, ks, kext;
@@ -411,16 +515,15 @@ vox_get_particle
}
static void
-vox_get_gas
+cloud_vox_get_gas
(const size_t xyz[3],
float vox[],
- const struct build_octree_context* ctx)
+ const struct build_tree_context* ctx)
{
struct htgop_layer layer;
struct htgop_layer_sw_spectral_interval band;
size_t ilayer;
size_t layer_range[2];
- size_t quad_range[2];
double x_h2o_range[2];
double vox_low[3], vox_upp[3]; /* Voxel AABB */
double ka_min, ka_max;
@@ -519,20 +622,20 @@ vox_get_gas
/* ... retrieve the considered spectral interval */
HTGOP(layer_get_sw_spectral_interval(&layer, ctx->iband, &band));
- quad_range[0] = 0;
- quad_range[1] = band.quadrature_length-1;
+ ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]);
+ ASSERT(ctx->quadrature_range[1] < band.quadrature_length);
/* ... and compute the radiative properties and upd their bounds */
HTGOP(layer_sw_spectral_interval_quadpoints_get_ka_bounds
- (&band, quad_range, x_h2o_range, k));
+ (&band, ctx->quadrature_range, x_h2o_range, k));
ka_min = MMIN(ka_min, k[0]);
ka_max = MMAX(ka_max, k[1]);
HTGOP(layer_sw_spectral_interval_quadpoints_get_ks_bounds
- (&band, quad_range, x_h2o_range, k));
+ (&band, ctx->quadrature_range, x_h2o_range, k));
ks_min = MMIN(ks_min, k[0]);
ks_max = MMAX(ks_max, k[1]);
HTGOP(layer_sw_spectral_interval_quadpoints_get_kext_bounds
- (&band, quad_range, x_h2o_range, k));
+ (&band, ctx->quadrature_range, x_h2o_range, k));
kext_min = MMIN(kext_min, k[0]);
kext_max = MMAX(kext_max, k[1]);
}
@@ -555,9 +658,9 @@ vox_get_gas
}
static void
-vox_get(const size_t xyz[3], void* dst, void* context)
+cloud_vox_get(const size_t xyz[3], void* dst, void* context)
{
- struct build_octree_context* ctx = context;
+ struct build_tree_context* ctx = context;
ASSERT(context);
if(ctx->grid) { /* Fetch voxel data from precomputed grid */
@@ -565,8 +668,8 @@ vox_get(const size_t xyz[3], void* dst, void* context)
memcpy(dst, vox_data, NFLOATS_PER_VOXEL * sizeof(float));
} else {
/* No precomputed grid. Compute the voxel data at runtime */
- vox_get_particle(xyz, dst, ctx);
- vox_get_gas(xyz, dst, ctx);
+ cloud_vox_get_particle(xyz, dst, ctx);
+ cloud_vox_get_gas(xyz, dst, ctx);
}
}
@@ -606,7 +709,7 @@ vox_merge_component
}
static void
-vox_merge(void* dst, const void* voxels[], const size_t nvoxs, void* context)
+cloud_vox_merge(void* dst, const void* voxels[], const size_t nvoxs, void* context)
{
ASSERT(dst && voxels && nvoxs);
(void)context;
@@ -617,31 +720,36 @@ vox_merge(void* dst, const void* voxels[], const size_t nvoxs, void* context)
static INLINE int
vox_challenge_merge_component
(const enum htrdr_sky_component comp,
- const float* voxels[],
+ const struct svx_voxel voxels[],
const size_t nvoxs,
- struct build_octree_context* ctx)
+ struct build_tree_context* ctx)
{
+ double lower_z = DBL_MAX;
+ double upper_z =-DBL_MAX;
+ double dst;
float kext_min = FLT_MAX;
float kext_max =-FLT_MAX;
size_t ivox;
ASSERT(voxels && nvoxs && ctx);
FOR_EACH(ivox, 0, nvoxs) {
- const float* vox = voxels[ivox];
+ const float* vox = voxels[ivox].data;
kext_min = MMIN(kext_min, voxel_get(vox, comp, HTRDR_Kext, HTRDR_SVX_MIN));
kext_max = MMAX(kext_max, voxel_get(vox, comp, HTRDR_Kext, HTRDR_SVX_MAX));
+ lower_z = MMIN(voxels[ivox].lower[2], lower_z);
+ upper_z = MMAX(voxels[ivox].upper[2], upper_z);
}
- return (kext_max - kext_min)*ctx->dst_max <= ctx->tau_threshold;
+ dst = upper_z - lower_z;
+ return (kext_max - kext_min)*dst <= ctx->tau_threshold;
}
static int
-vox_challenge_merge
- (const void* voxels[], const size_t nvoxs, void* ctx)
+cloud_vox_challenge_merge
+ (const struct svx_voxel voxels[], const size_t nvoxs, void* ctx)
{
- const float** voxs = (const float**)voxels;
ASSERT(voxels);
- return vox_challenge_merge_component(HTRDR_PARTICLES__, voxs, nvoxs, ctx)
- && vox_challenge_merge_component(HTRDR_GAS__, voxs, nvoxs, ctx);
+ return vox_challenge_merge_component(HTRDR_PARTICLES__, voxels, nvoxs, ctx)
+ && vox_challenge_merge_component(HTRDR_GAS__, voxels, nvoxs, ctx);
}
static res_T
@@ -702,6 +810,7 @@ setup_cloud_grid
(struct htrdr_sky* sky,
const size_t definition[3],
const size_t iband,
+ const size_t iquad,
const char* htcp_filename,
const char* htgop_filename,
const char* htmie_filename,
@@ -711,7 +820,8 @@ setup_cloud_grid
struct htrdr_grid* grid = NULL;
struct str path;
struct str str;
- struct build_octree_context ctx = BUILD_OCTREE_CONTEXT_NULL;
+ struct build_tree_context ctx = BUILD_TREE_CONTEXT_NULL;
+ struct htgop_spectral_interval band;
size_t sizeof_cell;
size_t ncells;
uint64_t mcode;
@@ -727,7 +837,7 @@ setup_cloud_grid
str_init(sky->htrdr->allocator, &str);
str_init(sky->htrdr->allocator, &path);
- CHK((size_t)snprintf(buf, sizeof(buf), ".%lu", iband) < sizeof(buf));
+ CHK((size_t)snprintf(buf, sizeof(buf), ".%lu.%lu", iband, iquad) < sizeof(buf));
res = setup_temp_dir(sky, htgop_filename, htmie_filename, &path);
if(res != RES_OK) goto error;
@@ -738,7 +848,7 @@ setup_cloud_grid
CHK(RES_OK == str_append(&path, str_cget(&str)));
CHK(RES_OK == str_append(&path, ".grid"));
CHK(RES_OK == str_append(&path, buf));
-
+
if(!force_update) {
/* Try to open an already saved grid */
res = htrdr_grid_open(sky->htrdr, str_cget(&path), &grid);
@@ -773,10 +883,13 @@ setup_cloud_grid
if(res != RES_OK) goto error;
ctx.sky = sky;
- ctx.dst_max = DBL_MAX; /* Unused for grid construction */
ctx.tau_threshold = DBL_MAX; /* Unused for grid construction */
ctx.iband = iband;
+ HTGOP(get_sw_spectral_interval(sky->htgop, ctx.iband, &band));
+ ctx.quadrature_range[0] = iquad;
+ ctx.quadrature_range[1] = iquad;
+
/* Compute the size of a SVX voxel */
ctx.vxsz[0] = sky->htcp_desc.upper[0] - sky->htcp_desc.lower[0];
ctx.vxsz[1] = sky->htcp_desc.upper[1] - sky->htcp_desc.lower[1];
@@ -797,7 +910,7 @@ setup_cloud_grid
/* Compute the overall number of voxels in the htrdr_grid */
ncells = definition[0] * definition[1] * definition[2];
- fprintf(stderr, "Generating cloud grid %lu: %3u%%", iband, 0);
+ fprintf(stderr, "Generating cloud grid %lu.%lu: %3u%%", iband, iquad, 0);
fflush(stderr);
omp_set_num_threads((int)sky->htrdr->nthreads);
@@ -813,7 +926,7 @@ setup_cloud_grid
if((xyz[2] = morton3D_decode_u21(mcode >> 0)) >= definition[2]) continue;
/* Compute the voxel data */
- vox_get(xyz, htrdr_grid_at_mcode(grid, mcode), &ctx);
+ cloud_vox_get(xyz, htrdr_grid_at_mcode(grid, mcode), &ctx);
/* Update the progress message */
n = (size_t)ATOMIC_INCR(&ncells_computed);
@@ -821,8 +934,8 @@ setup_cloud_grid
#pragma omp critical
if(pcent > progress) {
progress = pcent;
- fprintf(stderr, "%c[2K\rGenerating cloud grid %lu: %3u%%",
- 27, iband, (unsigned)pcent);
+ fprintf(stderr, "%c[2K\rGenerating cloud grid %lu.%lu: %3u%%",
+ 27, iband, iquad, (unsigned)pcent);
fflush(stderr);
}
}
@@ -848,18 +961,18 @@ setup_clouds
const char* htcp_filename,
const char* htgop_filename,
const char* htmie_filename,
+ const double optical_thickness_threshold,
const int force_cache_update)
{
struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL;
- struct build_octree_context ctx = BUILD_OCTREE_CONTEXT_NULL;
+ struct build_tree_context ctx = BUILD_TREE_CONTEXT_NULL;
size_t nvoxs[3];
double low[3];
double upp[3];
- double sz[3];
size_t nbands;
size_t i;
res_T res = RES_OK;
- ASSERT(sky && sky->sw_bands);
+ ASSERT(sky && sky->sw_bands && optical_thickness_threshold >= 0);
res = htcp_get_desc(sky->htcp, &sky->htcp_desc);
if(res != RES_OK) {
@@ -919,20 +1032,20 @@ setup_clouds
ASSERT(eq_eps(upp[2] - low[2], len_z, 1.e-6));
}
- /* Compute the size of of the AABB */
- sz[0] = upp[0] - low[0];
- sz[1] = upp[1] - low[1];
- sz[2] = upp[2] - low[2];
-
/* Setup the build context */
ctx.sky = sky;
- ctx.dst_max = sz[2];
- ctx.tau_threshold = 0.1;
+ ctx.tau_threshold = optical_thickness_threshold;
+ ctx.vxsz[0] = sky->htcp_desc.upper[0] - sky->htcp_desc.lower[0];
+ ctx.vxsz[1] = sky->htcp_desc.upper[1] - sky->htcp_desc.lower[1];
+ ctx.vxsz[2] = sky->htcp_desc.upper[2] - sky->htcp_desc.lower[2];
+ ctx.vxsz[0] = ctx.vxsz[0] / (double)sky->htcp_desc.spatial_definition[0];
+ ctx.vxsz[1] = ctx.vxsz[1] / (double)sky->htcp_desc.spatial_definition[1];
+ ctx.vxsz[2] = ctx.vxsz[2] / (double)sky->htcp_desc.spatial_definition[2];
/* Setup the voxel descriptor */
- vox_desc.get = vox_get;
- vox_desc.merge = vox_merge;
- vox_desc.challenge_merge = vox_challenge_merge;
+ vox_desc.get = cloud_vox_get;
+ vox_desc.merge = cloud_vox_merge;
+ vox_desc.challenge_merge = cloud_vox_challenge_merge;
vox_desc.context = &ctx;
vox_desc.size = sizeof(float) * NFLOATS_PER_VOXEL;
@@ -947,29 +1060,50 @@ setup_clouds
}
FOR_EACH(i, 0, nbands) {
+ size_t iquad;
+ struct htgop_spectral_interval band;
ctx.iband = i + sky->sw_bands_range[0];
- /* Compute grid of voxels data */
- res = setup_cloud_grid(sky, nvoxs, ctx.iband, htcp_filename,
- htgop_filename, htmie_filename, force_cache_update, &ctx.grid);
- if(res != RES_OK) goto error;
+ HTGOP(get_sw_spectral_interval(sky->htgop, ctx.iband, &band));
- /* Create the octree */
- res = svx_octree_create
- (sky->htrdr->svx, low, upp, nvoxs, &vox_desc, &sky->clouds[i].octree);
- if(res != RES_OK) {
- htrdr_log_err(sky->htrdr, "could not create the octree of the cloud "
- "properties for the band %lu.\n", (unsigned long)ctx.iband);
+ sky->clouds[i] = MEM_CALLOC(sky->htrdr->allocator,
+ band.quadrature_length, sizeof(*sky->clouds[i]));
+ if(!sky->clouds[i]) {
+ htrdr_log_err(sky->htrdr,
+ "could not create the list of per quadrature point cloud data "
+ "for the band %lu.\n", (unsigned long)ctx.iband);
+ res = RES_MEM_ERR;
goto error;
}
- /* Fetch the octree descriptor for future use */
- SVX(tree_get_desc
- (sky->clouds[i].octree, &sky->clouds[i].octree_desc));
+ /* Build a cloud octree for each quadrature point of the considered
+ * spectral band */
+ FOR_EACH(iquad, 0, band.quadrature_length) {
+ ctx.quadrature_range[0] = iquad;
+ ctx.quadrature_range[1] = iquad;
+
+ /* Compute grid of voxels data */
+ res = setup_cloud_grid(sky, nvoxs, ctx.iband, iquad, htcp_filename,
+ htgop_filename, htmie_filename, force_cache_update, &ctx.grid);
+ if(res != RES_OK) goto error;
+
+ /* Create the octree */
+ res = svx_octree_create(sky->htrdr->svx, low, upp, nvoxs, &vox_desc,
+ &sky->clouds[i][iquad].octree);
+ if(res != RES_OK) {
+ htrdr_log_err(sky->htrdr, "could not create the octree of the cloud "
+ "properties for the band %lu.\n", (unsigned long)ctx.iband);
+ goto error;
+ }
+
+ /* Fetch the octree descriptor for future use */
+ SVX(tree_get_desc
+ (sky->clouds[i][iquad].octree, &sky->clouds[i][iquad].octree_desc));
- if(ctx.grid) {
- htrdr_grid_ref_put(ctx.grid);
- ctx.grid = NULL;
+ if(ctx.grid) {
+ htrdr_grid_ref_put(ctx.grid);
+ ctx.grid = NULL;
+ }
}
}
@@ -977,16 +1111,197 @@ exit:
return res;
error:
if(ctx.grid) htrdr_grid_ref_put(ctx.grid);
- if(sky->clouds) {
- FOR_EACH(i, 0, nbands) {
- if(sky->clouds[i].octree) {
- SVX(tree_ref_put(sky->clouds[i].octree));
- sky->clouds[i].octree = NULL;
+ clean_clouds(sky);
+ darray_split_clear(&sky->svx2htcp_z);
+ goto exit;
+}
+
+static void
+atmosphere_vox_get(const size_t xyz[3], void* dst, void* context)
+{
+ float* vox = dst;
+ struct build_tree_context* ctx = context;
+ struct htgop_level level;
+ size_t layer_range[2];
+ size_t nlevels;
+ double vox_low, vox_upp;
+ double ka_min, ks_min, kext_min;
+ double ka_max, ks_max, kext_max;
+ size_t ilayer;
+ ASSERT(xyz && dst && context);
+
+ /* Compute the boundaries of the SVX voxel */
+ HTGOP(get_level(ctx->sky->htgop, 0, &level));
+ vox_low = (double)xyz[2] * ctx->vxsz[2] + level.height;
+ HTGOP(get_levels_count(ctx->sky->htgop, &nlevels));
+ HTGOP(get_level(ctx->sky->htgop, nlevels-1, &level));
+ vox_upp = MMIN(vox_low + ctx->vxsz[2], level.height); /* Handle numerical issues */
+
+ /* Define the atmospheric layers overlapped by the SVX voxel */
+ HTGOP(position_to_layer_id(ctx->sky->htgop, vox_low, &layer_range[0]));
+ HTGOP(position_to_layer_id(ctx->sky->htgop, vox_upp, &layer_range[1]));
+
+ ka_min = ks_min = kext_min = DBL_MAX;
+ ka_max = ks_max = kext_max =-DBL_MAX;
+
+ /* For each atmospheric layer that overlaps the SVX voxel ... */
+ FOR_EACH(ilayer, layer_range[0], layer_range[1]+1) {
+ struct htgop_layer layer;
+ struct htgop_layer_sw_spectral_interval band;
+ size_t iquad;
+
+ HTGOP(get_layer(ctx->sky->htgop, ilayer, &layer));
+
+ /* ... retrieve the considered spectral interval */
+ HTGOP(layer_get_sw_spectral_interval(&layer, ctx->iband, &band));
+
+ /* ... and update the radiative properties bound with the per quadrature
+ * point nominal values */
+ ASSERT(ctx->quadrature_range[0] <= ctx->quadrature_range[1]);
+ ASSERT(ctx->quadrature_range[1] < band.quadrature_length);
+ FOR_EACH(iquad, ctx->quadrature_range[0], ctx->quadrature_range[1]+1) {
+ ka_min = MMIN(ka_min, band.ka_nominal[iquad]);
+ ka_max = MMAX(ka_max, band.ka_nominal[iquad]);
+ ks_min = MMIN(ks_min, band.ks_nominal[iquad]);
+ ks_max = MMAX(ks_max, band.ks_nominal[iquad]);
+ kext_min = MMIN(kext_min, band.ka_nominal[iquad]+band.ks_nominal[iquad]);
+ kext_max = MMAX(kext_max, band.ka_nominal[iquad]+band.ks_nominal[iquad]);
+ }
+ }
+
+ /* Ensure that the single precision bounds include their double precision
+ * version. */
+ if(ka_min != (float)ka_min) ka_min = nextafterf((float)ka_min,-FLT_MAX);
+ if(ka_max != (float)ka_max) ka_max = nextafterf((float)ka_max, FLT_MAX);
+ if(ks_min != (float)ks_min) ks_min = nextafterf((float)ks_min,-FLT_MAX);
+ if(ks_max != (float)ks_max) ks_max = nextafterf((float)ks_max, FLT_MAX);
+ if(kext_min != (float)kext_min) kext_min = nextafterf((float)kext_min,-FLT_MAX);
+ if(kext_max != (float)kext_max) kext_max = nextafterf((float)kext_max, FLT_MAX);
+
+ /* Setup gas optical properties of the voxel */
+ voxel_set(vox, HTRDR_GAS__, HTRDR_Ka, HTRDR_SVX_MIN, (float)ka_min);
+ voxel_set(vox, HTRDR_GAS__, HTRDR_Ka, HTRDR_SVX_MAX, (float)ka_max);
+ voxel_set(vox, HTRDR_GAS__, HTRDR_Ks, HTRDR_SVX_MIN, (float)ks_min);
+ voxel_set(vox, HTRDR_GAS__, HTRDR_Ks, HTRDR_SVX_MAX, (float)ks_max);
+ voxel_set(vox, HTRDR_GAS__, HTRDR_Kext, HTRDR_SVX_MIN, (float)kext_min);
+ voxel_set(vox, HTRDR_GAS__, HTRDR_Kext, HTRDR_SVX_MAX, (float)kext_max);
+}
+
+static void
+atmosphere_vox_merge
+ (void* dst, const void* voxels[], const size_t nvoxs, void* context)
+{
+ ASSERT(dst && voxels && nvoxs);
+ (void)context;
+ vox_merge_component(dst, HTRDR_GAS__, (const float**)voxels, nvoxs);
+}
+
+static int
+atmosphere_vox_challenge_merge
+ (const struct svx_voxel voxels[], const size_t nvoxs, void* ctx)
+{
+ ASSERT(voxels);
+ return vox_challenge_merge_component(HTRDR_GAS__, voxels, nvoxs, ctx);
+}
+
+static res_T
+setup_atmosphere
+ (struct htrdr_sky* sky, const double optical_thickness_threshold)
+{
+ struct build_tree_context ctx = BUILD_TREE_CONTEXT_NULL;
+ struct htgop_level lvl_low, lvl_upp;
+ struct svx_voxel_desc vox_desc = SVX_VOXEL_DESC_NULL;
+ double low, upp;
+ size_t nlayers, nlevels;
+ size_t definition;
+ size_t nbands;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(sky && optical_thickness_threshold >= 0);
+
+ HTGOP(get_layers_count(sky->htgop, &nlayers));
+ HTGOP(get_levels_count(sky->htgop, &nlevels));
+ HTGOP(get_level(sky->htgop, 0, &lvl_low));
+ HTGOP(get_level(sky->htgop, nlevels-1, &lvl_upp));
+ low = lvl_low.height;
+ upp = lvl_upp.height;
+ definition = nlayers;
+
+ /* Setup the build context */
+ ctx.sky = sky;
+ ctx.tau_threshold = optical_thickness_threshold;
+ ctx.vxsz[0] = INF;
+ ctx.vxsz[1] = INF;
+ ctx.vxsz[2] = (upp-low)/(double)definition;
+
+ /* Setup the voxel descriptor for the atmosphere. Note that in contrats with
+ * the clouds, the voxel contains only NFLOATS_PER_COMPONENT floats and not
+ * NFLOATS_PER_VOXEL. Indeed, the atmosphere has only 1 component: the gas.
+ * Anyway, we still rely on the memory layout of the cloud voxels excepted
+ * that we assume that the optical properties of the particles are never
+ * fetched. We thus have to ensure that the gas properties are arranged
+ * before the particles, i.e. HTRDR_GAS__ == 0 */
+ vox_desc.get = atmosphere_vox_get;
+ vox_desc.merge = atmosphere_vox_merge;
+ vox_desc.challenge_merge = atmosphere_vox_challenge_merge;
+ vox_desc.context = &ctx;
+ vox_desc.size = sizeof(float) * NFLOATS_PER_COMPONENT;
+ STATIC_ASSERT(HTRDR_GAS__ == 0, Unexpected_enum_value);
+
+ /* Create as many atmospheric data structure than considered SW spectral
+ * bands */
+ nbands = htrdr_sky_get_sw_spectral_bands_count(sky);
+ sky->atmosphere = MEM_CALLOC
+ (sky->htrdr->allocator, nbands, sizeof(*sky->atmosphere));
+ if(!sky->atmosphere) {
+ htrdr_log_err(sky->htrdr,
+ "could not create the list of per SW band atmospheric data structure.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ FOR_EACH(i, 0, nbands) {
+ size_t iquad;
+ struct htgop_spectral_interval band;
+ ctx.iband = i + sky->sw_bands_range[0];
+
+ HTGOP(get_sw_spectral_interval(sky->htgop, ctx.iband, &band));
+
+ sky->atmosphere[i] = MEM_CALLOC(sky->htrdr->allocator,
+ band.quadrature_length, sizeof(*sky->atmosphere[i]));
+ if(!sky->atmosphere[i]) {
+ htrdr_log_err(sky->htrdr,
+ "could not create the list of per quadrature point atmospheric data "
+ "for the band %lu.\n", (unsigned long)ctx.iband);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Build an atmospheric binary tree for each quadrature point of the
+ * considered spectral band */
+ FOR_EACH(iquad, 0, band.quadrature_length) {
+ ctx.quadrature_range[0] = iquad;
+ ctx.quadrature_range[1] = iquad;
+
+ /* Create the atmospheric binary tree */
+ res = svx_bintree_create(sky->htrdr->svx, low, upp, definition, SVX_AXIS_Z,
+ &vox_desc, &sky->atmosphere[i][iquad].bitree);
+ if(res != RES_OK) {
+ htrdr_log_err(sky->htrdr, "could not create the binary tree of the "
+ "atmospheric properties for the band %lu.\n", (unsigned long)ctx.iband);
+ goto error;
}
+
+ /* Fetch the binary tree descriptor for future use */
+ SVX(tree_get_desc(sky->atmosphere[i][iquad].bitree,
+ &sky->atmosphere[i][iquad].bitree_desc));
}
- MEM_RM(sky->htrdr->allocator, sky->clouds);
}
- darray_split_clear(&sky->svx2htcp_z);
+
+exit:
+ return res;
+error:
+ clean_atmosphere(sky);
goto exit;
}
@@ -1071,22 +1386,13 @@ release_sky(ref_T* ref)
struct htrdr_sky* sky;
ASSERT(ref);
sky = CONTAINER_OF(ref, struct htrdr_sky, ref);
+ clean_clouds(sky);
+ clean_atmosphere(sky);
if(sky->sun) htrdr_sun_ref_put(sky->sun);
if(sky->htcp) HTCP(ref_put(sky->htcp));
if(sky->htgop) HTGOP(ref_put(sky->htgop));
if(sky->htmie) HTMIE(ref_put(sky->htmie));
if(sky->sw_bands) MEM_RM(sky->htrdr->allocator, sky->sw_bands);
- if(sky->clouds) {
- const size_t nbands = htrdr_sky_get_sw_spectral_bands_count(sky);
- size_t i;
- FOR_EACH(i, 0, nbands) {
- if(sky->clouds[i].octree) {
- SVX(tree_ref_put(sky->clouds[i].octree));
- sky->clouds[i].octree = NULL;
- }
- }
- MEM_RM(sky->htrdr->allocator, sky->clouds);
- }
darray_split_release(&sky->svx2htcp_z);
MEM_RM(sky->htrdr->allocator, sky);
}
@@ -1101,15 +1407,19 @@ htrdr_sky_create
const char* htcp_filename,
const char* htgop_filename,
const char* htmie_filename,
+ const double optical_thickness_threshold,
struct htrdr_sky** out_sky)
{
+ struct time t0, t1;
struct htrdr_sky* sky = NULL;
+ char buf[128];
int htcp_upd = 1;
int htmie_upd = 1;
int htgop_upd = 1;
int force_upd = 1;
res_T res = RES_OK;
ASSERT(htrdr && sun && htcp_filename && htmie_filename && out_sky);
+ ASSERT(optical_thickness_threshold >= 0);
sky = MEM_CALLOC(htrdr->allocator, 1, sizeof(*sky));
if(!sky) {
@@ -1176,9 +1486,34 @@ htrdr_sky_create
if(res != RES_OK) goto error;
force_upd = htcp_upd || htmie_upd || htgop_upd;
- res = setup_clouds
- (sky, htcp_filename, htgop_filename, htmie_filename, force_upd);
+ time_current(&t0);
+ res = setup_clouds(sky, htcp_filename, htgop_filename, htmie_filename,
+ optical_thickness_threshold, force_upd);
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, "Setup clouds in %s\n", buf);
+
+ /* Update the file stamps */
+ if(htcp_upd) {
+ res = update_file_stamp(sky->htrdr, htcp_filename);
+ if(res != RES_OK) goto error;
+ }
+ if(htmie_upd) {
+ res = update_file_stamp(sky->htrdr, htmie_filename);
+ if(res != RES_OK) goto error;
+ }
+ if(htgop_upd) {
+ res = update_file_stamp(sky->htrdr, htgop_filename);
+ if(res != RES_OK) goto error;
+ }
+
+ time_current(&t0);
+ res = setup_atmosphere(sky, optical_thickness_threshold);
+ 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, "Setup atmosphere in %s\n", buf);
exit:
*out_sky = sky;
@@ -1234,7 +1569,10 @@ htrdr_sky_fetch_raw_property
size_t ivox[3];
size_t i;
const struct svx_tree_desc* cloud_desc;
+ const struct svx_tree_desc* atmosphere_desc;
int comp_mask = components_mask;
+ int in_clouds; /* Defines if `pos' lies in the clouds */
+ int in_atmosphere; /* Defines if `pos' lies in the atmosphere */
double k_particle = 0;
double k_gas = 0;
double k = 0;
@@ -1244,31 +1582,49 @@ htrdr_sky_fetch_raw_property
ASSERT(comp_mask & HTRDR_ALL_COMPONENTS);
i = iband - sky->sw_bands_range[0];
- cloud_desc = &sky->clouds[i].octree_desc;
-
- /* Is the position outside the clouds? */
- if(pos[0] < cloud_desc->lower[0]
- || pos[1] < cloud_desc->lower[1]
- || pos[2] < cloud_desc->lower[2]
- || pos[0] > cloud_desc->upper[0]
- || pos[1] > cloud_desc->upper[1]
- || pos[2] > cloud_desc->upper[2]) {
- comp_mask &= ~HTRDR_PARTICLES; /* No particle */
- }
-
- /* Compute the index of the voxel to fetch */
- ivox[0] = (size_t)((pos[0] - cloud_desc->lower[0])/sky->htcp_desc.vxsz_x);
- ivox[1] = (size_t)((pos[1] - cloud_desc->lower[1])/sky->htcp_desc.vxsz_y);
- if(!sky->htcp_desc.irregular_z) {
- /* The voxels along the Z dimension have the same size */
- ivox[2] = (size_t)((pos[2] - cloud_desc->lower[2])/sky->htcp_desc.vxsz_z[0]);
+ cloud_desc = &sky->clouds[i][iquad].octree_desc;
+ atmosphere_desc = &sky->atmosphere[i][iquad].bitree_desc;
+ ASSERT(atmosphere_desc->frame[0] == SVX_AXIS_Z);
+
+ /* Is the position inside the clouds? */
+ in_clouds =
+ pos[0] >= cloud_desc->lower[0]
+ && pos[1] >= cloud_desc->lower[1]
+ && pos[2] >= cloud_desc->lower[2]
+ && pos[0] <= cloud_desc->upper[0]
+ && pos[1] <= cloud_desc->upper[1]
+ && pos[2] <= cloud_desc->upper[2];
+
+ /* Is the position inside the atmosphere? */
+ in_atmosphere =
+ pos[2] >= atmosphere_desc->lower[2]
+ && pos[2] <= atmosphere_desc->upper[2];
+
+ if(!in_clouds) {
+ /* Make invalid the voxel index */
+ ivox[0] = SIZE_MAX;
+ ivox[1] = SIZE_MAX;
+ ivox[2] = SIZE_MAX;
+ /* Not in clouds => No particle */
+ comp_mask &= ~HTRDR_PARTICLES;
+ /* Not in atmopshere => No gas */
+ if(!in_atmosphere) comp_mask &= ~HTRDR_GAS;
} else {
- /* Irregular voxel size along the Z dimension. Compute the index of the Z
- * position in the svx2htcp_z Look Up Table and use the LUT to define the
- * voxel index into the HTCP descripptor */
- const struct split* splits = darray_split_cdata_get(&sky->svx2htcp_z);
- const size_t ilut = (size_t)((pos[2] - cloud_desc->lower[2]) / sky->lut_cell_sz);
- ivox[2] = splits[ilut].index + (pos[2] > splits[ilut].height);
+ /* Compute the index of the voxel to fetch */
+ ivox[0] = (size_t)((pos[0] - cloud_desc->lower[0])/sky->htcp_desc.vxsz_x);
+ ivox[1] = (size_t)((pos[1] - cloud_desc->lower[1])/sky->htcp_desc.vxsz_y);
+ if(!sky->htcp_desc.irregular_z) {
+ /* The voxels along the Z dimension have the same size */
+ ivox[2] = (size_t)((pos[2] - cloud_desc->lower[2])/sky->htcp_desc.vxsz_z[0]);
+ } else {
+ /* Irregular voxel size along the Z dimension. Compute the index of the Z
+ * position in the svx2htcp_z Look Up Table and use the LUT to define the
+ * voxel index into the HTCP descripptor */
+ const struct split* splits = darray_split_cdata_get(&sky->svx2htcp_z);
+ const size_t ilut = (size_t)
+ ((pos[2] - cloud_desc->lower[2]) / sky->lut_cell_sz);
+ ivox[2] = splits[ilut].index + (pos[2] > splits[ilut].height);
+ }
}
if(comp_mask & HTRDR_PARTICLES) {
@@ -1277,6 +1633,7 @@ htrdr_sky_fetch_raw_property
double ql = 0; /* Droplet density In kg.m^-3 */
double Ca = 0; /* Massic absorption cross section in m^2.kg^-1 */
double Cs = 0; /* Massic scattering cross section in m^2.kg^-1 */
+ ASSERT(in_clouds);
/* Compute he dry air density */
rho_da = cloud_dry_air_density(&sky->htcp_desc, ivox);
@@ -1295,29 +1652,49 @@ htrdr_sky_fetch_raw_property
if(comp_mask & HTRDR_GAS) {
struct htgop_layer layer;
struct htgop_layer_sw_spectral_interval band;
- struct htgop_layer_sw_spectral_interval_tab tab;
- const double x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
size_t ilayer = 0;
+ ASSERT(in_atmosphere);
/* Retrieve the quadrature point into the spectral band of the layer into
* which `pos' lies */
HTGOP(position_to_layer_id(sky->htgop, pos[2], &ilayer));
HTGOP(get_layer(sky->htgop, ilayer, &layer));
HTGOP(layer_get_sw_spectral_interval(&layer, iband, &band));
- HTGOP(layer_sw_spectral_interval_get_tab(&band, iquad, &tab));
-
- /* Fetch the optical properties wrt x_h2o */
- switch(prop) {
- case HTRDR_Ka:
- HTGOP(layer_sw_spectral_interval_tab_fetch_ka(&tab, x_h2o, &k_gas));
- break;
- case HTRDR_Ks:
- HTGOP(layer_sw_spectral_interval_tab_fetch_ks(&tab, x_h2o, &k_gas));
- break;
- case HTRDR_Kext:
- HTGOP(layer_sw_spectral_interval_tab_fetch_kext(&tab, x_h2o, &k_gas));
- break;
- default: FATAL("Unreachable code.\n"); break;
+
+ if(!in_clouds) {
+ /* Pos is outside the clouds. Directly fetch the nominal optical
+ * properties */
+ ASSERT(iquad < band.quadrature_length);
+ switch(prop) {
+ case HTRDR_Ka: k_gas = band.ka_nominal[iquad]; break;
+ case HTRDR_Ks: k_gas = band.ks_nominal[iquad]; break;
+ case HTRDR_Kext:
+ k_gas = band.ka_nominal[iquad] + band.ks_nominal[iquad];
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ } else {
+ /* Pos is inside the clouds. Compute the water vapor molar fraction at
+ * the current voxel */
+ const double x_h2o = cloud_water_vapor_molar_fraction(&sky->htcp_desc, ivox);
+ struct htgop_layer_sw_spectral_interval_tab tab;
+
+ /* Retrieve the tabulated data for the quadrature point */
+ HTGOP(layer_sw_spectral_interval_get_tab(&band, iquad, &tab));
+
+ /* Fetch the optical properties wrt x_h2o */
+ switch(prop) {
+ case HTRDR_Ka:
+ HTGOP(layer_sw_spectral_interval_tab_fetch_ka(&tab, x_h2o, &k_gas));
+ break;
+ case HTRDR_Ks:
+ HTGOP(layer_sw_spectral_interval_tab_fetch_ks(&tab, x_h2o, &k_gas));
+ break;
+ case HTRDR_Kext:
+ HTGOP(layer_sw_spectral_interval_tab_fetch_kext(&tab, x_h2o, &k_gas));
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
}
}
@@ -1327,21 +1704,6 @@ htrdr_sky_fetch_raw_property
return k;
}
-struct svx_tree*
-htrdr_sky_get_svx_tree
- (struct htrdr_sky* sky,
- const size_t ispectral_band,
- const size_t iquadrature_pt)
-{
- size_t i;
- ASSERT(sky);
- ASSERT(ispectral_band >= sky->sw_bands_range[0]);
- ASSERT(ispectral_band <= sky->sw_bands_range[1]);
- (void)iquadrature_pt;
- i = ispectral_band - sky->sw_bands_range[0];
- return sky->clouds[i].octree;
-}
-
res_T
htrdr_sky_dump_clouds_vtk
(const struct htrdr_sky* sky,
@@ -1349,8 +1711,8 @@ htrdr_sky_dump_clouds_vtk
const size_t iquad, /* Index of the quadrature point */
FILE* stream)
{
+ const struct cloud* cloud;
struct htgop_spectral_interval specint;
- struct svx_tree_desc desc;
struct octree_data data;
const double* leaf_data;
size_t nvertices;
@@ -1363,14 +1725,15 @@ htrdr_sky_dump_clouds_vtk
i = iband - sky->sw_bands_range[0];
octree_data_init(sky, iband, iquad, &data);
- SVX(tree_get_desc(sky->clouds[i].octree, &desc));
- ASSERT(desc.type == SVX_OCTREE);
+ cloud = &sky->clouds[i][iquad];
+
+ ASSERT(cloud->octree_desc.type == SVX_OCTREE);
/* Register leaf data */
- SVX(tree_for_each_leaf(sky->clouds[i].octree, register_leaf, &data));
+ SVX(tree_for_each_leaf(cloud->octree, register_leaf, &data));
nvertices = darray_double_size_get(&data.vertices) / 3/*#coords per vertex*/;
ncells = darray_size_t_size_get(&data.cells)/8/*#ids per cell*/;
- ASSERT(ncells == desc.nleaves);
+ ASSERT(ncells == cloud->octree_desc.nleaves);
/* Fetch the spectral interval descriptor */
HTGOP(get_sw_spectral_interval(sky->htgop, iband, &specint));
@@ -1433,7 +1796,11 @@ htrdr_sky_fetch_svx_property
const double pos[3])
{
struct svx_voxel voxel = SVX_VOXEL_NULL;
+ struct atmosphere* atmosphere;
+ struct cloud* cloud;
size_t i;
+ int in_clouds; /* Defines if `pos' lies in the clouds */
+ int in_atmosphere; /* Defines if `pos' lies in the atmosphere */
int comp_mask = components_mask;
ASSERT(sky && pos);
ASSERT(comp_mask & HTRDR_ALL_COMPONENTS);
@@ -1441,18 +1808,39 @@ htrdr_sky_fetch_svx_property
ASSERT(iband <= sky->sw_bands_range[1]);
i = iband - sky->sw_bands_range[0];
-
- /* Is the position outside the clouds? */
- if(pos[0] < sky->clouds[i].octree_desc.lower[0]
- || pos[1] < sky->clouds[i].octree_desc.lower[1]
- || pos[2] < sky->clouds[i].octree_desc.lower[2]
- || pos[0] > sky->clouds[i].octree_desc.upper[0]
- || pos[1] > sky->clouds[i].octree_desc.upper[1]
- || pos[2] > sky->clouds[i].octree_desc.upper[2]) {
- comp_mask &= ~HTRDR_PARTICLES; /* No particle */
+ cloud = &sky->clouds[i][iquad];
+ atmosphere = &sky->atmosphere[i][iquad];
+
+ /* Is the position inside the clouds? */
+ in_clouds =
+ pos[0] >= cloud->octree_desc.lower[0]
+ && pos[1] >= cloud->octree_desc.lower[1]
+ && pos[2] >= cloud->octree_desc.lower[2]
+ && pos[0] <= cloud->octree_desc.upper[0]
+ && pos[1] <= cloud->octree_desc.upper[1]
+ && pos[2] <= cloud->octree_desc.upper[2];
+
+ ASSERT(atmosphere->bitree_desc.frame[0] = SVX_AXIS_Z);
+ in_atmosphere =
+ pos[2] >= atmosphere->bitree_desc.lower[2]
+ && pos[2] <= atmosphere->bitree_desc.upper[2];
+
+ if(!in_clouds) { /* Not in clouds => No particle */
+ comp_mask &= ~HTRDR_PARTICLES;
}
+ if(!in_atmosphere) { /* Not in atmosphere => No gas */
+ comp_mask &= ~HTRDR_GAS;
+ }
+
+ if(!in_clouds && !in_atmosphere) /* In vacuum */
+ return 0;
- SVX(tree_at(sky->clouds[i].octree, pos, NULL, NULL, &voxel));
+ if(in_clouds) {
+ SVX(tree_at(cloud->octree, pos, NULL, NULL, &voxel));
+ } else {
+ ASSERT(in_atmosphere);
+ SVX(tree_at(atmosphere->bitree, pos, NULL, NULL, &voxel));
+ }
return htrdr_sky_fetch_svx_voxel_property
(sky, prop, op, comp_mask, iband, iquad, &voxel);
}
@@ -1469,15 +1857,23 @@ htrdr_sky_fetch_svx_voxel_property
{
double gas = 0;
double par = 0;
+ int comp_mask = components_mask;
ASSERT(sky && voxel);
ASSERT((unsigned)prop < HTRDR_PROPERTIES_COUNT__);
ASSERT((unsigned)op < HTRDR_SVX_OPS_COUNT__);
(void)sky, (void)ispectral_band, (void)iquad;
- if(components_mask & HTRDR_PARTICLES) {
+ /* Check if the voxel has infinite bounds/degenerated. In such case it is
+ * atmospheric voxel with only gas properties */
+ if(IS_INF(voxel->upper[0]) || voxel->lower[0] > voxel->upper[0]) {
+ ASSERT(IS_INF(voxel->upper[1]) || voxel->lower[1] > voxel->upper[1]);
+ comp_mask &= ~HTRDR_PARTICLES;
+ }
+
+ if(comp_mask & HTRDR_PARTICLES) {
par = voxel_get(voxel->data, HTRDR_PARTICLES__, prop, op);
}
- if(components_mask & HTRDR_GAS) {
+ if(comp_mask & HTRDR_GAS) {
gas = voxel_get(voxel->data, HTRDR_GAS__, prop, op);
}
/* Interval arithmetic to ensure that the returned Min/Max includes the
@@ -1500,6 +1896,18 @@ htrdr_sky_get_sw_spectral_band_id
return sky->sw_bands_range[0] + i;
}
+size_t
+htrdr_sky_get_sw_spectral_band_quadrature_length
+ (const struct htrdr_sky* sky, const size_t iband)
+{
+ struct htgop_spectral_interval band;
+ ASSERT(sky);
+ ASSERT(iband >= sky->sw_bands_range[0]);
+ ASSERT(iband <= sky->sw_bands_range[1]);
+ HTGOP(get_sw_spectral_interval(sky->htgop, iband, &band));
+ return band.quadrature_length;
+}
+
res_T
htrdr_sky_get_sw_spectral_band_bounds
(const struct htrdr_sky* sky,
@@ -1555,3 +1963,98 @@ htrdr_sky_sample_sw_spectral_data_CIE_1931_Z
ispectral_band, iquadrature_pt);
}
+res_T
+htrdr_sky_trace_ray
+ (struct htrdr_sky* sky,
+ const double org[3],
+ const double dir[3], /* Must be normalized */
+ const double range[2],
+ const svx_hit_challenge_T challenge, /* NULL <=> Traversed up to the leaves */
+ const svx_hit_filter_T filter, /* NULL <=> Stop RT at the 1st hit voxel */
+ void* context, /* Data sent to the filter functor */
+ const size_t ispectral_band,
+ const size_t iquadrature_pt,
+ struct svx_hit* hit)
+{
+ double cloud_range[2];
+ struct svx_tree* clouds;
+ struct svx_tree* atmosphere;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(sky);
+ ASSERT(ispectral_band >= sky->sw_bands_range[0]);
+ ASSERT(ispectral_band <= sky->sw_bands_range[1]);
+ (void)iquadrature_pt;
+
+ /* Fetch the clouds/atmosphere corresponding to the submitted spectral data */
+ i = ispectral_band - sky->sw_bands_range[0];
+ clouds = sky->clouds[i][iquadrature_pt].octree;
+ atmosphere = sky->atmosphere[i][iquadrature_pt].bitree;
+
+ cloud_range[0] = DBL_MAX;
+ cloud_range[1] =-DBL_MAX;
+
+ /* Compute the ray range, intersecting the clouds AABB */
+ ray_intersect_aabb(org, dir, range, sky->htcp_desc.lower,
+ sky->htcp_desc.upper, cloud_range);
+
+ /* Reset the hit */
+ *hit = SVX_HIT_NULL;
+
+ if(cloud_range[0] > cloud_range[1]) { /* The ray does not traverse the clouds */
+ res = svx_tree_trace_ray(atmosphere, org, dir, range, challenge, filter,
+ context, hit);
+ if(res != RES_OK) {
+ htrdr_log_err(sky->htrdr,
+ "%s: could not trace the ray in the atmosphere.\n",
+ FUNC_NAME);
+ goto error;
+ }
+ } else { /* The ray may traverse the clouds */
+ double range_adjusted[2];
+
+ if(cloud_range[0] > range[0]) { /* The ray begins in the atmosphere */
+ /* Trace a ray in the atmosphere from range[0] to cloud_range[0] */
+ range_adjusted[0] = range[0];
+ range_adjusted[1] = nextafter(cloud_range[0], -DBL_MAX);
+ res = svx_tree_trace_ray(atmosphere, org, dir, range_adjusted, challenge,
+ filter, context, hit);
+ if(res != RES_OK) {
+ htrdr_log_err(sky->htrdr,
+ "%s: could not to trace the part that begins in the atmosphere.\n",
+ FUNC_NAME);
+ goto error;
+ }
+ if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */
+ }
+
+ /* Pursue ray traversal into the clouds */
+ res = svx_tree_trace_ray(clouds, org, dir, cloud_range, challenge, filter,
+ context, hit);
+ if(res != RES_OK) {
+ htrdr_log_err(sky->htrdr,
+ "%s: could not trace the ray part that intersects the clouds.\n",
+ FUNC_NAME);
+ goto error;
+ }
+ if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */
+
+ /* Pursue ray traversal into the atmosphere */
+ range_adjusted[0] = nextafter(cloud_range[1], DBL_MAX);
+ range_adjusted[1] = range[1];
+ res = svx_tree_trace_ray(atmosphere, org, dir, range_adjusted, challenge,
+ filter, context, hit);
+ if(res != RES_OK) {
+ htrdr_log_err(sky->htrdr,
+ "%s: could not trace the ray part that in the atmosphere.\n", FUNC_NAME);
+ goto error;
+ }
+ if(!SVX_HIT_NONE(hit)) goto exit; /* Collision */
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
diff --git a/src/htrdr_sky.h b/src/htrdr_sky.h
@@ -17,6 +17,7 @@
#define HTRDR_SKY_H
#include <rsys/rsys.h>
+#include <star/svx.h>
/* Properties of the sky */
enum htrdr_sky_property {
@@ -26,7 +27,9 @@ enum htrdr_sky_property {
HTRDR_PROPERTIES_COUNT__
};
-enum htrdr_sky_component { /* FIXME */
+/* FIXME Maybe rename this constant to avoid the confusion with the
+ * htrdr_sky_component_flag enumerate */
+enum htrdr_sky_component {
HTRDR_GAS__,
HTRDR_PARTICLES__,
HTRDR_COMPONENTS_COUNT__
@@ -60,6 +63,7 @@ htrdr_sky_create
const char* htcp_filename,
const char* htgop_filename,
const char* htmie_filename,
+ const double optical_thickness, /* Threshold used during octree building */
struct htrdr_sky** sky);
extern LOCAL_SYM void
@@ -87,12 +91,6 @@ htrdr_sky_fetch_raw_property
const double k_min, /* For debug only */
const double k_max); /* For debug only */
-extern LOCAL_SYM struct svx_tree*
-htrdr_sky_get_svx_tree
- (struct htrdr_sky* sky,
- const size_t ispectral_band,
- const size_t iquadrature_pt);
-
extern LOCAL_SYM double
htrdr_sky_fetch_svx_property
(const struct htrdr_sky* sky,
@@ -129,6 +127,11 @@ htrdr_sky_get_sw_spectral_band_id
(const struct htrdr_sky* sky,
const size_t i); /* in [0, htrdr_sky_get_sw_spectral_bands_count[ */
+extern LOCAL_SYM size_t
+htrdr_sky_get_sw_spectral_band_quadrature_length
+ (const struct htrdr_sky* sky,
+ const size_t iband);
+
extern LOCAL_SYM res_T
htrdr_sky_get_sw_spectral_band_bounds
(const struct htrdr_sky* sky,
@@ -156,5 +159,18 @@ htrdr_sky_sample_sw_spectral_data_CIE_1931_Z
size_t* ispectral_band,
size_t* iquadrature_pt);
+extern LOCAL_SYM res_T
+htrdr_sky_trace_ray
+ (struct htrdr_sky* sky,
+ const double ray_origin[3],
+ const double ray_direction[3], /* Must be normalized */
+ const double ray_range[2],
+ const svx_hit_challenge_T challenge, /* NULL <=> Traversed up to the leaves */
+ const svx_hit_filter_T filter, /* NULL <=> Stop RT at the 1st hit voxel */
+ void* context, /* Data sent to the filter functor */
+ const size_t ispectral_band,
+ const size_t iquadrature_pt,
+ struct svx_hit* hit);
+
#endif /* HTRDR_SKY_H */
