intern/cycles/kernel/light/point.h

@@ -13,27 +13,27 @@ ccl_device_inline bool point_light_sample(const ccl_global KernelLight *klight,
                                           const float3 P,
                                           ccl_private LightSample *ls)
 {
-  float3 center = klight->co;
-  float radius = klight->spot.radius;
-  /* disk oriented normal */
-  const float3 lightN = normalize(P - center);
-  ls->P = center;
+  float dist;
+  const float3 lightN = normalize_len(P - klight->co, &dist);
 
-  if (radius > 0.0f) {
-    ls->P += disk_light_sample(lightN, rand) * radius;
-  }
-  ls->pdf = klight->spot.invarea;
+  const float one_minus_cos_angle = tan_to_one_minus_cos(klight->spot.radius / dist);
 
-  ls->D = normalize_len(ls->P - P, &ls->t);
-  /* we set the light normal to the outgoing direction to support texturing */
-  ls->Ng = -ls->D;
+  /* We set the light normal to the outgoing direction to support texturing. */
+  sample_uniform_cone_concentric(lightN, one_minus_cos_angle, rand, &ls->Ng, &ls->pdf);
+  ls->D = -ls->Ng;
+  ls->t = dist / dot(ls->Ng, lightN);
+  ls->P = P + ls->D * ls->t;
 
+  /* TODO: change invarea to that of a sphere. */
   ls->eval_fac = M_1_PI_F * 0.25f * klight->spot.invarea;
-
-  float2 uv = map_to_sphere(ls->Ng);
+  if (one_minus_cos_angle == 0) {
+    /* Use intensity instead of radiance for point light. */
+    ls->eval_fac /= sqr(ls->t);
+  }
+  const float2 uv = map_to_sphere(ls->Ng);
   ls->u = uv.x;
   ls->v = uv.y;
-  ls->pdf *= lamp_light_pdf(lightN, -ls->D, ls->t);
+
   return true;
 }
 
