commit df29d4ea741ce6653129dfc6b32d9e955cd2e268
parent 5301884f22d0bb7a8730624f6f399eeb28051748
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 25 Sep 2015 09:19:54 +0200
Code refactoring
Split the source in 2 parts: the sgf_device implementation and the
sgf_estimator one. Add comments.
Diffstat:
7 files changed, 544 insertions(+), 463 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -40,7 +40,7 @@ find_package(RCMake 0.1 REQUIRED)
find_package(RSys 0.2 REQUIRED)
find_package(StarSP 0.1 REQUIRED)
find_package(Star3D 0.2 REQUIRED)
-find_package(OpenMP)
+find_package(OpenMP REQUIRED)
include_directories(
${RSys_INCLUDE_DIR}
@@ -59,16 +59,18 @@ set(VERSION_MINOR 0)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
-set(SGF_FILES_SRC sgf.c)
+set(SGF_FILES_SRC sgf_device.c sgf_estimator.c)
set(SGF_FILES_INC_API sgf.h)
+set(SGF_FILES_INC sgf_device_c.h)
set(SGF_FILES_DOC COPYING.en COPYING.fr README.md)
# Prepend each file in the `SGF_FILES_<SRC|INC>' list by `SGF_SOURCE_DIR'
rcmake_prepend_path(SGF_FILES_SRC ${SGF_SOURCE_DIR})
rcmake_prepend_path(SGF_FILES_INC_API ${SGF_SOURCE_DIR})
+rcmake_prepend_path(SGF_FILES_INC ${SGF_SOURCE_DIR})
rcmake_prepend_path(SGF_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
-add_library(sgf SHARED ${SGF_FILES_SRC} ${SGF_FILES_INC_API})
+add_library(sgf SHARED ${SGF_FILES_SRC} ${SGF_FILES_INC_API} ${SGF_FILES_INC})
set_target_properties(sgf PROPERTIES
DEFINE_SYMBOL SGF_SHARED_BUILD
VERSION ${VERSION}
diff --git a/src/sgf.c b/src/sgf.c
@@ -1,459 +0,0 @@
-/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
- *
- * This software is governed by the CeCILL license under French law and
- * abiding by the rules of distribution of free software. You can use,
- * modify and/or redistribute the software under the terms of the CeCILL
- * license as circulated by CEA, CNRS and INRIA at the following URL
- * "http://www.cecill.info".
- *
- * As a counterpart to the access to the source code and rights to copy,
- * modify and redistribute granted by the license, users are provided only
- * with a limited warranty and the software's author, the holder of the
- * economic rights, and the successive licensors have only limited
- * liability.
- *
- * In this respect, the user's attention is drawn to the risks associated
- * with loading, using, modifying and/or developing or reproducing the
- * software by the user in light of its specific status of free software,
- * that may mean that it is complicated to manipulate, and that also
- * therefore means that it is reserved for developers and experienced
- * professionals having in-depth computer knowledge. Users are therefore
- * encouraged to load and test the software's suitability as regards their
- * requirements in conditions enabling the security of their systems and/or
- * data to be ensured and, more generally, to use and operate it in the
- * same conditions as regards security.
- *
- * The fact that you are presently reading this means that you have had
- * knowledge of the CeCILL license and that you accept its terms. */
-
-#include "sgf.h"
-
-#include <star/s3d.h>
-#include <star/ssp.h>
-
-#include <rsys/dynamic_array.h>
-#include <rsys/logger.h>
-#include <rsys/mem_allocator.h>
-#include <rsys/ref_count.h>
-
-struct accum {
- double radiative_flux;
- double sqr_radiative_flux;
-};
-
-#define DARRAY_NAME accum
-#define DARRAY_DATA struct accum
-#include <rsys/dynamic_array.h>
-
-struct sgf_device {
- int verbose;
- struct logger* logger;
-
- struct mem_allocator* allocator;
- ref_T ref;
-};
-
-struct sgf_estimator {
- struct darray_accum buf;
- size_t nsteps;
-
- struct sgf_device* dev;
- ref_T ref;
-};
-
-/* Define how many failures are handled until an error occurs */
-#define MAX_FAILURES 10
-
-/*******************************************************************************
- * Helper function
- ******************************************************************************/
-static res_T
-estimator_create(struct sgf_device* dev, struct sgf_estimator** out_estimator)
-{
- struct sgf_estimator* estimator = NULL;
- res_T res = RES_OK;
- ASSERT(dev && out_estimator);
-
- estimator = MEM_CALLOC(dev->allocator, 1, sizeof(struct sgf_estimator));
- if(!estimator) {
- if(dev->verbose) {
- logger_print(dev->logger, LOG_ERROR,
-"Not enough memory space: couldn't allocato the Gebhart Factor estimator.\n");
- }
- res = RES_MEM_ERR;
- goto error;
- }
- ref_init(&estimator->ref);
- SGF(device_ref_get(dev));
- estimator->dev = dev;
- darray_accum_init(dev->allocator, &estimator->buf);
-
-exit:
- *out_estimator = estimator;
- return res;
-
-error:
- if(estimator) {
- SGF(estimator_ref_put(estimator));
- estimator = NULL;
- }
- goto exit;
-}
-
-static INLINE char
-check_scene_desc(struct sgf_scene_desc* desc)
-{
- ASSERT(desc);
- return desc->get_material_property && desc->scene;
-}
-
-static FINLINE void
-accum_weight
- (struct sgf_estimator* estimator,
- const size_t iface,
- const double weight)
-{
- struct accum* acc;
- ASSERT(estimator);
- acc = darray_accum_data_get(&estimator->buf) + iface;
- acc->radiative_flux += weight;
- acc->sqr_radiative_flux += weight * weight;
-}
-
-/* Return RES_BAD_OP on numerical issue. i.e. the radiative path is skipped */
-static res_T
-gebhart_radiative_path
- (struct sgf_estimator* estimator,
- struct ssp_rng* rng,
- const size_t primitive_id,
- struct sgf_scene_desc* desc)
-{
- struct s3d_hit hit;
- struct s3d_primitive prim_from, prim;
- struct s3d_attrib attrib;
- float st[2] = {0.f, 0.f};
- double sky_gebhart = 0.0;
- double transmissivity = 1.0;
- double proba_reflec_spec;
- double emissivity;
- double reflectivity;
- double specularity;
-#ifndef NDEBUG
- double sum_radiative_flux = 0.f;
-#endif
- const double trans_min = 1.e-8;
- size_t nprims;
- float vec0[3];
- float normal[3];
- float pos[3];
- float dir[4];
- float range[2];
-
- /* Discard faces with no emissivity */
- emissivity = desc->get_material_property(desc->material,
- SGF_MATERIAL_EMISSIVITY, primitive_id, 0);
- if(emissivity <= 0.0)
- return RES_OK;
-
- S3D(scene_primitives_count(desc->scene, &nprims));
- S3D(scene_get_primitive(desc->scene, (unsigned)primitive_id, &prim));
-
- /* Compute the geometric normal of the primitive */
- S3D(primitive_get_attrib(&prim, S3D_GEOMETRY_NORMAL, st, &attrib));
- f3_normalize(normal, attrib.value);
-
- /* Uniformly sample the primitive to define the starting position of the
- * radiative path */
- S3D(primitive_sample
- (&prim, (float)ssp_rng_canonical(rng), (float)ssp_rng_canonical(rng), st));
- S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &attrib));
- f3_set(pos, attrib.value);
-
- /* Cosine weighted sampling of the starting radiative path direction around
- * `normal' */
- ssp_ran_hemisphere_cos(rng, normal, dir);
-
- transmissivity = 1.0;
- sky_gebhart = 0.0;
-
- for(;;) { /* Here we go */
- prim_from = prim;
-
- range[0] = 1.e-6f, range[1] = FLT_MAX;
- for(;;) {
- S3D(scene_trace_ray(desc->scene, pos, dir, range, &hit));
- if(S3D_PRIMITIVE_EQ(&hit.prim, &prim_from)) range[0] += 1.e-6f; /* Self hit */
- else break;
- }
-
- if(S3D_HIT_NONE(&hit)) {
- sky_gebhart += transmissivity;
- break;
- }
- f3_normalize(normal, hit.normal);
-
- /* Retrieve the hit position and the hit primitive id */
- f3_mulf(vec0, dir, hit.distance);
- f3_add(pos, vec0, pos);
- prim = hit.prim;
-
- /* Fetch material property */
- emissivity = desc->get_material_property(desc->material,
- SGF_MATERIAL_EMISSIVITY, prim.scene_prim_id, 0);
- specularity = desc->get_material_property(desc->material,
- SGF_MATERIAL_SPECULARITY, prim.scene_prim_id, 0);
- reflectivity = desc->get_material_property(desc->material,
- SGF_MATERIAL_REFLECTIVITY, prim.scene_prim_id, 0);
-
- if(transmissivity > trans_min) {
- const double weight = transmissivity * emissivity;
- accum_weight(estimator, prim.scene_prim_id, weight);
- transmissivity = transmissivity * (1.0 - emissivity);
-#ifndef NDEBUG
- sum_radiative_flux += weight;
-#endif
- } else {
- /* Russian roulette */
- if(ssp_rng_canonical(rng) < emissivity) {
- const double weight = transmissivity;
- accum_weight(estimator, prim.scene_prim_id, weight);
-#ifndef NDEBUG
- sum_radiative_flux += weight;
-#endif
- break;
- }
- }
-
- if(reflectivity > 0.0) {
- proba_reflec_spec = specularity / reflectivity;
- } else {
- break;
- }
- if(ssp_rng_canonical(rng) >= proba_reflec_spec) { /* Diffuse reflection */
- ssp_ran_hemisphere_cos(rng, normal, dir);
- ASSERT(f3_dot(normal, dir) > 0);
- } else { /* Specular reflection */
- const float tmp = -2.f * f3_dot(dir, normal);
- ASSERT(0);
- f3_mulf(vec0, normal, tmp);
- f3_add(dir, dir, vec0);
- f3_normalize(dir, dir);
- }
- }
-#if !defined(NDEBUG)
- /* Check the energy conservation property */
- ASSERT(eq_eps(sum_radiative_flux + sky_gebhart, 1.0, 1.e6));
-#endif
-
- return RES_OK;
-}
-
-static void
-device_release(ref_T* ref)
-{
- struct sgf_device* dev;
- ASSERT(ref);
- dev = CONTAINER_OF(ref, struct sgf_device, ref);
- MEM_RM(dev->allocator, dev);
-}
-
-static void
-estimator_release(ref_T* ref)
-{
- struct sgf_device* dev;
- struct sgf_estimator* estimator;
- ASSERT(ref);
- estimator = CONTAINER_OF(ref, struct sgf_estimator, ref);
- darray_accum_release(&estimator->buf);
- dev = estimator->dev;
- MEM_RM(dev->allocator, estimator);
- SGF(device_ref_put(dev));
-}
-
-/*******************************************************************************
- * API functions
- ******************************************************************************/
-res_T
-sgf_device_create
- (struct logger* log,
- struct mem_allocator* allocator,
- const int verbose,
- struct sgf_device** out_dev)
-{
- struct mem_allocator* mem_allocator;
- struct sgf_device* dev = NULL;
- struct logger* logger = NULL;
- res_T res = RES_OK;
-
- if(!out_dev) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- mem_allocator = allocator ? allocator : &mem_default_allocator;
- logger = log ? log : LOGGER_DEFAULT;
-
- dev = MEM_CALLOC(mem_allocator, 1, sizeof(struct sgf_device));
- if(!dev) {
- if(verbose) {
- logger_print(logger, LOG_ERROR,
- "Couldn't allocate the Star-Gebhart-Factor device.\n");
- }
- res = RES_MEM_ERR;
- goto error;
- }
- ref_init(&dev->ref);
- dev->allocator = mem_allocator;
- dev->logger = logger;
- dev->verbose = verbose;
-
-exit:
- if(out_dev) *out_dev = dev;
- return res;
-error:
- if(dev) {
- SGF(device_ref_put(dev));
- dev = NULL;
- }
- goto exit;
-}
-
-res_T
-sgf_device_ref_get(struct sgf_device* dev)
-{
- if(!dev) return RES_BAD_ARG;
- ref_get(&dev->ref);
- return RES_OK;
-}
-
-res_T
-sgf_device_ref_put(struct sgf_device* dev)
-{
- if(!dev) return RES_BAD_ARG;
- ref_put(&dev->ref, device_release);
- return RES_OK;
-}
-
-res_T
-sgf_integrate
- (struct sgf_device* dev,
- struct ssp_rng* rng,
- const size_t primitive_id,
- struct sgf_scene_desc* desc,
- const size_t steps_count,
- struct sgf_estimator** out_estimator)
-{
- struct sgf_estimator* estimator = NULL;
- size_t istep = 0;
- size_t nprims;
- int mask;
- res_T res = RES_OK;
-
- if(!dev || !steps_count || !rng || !desc || !check_scene_desc(desc)
- || !out_estimator) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- res = estimator_create(dev, &estimator);
- if(res != RES_OK) goto error;
-
- /* Check scene active sessions */
- S3D(scene_get_session_mask(desc->scene, &mask));
- if(!(mask & S3D_TRACE) || !(mask & S3D_GET_PRIMITIVE)) {
- if(dev->verbose) {
- logger_print(dev->logger, LOG_ERROR,
- "No active trace/get_primitive session on the Star-3D scene\n");
- }
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Check submitted primitive_id */
- S3D(scene_primitives_count(desc->scene, &nprims));
- if(primitive_id >= nprims) {
- if(dev->verbose) {
- logger_print(dev->logger, LOG_ERROR,
- "Invalid primitive index `%lu'\n", (unsigned long)primitive_id);
- }
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Allocate internal memory space */
- res = darray_accum_resize(&estimator->buf, nprims*desc->spectral_bands_count);
- if(res != RES_OK) {
- if(dev->verbose) {
- logger_print(dev->logger, LOG_ERROR,
-"Couldn't allocate the Gebhart Factor integration data structure.\n");
- }
- goto error;
- }
- memset(darray_accum_data_get(&estimator->buf), 0,
- darray_accum_size_get(&estimator->buf)*sizeof(struct accum));
-
- /* Here we go ! Perform the Gebhart factor integration.
- * TODO interate per spectral bands */
- FOR_EACH(istep, 0, steps_count) {
- res = gebhart_radiative_path(estimator, rng, primitive_id, desc);
- if(res != RES_OK) goto error;
- }
- estimator->nsteps = steps_count;
-
-exit:
- if(out_estimator) *out_estimator = estimator;
- return res;
-error:
- if(estimator) SGF(estimator_ref_put(estimator));
- goto exit;
-}
-
-res_T
-sgf_estimator_ref_get(struct sgf_estimator* estimator)
-{
- if(!estimator) return RES_BAD_ARG;
- ref_get(&estimator->ref);
- return RES_OK;
-}
-
-res_T
-sgf_estimator_ref_put(struct sgf_estimator* estimator)
-{
- if(!estimator) return RES_BAD_ARG;
- ref_put(&estimator->ref, estimator_release);
- return RES_OK;
-}
-
-res_T
-sgf_estimator_get_status
- (struct sgf_estimator* estimator,
- const size_t iprimitive,
- struct sgf_status* status)
-{
- struct accum* acc;
- size_t nsteps;
-
- if(!estimator || !status)
- return RES_BAD_ARG;
-
- if(iprimitive >= darray_accum_size_get(&estimator->buf)) {
- if(estimator->dev->verbose) {
- logger_print(estimator->dev->logger, LOG_ERROR,
- "Invalid primitive index `%lu'\n", iprimitive);
- }
- return RES_BAD_ARG;
- }
-
- acc = darray_accum_data_get(&estimator->buf) + iprimitive;
- nsteps = estimator->nsteps;
-
- status->E = acc->radiative_flux / (double)nsteps;
- status->V =
- acc->sqr_radiative_flux / (double)nsteps
- - (acc->radiative_flux * acc->radiative_flux) / (double)(nsteps * nsteps);
- status->SE = sqrt(status->V / (double)nsteps);
- status->nsteps = nsteps;
-
- return RES_OK;
-}
-
diff --git a/src/sgf.h b/src/sgf.h
@@ -78,6 +78,7 @@ static const struct sgf_scene_desc SGF_SCENE_DESC_NULL = {
NULL, NULL, 0, 0
};
+/* Estimated Gebart Factor between 2 faces */
struct sgf_status {
double E; /* Expected value */
double V; /* Variance */
@@ -85,9 +86,13 @@ struct sgf_status {
size_t nsteps; /* #realisations */
};
+/* Opaque types */
struct sgf_device;
struct sgf_estimator;
+/*******************************************************************************
+ * Star Gebhart Factor API
+ ******************************************************************************/
BEGIN_DECLS
SGF_API res_T
@@ -105,6 +110,7 @@ SGF_API res_T
sgf_device_ref_put
(struct sgf_device* dev);
+/* TODO comment it */
SGF_API res_T
sgf_integrate
(struct sgf_device* dev,
diff --git a/src/sgf_device.c b/src/sgf_device.c
@@ -0,0 +1,110 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#include "sgf.h"
+#include "sgf_device_c.h"
+
+#include <rsys/logger.h>
+#include <rsys/mem_allocator.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static void
+device_release(ref_T* ref)
+{
+ struct sgf_device* dev;
+ ASSERT(ref);
+ dev = CONTAINER_OF(ref, struct sgf_device, ref);
+ MEM_RM(dev->allocator, dev);
+}
+
+/*******************************************************************************
+ * API functions
+ ******************************************************************************/
+res_T
+sgf_device_create
+ (struct logger* log,
+ struct mem_allocator* allocator,
+ const int verbose,
+ struct sgf_device** out_dev)
+{
+ struct mem_allocator* mem_allocator;
+ struct sgf_device* dev = NULL;
+ struct logger* logger = NULL;
+ res_T res = RES_OK;
+
+ if(!out_dev) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ mem_allocator = allocator ? allocator : &mem_default_allocator;
+ logger = log ? log : LOGGER_DEFAULT;
+
+ dev = MEM_CALLOC(mem_allocator, 1, sizeof(struct sgf_device));
+ if(!dev) {
+ if(verbose) {
+ logger_print(logger, LOG_ERROR,
+ "Couldn't allocate the Star-Gebhart-Factor device.\n");
+ }
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&dev->ref);
+ dev->allocator = mem_allocator;
+ dev->logger = logger;
+ dev->verbose = verbose;
+
+exit:
+ if(out_dev) *out_dev = dev;
+ return res;
+error:
+ if(dev) {
+ SGF(device_ref_put(dev));
+ dev = NULL;
+ }
+ goto exit;
+}
+
+res_T
+sgf_device_ref_get(struct sgf_device* dev)
+{
+ if(!dev) return RES_BAD_ARG;
+ ref_get(&dev->ref);
+ return RES_OK;
+}
+
+res_T
+sgf_device_ref_put(struct sgf_device* dev)
+{
+ if(!dev) return RES_BAD_ARG;
+ ref_put(&dev->ref, device_release);
+ return RES_OK;
+}
+
diff --git a/src/sgf_device_c.h b/src/sgf_device_c.h
@@ -0,0 +1,49 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#ifndef SGF_DEVICE_C_H
+#define SGF_DEVICE_C_H
+
+#include <rsys/ref_count.h>
+
+struct logger;
+struct mem_allocator;
+
+struct sgf_device {
+ int verbose;
+ struct logger* logger;
+
+ struct mem_allocator* allocator;
+ ref_T ref;
+};
+
+/* TODO implement a print function that print messages with respect to the
+ * device verbosity */
+
+#endif /* SGF_DEVICE_C_H */
+
diff --git a/src/sgf_estimator.c b/src/sgf_estimator.c
@@ -0,0 +1,373 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#include "sgf.h"
+#include "sgf_device_c.h"
+
+#include <star/s3d.h>
+#include <star/ssp.h>
+
+#include <rsys/dynamic_array.h>
+#include <rsys/logger.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+struct accum {
+ double radiative_flux;
+ double sqr_radiative_flux;
+};
+
+#define DARRAY_NAME accum
+#define DARRAY_DATA struct accum
+#include <rsys/dynamic_array.h>
+
+/* Estimator of the Gebhart Factor of a given face to all the other ones */
+struct sgf_estimator {
+ struct darray_accum buf; /* Integrated radiative flux of each primitive */
+ size_t nsteps; /* Number of realisations */
+
+ struct sgf_device* dev;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE char
+check_scene_desc(struct sgf_scene_desc* desc)
+{
+ ASSERT(desc);
+ return desc->get_material_property && desc->scene;
+}
+
+static FINLINE void
+accum_weight
+ (struct sgf_estimator* estimator,
+ const size_t iface,
+ const double weight)
+{
+ struct accum* acc;
+ ASSERT(estimator);
+ acc = darray_accum_data_get(&estimator->buf) + iface;
+ acc->radiative_flux += weight;
+ acc->sqr_radiative_flux += weight * weight;
+}
+
+/* Return RES_BAD_OP on numerical issue. i.e. the radiative path is skipped */
+static res_T
+gebhart_radiative_path
+ (struct sgf_estimator* estimator,
+ struct ssp_rng* rng,
+ const size_t primitive_id,
+ struct sgf_scene_desc* desc)
+{
+ struct s3d_hit hit;
+ struct s3d_primitive prim_from, prim;
+ struct s3d_attrib attrib;
+
+ const double trans_min = 1.e-8; /* Minimum transmissivity threshold */
+ double outside_radiative_flux = 0.0;
+ double proba_reflec_spec;
+ double transmissivity, emissivity, reflectivity, specularity;
+#ifndef NDEBUG
+ double sum_radiative_flux = 0.f;
+#endif
+
+ float st[2] = {0.f, 0.f}; /* Parametric coordinate of the sampled position */
+ float vec0[3]; /* Temporary vector */
+ float normal[3]; /* Geometric normal */
+ float pos[3]; /* Radiative path position */
+ float dir[4]; /* Radiative path direction. dir[3] <=> sampled dir pdf */
+ float range[2]; /* Traced ray range */
+
+ /* Discard faces with no emissivity */
+ emissivity = desc->get_material_property(desc->material,
+ SGF_MATERIAL_EMISSIVITY, primitive_id, 0);
+ if(emissivity <= 0.0)
+ return RES_OK;
+
+ /* Get the geometric normal of the input primitive */
+ S3D(scene_get_primitive(desc->scene, (unsigned)primitive_id, &prim));
+ S3D(primitive_get_attrib(&prim, S3D_GEOMETRY_NORMAL, st, &attrib));
+ f3_normalize(normal, attrib.value);
+
+ /* Uniformly sample the primitive to define the starting position of the
+ * radiative path */
+ S3D(primitive_sample(&prim,
+ (float)ssp_rng_canonical(rng), (float)ssp_rng_canonical(rng), st));
+ S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &attrib));
+ f3_set(pos, attrib.value);
+
+ /* Cosine weighted sampling of the radiative path direction around `normal' */
+ ssp_ran_hemisphere_cos(rng, normal, dir);
+
+ transmissivity = 1.0;
+ outside_radiative_flux = 0.0;
+ for(;;) { /* Here we go */
+
+ prim_from = prim;
+ range[0] = 1.e-6f, range[1] = FLT_MAX;
+
+ for(;;) {
+ S3D(scene_trace_ray(desc->scene, pos, dir, range, &hit));
+ if(S3D_PRIMITIVE_EQ(&hit.prim, &prim_from)) range[0] += 1.e-6f; /* Self hit */
+ else break;
+ }
+
+ if(S3D_HIT_NONE(&hit)) { /* The ray is outside the volume */
+ outside_radiative_flux += transmissivity;
+ break;
+ }
+
+ /* Retrieve the hit position and the hit primitive id */
+ f3_mulf(vec0, dir, hit.distance);
+ f3_add(pos, vec0, pos);
+ prim = hit.prim;
+
+ /* Fetch material property */
+ emissivity = desc->get_material_property(desc->material,
+ SGF_MATERIAL_EMISSIVITY, prim.scene_prim_id, 0);
+ specularity = desc->get_material_property(desc->material,
+ SGF_MATERIAL_SPECULARITY, prim.scene_prim_id, 0);
+ reflectivity = desc->get_material_property(desc->material,
+ SGF_MATERIAL_REFLECTIVITY, prim.scene_prim_id, 0);
+
+ if(transmissivity > trans_min) {
+ const double weight = transmissivity * emissivity;
+ accum_weight(estimator, prim.scene_prim_id, weight);
+ transmissivity = transmissivity * (1.0 - emissivity);
+#ifndef NDEBUG
+ sum_radiative_flux += weight;
+#endif
+ } else {
+ /* Russian roulette */
+ if(ssp_rng_canonical(rng) < emissivity) {
+ const double weight = transmissivity;
+ accum_weight(estimator, prim.scene_prim_id, weight);
+#ifndef NDEBUG
+ sum_radiative_flux += weight;
+#endif
+ break;
+ }
+ }
+
+ if(reflectivity <= 0.0)
+ break;
+
+ proba_reflec_spec = specularity / reflectivity;
+ f3_normalize(normal, hit.normal);
+
+ if(ssp_rng_canonical(rng) >= proba_reflec_spec) { /* Diffuse reflection */
+ ssp_ran_hemisphere_cos(rng, normal, dir);
+ ASSERT(f3_dot(normal, dir) > 0);
+ } else { /* Specular reflection */
+ const float tmp = -2.f * f3_dot(dir, normal);
+ ASSERT(0);
+ f3_mulf(vec0, normal, tmp);
+ f3_add(dir, dir, vec0);
+ f3_normalize(dir, dir);
+ }
+ }
+#if !defined(NDEBUG)
+ /* Check the energy conservation property */
+ ASSERT(eq_eps(sum_radiative_flux + outside_radiative_flux, 1.0, 1.e6));
+#endif
+ return RES_OK;
+}
+
+static res_T
+estimator_create(struct sgf_device* dev, struct sgf_estimator** out_estimator)
+{
+ struct sgf_estimator* estimator = NULL;
+ res_T res = RES_OK;
+ ASSERT(dev && out_estimator);
+
+ estimator = MEM_CALLOC(dev->allocator, 1, sizeof(struct sgf_estimator));
+ if(!estimator) {
+ if(dev->verbose) {
+ logger_print(dev->logger, LOG_ERROR,
+"Not enough memory space: couldn't allocato the Gebhart Factor estimator.\n");
+ }
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&estimator->ref);
+ SGF(device_ref_get(dev));
+ estimator->dev = dev;
+ darray_accum_init(dev->allocator, &estimator->buf);
+
+exit:
+ *out_estimator = estimator;
+ return res;
+
+error:
+ if(estimator) {
+ SGF(estimator_ref_put(estimator));
+ estimator = NULL;
+ }
+ goto exit;
+}
+
+static void
+estimator_release(ref_T* ref)
+{
+ struct sgf_device* dev;
+ struct sgf_estimator* estimator;
+ ASSERT(ref);
+ estimator = CONTAINER_OF(ref, struct sgf_estimator, ref);
+ darray_accum_release(&estimator->buf);
+ dev = estimator->dev;
+ MEM_RM(dev->allocator, estimator);
+ SGF(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sgf_integrate
+ (struct sgf_device* dev,
+ struct ssp_rng* rng,
+ const size_t primitive_id,
+ struct sgf_scene_desc* desc,
+ const size_t steps_count,
+ struct sgf_estimator** out_estimator)
+{
+ struct sgf_estimator* estimator = NULL;
+ size_t istep = 0;
+ size_t nprims;
+ int mask;
+ res_T res = RES_OK;
+
+ if(!dev || !steps_count || !rng || !desc || !check_scene_desc(desc)
+ || !out_estimator) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = estimator_create(dev, &estimator);
+ if(res != RES_OK) goto error;
+
+ /* Check scene active sessions */
+ S3D(scene_get_session_mask(desc->scene, &mask));
+ if(!(mask & S3D_TRACE) || !(mask & S3D_GET_PRIMITIVE)) {
+ if(dev->verbose) {
+ logger_print(dev->logger, LOG_ERROR,
+ "No active trace/get_primitive session on the Star-3D scene\n");
+ }
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Check submitted primitive_id */
+ S3D(scene_primitives_count(desc->scene, &nprims));
+ if(primitive_id >= nprims) {
+ if(dev->verbose) {
+ logger_print(dev->logger, LOG_ERROR,
+ "Invalid primitive index `%lu'\n", (unsigned long)primitive_id);
+ }
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Allocate internal memory space */
+ res = darray_accum_resize(&estimator->buf, nprims*desc->spectral_bands_count);
+ if(res != RES_OK) {
+ if(dev->verbose) {
+ logger_print(dev->logger, LOG_ERROR,
+"Couldn't allocate the Gebhart Factor integration data structure.\n");
+ }
+ goto error;
+ }
+ memset(darray_accum_data_get(&estimator->buf), 0,
+ darray_accum_size_get(&estimator->buf)*sizeof(struct accum));
+
+ /* Invoke the Gebhart factor integration. TODO interate per spectral bands */
+ FOR_EACH(istep, 0, steps_count) {
+ res = gebhart_radiative_path(estimator, rng, primitive_id, desc);
+ if(res != RES_OK) goto error;
+ }
+ estimator->nsteps = steps_count;
+
+exit:
+ if(out_estimator) *out_estimator = estimator;
+ return res;
+error:
+ if(estimator) SGF(estimator_ref_put(estimator));
+ goto exit;
+}
+
+res_T
+sgf_estimator_ref_get(struct sgf_estimator* estimator)
+{
+ if(!estimator) return RES_BAD_ARG;
+ ref_get(&estimator->ref);
+ return RES_OK;
+}
+
+res_T
+sgf_estimator_ref_put(struct sgf_estimator* estimator)
+{
+ if(!estimator) return RES_BAD_ARG;
+ ref_put(&estimator->ref, estimator_release);
+ return RES_OK;
+}
+
+res_T
+sgf_estimator_get_status
+ (struct sgf_estimator* estimator,
+ const size_t iprimitive,
+ struct sgf_status* status)
+{
+ struct accum* acc;
+ size_t nsteps;
+
+ if(!estimator || !status)
+ return RES_BAD_ARG;
+
+ if(iprimitive >= darray_accum_size_get(&estimator->buf)) {
+ if(estimator->dev->verbose) {
+ logger_print(estimator->dev->logger, LOG_ERROR,
+ "Invalid primitive index `%lu'\n", iprimitive);
+ }
+ return RES_BAD_ARG;
+ }
+
+ acc = darray_accum_data_get(&estimator->buf) + iprimitive;
+ nsteps = estimator->nsteps;
+
+ status->E = acc->radiative_flux / (double)nsteps; /* Expected value */
+ status->V = /* Variance */
+ acc->sqr_radiative_flux / (double)nsteps
+ - (acc->radiative_flux * acc->radiative_flux) / (double)(nsteps * nsteps);
+ status->SE = sqrt(status->V / (double)nsteps); /* Standard error */
+ status->nsteps = nsteps; /* # realisations */
+
+ return RES_OK;
+}
+
diff --git a/src/test_sgf_cube.c b/src/test_sgf_cube.c
@@ -37,7 +37,7 @@
#include <omp.h>
-#define NSTEPS 100000L
+#define NSTEPS 10000L
/*
* Analytic Gebhart Factor matrix
