star-schiff

commit 5943fc9b15460ffbcad3ad876a721f553faa62fb
parent 191dee04476221ae1f25a4878895af71b5437202
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 19 Oct 2015 15:00:42 +0200

Make the Schiff estimator more robust

Handle errors in sampled geometries, i.e. geometries that do not define
a closed volume.

Diffstat:
M
src/sschiff_estimator.c
 | 
31
+++++++++++++++++++++++--------

1 file changed, 23 insertions(+), 8 deletions(-)
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -44,6 +44,8 @@
 
 #include <math.h>
 
+#define MAX_FAILURES 10
+
 /* Type of computed Monte Carlo weights */
 enum weight_type {
   WEIGHT_EXTINCTION_CROSS_SECTION,
@@ -196,6 +198,7 @@ compute_monte_carlo_weight
   struct s3d_hit hit;
   double sample[2];
   const float range[2] = { 0, FLT_MAX };
+  size_t nfailures = 0;
   float axis_x[3], axis_y[3], axis_z[3];
   float plane_size[2];
   float ray_org[3];
@@ -236,7 +239,10 @@ compute_monte_carlo_weight
       double path_length;
 
       res = compute_path_length(dev, scn, &hit, ray_org, axis_z, &path_length);
-      if(res != RES_OK) continue;
+      if(res != RES_OK) {
+        ++nfailures;
+        continue;
+      }
 
       n_e = props->medium_refractive_index;
       n_r = props->relative_real_refractive_index;
@@ -257,9 +263,13 @@ compute_monte_carlo_weight
 
       weights[WEIGHT_AVERAGE_PROJECTED_AREA] = rcp_pdf;
     }
-  } while(res != RES_OK); /* TODO handle infinite loop */
-
-  return RES_OK;
+  } while(res != RES_OK && nfailures < MAX_FAILURES);
+  if(res != RES_OK) {
+    log_error(dev,
+"Too many failures in computing the radiative path length. The sampled geometry\n"
+"might not defined a closed volume.\n");
+  }
+  return res;
 }
 
 /* Cast rays into the RT volume */
@@ -353,6 +363,7 @@ radiative_properties
   float centroid[3];
   size_t idir;
   int i;
+  res_T res = RES_OK;
   ASSERT(dev && rng && ndirs && props && accums);
   (void)istep;
 
@@ -403,11 +414,11 @@ radiative_properties
       f3_max(upper, upper, pt[i]);
     }
 
-    compute_monte_carlo_weight(dev, scene, rng, props, wavelength, basis,
+    res = compute_monte_carlo_weight(dev, scene, rng, props, wavelength, basis,
       centroid, lower, upper, weights);
+    if(res != RES_OK) goto error;
 
-    FOR_EACH(i, 0, WEIGHTS_COUNT)
-      weights_dirs[i] += weights[i];
+    FOR_EACH(i, 0, WEIGHTS_COUNT) weights_dirs[i] += weights[i];
 #if 0
     /* Compute an image of the shape transformed in the footprint space */
     {
@@ -423,7 +434,11 @@ radiative_properties
     accums[i].sqr_weight += weights_dirs[i] * weights_dirs[i];
     accums[i].naccums += 1;
   }
-  return RES_OK;
+exit:
+  return res;
+error:
+  memset(accums, 0, sizeof(struct accum)*WEIGHTS_COUNT);
+  goto exit;
 }
 
 static void