Cycles: removed UV requirement for MNEE, along with fixes and cleanups

Remove need for shadow caustic caster geometry to have a UV layout. UVs were
useful to maintain a consistent tangent frame across the surface while
performing the walk. A consistent tangent frame is necessary for rough
surfaces where a normal offset encodes the sampled h, which should point
towards the same direction across the mesh.

In order to get a continuous surface parametrization without UVs, the
technique described in this paper was implemented:

"The Natural-Constraint Representation of the Path Space for Efficient
 Light Transport Simulation" (Supplementary Material), SIGGRAPH 2014.

In addition to implementing this feature:
* Shadow caustic casters without smooth normals are now ignored (triggered
  some refactoring and cleaning).
* Hit point calculation was refactored using existing utils functions,
  simplifying the code.
* The max number of solver iterations was reduced to 32, a solution is
  usually found by then.
* Added generalized geometry term clamping (transfer matrix calculation can
  sometimes get unstable).
* Add stop condition to Newton solver for more consistent CPU and GPU result.
* Add support for multi scatter GGX refraction.

Fixes T96990, T96991

Ref T94120

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

1 files changed, 115 insertions, 157 deletions

diff --git a/intern/cycles/kernel/integrator/mnee.h b/intern/cycles/kernel/integrator/mnee.h
index 7218165c67c..6eee597a0a1 100644
--- a/intern/cycles/kernel/integrator/mnee.h
+++ b/intern/cycles/kernel/integrator/mnee.h
@@ -36,11 +36,13 @@
  *  https://cg.ivd.kit.edu/english/HSLT.php
  */
 
-#  define MNEE_MAX_ITERATIONS 50
+#  define MNEE_MAX_ITERATIONS 32
 #  define MNEE_MAX_INTERSECTION_COUNT 10
 #  define MNEE_SOLVER_THRESHOLD 0.001f
+#  define MNEE_MINIMUM_STEP_SIZE 0.0001f
 #  define MNEE_MAX_CAUSTIC_CASTERS 6
 #  define MNEE_MIN_DISTANCE 0.001f
+#  define MNEE_MIN_PROGRESS_DISTANCE 0.0001f
 #  define MNEE_MIN_DETERMINANT 0.0001f
 #  define MNEE_PROJECTION_DISTANCE_MULTIPLIER 2.f
 
