commit c30ea2d0f67698513c8213dec6f33335e1785d40
parent c22998c6d97f1fc27a9976224b20346056ed52a0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 24 Jul 2015 17:25:55 +0200
Change the API of the sgf_context
Replace the spectral band identifier by the number of spectral bands.
The sgf_integrate use this data to compute the Gebhart Factor for *all*
spectral bands. Up to now, only one spectral band was computed per
sgf_integrate call. With this update, the radiative flux of each
spectral band of a primitive are computed and stored sequentially.
Diffstat:
4 files changed, 64 insertions(+), 34 deletions(-)
diff --git a/src/sgf.c b/src/sgf.c
@@ -55,9 +55,11 @@ static res_T
gebhart_radiative_path
(struct sgf_accum* result,
struct ssp_rng* rng,
- struct sgf_context* ctx)
+ struct sgf_context* ctx,
+ const size_t ispectral_band)
{
struct s3d_hit hit;
+ struct sgf_accum* gebfac; /* Gebhart factor */
double transmissivity = 1.0;
double proba_reflec_spec;
double emissivity;
@@ -65,18 +67,17 @@ gebhart_radiative_path
double specularity;
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[4];
float range[2];
- float dist;
+ ASSERT(ispectral_band < ctx->spectral_bands_count);
/* Discard faces with no emissivity */
emissivity = ctx->get_material_property(ctx->material,
- SGF_MATERIAL_EMISSIVITY, ctx->primitive_id, ctx->spectral_band_id);
+ SGF_MATERIAL_EMISSIVITY, ctx->primitive_id, ispectral_band);
if(emissivity <= 0.0)
return RES_OK;
@@ -122,31 +123,38 @@ gebhart_radiative_path
f3_add(pos, vec0, pos);
iprim = hit.prim.prim_id;
ASSERT(iprim < ctx->primitives_count);
-#ifndef NDEBUG
+
+ /* FIXME Why ? */
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));
+ ssp_ran_triangle_uniform(rng, triangle[0], triangle[1], triangle[2], pos);
+#ifndef NDEBUG
+ {
+ /*Ax + By + Cz + D*/
+ const float dist = f3_dot(pos, normal) - f3_dot(normal, triangle[0]);
+ ASSERT(eq_eps(dist, 0, 1.e-6f));
+ }
#endif
/* Fetch material property */
emissivity = ctx->get_material_property
- (ctx->material, SGF_MATERIAL_EMISSIVITY, iprim, ctx->spectral_band_id);
+ (ctx->material, SGF_MATERIAL_EMISSIVITY, iprim, ispectral_band);
specularity = ctx->get_material_property
- (ctx->material, SGF_MATERIAL_SPECULARITY, iprim, ctx->spectral_band_id);
+ (ctx->material, SGF_MATERIAL_SPECULARITY, iprim, ispectral_band);
reflectivity = ctx->get_material_property
- (ctx->material, SGF_MATERIAL_REFLECTIVITY, iprim, ctx->spectral_band_id);
+ (ctx->material, SGF_MATERIAL_REFLECTIVITY, iprim, ispectral_band);
+ gebfac = result + iprim * ctx->spectral_bands_count;
if(transmissivity > trans_min) {
const double weight = transmissivity * emissivity;
- result[iprim].radiative_flux += weight;
- result[iprim].sqr_radiative_flux += weight * weight;
+ 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;
- result[iprim].radiative_flux += weight;
- result[iprim].sqr_radiative_flux += weight * weight;
+ gebfac[ispectral_band].radiative_flux += weight;
+ gebfac[ispectral_band].sqr_radiative_flux += weight * weight;
break;
}
}
@@ -166,15 +174,15 @@ gebhart_radiative_path
f3_add(dir, dir, vec0);
f3_normalize(dir, dir);
}
- ++nbounce;
ctx->get_triangle(ctx->geometry, iprim, triangle);
}
-/* ASSERT(nbounce <= 1);*/
#if 1 && !defined(NDEBUG)
{ /* Ensure the energy conservation property */
double sum_radiative_flux = 0.0;
+ gebfac = result;
FOR_EACH(iprim, 0, ctx->primitives_count) {
- sum_radiative_flux += result[iprim].radiative_flux;
+ sum_radiative_flux += gebfac[ispectral_band].radiative_flux;
+ gebfac += ctx->spectral_bands_count;
}
ASSERT(eq_eps(sum_radiative_flux, 1.0, 1.e-6));
}
@@ -189,22 +197,29 @@ gebhart_factor_integrand
void* context)
{
struct sgf_context* ctx = context;
+ struct sgf_accum* buf;
size_t nfailures = 0;
+ size_t iband = 0;
res_T res = RES_OK;
- do {
- accum_buffer_clear(result);
- res = gebhart_radiative_path(darray_gfacc_data_get(result), rng, ctx);
- } while(res == RES_BAD_OP && (++nfailures < MAX_FAILURES));
+ ASSERT(result && rng && context);
+
+ buf = darray_gfacc_data_get(result);
+ FOR_EACH(iband, 0, ctx->spectral_bands_count) {
+ do {
+ accum_buffer_clear(result);
+ res = gebhart_radiative_path(buf, rng, ctx, iband);
+ } while(res == RES_BAD_OP && (++nfailures < MAX_FAILURES));
- if(ctx->verbose) {
- if(nfailures >= MAX_FAILURES) {
- logger_print(ctx->logger, LOG_ERROR,
+ if(ctx->verbose) {
+ if(nfailures >= MAX_FAILURES) {
+ logger_print(ctx->logger, LOG_ERROR,
"Too many radiative path failures. The geometry might be wrongly setuped. The\n"
"geometry should be closed and 2D-manifold. In addition its normals must point\n"
"*into* the volume with respect to a Clock Wise vertex ordering\n");
- } else if(res != RES_OK) {
- logger_print(ctx->logger, LOG_ERROR,
- "Critical error in the computation of the gebhart factor\n");
+ } else if(res != RES_OK) {
+ logger_print(ctx->logger, LOG_ERROR,
+ "Critical error in the computation of the Gebhart factor\n");
+ }
}
}
}
diff --git a/src/sgf.h b/src/sgf.h
@@ -69,7 +69,7 @@ struct sgf_context {
size_t primitives_count; /* Total number of primitives */
size_t primitive_id; /* Triangle id on which the integration is performed */
- size_t spectral_band_id; /* Identifier of the spectral band */
+ 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 */
@@ -88,10 +88,17 @@ struct sgf_accum {
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_integrate
(struct smc_device* smc,
- struct sgf_context* contex,
+ struct sgf_context* context,
const size_t steps_count,
struct sgf_accum* accum_buffer);
diff --git a/src/test_sgf_cube.c b/src/test_sgf_cube.c
@@ -135,7 +135,7 @@ main(int argc, char** argv)
ctx.get_material_property = get_material_property;
ctx.material = &mtr;
ctx.primitives_count = nprims;
- ctx.spectral_band_id = 0;
+ ctx.spectral_bands_count = 1;
ctx.scene = scn;
ctx.logger = LOGGER_DEFAULT;
ctx.verbose = 1;
@@ -214,17 +214,21 @@ main(int argc, char** argv)
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 * 12;
+ 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)
diff --git a/src/test_sgf_tetrahedron.c b/src/test_sgf_tetrahedron.c
@@ -35,7 +35,7 @@
#include <star/s3d.h>
#include <star/smc.h>
-#define NSTEPS 10000L
+#define NSTEPS 100000L
static const float vertices[] = {
0.57735026918962576450f, 0.f, 0.f,
@@ -79,7 +79,7 @@ main(int argc, char** argv)
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);
+ CHECK(smc_device_create(NULL, NULL, 1, &smc), RES_OK);
attribs[0].type = S3D_FLOAT3;
attribs[0].usage = S3D_POSITION;
@@ -102,7 +102,7 @@ main(int argc, char** argv)
ctx.get_material_property = get_material_property;
ctx.material = &mtr;
ctx.primitives_count = nprims;
- ctx.spectral_band_id = 0;
+ ctx.spectral_bands_count = 1;
ctx.scene = scn;
ctx.logger = LOGGER_DEFAULT;
ctx.verbose = 1;
@@ -121,11 +121,15 @@ main(int argc, char** argv)
FOR_EACH(iprim, 0, nprims) {
const struct sgf_accum* row = gfacc + iprim * nprims;
unsigned icol;
+ double sum = 0.0;
FOR_EACH(icol, 0, nprims) {
const double E = row[icol].radiative_flux / NSTEPS;
+ sum += E;
printf("%.6f ", E);
}
printf("\n");
+ /* Ensure the energy conservation property */
+ CHECK(eq_eps(sum, 1.0, 1.e-6), 1);
}
printf("\n");
