star-phor

commit bc2abf06454e245f12a300a471eb982b61e7e542
parent 8b22ce94b0c23100006ffff1a25ec132ef8b92ed
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date:   Fri, 23 Jan 2026 11:33:43 +0100

Implement a function to sample a direction according to cos^n radiance

Sample a source direction for a ray assuming a collimation degree n of
the source. The direction is sampled according to the solid angle PDF:

  pdf(Omega) = (n+2) / (2 * pi) * cos^(n+1)(theta)

using the inversion method on the cumulative distribution function
(CDF). The model assumes azimuthal independence: phi is sampled
uniformly in [0, 2 * pi).

The function returns the sampled direction in the global coordinate
system and optionally the corresponding probability density function
value.

Diffstat:
M
src/sphor_ran_source.c
 | 
42
+++++++++++++++++++++++++++++++++++++++++-

1 file changed, 41 insertions(+), 1 deletion(-)
diff --git a/src/sphor_ran_source.c b/src/sphor_ran_source.c
@@ -30,11 +30,50 @@
 
 #include <rsys/double3.h>
 #include <rsys/float3.h>
+#include <rsys/math.h>
 #include <rsys/rsys.h>
 #include <star/s3d.h>
 #include <star/sphin.h>
 #include <star/ssp.h>
 
+#include <math.h>
+
+static double*
+ran_hemisphere_cos_pow_n
+  (struct ssp_rng* rng,
+   double n,
+   double normal[3],
+   double sample[3],
+   double* pdf)
+{
+  double phi = 0;
+  double tmp = 0;
+  double cos_theta = 0;
+  double sin_theta = 0;
+  double basis[9] = {0};
+  double sample_local[3] = {0};
+
+  ASSERT(NULL != rng);
+  ASSERT(NULL != normal);
+  ASSERT(NULL != sample);
+  ASSERT(d3_is_normalized(normal));
+  ASSERT(0 <= n);
+
+  phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
+  tmp = ssp_rng_canonical(rng);
+
+  cos_theta = pow(tmp, 1 / (n + 2));
+  sin_theta = sqrt(1- cos_theta * cos_theta);
+
+  if(pdf) *pdf = (n+2) * pow(cos_theta, n + 1) / (2 * PI);
+
+  sample_local[0] = sin_theta * cos(phi);
+  sample_local[1] = sin_theta * sin(phi);
+  sample_local[2] = cos_theta;
+
+  return d33_muld3(sample, d33_basis(basis, normal), sample_local);
+}
+
 /*******************************************************************************
  * Local functions
  ******************************************************************************/
@@ -67,7 +106,8 @@ source_surface_sample_direction
       d3_set(sample, normal);
       break;
     case SPHIN_SOURCE_DIRECTION_COS_POW_N:
-      /* TODO Sample direction with collimation degree */
+      ran_hemisphere_cos_pow_n
+        (rng, src_dir_dist.cos_pow_n.collimation_degree, normal, sample, NULL);
       break;
     case SPHIN_SOURCE_DIRECTION_ISOTROPIC:
       ssp_ran_hemisphere_cos(rng, normal, sample, NULL);