commit 9463de1105aeda9fcd0546bf3ea28d1317c6bba0
parent d690a610a8c644778cb2fe13ccaf89d6ce02f70d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 28 May 2019 15:57:36 +0200
Merge branch 'feature_improve_camera_solver' into develop
Diffstat:
17 files changed, 608 insertions(+), 364 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -16,9 +16,6 @@
cmake_minimum_required(VERSION 3.0)
project(stardis C)
enable_testing()
-if(POLICY CMP0054)
- cmake_policy(SET CMP0054 OLD) # Disable OpenMP CMake warning
-endif()
include(CMakeDependentOption)
@@ -65,11 +62,11 @@ set(VERSION_PATCH 1)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SDIS_FILES_SRC
- sdis_accum_buffer.c
sdis_camera.c
sdis_data.c
sdis_device.c
sdis_estimator.c
+ sdis_estimator_buffer.c
sdis_green.c
sdis_heat_path.c
sdis_interface.c
@@ -171,7 +168,6 @@ if(NOT NO_TEST)
register_test(${_name} ${_name})
endfunction()
- new_test(test_sdis_accum_buffer)
new_test(test_sdis_camera)
new_test(test_sdis_conducto_radiative)
new_test(test_sdis_conducto_radiative_2d)
diff --git a/src/sdis.h b/src/sdis.h
@@ -56,11 +56,11 @@ struct senc_descriptor;
* get or release a reference on the data, i.e. they increment or decrement the
* reference counter, respectively. When this counter reaches 0, the object is
* silently destroyed and cannot be used anymore. */
-struct sdis_accum_buffer;
struct sdis_camera;
struct sdis_data;
struct sdis_device;
struct sdis_estimator;
+struct sdis_estimator_buffer;
struct sdis_green_function;
struct sdis_interface;
struct sdis_medium;
@@ -137,14 +137,6 @@ struct sdis_interface_fragment {
static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL =
SDIS_INTERFACE_FRAGMENT_NULL__;
-/* Monte-Carlo accumulator */
-struct sdis_accum {
- double sum_weights; /* Sum of Monte-Carlo weights */
- double sum_weights_sqr; /* Sum of Monte-Carlo square weights */
- size_t nweights; /* #accumulated weights */
- size_t nfailures; /* #failures */
-};
-
/* Monte-Carlo estimation */
struct sdis_mc {
double E; /* Expected value */
@@ -244,22 +236,6 @@ struct sdis_interface_shader {
static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL =
SDIS_INTERFACE_SHADER_NULL__;
-struct sdis_accum_buffer_layout {
- size_t width;
- size_t height;
-};
-#define SDIS_ACCUM_BUFFER_LAYOUT_NULL__ {0, 0}
-static const struct sdis_accum_buffer_layout SDIS_ACCUM_BUFFER_LAYOUT_NULL =
- SDIS_ACCUM_BUFFER_LAYOUT_NULL__;
-
-/* Functor use to write accumulations performed by sdis_solve_camera */
-typedef res_T
-(*sdis_write_accums_T)
- (void* context, /* User data */
- const size_t origin[2], /* Coordinates of the 1st accumulation in image plane */
- const size_t naccums[2], /* #accumulations in X and Y */
- const struct sdis_accum* accums); /* List of row ordered accumulations */
-
/* Vertex of heat path v*/
struct sdis_heat_vertex {
double P[3];
@@ -418,51 +394,47 @@ sdis_camera_look_at
const double up[3]);
/*******************************************************************************
- * A buffer of accumulations is a 2D array whose each cell stores an
- * Monte-Carlo accumulation, i.e. a sum of MC weights, the sum of their square
- * and the overall number of summed weights (see struct sdis_accum)
- ******************************************************************************/
+ * An estimator buffer is 2D array of estimators
+******************************************************************************/
SDIS_API res_T
-sdis_accum_buffer_create
- (struct sdis_device* dev,
- const size_t width, /* #cells in X */
- const size_t height, /* #cells in Y */
- struct sdis_accum_buffer** buf);
+sdis_estimator_buffer_ref_get
+ (struct sdis_estimator_buffer* buf);
SDIS_API res_T
-sdis_accum_buffer_ref_get
- (struct sdis_accum_buffer* buf);
+sdis_estimator_buffer_ref_put
+ (struct sdis_estimator_buffer* buf);
SDIS_API res_T
-sdis_accum_buffer_ref_put
- (struct sdis_accum_buffer* buf);
+sdis_estimator_buffer_get_definition
+ (const struct sdis_estimator_buffer* buf,
+ size_t definition[2]);
SDIS_API res_T
-sdis_accum_buffer_get_layout
- (const struct sdis_accum_buffer* buf,
- struct sdis_accum_buffer_layout* layout);
+sdis_estimator_buffer_at
+ (const struct sdis_estimator_buffer* buf,
+ const size_t x,
+ const size_t y,
+ const struct sdis_estimator** estimator);
-/* Get a read only pointer toward the memory space of the accum buffer */
SDIS_API res_T
-sdis_accum_buffer_map
- (const struct sdis_accum_buffer* buf,
- const struct sdis_accum** accums);
+sdis_estimator_buffer_get_realisation_count
+ (const struct sdis_estimator_buffer* buf,
+ size_t* nrealisations); /* Successful ones */
+
+SDIS_API res_T
+sdis_estimator_buffer_get_failure_count
+ (const struct sdis_estimator_buffer* buf,
+ size_t* nfailures);
SDIS_API res_T
-sdis_accum_buffer_unmap
- (const struct sdis_accum_buffer* buf);
+sdis_estimator_buffer_get_temperature
+ (const struct sdis_estimator_buffer* buf,
+ struct sdis_mc* temperature);
-/* Helper function that matches the `sdis_write_accums_T' functor type. One can
- * send this function directly to the sdis_solve_camera function, to fill the
- * accum buffer with the estimation of the radiative temperature that reaches
- * each pixel of an image whose definition matches the definition of the accum
- * buffer. */
SDIS_API res_T
-sdis_accum_buffer_write
- (void* buf, /* User data */
- const size_t origin[2], /* Coordinates of the 1st accum in image plane */
- const size_t naccum[2], /* #accum in X and Y */
- const struct sdis_accum* accums); /* List of row ordered accum */
+sdis_estimator_buffer_get_realisation_time
+ (const struct sdis_estimator_buffer* buf,
+ struct sdis_mc* time);
/*******************************************************************************
* A medium encapsulates the properties of either a fluid or a solid.
@@ -931,15 +903,15 @@ SDIS_API res_T
sdis_solve_camera
(struct sdis_scene* scn,
const struct sdis_camera* cam, /* Point of view */
- const double time, /* Observation time */
+ const double time_range[2], /* Observation time */
const double fp_to_meter, /* Scale from floating point units to meters */
const double ambient_radiative_temperature, /* In Kelvin */
const double reference_temperature, /* In Kelvin */
const size_t width, /* Image definition in in X */
const size_t height, /* Image definition in Y */
const size_t spp, /* #samples per pixel */
- sdis_write_accums_T writer,
- void* writer_data);
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
+ struct sdis_estimator_buffer** buf);
SDIS_API res_T
sdis_solve_medium
diff --git a/src/sdis_accum_buffer.c b/src/sdis_accum_buffer.c
@@ -1,160 +0,0 @@
-/* Copyright (C) 2016-2019 |Meso|Star> (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "sdis.h"
-#include "sdis_device_c.h"
-
-struct sdis_accum_buffer {
- struct sdis_accum* accums;
- size_t width;
- size_t height;
-
- ref_T ref;
- struct sdis_device* dev;
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static void
-accum_buffer_release(ref_T* ref)
-{
- struct sdis_accum_buffer* buf;
- struct sdis_device* dev;
- ASSERT(ref);
- buf = CONTAINER_OF(ref, struct sdis_accum_buffer, ref);
- dev = buf->dev;
- if(buf->accums) MEM_RM(dev->allocator, buf->accums);
- MEM_RM(dev->allocator, buf);
- SDIS(device_ref_put(dev));
-}
-
-/*******************************************************************************
- * Export functions
- ******************************************************************************/
-res_T
-sdis_accum_buffer_create
- (struct sdis_device* dev,
- const size_t width,
- const size_t height,
- struct sdis_accum_buffer** out_buf)
-{
- struct sdis_accum_buffer* buf = NULL;
- size_t sz;
- res_T res = RES_OK;
-
- if(!dev || !width || !height || !out_buf) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- buf = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_accum_buffer));
- if(!buf) {
- res = RES_MEM_ERR;
- goto error;
- }
- ref_init(&buf->ref);
- SDIS(device_ref_get(dev));
- buf->dev = dev;
- buf->width = width;
- buf->height = height;
-
- sz = buf->width * buf->height * sizeof(struct sdis_accum);
- buf->accums = MEM_ALLOC_ALIGNED(dev->allocator, sz, 16);
- if(!buf->accums) {
- log_err(dev, "%s: could not allocate a buffer accumulations of %lux%lu.\n",
- FUNC_NAME, (unsigned long)width, (unsigned long)height);
- res = RES_MEM_ERR;
- goto error;
- }
-
-exit:
- if(out_buf) *out_buf = buf;
- return res;
-error:
- if(buf) {
- SDIS(accum_buffer_ref_put(buf));
- buf = NULL;
- }
- goto exit;
-}
-
-res_T
-sdis_accum_buffer_ref_get(struct sdis_accum_buffer* buf)
-{
- if(!buf) return RES_BAD_ARG;
- ref_get(&buf->ref);
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_ref_put(struct sdis_accum_buffer* buf)
-{
- if(!buf) return RES_BAD_ARG;
- ref_put(&buf->ref, accum_buffer_release);
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_get_layout
- (const struct sdis_accum_buffer* buf,
- struct sdis_accum_buffer_layout* layout)
-{
- if(!buf || !layout) return RES_BAD_ARG;
- layout->width = buf->width;
- layout->height = buf->height;
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_map
- (const struct sdis_accum_buffer* buf, const struct sdis_accum** accums)
-{
- if(!buf || !accums) return RES_BAD_ARG;
- *accums = buf->accums;
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_unmap(const struct sdis_accum_buffer* buf)
-{
- if(!buf) return RES_BAD_ARG;
- /* Do nothing */
- return RES_OK;
-}
-
-res_T
-sdis_accum_buffer_write
- (void* buffer,
- const size_t origin[2],
- const size_t size[2],
- const struct sdis_accum* accums)
-{
- struct sdis_accum_buffer* buf = buffer;
- const struct sdis_accum* src_row = accums;
- size_t dst_iy;
-
- if(UNLIKELY(!buffer || !origin || !size || !accums)) return RES_BAD_ARG;
- if(UNLIKELY((origin[0] + size[0]) > buf->width)) return RES_BAD_ARG;
- if(UNLIKELY((origin[1] + size[1]) > buf->height)) return RES_BAD_ARG;
-
- FOR_EACH(dst_iy, origin[1], origin[1] + size[1]) {
- const size_t dst_irow = dst_iy*buf->width + origin[0];
- memcpy(buf->accums + dst_irow, src_row, size[0] * sizeof(struct sdis_accum));
- src_row += size[0];
- }
- return RES_OK;
-}
-
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -54,7 +54,6 @@ device_release(ref_T* ref)
ASSERT(flist_name_is_empty(&dev->media_names));
flist_name_release(&dev->interfaces_names);
flist_name_release(&dev->media_names);
- darray_tile_release(&dev->tiles);
MEM_RM(dev->allocator, dev);
}
@@ -98,14 +97,6 @@ sdis_device_create
ref_init(&dev->ref);
flist_name_init(allocator, &dev->interfaces_names);
flist_name_init(allocator, &dev->media_names);
- darray_tile_init(allocator, &dev->tiles);
-
- res = darray_tile_resize(&dev->tiles, dev->nthreads);
- if(res != RES_OK) {
- log_err(dev,
- "%s: could not allocate the per thread buffer of estimations.\n", FUNC_NAME);
- goto error;
- }
res = s2d_device_create(log, allocator, 0, &dev->s2d);
if(res != RES_OK) {
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -26,18 +26,6 @@ struct name { FITEM; };
#define FITEM_TYPE name
#include <rsys/free_list.h>
-#define DARRAY_NAME accum
-#define DARRAY_DATA struct sdis_accum
-#include <rsys/dynamic_array.h>
-
-#define DARRAY_NAME tile
-#define DARRAY_DATA struct darray_accum
-#define DARRAY_FUNCTOR_INIT darray_accum_init
-#define DARRAY_FUNCTOR_RELEASE darray_accum_release
-#define DARRAY_FUNCTOR_COPY darray_accum_copy
-#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_accum_copy_and_release
-#include <rsys/dynamic_array.h>
-
struct sdis_device {
struct logger* logger;
struct mem_allocator* allocator;
@@ -46,7 +34,6 @@ struct sdis_device {
struct flist_name interfaces_names;
struct flist_name media_names;
- struct darray_tile tiles;
struct s2d_device* s2d;
struct s3d_device* s3d;
diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c
@@ -82,12 +82,19 @@ sdis_estimator_get_failure_count
return RES_OK;
}
+#define SETUP_MC(Mc, Acc) { \
+ (Mc)->E = (Acc)->sum / (double)(Acc)->count; \
+ (Mc)->V = (Acc)->sum2 / (double)(Acc)->count - (Mc)->E*(Mc)->E; \
+ (Mc)->V = MMAX((Mc)->V, 0); \
+ (Mc)->SE = sqrt((Mc)->V / (double)(Acc)->count); \
+} (void)0
+
res_T
sdis_estimator_get_temperature
(const struct sdis_estimator* estimator, struct sdis_mc* mc)
{
if(!estimator || !mc) return RES_BAD_ARG;
- *mc = estimator->temperature;
+ SETUP_MC(mc, &estimator->temperature);
return RES_OK;
}
@@ -96,7 +103,7 @@ sdis_estimator_get_realisation_time
(const struct sdis_estimator* estimator, struct sdis_mc* mc)
{
if(!estimator || !mc) return RES_BAD_ARG;
- *mc = estimator->realisation_time;
+ SETUP_MC(mc, &estimator->realisation_time);
return RES_OK;
}
@@ -107,7 +114,7 @@ sdis_estimator_get_convective_flux
if(!estimator || !flux ||estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
ASSERT(estimator->fluxes);
- *flux = estimator->fluxes[FLUX_CONVECTIVE];
+ SETUP_MC(flux, &estimator->fluxes[FLUX_CONVECTIVE]);
return RES_OK;
}
@@ -118,7 +125,7 @@ sdis_estimator_get_radiative_flux
if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
ASSERT(estimator->fluxes);
- *flux = estimator->fluxes[FLUX_RADIATIVE];
+ SETUP_MC(flux, &estimator->fluxes[FLUX_RADIATIVE]);
return RES_OK;
}
@@ -129,10 +136,12 @@ sdis_estimator_get_total_flux
if(!estimator || !flux || estimator->type != SDIS_ESTIMATOR_FLUX)
return RES_BAD_ARG;
ASSERT(estimator->fluxes);
- *flux = estimator->fluxes[FLUX_TOTAL];
+ SETUP_MC(flux, &estimator->fluxes[FLUX_TOTAL]);
return RES_OK;
}
+#undef SETUP_MC
+
res_T
sdis_estimator_get_paths_count
(const struct sdis_estimator* estimator, size_t* npaths)
diff --git a/src/sdis_estimator_buffer.c b/src/sdis_estimator_buffer.c
@@ -0,0 +1,244 @@
+/* Copyright (C) 2016-2019 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "sdis_device_c.h"
+#include "sdis_estimator_c.h"
+#include "sdis_estimator_buffer_c.h"
+
+struct sdis_estimator_buffer {
+ struct sdis_estimator** estimators; /* Row major per pixe lestimators */
+ size_t width;
+ size_t height;
+
+ struct accum temperature;
+ struct accum realisation_time;
+ size_t nrealisations; /* #successes */
+ size_t nfailures;
+
+ ref_T ref;
+ struct sdis_device* dev;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+estimator_buffer_release(ref_T* ref)
+{
+ struct sdis_estimator_buffer* buf;
+ struct sdis_device* dev;
+ size_t i;
+ ASSERT(ref);
+
+ buf = CONTAINER_OF(ref, struct sdis_estimator_buffer, ref);
+ dev = buf->dev;
+
+ if(buf->estimators) {
+ FOR_EACH(i, 0, buf->width*buf->height) {
+ if(buf->estimators[i]) SDIS(estimator_ref_put(buf->estimators[i]));
+ }
+ MEM_RM(dev->allocator, buf->estimators);
+ }
+
+ MEM_RM(dev->allocator, buf);
+ SDIS(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sdis_estimator_buffer_ref_get(struct sdis_estimator_buffer* buf)
+{
+ if(!buf) return RES_BAD_ARG;
+ ref_get(&buf->ref);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_ref_put(struct sdis_estimator_buffer* buf)
+{
+ if(!buf) return RES_BAD_ARG;
+ ref_put(&buf->ref, estimator_buffer_release);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_get_definition
+ (const struct sdis_estimator_buffer* buf, size_t definition[2])
+{
+ if(!buf || !definition) return RES_BAD_ARG;
+ definition[0] = buf->width;
+ definition[1] = buf->height;
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_at
+ (const struct sdis_estimator_buffer* buf,
+ const size_t x,
+ const size_t y,
+ const struct sdis_estimator** estimator)
+{
+ if(!buf || x >= buf->width || y >= buf->height || !estimator)
+ return RES_BAD_ARG;
+ *estimator = estimator_buffer_grab(buf, x, y);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_get_realisation_count
+ (const struct sdis_estimator_buffer* buf, size_t* nrealisations)
+{
+ if(!buf || !nrealisations) return RES_BAD_ARG;
+ *nrealisations = buf->nrealisations;
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_get_failure_count
+ (const struct sdis_estimator_buffer* buf, size_t* nfailures)
+{
+ if(!buf || !nfailures) return RES_BAD_ARG;
+ *nfailures = buf->nfailures;
+ return RES_OK;
+}
+
+#define SETUP_MC(Mc, Acc) { \
+ (Mc)->E = (Acc)->sum / (double)(Acc)->count; \
+ (Mc)->V = (Acc)->sum2 / (double)(Acc)->count - (Mc)->E*(Mc)->E; \
+ (Mc)->V = MMAX((Mc)->V, 0); \
+ (Mc)->SE = sqrt((Mc)->V / (double)(Acc)->count); \
+} (void)0
+
+res_T
+sdis_estimator_buffer_get_temperature
+ (const struct sdis_estimator_buffer* buf, struct sdis_mc* mc)
+{
+ if(!buf || !mc) return RES_BAD_ARG;
+ SETUP_MC(mc, &buf->temperature);
+ return RES_OK;
+}
+
+res_T
+sdis_estimator_buffer_get_realisation_time
+(const struct sdis_estimator_buffer* buf, struct sdis_mc* mc)
+{
+ if(!buf || !mc) return RES_BAD_ARG;
+ SETUP_MC(mc, &buf->realisation_time);
+ return RES_OK;
+}
+
+#undef SETUP_MC
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+estimator_buffer_create
+ (struct sdis_device* dev,
+ const size_t width,
+ const size_t height,
+ struct sdis_estimator_buffer** out_buf)
+{
+ struct sdis_estimator_buffer* buf = NULL;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!dev || !width || !height || !out_buf) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ buf = MEM_CALLOC(dev->allocator, 1, sizeof(*buf));
+ if(!buf) { res = RES_MEM_ERR; goto error; }
+
+ ref_init(&buf->ref);
+ SDIS(device_ref_get(dev));
+ buf->dev = dev;
+ buf->width = width;
+ buf->height = height;
+
+ buf->estimators = MEM_CALLOC
+ (dev->allocator, width*height, sizeof(*buf->estimators));
+ if(!buf->estimators) {
+ log_err(dev, "%s: could not allocate a buffer of estimators of %lux%lu.\n",
+ FUNC_NAME, (unsigned long)width, (unsigned long)height);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ FOR_EACH(i, 0, width*height) {
+ res = estimator_create(dev, SDIS_ESTIMATOR_TEMPERATURE, buf->estimators+i);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ if(out_buf) *out_buf = buf;
+ return res;
+error:
+ if(buf) {
+ SDIS(estimator_buffer_ref_put(buf));
+ buf = NULL;
+ }
+ goto exit;
+}
+
+struct sdis_estimator*
+estimator_buffer_grab
+ (const struct sdis_estimator_buffer* buf,
+ const size_t x,
+ const size_t y)
+{
+ ASSERT(x < buf->width && y < buf->height);
+ return buf->estimators[y*buf->width + x];
+}
+
+void
+estimator_buffer_setup_realisations_count
+ (struct sdis_estimator_buffer* buf,
+ const size_t nrealisations,
+ const size_t nsuccesses)
+{
+ ASSERT(buf && nrealisations && nsuccesses && nsuccesses<=nrealisations);
+ buf->nrealisations = nsuccesses;
+ buf->nfailures = nrealisations - nsuccesses;
+}
+
+void
+estimator_buffer_setup_temperature
+ (struct sdis_estimator_buffer* buf,
+ const double sum,
+ const double sum2)
+{
+ ASSERT(buf && buf->nrealisations);
+ buf->temperature.sum = sum;
+ buf->temperature.sum2 = sum2;
+ buf->temperature.count = buf->nrealisations;
+}
+
+void
+estimator_buffer_setup_realisation_time
+ (struct sdis_estimator_buffer* buf,
+ const double sum,
+ const double sum2)
+{
+ ASSERT(buf && buf->nrealisations);
+ buf->realisation_time.sum = sum;
+ buf->realisation_time.sum2 = sum2;
+ buf->realisation_time.count = buf->nrealisations;
+}
+
diff --git a/src/sdis_estimator_buffer_c.h b/src/sdis_estimator_buffer_c.h
@@ -0,0 +1,58 @@
+/* Copyright (C) 2016-2019 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SDIS_ESTIMATOR_BUFFER_C_H
+#define SDIS_ESTIMATOR_BUFFER_C_H
+
+#include <rsys/rsys.h>
+
+/* Forward declaration */
+struct sdis_device;
+struct sdis_estimator;
+struct sdis_estimator_buffer;
+
+extern LOCAL_SYM res_T
+estimator_buffer_create
+ (struct sdis_device* dev,
+ const size_t width,
+ const size_t height,
+ struct sdis_estimator_buffer** buf);
+
+extern LOCAL_SYM struct sdis_estimator*
+estimator_buffer_grab
+ (const struct sdis_estimator_buffer* buf,
+ const size_t x,
+ const size_t y);
+
+extern LOCAL_SYM void
+estimator_buffer_setup_realisations_count
+ (struct sdis_estimator_buffer* buf,
+ const size_t nrealisations,
+ const size_t nsuccesses);
+
+extern LOCAL_SYM void
+estimator_buffer_setup_temperature
+ (struct sdis_estimator_buffer* buf,
+ const double sum,
+ const double sum2);
+
+extern LOCAL_SYM void
+estimator_buffer_setup_realisation_time
+ (struct sdis_estimator_buffer* buf,
+ const double sum,
+ const double sum2);
+
+#endif /* SDIS_ESTIMATOR_BUFFER_C_H */
+
diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h
@@ -17,6 +17,7 @@
#define SDIS_ESTIMATOR_C_H
#include "sdis_heat_path.h"
+#include "sdis_misc.h"
#include <rsys/math.h>
#include <rsys/ref_count.h>
@@ -34,10 +35,10 @@ enum flux_name {
};
struct sdis_estimator {
- struct sdis_mc temperature;
- struct sdis_mc realisation_time;
- struct sdis_mc fluxes[FLUX_NAMES_COUNT__];
- size_t nrealisations;
+ struct accum temperature;
+ struct accum realisation_time;
+ struct accum fluxes[FLUX_NAMES_COUNT__];
+ size_t nrealisations; /* #successes */
size_t nfailures;
struct mutex* mutex;
@@ -77,13 +78,6 @@ estimator_setup_realisations_count
estimator->nfailures = nrealisations - nsuccesses;
}
-#define SETUP_MC(McName, Sum, Sum2, Count) { \
- (McName).E = (Sum) / (double)(Count); \
- (McName).V = (Sum2) / (double)(Count) - (McName).E*(McName).E; \
- (McName).V = MMAX((McName).V, 0); \
- (McName).SE = sqrt((McName).V / (double)(Count)); \
-} (void)0
-
static INLINE void
estimator_setup_temperature
(struct sdis_estimator* estim,
@@ -91,7 +85,9 @@ estimator_setup_temperature
const double sum2)
{
ASSERT(estim && estim->nrealisations);
- SETUP_MC(estim->temperature, sum, sum2, estim->nrealisations);
+ estim->temperature.sum = sum;
+ estim->temperature.sum2 = sum2;
+ estim->temperature.count = estim->nrealisations;
}
static INLINE void
@@ -101,7 +97,9 @@ estimator_setup_realisation_time
const double sum2)
{
ASSERT(estim && estim->nrealisations);
- SETUP_MC(estim->realisation_time, sum, sum2, estim->nrealisations);
+ estim->realisation_time.sum = sum;
+ estim->realisation_time.sum2 = sum2;
+ estim->realisation_time.count = estim->nrealisations;
}
static INLINE void
@@ -112,10 +110,10 @@ estimator_setup_flux
const double sum2)
{
ASSERT(estim && (unsigned)name < FLUX_NAMES_COUNT__ && estim->nrealisations);
- SETUP_MC(estim->fluxes[name], sum, sum2, estim->nrealisations);
+ estim->fluxes[name].sum = sum;
+ estim->fluxes[name].sum2 = sum2;
+ estim->fluxes[name].count = estim->nrealisations;
}
-#undef SETUP_MC
-
#endif /* SDIS_PROBE_ESTIMATOR_C_H */
diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c
@@ -32,6 +32,7 @@ ray_realisation_3d
const double fp_to_meter,
const double Tarad,
const double Tref,
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -50,9 +51,14 @@ ray_realisation_3d
ctx.Tarad = Tarad;
ctx.Tref3 = Tref*Tref*Tref;
+ ctx.heat_path = heat_path;
f3_set_d3(dir, direction);
+ /* Register the starting position against the heat path */
+ res = register_heat_vertex(heat_path, &rwalk.vtx, 0, SDIS_HEAT_VERTEX_RADIATIVE);
+ if(res != RES_OK) goto error;
+
res = trace_radiative_path_3d(scn, dir, fp_to_meter, &ctx, &rwalk, rng, &T);
if(res != RES_OK) goto error;
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -46,8 +46,8 @@ probe_realisation_2d
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -61,8 +61,8 @@ probe_realisation_3d
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
/*******************************************************************************
@@ -79,8 +79,8 @@ boundary_realisation_2d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -94,8 +94,8 @@ boundary_realisation_3d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
extern LOCAL_SYM res_T
@@ -140,6 +140,7 @@ ray_realisation_3d
const double fp_to_meter,
const double ambient_radiative_temperature,
const double reference_temperature,
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight);
#endif /* SDIS_REALISATION_H */
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -121,8 +121,8 @@ XD(probe_realisation)
const double fp_to_meter,/* Scale factor from floating point unit to meter */
const double ambient_radiative_temperature,
const double reference_temperature,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
@@ -213,8 +213,8 @@ XD(boundary_realisation)
const double fp_to_meter,
const double Tarad,
const double Tref,
- struct green_path_handle* green_path,
- struct sdis_heat_path* heat_path,
+ struct green_path_handle* green_path, /* May be NULL */
+ struct sdis_heat_path* heat_path, /* May be NULL */
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -17,6 +17,7 @@
#include "sdis_camera.h"
#include "sdis_device_c.h"
#include "sdis_estimator_c.h"
+#include "sdis_estimator_buffer_c.h"
#include "sdis_interface_c.h"
#include <star/ssp.h>
@@ -66,28 +67,43 @@ solve_pixel
struct ssp_rng* rng,
struct sdis_medium* mdm,
const struct sdis_camera* cam,
- const double time, /* Observation time */
+ const double time_range[2], /* Observation time */
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
const size_t ipix[2], /* Pixel coordinate in the image plane */
const size_t nrealisations,
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct sdis_accum* accum)
+ struct sdis_estimator* estimator)
{
- double sum_weights = 0;
- double sum_weights_sqr = 0;
- size_t N = 0; /* #realisations that do not fail */
+ struct accum acc_temp = ACCUM_NULL;
+ struct accum acc_time = ACCUM_NULL;
size_t irealisation;
res_T res = RES_OK;
ASSERT(scn && mdm && rng && cam && ipix && nrealisations && Tref >= 0);
ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+ ASSERT(estimator && time_range);
FOR_EACH(irealisation, 0, nrealisations) {
+ struct time t0, t1;
double samp[2]; /* Pixel sample */
double ray_pos[3];
double ray_dir[3];
double w = 0;
+ struct sdis_heat_path* pheat_path = NULL;
+ struct sdis_heat_path heat_path;
+ double time;
+ res_T res_simul = RES_OK;
+
+ /* Begin time registration */
+ time_current(&t0);
+
+ time = sample_time(rng, time_range);
+ if(register_paths) {
+ heat_path_init(scn->dev->allocator, &heat_path);
+ pheat_path = &heat_path;
+ }
/* Generate a sample into the pixel to estimate */
samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
@@ -98,23 +114,45 @@ solve_pixel
camera_ray(cam, samp, ray_pos, ray_dir);
/* Launch the realisation */
- res = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
- time, fp_to_meter, Tarad, Tref, &w);
- if(res == RES_OK) {
- sum_weights += w;
- sum_weights_sqr += w*w;
- ++N;
- } else if(res == RES_BAD_OP) {
- res = RES_OK;
- } else {
+ res_simul = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
+ time, fp_to_meter, Tarad, Tref, pheat_path, &w);
+
+ /* Handle fatal error */
+ if(res_simul != RES_OK && res_simul != RES_BAD_OP) {
+ res = res_simul;
goto error;
}
+
+ if(pheat_path) {
+ pheat_path->status = res_simul == RES_OK
+ ? SDIS_HEAT_PATH_SUCCEED
+ : SDIS_HEAT_PATH_FAILED;
+
+ /* Check if the path must be saved regarding the register_paths mask */
+ if(!(register_paths & (int)pheat_path->status)) {
+ heat_path_release(pheat_path);
+ } else { /* Register the sampled path */
+ res = estimator_add_and_release_heat_path(estimator, pheat_path);
+ if(res != RES_OK) goto error;
+ }
+ }
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+
+ if(res_simul == RES_OK) {
+ /* Update global accumulators */
+ const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+ acc_temp.sum += w; acc_temp.sum2 += w*w; ++acc_temp.count;
+ acc_time.sum += usec; acc_time.sum2 += usec*usec; ++acc_time.count;
+ }
}
- accum->sum_weights = sum_weights;
- accum->sum_weights_sqr = sum_weights_sqr;
- accum->nweights = N;
- accum->nfailures = nrealisations - N;
+ /* Setup the pixel estimator */
+ ASSERT(acc_temp.count == acc_time.count);
+ estimator_setup_realisations_count(estimator, nrealisations, acc_temp.count);
+ estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2);
+ estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2);
exit:
return res;
@@ -128,22 +166,23 @@ solve_tile
struct ssp_rng* rng,
struct sdis_medium* mdm,
const struct sdis_camera* cam,
- const double time,
+ const double time_range[2],
const double fp_to_meter,
const double Tarad,
const double Tref,
const size_t origin[2], /* Tile origin in image plane */
const size_t size[2], /* #pixels in X and Y */
const size_t spp, /* #samples per pixel */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct sdis_accum* accums)
+ struct sdis_estimator_buffer* buf)
{
size_t mcode; /* Morton code of the tile pixel */
size_t npixels;
res_T res = RES_OK;
- ASSERT(scn && rng && mdm && cam && spp && origin && accums && Tref >= 0);
- ASSERT(size &&size[0] && size[1]);
- ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+ ASSERT(scn && rng && mdm && cam && spp && origin && Tref >= 0);
+ ASSERT(size &&size[0] && size[1] && buf);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0 && time_range);
/* Adjust the #pixels to process them wrt a morton order */
npixels = round_up_pow2(MMAX(size[0], size[1]));
@@ -151,19 +190,21 @@ solve_tile
FOR_EACH(mcode, 0, npixels) {
size_t ipix[2];
- struct sdis_accum* accum;
+ struct sdis_estimator* estimator;
ipix[0] = morton2D_decode((uint32_t)(mcode>>0));
if(ipix[0] >= size[0]) continue;
ipix[1] = morton2D_decode((uint32_t)(mcode>>1));
if(ipix[1] >= size[1]) continue;
- accum = accums + ipix[1]*size[0] + ipix[0];
ipix[0] = ipix[0] + origin[0];
ipix[1] = ipix[1] + origin[1];
- res = solve_pixel(scn, rng, mdm, cam, time, fp_to_meter, Tarad, Tref, ipix,
- spp, pix_sz, accum);
+ /* Fetch the pixel estimator */
+ estimator = estimator_buffer_grab(buf, ipix[0], ipix[1]);
+
+ res = solve_pixel(scn, rng, mdm, cam, time_range, fp_to_meter, Tarad, Tref,
+ ipix, spp, register_paths, pix_sz, estimator);
if(res != RES_OK) goto error;
}
@@ -360,39 +401,48 @@ res_T
sdis_solve_camera
(struct sdis_scene* scn,
const struct sdis_camera* cam,
- const double time,
+ const double time_range[2],
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
const size_t width, /* #pixels in X */
const size_t height, /* #pixels in Y */
const size_t spp, /* #samples per pixel */
- sdis_write_accums_T writer,
- void* writer_data)
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
+ struct sdis_estimator_buffer** out_buf)
{
#define TILE_SIZE 32 /* definition in X & Y of a tile */
STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
+ struct sdis_estimator_buffer* buf = NULL;
struct sdis_medium* medium = NULL;
- struct darray_accum* tiles = NULL;
+ struct accum acc_temp = ACCUM_NULL;
+ struct accum acc_time = ACCUM_NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** rngs = NULL;
size_t ntiles_x, ntiles_y, ntiles;
double pix_sz[2]; /* Size of a pixel in the normalized image plane */
int64_t mcode; /* Morton code of a tile */
+ size_t nrealisations;
+ size_t nsuccesses;
+ size_t ix, iy;
size_t i;
ATOMIC res = RES_OK;
- if(!scn || !cam || time < 0 || fp_to_meter <= 0 || Tref < 0 || !width
- || !height || !spp || !writer) {
+ if(!scn || !cam || fp_to_meter <= 0 || Tref < 0 || !width || !height || !spp
+ || !out_buf) {
res = RES_BAD_ARG;
goto error;
}
-
if(scene_is_2d(scn)) {
log_err(scn->dev, "%s: 2D scene are not supported.\n", FUNC_NAME);
goto error;
}
+ if(!time_range || time_range[0] < 0 || time_range[1] < time_range[0]
+ || (time_range[1] > DBL_MAX && time_range[0] != time_range[1])) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
/* Retrieve the medium in which the submitted position lies */
res = scene_get_medium(scn, cam->position, NULL, &medium);
@@ -422,15 +472,6 @@ sdis_solve_camera
if(res != RES_OK) goto error;
}
- /* Allocate per thread buffer of accumulations */
- tiles = darray_tile_data_get(&scn->dev->tiles);
- ASSERT(darray_tile_size_get(&scn->dev->tiles) == scn->dev->nthreads);
- FOR_EACH(i, 0, scn->dev->nthreads) {
- const size_t naccums = TILE_SIZE * TILE_SIZE;
- res = darray_accum_resize(tiles+i, naccums);
- if(res != RES_OK) goto error;
- }
-
ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
ntiles_y = (height + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y));
@@ -439,13 +480,16 @@ sdis_solve_camera
pix_sz[0] = 1.0 / (double)width;
pix_sz[1] = 1.0 / (double)height;
+ /* Create the global estimator */
+ res = estimator_buffer_create(scn->dev, width, height, &buf);
+ if(res != RES_OK) goto error;
+
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static, 1/*chunk size*/)
for(mcode = 0; mcode < (int64_t)ntiles; ++mcode) {
size_t tile_org[2] = {0, 0};
size_t tile_sz[2] = {0, 0};
const int ithread = omp_get_thread_num();
- struct sdis_accum* accums = NULL;
struct ssp_rng* rng = rngs[ithread];
res_T res_local = RES_OK;
@@ -462,25 +506,40 @@ sdis_solve_camera
tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]);
tile_sz[1] = MMIN(TILE_SIZE, height - tile_org[1]);
- /* Fetch the accumulations buffer */
- accums = darray_accum_data_get(tiles+ithread);
-
/* Draw the tile */
- res_local = solve_tile(scn, rng, medium, cam, time, fp_to_meter, Tarad,
- Tref, tile_org, tile_sz, spp, pix_sz, accums);
+ res_local = solve_tile(scn, rng, medium, cam, time_range, fp_to_meter,
+ Tarad, Tref, tile_org, tile_sz, spp, register_paths, pix_sz, buf);
if(res_local != RES_OK) {
ATOMIC_SET(&res, res_local);
continue;
}
+ }
- /* Write the accumulations */
- res_local = writer(writer_data, tile_org, tile_sz, accums);
- if(res_local != RES_OK) {
- ATOMIC_SET(&res, res_local);
- continue;
+ /* Setup the accumulators of the whole estimator buffer */
+ acc_temp = ACCUM_NULL;
+ acc_time = ACCUM_NULL;
+ nsuccesses = 0;
+ FOR_EACH(iy, 0, height) {
+ FOR_EACH(ix, 0, width) {
+ const struct sdis_estimator* estimator;
+ SDIS(estimator_buffer_at(buf, ix, iy, &estimator));
+ acc_temp.sum += estimator->temperature.sum;
+ acc_temp.sum2 += estimator->temperature.sum2;
+ acc_temp.count += estimator->temperature.count;
+ acc_time.sum += estimator->realisation_time.sum;
+ acc_time.sum2 += estimator->realisation_time.sum2;
+ acc_time.count += estimator->realisation_time.count;
+ nsuccesses += estimator->nrealisations;
}
}
+ nrealisations = width*height*spp;
+ ASSERT(acc_temp.count == acc_time.count);
+ ASSERT(acc_temp.count == nsuccesses);
+ estimator_buffer_setup_realisations_count(buf, nrealisations, nsuccesses);
+ estimator_buffer_setup_temperature(buf, acc_temp.sum, acc_temp.sum2);
+ estimator_buffer_setup_realisation_time(buf, acc_time.sum, acc_time.sum2);
+
exit:
if(rngs) {
FOR_EACH(i, 0, scn->dev->nthreads) {
@@ -489,6 +548,7 @@ exit:
MEM_RM(scn->dev->allocator, rngs);
}
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_buf) *out_buf = buf;
return (res_T)res;
error:
goto exit;
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -90,7 +90,7 @@ XD(solve_probe)
if(res != RES_OK) goto error;
}
- /* Create the per thread accumulator */
+ /* Create the per thread accumulators */
acc_temps = MEM_CALLOC
(scn->dev->allocator, scn->dev->nthreads, sizeof(*acc_temps));
if(!acc_temps) { res = RES_MEM_ERR; goto error; }
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -433,39 +433,59 @@ geometry_add_shape
}
static void
-dump_image(const struct sdis_accum_buffer* buf)
+dump_image(const struct sdis_estimator_buffer* buf)
{
- struct sdis_accum_buffer_layout layout = SDIS_ACCUM_BUFFER_LAYOUT_NULL;
struct image img;
double* temps = NULL;
- const struct sdis_accum* accums = NULL;
double Tmax = -DBL_MAX;
double Tmin = DBL_MAX;
double norm;
+ size_t definition[2];
size_t i, ix, iy;
CHK(buf != NULL);
- OK(sdis_accum_buffer_get_layout(buf, &layout));
+ OK(sdis_estimator_buffer_get_definition(buf, definition));
- temps = mem_alloc(layout.width*layout.height*sizeof(double));
+ temps = mem_alloc(definition[0]*definition[1]*sizeof(double));
CHK(temps != NULL);
- OK(sdis_accum_buffer_map(buf, &accums));
-
/* Check the results validity */
- FOR_EACH(i, 0, layout.height * layout.width) {
- CHK(accums[i].nweights + accums[i].nfailures == SPP);
- CHK(accums[i].nfailures <= SPP/100);
- CHK(accums[i].sum_weights >= 0);
- CHK(accums[i].sum_weights_sqr >= 0);
+ FOR_EACH(iy, 0, definition[1]) {
+ FOR_EACH(ix, 0, definition[0]) {
+ const struct sdis_estimator* estimator;
+ struct sdis_mc T;
+ struct sdis_mc time;
+ size_t nreals;
+ size_t nfails;
+
+ BA(sdis_estimator_buffer_at(NULL, ix, iy, &estimator));
+ BA(sdis_estimator_buffer_at(buf, definition[0]+1, iy, &estimator));
+ BA(sdis_estimator_buffer_at(buf, ix, definition[1]+1, &estimator));
+ BA(sdis_estimator_buffer_at(buf, ix, iy, NULL));
+ OK(sdis_estimator_buffer_at(buf, ix, iy, &estimator));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ CHK(nreals + nfails == SPP);
+ CHK(T.E > 0);
+ CHK(time.E > 0);
+ }
}
/* Compute the per pixel temperature */
- FOR_EACH(iy, 0, layout.height) {
- const struct sdis_accum* row_accums = accums + iy * layout.width;
- double* row = temps + iy * layout.width;
- FOR_EACH(ix, 0, layout.width) {
- row[ix] = row_accums[ix].sum_weights / (double)row_accums[ix].nweights;
+ FOR_EACH(iy, 0, definition[1]) {
+ double* row = temps + iy * definition[0];
+ FOR_EACH(ix, 0, definition[0]) {
+ const struct sdis_estimator* estimator;
+ struct sdis_mc T;
+
+ OK(sdis_estimator_buffer_at(buf, ix, iy, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ row[ix] = T.E;
Tmax = MMAX(row[ix], Tmax);
Tmin = MMIN(row[ix], Tmin);
}
@@ -481,10 +501,10 @@ dump_image(const struct sdis_accum_buffer* buf)
OK(image_init(NULL, &img));
OK(image_setup(&img, IMG_WIDTH, IMG_HEIGHT, IMG_WIDTH*3, IMAGE_RGB8, NULL));
- FOR_EACH(iy, 0, layout.height) {
- const double* src_row = temps + iy*layout.width;
+ FOR_EACH(iy, 0, definition[1]) {
+ const double* src_row = temps + iy*definition[0];
char* dst_row = img.pixels + iy*img.pitch;
- FOR_EACH(ix, 0, layout.width) {
+ FOR_EACH(ix, 0, definition[0]) {
unsigned char* pixels = (unsigned char*)
(dst_row + ix * sizeof_image_format(img.format));
const unsigned char T = (unsigned char)
@@ -509,9 +529,11 @@ main(int argc, char** argv)
struct geometry geom = GEOMETRY_NULL;
struct s3dut_mesh* msh = NULL;
struct s3dut_mesh_data msh_data;
- struct sdis_accum_buffer* buf = NULL;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
struct sdis_camera* cam = NULL;
struct sdis_device* dev = NULL;
+ struct sdis_estimator_buffer* buf = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid0 = NULL;
struct sdis_medium* fluid1 = NULL;
@@ -522,9 +544,12 @@ main(int argc, char** argv)
struct solid solid_param = SOLID_NULL;
struct interf interface_param = INTERF_NULL;
size_t ntris, npos;
+ size_t nreals, nfails;
+ size_t definition[2];
double pos[3];
double tgt[3];
double up[3];
+ double trange[2] = {INF, INF};
(void)argc, (void)argv;
OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
@@ -590,6 +615,11 @@ main(int argc, char** argv)
geometry_get_interface, npos, geometry_get_position,
&geom, &scn));
+#if 0
+ dump_mesh(stdout, geom.positions, sa_size(geom.positions)/3, geom.indices,
+ sa_size(geom.indices)/3);
+#endif
+
/* Setup the camera */
d3(pos, 3, 3, 3);
d3(tgt, 0, 0, 0);
@@ -599,28 +629,80 @@ main(int argc, char** argv)
OK(sdis_camera_set_fov(cam, MDEG2RAD(70)));
OK(sdis_camera_look_at(cam, pos, tgt, up));
- /* Create the accum buffer */
- OK(sdis_accum_buffer_create(dev, IMG_WIDTH, IMG_HEIGHT, &buf));
-
#if 0
dump_mesh(stdout, geom.positions, npos, geom.indices, ntris);
exit(0);
#endif
+ BA(sdis_solve_camera(NULL, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, NULL, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, NULL, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 0, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, -1, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, 0, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, 0,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ 0, SDIS_HEAT_PATH_NONE, &buf));
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, NULL));
+
+ trange[0] = -1;
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ trange[0] = 10; trange[1] = 1;
+ BA(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+ trange[0] = trange[1] = INF;
/* Launch the simulation */
- OK(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
- sdis_accum_buffer_write, buf));
+ OK(sdis_solve_camera(scn, cam, trange, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT,
+ SPP, SDIS_HEAT_PATH_NONE, &buf));
+
+ BA(sdis_estimator_buffer_get_realisation_count(NULL, &nreals));
+ BA(sdis_estimator_buffer_get_realisation_count(buf, NULL));
+ OK(sdis_estimator_buffer_get_realisation_count(buf, &nreals));
+
+ BA(sdis_estimator_buffer_get_failure_count(NULL, &nfails));
+ BA(sdis_estimator_buffer_get_failure_count(buf, NULL));
+ OK(sdis_estimator_buffer_get_failure_count(buf, &nfails));
+
+ BA(sdis_estimator_buffer_get_temperature(NULL, &T));
+ BA(sdis_estimator_buffer_get_temperature(buf, NULL));
+ OK(sdis_estimator_buffer_get_temperature(buf, &T));
+
+ BA(sdis_estimator_buffer_get_realisation_time(NULL, &time));
+ BA(sdis_estimator_buffer_get_realisation_time(buf, NULL));
+ OK(sdis_estimator_buffer_get_realisation_time(buf, &time));
+
+ CHK(nreals + nfails == IMG_WIDTH*IMG_HEIGHT*SPP);
+
+ fprintf(stderr, "Overall temperature ~ %g +/- %g\n", T.E, T.SE);
+ fprintf(stderr, "Time per realisation (in usec) ~ %g +/- %g\n", time.E, time.SE);
+ fprintf(stderr, "#failures = %lu/%lu\n",
+ (unsigned long)nfails, (unsigned long)(IMG_WIDTH*IMG_HEIGHT*SPP));
+
+ BA(sdis_estimator_buffer_get_definition(NULL, definition));
+ BA(sdis_estimator_buffer_get_definition(buf, NULL));
+ OK(sdis_estimator_buffer_get_definition(buf, definition));
+ CHK(definition[0] == IMG_WIDTH);
+ CHK(definition[1] == IMG_HEIGHT);
/* Write the image */
dump_image(buf);
/* Release memory */
+ OK(sdis_estimator_buffer_ref_put(buf));
OK(sdis_scene_ref_put(scn));
OK(sdis_camera_ref_put(cam));
OK(sdis_interface_ref_put(interf0));
OK(sdis_interface_ref_put(interf1));
OK(sdis_device_ref_put(dev));
- OK(sdis_accum_buffer_ref_put(buf));
geometry_release(&geom);
check_memory_allocator(&allocator);
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -270,7 +270,7 @@ check_green_function(struct sdis_green_function* green)
}
void
-dump_heat_paths(FILE* stream, struct sdis_estimator* estimator)
+dump_heat_paths(FILE* stream, const struct sdis_estimator* estimator)
{
const struct sdis_heat_path* path;
size_t ipath;
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -292,7 +292,7 @@ check_green_function
extern LOCAL_SYM void
dump_heat_paths
(FILE* stream,
- struct sdis_estimator* estimator);
+ const struct sdis_estimator* estimator);
#endif /* TEST_SDIS_UTILS_H */
