Cycles: fix bugs in point and spot light multiple importance sampling

* Spot lights are now handled as disks aligned with the direction of the
  spotlight instead of view aligned disks.

* Point light is now handled separately from the spot light, to fix a case
  where multiple lights are intersected in a row. Before the origin of the
  ray was the previously intersected light and not the origin of the initial
  ray traced from the last surface/volume interaction.

This makes both strategies in multiple importance sampling converge to the same
result. It changes the render results in some scenes, for example the junkshop
scene where there are large point lights overlapping scene geometry and each
other.

Differential Revision: https://developer.blender.org/D13233

3 files changed, 119 insertions, 26 deletions

diff --git a/intern/cycles/kernel/integrator/intersect_closest.h b/intern/cycles/kernel/integrator/intersect_closest.h
index 366bfba7aca..df710dc1d82 100644
--- a/intern/cycles/kernel/integrator/intersect_closest.h
+++ b/intern/cycles/kernel/integrator/intersect_closest.h
@@ -341,9 +341,8 @@ ccl_device void integrator_intersect_closest(KernelGlobals kg,
      * these in the path_state_init. */
     const int last_type = INTEGRATOR_STATE(state, isect, type);
     const uint32_t path_flag = INTEGRATOR_STATE(state, path, flag);
-
     hit = lights_intersect(
-              kg, &ray, &isect, last_isect_prim, last_isect_object, last_type, path_flag) ||
+              kg, state, &ray, &isect, last_isect_prim, last_isect_object, last_type, path_flag) ||
           hit;
   }
 
