star-sp

commit f79e41908caa31a4d83af8095dced7762f05627d
parent 8b7f4dddb43c85722b51cf454d32cc157ee5c0a4
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 10 Jul 2018 10:00:36 +0200

Add the uniform sphere cap random variate

Diffstat:
M
src/ssp.h
 | 
70
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
M
src/ssp_ran.c
 | 
58
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++-

2 files changed, 126 insertions(+), 2 deletions(-)
diff --git a/src/ssp.h b/src/ssp.h
@@ -74,7 +74,7 @@ struct ssp_rng_type {
   res_T (*set)(void* state, const uint64_t seed);
 
   uint64_t (*get)(void* state);
-  res_T(*discard)(void* state, uint64_t n);
+  res_T (*discard)(void* state, uint64_t n);
   res_T (*read)(void* state, FILE* file);
   res_T (*write)(const void* state, FILE* file);
   double (*entropy)(const void* state);
@@ -927,6 +927,74 @@ ssp_ran_uniform_disk_float
   return f33_mulf3(pt, f33_basis(basis, up), sample_local);
 }
 
+/*******************************************************************************
+ * Uniform distribution of a point ont a sphere cap
+ * pdf = 1/(2*PI*(1-cos(theta_max))
+ ******************************************************************************/
+/* Uniform sampling of unit sphere cap centered in zero whose up direction
+ * is implicity the Z axis. */
+SSP_API double*
+ssp_ran_sphere_cap_uniform_local
+  (struct ssp_rng* rng,
+   /* In [-1, 1]. Height where the sphere cap begins. It equal cos(theta_max)
+    * where theta_max is the aperture angle from the up axis */
+   const double height,
+   double sample[3],
+   double* pdf); /* Can be NULL => pdf not computed */
+
+SSP_API float*
+ssp_ran_sphere_cap_uniform_float_local
+  (struct ssp_rng* rng,
+   const float height, /* In [-1, 1]. Height of the sphere cap. */
+   float sample[3],
+   float* pdf); /* Can be NULL => pdf not computed */
+
+static INLINE double
+ssp_ran_sphere_cap_uniform_pdf(const double height)
+{
+  ASSERT(height <= 1.0 && height >= -1.0);
+  return height == 1.0 ? INF : 1.0/(2.0*PI*(1-height));
+}
+
+static INLINE float
+ssp_ran_sphere_cap_uniform_float_pdf(const float height)
+{
+  ASSERT(height <= 1.f && height >= -1.f);
+  return height == 1.f ? (float)INF : 1.f/(2.f*(float)PI*(1.f-height));
+}
+
+/* Uniform sampling of unit sphere cap centered in zero whose up direction
+ * is explicitly defined */
+static INLINE double*
+ssp_ran_sphere_cap_uniform
+  (struct ssp_rng* rng,
+   const double up[3],
+   const double height,
+   double sample[3],
+   double* pdf) /* Can be NULL */
+{
+  double sample_local[3];
+  double basis[9];
+  ASSERT(up);
+  ssp_ran_sphere_cap_uniform_local(rng, height, sample_local, pdf);
+  return d33_muld3(sample, d33_basis(basis, up), sample_local);
+}
+
+static  INLINE float*
+ssp_ran_sphere_cap_uniform_float
+  (struct ssp_rng* rng,
+   const float up[3],
+   const float height,
+   float sample[3],
+   float* pdf)
+{
+  float sample_local[3];
+  float basis[9];
+  ASSERT(up);
+  ssp_ran_sphere_cap_uniform_float_local(rng, height, sample_local, pdf);
+  return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+}
+
 END_DECLS
 
 #endif /* SSP_H */
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -219,7 +219,7 @@ ssp_ran_circle_uniform
   return sample;
 }
 
-float* 
+float*
 ssp_ran_circle_uniform_float
   (struct ssp_rng* rng,
    float sample[2],
@@ -521,3 +521,59 @@ ssp_ran_uniform_disk_float_local
   if (pdf) *pdf = 1 / (radius * radius);
   return pt;
 }
+
+double*
+ssp_ran_sphere_cap_uniform_local
+  (struct ssp_rng* rng,
+   const double height,
+   double pt[3],
+   double* pdf)
+{
+  ASSERT(rng && pt && (height >= -1 || height <= 1));
+
+  if(height == 1) {
+    d3(pt, 0, 0, 1);
+    if(pdf) *pdf = INF;
+  } else if(height == -1) {
+    ssp_ran_sphere_uniform(rng, pt, pdf);
+  } else {
+    const double cos_aperture = height;
+    double phi, cos_theta, sin_theta;
+    phi = ssp_rng_uniform_double(rng, 0, 2.0*PI);
+    cos_theta = ssp_rng_uniform_double(rng, cos_aperture, 1);
+    sin_theta = cos2sin(cos_theta);
+    pt[0] = cos(phi) * sin_theta;
+    pt[1] = sin(phi) * sin_theta;
+    pt[2] = cos_theta;
+    if(pdf) *pdf = 1.0/(2.0*PI*(1.0-cos_aperture));
+  }
+  return pt;
+}
+
+float*
+ssp_ran_sphere_cap_uniform_float_local
+  (struct ssp_rng* rng,
+   const float height,
+   float pt[3],
+   float* pdf)
+{
+  ASSERT(rng && pt && (height >= -1 || height <= 1));
+  if(height == 1) {
+    f3(pt, 0, 0, 1);
+    if(pdf) *pdf = (float)INF;
+  } else if(height == -1) {
+    ssp_ran_sphere_uniform_float(rng, pt, pdf);
+  } else {
+    const float cos_aperture = height;
+    float phi, cos_theta, sin_theta;
+    phi = ssp_rng_uniform_float(rng, 0, 2.f*(float)PI);
+    cos_theta = ssp_rng_uniform_float(rng, cos_aperture, 1);
+    sin_theta = (float)cos2sin(cos_theta);
+    pt[0] = cosf(phi) * sin_theta;
+    pt[1] = sinf(phi) * sin_theta;
+    pt[2] = cos_theta;
+    if(pdf) *pdf = 1.f/(2.f*(float)PI*(1.f-cos_aperture));
+  }
+  return pt;
+}
+