star-gf

commit 260323f4d715bf7e811ff94dd38d2bbd7e1601f9
parent f1763a37e41c06568935650bb4083ca47a89f5ef
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 18 Oct 2016 11:13:52 +0200

Fix the accumulation of the MC weights

Diffstat:
M
src/sgf_estimator.c
 | 
6
+++---
M
src/sgf_realisation.h
 | 
66
+++++++++++++++++++++++++++++++++++-------------------------------

2 files changed, 38 insertions(+), 34 deletions(-)
diff --git a/src/sgf_estimator.c b/src/sgf_estimator.c
@@ -150,7 +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 htable_bounce path;
   struct sgf_estimator* estimator = NULL;
   size_t istep;
   size_t iband;
@@ -163,7 +163,7 @@ sgf_integrate
   gebhart_radiative_path_T gebhart_radiative_path;
 
   if(!dev) return RES_BAD_ARG;
-  darray_bounce_init(dev->allocator, &path);
+  htable_bounce_init(dev->allocator, &path);
 
   if(!steps_count || !rng || !desc || !check_scene_desc(desc) || !out_estimator) {
     res = RES_BAD_ARG;
@@ -295,7 +295,7 @@ sgf_integrate
 
 exit:
   if(out_estimator) *out_estimator = estimator;
-  darray_bounce_release(&path);
+  htable_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,9 +23,9 @@
 #ifndef SGF_REALISATION_H
 #define SGF_REALISATION_H
 
-#include <rsys/dynamic_array.h>
 #include <rsys/float2.h>
 #include <rsys/float3.h>
+#include <rsys/hash_table.h>
 
 #include <star/s2d.h>
 #include <star/s3d.h>
@@ -41,16 +41,11 @@ struct accum {
   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>
+/* Declare the htable_bounce data structure */
+#define HTABLE_NAME bounce
+#define HTABLE_KEY unsigned /* Primitive ID */
+#define HTABLE_DATA double /* Weight */
+#include <rsys/hash_table.h>
 
 /* Type of the realisation function. Return RES_BAD_OP on numerical issue. i.e.
  * the radiative path is skipped */
@@ -61,7 +56,7 @@ typedef res_T
    struct accum* accum_infinity,
    struct accum* accum_medium,
    struct ssp_rng* rng,
-   struct darray_bounce* path, /* Store the path bounces */
+   struct htable_bounce* path, /* Store the path */
    const double absorption_coef, /* In m^-1  */
    const size_t ispectral_band,
    const float epsilon,
@@ -219,7 +214,7 @@ GEBHART_RADIATIVE_PATH
    struct accum* accum_infinity,
    struct accum* accum_medium,
    struct ssp_rng* rng,
-   struct darray_bounce* path,
+   struct htable_bounce* path,
    const double absorption_coef, /* m^-1 */
    const size_t ispectral_band,
    const float epsilon,
@@ -229,8 +224,8 @@ GEBHART_RADIATIVE_PATH
   sXd_scene_T* scene;
   sXd_hit_T hit;
   sXd_primitive_T prim_from, prim;
+  struct htable_bounce_iterator it, end;
   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;
@@ -251,7 +246,7 @@ GEBHART_RADIATIVE_PATH
   ASSERT(absorption_coef >= 0 || accum_infinity);
 
   scene = sgf_scene_desc_get_sXd_scene(desc);
-  darray_bounce_clear(path);
+  htable_bounce_clear(path);
 
   /* Discard faces with no emissivity */
   emissivity = desc->get_material_property(desc->material,
@@ -334,25 +329,30 @@ GEBHART_RADIATIVE_PATH
       SGF_MATERIAL_REFLECTIVITY, prim.scene_prim_id, ispectral_band);
 
     if(transmissivity > trans_min) {
-      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;
+      const double weight = transmissivity * emissivity;
+      double* pweight = htable_bounce_find(path, &prim.scene_prim_id);
+      if(pweight) {
+        *pweight += weight;
+      } else {
+        res = htable_bounce_set(path, &prim.scene_prim_id, &weight);
+        if(res != RES_OK) return res;
+      }
       transmissivity = transmissivity * (1.0 - emissivity);
 #ifndef NDEBUG
-      radiative_flux += bounce.weight;
+      radiative_flux += weight;
 #endif
     } else {
       /* Russian roulette */
       if(ssp_rng_canonical(rng) < emissivity) {
-        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;
+        double* pweight = htable_bounce_find(path, &prim.scene_prim_id);
+        if(pweight) {
+          *pweight += transmissivity;
+        } else {
+          res = htable_bounce_set(path, &prim.scene_prim_id, &transmissivity);
+          if(res != RES_OK) return res;
+        }
 #ifndef NDEBUG
-        radiative_flux += bounce.weight;
+        radiative_flux += transmissivity;
 #endif
         break;
       }
@@ -374,10 +374,14 @@ GEBHART_RADIATIVE_PATH
     }
   }
 
-  /* 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);
+  /* Flush MC weights */
+  htable_bounce_begin(path, &it);
+  htable_bounce_end(path, &end);
+  while(!htable_bounce_iterator_eq(&it, &end)) {
+    const size_t iprim = *htable_bounce_iterator_key_get(&it);
+    const double weight = *htable_bounce_iterator_data_get(&it);
+    accum_weight(accums + iprim, weight);
+    htable_bounce_iterator_next(&it);
   }
   if(medium_radiative_flux) {
     accum_weight(accum_medium, medium_radiative_flux);