commit 957bbfc9dd46d0536a99081177c8832941ac361c
parent 9fb45f8fb33244c05c12f37a110a09249b790b48
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 17 Apr 2019 11:18:39 +0200
Merge branch 'feature_time_per_realisation' into develop
Diffstat:
5 files changed, 120 insertions(+), 56 deletions(-)
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/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");
@@ -527,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;
}
@@ -629,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_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_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_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 */