@@ -168,7 +170,8 @@ ccl_device_forceinline void mnee_setup_manifold_vertex(KernelGlobals kg,
                                                        const float2 n_offset,
                                                        ccl_private const Ray *ray,
                                                        ccl_private const Intersection *isect,
-                                                       ccl_private ShaderData *sd_vtx)
+                                                       ccl_private ShaderData *sd_vtx,
+                                                       bool seed)
 {
   sd_vtx->object = (isect->object == OBJECT_NONE) ? kernel_tex_fetch(__prim_object, isect->prim) :
                                                     isect->object;
@@ -177,7 +180,7 @@ ccl_device_forceinline void mnee_setup_manifold_vertex(KernelGlobals kg,
   sd_vtx->flag = 0;
   sd_vtx->object_flag = kernel_tex_fetch(__object_flag, sd_vtx->object);
 
-  /* matrices and time */
+  /* Matrices and time. */
   shader_setup_object_transforms(kg, sd_vtx, ray->time);
   sd_vtx->time = ray->time;
 
@@ -191,83 +194,29 @@ ccl_device_forceinline void mnee_setup_manifold_vertex(KernelGlobals kg,
 
   float3 verts[3];
   float3 normals[3];
-  uint4 tri_vindex = kernel_tex_fetch(__tri_vindex, sd_vtx->prim);
-
   if (sd_vtx->type & PRIMITIVE_TRIANGLE) {
-    /* Static triangle. */
-
-    /* Load triangle vertices. */
-    verts[0] = kernel_tex_fetch(__tri_verts, tri_vindex.w + 0);
-    verts[1] = kernel_tex_fetch(__tri_verts, tri_vindex.w + 1);
-    verts[2] = kernel_tex_fetch(__tri_verts, tri_vindex.w + 2);
-
-    /* Vectors. */
-    sd_vtx->P = triangle_point_from_uv(kg, sd_vtx, isect->object, isect->prim, isect->u, isect->v);
-
-    /* Smooth normal. */
-    if (sd_vtx->shader & SHADER_SMOOTH_NORMAL) {
-      /* Load triangle vertices. */
-      normals[0] = kernel_tex_fetch(__tri_vnormal, tri_vindex.x);
-      normals[1] = kernel_tex_fetch(__tri_vnormal, tri_vindex.y);
-      normals[2] = kernel_tex_fetch(__tri_vnormal, tri_vindex.z);
+    /* Load triangle vertices and normals. */
+    triangle_vertices_and_normals(kg, sd_vtx->prim, verts, normals);
+
+    /* Compute refined position (same code as in triangle_point_from_uv). */
+    sd_vtx->P = isect->u * verts[0] + isect->v * verts[1] + (1.f - isect->u - isect->v) * verts[2];
+    if (!(sd_vtx->object_flag & SD_OBJECT_TRANSFORM_APPLIED)) {
+      const Transform tfm = object_get_transform(kg, sd_vtx);
+      sd_vtx->P = transform_point(&tfm, sd_vtx->P);
     }
   }
   else { /* if (sd_vtx->type & PRIMITIVE_MOTION_TRIANGLE) */
-    /* Motion triangle. */
-
-    /* Get motion info. */
-    int numsteps, numverts;
-    object_motion_info(kg, sd_vtx->object, &numsteps, &numverts, NULL);
-
-    /* Figure out which steps we need to fetch and their interpolation factor. */
-    int maxstep = numsteps * 2;
-    int step = min((int)(sd_vtx->time * maxstep), maxstep - 1);
-    float t = sd_vtx->time * maxstep - step;
-
-    /* Find attribute. */
-    int offset = intersection_find_attribute(kg, sd_vtx->object, ATTR_STD_MOTION_VERTEX_POSITION);
-    kernel_assert(offset != ATTR_STD_NOT_FOUND);
-
-    /* Fetch vertex coordinates. */
-    float3 next_verts[3];
-    uint4 tri_vindex = kernel_tex_fetch(__tri_vindex, sd_vtx->prim);
-    motion_triangle_verts_for_step(kg, tri_vindex, offset, numverts, numsteps, step, verts);
-    motion_triangle_verts_for_step(
-        kg, tri_vindex, offset, numverts, numsteps, step + 1, next_verts);
-
-    /* Interpolate between steps. */
-    verts[0] = (1.0f - t) * verts[0] + t * next_verts[0];
-    verts[1] = (1.0f - t) * verts[1] + t * next_verts[1];
-    verts[2] = (1.0f - t) * verts[2] + t * next_verts[2];
+    /* Load triangle vertices and normals. */
+    motion_triangle_vertices_and_normals(
+        kg, sd_vtx->object, sd_vtx->prim, sd_vtx->time, verts, normals);
 
     /* Compute refined position. */
     sd_vtx->P = motion_triangle_point_from_uv(
         kg, sd_vtx, isect->object, isect->prim, isect->u, isect->v, verts);
-
-    /* Compute smooth normal. */
-    if (sd_vtx->shader & SHADER_SMOOTH_NORMAL) {
-      /* Find attribute. */
-      int offset = intersection_find_attribute(kg, sd_vtx->object, ATTR_STD_MOTION_VERTEX_NORMAL);
-      kernel_assert(offset != ATTR_STD_NOT_FOUND);
-
-      /* Fetch vertex coordinates. */
-      float3 next_normals[3];
-      motion_triangle_normals_for_step(kg, tri_vindex, offset, numverts, numsteps, step, normals);
-      motion_triangle_normals_for_step(
-          kg, tri_vindex, offset, numverts, numsteps, step + 1, next_normals);
-
-      /* Interpolate between steps. */
-      normals[0] = (1.0f - t) * normals[0] + t * next_normals[0];
-      normals[1] = (1.0f - t) * normals[1] + t * next_normals[1];
-      normals[2] = (1.0f - t) * normals[2] + t * next_normals[2];
-    }
   }
 
-  /* manifold vertex position */
-  vtx->p = sd_vtx->P;
-
+  /* Instance transform. */
   if (!(sd_vtx->object_flag & SD_OBJECT_TRANSFORM_APPLIED)) {
-    /* Instance transform. */
     object_position_transform_auto(kg, sd_vtx, &verts[0]);
     object_position_transform_auto(kg, sd_vtx, &verts[1]);
     object_position_transform_auto(kg, sd_vtx, &verts[2]);
@@ -277,94 +226,70 @@ ccl_device_forceinline void mnee_setup_manifold_vertex(KernelGlobals kg,
   }
 
   /* Tangent space (position derivatives) WRT barycentric (u, v). */
-  vtx->dp_du = verts[0] - verts[2];
-  vtx->dp_dv = verts[1] - verts[2];
+  float3 dp_du = verts[0] - verts[2];
+  float3 dp_dv = verts[1] - verts[2];
 
   /* Geometric normal. */
-  vtx->ng = normalize(cross(vtx->dp_du, vtx->dp_dv));
+  vtx->ng = normalize(cross(dp_du, dp_dv));
   if (sd_vtx->object_flag & SD_OBJECT_NEGATIVE_SCALE_APPLIED)
     vtx->ng = -vtx->ng;
 
-  /* Shading normal. */
-  if (!(sd_vtx->shader & SHADER_SMOOTH_NORMAL)) {
-    vtx->n = vtx->ng;
-    vtx->dn_du = vtx->dn_dv = zero_float3();
-  }
-  else {
-    /* Interpolate normals between vertices. */
-    float n_len;
-    vtx->n = normalize_len(normals[0] * sd_vtx->u + normals[1] * sd_vtx->v +
-                               normals[2] * (1.0f - sd_vtx->u - sd_vtx->v),
-                           &n_len);
-
-    /* Shading normal derivatives WRT barycentric (u, v)
-     * we calculate the derivative of n = |u*n0 + v*n1 + (1-u-v)*n2| using:
-     * d/du [f(u)/|f(u)|] = [d/du f(u)]/|f(u)| - f(u)/|f(u)|^3 <f(u), d/du f(u)>. */
-    const float inv_n_len = 1.f / n_len;
-    vtx->dn_du = inv_n_len * (normals[0] - normals[2]);
-    vtx->dn_dv = inv_n_len * (normals[1] - normals[2]);
-    vtx->dn_du -= vtx->n * dot(vtx->n, vtx->dn_du);
-    vtx->dn_dv -= vtx->n * dot(vtx->n, vtx->dn_dv);
-  }
+  /* Shading normals: Interpolate normals between vertices. */
+  float n_len;
+  vtx->n = normalize_len(normals[0] * sd_vtx->u + normals[1] * sd_vtx->v +
+                             normals[2] * (1.0f - sd_vtx->u - sd_vtx->v),
+                         &n_len);
+
+  /* Shading normal derivatives WRT barycentric (u, v)
+   * we calculate the derivative of n = |u*n0 + v*n1 + (1-u-v)*n2| using:
+   * d/du [f(u)/|f(u)|] = [d/du f(u)]/|f(u)| - f(u)/|f(u)|^3 <f(u), d/du f(u)>. */
+  const float inv_n_len = 1.f / n_len;
+  float3 dn_du = inv_n_len * (normals[0] - normals[2]);
+  float3 dn_dv = inv_n_len * (normals[1] - normals[2]);
+  dn_du -= vtx->n * dot(vtx->n, dn_du);
+  dn_dv -= vtx->n * dot(vtx->n, dn_dv);
 
-  /* dp_du and dp_dv need to be continuous across triangles for the h normal
-   * offset to yield a consistent halfvector while walking on the manifold.
-   * It's usually best to rely on the mesh uv layout, which is assumed to be
-   * continuous across the mesh. */
-  float2 duv0, duv1;
-  bool found_uv = false;
-  AttributeDescriptor uv_desc = find_attribute(kg, sd_vtx, ATTR_STD_GENERATED);
-  if (uv_desc.offset != ATTR_STD_NOT_FOUND) {
-    float3 uvs[3];
-    uvs[0] = kernel_tex_fetch(__attributes_float3, uv_desc.offset + tri_vindex.x);
-    uvs[1] = kernel_tex_fetch(__attributes_float3, uv_desc.offset + tri_vindex.y);
-    uvs[2] = kernel_tex_fetch(__attributes_float3, uv_desc.offset + tri_vindex.z);
-    duv0 = make_float2(uvs[0].x - uvs[2].x, uvs[0].y - uvs[2].y);
-    duv1 = make_float2(uvs[1].x - uvs[2].x, uvs[1].y - uvs[2].y);
-    found_uv = true;
+  /* Orthonormalize (dp_du,dp_dv) using a linear transformation, which
+   * we use on (dn_du,dn_dv) as well so the new (u,v) are consistent. */
+  const float inv_len_dp_du = 1.f / len(dp_du);
+  dp_du *= inv_len_dp_du;
+  dn_du *= inv_len_dp_du;
+
+  const float dpdu_dot_dpdv = dot(dp_du, dp_dv);
+  dp_dv -= dpdu_dot_dpdv * dp_du;
+  dn_dv -= dpdu_dot_dpdv * dn_du;
+
+  const float inv_len_dp_dv = 1.f / len(dp_dv);
+  dp_dv *= inv_len_dp_dv;
+  dn_dv *= inv_len_dp_dv;
+
+  /* Final local differential geometry. */
+  if (seed) {
+    vtx->dp_du = dp_du;
+    vtx->dp_dv = dp_dv;
+    vtx->dn_du = dn_du;
+    vtx->dn_dv = dn_dv;
   }
   else {
-    uv_desc = find_attribute(kg, sd_vtx, ATTR_STD_UV);
-    if (uv_desc.offset != ATTR_STD_NOT_FOUND) {
-      float2 uvs[3];
-      uvs[0] = kernel_tex_fetch(__attributes_float2, uv_desc.offset + tri_vindex.x);
-      uvs[1] = kernel_tex_fetch(__attributes_float2, uv_desc.offset + tri_vindex.y);
-      uvs[2] = kernel_tex_fetch(__attributes_float2, uv_desc.offset + tri_vindex.z);
-      duv0 = make_float2(uvs[0].x - uvs[2].x, uvs[0].y - uvs[2].y);
-      duv1 = make_float2(uvs[1].x - uvs[2].x, uvs[1].y - uvs[2].y);
-      found_uv = true;
-    }
-  }
-  if (found_uv) {
-    const float det = duv0.x * duv1.y - duv0.y * duv1.x;
-    if (det != 0.f) {
-      const float inv_det = 1.f / det;
-
-      /* Tangent space (position derivatives) WRT texture (u, v). */
-      const float3 dp_du = vtx->dp_du;
-      const float3 dp_dv = vtx->dp_dv;
-      vtx->dp_du = (duv1.y * dp_du - duv0.y * dp_dv) * inv_det;
-      vtx->dp_dv = (-duv1.x * dp_du + duv0.x * dp_dv) * inv_det;
-
-      /* Shading normal derivatives WRT texture (u, v). */
-      const float3 dn_du = vtx->dn_du;
-      const float3 dn_dv = vtx->dn_dv;
-      vtx->dn_du = (duv1.y * dn_du - duv0.y * dn_dv) * inv_det;
-      vtx->dn_dv = (-duv1.x * dn_du + duv0.x * dn_dv) * inv_det;
-    }
+    /* Find angle subtended by reference direction (travel direction). */
+    const float3 reference_direction = normalize(sd_vtx->P - vtx->p);
+    const float reference_theta = atan2(dot(reference_direction, vtx->dp_dv),
+                                        dot(reference_direction, vtx->dp_du));
+    const float current_theta = atan2(dot(reference_direction, dp_dv),
+                                      dot(reference_direction, dp_du));
+    const float theta = reference_theta - current_theta;
+
+    /* Rotate (dp_du,dp_dv) to be consistent with previous tangent frame. */
+    float cos_theta, sin_theta;
+    fast_sincosf(theta, &sin_theta, &cos_theta);
+    vtx->dp_du = cos_theta * dp_du - sin_theta * dp_dv;
+    vtx->dp_dv = sin_theta * dp_du + cos_theta * dp_dv;
+    vtx->dn_du = cos_theta * dn_du - sin_theta * dn_dv;
+    vtx->dn_dv = sin_theta * dn_du + cos_theta * dn_dv;
   }
 
-  /* Orthonormalize (dp_du,dp_dv) using a linear transformation, which
-   * we use on (dn_du,dn_dv) as well so the new (u,v) are consistent. */
-  const float inv_len_dp_du = 1.f / len(vtx->dp_du);
-  vtx->dp_du *= inv_len_dp_du;
-  vtx->dn_du *= inv_len_dp_du;
-  const float dpdu_dot_dpdv = dot(vtx->dp_du, vtx->dp_dv);
-  const float3 dp_dv = vtx->dp_dv - dpdu_dot_dpdv * vtx->dp_du;
-  const float3 dn_dv = vtx->dn_dv - dpdu_dot_dpdv * vtx->dn_du;
-  const float inv_len_dp_dv = 1.f / len(dp_dv);
-  vtx->dp_dv = dp_dv * inv_len_dp_dv;
-  vtx->dn_dv = dn_dv * inv_len_dp_dv;
+  /* Manifold vertex position. */
+  vtx->p = sd_vtx->P;
 
   /* Initialize constraint and its derivates. */
   vtx->a = vtx->c = zero_float4();
@@ -638,15 +563,28 @@ ccl_device_forceinline bool mnee_newton_solver(KernelGlobals kg,
         break;
       }
 
+      /* Initialize tangent frame, which will be used as reference. */
+      ccl_private ManifoldVertex &tv = tentative[vi];
+      tv.p = mv.p;
+      tv.dp_du = mv.dp_du;
+      tv.dp_dv = mv.dp_dv;
+
       /* Setup corrected manifold vertex. */
       mnee_setup_manifold_vertex(kg,
-                                 &tentative[vi],
+                                 &tv,
                                  mv.bsdf,
                                  mv.eta,
                                  mv.n_offset,
                                  &projection_ray,
                                  &projection_isect,
-                                 sd_vtx);
+                                 sd_vtx,
+                                 false);
+
+      /* Fail newton solve if we are not making progress, probably stuck trying to move off the
+       * edge of the mesh. */
+      const float distance = len(tv.p - mv.p);
+      if (distance < MNEE_MIN_PROGRESS_DISTANCE)
+        return false;
     }
 
     /* Check that tentative path is still transmissive. */
@@ -672,6 +610,11 @@ ccl_device_forceinline bool mnee_newton_solver(KernelGlobals kg,
       reduce_stepsize = false;
       resolve_constraint = false;
       beta *= .5f;
+
+      /* Fail newton solve if the stepsize is too small. */
+      if (beta < MNEE_MINIMUM_STEP_SIZE)
+        return false;
+
       continue;
     }
 
@@ -714,8 +657,7 @@ mnee_sample_bsdf_dh(ClosureType type, float alpha_x, float alpha_y, float sample
   if (type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID) {
     tan2_theta *= -logf(1.0f - sample_u);
   }
-  else { /* if (type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID || type ==
-            CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_FRESNEL_ID) */
+  else { /* type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID assumed */
     tan2_theta *= sample_u / (1.0f - sample_u);
   }
   float cos2_theta = 1.0f / (1.0f + tan2_theta);
@@ -749,8 +691,7 @@ ccl_device_forceinline float3 mnee_eval_bsdf_contribution(ccl_private ShaderClos
     /* Eq. 26, 27: now calculate G1(i,m) and G1(o,m). */
     G = bsdf_beckmann_G1(bsdf->alpha_x, cosNO) * bsdf_beckmann_G1(bsdf->alpha_x, cosNI);
   }
-  else { /* if (type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID || type ==
-             CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_FRESNEL_ID) */
+  else { /* bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID assumed */
     /* Eq. 34: now calculate G1(i,m) and G1(o,m). */
     G = (2.f / (1.f + safe_sqrtf(1.f + alpha2 * (1.f - cosNO * cosNO) / (cosNO * cosNO)))) *
         (2.f / (1.f + safe_sqrtf(1.f + alpha2 * (1.f - cosNI * cosNI) / (cosNI * cosNI))));
@@ -929,7 +870,8 @@ ccl_device_forceinline bool mnee_path_contribution(KernelGlobals kg,
   /* Receiver bsdf eval above already contains |n.wo|. */
   const float dw0_dx1 = fabsf(dot(wo, vertices[0].n)) / sqr(wo_len);
 
-  const float G = dw0_dx1 * dx1_dxlight;
+  /* Clamp since it has a tendency to be unstable. */
+  const float G = fminf(dw0_dx1 * dx1_dxlight, 2.f);
   bsdf_eval_mul(throughput, G);
 
   /* Specular reflectance. */
@@ -1058,21 +1000,36 @@ ccl_device_forceinline bool kernel_path_mnee_sample(KernelGlobals kg,
       if (vertex_count >= MNEE_MAX_CAUSTIC_CASTERS)
         return false;
 
+      /* Reject caster if it is not a triangles mesh. */
+      if (!(probe_isect.type & PRIMITIVE_TRIANGLE))
+        return false;
+
       ccl_private ManifoldVertex &mv = vertices[vertex_count++];
 
       /* Setup shader data on caustic caster and evaluate context. */
       shader_setup_from_ray(kg, sd_mnee, &probe_ray, &probe_isect);
 
+      /* Reject caster if smooth normals are not available: Manifold exploration assumes local
+       * differential geometry can be created at any point on the surface which is not possible if
+       * normals are not smooth. */
+      if (!(sd_mnee->shader & SHADER_SMOOTH_NORMAL))
+        return false;
+
       /* Last bool argument is the MNEE flag (for TINY_MAX_CLOSURE cap in kernel_shader.h). */
       shader_eval_surface<KERNEL_FEATURE_NODE_MASK_SURFACE_SHADOW>(
           kg, state, sd_mnee, NULL, PATH_RAY_DIFFUSE, true);
 
-      /* get and sample refraction bsdf. */
+      /* Get and sample refraction bsdf */
       bool found_transimissive_microfacet_bsdf = false;
       for (int ci = 0; ci < sd_mnee->num_closure; ci++) {
         ccl_private ShaderClosure *bsdf = &sd_mnee->closure[ci];
         if (bsdf->type == CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID ||
-            bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID) {
+            bsdf->type == CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID ||
+            bsdf->type == CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID ||
+            bsdf->type == CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_FRESNEL_ID) {
+          /* Note that CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID and
+             CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_FRESNEL_ID are treated as
+             CLOSURE_BSDF_MICROFACET_GGX_REFRACTION_ID further below */
 
           found_transimissive_microfacet_bsdf = true;
           ccl_private MicrofacetBsdf *microfacet_bsdf = (ccl_private MicrofacetBsdf *)bsdf;
@@ -1088,7 +1045,8 @@ ccl_device_forceinline bool kernel_path_mnee_sample(KernelGlobals kg,
               bsdf->type, microfacet_bsdf->alpha_x, microfacet_bsdf->alpha_y, bsdf_u, bsdf_v);
 
           /* Setup differential geometry on vertex. */
-          mnee_setup_manifold_vertex(kg, &mv, bsdf, eta, h, &probe_ray, &probe_isect, sd_mnee);
+          mnee_setup_manifold_vertex(
+              kg, &mv, bsdf, eta, h, &probe_ray, &probe_isect, sd_mnee, true);
           break;
         }
       }