@@ -44,14 +44,22 @@ ccl_device_forceinline void point_light_mnee_sample_update(const ccl_global Kern
   ls->D = normalize_len(ls->P - P, &ls->t);
   ls->Ng = -ls->D;
 
-  float2 uv = map_to_sphere(ls->Ng);
+  const float2 uv = map_to_sphere(ls->Ng);
   ls->u = uv.x;
   ls->v = uv.y;
 
-  float invarea = klight->spot.invarea;
-  ls->eval_fac = (0.25f * M_1_PI_F) * invarea;
-  /* NOTE : preserve pdf in area measure. */
-  ls->pdf = invarea;
+  const float radius = klight->spot.radius;
+
+  if (radius > 0) {
+    const float dist = len(P - klight->co);
+    const float one_minus_cos_angle = tan_to_one_minus_cos(radius / dist);
+    /* NOTE : preserve pdf in area measure. */
+    ls->pdf = M_1_2PI_F / one_minus_cos_angle * dist / powf(ls->t, 3.0f);
+  }
+  else {
+    ls->eval_fac = M_1_PI_F * 0.25f * klight->spot.invarea;
+    /* PDF does not change. */
+  }
 }
 
 ccl_device_inline bool point_light_intersect(const ccl_global KernelLight *klight,
@@ -65,7 +73,7 @@ ccl_device_inline bool point_light_intersect(const ccl_global KernelLight *kligh
     return false;
   }
 
-  /* disk oriented normal */
+  /* Disk oriented normal. */
   const float3 lightN = normalize(ray->P - lightP);
   float3 P;
   return ray_disk_intersect(ray->P, ray->D, ray->tmin, ray->tmax, lightP, lightN, radius, &P, t);
@@ -78,22 +86,24 @@ ccl_device_inline bool point_light_sample_from_intersection(
     const float3 ray_D,
     ccl_private LightSample *ccl_restrict ls)
 {
-  const float3 lighN = normalize(ray_P - klight->co);
-
   /* We set the light normal to the outgoing direction to support texturing. */
   ls->Ng = -ls->D;
 
-  float invarea = klight->spot.invarea;
-  ls->eval_fac = (0.25f * M_1_PI_F) * invarea;
-  ls->pdf = invarea;
+  ls->eval_fac = (0.25f * M_1_PI_F) * klight->spot.invarea;
 
-  float2 uv = map_to_sphere(ls->Ng);
+  const float2 uv = map_to_sphere(ls->Ng);
   ls->u = uv.x;
   ls->v = uv.y;
 
-  /* compute pdf */
   if (ls->t != FLT_MAX) {
-    ls->pdf *= lamp_light_pdf(lighN, -ls->D, ls->t);
+    const float radius = klight->spot.radius;
+    if (radius > 0) {
+      ls->pdf = M_1_2PI_F / sin_to_one_minus_cos(radius / ls->t);
+    }
+    else {
+      ls->eval_fac /= sqr(ls->t);
+      ls->pdf = 1.0f;
+    }
   }
   else {
     ls->pdf = 0.f;

intern/cycles/kernel/sample/mapping.h

@@ -20,6 +20,29 @@ ccl_device void to_unit_disk(ccl_private float2 *rand)
   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)
+{
+  float phi, r;
+  float a = 2.0f * rand.x - 1.0f;
+  float b = 2.0f * rand.y - 1.0f;
+
+  if (a == 0.0f && b == 0.0f) {
+    return zero_float2();
+  }
+  else if (a * a > b * b) {
+    r = a;
+    phi = M_PI_4_F * (b / a);
+  }
+  else {
+    r = b;
+    phi = M_PI_2_F - M_PI_4_F * (a / b);
+  }
+
+  return make_float2(r * cosf(phi), r * sinf(phi));
+}
+
 /* return an orthogonal tangent and bitangent given a normal and tangent that
  * may not be exactly orthogonal */
 ccl_device void make_orthonormals_tangent(const float3 N,
@@ -92,6 +115,37 @@ ccl_device_inline float pdf_uniform_cone(const float3 N, float3 D, float angle)
   return 0.0f;
 }
 
+/* Sample direction uniformly distributed in a cone around `N`. Using concentric mapping to better
+ * preserve stratification. */
+ccl_device_inline void sample_uniform_cone_concentric(const float3 N,
+                                                      const float one_minus_cos_angle,
+                                                      const float2 rand,
+                                                      ccl_private float3 *wo,
+                                                      ccl_private float *pdf)
+{
+  if (one_minus_cos_angle > 0) {
+    /* Map random number from 2D to 1D. */
+    float2 xy = concentric_sample_disk(rand);
+    const float r2 = len_squared(xy);
+
+    /* Equivalent to `mix(cos_angle, 1.0f, 1.0f - r2)` */
+    const float cos_theta = 1.0f - r2 * one_minus_cos_angle;
+
+    /* Equivalent to `xy *= sin_theta / sqrt(r2); */
+    xy *= safe_sqrtf(one_minus_cos_angle * (2.0f - one_minus_cos_angle * r2));
+
+    float3 T, B;
+    make_orthonormals(N, &T, &B);
+
+    *wo = xy.x * T + xy.y * B + cos_theta * N;
+    *pdf = M_1_2PI_F / one_minus_cos_angle;
+  }
+  else {
+    *wo = N;
+    *pdf = 1.0f;
+  }
+}
+
 /* sample uniform point on the surface of a sphere */
 ccl_device float3 sample_uniform_sphere(const float2 rand)
 {
@@ -104,29 +158,6 @@ ccl_device float3 sample_uniform_sphere(const float2 rand)
   return make_float3(x, y, z);
 }
 
-/* distribute uniform xy 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)
-{
-  float phi, r;
-  float a = 2.0f * rand.x - 1.0f;
-  float b = 2.0f * rand.y - 1.0f;
-
-  if (a == 0.0f && b == 0.0f) {
-    return zero_float2();
-  }
-  else if (a * a > b * b) {
-    r = a;
-    phi = M_PI_4_F * (b / a);
-  }
-  else {
-    r = b;
-    phi = M_PI_2_F - M_PI_4_F * (a / b);
-  }
-
-  return make_float2(r * cosf(phi), r * sinf(phi));
-}
-
 /* sample point in unit polygon with given number of corners and rotation */
 ccl_device float2 regular_polygon_sample(float corners, float rotation, const float2 rand)
 {

intern/cycles/util/math.h

@@ -747,6 +747,20 @@ ccl_device_inline float cos_from_sin(const float s)
   return safe_sqrtf(1.0f - sqr(s));
 }
 
+ccl_device_inline float tan_to_one_minus_cos(const float t)
+{
+  /* Using second-order Taylor expansion at small angles for better accuracy. */
+  const float t2 = sqr(t);
+  return t > 0.02f ? 1.0f - inversesqrtf(t2 + 1.0f) : 0.5f * t2;
+}
+
+ccl_device_inline float sin_to_one_minus_cos(const float s)
+{
+  /* Using second-order Taylor expansion at small angles for better accuracy. */
+  const float s2 = sqr(s);
+  return s > 0.02f ? 1.0f - safe_sqrtf(1.0f - s2) : 0.5f * s2;
+}
+
 ccl_device_inline float pow20(float a)
 {
   return sqr(sqr(sqr(sqr(a)) * a));