Cycles: use low-distortion mapping when sampling cone and hemisphere

based on concentric disk mapping.
Concentric disk mapping was already present, but not used everywhere.
Now `sample_cos_hemisphere()`, `sample_uniform_hemisphere()`, and
`sample_uniform_cone()` use concentric disk mapping.
This changes the noise in many test images.

Pull Request: https://projects.blender.org/blender/blender/pulls/109774

intern/cycles/kernel/camera/camera.h

@@ -20,7 +20,7 @@ ccl_device float2 camera_sample_aperture(ccl_constant KernelCamera *cam, const f
 
   if (blades == 0.0f) {
     /* sample disk */
-    bokeh = concentric_sample_disk(rand);
+    bokeh = sample_uniform_disk(rand);
   }
   else {
     /* sample polygon */

intern/cycles/kernel/closure/bsdf_microfacet.h

@@ -195,7 +195,7 @@ ccl_device_forceinline float3 microfacet_ggx_sample_vndf(const float3 wi,
   }
 
   /* Section 4.2: Parameterization of the projected area. */
-  float2 t = concentric_sample_disk(rand);
+  float2 t = sample_uniform_disk(rand);
   t.y = mix(safe_sqrtf(1.0f - sqr(t.x)), t.y, 0.5f * (1.0f + wi_.z));
 
   /* Section 4.3: Reprojection onto hemisphere. */

intern/cycles/kernel/closure/bsdf_sheen.h

@@ -82,7 +82,7 @@ ccl_device int bsdf_sheen_sample(ccl_private const ShaderClosure *sc,
   const float3 N = bsdf->N, T = bsdf->T, B = bsdf->B;
   float a = bsdf->transformA, b = bsdf->transformB;
 
-  float2 disk = concentric_sample_disk(rand);
+  float2 disk = sample_uniform_disk(rand);
   float diskZ = safe_sqrtf(1.0f - dot(disk, disk));
   float3 localO = normalize(make_float3((disk.x - diskZ * b) / a, disk.y / a, diskZ));
 

intern/cycles/kernel/closure/bsdf_toon.h

@@ -90,7 +90,8 @@ ccl_device int bsdf_diffuse_toon_sample(ccl_private const ShaderClosure *sc,
   float angle = sample_angle * rand.x;
 
   if (sample_angle > 0.0f) {
-    sample_uniform_cone(bsdf->N, sample_angle, rand, wo, pdf);
+    float unused;
+    *wo = sample_uniform_cone(bsdf->N, one_minus_cos(sample_angle), rand, &unused, pdf);
 
     if (dot(Ng, *wo) > 0.0f) {
       *eval = make_spectrum(*pdf * bsdf_toon_get_intensity(max_angle, smooth, angle));
@@ -167,7 +168,8 @@ ccl_device int bsdf_glossy_toon_sample(ccl_private const ShaderClosure *sc,
     float sample_angle = bsdf_toon_get_sample_angle(max_angle, smooth);
     float angle = sample_angle * rand.x;
 
-    sample_uniform_cone(R, sample_angle, rand, wo, pdf);
+    float unused;
+    *wo = sample_uniform_cone(R, one_minus_cos(sample_angle), rand, &unused, pdf);
 
     if (dot(Ng, *wo) > 0.0f) {
       float cosNO = dot(bsdf->N, *wo);

intern/cycles/kernel/light/background.h

@@ -277,11 +277,10 @@ ccl_device_inline float3 background_sun_sample(KernelGlobals kg,
                                                float2 rand,
                                                ccl_private float *pdf)
 {
-  float3 D;
   const float3 N = float4_to_float3(kernel_data.background.sun);
   const float angle = kernel_data.background.sun.w;
-  sample_uniform_cone(N, angle, rand, &D, pdf);
-  return D;
+  float unused;
+  return sample_uniform_cone(N, one_minus_cos(angle), rand, &unused, pdf);
 }
 
 ccl_device_inline float background_sun_pdf(KernelGlobals kg, float3 D)

intern/cycles/kernel/light/common.h

@@ -31,7 +31,7 @@ typedef struct LightSample {
 
 ccl_device_inline float3 ellipse_sample(float3 ru, float3 rv, float2 rand)
 {
-  const float2 uv = concentric_sample_disk(rand);
+  const float2 uv = sample_uniform_disk(rand);
   return ru * uv.x + rv * uv.y;
 }
 

intern/cycles/kernel/light/distant.h

@@ -36,8 +36,8 @@ ccl_device_inline bool distant_light_sample(const ccl_global KernelLight *klight
                                             ccl_private LightSample *ls)
 {
   float unused;
-  sample_uniform_cone_concentric(
-      klight->co, klight->distant.one_minus_cosangle, rand, &unused, &ls->Ng, &ls->pdf);
+  ls->Ng = sample_uniform_cone(
+      klight->co, klight->distant.one_minus_cosangle, rand, &unused, &ls->pdf);
 
   ls->P = ls->Ng;
   ls->D = -ls->Ng;

intern/cycles/kernel/light/point.h

@@ -25,7 +25,7 @@ ccl_device_inline bool point_light_sample(const ccl_global KernelLight *klight,
   float cos_theta;
   if (d_sq > r_sq) {
     const float one_minus_cos = sin_sqr_to_one_minus_cos(r_sq / d_sq);
-    sample_uniform_cone_concentric(-lightN, one_minus_cos, rand, &cos_theta, &ls->D, &ls->pdf);
+    ls->D = sample_uniform_cone(-lightN, one_minus_cos, rand, &cos_theta, &ls->pdf);
   }
   else {
     const bool has_transmission = (shader_flags & SD_BSDF_HAS_TRANSMISSION);

intern/cycles/kernel/light/spot.h

@@ -61,13 +61,12 @@ ccl_device_inline bool spot_light_sample(const ccl_global KernelLight *klight,
 
     if (in_volume_segment || one_minus_cos_half_angle < one_minus_cos_half_spot_spread) {
       /* Sample visible part of the sphere. */
-      sample_uniform_cone_concentric(
-          -lightN, one_minus_cos_half_angle, rand, &cos_theta, &ls->D, &ls->pdf);
+      ls->D = sample_uniform_cone(-lightN, one_minus_cos_half_angle, rand, &cos_theta, &ls->pdf);
     }
     else {
       /* Sample spread cone. */
-      sample_uniform_cone_concentric(
-          -klight->spot.dir, one_minus_cos_half_spot_spread, rand, &cos_theta, &ls->D, &ls->pdf);
+      ls->D = sample_uniform_cone(
+          -klight->spot.dir, one_minus_cos_half_spot_spread, rand, &cos_theta, &ls->pdf);
 
       if (!ray_sphere_intersect(P, ls->D, 0.0f, FLT_MAX, center, radius, &ls->P, &ls->t)) {
         /* Sampled direction does not intersect with the light. */

intern/cycles/kernel/sample/mapping.h

@@ -9,19 +9,9 @@
 
 CCL_NAMESPACE_BEGIN
 
-/* distribute uniform xy on [0,1] over unit disk [-1,1] */
-ccl_device void to_unit_disk(ccl_private float2 *rand)
-{
-  float phi = M_2PI_F * rand->x;
-  float r = sqrtf(rand->y);
-
-  rand->x = r * cosf(phi);
-  rand->y = r * sinf(phi);
-}
-
 /* Distribute 2D uniform random samples on [0, 1] over unit disk [-1, 1], with concentric mapping
  * to better preserve stratification for some RNG sequences. */
-ccl_device float2 concentric_sample_disk(const float2 rand)
+ccl_device float2 sample_uniform_disk(const float2 rand)
 {
   float phi, r;
   float a = 2.0f * rand.x - 1.0f;
@@ -55,11 +45,11 @@ ccl_device void make_orthonormals_tangent(const float3 N,
 
 /* sample direction with cosine weighted distributed in hemisphere */
 ccl_device_inline void sample_cos_hemisphere(const float3 N,
-                                             float2 rand,
+                                             float2 rand_in,
                                              ccl_private float3 *wo,
                                              ccl_private float *pdf)
 {
-  to_unit_disk(&rand);
+  float2 rand = sample_uniform_disk(rand_in);
   float costheta = safe_sqrtf(1.0f - len_squared(rand));
 
   float3 T, B;
@@ -80,42 +70,23 @@ ccl_device_inline void sample_uniform_hemisphere(const float3 N,
                                                  ccl_private float3 *wo,
                                                  ccl_private float *pdf)
 {
-  float z = rand.x;
-  float r = sin_from_cos(z);
-  float phi = M_2PI_F * rand.y;
-  float x = r * cosf(phi);
-  float y = r * sinf(phi);
+  float2 xy = sample_uniform_disk(rand);
+  float z = 1.0f - len_squared(xy);
+
+  xy *= safe_sqrtf(z + 1.0f);
 
   float3 T, B;
   make_orthonormals(N, &T, &B);
-  *wo = x * T + y * B + z * N;
-  *pdf = 0.5f * M_1_PI_F;
-}
 
-/* sample direction uniformly distributed in cone */
-ccl_device_inline void sample_uniform_cone(
-    const float3 N, float angle, const float2 rand, ccl_private float3 *wo, ccl_private float *pdf)
-{
-  const float cosThetaMin = cosf(angle);
-  const float cosTheta = mix(cosThetaMin, 1.0f, rand.x);
-  const float sinTheta = sin_from_cos(cosTheta);
-  const float phi = M_2PI_F * rand.y;
-  const float x = sinTheta * cosf(phi);
-  const float y = sinTheta * sinf(phi);
-  const float z = cosTheta;
-
-  float3 T, B;
-  make_orthonormals(N, &T, &B);
-  *wo = x * T + y * B + z * N;
-  *pdf = M_1_2PI_F / (1.0f - cosThetaMin);
+  *wo = xy.x * T + xy.y * B + z * N;
+  *pdf = M_1_2PI_F;
 }
 
 ccl_device_inline float pdf_uniform_cone(const float3 N, float3 D, float angle)
 {
-  float zMin = cosf(angle);
   float z = precise_angle(N, D);
   if (z < angle) {
-    return M_1_2PI_F / (1.0f - zMin);
+    return M_1_2PI_F / one_minus_cos(angle);
   }
   return 0.0f;
 }
@@ -125,12 +96,11 @@ ccl_device_inline float pdf_uniform_cone(const float3 N, float3 D, float angle)
  * `cos_theta`.
  * Pass `1 - cos(angle)` as argument instead of `angle` to alleviate precision issues at small
  * angles (see sphere light for reference). */
-ccl_device_inline void sample_uniform_cone_concentric(const float3 N,
-                                                      const float one_minus_cos_angle,
-                                                      const float2 rand,
-                                                      ccl_private float *cos_theta,
-                                                      ccl_private float3 *wo,
-                                                      ccl_private float *pdf)
+ccl_device_inline float3 sample_uniform_cone(const float3 N,
+                                             const float one_minus_cos_angle,
+                                             const float2 rand,
+                                             ccl_private float *cos_theta,
+                                             ccl_private float *pdf)
 {
   if (one_minus_cos_angle > 0) {
     /* Remap radius to get a uniform distribution w.r.t. solid angle on the cone.
@@ -155,7 +125,7 @@ ccl_device_inline void sample_uniform_cone_concentric(const float3 N,
      * r_cone: `r_cone(r_disk) = r_cone(rand(r_disk)) = sin_from_cos(mix(cos_angle, 1, r_disk^2))`.
      * In practice, we need to replace `rand` with `1 - rand` to preserve the stratification,
      * but since it's uniform, that's fine. */
-    float2 xy = concentric_sample_disk(rand);
+    float2 xy = sample_uniform_disk(rand);
     const float r2 = len_squared(xy);
 
     /* Equivalent to `mix(cos_angle, 1.0f, 1.0f - r2)` */
@@ -164,17 +134,17 @@ ccl_device_inline void sample_uniform_cone_concentric(const float3 N,
     /* Remap disk radius to cone radius, equivalent to `xy *= sin_theta / sqrt(r2); */
     xy *= safe_sqrtf(one_minus_cos_angle * (2.0f - one_minus_cos_angle * r2));
 
+    *pdf = M_1_2PI_F / one_minus_cos_angle;
+
     float3 T, B;
     make_orthonormals(N, &T, &B);
+    return xy.x * T + xy.y * B + *cos_theta * N;
+  }
 
-    *wo = xy.x * T + xy.y * B + *cos_theta * N;
-    *pdf = M_1_2PI_F / one_minus_cos_angle;
-  }
-  else {
-    *cos_theta = 1.0f;
-    *wo = N;
-    *pdf = 1.0f;
-  }
+  *cos_theta = 1.0f;
+  *pdf = 1.0f;
+
+  return N;
 }
 
 /* sample uniform point on the surface of a sphere */

intern/cycles/kernel/svm/ao.h

@@ -53,7 +53,7 @@ ccl_device float svm_ao(
     const float2 rand_disk = path_branched_rng_2D(
         kg, &rng_state, sample, num_samples, PRNG_SURFACE_AO);
 
-    float2 d = concentric_sample_disk(rand_disk);
+    float2 d = sample_uniform_disk(rand_disk);
     float3 D = make_float3(d.x, d.y, safe_sqrtf(1.0f - dot(d, d)));
 
     /* Create ray. */

intern/cycles/util/math.h

@@ -758,6 +758,12 @@ ccl_device_inline float sin_sqr_to_one_minus_cos(const float s_sq)
   return s_sq > 0.0004f ? 1.0f - safe_sqrtf(1.0f - s_sq) : 0.5f * s_sq;
 }
 
+ccl_device_inline float one_minus_cos(const float angle)
+{
+  /* Using second-order Taylor expansion at small angles for better accuracy. */
+  return angle > 0.02f ? 1.0f - cosf(angle) : 0.5f * sqr(angle);
+}
+
 ccl_device_inline float pow20(float a)
 {
   return sqr(sqr(sqr(sqr(a)) * a));