star-sf

commit 9a86fd81ba0a61fbb65e0bf0e757913fc7bc1b85
parent 5b1c996b2ff3ed94cf48041b0f6f34d7020b3057
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 14 Sep 2016 16:15:59 +0200

Rewrite the ssf_bsdf_sample function

Diffstat:
M
src/ssf_bsdf.c
 | 
101
+++++++++++++++++++++++++++++++++++++++++--------------------------------------

1 file changed, 52 insertions(+), 49 deletions(-)
diff --git a/src/ssf_bsdf.c b/src/ssf_bsdf.c
@@ -21,12 +21,12 @@
 #include <rsys/mem_allocator.h>
 #include <rsys/ref_count.h>
 
-#define MAX_BxDFs 8
+#define MAX_COMPONENTS 8
 
 struct ssf_bsdf {
-  struct ssf_bxdf* bxdfs[MAX_BxDFs];
-  double weights[MAX_BxDFs];
-  size_t nbxdfs;
+  struct ssf_bxdf* components[MAX_COMPONENTS];
+  double weights[MAX_COMPONENTS];
+  size_t ncomponents;
 
   ref_T ref;
   struct mem_allocator* allocator;
@@ -99,11 +99,12 @@ res_T
 ssf_bsdf_add(struct ssf_bsdf* bsdf, struct ssf_bxdf* bxdf, const double weight)
 {
   if(!bsdf || !bxdf || weight < 0) return RES_BAD_ARG;
-  if(bsdf->nbxdfs >= MAX_BxDFs) return RES_MEM_ERR;
+  if(bsdf->ncomponents >= MAX_COMPONENTS) return RES_MEM_ERR;
+  if(weight == 0) return RES_OK;
   SSF(bxdf_ref_get(bxdf));
-  bsdf->bxdfs[bsdf->nbxdfs] = bxdf;
-  bsdf->weights[bsdf->nbxdfs] = weight;
-  ++bsdf->nbxdfs;
+  bsdf->components[bsdf->ncomponents] = bxdf;
+  bsdf->weights[bsdf->ncomponents] = weight;
+  ++bsdf->ncomponents;
   return RES_OK;
 }
 
@@ -113,16 +114,13 @@ ssf_bsdf_clear(struct ssf_bsdf* bsdf)
 {
   size_t i;
   if(!bsdf) return RES_BAD_ARG;
-  FOR_EACH(i, 0, bsdf->nbxdfs) {
-    SSF(bxdf_ref_put(bsdf->bxdfs[i]));
+  FOR_EACH(i, 0, bsdf->ncomponents) {
+    SSF(bxdf_ref_put(bsdf->components[i]));
   }
-  bsdf->nbxdfs = 0;
+  bsdf->ncomponents = 0;
   return RES_OK;
 }
 
-/* TODO fix the use of the u & v canonical random number. Currently u is used 2
- * times
- * TODO fix the sampling that is actually totally wrong */
 double
 ssf_bsdf_sample
   (struct ssf_bsdf* bsdf,
@@ -132,48 +130,53 @@ ssf_bsdf_sample
    const double N[3],
    double wo[4])
 {
-  const size_t PDF = 3;
-  double reflectivities[MAX_BxDFs];
-  double dirs[MAX_BxDFs][4];
-  double probas[MAX_BxDFs];
-  double cumul[MAX_BxDFs];
+  double probas[MAX_COMPONENTS];
+  double cumul[MAX_COMPONENTS];
   double probas_sum;
-  size_t i, n;
-
-  ASSERT(bsdf && u>=0 && u<1 && v>=0 && v<1 && wi && N && wo);
-  ASSERT(d3_is_normalized(wi) && d3_is_normalized(N));
-
-  /* Build the probability distribution by sampling each BRDF */
-  n = 0;
-  probas_sum = 0.0f;
-  FOR_EACH(i, 0, bsdf->nbxdfs) {
-    struct ssf_bxdf* bxdf = bsdf->bxdfs[i];
-
-    reflectivities[n] = ssf_bxdf_sample(bxdf, u, v, wi, N, dirs[n]);
-    if(reflectivities[n] <= 0 || dirs[n][PDF] <= 0)
-      continue; /* Discard BxDF */
+  double reflectivity;
+  double w;
+  double pdf;
+  size_t i, icomponent;
+
+  if(!bsdf->ncomponents) return 0;
+  if(bsdf->ncomponents == 1) {
+    return ssf_bxdf_sample(bsdf->components[0], u, v, wi, N, wo);
+  }
 
-    probas[n] = reflectivities[n] / dirs[n][PDF];
-    probas_sum += probas[n];
-    ++n;
+  /* Compute the normalization factor of the BxDF weights */
+  probas_sum = 0;
+  FOR_EACH(icomponent, 0, bsdf->ncomponents) {
+    probas_sum += bsdf->weights[icomponent];
   }
 
-  if(!n) { /* No valid BxDF to sample */
-    d4_splat(wo, 0.0);
-    return 0.0;
+  /* Compute the components probabilities and the cumulative */
+  cumul[0] = probas[0] = bsdf->weights[0] / probas_sum;
+  FOR_EACH(icomponent, 0, bsdf->ncomponents) {
+    probas[icomponent] = bsdf->weights[icomponent] / probas_sum;
+    cumul[icomponent] = cumul[icomponent-1] + probas[icomponent];
   }
+  /* Sample a component */
+  FOR_EACH(icomponent, 0, bsdf->ncomponents-1) {
+    if(u >= cumul[icomponent]) break;
+  }
+
+  /* Rescale the random number to reuse it (NOTE why?) */
+  w = (u - cumul[icomponent] ) / (cumul[icomponent+1] - cumul[icomponent]);
 
-  /* Normalize the probability distribution */
-  FOR_EACH(i, 0, n) probas[i] /= probas_sum;
+  /* Sample the component */
+  reflectivity = ssf_bxdf_sample(bsdf->components[icomponent], w, v, wi, N, wo);
+  if(reflectivity == 0) return 0;
 
-  /* Compute the cumulative */
-  cumul[0] = probas[0];
-  cumul[n-1] = 1.f;
-  FOR_EACH(i, 1, n-1) cumul[i] = cumul[i-1] + probas[i];
+  pdf = wo[3] * probas[icomponent];
+  reflectivity *= pdf;
 
-  /* Finally sample the distribution */
-  FOR_EACH(i, 0, n-1) if(u <= cumul[i]) break;
-  d4_set(wo, dirs[i]);
-  return reflectivities[i];
+  /* Add the contribution of the others components */
+  FOR_EACH(i, 0, bsdf->ncomponents) {
+    if(i == icomponent) continue;
+    pdf += ssf_bxdf_pdf(bsdf->components[i], wi, wo) * probas[i];
+    reflectivity += ssf_bxdf_eval(bsdf->components[i], wi, wo) * probas[i];
+  }
+  wo[3] = pdf;
+  return reflectivity / pdf;
 }