star-gf

Compute Gebhart factors
git clone git://git.meso-star.fr/star-gf.git
Log | Files | Refs | README | LICENSE

commit 34368761dc359da86f418c7b2b52edefab6a9212
parent 385712a6f0fe3098deb9c78217eb479f7b962908
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 24 Sep 2015 09:36:32 +0200

Major refactoring of the SGF API

Diffstat:
Mcmake/CMakeLists.txt | 17++++++++++-------
Msrc/sgf.c | 398++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------
Msrc/sgf.h | 85++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------
Dsrc/sgf_accum.h | 122-------------------------------------------------------------------------------
Msrc/test_sgf_cube.c | 190+++++++++++++++++++++++++++++++++++++------------------------------------------
Asrc/test_sgf_device.c | 86+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/test_sgf_utils.h | 6+++---
7 files changed, 554 insertions(+), 350 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -39,13 +39,12 @@ option(NO_TEST "Disable the test" OFF) find_package(RCMake 0.1 REQUIRED) find_package(RSys 0.2 REQUIRED) find_package(StarSP 0.1 REQUIRED) -find_package(StarMC 0.2 REQUIRED) find_package(Star3D 0.2 REQUIRED) +find_package(OpenMP) include_directories( ${RSys_INCLUDE_DIR} ${StarSP_INCLUDE_DIR} - ${StarMC_INCLUDE_DIR} ${Star3D_INCLUDE_DIR}) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR}) @@ -62,22 +61,20 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) set(SGF_FILES_SRC sgf.c) set(SGF_FILES_INC_API sgf.h) -set(SGF_FILES_INC sgf_accum.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 ${SGF_SOURCE_DIR}) rcmake_prepend_path(SGF_FILES_INC_API ${SGF_SOURCE_DIR}) rcmake_prepend_path(SGF_FILES_DOC ${PROJECT_SOURCE_DIR}/../) -add_library(sgf SHARED ${SGF_FILES_SRC} ${SGF_FILES_INC} ${SGF_FILES_INC_API}) +add_library(sgf SHARED ${SGF_FILES_SRC} ${SGF_FILES_INC_API}) set_target_properties(sgf PROPERTIES DEFINE_SYMBOL SGF_SHARED_BUILD VERSION ${VERSION} SOVERSION ${VERSION_MAJOR}) -target_link_libraries(sgf RSys Star3D StarMC StarSP) +target_link_libraries(sgf RSys Star3D StarSP) if(CMAKE_COMPILER_IS_GNUCC) target_link_libraries(sgf m) endif() @@ -93,8 +90,14 @@ if(NOT NO_TEST) target_link_libraries(${_name} sgf RSys) add_test(${_name} ${_name}) endfunction(new_test) + new_test(test_sgf_device) new_test(test_sgf_cube) - new_test(test_sgf_tetrahedron) + # new_test(test_sgf_tetrahedron) + + set_target_properties(test_sgf_cube PROPERTIES COMPILE_FLAGS ${OpenMP_C_FLAGS}) + if(CMAKE_COMPILER_IS_GNUCC) + set_target_properties(test_sgf_cube PROPERTIES LINK_FLAGS ${OpenMP_C_FLAGS}) + endif() endif() ################################################################################ diff --git a/src/sgf.c b/src/sgf.c @@ -27,11 +27,43 @@ * knowledge of the CeCILL license and that you accept its terms. */ #include "sgf.h" -#include "sgf_accum.h" #include <star/s3d.h> #include <star/ssp.h> +#include <rsys/dynamic_array_size_t.h> +#include <rsys/hash_table.h> +#include <rsys/logger.h> +#include <rsys/mem_allocator.h> +#include <rsys/mutex.h> +#include <rsys/ref_count.h> + +struct accum { + double radiative_flux; + double sqr_radiative_flux; +}; + +#define HTABLE_NAME accum +#define HTABLE_DATA struct accum +#define HTABLE_KEY uint64_t +#include <rsys/hash_table.h> + +#define CELL_KEY(From, To) (((uint64_t)From<<32) | (uint64_t)(To)) + +struct sgf_device { + struct htable_accum matrix; + struct darray_size_t nsteps; + struct mutex* mutex; + + int verbose; + struct logger* logger; + struct mem_allocator* allocator; + + struct sgf_scene_desc scn_desc; + int is_integrating; + ref_T ref; +}; + /* Define how many failures are handled until an error occurs */ #define MAX_FAILURES 10 @@ -39,24 +71,48 @@ * Helper function ******************************************************************************/ static INLINE char -check_context(struct sgf_context* ctx) +check_scene_desc(struct sgf_scene_desc* desc) { - ASSERT(ctx); - return ctx->get_material_property - && ctx->scene - && ctx->logger; + ASSERT(desc); + return desc->get_material_property && desc->scene; +} + +static FINLINE res_T +accum_weight + (struct sgf_device* dev, + const size_t ifrom, + const size_t ito, + const double weight) +{ + struct accum* acc; + uint64_t key; + res_T res = RES_OK; + ASSERT(dev); + + key = CELL_KEY(ifrom, ito); + mutex_lock(dev->mutex); + acc = htable_accum_find(&dev->matrix, &key); + if(acc) { + acc->radiative_flux += weight; + acc->sqr_radiative_flux += weight * weight; + } else { + struct accum tmp; + tmp.radiative_flux = weight; + tmp.sqr_radiative_flux = weight * weight; + res = htable_accum_set(&dev->matrix, &key, &tmp); + } + mutex_unlock(dev->mutex); + return res; } /* Return RES_BAD_OP on numerical issue. i.e. the radiative path is skipped */ static res_T gebhart_radiative_path - (struct sgf_accum* result, - struct ssp_rng* rng, - struct sgf_context* ctx, - const size_t ispectral_band) + (struct sgf_device* dev, + const size_t primitive_id, + struct ssp_rng* rng) { struct s3d_hit hit; - struct sgf_accum* gebfac; /* Gebhart factor */ struct s3d_primitive prim_from, prim; struct s3d_attrib attrib; float st[2] = {0.f, 0.f}; @@ -76,16 +132,16 @@ gebhart_radiative_path float pos[3]; float dir[4]; float range[2]; - ASSERT(ispectral_band < ctx->spectral_bands_count); + res_T res = RES_OK; /* Discard faces with no emissivity */ - emissivity = ctx->get_material_property(ctx->material, - SGF_MATERIAL_EMISSIVITY, ctx->primitive_id, ispectral_band); + emissivity = dev->scn_desc.get_material_property(dev->scn_desc.material, + SGF_MATERIAL_EMISSIVITY, primitive_id, 0); if(emissivity <= 0.0) return RES_OK; - S3D(scene_primitives_count(ctx->scene, &nprims)); - S3D(scene_get_primitive(ctx->scene, (unsigned)ctx->primitive_id, &prim)); + S3D(scene_primitives_count(dev->scn_desc.scene, &nprims)); + S3D(scene_get_primitive(dev->scn_desc.scene, (unsigned)primitive_id, &prim)); /* Compute the geometric normal of the primitive */ S3D(primitive_get_attrib(&prim, S3D_GEOMETRY_NORMAL, st, &attrib)); @@ -110,10 +166,9 @@ gebhart_radiative_path range[0] = 1.e-6f, range[1] = FLT_MAX; for(;;) { - S3D(scene_trace_ray(ctx->scene, pos, dir, range, &hit)); + S3D(scene_trace_ray(dev->scn_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)) { @@ -128,28 +183,27 @@ gebhart_radiative_path prim = hit.prim; /* Fetch material property */ - emissivity = ctx->get_material_property(ctx->material, - SGF_MATERIAL_EMISSIVITY, prim.scene_prim_id, ispectral_band); - specularity = ctx->get_material_property(ctx->material, - SGF_MATERIAL_SPECULARITY, prim.scene_prim_id, ispectral_band); - reflectivity = ctx->get_material_property(ctx->material, - SGF_MATERIAL_REFLECTIVITY, prim.scene_prim_id, ispectral_band); - - gebfac = result + prim.scene_prim_id * ctx->spectral_bands_count; + emissivity = dev->scn_desc.get_material_property(dev->scn_desc.material, + SGF_MATERIAL_EMISSIVITY, prim.scene_prim_id, 0); + specularity = dev->scn_desc.get_material_property(dev->scn_desc.material, + SGF_MATERIAL_SPECULARITY, prim.scene_prim_id, 0); + reflectivity = dev->scn_desc.get_material_property(dev->scn_desc.material, + SGF_MATERIAL_REFLECTIVITY, prim.scene_prim_id, 0); + if(transmissivity > trans_min) { const double weight = transmissivity * emissivity; + res = accum_weight(dev, primitive_id, prim.scene_prim_id, weight); + if(res != RES_OK) goto error; + transmissivity = transmissivity * (1.0 - emissivity); #ifndef NDEBUG sum_radiative_flux += weight; #endif - gebfac[ispectral_band].radiative_flux += weight; - gebfac[ispectral_band].sqr_radiative_flux += weight * weight; - transmissivity = transmissivity * (1.0 - emissivity); } else { /* Russian roulette */ if(ssp_rng_canonical(rng) < emissivity) { const double weight = transmissivity; - gebfac[ispectral_band].radiative_flux += weight; - gebfac[ispectral_band].sqr_radiative_flux += weight * weight; + res = accum_weight(dev, primitive_id, prim.scene_prim_id, weight); + if(res != RES_OK) goto error; #ifndef NDEBUG sum_radiative_flux += weight; #endif @@ -177,108 +231,266 @@ gebhart_radiative_path /* Check the energy conservation property */ ASSERT(eq_eps(sum_radiative_flux + sky_gebhart, 1.0, 1.e6)); #endif - return RES_OK; + +exit: + return res; +error: + goto exit; } static void -gebhart_factor_integrand - (void* result, - struct ssp_rng* rng, - void* context) +device_release(ref_T* ref) { - struct sgf_context* ctx = context; - struct sgf_accum* buf; - size_t iband = 0; - res_T res = RES_OK; (void)res; - ASSERT(result && rng && context); - - buf = darray_gfacc_data_get(result); - FOR_EACH(iband, 0, ctx->spectral_bands_count) { - res = gebhart_radiative_path(buf, rng, ctx, iband); - ASSERT(res == RES_OK); - } + struct sgf_device* dev; + ASSERT(ref); + dev = CONTAINER_OF(ref, struct sgf_device, ref); + htable_accum_release(&dev->matrix); + darray_size_t_release(&dev->nsteps); + if(dev->mutex) mutex_destroy(dev->mutex); + MEM_RM(dev->allocator, dev); } /******************************************************************************* * API functions ******************************************************************************/ res_T -sgf_integrate - (struct smc_device* smc, - struct sgf_context* context, - const size_t steps_count, - struct sgf_accum* accbuf) +sgf_device_create + (struct logger* log, + struct mem_allocator* allocator, + const int verbose, + struct sgf_device** out_dev) { - struct smc_integrator_accum integrator; - struct smc_accumulator* accum = NULL; - struct smc_accumulator_status status; - size_t i; - size_t nprims; - int mask; + struct mem_allocator* mem_allocator; + struct sgf_device* dev = NULL; + struct logger* logger = NULL; res_T res = RES_OK; - if(!smc || !context || !accbuf || !check_context(context) || !steps_count) { + if(!out_dev) { res = RES_BAD_ARG; goto error; } - integrator.integrand = gebhart_factor_integrand; - integrator.type = &accum_buffer_type; - integrator.max_steps = steps_count; + mem_allocator = allocator ? allocator : &mem_default_allocator; + logger = log ? log : LOGGER_DEFAULT; - S3D(scene_get_session_mask(context->scene, &mask)); - if((mask & S3D_TRACE) == 0) { - if(context->verbose) { - logger_print(context->logger, LOG_ERROR, - "No active \"trace\" session on the Star-3D scene\n"); + 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_BAD_OP; + res = RES_MEM_ERR; goto error; } - if((mask & S3D_GET_PRIMITIVE) == 0) { - if(context->verbose) { - logger_print(context->logger, LOG_ERROR, - "No active \"get primitive\" session on the Star-3D scene\n"); + ref_init(&dev->ref); + dev->allocator = mem_allocator; + dev->logger = logger; + dev->verbose = verbose; + htable_accum_init(dev->allocator, &dev->matrix); + darray_size_t_init(dev->allocator, &dev->nsteps); + + dev->mutex = mutex_create(); + if(!dev->mutex) { + if(verbose) { + logger_print(logger, LOG_ERROR, "Couldn't create the SGF mutex.\n"); } - res = RES_BAD_OP; + res = RES_MEM_ERR; goto error; } - S3D(scene_primitives_count(context->scene, &nprims)); - if(context->primitive_id >= nprims) { - if(context->verbose) { - logger_print(context->logger, LOG_ERROR, - "Invalid primitive index `%lu'\n", - (unsigned long)context->primitive_id); +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_begin_integrate(struct sgf_device* dev, struct sgf_scene_desc* desc) +{ + res_T res = RES_OK; + size_t nprims; + + if(!dev || !desc || !check_scene_desc(desc)) return RES_BAD_ARG; + + if(dev->is_integrating) { + if(dev->verbose) { + logger_print(dev->logger, LOG_ERROR, + "A Gebhart Factor integration is already running.\n"); } - res = RES_BAD_ARG; - goto error; + return RES_BAD_OP; } - res = smc_integrate(smc, &integrator, context, &accum); + res = s3d_scene_begin_session(desc->scene, S3D_TRACE|S3D_GET_PRIMITIVE); if(res != RES_OK) { - if(context->verbose) { - logger_print(context->logger, LOG_ERROR, - "Couldn't integrate the radiative flux\n"); + if(dev->verbose) { + logger_print(dev->logger, LOG_ERROR, + "Couldn't begin a Star-3D sesion on the submitted scene\n"); } - goto error; + return res; } - res = smc_accumulator_get_status(accum, &status); + dev->scn_desc = *desc; + S3D(scene_primitives_count(dev->scn_desc.scene, &nprims)); + + res = darray_size_t_resize(&dev->nsteps, nprims); if(res != RES_OK) { - if(context->verbose) { - logger_print(context->logger, LOG_ERROR, - "Couldn't get the radiative flux integration status\n"); + if(dev->verbose) { + logger_print(dev->logger, LOG_ERROR, +"Couldn't allocate the Gebhart Factor integration data structure.\n"); + } + S3D(scene_end_session(desc->scene)); + return res; + } + + memset(darray_size_t_data_get(&dev->nsteps), 0, nprims * sizeof(size_t)); + htable_accum_clear(&dev->matrix); + S3D(scene_ref_get(dev->scn_desc.scene)); + dev->is_integrating = 1; + + return RES_OK; +} + +res_T +sgf_end_integrate(struct sgf_device* dev) +{ + res_T res = RES_OK; + + if(!dev) return RES_BAD_ARG; + if(!dev->is_integrating) { + if(dev->verbose) { + logger_print(dev->logger, LOG_ERROR, + "No Gebhart Factor integration is running\n"); + } + return RES_BAD_OP; + } + res = s3d_scene_end_session(dev->scn_desc.scene); + if(res != RES_OK) { + if(dev->verbose) { + logger_print(dev->logger, LOG_ERROR, + "Couldn't stop the Star-3D session.\n"); + } + return res; + } + + S3D(scene_ref_put(dev->scn_desc.scene)); + dev->is_integrating = 0; + htable_accum_clear(&dev->matrix); + return RES_OK; +} + +res_T +sgf_integrate + (struct sgf_device* dev, + struct ssp_rng* rng, + const size_t primitive_id, + const size_t steps_count) +{ + size_t nprims; + size_t i; + res_T res = RES_OK; + + if(!dev || !steps_count || !rng) { + res = RES_BAD_ARG; + goto error; + } + if(!dev->is_integrating) { + if(dev->verbose) { + logger_print(dev->logger, LOG_ERROR, +"The sgf_integrate function can be called only if sgf_begin_integrate is active" +"on the SGF device.\n"); } + res = RES_BAD_OP; goto error; } - FOR_EACH(i, 0, nprims) { - accbuf[i] = darray_gfacc_cdata_get((struct darray_gfacc*)status.value)[i]; + S3D(scene_primitives_count(dev->scn_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; } + + darray_size_t_data_get(&dev->nsteps)[primitive_id] += steps_count; + FOR_EACH(i, 0, steps_count) { + res = gebhart_radiative_path(dev, primitive_id, rng); + ASSERT(res == RES_OK); + } + exit: - if(accum) SMC(accumulator_ref_put(accum)); return res; error: goto exit; } +res_T +sgf_get + (struct sgf_device* dev, + const size_t ifrom, + const size_t ito, + struct sgf_estimator* estimator) +{ + uint64_t key; + struct accum* acc; + size_t nprims; + + if(!dev || !estimator) + return RES_BAD_ARG; + + if(!dev->is_integrating) { + if(dev->verbose) { + logger_print(dev->logger, LOG_ERROR, +"The sgf_get function can be called only if sgf_begin_integrate is active" +"on the SGF device.\n"); + } + return RES_BAD_OP; + } + + S3D(scene_primitives_count(dev->scn_desc.scene, &nprims)); + if(ifrom >= nprims || ito >= nprims) { + if(dev->verbose) { + logger_print(dev->logger, LOG_ERROR, + "Invalid primitive index `%lu'\n", + (unsigned long)(ifrom >= nprims ? ifrom : ito)); + } + return RES_BAD_ARG; + } + + key = CELL_KEY(ifrom, ito); + acc = htable_accum_find(&dev->matrix, &key); + if(!acc) { + estimator->E = estimator->V = estimator->SE = 0.0; + estimator->nsteps = 0; + } else { + size_t nsteps = darray_size_t_data_get(&dev->nsteps)[ifrom]; + estimator->E = acc->radiative_flux / (double)nsteps; + estimator->V = + acc->sqr_radiative_flux / (double)nsteps + - (acc->radiative_flux * acc->radiative_flux) / (double)(nsteps * nsteps); + estimator->SE = sqrt(estimator->V / (double)nsteps); + estimator->nsteps = nsteps; + } + return RES_OK; +} diff --git a/src/sgf.h b/src/sgf.h @@ -40,9 +40,20 @@ #define SGF_API extern IMPORT_SYM #endif +/* Helper macro that asserts if the invocation of the sgf function `Func' + * returns an error. One should use this macro on sgf function calls for which + * no explicit error checking is performed */ +#ifndef NDEBUG + #define SGF(Func) ASSERT(sgf_ ## Func == RES_OK) +#else + #define SGF(Func) sgf_ ## Func +#endif + /* Forward declaration of external types */ -struct smc_device; +struct logger; +struct mem_allocator; struct s3d_scene; +struct ssp_rng; enum sgf_material_property { SGF_MATERIAL_EMISSIVITY, @@ -50,8 +61,8 @@ enum sgf_material_property { SGF_MATERIAL_SPECULARITY }; -/* Descriptor of a gebhart factor integration */ -struct sgf_context { +/* Descriptor of the scene */ +struct sgf_scene_desc { /* Retrieve material properties from client side memory */ double (*get_material_property) (void* material, /* Client side material */ @@ -59,40 +70,64 @@ struct sgf_context { const size_t primitive_id, /* Triangle identifier */ const size_t spectral_band_id); /* Spectral band identifier */ void* material; /* Client side material */ - - size_t primitive_id; /* Triangle id on which the integration is performed */ size_t spectral_bands_count; /* Total number of spectral bands */ - struct s3d_scene* scene; /* Star-3D encapsulation of the geometry */ - struct logger* logger; /* Message logger */ - int verbose; /* Make the integration more verbose */ }; -static const struct sgf_context SGF_CONTEXT_NULL = { - NULL, NULL, 0, 0, NULL, NULL, 0 +static const struct sgf_scene_desc SGF_SCENE_DESC_NULL = { + NULL, NULL, 0, 0 }; -/* Integrated value */ -struct sgf_accum { - double radiative_flux; /* in W/m^2 */ - double sqr_radiative_flux; /* radiative_flux * radiative_flux */ +struct sgf_estimator { + double E; /* Expected value */ + double V; /* Variance */ + double SE; /* Standard error */ + size_t nsteps; /* #realisations */ }; +struct sgf_device; + BEGIN_DECLS -/* Integrate the Gebhart Factor. The output `accum_buffer' stores the - * accumulated radiative flux emitted from `ctx->primitive_id' to each - * primitive and for each spectral band. The size of accum_buffer must be at - * least equal to ctx->primitives_count * ctx->spectral_bands_count. The - * spectral bands of a primitive are stored sequentially. For instance, let the - * face F and N spectral bands, the accumulated radiative flux of F are saved - * in accum_buffer[F*N .. F*N + N). */ +SGF_API res_T +sgf_device_create + (struct logger* logger, /* May be NULL <=> use default logger */ + struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ + const int verbose, /* Verbosity level */ + struct sgf_device** dev); + +SGF_API res_T +sgf_device_ref_get + (struct sgf_device* dev); + +SGF_API res_T +sgf_device_ref_put + (struct sgf_device* dev); + +SGF_API res_T +sgf_begin_integrate + (struct sgf_device* dev, + struct sgf_scene_desc* scene); + +SGF_API res_T +sgf_end_integrate + (struct sgf_device* dev); + +/* Can be called only if sgf_begin_integrate is active on `dev' */ SGF_API res_T sgf_integrate - (struct smc_device* smc, - struct sgf_context* context, - const size_t steps_count, - struct sgf_accum* accum_buffer); + (struct sgf_device* dev, + struct ssp_rng* rng, + const size_t primitive_id, + const size_t steps_count); + +/* Can be called only if sgf_begin_integrate is active on `dev' */ +SGF_API res_T +sgf_get + (struct sgf_device* dev, + const size_t primitive_id_from, + const size_t primitive_id_to, + struct sgf_estimator* estimator); END_DECLS diff --git a/src/sgf_accum.h b/src/sgf_accum.h @@ -1,122 +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. */ - -#ifndef SGF_ACCUM_H -#define SGF_ACCUM_H - -#include <rsys/dynamic_array.h> -#include <rsys/logger.h> - -#include <star/s3d.h> -#include <star/smc.h> - -#define DARRAY_DATA struct sgf_accum -#define DARRAY_NAME gfacc -#include <rsys/dynamic_array.h> - -static void* -accum_buffer_create(struct mem_allocator* allocator, void* context) -{ - struct sgf_context* ctx = context; - struct darray_gfacc* buf = NULL; - size_t nprims; - res_T res = RES_OK; - ASSERT(allocator && ctx); - - buf = MEM_ALLOC(allocator, sizeof(struct darray_gfacc)); - if(!buf) goto error; - - S3D(scene_primitives_count(ctx->scene, &nprims)); - - darray_gfacc_init(allocator, buf); - res = darray_gfacc_resize(buf, nprims); - if(res != RES_OK) goto error; - -exit: - return buf; -error: - if(ctx->verbose) { - logger_print(ctx->logger, LOG_ERROR, - "Couldn't allocate the gebhart factor accum buffer\n"); - } - if(buf) { - MEM_RM(allocator, buf); - buf = NULL; - } - goto exit; -} - -static void -accum_buffer_destroy(struct mem_allocator* allocator, void* data) -{ - struct darray_gfacc* buf = data; - ASSERT(allocator && data); - darray_gfacc_release(buf); - MEM_RM(allocator, buf); -} - -static void -accum_buffer_clear(void* data) -{ - struct darray_gfacc* buf = data; - ASSERT(buf); - memset(darray_gfacc_data_get(buf), 0, - darray_gfacc_size_get(buf) * sizeof(struct sgf_accum)); -} - -static void -accum_buffer_add(void* result, const void* op0, const void* op1) -{ - struct darray_gfacc* dst = result; - const struct darray_gfacc* a = op0; - const struct darray_gfacc* b = op1; - const size_t ngebaccs = darray_gfacc_size_get(dst); - size_t i; - ASSERT(result && op0 && op1); - ASSERT(ngebaccs==darray_gfacc_size_get(a)); - ASSERT(ngebaccs==darray_gfacc_size_get(b)); - - FOR_EACH(i, 0, ngebaccs) { - darray_gfacc_data_get(dst)[i].radiative_flux = - darray_gfacc_cdata_get(a)[i].radiative_flux - + darray_gfacc_cdata_get(b)[i].radiative_flux; - darray_gfacc_data_get(dst)[i].sqr_radiative_flux = - darray_gfacc_cdata_get(a)[i].sqr_radiative_flux - + darray_gfacc_cdata_get(b)[i].sqr_radiative_flux; - } -} - -static const struct smc_type_accum accum_buffer_type = { - accum_buffer_create, - accum_buffer_destroy, - accum_buffer_clear, - accum_buffer_add -}; - -#endif /* SGF_ACCUM_H */ - diff --git a/src/test_sgf_cube.c b/src/test_sgf_cube.c @@ -33,7 +33,7 @@ #include <rsys/stretchy_array.h> #include <star/s3d.h> -#include <star/smc.h> +#include <star/ssp.h> #define NSTEPS 10000L @@ -99,20 +99,21 @@ main(int argc, char** argv) struct s3d_vertex_data attribs[2]; struct triangle_mesh mesh; struct material mtr; - struct smc_device* smc; - struct sgf_context ctx = SGF_CONTEXT_NULL; - struct sgf_accum* gfacc = NULL; + struct sgf_device* sgf = NULL; + struct sgf_scene_desc scn_desc = SGF_SCENE_DESC_NULL; + struct sgf_estimator* estimators = NULL; + struct ssp_rng* rng; size_t iprim; - size_t nsteps_adjusted; (void)argc, (void)argv; mem_init_proxy_allocator(&allocator, &mem_default_allocator); - CHECK(s3d_device_create(NULL, &allocator, 1, &s3d), RES_OK); + CHECK(s3d_device_create(NULL, &allocator, 0, &s3d), RES_OK); CHECK(s3d_shape_create_mesh(s3d, &shape), RES_OK); CHECK(s3d_scene_create(s3d, &scn), RES_OK); CHECK(s3d_scene_attach_shape(scn, shape), RES_OK); - CHECK(smc_device_create(NULL, NULL, 1, NULL, &smc), RES_OK); + CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK); + CHECK(sgf_device_create(NULL, NULL, 1, &sgf), RES_OK); attribs[0].type = S3D_FLOAT3; attribs[0].usage = S3D_POSITION; @@ -130,118 +131,107 @@ main(int argc, char** argv) CHECK(s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, get_ids, (unsigned)nvertices, attribs, &mesh), RES_OK); - ctx.get_material_property = get_material_property; - ctx.material = &mtr; - ctx.spectral_bands_count = 1; - ctx.scene = scn; - ctx.logger = LOGGER_DEFAULT; - ctx.verbose = 1; - - gfacc = sa_add(gfacc, nprims*nprims); - - CHECK(sgf_integrate(NULL, NULL, 0, NULL), RES_BAD_ARG); - CHECK(sgf_integrate(smc, NULL, 0, NULL), RES_BAD_ARG); - CHECK(sgf_integrate(NULL, &ctx, 0, NULL), RES_BAD_ARG); - CHECK(sgf_integrate(smc, &ctx, 0, NULL), RES_BAD_ARG); - CHECK(sgf_integrate(NULL, NULL, NSTEPS, NULL), RES_BAD_ARG); - CHECK(sgf_integrate(smc, NULL, NSTEPS, NULL), RES_BAD_ARG); - CHECK(sgf_integrate(NULL, &ctx, NSTEPS, NULL), RES_BAD_ARG); - CHECK(sgf_integrate(smc, &ctx, NSTEPS, NULL), RES_BAD_ARG); - CHECK(sgf_integrate(NULL, NULL, 0, gfacc), RES_BAD_ARG); - CHECK(sgf_integrate(smc, NULL, 0, gfacc), RES_BAD_ARG); - CHECK(sgf_integrate(NULL, &ctx, 0, gfacc), RES_BAD_ARG); - CHECK(sgf_integrate(smc, &ctx, 0, gfacc), RES_BAD_ARG); - CHECK(sgf_integrate(NULL, NULL, NSTEPS, gfacc), RES_BAD_ARG); - CHECK(sgf_integrate(smc, NULL, NSTEPS, gfacc), RES_BAD_ARG); - CHECK(sgf_integrate(NULL, &ctx, NSTEPS, gfacc), RES_BAD_ARG); - CHECK(sgf_integrate(smc, &ctx, NSTEPS, gfacc), RES_BAD_OP); - - CHECK(s3d_scene_begin_session(scn, S3D_TRACE|S3D_GET_PRIMITIVE), RES_OK); - FOR_EACH(iprim, 0, nprims) { - struct sgf_accum* row = gfacc + iprim*nprims; - ctx.primitive_id = iprim; - CHECK(sgf_integrate(smc, &ctx, NSTEPS, row), RES_OK); - } - CHECK(s3d_scene_end_session(scn), RES_OK); - -#if 0 - FOR_EACH(iprim, 0, nprims) { - const struct sgf_accum* row = gfacc + iprim * nprims; - unsigned icol; - FOR_EACH(icol, 0, nprims) { - const double E = row[icol].radiative_flux / NSTEPS; - printf("%.6f ", E); + scn_desc.get_material_property = get_material_property; + scn_desc.material = &mtr; + scn_desc.spectral_bands_count = 1; + scn_desc.scene = scn; + + estimators = sa_add(estimators, nprims * nprims); + + CHECK(sgf_begin_integrate(NULL, NULL), RES_BAD_ARG); + CHECK(sgf_begin_integrate(sgf, NULL), RES_BAD_ARG); + CHECK(sgf_begin_integrate(NULL, &scn_desc), RES_BAD_ARG); + CHECK(sgf_begin_integrate(sgf, &scn_desc), RES_OK); + CHECK(sgf_begin_integrate(sgf, &scn_desc), RES_BAD_OP); + + CHECK(sgf_integrate(NULL, NULL, SIZE_MAX, 0), RES_BAD_ARG); + CHECK(sgf_integrate(sgf, NULL, SIZE_MAX, 0), RES_BAD_ARG); + CHECK(sgf_integrate(NULL, rng, SIZE_MAX, 0), RES_BAD_ARG); + CHECK(sgf_integrate(sgf, rng, SIZE_MAX, 0), RES_BAD_ARG); + CHECK(sgf_integrate(NULL, NULL, 0, 0), RES_BAD_ARG); + CHECK(sgf_integrate(sgf, NULL, 0, 0), RES_BAD_ARG); + CHECK(sgf_integrate(NULL, rng, 0, 0), RES_BAD_ARG); + CHECK(sgf_integrate(sgf, rng, 0, 0), RES_BAD_ARG); + CHECK(sgf_integrate(NULL, NULL, SIZE_MAX, NSTEPS), RES_BAD_ARG); + CHECK(sgf_integrate(sgf, NULL, SIZE_MAX, NSTEPS), RES_BAD_ARG); + CHECK(sgf_integrate(NULL, rng, SIZE_MAX, NSTEPS), RES_BAD_ARG); + CHECK(sgf_integrate(sgf, rng, SIZE_MAX, NSTEPS), RES_BAD_ARG); + CHECK(sgf_integrate(NULL, NULL, 0, NSTEPS), RES_BAD_ARG); + CHECK(sgf_integrate(sgf, NULL, 0, NSTEPS), RES_BAD_ARG); + CHECK(sgf_integrate(NULL, rng, 0, NSTEPS), RES_BAD_ARG); + CHECK(sgf_integrate(sgf, rng, 0, NSTEPS), RES_OK); + + CHECK(sgf_get(NULL, SIZE_MAX, SIZE_MAX, NULL), RES_BAD_ARG); + CHECK(sgf_get(sgf, SIZE_MAX, SIZE_MAX, NULL), RES_BAD_ARG); + CHECK(sgf_get(NULL, 0, SIZE_MAX, NULL), RES_BAD_ARG); + CHECK(sgf_get(sgf, 0, SIZE_MAX, NULL), RES_BAD_ARG); + CHECK(sgf_get(NULL, SIZE_MAX, 0, NULL), RES_BAD_ARG); + CHECK(sgf_get(sgf, SIZE_MAX, 0, NULL), RES_BAD_ARG); + CHECK(sgf_get(NULL, 0, 0, NULL), RES_BAD_ARG); + CHECK(sgf_get(sgf, 0, 0, NULL), RES_BAD_ARG); + CHECK(sgf_get(NULL, SIZE_MAX, SIZE_MAX, estimators + 0), RES_BAD_ARG); + CHECK(sgf_get(sgf, SIZE_MAX, SIZE_MAX, estimators + 0), RES_BAD_ARG); + CHECK(sgf_get(NULL, 0, SIZE_MAX, estimators + 0), RES_BAD_ARG); + CHECK(sgf_get(sgf, 0, SIZE_MAX, estimators + 0), RES_BAD_ARG); + CHECK(sgf_get(NULL, SIZE_MAX, 0, estimators + 0), RES_BAD_ARG); + CHECK(sgf_get(sgf, SIZE_MAX, 0, estimators + 0), RES_BAD_ARG); + CHECK(sgf_get(NULL, 0, 0, estimators + 0), RES_BAD_ARG); + CHECK(sgf_get(sgf, 0, 0, estimators + 0), RES_OK); + + CHECK(estimators[0].nsteps, NSTEPS); + + CHECK(sgf_end_integrate(NULL), RES_BAD_ARG); + CHECK(sgf_end_integrate(sgf), RES_OK); + CHECK(sgf_end_integrate(sgf), RES_BAD_OP); + + CHECK(sgf_integrate(sgf, rng, 0, NSTEPS), RES_BAD_OP); + CHECK(sgf_get(sgf, 0, 0, estimators + 0), RES_BAD_OP); + + /* Integrate the Gebhart Factors */ + CHECK(sgf_begin_integrate(sgf, &scn_desc), RES_OK); + +/* #pragma omp parallel for */ + for(iprim = 0; iprim < nprims; ++iprim) { + size_t iprim2; + struct sgf_estimator* row = estimators + iprim * nprims; + + CHECK(sgf_integrate(sgf, rng, iprim, NSTEPS), RES_OK); + + FOR_EACH(iprim2, 0, nprims) { + CHECK(sgf_get(sgf, iprim, iprim2, row + iprim2), RES_OK); + CHECK(row[iprim2].nsteps, NSTEPS); } - printf("\n"); } - printf("\n"); -#endif + CHECK(sgf_end_integrate(sgf), RES_OK); /* Merge the radiative flux of coplanar primitives */ FOR_EACH(iprim, 0, nprims/2) { - const struct sgf_accum* row_src0 = gfacc + iprim * 2 * nprims; - const struct sgf_accum* row_src1 = row_src0 + nprims; - struct sgf_accum* row_dst = gfacc + iprim * nprims; + const struct sgf_estimator* row_src0 = estimators + iprim * 2 * nprims; + const struct sgf_estimator* row_src1 = row_src0 + nprims; size_t icol; + double sum = 0; FOR_EACH(icol, 0, nprims/2) { - const struct sgf_accum* src0 = row_src0 + icol * 2; - const struct sgf_accum* src1 = src0 + 1; - const struct sgf_accum* src2 = row_src1 + icol * 2; - const struct sgf_accum* src3 = src2 + 1; - struct sgf_accum* dst = row_dst + icol; - - dst->radiative_flux = - src0->radiative_flux - + src1->radiative_flux - + src2->radiative_flux - + src3->radiative_flux; - - dst->sqr_radiative_flux = - src0->sqr_radiative_flux - + src1->sqr_radiative_flux - + src2->sqr_radiative_flux - + src3->sqr_radiative_flux; - } - } - nsteps_adjusted = NSTEPS * 2; + const struct sgf_estimator* src0 = row_src0 + icol * 2; + const struct sgf_estimator* src1 = src0 + 1; + const struct sgf_estimator* src2 = row_src1 + icol * 2; + const struct sgf_estimator* src3 = src2 + 1; + double E = (src0->E + src1->E + src2->E + src3->E) / 2; - /* Estimated value */ - FOR_EACH(iprim, 0, nprims/2) { - const struct sgf_accum* row = gfacc + iprim * nprims; - size_t icol; - double sum = 0.0; - FOR_EACH(icol, 0, nprims/2) { - const double E = row[icol].radiative_flux / ((double)nsteps_adjusted); sum += E; printf("%.6f ", E); } printf("\n"); - /* Ensure the energy conservation property */ - CHECK(eq_eps(sum, 1.0, 1.e-6), 1); - } - printf("\n"); - - /* Standard error */ - FOR_EACH(iprim, 0, nprims/2) { - const struct sgf_accum* row = gfacc + iprim * nprims; - size_t icol; - FOR_EACH(icol, 0, nprims/2) { - const double V = row[icol].sqr_radiative_flux / ((double)nsteps_adjusted) - - (row[icol].radiative_flux*row[icol].radiative_flux) - / (double)(nsteps_adjusted*nsteps_adjusted); - const double SE = sqrt(V / NSTEPS); - printf("%.6f ", SE); - } - printf("\n"); + CHECK(eq_eps(sum, 1.0, 1.e-6), 1); /* Ensure the energy conservation */ } - sa_release(gfacc); + sa_release(estimators); CHECK(s3d_shape_ref_put(shape), RES_OK); CHECK(s3d_scene_ref_put(scn), RES_OK); CHECK(s3d_device_ref_put(s3d), RES_OK); - CHECK(smc_device_ref_put(smc), RES_OK); + CHECK(ssp_rng_ref_put(rng), RES_OK); + CHECK(sgf_device_ref_put(sgf), RES_OK); check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); diff --git a/src/test_sgf_device.c b/src/test_sgf_device.c @@ -0,0 +1,86 @@ +/* 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 "test_sgf_utils.h" + +#include <rsys/logger.h> +#include <rsys/rsys.h> + +static void +log_stream(const char* msg, void* ctx) +{ + ASSERT(msg); + (void)msg, (void)ctx; + printf("%s\n", msg); +} + +int +main(int argc, char** argv) +{ + struct logger logger; + struct mem_allocator allocator; + struct sgf_device* dev; + (void)argc, (void)argv; + + CHECK(sgf_device_create(NULL, NULL, 0, NULL), RES_BAD_ARG); + CHECK(sgf_device_create(NULL, NULL, 0, &dev), RES_OK); + + CHECK(sgf_device_ref_get(NULL), RES_BAD_ARG); + CHECK(sgf_device_ref_get(dev), RES_OK); + CHECK(sgf_device_ref_put(NULL), RES_BAD_ARG); + CHECK(sgf_device_ref_put(dev), RES_OK); + CHECK(sgf_device_ref_put(dev), RES_OK); + + mem_init_proxy_allocator(&allocator, &mem_default_allocator); + + CHECK(MEM_ALLOCATED_SIZE(&allocator), 0); + CHECK(sgf_device_create(NULL, &allocator, 1, NULL), RES_BAD_ARG); + CHECK(sgf_device_create(NULL, &allocator, 1, &dev), RES_OK); + CHECK(sgf_device_ref_put(dev), RES_OK); + CHECK(MEM_ALLOCATED_SIZE(&allocator), 0); + + CHECK(logger_init(&allocator, &logger), RES_OK); + logger_set_stream(&logger, LOG_OUTPUT, log_stream, NULL); + logger_set_stream(&logger, LOG_ERROR, log_stream, NULL); + logger_set_stream(&logger, LOG_WARNING, log_stream, NULL); + + CHECK(sgf_device_create(&logger, NULL, 0, NULL), RES_BAD_ARG); + CHECK(sgf_device_create(&logger, NULL, 0, &dev), RES_OK); + CHECK(sgf_device_ref_put(dev), RES_OK); + + CHECK(sgf_device_create(&logger, &allocator, 4, NULL), RES_BAD_ARG); + CHECK(sgf_device_create(&logger, &allocator, 4, &dev), RES_OK); + CHECK(sgf_device_ref_put(dev), RES_OK); + + logger_release(&logger); + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHECK(mem_allocated_size(), 0); + return 0; +} diff --git a/src/test_sgf_utils.h b/src/test_sgf_utils.h @@ -44,7 +44,7 @@ struct material { const double* specularity; }; -static void +static INLINE void get_ids(const unsigned itri, unsigned ids[3], void* data) { const struct triangle_mesh* mesh = data; @@ -57,7 +57,7 @@ get_ids(const unsigned itri, unsigned ids[3], void* data) ids[2] = mesh->indices[id + 2]; } -static void +static INLINE void get_pos(const unsigned ivert, float pos[3], void* data) { const struct triangle_mesh* mesh = data; @@ -70,7 +70,7 @@ get_pos(const unsigned ivert, float pos[3], void* data) pos[2] = mesh->vertices[i + 2]; } -static double +static INLINE double get_material_property (void* material, const enum sgf_material_property prop,