commit c22998c6d97f1fc27a9976224b20346056ed52a0
parent edcb49d51b3fc7ec8d9f400e9b2c8994467c7014
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 23 Jul 2015 22:05:28 +0200
Fix an integration issue and add a new test
The transmissivity was updated even though the russian roulette failed.
Diffstat:
6 files changed, 581 insertions(+), 325 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -93,7 +93,8 @@ if(NOT NO_TEST)
target_link_libraries(${_name} sgf RSys)
add_test(${_name} ${_name})
endfunction(new_test)
- new_test(test_sgf)
+ new_test(test_sgf_cube)
+ new_test(test_sgf_tetrahedron)
endif()
################################################################################
diff --git a/src/sgf.c b/src/sgf.c
@@ -63,16 +63,23 @@ gebhart_radiative_path
double emissivity;
double reflectivity;
double specularity;
- const double trans_min = 1.e-6;
+ const double trans_min = 1.e-8;
size_t iprim, iprim_from;
+ size_t nbounce = 0;
float vec0[3], vec1[3];
float triangle[3][3];
float normal[3];
float pos[3];
- float dir[3];
+ float dir[4];
float range[2];
float dist;
+ /* Discard faces with no emissivity */
+ emissivity = ctx->get_material_property(ctx->material,
+ SGF_MATERIAL_EMISSIVITY, ctx->primitive_id, ctx->spectral_band_id);
+ if(emissivity <= 0.0)
+ return RES_OK;
+
/* Compute the geometric normal of the triangle */
ctx->get_triangle(ctx->geometry, ctx->primitive_id, triangle);
f3_sub(vec0, triangle[1], triangle[0]);
@@ -89,18 +96,17 @@ gebhart_radiative_path
ssp_ran_hemisphere_cos(rng, normal, dir);
iprim = ctx->primitive_id;
+ transmissivity = 1.0;
+
for(;;) { /* Here we go */
iprim_from = iprim;
- /* Ensure that the point lies on or slightly forward the triangle plane */
- range[1] = FLT_MAX;
- dist = f3_dot(pos, normal) - f3_dot(normal, triangle[0]);/*Ax + By + Cz + D*/
- for(range[0] = FLT_MIN; dist + range[0] < 0.f; range[0] += FLT_EPSILON);
-
+ range[0] = FLT_MIN, range[1] = FLT_MAX;
for(;;) {
S3D(scene_trace_ray(ctx->scene, pos, dir, range, &hit));
if(hit.prim.prim_id == iprim_from) range[0] += FLT_EPSILON; /* Self hit */
else break;
+
}
if(S3D_HIT_NONE(&hit)) {
@@ -116,6 +122,11 @@ gebhart_radiative_path
f3_add(pos, vec0, pos);
iprim = hit.prim.prim_id;
ASSERT(iprim < ctx->primitives_count);
+#ifndef NDEBUG
+ ctx->get_triangle(ctx->geometry, iprim, triangle);
+ dist = f3_dot(pos, normal) - f3_dot(normal, triangle[0]);/*Ax + By + Cz + D*/
+ ASSERT(eq_eps(dist, 0, 1.e-6f));
+#endif
/* Fetch material property */
emissivity = ctx->get_material_property
@@ -129,29 +140,37 @@ gebhart_radiative_path
const double weight = transmissivity * emissivity;
result[iprim].radiative_flux += weight;
result[iprim].sqr_radiative_flux += weight * weight;
+ transmissivity = transmissivity * (1.0 - emissivity);
} else {
/* Russian roulette */
if(ssp_rng_canonical(rng) < emissivity) {
- result[iprim].radiative_flux += transmissivity;
- result[iprim].sqr_radiative_flux += transmissivity * transmissivity;
+ const double weight = transmissivity;
+ result[iprim].radiative_flux += weight;
+ result[iprim].sqr_radiative_flux += weight * weight;
break;
}
}
- transmissivity = transmissivity * (1.0 - emissivity);
- proba_reflec_spec = specularity / reflectivity;
+ 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);
}
-
+ ++nbounce;
ctx->get_triangle(ctx->geometry, iprim, triangle);
}
-#ifndef NDEBUG
+/* ASSERT(nbounce <= 1);*/
+#if 1 && !defined(NDEBUG)
{ /* Ensure the energy conservation property */
double sum_radiative_flux = 0.0;
FOR_EACH(iprim, 0, ctx->primitives_count) {
diff --git a/src/test_sgf.c b/src/test_sgf.c
@@ -1,311 +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 <rsys/logger.h>
-#include <rsys/mem_allocator.h>
-#include <rsys/stretchy_array.h>
-
-#include <star/s3d.h>
-#include <star/smc.h>
-
-#define NSTEPS 10000L
-
-/*
- * Matrix of Gebhart factors
- *
- * 0.065871765 0.245053213 0.281090450 0.210375872 0.0838339939 0.113774706
- * 0.183789910 0.100632183 0.291901621 0.218467251 0.0870583783 0.118150656
- * 0.187393634 0.259468108 0.121154594 0.222750923 0.0887654054 0.120467336
- * 0.180322176 0.249676859 0.286394044 0.082269756 0.0854157674 0.115921399
- * 0.167667988 0.232155676 0.266296216 0.199303457 0.0267900995 0.107786564
- * 0.170662059 0.236301313 0.271051506 0.202862448 0.0808399227 0.038282752
- */
-
-static const float vertices[] = {
- 0.f, 0.f, 0.f,
- 1.f, 0.f, 0.f,
- 0.f, 1.f, 0.f,
- 1.f, 1.f, 0.f,
- 0.f, 0.f, 1.f,
- 1.f, 0.f, 1.f,
- 0.f, 1.f, 1.f,
- 1.f, 1.f, 1.f
-};
-
-/* Front faces are CW. The normals point into the cube */
-static const unsigned indices[] = {
- 0, 2, 1, 1, 2, 3, /* Front */
- 0, 4, 2, 2, 4, 6, /* Left */
- 4, 5, 6, 6, 5, 7, /* Back */
- 3, 7, 1, 1, 7, 5, /* Right */
- 2, 6, 3, 3, 6, 7, /* Top */
- 0, 1, 4, 4, 1, 5 /* Bottom */
-};
-
-static const double emissivity[] = {
- 0.6, 0.6, /* Front */
- 0.8, 0.8, /* Left */
- 0.9, 0.9, /* Back */
- 0.7, 0.7, /* Right */
- 0.3, 0.3, /* Top */
- 0.4, 0.4 /* Bottom */
-};
-
-static const double reflectivity[] = {
- 0.4, 0.4, /* Front */
- 0.2, 0.2, /* Left */
- 0.1, 0.1, /* Back */
- 0.3, 0.3, /* Right */
- 0.7, 0.7, /* Top */
- 0.6, 0.6 /* Bottom */
-};
-
-static const double specularity[] = {
- 0.0, 0.0, /* Front */
- 0.0, 0.0, /* Left */
- 0.0, 0.0, /* Back */
- 0.0, 0.0, /* Right */
- 0.0, 0.0, /* Top */
- 0.0, 0.0, /* Bottom */
-};
-
-#if 0
-static void
-dump_mesh(void)
-{
- int i;
- FOR_EACH(i, 0, 8)
- printf("v %f %f %f\n", SPLIT3(vertices + i*3));
- FOR_EACH(i, 0, 12) {
- printf("f %d %d %d\n",
- indices[i*3 + 0] + 1,
- indices[i*3 + 1] + 1,
- indices[i*3 + 2] + 1);
- }
-}
-#endif
-
-static void
-get_ids(const unsigned itri, unsigned ids[3], void* data)
-{
- const unsigned id = itri * 3;
- (void)data;
- NCHECK(ids, NULL);
- ids[0] = indices[id + 0];
- ids[1] = indices[id + 1];
- ids[2] = indices[id + 2];
-}
-
-static void
-get_pos(const unsigned ivert, float pos[3], void* data)
-{
- const unsigned i = ivert*3;
- (void)data;
- NCHECK(pos, NULL);
- pos[0] = vertices[i + 0];
- pos[1] = vertices[i + 1];
- pos[2] = vertices[i + 2];
-}
-
-static void
-get_triangle(void* geometry, const size_t itri, float triangle[3][3])
-{
- unsigned ids[3];
- int ivert;
- (void)geometry;
-
- get_ids((unsigned)itri, ids, NULL);
- FOR_EACH(ivert, 0, 3) {
- float pos[3];
- int icoord;
- get_pos(ids[ivert], pos, NULL);
- FOR_EACH(icoord, 0, 3) {
- triangle[ivert][icoord] = pos[icoord];
- }
- }
-}
-
-static double
-get_material_property
- (void* material,
- const enum sgf_material_property prop,
- const size_t itri,
- const size_t ispec_band)
-{
- (void)material, (void)ispec_band;
- switch(prop) {
- case SGF_MATERIAL_EMISSIVITY: return emissivity[itri]; break;
- case SGF_MATERIAL_REFLECTIVITY: return reflectivity[itri]; break;
- case SGF_MATERIAL_SPECULARITY: return specularity[itri]; break;
- default: FATAL("Unreachable code\n"); break;
- }
-}
-
-int
-main(int argc, char** argv)
-{
- struct mem_allocator allocator;
- struct s3d_device* s3d;
- struct s3d_scene* scn;
- struct s3d_shape* shape;
- struct s3d_vertex_data attribs[2];
- struct smc_device* smc;
- struct sgf_context ctx = SGF_CONTEXT_NULL;
- struct sgf_accum* gfacc = NULL;
- size_t iprim;
- (void)argc, (void)argv;
-
- mem_init_proxy_allocator(&allocator, &mem_default_allocator);
-
- /*dump_mesh();
- exit(0);*/
-
- CHECK(s3d_device_create(NULL, &allocator, 1, &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, SMC_NTHREADS_DEFAULT, &smc), RES_OK);
-
- attribs[0].type = S3D_FLOAT3;
- attribs[0].usage = S3D_POSITION;
- attribs[0].get = get_pos;
- attribs[1] = S3D_VERTEX_DATA_NULL;
-
- CHECK(s3d_mesh_setup_indexed_vertices
- (shape, 12/*#tris*/, get_ids, 8/*#verts*/, attribs, NULL), RES_OK);
-
- ctx.get_triangle = get_triangle;
- ctx.geometry = NULL;
- ctx.get_material_property = get_material_property;
- ctx.material = NULL;
- ctx.primitives_count = 12;
- ctx.spectral_band_id = 0;
- ctx.scene = scn;
- ctx.logger = LOGGER_DEFAULT;
- ctx.verbose = 1;
- gfacc = sa_add(gfacc, 12*12);
-
- 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_trace(scn), RES_OK);
- FOR_EACH(iprim, 0, 12) {
- struct sgf_accum* row = gfacc + iprim*12;
- ctx.primitive_id = iprim;
- CHECK(sgf_integrate(smc, &ctx, NSTEPS, row), RES_OK);
- }
-
- CHECK(s3d_scene_end_trace(scn), RES_OK);
-
- /* Merge the radiative flux of coplanar primitives */
- FOR_EACH(iprim, 0, 6) {
- const struct sgf_accum* row_src0 = gfacc + iprim * 2 * 12/* #prims per row */;
- const struct sgf_accum* row_src1 = row_src0 + 12/* #prims per row */;
- struct sgf_accum* row_dst = gfacc + iprim * 12/* #prims per row */;
-
- int icol;
- FOR_EACH(icol, 0, 6) {
- 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;
- }
- }
-
- /* Estimated value */
- FOR_EACH(iprim, 0, 6) {
- const struct sgf_accum* row = gfacc + iprim * 12;
- int icol;
- FOR_EACH(icol, 0, 6) {
- const double E = row[icol].radiative_flux / NSTEPS;
- printf("%.6f ", E);
- }
- printf("\n");
- }
-
- printf("\n");
- /* Standard error */
- FOR_EACH(iprim, 0, 6) {
- const struct sgf_accum* row = gfacc + iprim * 12;
- int icol;
- FOR_EACH(icol, 0, 6) {
- const double V = row[icol].sqr_radiative_flux / NSTEPS
- - (row[icol].radiative_flux*row[icol].radiative_flux) / (NSTEPS*NSTEPS);
- const double SE = sqrt(V / NSTEPS);
- printf("%.6f ", SE);
- }
- printf("\n");
- }
-
- sa_release(gfacc);
-
- 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);
-
- if(MEM_ALLOCATED_SIZE(&allocator)) {
- char dump[512];
- MEM_DUMP(&allocator, dump, sizeof(dump)/sizeof(char));
- fprintf(stderr, "%s\n", dump);
- FATAL("Memory leaks\n");
- }
- mem_shutdown_proxy_allocator(&allocator);
- CHECK(mem_allocated_size(), 0);
- return 0;
-}
-
diff --git a/src/test_sgf_cube.c b/src/test_sgf_cube.c
@@ -0,0 +1,251 @@
+/* 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/stretchy_array.h>
+
+#include <star/s3d.h>
+#include <star/smc.h>
+
+#define NSTEPS 10000L
+
+/*
+ * Matrix of Gebhart factors
+ *
+ * 0.065871765 0.245053213 0.281090450 0.210375872 0.0838339939 0.113774706
+ * 0.183789910 0.100632183 0.291901621 0.218467251 0.0870583783 0.118150656
+ * 0.187393634 0.259468108 0.121154594 0.222750923 0.0887654054 0.120467336
+ * 0.180322176 0.249676859 0.286394044 0.082269756 0.0854157674 0.115921399
+ * 0.167667988 0.232155676 0.266296216 0.199303457 0.0267900995 0.107786564
+ * 0.170662059 0.236301313 0.271051506 0.202862448 0.0808399227 0.038282752
+ */
+
+static const float vertices[] = {
+ 0.f, 0.f, 0.f,
+ 1.f, 0.f, 0.f,
+ 0.f, 1.f, 0.f,
+ 1.f, 1.f, 0.f,
+ 0.f, 0.f, 1.f,
+ 1.f, 0.f, 1.f,
+ 0.f, 1.f, 1.f,
+ 1.f, 1.f, 1.f
+};
+static const size_t nvertices = sizeof(vertices) / sizeof(float[3]);
+
+/* Front faces are CW. The normals point into the cube */
+static const unsigned indices[] = {
+ 0, 2, 1, 1, 2, 3, /* Front */
+ 0, 4, 2, 2, 4, 6, /* Left */
+ 4, 5, 6, 6, 5, 7, /* Back */
+ 3, 7, 1, 1, 7, 5, /* Right */
+ 2, 6, 3, 3, 6, 7, /* Top */
+ 0, 1, 4, 4, 1, 5 /* Bottom */
+};
+static const size_t nprims = sizeof(indices) / sizeof(unsigned[3]);
+
+static const double emissivity[] = {
+ 0.6, 0.6, /* Front */
+ 0.8, 0.8, /* Left */
+ 0.9, 0.9, /* Back */
+ 0.7, 0.7, /* Right */
+ 0.3, 0.3, /* Top */
+ 0.4, 0.4 /* Bottom */
+};
+
+static const double specularity[] = {
+ 0.0, 0.0, /* Front */
+ 0.0, 0.0, /* Left */
+ 0.0, 0.0, /* Back */
+ 0.0, 0.0, /* Right */
+ 0.0, 0.0, /* Top */
+ 0.0, 0.0 /* Bottom */
+};
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct s3d_device* s3d;
+ struct s3d_scene* scn;
+ struct s3d_shape* shape;
+ 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;
+ 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_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, &smc), RES_OK);
+
+ attribs[0].type = S3D_FLOAT3;
+ attribs[0].usage = S3D_POSITION;
+ attribs[0].get = get_pos;
+ attribs[1] = S3D_VERTEX_DATA_NULL;
+
+ mesh.vertices = vertices;
+ mesh.nvertices = nvertices;
+ mesh.indices = indices;
+ mesh.ntriangles = nprims;
+
+ mtr.emissivity = emissivity;
+ mtr.specularity = specularity;
+
+ CHECK(s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, get_ids,
+ (unsigned)nvertices, attribs, &mesh), RES_OK);
+
+ ctx.get_triangle = get_triangle;
+ ctx.geometry = &mesh;
+ ctx.get_material_property = get_material_property;
+ ctx.material = &mtr;
+ ctx.primitives_count = nprims;
+ ctx.spectral_band_id = 0;
+ 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_trace(scn), 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_trace(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);
+ }
+ printf("\n");
+ }
+ printf("\n");
+#endif
+
+ /* 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;
+ size_t icol;
+
+ 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;
+
+ /* Estimated value */
+ 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 E = row[icol].radiative_flux / ((double)nsteps_adjusted);
+ printf("%.6f ", E);
+ }
+ printf("\n");
+ }
+ printf("\n");
+
+ /* Standard error */
+ FOR_EACH(iprim, 0, nprims/2) {
+ const struct sgf_accum* row = gfacc + iprim * 12;
+ 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");
+ }
+
+ sa_release(gfacc);
+
+ 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_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+ return 0;
+}
+
diff --git a/src/test_sgf_tetrahedron.c b/src/test_sgf_tetrahedron.c
@@ -0,0 +1,158 @@
+/* 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/stretchy_array.h>
+
+#include <star/s3d.h>
+#include <star/smc.h>
+
+#define NSTEPS 10000L
+
+static const float vertices[] = {
+ 0.57735026918962576450f, 0.f, 0.f,
+ -0.28867513459481288225f, 0.5f, 0.f,
+ -0.28867513459481288225f, -0.5f, 0.f,
+ 0.f, 0.f, 0.81649658092772603273f
+};
+static const size_t nvertices = sizeof(vertices) / sizeof(float[3]);
+
+/* Front faces are CW. The normals point into the cube */
+static const unsigned indices[] = {
+ 0, 1, 3,
+ 0, 2, 1,
+ 1, 2, 3,
+ 0, 3, 2
+};
+static const size_t nprims = sizeof(indices) / sizeof(unsigned[3]);
+
+static const double emissivity[] = { 0.5, 1.0, 1.0, 1.0 };
+static const double specularity[] = { 0.0, 0.0, 0.0, 0.0 };
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct s3d_device* s3d;
+ struct s3d_scene* scn;
+ struct s3d_shape* shape;
+ 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;
+ size_t iprim;
+ (void)argc, (void)argv;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(s3d_device_create(NULL, &allocator, 1, &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, SMC_NTHREADS_DEFAULT, &smc), RES_OK);
+
+ attribs[0].type = S3D_FLOAT3;
+ attribs[0].usage = S3D_POSITION;
+ attribs[0].get = get_pos;
+ attribs[1] = S3D_VERTEX_DATA_NULL;
+
+ mesh.vertices = vertices;
+ mesh.nvertices = nvertices;
+ mesh.indices = indices;
+ mesh.ntriangles = nprims;
+
+ mtr.emissivity = emissivity;
+ mtr.specularity = specularity;
+
+ CHECK(s3d_mesh_setup_indexed_vertices(shape, (unsigned)nprims, get_ids,
+ (unsigned)nvertices, attribs, &mesh), RES_OK);
+
+ ctx.get_triangle = get_triangle;
+ ctx.geometry = &mesh;
+ ctx.get_material_property = get_material_property;
+ ctx.material = &mtr;
+ ctx.primitives_count = nprims;
+ ctx.spectral_band_id = 0;
+ ctx.scene = scn;
+ ctx.logger = LOGGER_DEFAULT;
+ ctx.verbose = 1;
+
+ gfacc = sa_add(gfacc, nprims*nprims);
+
+ CHECK(s3d_scene_begin_trace(scn), 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_trace(scn), RES_OK);
+
+ /* Expected value */
+ 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);
+ }
+ printf("\n");
+ }
+ printf("\n");
+
+ /* Standard error */
+ FOR_EACH(iprim, 0, nprims) {
+ const struct sgf_accum* row = gfacc + iprim * nprims;
+ unsigned icol;
+ FOR_EACH(icol, 0, nprims) {
+ const double V = row[icol].sqr_radiative_flux / ((double)NSTEPS)
+ - (row[icol].radiative_flux*row[icol].radiative_flux)
+ / (double)(NSTEPS*NSTEPS);
+ const double SE = sqrt(V / NSTEPS);
+ printf("%.6f ", SE);
+ }
+ printf("\n");
+ }
+
+ sa_release(gfacc);
+
+ 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_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
@@ -0,0 +1,138 @@
+/* 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 TEST_SGF_UTILS_H
+#define TEST_SGF_UTILS_H
+
+#include <rsys/mem_allocator.h>
+#include <stdio.h>
+
+struct triangle_mesh {
+ const float* vertices;
+ size_t nvertices;
+ const unsigned* indices;
+ size_t ntriangles;
+};
+
+struct material {
+ const double* emissivity;
+ const double* specularity;
+};
+
+static void
+get_ids(const unsigned itri, unsigned ids[3], void* data)
+{
+ const struct triangle_mesh* mesh = data;
+ const unsigned id = itri * 3;
+ NCHECK(mesh, NULL);
+ NCHECK(ids, NULL);
+ CHECK(itri < mesh->ntriangles, 1);
+ ids[0] = mesh->indices[id + 0];
+ ids[1] = mesh->indices[id + 1];
+ ids[2] = mesh->indices[id + 2];
+}
+
+static void
+get_pos(const unsigned ivert, float pos[3], void* data)
+{
+ const struct triangle_mesh* mesh = data;
+ const unsigned i = ivert*3;
+ NCHECK(mesh, NULL);
+ NCHECK(pos, NULL);
+ CHECK(ivert < mesh->nvertices, 1);
+ pos[0] = mesh->vertices[i + 0];
+ pos[1] = mesh->vertices[i + 1];
+ pos[2] = mesh->vertices[i + 2];
+}
+
+static void
+get_triangle(void* geometry, const size_t itri, float triangle[3][3])
+{
+ struct triangle_mesh* mesh = geometry;
+ unsigned ids[3];
+ int ivert;
+ NCHECK(mesh, NULL);
+ NCHECK(triangle, NULL);
+ CHECK(itri < mesh->ntriangles, 1);
+ get_ids((unsigned)itri, ids, mesh);
+ FOR_EACH(ivert, 0, 3) {
+ float pos[3];
+ int icoord;
+ get_pos(ids[ivert], pos, mesh);
+ FOR_EACH(icoord, 0, 3) {
+ triangle[ivert][icoord] = pos[icoord];
+ }
+ }
+}
+
+static double
+get_material_property
+ (void* material,
+ const enum sgf_material_property prop,
+ const size_t itri,
+ const size_t ispec_band)
+{
+ struct material* mtr = material;
+ NCHECK(material, NULL);
+ (void)ispec_band;
+ switch(prop) {
+ case SGF_MATERIAL_EMISSIVITY: return mtr->emissivity[itri]; break;
+ case SGF_MATERIAL_REFLECTIVITY: return 1.0 - mtr->emissivity[itri]; break;
+ case SGF_MATERIAL_SPECULARITY: return mtr->specularity[itri]; break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+}
+
+static INLINE void
+dump_triangle_mesh(struct triangle_mesh* mesh)
+{
+ unsigned i;
+ NCHECK(mesh, NULL);
+ FOR_EACH(i, 0, mesh->nvertices)
+ printf("v %f %f %f\n", SPLIT3(mesh->vertices + i*3));
+ FOR_EACH(i, 0, mesh->ntriangles) {
+ printf("f %d %d %d\n",
+ mesh->indices[i*3 + 0] + 1,
+ mesh->indices[i*3 + 1] + 1,
+ mesh->indices[i*3 + 2] + 1);
+ }
+}
+
+static void
+check_memory_allocator(struct mem_allocator* allocator)
+{
+ if(MEM_ALLOCATED_SIZE(allocator)) {
+ char dump[512];
+ MEM_DUMP(allocator, dump, sizeof(dump)/sizeof(char));
+ fprintf(stderr, "%s\n", dump);
+ FATAL("Memory leaks\n");
+ }
+}
+
+#endif /* TEST_SGF_UTILS_H */
+
