star-sp

commit a1eecf0fa952a1af2f98d4464e7a2b217d7cce6b
parent 6294e35f20351c7f58cc442cfdce8eedb0774000
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Tue,  9 Jul 2019 15:56:52 +0200

Fix HG pdf when |g|==1

Diffstat:
M
src/ssp.h
 | 
20
++++++++++++++++----
M
src/ssp_ran.c
 | 
48
++++++++++++++++++++++++++++++------------------

2 files changed, 46 insertions(+), 22 deletions(-)
diff --git a/src/ssp.h b/src/ssp.h
@@ -674,8 +674,14 @@ ssp_ran_sphere_hg
   double sample_local[3];
   double basis[9];
   ASSERT(-1 <= g && g <= +1 && rng && up && sample && d3_is_normalized(up));
-  ssp_ran_sphere_hg_local(rng, g, sample_local, pdf);
-  return d33_muld3(sample, d33_basis(basis, up), sample_local);
+  if (fabs(g) == 1) {
+    d3_muld(sample, up, g);
+    if(pdf) *pdf=INF;
+  } else {
+    ssp_ran_sphere_hg_local(rng, g, sample_local, pdf);
+    d33_muld3(sample, d33_basis(basis, up), sample_local);
+  }
+  return sample;
 }
 
 static INLINE float*
@@ -689,8 +695,14 @@ ssp_ran_sphere_hg_float
   float sample_local[3];
   float basis[9];
   ASSERT(-1 <= g && g <= +1 && rng && up && sample && f3_is_normalized(up));
-  ssp_ran_sphere_hg_float_local(rng, g, sample_local, pdf);
-  return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+  if(fabsf(g) == 1) {
+    f3_mulf(sample, up, g);
+    if(pdf) *pdf = (float)INF;
+  } else {
+    ssp_ran_sphere_hg_float_local(rng, g, sample_local, pdf);
+    f33_mulf3(sample, f33_basis(basis, up), sample_local);
+  }
+  return sample;
 }
 
 /* Return the probability distribution for a Henyey-Greenstein random
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -411,15 +411,20 @@ ssp_ran_sphere_hg_local
 {
   double phi, R, cos_theta, sin_theta;
   ASSERT(fabs(g) <= 1 && rng && sample);
-  phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
-  R = ssp_rng_canonical(rng);
-  cos_theta = 2 * R*(1 + g)*(1 + g)*(g*(R - 1) + 1)
-    / ((1 - g*(1 - 2 * R))*(1 - g*(1 - 2 * R))) - 1;
-  sin_theta = cos2sin(cos_theta);
-  sample[0] = cos(phi) * sin_theta;
-  sample[1] = sin(phi) * sin_theta;
-  sample[2] = cos_theta;
-  if (pdf) *pdf = ssp_ran_sphere_hg_local_pdf(g, sample);
+  if(fabs(g) == 1) {
+    d3(sample, 0, 0, g);
+    if(pdf) *pdf = INF;
+  } else {
+    phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
+    R = ssp_rng_canonical(rng);
+    cos_theta = 2 * R*(1 + g)*(1 + g)*(g*(R - 1) + 1)
+      / ((1 - g * (1 - 2 * R))*(1 - g * (1 - 2 * R))) - 1;
+    sin_theta = cos2sin(cos_theta);
+    sample[0] = cos(phi) * sin_theta;
+    sample[1] = sin(phi) * sin_theta;
+    sample[2] = cos_theta;
+    if(pdf) *pdf = ssp_ran_sphere_hg_local_pdf(g, sample);
+  }
   return sample;
 }
 
@@ -432,15 +437,20 @@ ssp_ran_sphere_hg_float_local
 {
   float phi, R, cos_theta, sin_theta;
   ASSERT(fabsf(g) <= 1 && rng && sample);
-  phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI);
-  R = ssp_rng_canonical_float(rng);
-  cos_theta = 2 * R*(1 + g)*(1 + g)*(g*(R - 1) + 1)
-    / ((1 - g*(1 - 2 * R))*(1 - g*(1 - 2 * R))) - 1;
-  sin_theta = cosf2sinf(cos_theta);
-  sample[0] = cosf(phi) * sin_theta;
-  sample[1] = sinf(phi) * sin_theta;
-  sample[2] = cos_theta;
-  if (pdf) *pdf = ssp_ran_sphere_hg_float_local_pdf(g, sample);
+  if(fabsf(g) == 1) {
+    f3(sample, 0, 0, g);
+    if(pdf) *pdf = (float)INF;
+  } else {
+    phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI);
+    R = ssp_rng_canonical_float(rng);
+    cos_theta = 2 * R*(1 + g)*(1 + g)*(g*(R - 1) + 1)
+      / ((1 - g * (1 - 2 * R))*(1 - g * (1 - 2 * R))) - 1;
+    sin_theta = cosf2sinf(cos_theta);
+    sample[0] = cosf(phi) * sin_theta;
+    sample[1] = sinf(phi) * sin_theta;
+    sample[2] = cos_theta;
+    if(pdf) *pdf = ssp_ran_sphere_hg_float_local_pdf(g, sample);
+  }
   return sample;
 }
 
@@ -453,6 +463,7 @@ ssp_ran_sphere_hg_pdf
   double epsilon_g, epsilon_mu;
   ASSERT(-1 <= g && g <= +1 &&
     up && sample && d3_is_normalized(up) && d3_is_normalized(sample));
+  if(fabs(g) == 1) return INF;
   if(g>0) {
     epsilon_g = 1 - g;
     epsilon_mu = 1 - d3_dot(sample, up);
@@ -474,6 +485,7 @@ ssp_ran_sphere_hg_float_pdf
   float epsilon_g, epsilon_mu;
   ASSERT(-1 <= g && g <= +1 &&
     up && sample && f3_is_normalized(up) && f3_is_normalized(sample));
+  if(fabsf(g) == 1) return (float)INF;
   if(g>0) {
     epsilon_g = 1 - g;
     epsilon_mu = 1 - f3_dot(sample, up);