commit cc419c395bf5eae2cef4b6805b8255e9bb52b6fd
parent bb73632b6f8ed6f92630c66f57589c12a0ac7b07
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 17 Apr 2019 14:02:25 +0200
Merge branch 'release_0.1'
Diffstat:
13 files changed, 374 insertions(+), 146 deletions(-)
diff --git a/README.md b/README.md
@@ -60,6 +60,12 @@ informations on CMake.
## Release notes
+### Version 0.1
+
+- Add the `-V` option that fixes the maximum definition of the octrees used to
+ partitioned the radiative properties of the clouds.
+- Add a per pixel estimation of the per radiative path computation time.
+
### Version 0.0.4
- Fix the computation of the surface scattering: there was a bug in how Russian
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -59,8 +59,8 @@ include_directories(
# Generate files
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 0)
-set(VERSION_PATCH 4)
+set(VERSION_MINOR 1)
+set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(HTRDR_ARGS_DEFAULT_CAMERA_POS "0,0,0")
@@ -87,7 +87,6 @@ add_subdirectory(doc)
################################################################################
# Configure and define targets
################################################################################
-
set(HTRDR_FILES_SRC
htrdr.c
htrdr_args.c
@@ -117,8 +116,7 @@ set(HTRDR_FILES_INC
set(HTRDR_FILES_INC2 # Generated files
htrdr_args.h
htrdr_version.h)
-
-set(HTDRD_FILES_DOC COPYING README.md)
+set(HTRDR_FILES_DOC COPYING README.md)
# Prepend each file in the `HTRDR_FILES_<SRC|INC>' list by `HTRDR_SOURCE_DIR'
rcmake_prepend_path(HTRDR_FILES_SRC ${HTRDR_SOURCE_DIR})
diff --git a/doc/htrdr-image.5.txt b/doc/htrdr-image.5.txt
@@ -27,12 +27,14 @@ The *htrdr-image* is a raw image file format where data are stored in plain
text. Characters after the '#' character are considered as comments and are
thus ignored as well as empty lines. The first valid line stores 2 unsigned
integers that represent the image definition, i.e. the number of pixels per
-line and per column. Then each line stores 3 pairs of floating points data
-representing the pixel color encoded in the CIE 1931 XYZ color space. The
-first, second and third pair encodes the estimated radiance of the X, Y and Z
-pixel component, respectively. The first value of each pair is the expected
-value of the estimated radiance while the second one is its associated
-standard deviation.
+line and per column. Then each line stores 4 pairs of floating points data
+representing the pixel color encoded in the CIE 1931 XYZ color space and the
+time spent per realisation to estimate the pixel. The first, second and third
+pair encodes the estimated radiance in W.sr^-1.m^-2 of the X, Y and Z pixel
+components, respectively. The first value of each pair is the expected value of
+the estimated radiance while the second one is its associated standard
+deviation. The fourth pair saves the estimate in microseconds of the per
+radiative path computation time and its standard error.
Pixels are sorted line by line, with the origin defined at the top left corner
of the image. With an image definition of N by M pixels, with N the number of
@@ -56,12 +58,13 @@ GRAMMAR
<width> ::= INTEGER
<height> ::= INTEGER
-<pixel> ::= <X> <Y> <Z>
+<pixel> ::= <X> <Y> <Z> <time>
<X> ::= <estimate>
<Y> ::= <estimate>
<Z> ::= <estimate>
+<time> ::= <estimate>
-<estimate> ::= <exoected-value> <standard-error>
+<estimate> ::= <expected-value> <standard-error>
<expected-value> ::= REAL
<standard-error> ::= REAL
-------
@@ -80,30 +83,30 @@ the index of the line and the column of the pixel in the image, respectively.
800 600 # Image definition
# Pixels of the 1st line
-0.0002559 2.90278e-05 0.0003754 4.48780e-05 0.0003200 3.16818e-05 # (1,1)
-0.0002956 3.37333e-05 0.0003395 4.16836e-05 0.0003384 4.60855e-05 # (2,1)
-0.0003761 5.43206e-05 0.0003130 3.48458e-05 0.0003380 3.32255e-05 # (3,1)
+2.55e-4 2.90e-5 3.75e-4 4.48e-5 3.20e-4 3.16e-5 306.484 259.723 # (1,1)
+2.95e-4 3.37e-5 3.39e-4 4.16e-5 3.38e-4 4.60e-5 18.3633 2.66317 # (2,1)
+3.76e-4 5.43e-5 3.13e-4 3.48e-5 3.38e-4 3.32e-5 19.6252 2.67015 # (3,1)
...
-0.0007135 0.000114694 0.0007661 0.000135251 0.0007978 0.000126122 # (799,1)
-0.0006591 0.000114968 0.0007473 0.000141289 0.0004398 7.33230e-05 # (800,1)
+7.13e-4 1.14e-4 7.66e-4 1.35e-4 7.97e-4 1.26e-4 119.820 93.7820 # (799,1)
+6.59e-4 1.14e-4 7.47e-4 1.41e-4 4.39e-4 7.33e-5 24.8655 2.46348 # (800,1)
# Pixels of the 2nd line
-0.0003330 6.02808e-05 0.0004218 7.66163e-05 0.0003446 3.81472e-05 # (1,2)
-0.0003504 4.93154e-05 0.0003232 2.52040e-05 0.0003036 2.42803e-05 # (2,2)
-0.0002720 4.69853e-05 0.0003413 4.12614e-05 0.0002525 2.06930e-05 # (3,2)
+3.33e-4 6.02e-5 4.21e-4 7.66e-5 3.44e-4 3.81e-5 19.4580 2.50692 # (1,2)
+3.50e-4 4.93e-5 3.23e-4 2.52e-5 3.03e-4 2.42e-5 102.566 81.2936 # (2,2)
+2.72e-4 4.69e-5 3.41e-4 4.12e-5 2.52e-4 2.06e-5 25.5801 5.37736 # (3,2)
...
-0.0007528 0.000131807 0.0008919 0.000184047 0.0005487 0.000130065 # (799,2)
-0.0006822 0.000142027 0.0006618 7.85956e-05 0.0004443 5.99220e-05 # (800,2)
+7.52e-4 1.31e-4 8.91e-4 1.84e-4 5.48e-4 1.30e-4 46.5418 12.4728 # (799,2)
+6.82e-4 1.42e-4 6.61e-4 7.85e-5 4.44e-4 5.99e-5 59.8728 32.1468 # (800,2)
...
# Pixels of the 600th line
-0.0002697 7.44208e-05 0.0002310 2.56492e-05 0.0001958 2.30952e-05 # (1,600)
-0.0004325 0.000125897 0.0002222 2.22320e-05 0.0002047 2.60948e-05 # (2,600)
-0.0002782 5.81477e-05 0.0002756 4.99110e-05 0.0002172 3.30041e-05 # (3,600)
+2.69e-4 7.44e-5 2.31e-4 2.56e-5 1.95e-4 2.30e-5 43.8242 15.0047 # (1,600)
+4.32e-4 1.25e-4 2.22e-4 2.22e-5 2.04e-4 2.60e-5 25.5498 1.73942 # (2,600)
+2.78e-4 5.81e-5 2.75e-4 4.99e-5 2.17e-4 3.30e-5 38.4448 7.16199 # (3,600)
...
-0.0003548 4.32117e-05 0.0003070 3.80437e-05 0.0002387 2.49820e-05 # (799,600)
-0.0003078 2.61915e-05 0.0004609 0.000113663 0.0002699 4.29776e-05 # (800,600)
+3.54e-4 4.32e-5 3.07e-4 3.80e-5 2.38e-4 2.49e-5 102.893 36.9865 # (799,600)
+3.07e-4 2.61e-5 4.60e-4 1.13e-4 2.69e-4 4.29e-5 42.75070 11.913 # (800,600)
------
NOTES
diff --git a/doc/htrdr.1.txt.in b/doc/htrdr.1.txt.in
@@ -143,7 +143,7 @@ OPTIONS
*-i* <__image-parameter__:...>::
Define the image to render. Available image parameters are:
- **def**=**_width_**,**_height_**;;
+ **def**=**_width_**x**_height_**;;
Definition of the image. By default the image definition is
@HTRDR_ARGS_DEFAULT_IMG_WIDTH@x@HTRDR_ARGS_DEFAULT_IMG_HEIGHT@.
@@ -175,6 +175,11 @@ OPTIONS
Hint on the number of threads to use. By default use as many threads as CPU
cores.
+*-V* **_x_**,**_y_**,**_z_**::
+ Define the maximum definition of the acceleration data structure that
+ partitions the cloud properties. By default the finest definition is the
+ definition of the submitted *htcp*(5) file.
+
*-v*::
Make *htrdr* verbose.
diff --git a/src/htrdr.c b/src/htrdr.c
@@ -102,6 +102,7 @@ static res_T
dump_accum_buffer
(struct htrdr* htrdr,
struct htrdr_buffer* buf,
+ struct htrdr_accum* time_acc, /* May be NULL */
const char* stream_name,
FILE* stream)
{
@@ -112,34 +113,34 @@ dump_accum_buffer
(void)stream_name;
htrdr_buffer_get_layout(buf, &layout);
- if(layout.elmt_size != sizeof(struct htrdr_accum[3])/*#channels*/
- || layout.alignment < ALIGNOF(struct htrdr_accum[3])) {
+ if(layout.elmt_size != sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__])
+ || layout.alignment < ALIGNOF(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__])) {
htrdr_log_err(htrdr,
"%s: invalid buffer layout. "
- "The pixel size must be the size of an accumulator.\n",
- FUNC_NAME);
+ "The pixel size must be the size of %lu accumulators.\n",
+ FUNC_NAME, (unsigned long)HTRDR_ESTIMATES_COUNT__);
res = RES_BAD_ARG;
goto error;
}
fprintf(stream, "%lu %lu\n", layout.width, layout.height);
+
+ if(time_acc) *time_acc = HTRDR_ACCUM_NULL;
FOR_EACH(y, 0, layout.height) {
FOR_EACH(x, 0, layout.width) {
const struct htrdr_accum* accums = htrdr_buffer_at(buf, x, y);
int i;
- FOR_EACH(i, 0, 3) {
- const double N = (double)accums[i].nweights;
- double E = 0;
- double V = 0;
- double SE = 0;
-
- if(accums[i].nweights) {
- E = accums[i].sum_weights / N;
- V = MMAX(accums[i].sum_weights_sqr / N - E*E, 0);
- SE = sqrt(V/N);
- }
+ FOR_EACH(i, 0, HTRDR_ESTIMATES_COUNT__) {
+ double E, SE;
+
+ htrdr_accum_get_estimation(&accums[i], &E, &SE);
fprintf(stream, "%g %g ", E, SE);
}
+ if(time_acc) {
+ time_acc->sum_weights += accums[HTRDR_ESTIMATE_TIME].sum_weights;
+ time_acc->sum_weights_sqr += accums[HTRDR_ESTIMATE_TIME].sum_weights_sqr;
+ time_acc->nweights += accums[HTRDR_ESTIMATE_TIME].nweights;
+ }
fprintf(stream, "\n");
}
fprintf(stream, "\n");
@@ -460,6 +461,9 @@ htrdr_init
htrdr->spp = args->image.spp;
htrdr->width = args->image.definition[0];
htrdr->height = args->image.definition[1];
+ htrdr->grid_max_definition[0] = args->grid_max_definition[0];
+ htrdr->grid_max_definition[1] = args->grid_max_definition[1];
+ htrdr->grid_max_definition[2] = args->grid_max_definition[2];
res = mem_init_regular_allocator(&htrdr->svx_allocator);
if(res != RES_OK) goto error;
@@ -524,12 +528,15 @@ htrdr_init
/* Create the image buffer only on the master process; the image parts
* rendered by the processes are gathered onto the master process. */
if(!htrdr->dump_vtk && htrdr->mpi_rank == 0) {
+ /* 4 accums: X, Y, Z components and one more for the per realisation time */
+ const size_t pixsz = sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]);
+ const size_t pixal = ALIGNOF(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]);
res = htrdr_buffer_create(htrdr,
args->image.definition[0], /* Width */
args->image.definition[1], /* Height */
- args->image.definition[0]*sizeof(struct htrdr_accum[3]), /* Pitch */
- sizeof(struct htrdr_accum[3]),
- ALIGNOF(struct htrdr_accum[3]), /* Alignment */
+ args->image.definition[0] * pixsz, /* Pitch */
+ pixsz, /* Size of a pixel */
+ pixal, /* Alignment of a pixel */
&htrdr->buf);
if(res != RES_OK) goto error;
}
@@ -626,9 +633,16 @@ htrdr_run(struct htrdr* htrdr)
htrdr->height, htrdr->spp, htrdr->buf);
if(res != RES_OK) goto error;
if(htrdr->mpi_rank == 0) {
- res = dump_accum_buffer
- (htrdr, htrdr->buf, str_cget(&htrdr->output_name), htrdr->output);
+ struct htrdr_accum path_time_acc = HTRDR_ACCUM_NULL;
+ double E, SE;
+
+ res = dump_accum_buffer(htrdr, htrdr->buf, &path_time_acc,
+ str_cget(&htrdr->output_name), htrdr->output);
if(res != RES_OK) goto error;
+
+ htrdr_accum_get_estimation(&path_time_acc, &E, &SE);
+ htrdr_log(htrdr,
+ "Time per radiative path (in micro seconds): %g +/- %g\n", E, SE);
}
}
exit:
diff --git a/src/htrdr.h b/src/htrdr.h
@@ -59,6 +59,7 @@ struct htrdr {
FILE* output;
struct str output_name;
+ unsigned grid_max_definition[3]; /* Max definition of the acceleration grids */
unsigned nthreads; /* #threads of the process */
int dump_vtk; /* Dump octree VTK */
int cache_grids; /* Use/Precompute grid caches */
diff --git a/src/htrdr_args.c b/src/htrdr_args.c
@@ -79,6 +79,9 @@ 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(
+" -V X,Y,Z maximum definition of the cloud acceleration grids along\n"
+" the 3 axis. By default use the definition of the clouds\n");
+ printf(
" -v make the program verbose.\n");
printf(
" --version display version information and exit.\n");
@@ -319,6 +322,38 @@ error:
goto exit;
}
+static res_T
+parse_grid_definition(struct htrdr_args* args, const char* str)
+{
+ unsigned def[3];
+ size_t len;
+ res_T res = RES_OK;
+ ASSERT(args && str);
+
+ res = cstr_to_list_uint(str, ',', def, &len, 3);
+ if(res == RES_OK && len != 3) res = RES_BAD_ARG;
+ if(res != RES_OK) {
+ fprintf(stderr, "Invalid grid definition `%s'.\n", str);
+ goto error;
+ }
+
+ if(!def[0] || !def[1] || !def[2]) {
+ fprintf(stderr,
+ "Invalid null grid definition {%u, %u, %u}.\n", SPLIT3(def));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ args->grid_max_definition[0] = def[0];
+ args->grid_max_definition[1] = def[1];
+ args->grid_max_definition[2] = def[2];
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -343,7 +378,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
}
}
- while((opt = getopt(argc, argv, "a:C:c:D:de:fGg:hi:m:o:RrT:t:v")) != -1) {
+ while((opt = getopt(argc, argv, "a:C:c:D:de:fGg:hi:m:o:RrT:t:V:v")) != -1) {
switch(opt) {
case 'a': args->filename_gas = optarg; break;
case 'C':
@@ -383,6 +418,7 @@ htrdr_args_init(struct htrdr_args* args, int argc, char** argv)
res = cstr_to_uint(optarg, &args->nthreads);
if(res == RES_OK && !args->nthreads) res = RES_BAD_ARG;
break;
+ case 'V': res = parse_grid_definition(args, optarg); break;
case 'v': args->verbose = 1; break;
default: res = RES_BAD_ARG; break;
}
diff --git a/src/htrdr_args.h.in b/src/htrdr_args.h.in
@@ -16,6 +16,7 @@
#ifndef HTRDR_ARGS_H
#define HTRDR_ARGS_H
+#include <limits.h>
#include <rsys/rsys.h>
struct htrdr_args {
@@ -48,6 +49,7 @@ struct htrdr_args {
double sun_elevation; /* In degrees */
double optical_thickness; /* Threshold used during octree building */
double ground_reflectivity; /* Reflectivity of the ground */
+ unsigned grid_max_definition[3]; /* Maximum definition of the grid */
unsigned nthreads; /* Hint on the number of threads to use */
int force_overwriting;
@@ -83,6 +85,7 @@ struct htrdr_args {
90, /* Sun elevation */ \
@HTRDR_ARGS_DEFAULT_OPTICAL_THICKNESS_THRESHOLD@, /* Optical thickness */ \
@HTRDR_ARGS_DEFAULT_GROUND_REFLECTIVITY@, /* Ground reflectivity */ \
+ {UINT_MAX, UINT_MAX, UINT_MAX}, /* Maximum definition of the grid */ \
(unsigned)~0, /* #threads */ \
0, /* Force overwriting */ \
0, /* dump VTK */ \
diff --git a/src/htrdr_c.h b/src/htrdr_c.h
@@ -32,6 +32,14 @@ enum htrdr_mpi_message {
HTRDR_MPI_TILE_DATA
};
+enum htrdr_estimate {
+ HTRDR_ESTIMATE_X,
+ HTRDR_ESTIMATE_Y,
+ HTRDR_ESTIMATE_Z,
+ HTRDR_ESTIMATE_TIME, /* Time per realisation */
+ HTRDR_ESTIMATES_COUNT__
+};
+
struct htrdr;
#define SW_WAVELENGTH_MIN 380 /* In nanometer */
diff --git a/src/htrdr_camera.c b/src/htrdr_camera.c
@@ -72,7 +72,7 @@ htrdr_camera_create
ref_init(&cam->ref);
cam->htrdr = htrdr;
- if(fov <= 0) {
+ if(fov <= 0) {
htrdr_log_err(htrdr, "invalid horizontal camera field of view `%g'\n", fov);
res = RES_BAD_ARG;
goto error;
diff --git a/src/htrdr_draw_radiance_sw.c b/src/htrdr_draw_radiance_sw.c
@@ -92,7 +92,7 @@ tile_create(struct mem_allocator* allocator)
const size_t tile_sz =
sizeof(struct tile) - sizeof(struct htrdr_accum)/*rm dummy accum*/;
const size_t buf_sz = /* Flexiblbe array element */
- TILE_SIZE*TILE_SIZE*sizeof(struct htrdr_accum)*3/*#channels*/;
+ TILE_SIZE*TILE_SIZE*sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]);
ASSERT(allocator);
tile = MEM_ALLOC(allocator, tile_sz+buf_sz);
@@ -135,7 +135,7 @@ tile_at
const size_t y) /* In tile space */
{
ASSERT(tile && x < TILE_SIZE && y < TILE_SIZE);
- return tile->data.accums + (y*TILE_SIZE + x)*3/*#channels*/;
+ return tile->data.accums + (y*TILE_SIZE + x) * HTRDR_ESTIMATES_COUNT__;
}
static void
@@ -150,6 +150,7 @@ write_tile_data(struct htrdr_buffer* buf, const struct tile_data* tile_data)
htrdr_buffer_get_layout(buf, &layout);
buf_mem = htrdr_buffer_get_data(buf);
+ ASSERT(layout.elmt_size == sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]));
/* Compute the row/column of the tile origin into the buffer */
icol = tile_data->x * (size_t)TILE_SIZE;
@@ -163,9 +164,9 @@ write_tile_data(struct htrdr_buffer* buf, const struct tile_data* tile_data)
FOR_EACH(irow_tile, 0, nrows_tile) {
char* buf_row = buf_mem + (irow + irow_tile) * layout.pitch;
const struct htrdr_accum* tile_row =
- tile_data->accums + irow_tile*TILE_SIZE*3/*#channels*/;
- memcpy(buf_row + icol*sizeof(struct htrdr_accum)*3, tile_row,
- ncols_tile*sizeof(struct htrdr_accum)*3/*#channels*/);
+ tile_data->accums + irow_tile*TILE_SIZE*HTRDR_ESTIMATES_COUNT__;
+ memcpy(buf_row + icol*sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]),
+ tile_row, ncols_tile*sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]));
}
}
@@ -391,7 +392,7 @@ mpi_gather_tiles
/* Compute the size of the tile_data */
const size_t msg_sz =
sizeof(struct tile_data) - sizeof(struct htrdr_accum)/*dummy*/
- + TILE_SIZE*TILE_SIZE*sizeof(struct htrdr_accum)*3/*#channels*/;
+ + TILE_SIZE*TILE_SIZE*sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__]);
struct list_node* node = NULL;
struct tile* tile = NULL;
@@ -480,22 +481,32 @@ draw_tile
/* Fetch and reset the pixel accumulator */
pix_accums = tile_at(tile, ipix_tile[0], ipix_tile[1]);
+ pix_accums[HTRDR_ESTIMATE_TIME] = HTRDR_ACCUM_NULL; /* Reset time per radiative path */
/* Compute the pixel coordinate */
ipix[0] = tile_org[0] + ipix_tile[0];
ipix[1] = tile_org[1] + ipix_tile[1];
FOR_EACH(ichannel, 0, 3) {
+ /* Check that the X, Y and Z estimates are stored in accumulators 0, 1 et
+ * 2, respectively */
+ STATIC_ASSERT
+ ( HTRDR_ESTIMATE_X == 0
+ && HTRDR_ESTIMATE_Y == 1
+ && HTRDR_ESTIMATE_Z == 2,
+ Unexpected_htrdr_estimate_enumerate);
size_t isamp;
- pix_accums[ichannel] = HTRDR_ACCUM_NULL;
+ pix_accums[ichannel] = HTRDR_ACCUM_NULL;
FOR_EACH(isamp, 0, spp) {
+ struct time t0, t1;
double pix_samp[2];
double ray_org[3];
double ray_dir[3];
double weight;
size_t iband;
size_t iquad;
+ double usec;
/* Sample a position into the pixel, in the normalized image plane */
pix_samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
@@ -523,14 +534,22 @@ draw_tile
}
/* Compute the radiance that reach the pixel through the ray */
+ time_current(&t0);
weight = htrdr_compute_radiance_sw
(htrdr, ithread, rng, ray_org, ray_dir, iband, iquad);
ASSERT(weight >= 0);
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
- /* Update the pixel accumulator */
+ /* Update the pixel accumulator of the current channel */
pix_accums[ichannel].sum_weights += weight;
pix_accums[ichannel].sum_weights_sqr += weight*weight;
pix_accums[ichannel].nweights += 1;
+
+ /* Update the pixel accumulator of per realisation time */
+ pix_accums[HTRDR_ESTIMATE_TIME].sum_weights += usec;
+ pix_accums[HTRDR_ESTIMATE_TIME].sum_weights_sqr += usec*usec;
+ pix_accums[HTRDR_ESTIMATE_TIME].nweights += 1;
}
}
}
@@ -724,12 +743,12 @@ htrdr_draw_radiance_sw
ASSERT(layout.width || layout.height || layout.elmt_size);
ASSERT(layout.width == width && layout.height == height);
- if(layout.elmt_size != sizeof(struct htrdr_accum[3])/*#channels*/
- || layout.alignment < ALIGNOF(struct htrdr_accum[3])) {
+ if(layout.elmt_size != sizeof(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__])
+ || layout.alignment < ALIGNOF(struct htrdr_accum[HTRDR_ESTIMATES_COUNT__])) {
htrdr_log_err(htrdr,
"%s: invalid buffer layout. "
- "The pixel size must be the size of 3 * accumulators.\n",
- FUNC_NAME);
+ "The pixel size must be the size of %lu accumulators.\n",
+ FUNC_NAME, (unsigned long)HTRDR_ESTIMATES_COUNT__);
res = RES_BAD_ARG;
goto error;
}
diff --git a/src/htrdr_sky.c b/src/htrdr_sky.c
@@ -207,6 +207,21 @@ struct htrdr_sky {
/*******************************************************************************
* Helper function
******************************************************************************/
+static INLINE int
+aabb_intersect
+ (const double aabb0_low[3],
+ const double aabb0_upp[3],
+ const double aabb1_low[3],
+ const double aabb1_upp[3])
+{
+ ASSERT(aabb0_low[0] < aabb0_upp[0] && aabb1_low[0] < aabb1_upp[0]);
+ ASSERT(aabb0_low[1] < aabb0_upp[1] && aabb1_low[1] < aabb1_upp[1]);
+ ASSERT(aabb0_low[2] < aabb0_upp[2] && aabb1_low[2] < aabb1_upp[2]);
+ return !(aabb0_upp[0] < aabb1_low[0]) && !(aabb0_low[0] > aabb1_upp[0])
+ && !(aabb0_upp[1] < aabb1_low[1]) && !(aabb0_low[1] > aabb1_upp[1])
+ && !(aabb0_upp[2] < aabb1_low[2]) && !(aabb0_low[2] > aabb1_upp[2]);
+}
+
static void
log_svx_memory_usage(struct htrdr* htrdr)
{
@@ -222,7 +237,7 @@ log_svx_memory_usage(struct htrdr* htrdr)
memsz = MEM_ALLOCATED_SIZE(&htrdr->svx_allocator);
if((ngigas = memsz / GIGA_BYTE) != 0) {
- len = (size_t)snprintf(dst, available_space,
+ len = (size_t)snprintf(dst, available_space,
"%lu GB ", (unsigned long)ngigas);
CHK(len < available_space);
dst += len;
@@ -230,7 +245,7 @@ log_svx_memory_usage(struct htrdr* htrdr)
memsz -= ngigas * GIGA_BYTE;
}
if((nmegas = memsz / MEGA_BYTE) != 0) {
- len = (size_t)snprintf(dst, available_space,
+ len = (size_t)snprintf(dst, available_space,
"%lu MB ", (unsigned long)nmegas);
CHK(len < available_space);
dst += len;
@@ -245,7 +260,7 @@ log_svx_memory_usage(struct htrdr* htrdr)
memsz -= nkilos * KILO_BYTE;
}
if(memsz) {
- len = (size_t)snprintf(dst, available_space,
+ len = (size_t)snprintf(dst, available_space,
"%lu Byte%s", (unsigned long)memsz, memsz > 1 ? "s" : "");
CHK(len < available_space);
}
@@ -559,6 +574,11 @@ cloud_vox_get_particle
float vox[],
const struct build_tree_context* ctx)
{
+ const struct htcp_desc* htcp_desc;
+ size_t ivox[3];
+ size_t igrid_low[3], igrid_upp[3];
+ double vox_low[3], vox_upp[3];
+ double low[3], upp[3];
double rct;
double ka, ks, kext;
double ka_min, ka_max;
@@ -567,63 +587,107 @@ cloud_vox_get_particle
double rho_da; /* Dry air density */
double Ca_avg;
double Cs_avg;
+ double ipart;
size_t i;
ASSERT(xyz && vox && ctx);
i = ctx->iband - ctx->sky->sw_bands_range[0];
+ htcp_desc = &ctx->sky->htcp_desc;
/* Fetch the optical properties of the spectral band */
Ca_avg = ctx->sky->sw_bands[i].Ca_avg;
Cs_avg = ctx->sky->sw_bands[i].Cs_avg;
+ /* Compute the AABB of the SVX voxel */
+ vox_low[0] = (double)xyz[0] * ctx->vxsz[0] + htcp_desc->lower[0];
+ vox_low[1] = (double)xyz[1] * ctx->vxsz[1] + htcp_desc->lower[1];
+ vox_low[2] = (double)xyz[2] * ctx->vxsz[2] + htcp_desc->lower[2];
+ vox_upp[0] = vox_low[0] + ctx->vxsz[0];
+ vox_upp[1] = vox_low[1] + ctx->vxsz[1];
+ vox_upp[2] = vox_low[2] + ctx->vxsz[2];
+
+ /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
+ * the SVX voxel in X and Y */
+ low[0] = (vox_low[0]-htcp_desc->lower[0])/htcp_desc->vxsz_x;
+ low[1] = (vox_low[1]-htcp_desc->lower[1])/htcp_desc->vxsz_x;
+ upp[0] = (vox_upp[0]-htcp_desc->lower[0])/htcp_desc->vxsz_y;
+ upp[1] = (vox_upp[1]-htcp_desc->lower[1])/htcp_desc->vxsz_y;
+ igrid_low[0] = (size_t)low[0];
+ igrid_low[1] = (size_t)low[1];
+ igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0);
+ igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0);
+ ASSERT(igrid_low[0] <= igrid_upp[0]);
+ ASSERT(igrid_low[1] <= igrid_upp[1]);
+
if(!ctx->sky->htcp_desc.irregular_z) { /* 1 LES voxel <=> 1 SVX voxel */
- 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 = Ca_avg * rho_da * rct;
- ks_min = ks_max = ks = Cs_avg * rho_da * rct;
- kext_min = kext_max = kext = ka + ks;
- } else { /* A SVX voxel can be overlapped by 2 LES voxels */
- double vox_low[3], vox_upp[3];
+ /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
+ * the SVX voxel along the Z axis */
+ low[2] = (vox_low[2]-htcp_desc->lower[2])/htcp_desc->vxsz_z[0];
+ upp[2] = (vox_upp[2]-htcp_desc->lower[2])/htcp_desc->vxsz_z[0];
+ igrid_low[2] = (size_t)low[2];
+ igrid_upp[2] = (size_t)upp[2] - (modf(upp[2], &ipart)==0);
+ ASSERT(igrid_low[2] <= igrid_upp[2]);
+
+ /* Prepare the iteration over the grid voxels overlapped by the SVX voxel */
+ ka_min = ks_min = kext_min = DBL_MAX;
+ ka_max = ks_max = kext_max =-DBL_MAX;
+
+ /* Iterate over the grid voxels overlapped by the SVX voxel */
+ FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
+ FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
+ FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) {
+ /* Compute the radiative properties */
+ rho_da = cloud_dry_air_density(htcp_desc, ivox);
+ rct = htcp_desc_RCT_at(htcp_desc, ivox[0], ivox[1], ivox[2], 0);
+ ka = Ca_avg * rho_da * rct;
+ ks = Cs_avg * rho_da * rct;
+ kext = ka + ks;
+ /* Update the boundaries of the radiative properties */
+ ka_min = MMIN(ka_min, ka);
+ ka_max = MMAX(ka_max, ka);
+ ks_min = MMIN(ks_min, ks);
+ ks_max = MMAX(ks_max, ks);
+ kext_min = MMIN(kext_min, kext);
+ kext_max = MMAX(kext_max, kext);
+ #ifndef NDEBUG
+ {
+ double tmp_low[3], tmp_upp[3];
+ htcp_desc_get_voxel_aabb
+ (&ctx->sky->htcp_desc, ivox[0], ivox[1], ivox[2], tmp_low, tmp_upp);
+ ASSERT(aabb_intersect(tmp_low, tmp_upp, vox_low, vox_upp));
+ }
+ #endif
+ }
+ }
+ }
+ } else {
double pos_z;
- size_t ivox_z_prev;
size_t ilut_low, ilut_upp;
size_t ilut;
-
- /* Compute the AABB of the SVX voxel */
- vox_low[0] = (double)xyz[0] * ctx->vxsz[0] + ctx->sky->htcp_desc.lower[0];
- vox_low[1] = (double)xyz[1] * ctx->vxsz[1] + ctx->sky->htcp_desc.lower[1];
- vox_low[2] = (double)xyz[2] * ctx->vxsz[2] + ctx->sky->htcp_desc.lower[2];
- vox_upp[0] = vox_low[0] + ctx->vxsz[0];
- vox_upp[1] = vox_low[1] + ctx->vxsz[1];
- vox_upp[2] = vox_low[2] + ctx->vxsz[2];
+ size_t ivox_z_prev = SIZE_MAX;
/* Compute the *inclusive* bounds of the indices of the LUT cells
* overlapped by the SVX voxel */
- ilut_low = (size_t)
- ((vox_low[2] - ctx->sky->htcp_desc.lower[2]) / ctx->sky->lut_cell_sz);
- ilut_upp = (size_t)
- ((vox_upp[2] - ctx->sky->htcp_desc.lower[2]) / ctx->sky->lut_cell_sz);
+ ilut_low = (size_t)((vox_low[2]-htcp_desc->lower[2])/ctx->sky->lut_cell_sz);
+ ilut_upp = (size_t)((vox_upp[2]-htcp_desc->lower[2])/ctx->sky->lut_cell_sz);
ASSERT(ilut_low <= ilut_upp);
/* Prepare the iteration over the LES voxels overlapped by the SVX voxel */
- ivox_z_prev = SIZE_MAX;
ka_min = ks_min = kext_min = DBL_MAX;
ka_max = ks_max = kext_max =-DBL_MAX;
+ ivox_z_prev = SIZE_MAX;
pos_z = vox_low[2];
ASSERT(pos_z >= (double)ilut_low * ctx->sky->lut_cell_sz);
ASSERT(pos_z <= (double)ilut_upp * ctx->sky->lut_cell_sz);
/* Iterate over the LUT cells overlapped by the voxel */
FOR_EACH(ilut, ilut_low, ilut_upp+1) {
- size_t ivox[3];
const struct split* split = darray_split_cdata_get(&ctx->sky->svx2htcp_z)+ilut;
ASSERT(ilut < darray_split_size_get(&ctx->sky->svx2htcp_z));
- ivox[0] = xyz[0];
- ivox[1] = xyz[1];
ivox[2] = pos_z <= split->height ? split->index : split->index + 1;
if(ivox[2] >= ctx->sky->htcp_desc.spatial_definition[2]
- && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
+ && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
ivox[2] = split->index;
}
@@ -638,20 +702,25 @@ cloud_vox_get_particle
if(ivox[2] == ivox_z_prev) continue;
ivox_z_prev = ivox[2];
- /* Compute the radiative properties */
- 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 = Ca_avg * rho_da * rct;
- ks = Cs_avg * rho_da * rct;
- kext = ka + ks;
-
- /* Update the boundaries of the radiative properties */
- ka_min = MMIN(ka_min, ka);
- ka_max = MMAX(ka_max, ka);
- ks_min = MMIN(ks_min, ks);
- ks_max = MMAX(ks_max, ks);
- kext_min = MMIN(kext_min, kext);
- kext_max = MMAX(kext_max, kext);
+ FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
+ FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
+
+ /* Compute the radiative properties */
+ rho_da = cloud_dry_air_density(htcp_desc, ivox);
+ rct = htcp_desc_RCT_at(htcp_desc, ivox[0], ivox[1], ivox[2], 0);
+ ka = Ca_avg * rho_da * rct;
+ ks = Cs_avg * rho_da * rct;
+ kext = ka + ks;
+
+ /* Update the boundaries of the radiative properties */
+ ka_min = MMIN(ka_min, ka);
+ ka_max = MMAX(ka_max, ka);
+ ks_min = MMIN(ks_min, ks);
+ ks_max = MMAX(ks_max, ks);
+ kext_min = MMIN(kext_min, kext);
+ kext_max = MMAX(kext_max, kext);
+ }
+ }
}
}
@@ -678,58 +747,98 @@ cloud_vox_get_gas
float vox[],
const struct build_tree_context* ctx)
{
+ const struct htcp_desc* htcp_desc;
+ size_t ivox[3];
+ size_t igrid_low[3], igrid_upp[3];
struct htgop_layer layer;
struct htgop_layer_sw_spectral_interval band;
size_t ilayer;
size_t layer_range[2];
double x_h2o_range[2];
+ double low[3], upp[3];
double vox_low[3], vox_upp[3]; /* Voxel AABB */
double ka_min, ka_max;
double ks_min, ks_max;
double kext_min, kext_max;
+ double x_h2o;
+ double ipart;
ASSERT(xyz && vox && ctx);
+ htcp_desc = &ctx->sky->htcp_desc;
+
/* Compute the AABB of the SVX voxel */
- vox_low[0] = (double)xyz[0] * ctx->vxsz[0] + ctx->sky->htcp_desc.lower[0];
- vox_low[1] = (double)xyz[1] * ctx->vxsz[1] + ctx->sky->htcp_desc.lower[1];
- vox_low[2] = (double)xyz[2] * ctx->vxsz[2] + ctx->sky->htcp_desc.lower[2];
+ vox_low[0] = (double)xyz[0] * ctx->vxsz[0] + htcp_desc->lower[0];
+ vox_low[1] = (double)xyz[1] * ctx->vxsz[1] + htcp_desc->lower[1];
+ vox_low[2] = (double)xyz[2] * ctx->vxsz[2] + htcp_desc->lower[2];
vox_upp[0] = vox_low[0] + ctx->vxsz[0];
vox_upp[1] = vox_low[1] + ctx->vxsz[1];
vox_upp[2] = vox_low[2] + ctx->vxsz[2];
+ /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
+ * the SVX voxel in X and Y */
+ low[0] = (vox_low[0]-htcp_desc->lower[0])/htcp_desc->vxsz_x;
+ low[1] = (vox_low[1]-htcp_desc->lower[1])/htcp_desc->vxsz_x;
+ upp[0] = (vox_upp[0]-htcp_desc->lower[0])/htcp_desc->vxsz_y;
+ upp[1] = (vox_upp[1]-htcp_desc->lower[1])/htcp_desc->vxsz_y;
+ igrid_low[0] = (size_t)low[0];
+ igrid_low[1] = (size_t)low[1];
+ igrid_upp[0] = (size_t)upp[0] - (modf(upp[0], &ipart)==0);
+ igrid_upp[1] = (size_t)upp[1] - (modf(upp[1], &ipart)==0);
+ ASSERT(igrid_low[0] <= igrid_upp[0]);
+ ASSERT(igrid_low[1] <= igrid_upp[1]);
+
/* Define the xH2O range from the LES data */
if(!ctx->sky->htcp_desc.irregular_z) { /* 1 LES voxel <=> 1 SVX voxel */
- double x_h2o;
- x_h2o = cloud_water_vapor_molar_fraction(&ctx->sky->htcp_desc, xyz);
- x_h2o_range[0] = x_h2o_range[1] = x_h2o;
-#ifndef NDEBUG
- {
- double low[3], upp[3];
- htcp_desc_get_voxel_aabb
- (&ctx->sky->htcp_desc, xyz[0], xyz[1], xyz[2], low, upp);
- ASSERT(d3_eq_eps(low, vox_low, 1.e-6));
- ASSERT(d3_eq_eps(upp, vox_upp, 1.e-6));
+ /* Compute the *inclusive* bounds of the indices of cloud grid overlapped by
+ * the SVX voxel in Z */
+ low[2] = (vox_low[2]-htcp_desc->lower[2])/htcp_desc->vxsz_z[0];
+ upp[2] = (vox_upp[2]-htcp_desc->lower[2])/htcp_desc->vxsz_z[0];
+ igrid_low[2] = (size_t)low[2];
+ igrid_upp[2] = (size_t)upp[2] - (modf(upp[2], &ipart)==0);
+ ASSERT(igrid_low[2] <= igrid_upp[2]);
+
+ /* Prepare the iteration overt the grid voxels overlapped by the SVX voxel */
+ x_h2o_range[0] = DBL_MAX;
+ x_h2o_range[1] =-DBL_MAX;
+
+ FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
+ FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
+ FOR_EACH(ivox[2], igrid_low[2], igrid_upp[2]+1) {
+
+ /* Compute the xH2O for the current LES voxel */
+ x_h2o = cloud_water_vapor_molar_fraction(htcp_desc, ivox);
+
+ /* Update the xH2O range of the SVX voxel */
+ x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
+ x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
+ #ifndef NDEBUG
+ {
+ double tmp_low[3], tmp_upp[3];
+ htcp_desc_get_voxel_aabb
+ (&ctx->sky->htcp_desc, ivox[0], ivox[1], ivox[2], tmp_low, tmp_upp);
+ ASSERT(aabb_intersect(tmp_low, tmp_upp, vox_low, vox_upp));
+ }
+ #endif
+ }
+ }
}
-#endif
} else { /* A SVX voxel can be overlapped by 2 LES voxels */
double pos_z;
size_t ilut_low, ilut_upp;
- size_t ivox_z_prev;
size_t ilut;
+ size_t ivox_z_prev;
ASSERT(xyz[2] < darray_split_size_get(&ctx->sky->svx2htcp_z));
/* Compute the *inclusive* bounds of the indices of the LUT cells
* overlapped by the SVX voxel */
- ilut_low = (size_t)
- ((vox_low[2] - ctx->sky->htcp_desc.lower[2]) / ctx->sky->lut_cell_sz);
- ilut_upp = (size_t)
- ((vox_upp[2] - ctx->sky->htcp_desc.lower[2]) / ctx->sky->lut_cell_sz);
+ ilut_low = (size_t)((vox_low[2] - htcp_desc->lower[2]) / ctx->sky->lut_cell_sz);
+ ilut_upp = (size_t)((vox_upp[2] - htcp_desc->lower[2]) / ctx->sky->lut_cell_sz);
ASSERT(ilut_low <= ilut_upp);
/* Prepare the iteration over the LES voxels overlapped by the SVX voxel */
- ivox_z_prev = SIZE_MAX;
x_h2o_range[0] = DBL_MAX;
x_h2o_range[1] =-DBL_MAX;
+ ivox_z_prev = SIZE_MAX;
pos_z = vox_low[2];
ASSERT(pos_z >= (double)ilut_low * ctx->sky->lut_cell_sz);
ASSERT(pos_z <= (double)ilut_upp * ctx->sky->lut_cell_sz);
@@ -737,15 +846,11 @@ cloud_vox_get_gas
/* Iterate over the LUT cells overlapped by the voxel */
FOR_EACH(ilut, ilut_low, ilut_upp+1) {
const struct split* split = darray_split_cdata_get(&ctx->sky->svx2htcp_z)+ilut;
- size_t ivox[3];
- double x_h2o;
ASSERT(ilut < darray_split_size_get(&ctx->sky->svx2htcp_z));
- ivox[0] = xyz[0];
- ivox[1] = xyz[1];
ivox[2] = pos_z <= split->height ? split->index : split->index + 1;
if(ivox[2] >= ctx->sky->htcp_desc.spatial_definition[2]
- && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
+ && eq_eps(pos_z, split->height, 1.e-6)) { /* Handle numerical inaccuracy */
ivox[2] = split->index;
}
@@ -760,12 +865,17 @@ cloud_vox_get_gas
if(ivox[2] == ivox_z_prev) continue;
ivox_z_prev = ivox[2];
- /* Compute the xH2O for the current LES voxel */
- x_h2o = cloud_water_vapor_molar_fraction(&ctx->sky->htcp_desc, ivox);
+ FOR_EACH(ivox[0], igrid_low[0], igrid_upp[0]+1) {
+ FOR_EACH(ivox[1], igrid_low[1], igrid_upp[1]+1) {
+
+ /* Compute the xH2O for the current LES voxel */
+ x_h2o = cloud_water_vapor_molar_fraction(&ctx->sky->htcp_desc, ivox);
- /* Update the xH2O range of the SVX voxel */
- x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
- x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
+ /* Update the xH2O range of the SVX voxel */
+ x_h2o_range[0] = MMIN(x_h2o, x_h2o_range[0]);
+ x_h2o_range[1] = MMAX(x_h2o, x_h2o_range[1]);
+ }
+ }
}
}
@@ -1129,6 +1239,7 @@ setup_clouds
const int force_cache_update)
{
struct darray_specdata specdata;
+ const size_t* raw_def;
size_t nvoxs[3];
double vxsz[3];
double low[3];
@@ -1152,21 +1263,22 @@ setup_clouds
SPLIT3(sky->htcp_desc.lower), SPLIT3(sky->htcp_desc.upper));
/* Define the number of voxels */
- nvoxs[0] = sky->htcp_desc.spatial_definition[0];
- nvoxs[1] = sky->htcp_desc.spatial_definition[1];
- nvoxs[2] = sky->htcp_desc.spatial_definition[2];
+ raw_def = sky->htcp_desc.spatial_definition;
+ nvoxs[0] = MMIN(raw_def[0], sky->htrdr->grid_max_definition[0]);
+ nvoxs[1] = MMIN(raw_def[1], sky->htrdr->grid_max_definition[1]);
+ nvoxs[2] = MMIN(raw_def[2], sky->htrdr->grid_max_definition[2]);
/* Define the octree AABB excepted for the Z dimension */
low[0] = sky->htcp_desc.lower[0];
low[1] = sky->htcp_desc.lower[1];
low[2] = sky->htcp_desc.lower[2];
- upp[0] = low[0] + (double)nvoxs[0] * sky->htcp_desc.vxsz_x;
- upp[1] = low[1] + (double)nvoxs[1] * sky->htcp_desc.vxsz_y;
+ upp[0] = low[0] + (double)raw_def[0] * sky->htcp_desc.vxsz_x;
+ upp[1] = low[1] + (double)raw_def[1] * sky->htcp_desc.vxsz_y;
if(!sky->htcp_desc.irregular_z) {
/* Regular voxel size along the Z dimension: compute its upper boundary as
* the others dimensions */
- upp[2] = low[2] + (double)nvoxs[2] * sky->htcp_desc.vxsz_z[0];
+ upp[2] = low[2] + (double)raw_def[2] * sky->htcp_desc.vxsz_z[0];
} else { /* Irregular voxel size along Z */
double len_z;
size_t nsplits;
@@ -1191,7 +1303,7 @@ setup_clouds
}
/* Setup the svx2htcp LUT. Each LUT entry stores the index of the current Z
* voxel in the htcp file that overlaps the entry lower bound as well as the
- * lower bound in Z of the the next htcp voxel. */
+ * lower bound in Z of the next htcp voxel. */
iz2 = 0;
upp[2] = low[2] + sky->htcp_desc.vxsz_z[iz2];
FOR_EACH(iz, 0, nsplits) {
@@ -1210,9 +1322,9 @@ setup_clouds
vxsz[0] = sky->htcp_desc.upper[0] - sky->htcp_desc.lower[0];
vxsz[1] = sky->htcp_desc.upper[1] - sky->htcp_desc.lower[1];
vxsz[2] = sky->htcp_desc.upper[2] - sky->htcp_desc.lower[2];
- vxsz[0] = vxsz[0] / (double)sky->htcp_desc.spatial_definition[0];
- vxsz[1] = vxsz[1] / (double)sky->htcp_desc.spatial_definition[1];
- vxsz[2] = vxsz[2] / (double)sky->htcp_desc.spatial_definition[2];
+ vxsz[0] = vxsz[0] / (double)nvoxs[0];
+ vxsz[1] = vxsz[1] / (double)nvoxs[1];
+ vxsz[2] = vxsz[2] / (double)nvoxs[2];
/* Create as many cloud data structure than considered SW spectral bands */
nbands = htrdr_sky_get_sw_spectral_bands_count(sky);
diff --git a/src/htrdr_solve.h b/src/htrdr_solve.h
@@ -51,7 +51,7 @@ htrdr_draw_radiance_sw
const size_t width, /* Image width */
const size_t height, /* Image height */
const size_t spp, /* #samples per pixel, i.e. #realisations */
- /* Buffer of struct htrdr_accum[3]. May be NULL on on non master processes */
+ /* Buffer of struct htrdr_accum[4]. May be NULL on non master processes */
struct htrdr_buffer* buf);
extern LOCAL_SYM int
@@ -62,4 +62,27 @@ htrdr_ground_filter
void* ray_data,
void* filter_data);
+static FINLINE void
+htrdr_accum_get_estimation
+ (const struct htrdr_accum* acc,
+ double* expected_value,
+ double* std_err)
+{
+ ASSERT(acc && expected_value && std_err);
+
+ if(!acc->nweights) {
+ *expected_value = 0;
+ *std_err = 0;
+ } else {
+ const double N = (double)acc->nweights;
+ double E, V, SE;
+ E = acc->sum_weights / N;
+ V = MMAX(acc->sum_weights_sqr / N - E*E, 0);
+ SE = sqrt(V/N);
+
+ *expected_value = E;
+ *std_err = SE;
+ }
+}
+
#endif /* HTRDR_SOLVE_H */
