star-sf

commit 0a569165b52292bd52dac89680da988386940ada
parent 5d6ed65f64f2cd87d92abb8953dc307606b3e349
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Fri, 25 Aug 2017 16:26:24 +0200

Fix thin dielectric non-conservation bug.

Diffstat:
M
src/ssf_optics.h
 | 
2
+-
M
src/ssf_thin_specular_dielectric.c
 | 
12
++++++++----
M
src/test_ssf_thin_specular_dielectric.c
 | 
36
++++++++++++++++++++++++++++++++++--

3 files changed, 43 insertions(+), 7 deletions(-)
diff --git a/src/ssf_optics.h b/src/ssf_optics.h
@@ -53,7 +53,7 @@ refract(double res[3], const double V[3], const double N[3], const double eta)
   cos_theta_i = d3_dot(V, N);
   sin2_theta_i = MMAX(0, 1.0 - cos_theta_i*cos_theta_i);
   sin2_theta_t = eta * eta * sin2_theta_i;
-  if(sin2_theta_t >= 1) return NULL; /* Total internal reflection */
+  if(sin2_theta_t >= 1) return NULL; /* Total reflection */
   cos_theta_t = sqrt(1 - sin2_theta_t);
 
   d3_muld(tmp0, V, eta);
diff --git a/src/ssf_thin_specular_dielectric.c b/src/ssf_thin_specular_dielectric.c
@@ -73,7 +73,7 @@ thin_specular_dielectric_sample
   (void)rng;
 
   eta = bsdf->eta_i / bsdf->eta_t;
-  if(!refract(wt, wo, N, eta)) { /* Total internal reflection */
+  if(!refract(wt, wo, N, eta)) { /* Total reflection */
     reflect(wi, wo, N);
     *pdf = INF;
     *type = SSF_SPECULAR | SSF_REFLECTION;
@@ -107,14 +107,18 @@ thin_specular_dielectric_sample
   *pdf = INF;
 
   /* Importance sample the BTDF wrt R */
-  if(ssp_rng_canonical(rng) < R) { /* Sample the reflective part */
+  if(ssp_rng_uniform_double(rng, 0, R + T) < R) {
     reflect(wi, wo, N);
     *type = SSF_SPECULAR | SSF_REFLECTION;
-    return 1;
   } else { /* Sample the transmissive part */
     d3_minus(wi, wo);
     *type = SSF_SPECULAR | SSF_TRANSMISSION;
-    return T;
+  }
+  if(bsdf->absorption == 0) {
+    /* avoid numerical loss of energy if no absorption */
+    return 1;
+  } else {
+    return R + T;
   }
 }
 
diff --git a/src/test_ssf_thin_specular_dielectric.c b/src/test_ssf_thin_specular_dielectric.c
@@ -21,7 +21,7 @@
 int
 main(int argc, char** argv)
 {
-  const size_t NSTEPS = 100000;
+  const size_t NSTEPS = 1000000;
   struct mem_allocator allocator;
   struct ssp_rng* rng;
   struct ssf_bxdf* bsdf;
@@ -33,6 +33,8 @@ main(int argc, char** argv)
   double refract[3];
   double pdf;
   double R; /* Directional reflectance */
+  double SR, ST;
+  double SR2 = 0, ST2 = 0;
   size_t i;
   int type;
   (void)argc, (void)argv;
@@ -109,10 +111,40 @@ main(int argc, char** argv)
   d3_normalize(wo, d3(wo, 1, 0.0000001, 0));
   CHECK(ssf_thin_specular_dielectric_setup(bsdf, 0, 1.4, 1.0, 1), RES_OK);
   R = ssf_bxdf_sample(bsdf, rng, wo, N, wi, &type, &pdf);
-  NCHECK(R, 0);
+  CHECK(R, 1);
   CHECK(IS_INF(pdf), 1);
   CHECK(d3_eq_eps(wi, d3(tmp, -wo[0], wo[1], 0), 1.e-6), 1);
   CHECK(type, SSF_SPECULAR | SSF_REFLECTION);
+  FOR_EACH(i, 0, NSTEPS) {
+    R = ssf_bxdf_sample(bsdf, rng, wo, N, wi, &type, &pdf);
+    CHECK(type & SSF_REFLECTION, 1);
+    CHECK(R, 1);
+  }
+
+  /* Check T VS R proportion and E conservation */
+  SR = ST = 0;
+  SR2 = ST2 = 0;
+  d3(wo, 0.0, 1.0, 0.0);
+  FOR_EACH(i, 0, NSTEPS) {
+    R = ssf_bxdf_sample(bsdf, rng, wo, N, wi, &type, &pdf);
+    if (type & SSF_TRANSMISSION) {
+      ST += R; ST2 += R*R;
+    }
+    else {
+      SR += R; SR2 += R*R;
+    }
+  }
+
+#define MEAN(x, n) ((x) / (n))
+#define VAR(x, x2, n) (MEAN((x2), (n)) - MEAN((x), (n))*MEAN((x), (n)))
+#define STD(x, x2, n) (VAR((x), (x2), (n)) > 0 ? sqrt(VAR((x), (x2), (n)) / (n)) : 0)
+  /* Check conservation of energy */
+  CHECK(MEAN(SR+ST, NSTEPS), 1);
+  /* Check T VS R proportion */
+  CHECK(eq_eps(MEAN(SR, NSTEPS), 0.0540540540, 3 * STD(SR,SR2,NSTEPS)), 1);
+#undef MEAN
+#undef VAR
+#undef STD
 
   wo[0] = ssp_rng_uniform_double(rng, -1, 1);
   wo[1] = ssp_rng_uniform_double(rng, -1, 1);