commit 6e50b71038e29f34cf3da1c9fda874daea997af1
parent c6d6772a4922c273c7b7dd05d3ec3b7251dcdf3d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 10 Jul 2018 16:50:48 +0200
Implement the htrdr_sky_fetch_raw_property function
Diffstat:
4 files changed, 140 insertions(+), 31 deletions(-)
diff --git a/src/htrdr_draw_sw.c b/src/htrdr_draw_sw.c
@@ -15,22 +15,30 @@
struct scattering_context {
struct ssp_rng* rng;
+ const struct htrdr_sky* sky
+ double wavelength;
double Ts; /* Sampled optical thickness */
double traversal_dst; /* Distance traversed along the ray */
};
static const struct scattering_context SCATTERING_CONTEXT_NULL = {
- NULL, 0, 0, 0, 0
+ NULL, NULL, 0, 0, 0, 0, 0
};
struct transmissivity_context {
struct ssp_rng* rng;
+ const struct htrdr_sky* sky;
+ double wavelength;
+
double Ts; /* Sampled optical thickness */
double Tmin; /* Minimal optical thickness */
double traversal_dst; /* Distance traversed along the ray */
+
+
+ int prop_mask; /* Mask of htrdr_sky_property_flag */
};
static const struct transmissivity_context TRANSMISSION_CONTEXT_NULL = {
- NULL, 0, 0;
+ NULL, NULL, 0, 0, 0, 0, 0
};
/*******************************************************************************
@@ -53,7 +61,7 @@ scattering_hit_filter
(void)range;
vox_data = hit->voxel.data;
- ks_max = vox_data[HTRDR_Ks_MAX];
+ ks_max = vox_data[HTRDR_Ks_MAX/*TODO*/];
for(;;) {
dst = hit->distance[1] - ctx->traversal_dst;
@@ -67,7 +75,7 @@ scattering_hit_filter
break;
/* A real/null collision occurs before `dst' */
- } else {
+ else {
double pos[3];
double proba_s;
double ks;
@@ -82,7 +90,14 @@ scattering_hit_filter
pos[1] = org[1] + dst * dir[1];
pos[2] = org[2] + dst * dir[2];
- ks = fetch_ks(ctx->medium, x); /* TODO */
+ /* TODO use a per wavelength getter that will precompute the interpolated
+ * cross sections for a wavelength to improve the fetch efficiency */
+ ks = htrdr_sky_fetch_raw_property
+ (ctx->sky,
+ HTRDR_SKY_Ks,
+ HTRDR_SKY_PARTICLE | HTRDR_SKY_GAZ,
+ ctx->wavelength,
+ pos);
/* Handle the case that ks_max is not *really* the max */
proba_s = ks / ks_max;
@@ -149,7 +164,12 @@ transmissivty_hit_filter
x[1] = org[1] + dst * dir[1];
x[2] = org[2] + dst * dir[2];
- k = fetch_k(ctx->medium, x); /* TODO */
+ k = htrdr_sky_fetch_raw_property
+ (ctx->sky,
+ ctx->prop_mask,
+ HTRDR_SKY_PARTICLE | HTRDR_SKY_GAZ,
+ ctx->wavelength,
+ x);
/* Handle the case that k_max is not *really* the max */
proba = (k - k_min) / (k_max - k_min);
@@ -170,6 +190,8 @@ static double
transmissivity
(struct htrdr* htrdr,
struct ssp_rng* rng,
+ const int prop_mask, /* Combination of htrdr_sky_property_flag */
+ const double wavelength,
const double pos[3],
const double dir[3],
const double range[2],
@@ -184,7 +206,8 @@ transmissivity
float ray_dir[3];
float ray_range[2];
- ASSERT(htrdr && rng && pos && dir && range);
+ ASSERT(htrdr && rng && prop_mask && pos && dir && range);
+ ASSERT(prop_mask & HTRDR_SKY_PROP_Kext);
/* Find the first intersection with a surface */
f3_set_d3(ray_pos, pos);
@@ -196,7 +219,10 @@ transmissivity
if(!S3D_HIT_NONE(&s3d_hit)) return 0;
transmissivity_ctx->rng = rng;
+ transmissivity_ctx->sky = htrdr->sky;
+ transmissivity_ctx->wavelength = wavelength;
transmissivity_ctx->Ts = ssp_ran_exp(rng, 1); /* Sample an optical thickness */
+ transmissivity_ctx->prop_mask = prop_mask;
/* Compute the transmissivity */
SVX(octree_trace_ray(htrdr->clouds, pos, dir, range, NULL,
@@ -246,23 +272,28 @@ radiance_sw
ASSERT(htrdr && rng && pos && dir);
+ /* Fetch sun properties */
sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
+ L_sun = htrdr_sun_get_radiance(htrdr->sun, wavelength);
+
d3_set(pos, pos_in);
d3_set(dir, dir_in);
- /* Check that the sun is directly visible */
- f3_set_d3(ray_pos, pos);
- f3_set_d3(ray_dir, sun_dir);
- f2(ray_range, 0, FLT_MAX);
- S3D(scene_view_trace
- (htrdr->s3d_scn, ray_pos, ray_dir, ray_range, NULL, &s3d_hit));
-
- /* Add the direct contribution of the sun */
- if(!S3D_HIT_NONE(&s3d_hit)) {
- d3(range, 0, FLT_MAX);
- L_sun = htrdr_sun_get_luminance(htrdr->sun, sun_dir, wavelength);
- Tr = transmissivity(htrdr, rng, pos, sun_dir, range);
- w = L_sun * Tr;
+ /* Add the directly contribution of the sun */
+ if(htrdr_sun_is_dir_in_solar_cone(sun, dir)) {
+ f3_set_d3(ray_pos, pos);
+ f3_set_d3(ray_dir, dir);
+ f2(ray_range, 0, FLT_MAX);
+ S3D(scene_view_trace
+ (htrdr->s3d_scn, ray_pos, ray_dir, ray_range, NULL, &s3d_hit));
+
+ /* Add the direct contribution of the sun */
+ if(!S3D_HIT_NONE(&s3d_hit)) {
+ d3(range, 0, FLT_MAX);
+ Tr = transmissivity
+ (htrdr, rng, HTRDR_SKY_PROP_Kext, wavelength , pos, dir, range);
+ w = L_sun * Tr;
+ }
}
/* Radiative random walk */
@@ -279,9 +310,13 @@ radiance_sw
&s3d_hit_prev, &s3d_hit));
/* Sample an optical thicknes */
- scattering_ctx->rng = rng;
scattering_ctx->Ts = ssp_ran_exp(rng, 1);
+ /* Setup the remaining scattering context fields */
+ scattering_ctx->rng = rng;
+ scattering_ctx->sky = htrdr->sky;
+ scattering_ctx->wavelength = wavelength;
+
/* Define if a scattering event occurs */
d2(range, 0, s3d_hit.distance);
SVX(octree_trace_ray(htrdr->clouds, pos, dir, range, NULL,
@@ -304,13 +339,15 @@ radiance_sw
/* Define the absorption transmissivity from the current position to the
* next position */
d2(range, 0, scattering_ctx->traversal_dst);
- Tr_abs = transmissivity_absorption(htrdr, rng, pos, dir, range);
+ Tr_abs = transmissivity
+ (htrdr, rng, HTRDR_SKY_PROP_Ka, wavelength, pos, dir, range);
if(Tr_abs <= 0) break;
/* Define the transmissivity from the new position to the sun */
d2(range, 0, FLT_MAX);
htrdr_sun_sample_dir(htrdr->sun, rng, pos_next, sun_dir);
- Tr = transmissivity(htrdr, rng, pos_next, sun_dir, range);
+ Tr = transmissivity
+ (htrdr, rng, HTRDR_SKY_PROP_Kext, wavelength, pos_next, sun_dir, range);
/* Scattering at a surface */
if(SVX_HIT_NONE(&svx_hit)) {
diff --git a/src/htrdr_sky.c b/src/htrdr_sky.c
@@ -71,6 +71,8 @@ struct octree_data {
struct htrdr_sky {
struct svx_tree* clouds;
+ struct svx_tree_desc cloud_desc;
+
struct htcp* htcp;
struct htmie* htmie;
@@ -354,6 +356,9 @@ setup_clouds(struct htrdr_sky* sky)
goto error;
}
+ /* Fetch the octree descriptor for future use */
+ SVX(tree_get_desc(sky->clouds, &sky->cloud_desc));
+
exit:
return res;
error:
@@ -453,6 +458,70 @@ htrdr_sky_ref_put(struct htrdr_sky* sky)
ref_put(&sky->ref, release_sky);
}
+double
+htrdr_sky_fetch_raw_property
+ (const struct htrdr_sky* sky,
+ const int prop_mask, /* Combination of htrdr_sky_property_flag */
+ const int comp_mask_in, /* Combination of htrdr_sky_component_flag */
+ const double wavelength,
+ const double pos[3])
+{
+ size_t ivox[3];
+ int comp_mask = comp_mask_in;
+ double k_particle = 0;
+ double k_gaz = 0;
+ ASSERT(sky && pos && wavelength > 0);
+ ASSERT(prop_mask & HTRDR_SKY_PROP_Kext);
+ ASSERT(components_mask & HTRDR_SKY_COMP_ALL);
+
+ /* Is the position outside the clouds? */
+ if(pos[0] < sky->cloud_desc.lower[0]
+ || pos[1] < sky->cloud_desc.lower[1]
+ || pos[2] < sky->cloud_desc.lower[2]
+ || pos[0] > sky->cloud_desc.upper[0]
+ || pos[1] > sky->cloud_desc.upper[1]
+ || pos[2] > sky->cloud_desc.upper[2]) {
+ comp_mask &= ~HTRDR_SKY_COMP_PARTICLE; /* No particle */
+ }
+
+ /* Compute the index of the voxel to fetch */
+ ivox[0] = (size_t)((pos[0] - sky->cloud_desc.lower[0])/sky->htcp_desc.vxsz_x);
+ ivox[1] = (size_t)((pos[1] - sky->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] - sky->cloud_desc.lower[2])/sky->htcp_desc.vxsz_z[0]);
+ } else {
+ /* Irregular voxel size along the Z dimension. Compute the index 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 n = darray_split_size_get(&sky->svx2htcp_z);
+ const double sz = sky->cloud_desc.upper[2] - sky->cloud_desc.lower[2];
+ const double vxsz_lut = sz / (double)n;
+ const size_t ilut = (size_t)((pos[2] - sky->cloud_desc.lower[2])/vxsz_lut);
+ ivox[2] = splits[ilut].index + (pos[2] >= splits[ilut].height);
+ }
+
+ if(comp_mask & HTRDR_SKY_COMP_PARTICLE) {
+ 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 */
+
+ ql = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
+ if(prop_mask & HTRDR_SKY_PROP_Ks) {
+ cs = htmie_fetch_xsection_scattering
+ (sky->htmie, wavelength, HTMIE_FILTER_LINEAR);
+ }
+ if(prop_mask & HTRDR_SKY_PROP_Ka) {
+ ca = htmie_fetch_xsection_absorption
+ (sky->htmie, wavelength, HTMIE_FILTER_LINEAR);
+ }
+ k_particle = ql*(ca + cs);
+ }
+ if(comp_mask & HTRDR_SKY_COMP_GAZ) { /* TODO not implemented yet */ }
+ return k_particle + k_gaz;
+}
+
struct svx_tree*
htrdr_sky_get_svx_tree(struct htrdr_sky* sky)
{
@@ -567,7 +636,7 @@ htrdr_sky_fetch_svx_voxel_property
ASSERT((unsigned)prop < HTRDR_SKY_SVX_PROPS_COUNT__);
(void)sky, (void)wavelength;
- if(components_mask != (HTRDR_SKY_GAZ|HTRDR_SKY_PARTICLE)) {
+ if(components_mask != HTRDR_SKY_COMP_ALL) {
FATAL("Unsupported sky component\n");
}
pdbl = voxel->data;
diff --git a/src/htrdr_sky.h b/src/htrdr_sky.h
@@ -19,9 +19,10 @@
#include <rsys/rsys.h>
/* Raw sky properties */
-enum htrdr_sky_property {
- HTRDR_SKY_Ks, /* Scattering coefficient */
- HTRDR_SKY_Ka /* Absorption coefficient */
+enum htrdr_sky_property_flag {
+ HTRDR_SKY_PROP_Ks = BIT(0), /* Scattering coefficient */
+ HTRDR_SKY_PROP_Ka = BIT(1), /* Absorption coefficient */
+ HTRDR_SKY_PROP_Kext = HTRDR_SKY_PROP_Ks | HTRDR_SKY_PROP_Ka
};
/* Property of the sky computed by region and managed by Star-VoXel */
@@ -29,12 +30,14 @@ enum htrdr_sky_svx_property {
HTRDR_SKY_SVX_Kext_MIN,
HTRDR_SKY_SVX_Kext_MAX,
HTRDR_SKY_SVX_PROPS_COUNT__
+
};
/* Component of the sky for which the properties are queried */
enum htrdr_sky_component_flag {
- HTRDR_SKY_GAZ = BIT(0),
- HTRDR_SKY_PARTICLE = BIT(1)
+ HTRDR_SKY_COMP_GAZ = BIT(0),
+ HTRDR_SKY_COMP_PARTICLE = BIT(1),
+ HTRDR_SKY_COMP_ALL = HTRDR_SKY_COMP_GAZ | HTRDR_SKY_COMP_PARTICLE
};
/* Forward declaration */
@@ -61,7 +64,7 @@ htrdr_sky_ref_put
extern LOCAL_SYM double
htrdr_sky_fetch_raw_property
(const struct htrdr_sky* sky,
- const enum htrdr_sky_property prop,
+ const int prop_mask, /* Combination of htrdr_sky_property_flag */
const int components_mask, /* Combination of htrdr_sky_component_flag */
const double wavelength,
const double pos[3]);
diff --git a/src/htrdr_transmission.c b/src/htrdr_transmission.c
@@ -55,7 +55,7 @@ discard_hit
const double range[2],
void* context)
{
- const int comp = HTRDR_SKY_GAZ | HTRDR_SKY_PARTICLE;
+ const int comp = HTRDR_SKY_COMP_GAZ | HTRDR_SKY_COMP_PARTICLE;
struct transmit_context* ctx = context;
double dst;
double k_ext_min;
