commit 558fdd4538ca59e5a8353a893b0a3a2c11a3a7e8
parent 1a8e8e8f205b573a0a7e4251a9c12f9a36a1e0f2
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 6 Aug 2018 16:48:03 +0200
Change how the spectral dimension is handled
Use the spectral data provided by the htgop. At the beginning of a
realisation, sample a spectral band and then a quadrature point into the
sampled spectral band.
Diffstat:
13 files changed, 534 insertions(+), 180 deletions(-)
diff --git a/README.md b/README.md
@@ -10,6 +10,7 @@ This program is compatible GNU/Linux 64-bits. It relies on the
[RCMake](https://gitlab.com/vaplv/rcmake/) packages to build. It also depends
on the
[HTCP](https://gitlab.com/meso-star/htcp/),
+[HTGOP](https://gitlab.com/meso-star/htgop/),
[HTMIE](https://gitlab.com/meso-star/htmie/),
[RSys](https://gitlab.com/vaplv/rsys/),
[Star-3D](https://gitlab.com/meso-star/star-3d/),
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -31,6 +31,7 @@ find_package(StarSF 0.6 REQUIRED)
find_package(StarSP 0.8 REQUIRED)
find_package(StarVX REQUIRED)
find_package(HTCP REQUIRED)
+find_package(HTGOP REQUIRED)
find_package(HTMIE REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
@@ -46,6 +47,7 @@ include_directories(
${StarSP_INCLUDE_DIR}
${StarVX_INCLUDE_DIR}
${HTCP_INCLUDE_DIR}
+ ${HTGOP_INCLUDE_DIR}
${HTMIE_INCLUDE_DIR})
################################################################################
@@ -85,7 +87,7 @@ rcmake_prepend_path(HTRDR_FILES_INC ${HTRDR_SOURCE_DIR})
rcmake_prepend_path(HTRDR_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_executable(htrdr ${HTRDR_FILES_SRC} ${HTRDR_FILES_INC})
-target_link_libraries(htrdr HTCP HTMIE RSys Star3D Star3DAW StarSF StarSP StarVX)
+target_link_libraries(htrdr HTCP HTGOP HTMIE RSys Star3D Star3DAW StarSF StarSP StarVX)
if(CMAKE_COMPILER_IS_GNUCC)
target_link_libraries(htrdr m)
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -240,6 +240,7 @@ htrdr_init
{
double proj_ratio;
const char* output_name = NULL;
+ size_t ithread;
res_T res = RES_OK;
ASSERT(args && htrdr);
@@ -312,31 +313,27 @@ htrdr_init
htrdr_sun_set_direction(htrdr->sun, args->main_dir);
res = htrdr_sky_create(htrdr, htrdr->sun, args->filename_les,
- args->filename_mie, &htrdr->sky);
+ args->filename_gas, args->filename_mie, &htrdr->sky);
if(res != RES_OK) goto error;
- res = ssf_bsdf_create
- (htrdr->allocator, &ssf_lambertian_reflection, &htrdr->bsdf);
- if(res != RES_OK) {
- htrdr_log_err(htrdr, "could not create the BSDF of the ground.\n");
- goto error;
- }
- SSF(lambertian_reflection_setup(htrdr->bsdf, 0.02));
-/* SSF(lambertian_reflection_setup(htrdr->bsdf, 1));*/
-
- res = ssf_phase_create(htrdr->allocator, &ssf_phase_hg, &htrdr->phase_hg);
- if(res != RES_OK) {
+ htrdr->lifo_allocators = MEM_CALLOC
+ (htrdr->allocator, htrdr->nthreads, sizeof(*htrdr->lifo_allocators));
+ if(!htrdr->lifo_allocators) {
htrdr_log_err(htrdr,
- "could not create the Henyey & Greenstein phase function.\n");
+ "could not allocate the list of per thread LIFO allocator.\n");
+ res = RES_MEM_ERR;
goto error;
}
- SSF(phase_hg_setup(htrdr->phase_hg, 0));
- res = ssf_phase_create
- (htrdr->allocator, &ssf_phase_rayleigh, &htrdr->phase_rayleigh);
- if(res != RES_OK) {
- htrdr_log_err(htrdr,"could not create the Rayleigh phase function.\n");
- goto error;
+ FOR_EACH(ithread, 0, htrdr->nthreads) {
+ res = mem_init_lifo_allocator
+ (&htrdr->lifo_allocators[ithread], htrdr->allocator, 4096);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr,
+ "could not initialise the LIFO allocator of the thread %lu.\n",
+ (unsigned long)ithread);
+ goto error;
+ }
}
res = setup_geometry(htrdr, args->filename_obj);
@@ -353,9 +350,6 @@ void
htrdr_release(struct htrdr* htrdr)
{
ASSERT(htrdr);
- if(htrdr->bsdf) SSF(bsdf_ref_put(htrdr->bsdf));
- if(htrdr->phase_hg) SSF(phase_ref_put(htrdr->phase_hg));
- if(htrdr->phase_rayleigh) SSF(phase_ref_put(htrdr->phase_rayleigh));
if(htrdr->s3d_scn_view) S3D(scene_view_ref_put(htrdr->s3d_scn_view));
if(htrdr->s3d) S3D(device_ref_put(htrdr->s3d));
if(htrdr->svx) SVX(device_ref_put(htrdr->svx));
@@ -363,6 +357,13 @@ htrdr_release(struct htrdr* htrdr)
if(htrdr->sun) htrdr_sun_ref_put(htrdr->sun);
if(htrdr->cam) htrdr_camera_ref_put(htrdr->cam);
if(htrdr->buf) htrdr_buffer_ref_put(htrdr->buf);
+ if(htrdr->lifo_allocators) {
+ size_t i;
+ FOR_EACH(i, 0, htrdr->nthreads) {
+ mem_shutdown_lifo_allocator(&htrdr->lifo_allocators[i]);
+ }
+ MEM_RM(htrdr->allocator, htrdr->lifo_allocators);
+ }
str_release(&htrdr->output_name);
logger_release(&htrdr->logger);
}
@@ -372,7 +373,8 @@ htrdr_run(struct htrdr* htrdr)
{
res_T res = RES_OK;
if(htrdr->dump_vtk) {
- res = htrdr_sky_dump_clouds_vtk(htrdr->sky, htrdr->output);
+ const size_t iband = htrdr_sky_get_sw_spectral_band_id(htrdr->sky, 0);
+ res = htrdr_sky_dump_clouds_vtk(htrdr->sky, iband, 0, htrdr->output);
if(res != RES_OK) goto error;
} else {
struct time t0, t1;
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -36,16 +36,11 @@ struct htrdr {
struct svx_device* svx;
struct s3d_device* s3d;
- struct s3d_scene_view* s3d_scn_view;
+ struct s3d_scene_view* s3d_scn_view;
struct htrdr_sky* sky;
struct htrdr_sun* sun;
- /* Scattering functions */
- struct ssf_bsdf* bsdf; /* BSDF of the 3D geometry */
- struct ssf_phase* phase_hg; /* Henyey & Greenstein phase function */
- struct ssf_phase* phase_rayleigh; /* Rayleigh phase function */
-
struct htrdr_camera* cam;
struct htrdr_buffer* buf;
size_t spp; /* #samples per pixel */
@@ -59,6 +54,7 @@ struct htrdr {
struct logger logger;
struct mem_allocator* allocator;
+ struct mem_allocator* lifo_allocators; /* Per thread lifo allocator */
};
extern LOCAL_SYM res_T
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -32,6 +32,8 @@ print_help(const char* cmd)
ASSERT(cmd);
printf("Usage: %s -i INPUT [OPIONS]\n\n", cmd);
printf(
+" -a FILENAME path of gas optical properties file.\n");
+ printf(
" -c FILENAME path of the HTCP cloud properties file.\n");
printf(
" -C <camera> define the rendering point of view.\n");
@@ -292,8 +294,9 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
*args = HTRDR_ARGS_DEFAULT;
- while((opt = getopt(argc, argv, "C:c:D:dfg:hi:m:o:t:v")) != -1) {
+ while((opt = getopt(argc, argv, "a:C:c:D:dfg:hi:m:o:t:v")) != -1) {
switch(opt) {
+ case 'a': args->filename_gas = optarg; break;
case 'C':
res = parse_multiple_parameters
(args, optarg, parse_camera_parameter);
@@ -329,6 +332,12 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
goto error;
}
}
+ if(!args->filename_gas) {
+ fprintf(stderr,
+ "Missing the path of the gas optical properties file -- option '-a'\n");
+ res = RES_BAD_ARG;
+ goto error;
+ }
if(!args->filename_les) {
fprintf(stderr,
"Missing the path of the cloud properties file -- option '-c'\n");
diff --git a/src/htrdr_args.h b/src/htrdr_args.h
@@ -19,9 +19,10 @@
#include <rsys/rsys.h>
struct htrdr_args {
- const char* filename_les; /* Path toward the HTCP file */
- const char* filename_mie; /* Path toward the Mie properties */
- const char* filename_obj; /* Path toward the 3D geometry */
+ const char* filename_gas; /* Path of the gas file */
+ const char* filename_les; /* Path of the HTCP file */
+ const char* filename_mie; /* Path of the Mie properties */
+ const char* filename_obj; /* Path of the 3D geometry */
const char* output;
struct {
@@ -53,6 +54,7 @@ struct htrdr_args {
};
#define HTRDR_ARGS_DEFAULT__ { \
+ NULL, /* Gas filename */ \
NULL, /* LES filename */ \
NULL, /* Mie filename */ \
NULL, /* Obj filename */ \
diff --git a/src/htrdr_compute_radiance_sw.c b/src/htrdr_compute_radiance_sw.c
@@ -31,19 +31,21 @@
struct scattering_context {
struct ssp_rng* rng;
const struct htrdr_sky* sky;
- double wavelength;
+ size_t iband; /* Index of the spectrald */
+ size_t iquad; /* Index of the quadrature point into the band */
double Ts; /* Sampled optical thickness */
double traversal_dst; /* Distance traversed along the ray */
};
static const struct scattering_context SCATTERING_CONTEXT_NULL = {
- NULL, NULL, 0, 0, 0
+ NULL, NULL, 0, 0, 0, 0
};
struct transmissivity_context {
struct ssp_rng* rng;
const struct htrdr_sky* sky;
- double wavelength;
+ size_t iband; /* Index of the spectrald */
+ size_t iquad; /* Index of the quadrature point into the band */
double Ts; /* Sampled optical thickness */
double Tmin; /* Minimal optical thickness */
@@ -52,7 +54,7 @@ struct transmissivity_context {
enum htrdr_sky_property prop;
};
static const struct transmissivity_context TRANSMISSION_CONTEXT_NULL = {
- NULL, NULL, 0, 0, 0, 0, 0
+ NULL, NULL, 0, 0, 0, 0, 0, 0
};
/*******************************************************************************
@@ -74,7 +76,7 @@ scattering_hit_filter
(void)range;
ks_max = htrdr_sky_fetch_svx_voxel_property(ctx->sky, HTRDR_Ks,
- HTRDR_SVX_MAX, HTRDR_ALL_COMPONENTS, ctx->wavelength, &hit->voxel);
+ HTRDR_SVX_MAX, HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, &hit->voxel);
/* Compute the distance to traverse into the voxel */
vox_dst = hit->distance[1] - hit->distance[0];
@@ -108,10 +110,8 @@ scattering_hit_filter
pos[1] = org[1] + ctx->traversal_dst * dir[1];
pos[2] = org[2] + ctx->traversal_dst * dir[2];
- /* 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_Ks, HTRDR_ALL_COMPONENTS, ctx->wavelength, pos);
+ ks = htrdr_sky_fetch_raw_property(ctx->sky, HTRDR_Ks,
+ HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, pos);
/* Handle the case that ks_max is not *really* the max */
proba = ks / ks_max;
@@ -146,9 +146,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->wavelength, &hit->voxel);
+ HTRDR_SVX_MIN, HTRDR_ALL_COMPONENTS, 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->wavelength, &hit->voxel);
+ HTRDR_SVX_MAX, HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, &hit->voxel);
ASSERT(k_min <= k_max);
/* Compute the distance to traverse into the voxel */
@@ -184,10 +184,8 @@ transmissivity_hit_filter
x[1] = org[1] + ctx->traversal_dst * dir[1];
x[2] = org[2] + ctx->traversal_dst * dir[2];
- /* TODO use a per wavelength getter that will precompute the interpolated
- * cross sections for a wavelength to improve the fetch efficiency */
- k = htrdr_sky_fetch_raw_property
- (ctx->sky, ctx->prop, HTRDR_ALL_COMPONENTS, ctx->wavelength, x);
+ k = htrdr_sky_fetch_raw_property(ctx->sky, ctx->prop,
+ HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, x);
ASSERT(k >= k_min && k <= k_max);
/* Handle the case that k_max is not *really* the max */
@@ -211,7 +209,8 @@ transmissivity
(struct htrdr* htrdr,
struct ssp_rng* rng,
const enum htrdr_sky_property prop,
- const double wavelength,
+ const size_t iband,
+ const size_t iquad,
const double pos[3],
const double dir[3],
const double range[2],
@@ -240,12 +239,13 @@ transmissivity
transmissivity_ctx.rng = rng;
transmissivity_ctx.sky = htrdr->sky;
- transmissivity_ctx.wavelength = wavelength;
+ transmissivity_ctx.iband = iband;
+ transmissivity_ctx.iquad = iquad;
transmissivity_ctx.Ts = ssp_ran_exp(rng, 1); /* Sample an optical thickness */
transmissivity_ctx.prop = prop;
/* Compute the transmissivity */
- svx_tree = htrdr_sky_get_svx_tree(htrdr->sky);
+ svx_tree = htrdr_sky_get_svx_tree(htrdr->sky, iband, iquad);
SVX(octree_trace_ray(svx_tree, pos, dir, range, NULL,
transmissivity_hit_filter, &transmissivity_ctx, &svx_hit));
@@ -256,29 +256,35 @@ transmissivity
}
}
-
/*******************************************************************************
* Local functions
******************************************************************************/
double
htrdr_compute_radiance_sw
(struct htrdr* htrdr,
+ const size_t ithread,
struct ssp_rng* rng,
const double pos_in[3],
const double dir_in[3],
- const double wavelength)
+ const size_t iband,
+ const size_t iquad)
{
struct s3d_hit s3d_hit = S3D_HIT_NULL;
struct svx_hit svx_hit = SVX_HIT_NULL;
struct s3d_hit s3d_hit_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;
double pos[3];
double dir[3];
double range[2];
double pos_next[3];
double dir_next[3];
+ double band_bounds[2]; /* In nanometers */
+ double wlen;
double R;
double r; /* Random number */
double wo[3]; /* -dir */
@@ -296,19 +302,33 @@ htrdr_compute_radiance_sw
float ray_dir[3];
float ray_range[2];
- ASSERT(wavelength >= SW_WAVELENGTH_MIN && wavelength <= SW_WAVELENGTH_MAX);
- ASSERT(htrdr && rng && pos_in && dir_in);
+ ASSERT(htrdr && rng && pos_in && dir_in && ithread < htrdr->nthreads);
- /* Setup the phase function for this wavelength */
- g = htrdr_sky_fetch_particle_phase_function_asymmetry_parameter
- (htrdr->sky, wavelength);
- SSF(phase_hg_setup(htrdr->phase_hg, g));
+ CHK(RES_OK == ssf_bsdf_create
+ (&htrdr->lifo_allocators[ithread], &ssf_lambertian_reflection, &bsdf));
+ CHK(RES_OK == ssf_phase_create
+ (&htrdr->lifo_allocators[ithread], &ssf_phase_hg, &phase_hg));
+ CHK(RES_OK == ssf_phase_create
+ (&htrdr->lifo_allocators[ithread], &ssf_phase_rayleigh, &phase_rayleigh));
- /* Fetch sun properties */
+ SSF(lambertian_reflection_setup(bsdf, 0.02));
+
+ /* Setup the phase function for this spectral band & quadrature point */
+ g = htrdr_sky_fetch_particle_phase_function_asymmetry_parameter
+ (htrdr->sky, iband, iquad);
+ SSF(phase_hg_setup(phase_hg, g));
+
+ /* Fetch sun properties. Arbitrarily use the wavelength at the center of the
+ * band to retrieve the sun radiance of the current band. Note that the sun
+ * spectral data are defined by bands that, actually are the same of the SW
+ * spectral bands defined in the default "ecrad_opt_prot.txt" file provided
+ * by the HTGOP project. */
+ htrdr_sky_get_sw_spectral_band_bounds(htrdr->sky, iband, band_bounds);
+ wlen = (band_bounds[0] + band_bounds[1]) * 0.5;
sun_solid_angle = htrdr_sun_get_solid_angle(htrdr->sun);
- L_sun = htrdr_sun_get_radiance(htrdr->sun, wavelength);
+ L_sun = htrdr_sun_get_radiance(htrdr->sun, wlen);
- svx_tree = htrdr_sky_get_svx_tree(htrdr->sky);
+ svx_tree = htrdr_sky_get_svx_tree(htrdr->sky, iband, iquad);
d3_set(pos, pos_in);
d3_set(dir, dir_in);
@@ -325,7 +345,7 @@ htrdr_compute_radiance_sw
if(S3D_HIT_NONE(&s3d_hit)) {
d2(range, 0, FLT_MAX);
Tr = transmissivity
- (htrdr, rng, HTRDR_Kext, wavelength , pos, dir, range, NULL);
+ (htrdr, rng, HTRDR_Kext, iband, iquad , pos, dir, range, NULL);
w = L_sun * Tr;
}
}
@@ -349,7 +369,8 @@ htrdr_compute_radiance_sw
/* Setup the remaining scattering context fields */
scattering_ctx.rng = rng;
scattering_ctx.sky = htrdr->sky;
- scattering_ctx.wavelength = wavelength;
+ scattering_ctx.iband = iband;
+ scattering_ctx.iquad = iquad;
/* Define if a scattering event occurs */
d2(range, 0, s3d_hit.distance);
@@ -377,14 +398,14 @@ htrdr_compute_radiance_sw
* next position */
d2(range, 0, scattering_ctx.traversal_dst);
Tr_abs = transmissivity
- (htrdr, rng, HTRDR_Ka, wavelength, pos, dir, range, &s3d_hit_prev);
+ (htrdr, rng, HTRDR_Ka, iband, iquad, pos, dir, range, &s3d_hit_prev);
if(Tr_abs <= 0) break;
/* Define the transmissivity from the new position to the sun */
d2(range, 0, FLT_MAX);
htrdr_sun_sample_direction(htrdr->sun, rng, sun_dir);
- Tr = transmissivity
- (htrdr, rng, HTRDR_Kext, wavelength, pos_next, sun_dir, range, &s3d_hit);
+ Tr = transmissivity(htrdr, rng, HTRDR_Kext, iband, iquad, pos_next,
+ sun_dir, range, &s3d_hit);
/* Scattering at a surface */
if(SVX_HIT_NONE(&svx_hit)) {
@@ -395,10 +416,10 @@ htrdr_compute_radiance_sw
d3_normalize(N, d3_set_f3(N, s3d_hit.normal));
if(d3_dot(N, wo) < 0) d3_minus(N, N);
reflectivity = ssf_bsdf_sample
- (htrdr->bsdf, rng, wo, N, dir_next, &type, &pdf);
+ (bsdf, rng, wo, N, dir_next, &type, &pdf);
if(ssp_rng_canonical(rng) > reflectivity) break; /* Russian roulette */
- R = ssf_bsdf_eval(htrdr->bsdf, wo, N, sun_dir);
+ R = ssf_bsdf_eval(bsdf, wo, N, sun_dir);
/* Scattering in the medium */
} else {
@@ -408,16 +429,16 @@ htrdr_compute_radiance_sw
double ks; /* Overall scattering coefficient */
ks_gas = htrdr_sky_fetch_raw_property
- (htrdr->sky, HTRDR_Ks, HTRDR_GAS, wavelength, pos_next);
+ (htrdr->sky, HTRDR_Ks, HTRDR_GAS, iband, iquad, pos_next);
ks_particle = htrdr_sky_fetch_raw_property
- (htrdr->sky, HTRDR_Ks, HTRDR_PARTICLES, wavelength, pos_next);
+ (htrdr->sky, HTRDR_Ks, HTRDR_PARTICLES, iband, iquad, pos_next);
ks = ks_particle + ks_gas;
r = ssp_rng_canonical(rng);
if(r < ks_gas / ks) { /* Gas scattering */
- phase = htrdr->phase_rayleigh;
+ phase = phase_rayleigh;
} else { /* Cloud scattering */
- phase = htrdr->phase_hg;
+ phase = phase_hg;
}
ssf_phase_sample(phase, rng, wo, dir_next, NULL);
@@ -434,6 +455,9 @@ htrdr_compute_radiance_sw
s3d_hit_prev = s3d_hit;
}
+ SSF(bsdf_ref_put(bsdf));
+ SSF(phase_ref_put(phase_hg));
+ SSF(phase_ref_put(phase_rayleigh));
return w;
}
diff --git a/src/htrdr_draw_radiance_sw.c b/src/htrdr_draw_radiance_sw.c
@@ -17,6 +17,7 @@
#include "htrdr_c.h"
#include "htrdr_buffer.h"
#include "htrdr_camera.h"
+#include "htrdr_sky.h"
#include "htrdr_solve.h"
#include <star/ssp.h>
@@ -44,7 +45,8 @@ morton2D_decode(const uint32_t u32)
static res_T
draw_tile
(struct htrdr* htrdr,
- int64_t tile_mcode, /* For debug only */
+ const size_t ithread,
+ const int64_t tile_mcode, /* For debug only */
const size_t tile_org[2], /* Origin of the tile in pixel space */
const size_t tile_sz[2], /* Definition of the tile */
const double pix_sz[2], /* Size of a pixel in the normalized image plane */
@@ -85,6 +87,8 @@ draw_tile
double ray_org[3];
double ray_dir[3];
double weight;
+ size_t iband;
+ size_t iquad;
/* Sample a position into the pixel, in the normalized image plane */
pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
@@ -94,9 +98,13 @@ draw_tile
* pixel sample */
htrdr_camera_ray(cam, pix_samp, ray_org, ray_dir);
+ /* Sample a spectral band and a quadrature point */
+ htrdr_sky_sample_sw_spectral_data_CIE_1931_X
+ (htrdr->sky, rng, &iband, &iquad);
+
/* Compute the radiance that reach the pixel through the ray */
weight = htrdr_compute_radiance_sw
- (htrdr, rng, ray_org, ray_dir, SW_WAVELENGTH_MIN/*FIXME wavelength*/);
+ (htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
ASSERT(weight >= 0);
/* Update the pixel accumulator */
@@ -207,8 +215,8 @@ htrdr_draw_radiance_sw
tile_sz[0] = MMIN(TILE_SIZE, layout.width - tile_org[0]);
tile_sz[1] = MMIN(TILE_SIZE, layout.height - tile_org[1]);
- res_local = draw_tile
- (htrdr, mcode, tile_org, tile_sz, pix_sz, cam, spp, rng, buf);
+ res_local = draw_tile(htrdr, (size_t)ithread, mcode, tile_org, tile_sz,
+ pix_sz, cam, spp, rng, buf);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
continue;
diff --git a/src/htrdr_sky.c b/src/htrdr_sky.c
@@ -20,8 +20,10 @@
#include "htrdr_sky.h"
#include "htrdr_sun.h"
+#include <star/ssp.h>
#include <star/svx.h>
#include <high_tune/htcp.h>
+#include <high_tune/htgop.h>
#include <high_tune/htmie.h>
#include <rsys/dynamic_array_double.h>
@@ -47,6 +49,7 @@ struct build_octree_context {
const struct htrdr_sky* sky;
double dst_max; /* Max traversal distance */
double tau_threshold; /* Threshold criteria for the merge process */
+ size_t iband; /* Index of the band that overlaps the CIE XYZ color space */
};
struct vertex {
@@ -70,34 +73,53 @@ vertex_eq(const struct vertex* v0, const struct vertex* v1)
#define HTABLE_KEY_FUNCTOR_EQ vertex_eq
#include <rsys/hash_table.h>
+/* Temporary data structure used to dump the octree data into a VTK file */
struct octree_data {
struct htable_vertex vertex2id; /* Map a coordinate to its vertex id */
struct darray_double vertices; /* Array of unique vertices */
struct darray_double data;
struct darray_size_t cells;
+ size_t iband; /* Index of the band that overlaps the CIE XYZ color space */
+ size_t iquad; /* Index of the quadrature point into the band */
const struct htrdr_sky* sky;
};
+/* Properties of a short wave spectral band */
+struct sw_band_prop {
+ /* Average cross section in the band */
+ double Ca_avg; /* Absorption cross section */
+ double Cs_avg; /* Scattering cross section */
+
+ /* Average asymmetry parameter the band */
+ double g_avg;
+};
+
+/* Encompass the hierarchical data structure of the cloud data and its
+ * associated descriptor */
+struct cloud {
+ struct svx_tree* octree;
+ struct svx_tree_desc octree_desc;
+};
+
struct htrdr_sky {
- struct svx_tree* clouds;
- struct svx_tree_desc cloud_desc;
+ struct cloud* clouds; /* Per sw_band cloud data structure */
- struct htrdr_sun* sun;
+ struct htrdr_sun* sun; /* Sun attached to the sky */
- struct htcp* htcp;
- struct htmie* htmie;
+ /* Loaders of... */
+ struct htcp* htcp; /* ... Cloud properties */
+ struct htgop* htgop; /* ... Gas optical properties */
+ struct htmie* htmie; /* ... Mie's data */
struct htcp_desc htcp_desc;
/* Map the index in Z from the regular SVX to the irregular HTCP data */
struct darray_split svx2htcp_z;
- /* Average cross section in Short Wave, i.e. in [380, 780] nanometers */
- double Ca_avg_sw;
- double Cs_avg_sw;
-
- /* Average asymmetry parameter in Short Wave, in [380, 780] nanometers */
- double g_avg_sw;
+ /* Ids and optical properties of the short wave spectral bands loaded by
+ * HTGOP and that overlap the CIE XYZ color space */
+ size_t sw_bands_range[2];
+ struct sw_band_prop* sw_bands;
ref_T ref;
struct htrdr* htrdr;
@@ -121,14 +143,23 @@ cloud_dry_air_density
}
static INLINE void
-octree_data_init(const struct htrdr_sky* sky, struct octree_data* data)
+octree_data_init
+ (const struct htrdr_sky* sky,
+ const size_t ispectral_band,
+ const size_t iquad,
+ struct octree_data* data)
{
ASSERT(data);
+ ASSERT(ispectral_band >= sky->sw_bands_range[0]);
+ ASSERT(ispectral_band <= sky->sw_bands_range[1]);
+ (void)iquad;
htable_vertex_init(sky->htrdr->allocator, &data->vertex2id);
darray_double_init(sky->htrdr->allocator, &data->vertices);
darray_double_init(sky->htrdr->allocator, &data->data);
darray_size_t_init(sky->htrdr->allocator, &data->cells);
data->sky = sky;
+ data->iband = ispectral_band - sky->sw_bands_range[0];
+ data->iquad = iquad;
}
static INLINE void
@@ -182,10 +213,10 @@ register_leaf
}
/* Register the leaf data */
- kext_max = htrdr_sky_fetch_svx_voxel_property
- (ctx->sky, HTRDR_Kext, HTRDR_SVX_MAX, HTRDR_ALL_COMPONENTS, -1/*FIXME*/, leaf);
- kext_min = htrdr_sky_fetch_svx_voxel_property
- (ctx->sky, HTRDR_Kext, HTRDR_SVX_MIN, HTRDR_ALL_COMPONENTS, -1/*FIXME*/, leaf);
+ kext_max = htrdr_sky_fetch_svx_voxel_property(ctx->sky, HTRDR_Kext,
+ HTRDR_SVX_MAX, HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, leaf);
+ kext_min = htrdr_sky_fetch_svx_voxel_property(ctx->sky, HTRDR_Kext,
+ HTRDR_SVX_MIN, HTRDR_ALL_COMPONENTS, ctx->iband, ctx->iquad, leaf);
CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min));
CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max));
}
@@ -200,14 +231,20 @@ vox_get(const size_t xyz[3], void* dst, void* context)
double ks_min, ks_max;
double kext_min, kext_max;
double rho_da; /* Dry air density */
+ double Ca_avg;
+ double Cs_avg;
float* pflt = dst;
ASSERT(xyz && dst && context);
+ /* Fetch the optical properties of the spectral band */
+ Ca_avg = ctx->sky->sw_bands[ctx->iband].Ca_avg;
+ Cs_avg = ctx->sky->sw_bands[ctx->iband].Cs_avg;
+
if(!ctx->sky->htcp_desc.irregular_z) {
rho_da = cloud_dry_air_density(&ctx->sky->htcp_desc, xyz);
rct = htcp_desc_RCT_at(&ctx->sky->htcp_desc, xyz[0], xyz[1], xyz[2], 0);
- ka_min = ka_max = ka = ctx->sky->Ca_avg_sw * rho_da * rct;
- ks_min = ks_max = ks = ctx->sky->Cs_avg_sw * rho_da * rct;
+ ka_min = ka_max = ka = Ca_avg * rho_da * rct;
+ ks_min = ks_max = ks = Cs_avg * rho_da * rct;
kext_min = kext_max = kext = ka + ks;
} else {
size_t ivox[3];
@@ -221,8 +258,8 @@ vox_get(const size_t xyz[3], void* dst, void* context)
rho_da = cloud_dry_air_density(&ctx->sky->htcp_desc, ivox);
rct = htcp_desc_RCT_at(&ctx->sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
- ka_min = ka_max = ka = ctx->sky->Ca_avg_sw * rho_da * rct;
- ks_min = ks_max = ks = ctx->sky->Cs_avg_sw * rho_da * rct;
+ ka_min = ka_max = ka = Ca_avg * rho_da * rct;
+ ks_min = ks_max = ks = Cs_avg * rho_da * rct;
kext_min = kext_max = kext = ka + ks;
/* Define if the SVX voxel is overlapped by 2 HTCP voxels */
@@ -236,8 +273,8 @@ vox_get(const size_t xyz[3], void* dst, void* context)
rho_da = cloud_dry_air_density(&ctx->sky->htcp_desc, ivox);
rct = htcp_desc_RCT_at(&ctx->sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
- ka = ctx->sky->Ca_avg_sw * rho_da * rct;
- ks = ctx->sky->Cs_avg_sw * rho_da * rct;
+ ka = Ca_avg * rho_da * rct;
+ ks = Cs_avg * rho_da * rct;
kext = ka + ks;
ka_min = MMIN(ka_min, ka);
@@ -332,8 +369,10 @@ setup_clouds(struct htrdr_sky* sky)
double low[3];
double upp[3];
double sz[3];
+ size_t nbands;
+ size_t iband;
res_T res = RES_OK;
- ASSERT(sky);
+ ASSERT(sky && sky->sw_bands);
res = htcp_get_desc(sky->htcp, &sky->htcp_desc);
if(res != RES_OK) {
@@ -412,26 +451,124 @@ setup_clouds(struct htrdr_sky* sky)
vox_desc.size = sizeof(float)
* HTRDR_SVX_OPS_COUNT__ * HTRDR_PROPERTIES_COUNT__;
- /* Create the octree */
- res = svx_octree_create
- (sky->htrdr->svx, low, upp, nvoxs, &vox_desc, &sky->clouds);
- if(res != RES_OK) {
+ /* Create as many cloud data structure than considered SW spectral bands */
+ nbands = htrdr_sky_get_sw_spectral_bands_count(sky);
+ sky->clouds = MEM_CALLOC(sky->htrdr->allocator, nbands, sizeof(*sky->clouds));
+ if(!sky->clouds) {
htrdr_log_err(sky->htrdr,
- "could not create the octree of the cloud properties.\n");
+ "could not create the list of per SW band cloud data structure.\n");
+ res = RES_MEM_ERR;
goto error;
}
- /* Fetch the octree descriptor for future use */
- SVX(tree_get_desc(sky->clouds, &sky->cloud_desc));
+ FOR_EACH(iband, 0, nbands) {
+ ctx.iband = iband;
+
+ /* Create the octree */
+ res = svx_octree_create
+ (sky->htrdr->svx, low, upp, nvoxs, &vox_desc, &sky->clouds[iband].octree);
+ if(res != RES_OK) {
+ htrdr_log_err(sky->htrdr, "could not create the octree of the cloud "
+ "properties for the %luth band.\n", (unsigned long)iband);
+ goto error;
+ }
+
+ /* Fetch the octree descriptor for future use */
+ SVX(tree_get_desc
+ (sky->clouds[iband].octree, &sky->clouds[iband].octree_desc));
+ }
exit:
return res;
error:
- if(sky->clouds) SVX(tree_ref_put(sky->clouds));
+ if(sky->clouds) {
+ FOR_EACH(iband, 0, nbands) {
+ if(sky->clouds[iband].octree) {
+ SVX(tree_ref_put(sky->clouds[iband].octree));
+ sky->clouds[iband].octree = NULL;
+ }
+ }
+ MEM_RM(sky->htrdr->allocator, sky->clouds);
+ }
darray_split_clear(&sky->svx2htcp_z);
goto exit;
}
+static res_T
+setup_sw_bands_properties(struct htrdr_sky* sky)
+{
+ res_T res = RES_OK;
+ size_t iband;
+ size_t nbands;
+ ASSERT(sky);
+
+ res = htgop_get_sw_spectral_intervals_CIE_XYZ(sky->htgop, sky->sw_bands_range);
+ if(res != RES_OK) goto error;
+
+ nbands = htrdr_sky_get_sw_spectral_bands_count(sky);
+ ASSERT(nbands);
+ sky->sw_bands = MEM_CALLOC
+ (sky->htrdr->allocator, nbands, sizeof(*sky->sw_bands));
+ if(!sky->sw_bands) {
+ htrdr_log_err(sky->htrdr,
+ "could not allocate the list of SW band properties.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ FOR_EACH(iband, 0, nbands) {
+ struct htgop_spectral_interval band;
+ double band_wlens[2];
+
+ HTGOP(get_sw_spectral_interval
+ (sky->htgop, iband + sky->sw_bands_range[0], &band));
+ band_wlens[0] = wavenumber_to_wavelength(band.wave_numbers[1]);
+ band_wlens[1] = wavenumber_to_wavelength(band.wave_numbers[0]);
+ ASSERT(band_wlens[0] < band_wlens[1]);
+
+ sky->sw_bands[iband].Ca_avg = htmie_compute_xsection_absorption_average
+ (sky->htmie, band_wlens, HTMIE_FILTER_LINEAR);
+ sky->sw_bands[iband].Cs_avg = htmie_compute_xsection_scattering_average
+ (sky->htmie, band_wlens, HTMIE_FILTER_LINEAR);
+ sky->sw_bands[iband].g_avg = htmie_compute_asymmetry_parameter_average
+ (sky->htmie, band_wlens, HTMIE_FILTER_LINEAR);
+ ASSERT(sky->sw_bands[iband].Ca_avg > 0);
+ ASSERT(sky->sw_bands[iband].Cs_avg > 0);
+ ASSERT(sky->sw_bands[iband].g_avg > 0);
+ }
+
+exit:
+ return res;
+error:
+ if(sky->sw_bands) {
+ MEM_RM(sky->htrdr->allocator, sky->sw_bands);
+ sky->sw_bands = NULL;
+ }
+ goto exit;
+}
+
+static void
+sample_sw_spectral_data
+ (struct htgop* htgop,
+ struct ssp_rng* rng,
+ res_T (*sample_sw_band)(const struct htgop*, const double, size_t*),
+ size_t* ispectral_band,
+ size_t* iquadrature_pt)
+{
+ struct htgop_spectral_interval specint;
+ double r1, r2;
+ res_T res = RES_OK;
+ ASSERT(htgop && rng && sample_sw_band && ispectral_band && iquadrature_pt);
+ ASSERT(ispectral_band && iquadrature_pt);
+ (void)res;
+ r1 = ssp_rng_canonical(rng);
+ r2 = ssp_rng_canonical(rng);
+ res = sample_sw_band(htgop, r1, ispectral_band);
+ ASSERT(res == RES_OK);
+ HTGOP(get_sw_spectral_interval(htgop, *ispectral_band, &specint));
+ HTGOP(spectral_interval_sample_quadrature(&specint, r2, iquadrature_pt));
+}
+
static void
release_sky(ref_T* ref)
{
@@ -439,9 +576,21 @@ release_sky(ref_T* ref)
ASSERT(ref);
sky = CONTAINER_OF(ref, struct htrdr_sky, ref);
if(sky->sun) htrdr_sun_ref_put(sky->sun);
- if(sky->clouds) SVX(tree_ref_put(sky->clouds));
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 iband;
+ FOR_EACH(iband, 0, nbands) {
+ if(sky->clouds[iband].octree) {
+ SVX(tree_ref_put(sky->clouds[iband].octree));
+ sky->clouds[iband].octree = NULL;
+ }
+ }
+ MEM_RM(sky->htrdr->allocator, sky->clouds);
+ }
darray_split_release(&sky->svx2htcp_z);
MEM_RM(sky->htrdr->allocator, sky);
}
@@ -454,10 +603,10 @@ htrdr_sky_create
(struct htrdr* htrdr,
struct htrdr_sun* sun,
const char* htcp_filename,
+ const char* htgop_filename,
const char* htmie_filename,
struct htrdr_sky** out_sky)
{
- const double band_sw[2] = {SW_WAVELENGTH_MIN, SW_WAVELENGTH_MAX};
struct htrdr_sky* sky = NULL;
res_T res = RES_OK;
ASSERT(htrdr && sun && htcp_filename && htmie_filename && out_sky);
@@ -474,12 +623,12 @@ htrdr_sky_create
sky->sun = sun;
darray_split_init(htrdr->allocator, &sky->svx2htcp_z);
+ /* Load clouds properties */
res = htcp_create(&htrdr->logger, htrdr->allocator, htrdr->verbose, &sky->htcp);
if(res != RES_OK) {
htrdr_log_err(htrdr, "could not create the loader of cloud properties.\n");
goto error;
}
-
res = htcp_load(sky->htcp, htcp_filename);
if(res != RES_OK) {
htrdr_log_err(htrdr, "error loading the cloud properties -- `%s'.\n",
@@ -487,12 +636,12 @@ htrdr_sky_create
goto error;
}
+ /* Load MIE data */
res = htmie_create(&htrdr->logger, htrdr->allocator, htrdr->verbose, &sky->htmie);
if(res != RES_OK) {
htrdr_log_err(htrdr, "could not create the loader of Mie's data.\n");
goto error;
}
-
res = htmie_load(sky->htmie, htmie_filename);
if(res != RES_OK) {
htrdr_log_err(htrdr, "error loading the Mie's data -- `%s'.\n",
@@ -500,13 +649,23 @@ htrdr_sky_create
goto error;
}
- sky->Ca_avg_sw = htmie_compute_xsection_absorption_average
- (sky->htmie, band_sw, HTMIE_FILTER_LINEAR);
- sky->Cs_avg_sw = htmie_compute_xsection_scattering_average
- (sky->htmie, band_sw, HTMIE_FILTER_LINEAR);
- sky->g_avg_sw = htmie_compute_asymmetry_parameter_average
- (sky->htmie, band_sw, HTMIE_FILTER_LINEAR);
- ASSERT(sky->Ca_avg_sw > 0 && sky->Cs_avg_sw > 0 && sky->g_avg_sw > 0);
+ /* Load the gas optical properties */
+ res = htgop_create
+ (&htrdr->logger, htrdr->allocator, htrdr->verbose, &sky->htgop);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr,
+ "could not create the loader of the gas optical properties.\n");
+ goto error;
+ }
+ res = htgop_load(sky->htgop, htgop_filename);
+ if(res != RES_OK) {
+ htrdr_log_err(htrdr, "error loading the gas optical properties -- `%s'.\n",
+ htgop_filename);
+ goto error;
+ }
+
+ res = setup_sw_bands_properties(sky);
+ if(res != RES_OK) goto error;
res = setup_clouds(sky);
if(res != RES_OK) goto error;
@@ -538,11 +697,17 @@ htrdr_sky_ref_put(struct htrdr_sky* sky)
double
htrdr_sky_fetch_particle_phase_function_asymmetry_parameter
- (const struct htrdr_sky* sky, const double wavelength)
+ (const struct htrdr_sky* sky,
+ const size_t ispectral_band,
+ const size_t iquad)
{
- ASSERT(sky && wavelength > 0);
- (void)wavelength;
- return sky->g_avg_sw;
+ size_t iband;
+ ASSERT(sky);
+ ASSERT(ispectral_band >= sky->sw_bands_range[0]);
+ ASSERT(ispectral_band <= sky->sw_bands_range[1]);
+ (void)iquad;
+ iband = ispectral_band - sky->sw_bands_range[0];
+ return sky->sw_bands[iband].g_avg;
}
double
@@ -550,42 +715,50 @@ htrdr_sky_fetch_raw_property
(const struct htrdr_sky* sky,
const enum htrdr_sky_property prop,
const int components_mask, /* Combination of htrdr_sky_component_flag */
- const double wavelength, /* FIXME Unused */
+ const size_t ispectral_band, /* Index of the spectral band */
+ const size_t iquad, /* Index of the quadrature point in the spectral band */
const double pos[3])
{
size_t ivox[3];
+ size_t iband;
+ const struct svx_tree_desc* cloud_desc;
int comp_mask = components_mask;
double k_particle = 0;
double k_gaz = 0;
ASSERT(sky && pos);
+ ASSERT(ispectral_band >= sky->sw_bands_range[0]);
+ ASSERT(ispectral_band <= sky->sw_bands_range[1]);
ASSERT(comp_mask & HTRDR_ALL_COMPONENTS);
- (void)wavelength; /* FIXME */
+ (void)iquad; /* TODO */
+
+ iband = ispectral_band - sky->sw_bands_range[0];
+ cloud_desc = &sky->clouds[iband].octree_desc;
/* 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]) {
+ 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] - 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);
+ 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] - sky->cloud_desc.lower[2])/sky->htcp_desc.vxsz_z[0]);
+ 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 n = darray_split_size_get(&sky->svx2htcp_z);
- const double sz = sky->cloud_desc.upper[2] - sky->cloud_desc.lower[2];
+ const double sz = cloud_desc->upper[2] - 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);
+ const size_t ilut = (size_t)((pos[2] - cloud_desc->lower[2])/vxsz_lut);
ivox[2] = splits[ilut].index + (pos[2] >= splits[ilut].height);
}
@@ -603,21 +776,10 @@ htrdr_sky_fetch_raw_property
rct = htcp_desc_RCT_at(&sky->htcp_desc, ivox[0], ivox[1], ivox[2], 0);
ql = rho_da * rct;
- /* FIXME do not use the wavelength yet. Simply use the average cross
- * section on the while short wave range */
-#if 1
- if(prop == HTRDR_Ka || prop == HTRDR_Kext) Ca = sky->Ca_avg_sw;
- if(prop == HTRDR_Ks || prop == HTRDR_Kext) Cs = sky->Cs_avg_sw;
-#else
- if(prop == HTRDR_Ks || prop == HTRDR_Kext) {
- Cs = htmie_fetch_xsection_scattering
- (sky->htmie, wavelength, HTMIE_FILTER_LINEAR);
- }
- if(prop == HTRDR_Ka || prop == HTRDR_Kext) {
- Ca = htmie_fetch_xsection_absorption
- (sky->htmie, wavelength, HTMIE_FILTER_LINEAR);
- }
-#endif
+ /* Use the average cross section of the current spectral band */
+ if(prop == HTRDR_Ka || prop == HTRDR_Kext) Ca = sky->sw_bands[iband].Ca_avg;
+ if(prop == HTRDR_Ks || prop == HTRDR_Kext) Cs = sky->sw_bands[iband].Cs_avg;
+
k_particle = ql*(Ca + Cs);
}
if(comp_mask & HTRDR_GAS) { /* TODO not implemented yet */ }
@@ -625,36 +787,59 @@ htrdr_sky_fetch_raw_property
}
struct svx_tree*
-htrdr_sky_get_svx_tree(struct htrdr_sky* sky)
+htrdr_sky_get_svx_tree
+ (struct htrdr_sky* sky,
+ const size_t ispectral_band,
+ const size_t iquadrature_pt)
{
+ size_t iband;
ASSERT(sky);
- return sky->clouds;
+ ASSERT(ispectral_band >= sky->sw_bands_range[0]);
+ ASSERT(ispectral_band <= sky->sw_bands_range[1]);
+ (void)iquadrature_pt;
+ iband = ispectral_band - sky->sw_bands_range[0];
+ return sky->clouds[iband].octree;
}
res_T
-htrdr_sky_dump_clouds_vtk(const struct htrdr_sky* sky, FILE* stream)
+htrdr_sky_dump_clouds_vtk
+ (const struct htrdr_sky* sky,
+ const size_t ispectral_band, /* Index of the spectral band */
+ const size_t iquad, /* Index of the quadrature point */
+ FILE* stream)
{
+ struct htgop_spectral_interval specint;
struct svx_tree_desc desc;
struct octree_data data;
const double* leaf_data;
+ size_t iband;
size_t nvertices;
size_t ncells;
size_t i;
ASSERT(sky && stream);
+ ASSERT(ispectral_band >= sky->sw_bands_range[0]);
+ ASSERT(ispectral_band <= sky->sw_bands_range[1]);
+
+ iband = ispectral_band - sky->sw_bands_range[0];
- octree_data_init(sky, &data);
- SVX(tree_get_desc(sky->clouds, &desc));
+ octree_data_init(sky, ispectral_band, iquad, &data);
+ SVX(tree_get_desc(sky->clouds[iband].octree, &desc));
ASSERT(desc.type == SVX_OCTREE);
/* Register leaf data */
- SVX(tree_for_each_leaf(sky->clouds, register_leaf, &data));
+ SVX(tree_for_each_leaf(sky->clouds[iband].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);
+ /* Fetch the spectral interval descriptor */
+ HTGOP(get_sw_spectral_interval(sky->htgop, ispectral_band, &specint));
+
/* Write headers */
fprintf(stream, "# vtk DataFile Version 2.0\n");
- fprintf(stream, "Volume\n");
+ fprintf(stream, "Clouds optical properties in [%g, %g] nanometers\n",
+ wavenumber_to_wavelength(specint.wave_numbers[1]),
+ wavenumber_to_wavelength(specint.wave_numbers[0]));
fprintf(stream, "ASCII\n");
fprintf(stream, "DATASET UNSTRUCTURED_GRID\n");
@@ -703,27 +888,33 @@ htrdr_sky_fetch_svx_property
const enum htrdr_sky_property prop,
const enum htrdr_svx_op op,
const int components_mask, /* Combination of htrdr_sky_component_flag */
- const double wavelength,
+ const size_t ispectral_band, /* Index of the spectral band */
+ const size_t iquad, /* Index of the quadrature point in the spectral band */
const double pos[3])
{
struct svx_voxel voxel = SVX_VOXEL_NULL;
+ size_t iband;
int comp_mask = components_mask;
ASSERT(sky && pos);
ASSERT(comp_mask & HTRDR_ALL_COMPONENTS);
+ ASSERT(ispectral_band >= sky->sw_bands_range[0]);
+ ASSERT(ispectral_band <= sky->sw_bands_range[1]);
+
+ iband = ispectral_band - sky->sw_bands_range[0];
/* 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]) {
+ if(pos[0] < sky->clouds[iband].octree_desc.lower[0]
+ || pos[1] < sky->clouds[iband].octree_desc.lower[1]
+ || pos[2] < sky->clouds[iband].octree_desc.lower[2]
+ || pos[0] > sky->clouds[iband].octree_desc.upper[0]
+ || pos[1] > sky->clouds[iband].octree_desc.upper[1]
+ || pos[2] > sky->clouds[iband].octree_desc.upper[2]) {
comp_mask &= ~HTRDR_PARTICLES; /* No particle */
}
- SVX(tree_at(sky->clouds, pos, NULL, NULL, &voxel));
+ SVX(tree_at(sky->clouds[iband].octree, pos, NULL, NULL, &voxel));
return htrdr_sky_fetch_svx_voxel_property
- (sky, prop, op, comp_mask, wavelength, &voxel);
+ (sky, prop, op, comp_mask, ispectral_band, iquad, &voxel);
}
double
@@ -732,7 +923,8 @@ htrdr_sky_fetch_svx_voxel_property
const enum htrdr_sky_property prop,
const enum htrdr_svx_op op,
const int components_mask,
- const double wavelength, /* FIXME Unused */
+ const size_t ispectral_band, /* Index of the spectral band */
+ const size_t iquad, /* Index of the quadrature point in the spectral band */
const struct svx_voxel* voxel)
{
const float* pflt = NULL;
@@ -743,7 +935,7 @@ htrdr_sky_fetch_svx_voxel_property
ASSERT(sky && voxel);
ASSERT((unsigned)prop < HTRDR_PROPERTIES_COUNT__);
ASSERT((unsigned)op < HTRDR_SVX_OPS_COUNT__);
- (void)sky, (void)wavelength;
+ (void)sky, (void)ispectral_band, (void)iquad;
pflt = voxel->data;
@@ -768,3 +960,73 @@ htrdr_sky_fetch_svx_voxel_property
return data;
}
+size_t
+htrdr_sky_get_sw_spectral_bands_count(const struct htrdr_sky* sky)
+{
+ ASSERT(sky);
+ return sky->sw_bands_range[1] - sky->sw_bands_range[0] + 1;
+}
+
+size_t
+htrdr_sky_get_sw_spectral_band_id
+ (const struct htrdr_sky* sky, const size_t iband)
+{
+ ASSERT(sky && iband < htrdr_sky_get_sw_spectral_bands_count(sky));
+ return sky->sw_bands_range[0] + iband;
+}
+
+res_T
+htrdr_sky_get_sw_spectral_band_bounds
+ (const struct htrdr_sky* sky,
+ const size_t iband,
+ double wavelengths[2])
+{
+ struct htgop_spectral_interval specint;
+ res_T res = RES_OK;
+ ASSERT(sky && wavelengths);
+
+ res = htgop_get_sw_spectral_interval(sky->htgop, iband, &specint);
+ if(res != RES_OK) return res;
+
+ wavelengths[0] = wavenumber_to_wavelength(specint.wave_numbers[1]);
+ wavelengths[1] = wavenumber_to_wavelength(specint.wave_numbers[0]);
+ ASSERT(wavelengths[0] < wavelengths[1]);
+ return RES_OK;
+}
+
+void
+htrdr_sky_sample_sw_spectral_data_CIE_1931_X
+ (const struct htrdr_sky* sky,
+ struct ssp_rng* rng,
+ size_t* ispectral_band,
+ size_t* iquadrature_pt)
+{
+ sample_sw_spectral_data
+ (sky->htgop, rng, htgop_sample_sw_spectral_interval_CIE_1931_X,
+ ispectral_band, iquadrature_pt);
+}
+
+void
+htrdr_sky_sample_sw_spectral_data_CIE_1931_Y
+ (const struct htrdr_sky* sky,
+ struct ssp_rng* rng,
+ size_t* ispectral_band,
+ size_t* iquadrature_pt)
+{
+ sample_sw_spectral_data
+ (sky->htgop, rng, htgop_sample_sw_spectral_interval_CIE_1931_Y,
+ ispectral_band, iquadrature_pt);
+}
+
+void
+htrdr_sky_sample_sw_spectral_data_CIE_1931_Z
+ (const struct htrdr_sky* sky,
+ struct ssp_rng* rng,
+ size_t* ispectral_band,
+ size_t* iquadrature_pt)
+{
+ sample_sw_spectral_data
+ (sky->htgop, rng, htgop_sample_sw_spectral_interval_CIE_1931_Z,
+ ispectral_band, iquadrature_pt);
+}
+
diff --git a/src/htrdr_sky.h b/src/htrdr_sky.h
@@ -43,6 +43,7 @@ enum htrdr_svx_op {
struct htrdr;
struct htrdr_sky;
struct htrdr_sun;
+struct ssp_rng;
struct svx_tree;
struct svx_voxel;
@@ -51,6 +52,7 @@ htrdr_sky_create
(struct htrdr* htrdr,
struct htrdr_sun* sun,
const char* htcp_filename,
+ const char* htgop_filename,
const char* htmie_filename,
struct htrdr_sky** sky);
@@ -65,19 +67,23 @@ htrdr_sky_ref_put
extern LOCAL_SYM double
htrdr_sky_fetch_particle_phase_function_asymmetry_parameter
(const struct htrdr_sky* sky,
- const double wavelength);
+ const size_t ispectral_band,
+ const size_t iquadrature_pt);
extern LOCAL_SYM double
htrdr_sky_fetch_raw_property
(const struct htrdr_sky* sky,
const enum htrdr_sky_property prop,
const int comp_mask, /* Combination of htrdr_sky_component_flag */
- const double wavelength,
+ const size_t ispectral_band,
+ const size_t iquadrature_pt,
const double pos[3]);
extern LOCAL_SYM struct svx_tree*
htrdr_sky_get_svx_tree
- (struct htrdr_sky* sky);
+ (struct htrdr_sky* sky,
+ const size_t ispectral_band,
+ const size_t iquadrature_pt);
extern LOCAL_SYM double
htrdr_sky_fetch_svx_property
@@ -85,7 +91,8 @@ htrdr_sky_fetch_svx_property
const enum htrdr_sky_property prop,
const enum htrdr_svx_op op,
const int comp_mask, /* Combination of htrdr_sky_component_flag */
- const double wavelength,
+ const size_t ispectral_band,
+ const size_t iquadrature_pt,
const double pos[3]);
extern LOCAL_SYM double
@@ -94,13 +101,52 @@ htrdr_sky_fetch_svx_voxel_property
const enum htrdr_sky_property prop,
const enum htrdr_svx_op op,
const int comp_mask, /* Combination of htrdr_sky_component_flag */
- const double wavelength,
+ const size_t ispectral_band,
+ const size_t iquadrature_pt,
const struct svx_voxel* voxel);
extern LOCAL_SYM res_T
htrdr_sky_dump_clouds_vtk
(const struct htrdr_sky* sky,
+ const size_t ispectral_band,
+ const size_t iquadrature_pt,
FILE* stream);
+extern LOCAL_SYM size_t
+htrdr_sky_get_sw_spectral_bands_count
+ (const struct htrdr_sky* sky);
+
+extern LOCAL_SYM size_t
+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 res_T
+htrdr_sky_get_sw_spectral_band_bounds
+ (const struct htrdr_sky* sky,
+ const size_t iband,
+ double wavelengths[2]);
+
+extern LOCAL_SYM void
+htrdr_sky_sample_sw_spectral_data_CIE_1931_X
+ (const struct htrdr_sky* sky,
+ struct ssp_rng* rng,
+ size_t* ispectral_band,
+ size_t* iquadrature_pt);
+
+extern LOCAL_SYM void
+htrdr_sky_sample_sw_spectral_data_CIE_1931_Y
+ (const struct htrdr_sky* sky,
+ struct ssp_rng* rng,
+ size_t* ispectral_band,
+ size_t* iquadrature_pt);
+
+extern LOCAL_SYM void
+htrdr_sky_sample_sw_spectral_data_CIE_1931_Z
+ (const struct htrdr_sky* sky,
+ struct ssp_rng* rng,
+ size_t* ispectral_band,
+ size_t* iquadrature_pt);
+
#endif /* HTRDR_SKY_H */
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -36,10 +36,12 @@ struct ssp_rng;
extern LOCAL_SYM double
htrdr_compute_radiance_sw
(struct htrdr* htrdr,
+ const size_t ithread,
struct ssp_rng* rng,
const double pos[3],
const double dir[3],
- const double wavelength);
+ const size_t iband, /* Index of the spectral band */
+ const size_t iquad); /* Index of the quadrature point into the band */
extern LOCAL_SYM res_T
htrdr_draw_radiance_sw
diff --git a/src/htrdr_sun.c b/src/htrdr_sun.c
@@ -250,7 +250,7 @@ htrdr_sun_get_spectral_band_bounds
}
void
-htrdr_sun_get_XYZ_spectral_bands_range
+htrdr_sun_get_CIE_XYZ_spectral_bands_range
(const struct htrdr_sun* sun, size_t band_range[2])
{
/* Bounds of the X, Y and Z functions as defined by the CIE */
diff --git a/src/htrdr_sun.h b/src/htrdr_sun.h
@@ -73,10 +73,10 @@ htrdr_sun_get_spectral_band_bounds
const size_t ispectral_band,
double bounds[2]); /* Lower and upper wavelength in nanometer */
-/* Return the ranges of the spectral bands where the XYZ color space is defined
- * XYZ in [band_min, band_max] */
+/* Return the ranges of the spectral bands where the CIE XYZ color space is
+ * defined. CIE XYZ in [band_range[0], band_range[1]] */
extern LOCAL_SYM void
-htrdr_sun_get_XYZ_spectral_bands_range
+htrdr_sun_get_CIE_XYZ_spectral_bands_range
(const struct htrdr_sun* sun,
size_t band_range[2]);
