commit 10fbe01e1892bb1ad8c9f0f96caed8a11930ffee
parent 0dfae75c6705adbb51caba57462236719996e1ee
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 6 Jul 2016 17:33:00 +0200
Add support of Gebhart Factor integration on 2D scene
Make the gebhart_factor realisation generic to the dimensionality of the
submitted scene and instantiate it for both 2D and 3D scene. The 2D
scenes are handled through the Star-2D library.
Diffstat:
6 files changed, 401 insertions(+), 167 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -27,6 +27,7 @@ find_package(RCMake 0.2.3 REQUIRED)
find_package(RSys 0.3 REQUIRED)
find_package(StarSP 0.3 REQUIRED)
find_package(Star3D 0.3 REQUIRED)
+find_package(Star2D REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
include_directories(
@@ -48,7 +49,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
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_INC sgf_device_c.h sgf_realisation.h)
set(SGF_FILES_DOC COPYING README.md)
# Prepend each file in the `SGF_FILES_<SRC|INC>' list by `SGF_SOURCE_DIR'
@@ -63,7 +64,7 @@ set_target_properties(sgf PROPERTIES
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
-target_link_libraries(sgf RSys Star3D StarSP)
+target_link_libraries(sgf RSys Star2D Star3D StarSP)
if(CMAKE_COMPILER_IS_GNUCC)
target_link_libraries(sgf m)
endif()
diff --git a/src/sgf.h b/src/sgf.h
@@ -48,6 +48,11 @@ enum sgf_material_property {
SGF_MATERIAL_SPECULARITY
};
+enum sgf_dimensionality {
+ SGF_2D,
+ SGF_3D
+};
+
/* Descriptor of the scene */
struct sgf_scene_desc {
/* Retrieve material properties from client side memory */
@@ -58,11 +63,17 @@ struct sgf_scene_desc {
const size_t spectral_band_id); /* Spectral band identifier */
void* material; /* Client side material */
size_t spectral_bands_count; /* Total number of spectral bands */
- struct s3d_scene* scene; /* Star-3D encapsulation of the geometry */
+
+ /* Star-3D encapsulation of the geometry */
+ enum sgf_dimensionality dimensionality;
+ union {
+ struct s3d_scene* scn3d;
+ struct s2d_scene* scn2d;
+ } geometry;
};
static const struct sgf_scene_desc SGF_SCENE_DESC_NULL = {
- NULL, NULL, 0, 0
+ NULL, NULL, 0, SGF_3D, { NULL }
};
/* Estimated Gebart Factor between 2 faces */
diff --git a/src/sgf_estimator.c b/src/sgf_estimator.c
@@ -17,6 +17,7 @@
#include "sgf.h"
#include "sgf_device_c.h"
+#include "sgf_realisation.h"
#include <star/s3d.h>
#include <star/ssp.h>
@@ -31,18 +32,18 @@ enum status_type {
INFINITY_RADIATIVE_FLUX
};
-struct accum {
- double radiative_flux;
- double sqr_radiative_flux;
-};
-
+/* Generate the accum dynamic array data type */
#define DARRAY_NAME accum
#define DARRAY_DATA struct accum
#include <rsys/dynamic_array.h>
-/* How many self intersections are tested on the same ray before an error
- * occurs */
-#define NSELF_HITS_MAX 10
+/* Generate the 2D realisation function */
+#define SGF_DIMENSIONALITY 2
+#include "sgf_realisation.h"
+
+/* Generate the 3D realisation function */
+#define SGF_DIMENSIONALITY 3
+#include "sgf_realisation.h"
/* Estimator of the Gebhart Factor of a given face to all the other ones */
struct sgf_estimator {
@@ -60,19 +61,16 @@ struct sgf_estimator {
/*******************************************************************************
* Helper functions
******************************************************************************/
-static INLINE char
+static INLINE int
check_scene_desc(struct sgf_scene_desc* desc)
{
ASSERT(desc);
- return desc->get_material_property && desc->scene;
-}
-
-static FINLINE void
-accum_weight(struct accum* accums, const size_t iface, const double weight)
-{
- ASSERT(accums);
- accums[iface].radiative_flux += weight;
- accums[iface].sqr_radiative_flux += weight * weight;
+ if(!desc->get_material_property) return 0;
+ switch(desc->dimensionality) {
+ case SGF_2D: return desc->geometry.scn2d != NULL;
+ case SGF_3D: return desc->geometry.scn3d != NULL;
+ default: return 0;
+ }
}
static FINLINE void
@@ -91,144 +89,6 @@ setup_status
status->nsteps = nsteps; /* # realisations */
}
-/* Return RES_BAD_OP on numerical issue. i.e. the radiative path is skipped */
-static res_T
-gebhart_radiative_path
- (struct sgf_device* dev,
- struct accum* accums,
- struct accum* accums_infinity,
- struct ssp_rng* rng,
- const float epsilon,
- const size_t primitive_id,
- struct sgf_scene_desc* desc)
-{
- struct s3d_hit hit;
- struct s3d_primitive prim_from, prim;
- struct s3d_attrib attrib;
- size_t nself_hits = 0;
- const double trans_min = 1.e-8; /* Minimum transmissivity threshold */
- double proba_reflec_spec;
- double transmissivity, emissivity, reflectivity, specularity;
-#ifndef NDEBUG
- double infinite_radiative_flux = 0.0;
- 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 */
- float u, v;
-
- ASSERT(accums && accums_infinity && epsilon > 0.f && rng && desc);
-
- /* 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. */
- u = ssp_rng_canonical_float(rng);
- v = ssp_rng_canonical_float(rng);
- S3D(primitive_sample(&prim, u, v, 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;
- for(;;) { /* Here we go */
- prim_from = prim;
- range[0] = epsilon, range[1] = FLT_MAX;
-
- nself_hits = 0;
- do { /* Handle auto intersection */
- S3D(scene_trace_ray(desc->scene, pos, dir, range, NULL, &hit));
- if(!S3D_HIT_NONE(&hit)) range[0] = MMAX(range[0], hit.distance);
- range[0] = nextafterf(range[0], range[1]);
- } while(S3D_PRIMITIVE_EQ(&hit.prim, &prim_from) /* Self hit */
- && ++nself_hits < NSELF_HITS_MAX); /* Too many self hits */
-
- if(nself_hits >= NSELF_HITS_MAX) {
- if(dev->verbose) {
- logger_print(dev->logger, LOG_ERROR,
- "Too many self hit on a given ray. Ray origin: %g, %g, %g\n",
- SPLIT3(pos));
- }
- return RES_BAD_OP;
- }
-
- if(S3D_HIT_NONE(&hit)) { /* The ray is outside the volume */
- accum_weight(accums_infinity, prim.scene_prim_id, transmissivity);
-#ifndef NDEBUG
- infinite_radiative_flux = transmissivity;
-#endif
- 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(accums, 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(accums, 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);
- 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 + infinite_radiative_flux, 1.0, 1.e6));
-#endif
- return RES_OK;
-}
-
static res_T
estimator_create(struct sgf_device* dev, struct sgf_estimator** out_estimator)
{
@@ -334,14 +194,18 @@ sgf_integrate
const size_t steps_count,
struct sgf_estimator** out_estimator)
{
+ #define SXD_FUNC(Func) (desc->dimensionality == SGF_2D ? S2D(Func) : S3D(Func))
+ #define SXD_ENUM(Enum) (desc->dimensionality == SGF_2D ? S2D_##Enum : S3D_##Enum)
struct sgf_estimator* estimator = NULL;
size_t istep;
size_t ispectral_band;
size_t nprims;
float epsilon;
float lower[3], upper[3];
+ void* scene;
int mask;
res_T res = RES_OK;
+ gebhart_radiative_path_T gebhart_radiative_path;
if(!dev || !steps_count || !rng || !desc || !check_scene_desc(desc)
|| !out_estimator) {
@@ -352,19 +216,32 @@ sgf_integrate
res = estimator_create(dev, &estimator);
if(res != RES_OK) goto error;
+ switch(desc->dimensionality) {
+ case SGF_2D:
+ scene = desc->geometry.scn2d;
+ gebhart_radiative_path = gebhart_radiative_path_2d;
+ break;
+ case SGF_3D:
+ scene = desc->geometry.scn3d;
+ gebhart_radiative_path = gebhart_radiative_path_3d;
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
/* Check scene active sessions */
- S3D(scene_get_session_mask(desc->scene, &mask));
- if(!(mask & S3D_TRACE) || !(mask & S3D_GET_PRIMITIVE)) {
+ SXD_FUNC(scene_get_session_mask(scene, &mask));
+ if(!(mask & SXD_ENUM(TRACE)) || !(mask & SXD_ENUM(GET_PRIMITIVE))) {
if(dev->verbose) {
logger_print(dev->logger, LOG_ERROR,
- "No active trace/get_primitive session on the Star-3D scene\n");
+ "No active trace/get_primitive session on the Star-%dD scene\n",
+ desc->dimensionality == SGF_2D ? 2 : 3);
}
res = RES_BAD_ARG;
goto error;
}
/* Check submitted primitive_id */
- S3D(scene_primitives_count(desc->scene, &nprims));
+ SXD_FUNC(scene_primitives_count(scene, &nprims));
if(primitive_id >= nprims) {
if(dev->verbose) {
logger_print(dev->logger, LOG_ERROR,
@@ -375,7 +252,7 @@ sgf_integrate
}
/* Compute an epsilon relatively to the scene size */
- S3D(scene_get_aabb(desc->scene, lower, upper));
+ SXD_FUNC(scene_get_aabb(scene, lower, upper));
epsilon = upper[0] - lower[0];
epsilon = MMIN(epsilon, upper[1] - lower[1]);
epsilon = MMIN(epsilon, upper[2] - lower[2]);
@@ -439,6 +316,9 @@ exit:
error:
if(estimator) SGF(estimator_ref_put(estimator));
goto exit;
+
+ #undef SXD_FUNC
+ #undef SXD_ENUM
}
res_T
diff --git a/src/sgf_realisation.h b/src/sgf_realisation.h
@@ -0,0 +1,340 @@
+/* Copyright (C) 2015-2016 EDF S.A., France (syrthes-support@edf.fr)
+ *
+ * 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/>. */
+
+/*
+ * Generate the gebhart_radiative_path_<2|3>d realisation functions with
+ * respect to the SGF_DIMENSIONALITY macro
+ */
+
+#ifndef SGF_DIMENSIONALITY
+
+#ifndef SGF_REALISATION_H
+#define SGF_REALISATION_H
+
+#include <star/s2d.h>
+#include <star/s3d.h>
+#include <star/ssp.h>
+#include <rsys/float2.h>
+
+/* How many self intersections are tested on the same ray before an error
+ * occurs */
+#define NSELF_HITS_MAX 10
+
+struct accum {
+ double radiative_flux;
+ double sqr_radiative_flux;
+};
+
+/* Type of the realisation function. Return RES_BAD_OP on numerical issue. i.e.
+ * the radiative path is skipped */
+typedef res_T
+(*gebhart_radiative_path_T)
+ (struct sgf_device* dev,
+ struct accum* accums,
+ struct accum* accums_infinity,
+ struct ssp_rng* rng,
+ const float epsilon,
+ const size_t primitive_id,
+ struct sgf_scene_desc* desc);
+
+static FINLINE void
+accum_weight(struct accum* accums, const size_t iface, const double weight)
+{
+ ASSERT(accums);
+ accums[iface].radiative_flux += weight;
+ accums[iface].sqr_radiative_flux += weight * weight;
+}
+
+/*******************************************************************************
+ * 2D helper functions
+ ******************************************************************************/
+static FINLINE void
+primitive2d_get_normal(struct s2d_primitive* prim, float normal[3])
+{
+ struct s2d_attrib attr;
+ const float s = 0;
+ ASSERT(prim && normal);
+ S2D(primitive_get_attrib(prim, S2D_GEOMETRY_NORMAL, s, &attr));
+ ASSERT(attr.type == S2D_FLOAT2);
+ f2_normalize(normal, attr.value);
+ normal[2] = 0.f;
+}
+
+static FINLINE void
+primitive2d_sample_position
+ (struct s2d_primitive* prim,
+ struct ssp_rng* rng,
+ float position[3])
+{
+ struct s2d_attrib attr;
+ float s;
+ ASSERT(prim && position);
+ S2D(primitive_sample(prim, ssp_rng_canonical_float(rng), &s));
+ S2D(primitive_get_attrib(prim, S2D_POSITION, s, &attr));
+ ASSERT(attr.type == S2D_FLOAT2);
+ f2_set(position, attr.value);
+ position[2] = 0.f;
+}
+
+static FINLINE void
+hit2d_get_normal(struct s2d_hit* hit, float normal[3])
+{
+ ASSERT(hit && normal);
+ f2_normalize(normal, hit->normal);
+ normal[2] = 0.f;
+}
+
+/*******************************************************************************
+ * 3D helper functions
+ ******************************************************************************/
+static FINLINE void
+primitive3d_get_normal(struct s3d_primitive* prim, float normal[3])
+{
+ struct s3d_attrib attr;
+ const float st[2] = { 0.f, 0.f };
+ ASSERT(prim && normal);
+ S3D(primitive_get_attrib(prim, S3D_GEOMETRY_NORMAL, st, &attr));
+ ASSERT(attr.type == S3D_FLOAT3);
+ f3_normalize(normal, attr.value);
+}
+
+static FINLINE void
+primitive3d_sample_position
+ (struct s3d_primitive* prim,
+ struct ssp_rng* rng,
+ float position[3])
+{
+ struct s3d_attrib attr;
+ float st[2];
+ float u, v;
+ ASSERT(prim && position);
+
+ u = ssp_rng_canonical_float(rng);
+ v = ssp_rng_canonical_float(rng);
+ S3D(primitive_sample(prim, u, v, st));
+ S3D(primitive_get_attrib(prim, S3D_POSITION, st, &attr));
+ ASSERT(attr.type == S3D_FLOAT3);
+ f3_set(position, attr.value);
+}
+
+static FINLINE void
+hit3d_get_normal(struct s3d_hit* hit, float normal[3])
+{
+ ASSERT(hit && normal);
+ f3_normalize(normal, hit->normal);
+}
+
+#endif /* SGF_REALISATION_H */
+
+#else /* !SGF_DIMENSIONALITY */
+
+#if SGF_DIMENSIONALITY == 2
+ #define GEBHART_RADIATIVE_PATH gebhart_radiative_path_2d
+
+ /* Wrap type */
+ #define sXd_scene_T struct s2d_scene
+ #define sXd_hit_T struct s2d_hit
+ #define sXd_primitive_T struct s2d_primitive
+
+ /* Wrap macros */
+ #define SXD_HIT_NONE(Hit) S2D_HIT_NONE(Hit)
+ #define SXD_PRIMITIVE_EQ(A, B) S2D_PRIMITIVE_EQ(A, B)
+
+ /* Wrap Functions */
+ #define sXd_scene_get_primitive(Scn, IPrim, prim) \
+ S2D(scene_get_primitive(Scn, IPrim, prim))
+ #define sXd_primitive_get_normal(Prim, Normal) \
+ primitive2d_get_normal(Prim, Normal)
+ #define sXd_primitive_sample_position(Prim, Rng, Position) \
+ primitive2d_sample_position(Prim, Rng, Position)
+ #define sXd_scene_trace_ray(Scn, Org, Dir, Range, Hit) \
+ S2D(scene_trace_ray_3d(Scn, Org, Dir, Range, NULL, Hit))
+ #define sXd_hit_get_normal(Hit, Normal) \
+ hit2d_get_normal(Hit, Normal)
+ #define sgf_scene_desc_get_sXd_scene(Desc) (Desc)->geometry.scn2d
+
+#elif SGF_DIMENSIONALITY == 3
+ #define GEBHART_RADIATIVE_PATH gebhart_radiative_path_3d
+
+ /* Wrap type */
+ #define sXd_scene_T struct s3d_scene
+ #define sXd_hit_T struct s3d_hit
+ #define sXd_primitive_T struct s3d_primitive
+
+ /* Wrap macros */
+ #define SXD_HIT_NONE(Hit) S3D_HIT_NONE(Hit)
+ #define SXD_PRIMITIVE_EQ(A, B) S3D_PRIMITIVE_EQ(A, B)
+
+ /* Wrap Functions */
+ #define sXd_scene_get_primitive(Scn, IPrim, prim) \
+ S3D(scene_get_primitive(Scn, IPrim, prim))
+ #define sXd_primitive_get_normal(Prim, Normal) \
+ primitive3d_get_normal(Prim, Normal)
+ #define sXd_primitive_sample_position(Prim, Rng, Position) \
+ primitive3d_sample_position(Prim, Rng, Position)
+ #define sXd_scene_trace_ray(Scn, Org, Dir, Range, Hit) \
+ S3D(scene_trace_ray(Scn, Org, Dir, Range, NULL, Hit))
+ #define sXd_hit_get_normal(Hit, Normal) \
+ hit3d_get_normal(Hit, Normal)
+ #define sgf_scene_desc_get_sXd_scene(Desc) (Desc)->geometry.scn3d
+#else
+ #error Unexpected dimensionility
+#endif
+
+static res_T
+GEBHART_RADIATIVE_PATH
+ (struct sgf_device* dev,
+ struct accum* accums,
+ struct accum* accums_infinity,
+ struct ssp_rng* rng,
+ const float epsilon,
+ const size_t primitive_id,
+ struct sgf_scene_desc* desc)
+{
+ sXd_scene_T* scene;
+ sXd_hit_T hit;
+ sXd_primitive_T prim_from, prim;
+ size_t nself_hits = 0;
+ const double trans_min = 1.e-8; /* Minimum transmissivity threshold */
+ double proba_reflec_spec;
+ double transmissivity, emissivity, reflectivity, specularity;
+#ifndef NDEBUG
+ double infinite_radiative_flux = 0.0;
+ double sum_radiative_flux = 0.f;
+#endif
+ 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 */
+ ASSERT(accums && accums_infinity && epsilon > 0.f && rng && desc);
+
+ scene = sgf_scene_desc_get_sXd_scene(desc);
+
+ /* 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;
+
+ /* Retrieve the scene primitive */
+ sXd_scene_get_primitive(scene, (unsigned)primitive_id, &prim);
+ /* Get the geometric normal of the input primitive */
+ sXd_primitive_get_normal(&prim, normal);
+ /* Uniformly sample prim to define the origin of the radiative path */
+ sXd_primitive_sample_position(&prim, rng, pos);
+ /* Cosine weighted sampling of the radiative path direction around `normal' */
+ ssp_ran_hemisphere_cos(rng, normal, dir);
+
+ transmissivity = 1.0;
+ for(;;) { /* Here we go */
+ prim_from = prim;
+ range[0] = epsilon, range[1] = FLT_MAX;
+
+ nself_hits = 0;
+ do { /* Handle auto intersection */
+ sXd_scene_trace_ray(scene, pos, dir, range, &hit);
+ if(!SXD_HIT_NONE(&hit)) range[0] = MMAX(range[0], hit.distance);
+ range[0] = nextafterf(range[0], range[1]);
+ } while(SXD_PRIMITIVE_EQ(&hit.prim, &prim_from) /* Self hit */
+ && ++nself_hits < NSELF_HITS_MAX); /* Too many self hits */
+
+ if(nself_hits >= NSELF_HITS_MAX) {
+ if(dev->verbose) {
+ logger_print(dev->logger, LOG_ERROR,
+ "Too many self hit on a given ray. Ray origin: %g, %g, %g\n",
+ SPLIT3(pos));
+ }
+ return RES_BAD_OP;
+ }
+
+ if(SXD_HIT_NONE(&hit)) { /* The ray is outside the volume */
+ accum_weight(accums_infinity, prim.scene_prim_id, transmissivity);
+#ifndef NDEBUG
+ infinite_radiative_flux = transmissivity;
+#endif
+ 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(accums, 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(accums, prim.scene_prim_id, weight);
+#ifndef NDEBUG
+ sum_radiative_flux += weight;
+#endif
+ break;
+ }
+ }
+
+ if(reflectivity <= 0.0) break;
+
+ proba_reflec_spec = specularity / reflectivity;
+ sXd_hit_get_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);
+ 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 + infinite_radiative_flux, 1.0, 1.e6));
+#endif
+ return RES_OK;
+}
+
+#undef GEBHART_RADIATIVE_PATH
+#undef sXd_scene_T
+#undef sXd_hit_T
+#undef sXd_primitive_T
+#undef SXD_HIT_NONE
+#undef SXD_PRIMITIVE_EQ
+#undef sXd_scene_get_primitive
+#undef sXd_primitive_get_normal
+#undef sXd_primitive_sample_position
+#undef sXd_scene_trace_ray
+#undef sXd_hit_get_normal
+#undef sgf_scene_desc_get_sXd_scene
+
+#undef SGF_DIMENSIONALITY
+
+#endif /* !SGF_DIMENSIONALITY */
diff --git a/src/test_sgf_cube.c b/src/test_sgf_cube.c
@@ -127,7 +127,8 @@ main(int argc, char** argv)
scn_desc.get_material_property = get_material_property;
scn_desc.material = &mtr;
scn_desc.spectral_bands_count = 1;
- scn_desc.scene = scn;
+ scn_desc.dimensionality = SGF_3D;
+ scn_desc.geometry.scn3d = scn;
status = sa_add(status, nprims*nprims);
rngs = sa_add(rngs, nbuckets);
diff --git a/src/test_sgf_tetrahedron.c b/src/test_sgf_tetrahedron.c
@@ -94,7 +94,8 @@ main(int argc, char** argv)
desc.get_material_property = get_material_property;
desc.material = &mtr;
desc.spectral_bands_count = 1;
- desc.scene = scn;
+ desc.dimensionality = SGF_3D;
+ desc.geometry.scn3d = scn;
status = sa_add(status, nprims*nprims);
rngs = sa_add(rngs, nbuckets);