diff --git a/intern/cycles/kernel/light/light.h b/intern/cycles/kernel/light/light.h
index e0a9f1c57f5..f47e8fbcf82 100644
--- a/intern/cycles/kernel/light/light.h
+++ b/intern/cycles/kernel/light/light.h
@@ -112,35 +112,58 @@ ccl_device_inline bool light_sample(KernelGlobals kg,
   else {
     ls->P = make_float3(klight->co[0], klight->co[1], klight->co[2]);
 
-    if (type == LIGHT_POINT || type == LIGHT_SPOT) {
+    if (type == LIGHT_SPOT) {
+      ls->Ng = make_float3(klight->spot.dir[0], klight->spot.dir[1], klight->spot.dir[2]);
       float radius = klight->spot.radius;
 
       if (radius > 0.0f)
         /* sphere light */
-        ls->P += sphere_light_sample(P, ls->P, radius, randu, randv);
+        ls->P += disk_light_sample(ls->Ng, randu, randv) * radius;
 
       ls->D = normalize_len(ls->P - P, &ls->t);
-      ls->Ng = -ls->D;
 
       float invarea = klight->spot.invarea;
       ls->eval_fac = (0.25f * M_1_PI_F) * invarea;
       ls->pdf = invarea;
 
-      if (type == LIGHT_SPOT) {
-        /* spot light attenuation */
-        float3 dir = make_float3(klight->spot.dir[0], klight->spot.dir[1], klight->spot.dir[2]);
-        ls->eval_fac *= spot_light_attenuation(
-            dir, klight->spot.spot_angle, klight->spot.spot_smooth, ls->Ng);
-        if (!in_volume_segment && ls->eval_fac == 0.0f) {
-          return false;
-        }
+      /* spot light attenuation */
+      ls->eval_fac *= spot_light_attenuation(
+          ls->Ng, klight->spot.spot_angle, klight->spot.spot_smooth, -ls->D);
+      if (!in_volume_segment && ls->eval_fac == 0.0f) {
+        return false;
       }
+
       float2 uv = map_to_sphere(ls->Ng);
       ls->u = uv.x;
       ls->v = uv.y;
 
       ls->pdf *= lamp_light_pdf(kg, ls->Ng, -ls->D, ls->t);
     }
+    else if (type == LIGHT_POINT) {
+      float3 center = make_float3(klight->co[0], klight->co[1], klight->co[2]);
+      float radius = klight->spot.radius;
+      ls->P = center;
+      float pdf = 1.0;
+
+      if (radius > 0.0f) {
+        ls->Ng = normalize(P - center);
+        ls->P += disk_light_sample(ls->Ng, randu, randv) * radius;
+        pdf = klight->spot.invarea;
+        ls->D = normalize_len(ls->P - P, &ls->t);
+      }
+      else {
+        ls->Ng = normalize(P - center);
+      }
+
+      ls->D = normalize_len(ls->P - P, &ls->t);
+      ls->pdf = pdf;
+      ls->eval_fac = M_1_PI_F * 0.25f * klight->spot.invarea;
+
+      float2 uv = map_to_sphere(ls->Ng);
+      ls->u = uv.x;
+      ls->v = uv.y;
+      ls->pdf *= lamp_light_pdf(kg, ls->Ng, -ls->D, ls->t);
+    }
     else {
       /* area light */
       float3 axisu = make_float3(
@@ -207,6 +230,7 @@ ccl_device_inline bool light_sample(KernelGlobals kg,
 }
 
 ccl_device bool lights_intersect(KernelGlobals kg,
+                                 IntegratorState state,
                                  ccl_private const Ray *ccl_restrict ray,
                                  ccl_private Intersection *ccl_restrict isect,
                                  const int last_prim,
@@ -237,8 +261,31 @@ ccl_device bool lights_intersect(KernelGlobals kg,
     LightType type = (LightType)klight->type;
     float t = 0.0f, u = 0.0f, v = 0.0f;
 
-    if (type == LIGHT_POINT || type == LIGHT_SPOT) {
-      /* Sphere light. */
+    if (type == LIGHT_SPOT) {
+      /* Spot/Disk light. */
+      const float3 lightP = make_float3(klight->co[0], klight->co[1], klight->co[2]);
+      const float3 lightN = make_float3(
+          klight->spot.dir[0], klight->spot.dir[1], klight->spot.dir[2]);
+      const float radius = klight->spot.radius;
+      if (radius == 0.0f) {
+        continue;
+      }
+
+      /* One sided. */
+      if (dot(ray->D, lightN) >= 0.0f) {
+        continue;
+      }
+
+      float3 P;
+      if (!ray_disk_intersect(ray->P, ray->D, ray->t, lightP, lightN, radius, &P, &t)) {
+        continue;
+      }
+    }
+    else if (type == LIGHT_POINT) {
+      /* Sphere light (aka, aligned disk light). */
+      const float mis_ray_t = INTEGRATOR_STATE(state, path, mis_ray_t);
+      const float3 ray_P = ray->P - ray->D * mis_ray_t;
+
       const float3 lightP = make_float3(klight->co[0], klight->co[1], klight->co[2]);
       const float radius = klight->spot.radius;
       if (radius == 0.0f) {
@@ -246,7 +293,8 @@ ccl_device bool lights_intersect(KernelGlobals kg,
       }
 
       float3 P;
-      if (!ray_aligned_disk_intersect(ray->P, ray->D, ray->t, lightP, radius, &P, &t)) {
+      const float3 lsN = normalize(ray_P - lightP);
+      if (!ray_disk_intersect(ray->P, ray->D, ray->t, lightP, lsN, radius, &P, &t)) {
         continue;
       }
     }
@@ -378,23 +426,38 @@ ccl_device bool light_sample_from_intersection(KernelGlobals kg,
   ls->P = ray_P + ray_D * ls->t;
   ls->D = ray_D;
 
-  if (type == LIGHT_POINT || type == LIGHT_SPOT) {
-    ls->Ng = -ray_D;
+  if (type == LIGHT_SPOT) {
+    ls->Ng = make_float3(klight->spot.dir[0], klight->spot.dir[1], klight->spot.dir[2]);
 
     float invarea = klight->spot.invarea;
     ls->eval_fac = (0.25f * M_1_PI_F) * invarea;
     ls->pdf = invarea;
 
-    if (type == LIGHT_SPOT) {
-      /* spot light attenuation */
-      float3 dir = make_float3(klight->spot.dir[0], klight->spot.dir[1], klight->spot.dir[2]);
-      ls->eval_fac *= spot_light_attenuation(
-          dir, klight->spot.spot_angle, klight->spot.spot_smooth, ls->Ng);
+    /* spot light attenuation */
+    float3 dir = make_float3(klight->spot.dir[0], klight->spot.dir[1], klight->spot.dir[2]);
+    ls->eval_fac *= spot_light_attenuation(
+        ls->Ng, klight->spot.spot_angle, klight->spot.spot_smooth, -ls->D);
 
-      if (ls->eval_fac == 0.0f) {
-        return false;
-      }
+    if (ls->eval_fac == 0.0f) {
+      return false;
     }
+
+    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(kg, ls->Ng, -ls->D, ls->t);
+  }
+  else if (type == LIGHT_POINT) {
+    float3 center = make_float3(klight->co[0], klight->co[1], klight->co[2]);
+
+    ls->Ng = normalize(ray_P - center);
+    float invarea = klight->spot.invarea;
+    ls->eval_fac = (0.25f * M_1_PI_F) * invarea;
+    ls->pdf = invarea;
+
     float2 uv = map_to_sphere(ls->Ng);
     ls->u = uv.x;
     ls->v = uv.y;
@@ -402,6 +465,8 @@ ccl_device bool light_sample_from_intersection(KernelGlobals kg,
     /* compute pdf */
     if (ls->t != FLT_MAX)
       ls->pdf *= lamp_light_pdf(kg, ls->Ng, -ls->D, ls->t);
+    else
+      ls->pdf = 0.f;
   }
   else if (type == LIGHT_AREA) {
     /* area light */
diff --git a/intern/cycles/util/math_intersect.h b/intern/cycles/util/math_intersect.h
index 54ce3ab4b66..0fce9ff24fd 100644
--- a/intern/cycles/util/math_intersect.h
+++ b/intern/cycles/util/math_intersect.h
@@ -85,6 +85,35 @@ ccl_device bool ray_aligned_disk_intersect(float3 ray_P,
   return true;
 }
 
+ccl_device bool ray_disk_intersect(float3 ray_P,
+                                   float3 ray_D,
+                                   float ray_t,
+                                   float3 disk_P,
+                                   float3 disk_N,
+                                   float disk_radius,
+                                   ccl_private float3 *isect_P,
+                                   ccl_private float *isect_t)
+{
+  const float3 vp = ray_P - disk_P;
+  const float dp = dot(vp, disk_N);
+  const float cos_angle = dot(disk_N, -ray_D);
+  if (dp * cos_angle > 0.f)  // front of light
+  {
+    float t = dp / cos_angle;
+    if (t < 0.f) { /* Ray points away from the light. */
+      return false;
+    }
+    float3 P = ray_P + t * ray_D;
+    float3 T = P - disk_P;
+    if (dot(T, T) < sqr(disk_radius) /*&& t > 0.f*/ && t <= ray_t) {
+      *isect_P = ray_P + t * ray_D;
+      *isect_t = t;
+      return true;
+    }
+  }
+  return false;
+}
+
 ccl_device_forceinline bool ray_triangle_intersect(float3 ray_P,
                                                    float3 ray_dir,
                                                    float ray_t,