commit a9633e482b551d456a5dfa959991d93b02edffb6
parent ba0cef06873f1ffa058b68f3ab9828e5389897df
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 17 Oct 2016 16:34:53 +0200
Fix a bug in the accumulation of the MC weight
The expected value was correctly estimated but the variance was totally
wrong.
Diffstat:
2 files changed, 56 insertions(+), 19 deletions(-)
diff --git a/src/sgf_estimator.c b/src/sgf_estimator.c
@@ -65,6 +65,7 @@ check_scene_desc(struct sgf_scene_desc* desc)
{
ASSERT(desc);
if(!desc->get_material_property) return 0;
+ if(!desc->spectral_bands_count) return 0;
switch(desc->dimensionality) {
case SGF_2D: return desc->geometry.scn2d != NULL;
case SGF_3D: return desc->geometry.scn3d != NULL;
@@ -149,6 +150,7 @@ sgf_integrate
{
#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 darray_bounce path;
struct sgf_estimator* estimator = NULL;
size_t istep;
size_t iband;
@@ -160,8 +162,10 @@ sgf_integrate
res_T res = RES_OK;
gebhart_radiative_path_T gebhart_radiative_path;
- if(!dev || !steps_count || !rng || !desc || !check_scene_desc(desc)
- || !out_estimator) {
+ if(!dev) return RES_BAD_ARG;
+ darray_bounce_init(dev->allocator, &path);
+
+ if(!steps_count || !rng || !desc || !check_scene_desc(desc) || !out_estimator) {
res = RES_BAD_ARG;
goto error;
}
@@ -281,7 +285,7 @@ sgf_integrate
FOR_EACH(istep, 0, steps_count) {
res = gebhart_radiative_path(dev, accums, accum_infinity, accum_medium,
- rng, ka, iband, epsilon, iprim, desc);
+ rng, &path, ka, iband, epsilon, iprim, desc);
if(res != RES_OK) goto error;
}
}
@@ -291,6 +295,7 @@ sgf_integrate
exit:
if(out_estimator) *out_estimator = estimator;
+ darray_bounce_release(&path);
return res;
error:
if(estimator) SGF(estimator_ref_put(estimator));
diff --git a/src/sgf_realisation.h b/src/sgf_realisation.h
@@ -23,21 +23,35 @@
#ifndef SGF_REALISATION_H
#define SGF_REALISATION_H
+#include <rsys/dynamic_array.h>
+#include <rsys/float2.h>
+#include <rsys/float3.h>
+
#include <star/s2d.h>
#include <star/s3d.h>
#include <star/ssp.h>
-#include <rsys/float2.h>
-#include <rsys/float3.h>
/* How many self intersections are tested on the same ray before an error
* occurs */
#define NSELF_HITS_MAX 10
+/* Monte Carlo estimator */
struct accum {
double radiative_flux;
double sqr_radiative_flux;
};
+/* Define a radiative path bounce */
+struct bounce {
+ size_t iprim; /* Current primitive */
+ double weight; /* Current MC weight */
+};
+
+/* Declare the darray_bounce data structure */
+#define DARRAY_NAME bounce
+#define DARRAY_DATA struct bounce
+#include <rsys/dynamic_array.h>
+
/* Type of the realisation function. Return RES_BAD_OP on numerical issue. i.e.
* the radiative path is skipped */
typedef res_T
@@ -47,7 +61,8 @@ typedef res_T
struct accum* accum_infinity,
struct accum* accum_medium,
struct ssp_rng* rng,
- const double absorption_coef,
+ struct darray_bounce* path, /* Store the path bounces */
+ const double absorption_coef, /* In m^-1 */
const size_t ispectral_band,
const float epsilon,
const size_t primitive_id,
@@ -204,6 +219,7 @@ GEBHART_RADIATIVE_PATH
struct accum* accum_infinity,
struct accum* accum_medium,
struct ssp_rng* rng,
+ struct darray_bounce* path,
const double absorption_coef, /* m^-1 */
const size_t ispectral_band,
const float epsilon,
@@ -214,25 +230,28 @@ GEBHART_RADIATIVE_PATH
sXd_hit_T hit;
sXd_primitive_T prim_from, prim;
size_t nself_hits = 0;
+ size_t i;
const double trans_min = 1.e-8; /* Minimum transmissivity threshold */
double proba_reflec_spec;
double transmissivity, emissivity, reflectivity, specularity;
double medium_transmissivity;
#ifndef NDEBUG
double radiative_flux = 0.0;
+#endif
double infinite_radiative_flux = 0.0;
double medium_radiative_flux = 0.0;
-#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 */
+ res_T res = RES_OK;
ASSERT(accums && epsilon > 0.f && rng && desc);
ASSERT(absorption_coef < 0 || accum_medium);
ASSERT(absorption_coef >= 0 || accum_infinity);
scene = sgf_scene_desc_get_sXd_scene(desc);
+ darray_bounce_clear(path);
/* Discard faces with no emissivity */
emissivity = desc->get_material_property(desc->material,
@@ -291,19 +310,13 @@ GEBHART_RADIATIVE_PATH
medium_transmissivity = exp(-absorption_coef*hit.distance);
weight = transmissivity * (1-medium_transmissivity);
- accum_weight(accum_medium, weight);
transmissivity *= medium_transmissivity;
-#ifndef NDEBUG
medium_radiative_flux += weight;
-#endif
}
}
if(SXD_HIT_NONE(&hit)) { /* The ray is outside the volume */
- accum_weight(accum_infinity, transmissivity);
-#ifndef NDEBUG
infinite_radiative_flux = transmissivity;
-#endif
break;
}
@@ -321,19 +334,25 @@ GEBHART_RADIATIVE_PATH
SGF_MATERIAL_REFLECTIVITY, prim.scene_prim_id, ispectral_band);
if(transmissivity > trans_min) {
- const double weight = transmissivity * emissivity;
- accum_weight(accums + prim.scene_prim_id, weight);
+ struct bounce bounce;
+ bounce.weight = transmissivity * emissivity;
+ bounce.iprim = prim.scene_prim_id;
+ res = darray_bounce_push_back(path, &bounce);
+ if(res != RES_OK) return res;
transmissivity = transmissivity * (1.0 - emissivity);
#ifndef NDEBUG
- radiative_flux += weight;
+ radiative_flux += bounce.weight;
#endif
} else {
/* Russian roulette */
if(ssp_rng_canonical(rng) < emissivity) {
- const double weight = transmissivity;
- accum_weight(accums + prim.scene_prim_id, weight);
+ struct bounce bounce;
+ bounce.weight = transmissivity;
+ bounce.iprim = prim.scene_prim_id;
+ res = darray_bounce_push_back(path, &bounce);
+ if(res != RES_OK) return res;
#ifndef NDEBUG
- radiative_flux += weight;
+ radiative_flux += bounce.weight;
#endif
break;
}
@@ -354,6 +373,19 @@ GEBHART_RADIATIVE_PATH
f3_normalize(dir, dir);
}
}
+
+ /* Store radiative path weights */
+ FOR_EACH(i, 0, darray_bounce_size_get(path)) {
+ const struct bounce* bounce = darray_bounce_cdata_get(path) + i;
+ accum_weight(accums + bounce->iprim, bounce->weight);
+ }
+ if(medium_radiative_flux) {
+ accum_weight(accum_medium, medium_radiative_flux);
+ }
+ if(infinite_radiative_flux) {
+ accum_weight(accum_infinity, infinite_radiative_flux);
+ }
+
#if !defined(NDEBUG)
/* Check the energy conservation property */
ASSERT(eq_eps
