Cycles: approximate shadow caustics using manifold next event estimation

This adds support for selective rendering of caustics in shadows of refractive
objects. Example uses are rendering of underwater caustics and eye caustics.

This is based on "Manifold Next Event Estimation", a method developed for
production rendering. The idea is to selectively enable shadow caustics on a
few objects in the scene where they have a big visual impact, without impacting
render performance for the rest of the scene.

The Shadow Caustic option must be manually enabled on light, caustic receiver
and caster objects. For such light paths, the Filter Glossy option will be
ignored and replaced by sharp caustics.

Currently this method has a various limitations:

* Only caustics in shadows of refractive objects work, which means no caustics
  from reflection or caustics that outside shadows. Only up to 4 refractive
  caustic bounces are supported.
* Caustic caster objects should have smooth normals.
* Not currently support for Metal GPU rendering.

In the future this method may be extended for more general caustics.

TECHNICAL DETAILS

This code adds manifold next event estimation through refractive surface(s) as a
new sampling technique for direct lighting, i.e. finding the point on the
refractive surface(s) along the path to a light sample, which satisfies Fermat's
principle for a given microfacet normal and the path's end points. This
technique involves walking on the "specular manifold" using a pseudo newton
solver. Such a manifold is defined by the specular constraint matrix from the
manifold exploration framework [2]. For each refractive interface, this
constraint is defined by enforcing that the generalized half-vector projection
onto the interface local tangent plane is null. The newton solver guides the
walk by linearizing the manifold locally before reprojecting the linear solution
onto the refractive surface. See paper [1] for more details about the technique
itself and [3] for the half-vector light transport formulation, from which it is
derived.

[1] Manifold Next Event Estimation
Johannes Hanika, Marc Droske, and Luca Fascione. 2015.
Comput. Graph. Forum 34, 4 (July 2015), 87–97.
https://jo.dreggn.org/home/2015_mnee.pdf

[2] Manifold exploration: a Markov Chain Monte Carlo technique for rendering
scenes with difficult specular transport Wenzel Jakob and Steve Marschner.
2012. ACM Trans. Graph. 31, 4, Article 58 (July 2012), 13 pages.
https://www.cs.cornell.edu/projects/manifolds-sg12/

[3] The Natural-Constraint Representation of the Path Space for Efficient
Light Transport Simulation. Anton S. Kaplanyan, Johannes Hanika, and Carsten
Dachsbacher. 2014. ACM Trans. Graph. 33, 4, Article 102 (July 2014), 13 pages.
https://cg.ivd.kit.edu/english/HSLT.php

The code for this samping technique was inserted at the light sampling stage
(direct lighting). If the walk is successful, it turns off path regularization
using a specialized flag in the path state (PATH_MNEE_SUCCESS). This flag tells
the integrator not to blur the brdf roughness further down the path (in a child
ray created from BSDF sampling). In addition, using a cascading mechanism of
flag values, we cull connections to caustic lights for this and children rays,
which should be resolved through MNEE.

This mechanism also cancels the MIS bsdf counter part at the casutic receiver
depth, in essence leaving MNEE as the only sampling technique from receivers
through refractive casters to caustic lights. This choice might not be optimal
when the light gets large wrt to the receiver, though this is usually not when
you want to use MNEE.

This connection culling strategy removes a fair amount of fireflies, at the cost
of introducing a slight bias. Because of the selective nature of the culling
mechanism, reflective caustics still benefit from the native path
regularization, which further removes fireflies on other surfaces (bouncing
light off casters).

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

intern/cycles/blender/addon/properties.py

@@ -1012,6 +1012,12 @@ class CyclesLightSettings(bpy.types.PropertyGroup):
         "note that this will make the light invisible",
         default=False,
     )
+    is_caustics_light: BoolProperty(
+        name="Shadow Caustics",
+        description="Generate approximate caustics in shadows of refractive surfaces. "
+        "Lights, caster and receiver objects must have shadow caustics options set to enable this",
+        default=False,
+    )
 
     @classmethod
     def register(cls):
@@ -1028,6 +1034,12 @@ class CyclesLightSettings(bpy.types.PropertyGroup):
 
 class CyclesWorldSettings(bpy.types.PropertyGroup):
 
+    is_caustics_light: BoolProperty(
+        name="Shadow Caustics",
+        description="Generate approximate caustics in shadows of refractive surfaces. "
+        "Lights, caster and receiver objects must have shadow caustics options set to enable this",
+        default=False,
+    )
     sampling_method: EnumProperty(
         name="Sampling Method",
         description="How to sample the background light",
@@ -1226,6 +1238,21 @@ class CyclesObjectSettings(bpy.types.PropertyGroup):
         subtype='DISTANCE',
     )
 
+    is_caustics_caster: BoolProperty(
+        name="Cast Shadow Caustics",
+        description="With refractive materials, generate approximate caustics in shadows of this object. "
+        "Up to 10 bounces inside this object are taken into account. Lights, caster and receiver objects "
+        "must have shadow caustics options set to enable this",
+        default=False,
+    )
+
+    is_caustics_receiver: BoolProperty(
+        name="Receive Shadow Caustics",
+        description="Receive approximate caustics from refractive materials in shadows on this object. "
+        "Lights, caster and receiver objects must have shadow caustics options set to enable this",
+        default=False,
+    )
+
     @classmethod
     def register(cls):
         bpy.types.Object.cycles = PointerProperty(

intern/cycles/blender/addon/ui.py

@@ -1082,7 +1082,7 @@ class CYCLES_OBJECT_PT_shading(CyclesButtonsPanel, Panel):
             return False
 
         ob = context.object
-        return ob and has_geometry_visibility(ob)
+        return ob and has_geometry_visibility(ob) and ob.type != 'LIGHT'
 
     def draw(self, context):
         pass
@@ -1125,6 +1125,28 @@ class CYCLES_OBJECT_PT_shading_gi_approximation(CyclesButtonsPanel, Panel):
         col.prop(cob, "ao_distance")
 
 
+class CYCLES_OBJECT_PT_shading_caustics(CyclesButtonsPanel, Panel):
+    bl_label = "Caustics"
+    bl_parent_id = "CYCLES_OBJECT_PT_shading"
+    bl_context = "object"
+
+    @classmethod
+    def poll(cls, context):
+        return CyclesButtonsPanel.poll(context) and not use_metal(context)
+
+    def draw(self, context):
+        layout = self.layout
+        layout.use_property_split = True
+        layout.use_property_decorate = False
+
+        col = layout.column()
+
+        ob = context.object
+        cob = ob.cycles
+        col.prop(cob, "is_caustics_caster")
+        col.prop(cob, "is_caustics_receiver")
+
+
 class CYCLES_OBJECT_PT_visibility(CyclesButtonsPanel, Panel):
     bl_label = "Visibility"
     bl_context = "object"
@@ -1300,6 +1322,8 @@ class CYCLES_LIGHT_PT_light(CyclesButtonsPanel, Panel):
         sub.active = not (light.type == 'AREA' and clamp.is_portal)
         sub.prop(clamp, "cast_shadow")
         sub.prop(clamp, "use_multiple_importance_sampling", text="Multiple Importance")
+        if not use_metal(context):
+            sub.prop(clamp, "is_caustics_light", text="Shadow Caustics")
 
         if light.type == 'AREA':
             col.prop(clamp, "is_portal", text="Portal")
@@ -1496,6 +1520,8 @@ class CYCLES_WORLD_PT_settings_surface(CyclesButtonsPanel, Panel):
         subsub.active = cworld.sampling_method == 'MANUAL'
         subsub.prop(cworld, "sample_map_resolution")
         sub.prop(cworld, "max_bounces")
+        sub.prop(cworld, "is_caustics_light", text="Shadow Caustics")
+
 
 
 class CYCLES_WORLD_PT_settings_volume(CyclesButtonsPanel, Panel):
@@ -2193,6 +2219,7 @@ classes = (
     CYCLES_OBJECT_PT_shading,
     CYCLES_OBJECT_PT_shading_shadow_terminator,
     CYCLES_OBJECT_PT_shading_gi_approximation,
+    CYCLES_OBJECT_PT_shading_caustics,
     CYCLES_OBJECT_PT_visibility,
     CYCLES_OBJECT_PT_visibility_ray_visibility,
     CYCLES_OBJECT_PT_visibility_culling,

intern/cycles/blender/light.cpp

@@ -114,6 +114,9 @@ void BlenderSync::sync_light(BL::Object &b_parent,
   light->set_cast_shadow(get_boolean(clight, "cast_shadow"));
   light->set_use_mis(get_boolean(clight, "use_multiple_importance_sampling"));
 
+  /* caustics light */
+  light->set_use_caustics(get_boolean(clight, "is_caustics_light"));
+
   light->set_max_bounces(get_int(clight, "max_bounces"));
 
   if (b_ob_info.real_object != b_ob_info.iter_object) {
@@ -176,6 +179,9 @@ void BlenderSync::sync_background_light(BL::SpaceView3D &b_v3d, bool use_portal)
         /* force enable light again when world is resynced */
         light->set_is_enabled(true);
 
+        /* caustic light */
+        light->set_use_caustics(get_boolean(cworld, "is_caustics_light"));
+
         light->tag_update(scene);
         light_map.set_recalc(b_world);
       }

intern/cycles/blender/object.cpp

@@ -298,6 +298,12 @@ Object *BlenderSync::sync_object(BL::Depsgraph &b_depsgraph,
   }
   object->set_ao_distance(ao_distance);
 
+  bool is_caustics_caster = get_boolean(cobject, "is_caustics_caster");
+  object->set_is_caustics_caster(is_caustics_caster);
+
+  bool is_caustics_receiver = get_boolean(cobject, "is_caustics_receiver");
+  object->set_is_caustics_receiver(is_caustics_receiver);
+
   /* sync the asset name for Cryptomatte */
   BL::Object parent = b_ob.parent();
   ustring parent_name;

intern/cycles/kernel/CMakeLists.txt

@@ -223,6 +223,7 @@ set(SRC_KERNEL_INTEGRATOR_HEADERS
   integrator/intersect_subsurface.h
   integrator/intersect_volume_stack.h
   integrator/megakernel.h
+  integrator/mnee.h
   integrator/path_state.h
   integrator/shade_background.h
   integrator/shade_light.h

intern/cycles/kernel/integrator/init_from_bake.h

@@ -230,7 +230,11 @@ ccl_device bool integrator_init_from_bake(KernelGlobals kg,
     /* Setup next kernel to execute. */
     const int shader_index = shader & SHADER_MASK;
     const int shader_flags = kernel_tex_fetch(__shaders, shader_index).flags;
-    if (shader_flags & SD_HAS_RAYTRACE) {
+    const bool use_caustics = kernel_data.integrator.use_caustics &&
+                              (object_flag & SD_OBJECT_CAUSTICS);
+    const bool use_raytrace_kernel = (shader_flags & SD_HAS_RAYTRACE) || use_caustics;
+
+    if (use_raytrace_kernel) {
       INTEGRATOR_PATH_INIT_SORTED(DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE, shader_index);
     }
     else {

intern/cycles/kernel/integrator/intersect_closest.h

@@ -123,7 +123,9 @@ ccl_device_forceinline void integrator_split_shadow_catcher(
   /* Continue with shading shadow catcher surface. */
   const int shader = intersection_get_shader(kg, isect);
   const int flags = kernel_tex_fetch(__shaders, shader).flags;
-  const bool use_raytrace_kernel = (flags & SD_HAS_RAYTRACE);
+  const bool use_caustics = kernel_data.integrator.use_caustics &&
+                            (object_flags & SD_OBJECT_CAUSTICS);
+  const bool use_raytrace_kernel = (flags & SD_HAS_RAYTRACE) || use_caustics;
 
   if (use_raytrace_kernel) {
     INTEGRATOR_PATH_INIT_SORTED(DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE, shader);
@@ -145,7 +147,10 @@ ccl_device_forceinline void integrator_intersect_next_kernel_after_shadow_catche
 
   const int shader = intersection_get_shader(kg, &isect);
   const int flags = kernel_tex_fetch(__shaders, shader).flags;
-  const bool use_raytrace_kernel = (flags & SD_HAS_RAYTRACE);
+  const int object_flags = intersection_get_object_flags(kg, &isect);
+  const bool use_caustics = kernel_data.integrator.use_caustics &&
+                            (object_flags & SD_OBJECT_CAUSTICS);
+  const bool use_raytrace_kernel = (flags & SD_HAS_RAYTRACE) || use_caustics;
 
   if (use_raytrace_kernel) {
     INTEGRATOR_PATH_NEXT_SORTED(
@@ -214,7 +219,10 @@ ccl_device_forceinline void integrator_intersect_next_kernel(
       const int flags = kernel_tex_fetch(__shaders, shader).flags;
 
       if (!integrator_intersect_terminate(kg, state, flags)) {
-        const bool use_raytrace_kernel = (flags & SD_HAS_RAYTRACE);
+        const int object_flags = intersection_get_object_flags(kg, isect);
+        const bool use_caustics = kernel_data.integrator.use_caustics &&
+                                  (object_flags & SD_OBJECT_CAUSTICS);
+        const bool use_raytrace_kernel = (flags & SD_HAS_RAYTRACE) || use_caustics;
         if (use_raytrace_kernel) {
           INTEGRATOR_PATH_NEXT_SORTED(
               current_kernel, DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE, shader);
@@ -261,7 +269,10 @@ ccl_device_forceinline void integrator_intersect_next_kernel_after_volume(
       /* Hit a surface, continue with surface kernel unless terminated. */
       const int shader = intersection_get_shader(kg, isect);
       const int flags = kernel_tex_fetch(__shaders, shader).flags;
-      const bool use_raytrace_kernel = (flags & SD_HAS_RAYTRACE);
+      const int object_flags = intersection_get_object_flags(kg, isect);
+      const bool use_caustics = kernel_data.integrator.use_caustics &&
+                                (object_flags & SD_OBJECT_CAUSTICS);
+      const bool use_raytrace_kernel = (flags & SD_HAS_RAYTRACE) || use_caustics;
 
       if (use_raytrace_kernel) {
         INTEGRATOR_PATH_NEXT_SORTED(
@@ -328,12 +339,37 @@ ccl_device void integrator_intersect_closest(KernelGlobals kg,
     isect.prim = PRIM_NONE;
   }
 
+  /* Setup mnee flag to signal last intersection with a caster */
+  const uint32_t path_flag = INTEGRATOR_STATE(state, path, flag);
+
+#ifdef __MNEE__
+  /* Path culling logic for MNEE (removes fireflies at the cost of bias) */
+  if (kernel_data.integrator.use_caustics) {
+    /* The following firefly removal mechanism works by culling light connections when
+     * a ray comes from a caustic caster directly after bouncing off a different caustic
+     * receiver */
+    bool from_caustic_caster = false;
+    bool from_caustic_receiver = false;
+    if (!(path_flag & PATH_RAY_CAMERA) && last_isect_object != OBJECT_NONE) {
+      const int object_flags = kernel_tex_fetch(__object_flag, last_isect_object);
+      from_caustic_receiver = (object_flags & SD_OBJECT_CAUSTICS_RECEIVER);
+      from_caustic_caster = (object_flags & SD_OBJECT_CAUSTICS_CASTER);
+    }
+
+    bool has_receiver_ancestor = INTEGRATOR_STATE(state, path, mnee) & PATH_MNEE_RECEIVER_ANCESTOR;
+    INTEGRATOR_STATE_WRITE(state, path, mnee) &= ~PATH_MNEE_CULL_LIGHT_CONNECTION;
+    if (from_caustic_caster && has_receiver_ancestor)
+      INTEGRATOR_STATE_WRITE(state, path, mnee) |= PATH_MNEE_CULL_LIGHT_CONNECTION;
+    if (from_caustic_receiver)
+      INTEGRATOR_STATE_WRITE(state, path, mnee) |= PATH_MNEE_RECEIVER_ANCESTOR;
+  }
+#endif /* __MNEE__ */
+
   /* Light intersection for MIS. */
   if (kernel_data.integrator.use_lamp_mis) {
     /* NOTE: if we make lights visible to camera rays, we'll need to initialize
      * 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, state, &ray, &isect, last_isect_prim, last_isect_object, last_type, path_flag) ||
           hit;

intern/cycles/kernel/integrator/path_state.h

@@ -57,6 +57,10 @@ ccl_device_inline void path_state_init_integrator(KernelGlobals kg,
   INTEGRATOR_STATE_WRITE(state, path, continuation_probability) = 1.0f;
   INTEGRATOR_STATE_WRITE(state, path, throughput) = make_float3(1.0f, 1.0f, 1.0f);
 
+#ifdef __MNEE__
+  INTEGRATOR_STATE_WRITE(state, path, mnee) = 0;
+#endif
+
   INTEGRATOR_STATE_WRITE(state, isect, object) = OBJECT_NONE;
   INTEGRATOR_STATE_WRITE(state, isect, prim) = PRIM_NONE;
 

intern/cycles/kernel/integrator/shade_background.h

@@ -101,6 +101,22 @@ ccl_device_inline void integrate_background(KernelGlobals kg,
 #endif
   }
 
+#ifdef __MNEE__
+  if (INTEGRATOR_STATE(state, path, mnee) & PATH_MNEE_CULL_LIGHT_CONNECTION) {
+    if (kernel_data.background.use_mis) {
+      for (int lamp = 0; lamp < kernel_data.integrator.num_all_lights; lamp++) {
+        /* This path should have been resolved with mnee, it will
+         * generate a firefly for small lights since it is improbable. */
+        const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, lamp);
+        if (klight->type == LIGHT_BACKGROUND && klight->use_caustics) {
+          eval_background = false;
+          break;
+        }
+      }
+    }
+  }
+#endif /* __MNEE__ */
+
   /* Evaluate background shader. */
   float3 L = (eval_background) ? integrator_eval_background_shader(kg, state, render_buffer) :
                                  zero_float3();
@@ -140,6 +156,16 @@ ccl_device_inline void integrate_distant_lights(KernelGlobals kg,
       }
 #endif
 
+#ifdef __MNEE__
+      if (INTEGRATOR_STATE(state, path, mnee) & PATH_MNEE_CULL_LIGHT_CONNECTION) {
+        /* This path should have been resolved with mnee, it will
+         * generate a firefly for small lights since it is improbable. */
+        const ccl_global KernelLight *klight = &kernel_tex_fetch(__lights, lamp);
+        if (klight->use_caustics)
+          return;
+      }
+#endif /* __MNEE__ */
+
       /* Evaluate light shader. */
       /* TODO: does aliasing like this break automatic SoA in CUDA? */
       ShaderDataTinyStorage emission_sd_storage;

intern/cycles/kernel/integrator/shade_surface.h

@@ -6,6 +6,8 @@
 #include "kernel/film/accumulate.h"
 #include "kernel/film/passes.h"
 
+#include "kernel/integrator/mnee.h"
+
 #include "kernel/integrator/path_state.h"
 #include "kernel/integrator/shader_eval.h"
 #include "kernel/integrator/subsurface.h"
@@ -92,6 +94,7 @@ ccl_device_forceinline void integrate_surface_emission(KernelGlobals kg,
 #ifdef __EMISSION__
 /* Path tracing: sample point on light and evaluate light shader, then
  * queue shadow ray to be traced. */
+template<uint node_feature_mask>
 ccl_device_forceinline void integrate_surface_direct_light(KernelGlobals kg,
                                                            IntegratorState state,
                                                            ccl_private ShaderData *sd,
@@ -124,34 +127,65 @@ ccl_device_forceinline void integrate_surface_direct_light(KernelGlobals kg,
    * integrate_surface_bounce, evaluate the BSDF, and only then evaluate
    * the light shader. This could also move to its own kernel, for
    * non-constant light sources. */
-  ShaderDataTinyStorage emission_sd_storage;
+  ShaderDataCausticsStorage emission_sd_storage;
   ccl_private ShaderData *emission_sd = AS_SHADER_DATA(&emission_sd_storage);
-  const float3 light_eval = light_sample_shader_eval(kg, state, emission_sd, &ls, sd->time);
-  if (is_zero(light_eval)) {
-    return;
-  }
 
-  /* Evaluate BSDF. */
+  Ray ray ccl_optional_struct_init;
+  BsdfEval bsdf_eval ccl_optional_struct_init;
   const bool is_transmission = shader_bsdf_is_transmission(sd, ls.D);
 
-  BsdfEval bsdf_eval ccl_optional_struct_init;
-  const float bsdf_pdf = shader_bsdf_eval(kg, sd, ls.D, is_transmission, &bsdf_eval, ls.shader);
-  bsdf_eval_mul3(&bsdf_eval, light_eval / ls.pdf);
+#  ifdef __MNEE__
+  bool skip_nee = false;
+  IF_KERNEL_NODES_FEATURE(RAYTRACE)
+  {
+    if (ls.lamp != LAMP_NONE) {
+      /* Is this a caustic light? */
+      const bool use_caustics = kernel_tex_fetch(__lights, ls.lamp).use_caustics;
+      if (use_caustics) {
+        /* Are we on a caustic caster? */
+        if (is_transmission && (sd->object_flag & SD_OBJECT_CAUSTICS_CASTER))
+          return;
 
-  if (ls.shader & SHADER_USE_MIS) {
-    const float mis_weight = light_sample_mis_weight_nee(kg, ls.pdf, bsdf_pdf);
-    bsdf_eval_mul(&bsdf_eval, mis_weight);
+        /* Are we on a caustic receiver? */
+        if (!is_transmission && (sd->object_flag & SD_OBJECT_CAUSTICS_RECEIVER))
+          skip_nee = kernel_path_mnee_sample(
+              kg, state, sd, emission_sd, rng_state, &ls, &bsdf_eval);
+      }
+    }
+  }
+  if (skip_nee) {
+    /* Create shadow ray after successful manifold walk:
+     * emission_sd contains the last interface intersection and
+     * the light sample ls has been updated */
+    light_sample_to_surface_shadow_ray(kg, emission_sd, &ls, &ray);
+  }
+  else
+#  endif /* __MNEE__ */
+  {
+    const float3 light_eval = light_sample_shader_eval(kg, state, emission_sd, &ls, sd->time);
+    if (is_zero(light_eval)) {
+      return;
+    }
+
+    /* Evaluate BSDF. */
+    const float bsdf_pdf = shader_bsdf_eval(kg, sd, ls.D, is_transmission, &bsdf_eval, ls.shader);
+    bsdf_eval_mul3(&bsdf_eval, light_eval / ls.pdf);
+
+    if (ls.shader & SHADER_USE_MIS) {
+      const float mis_weight = light_sample_mis_weight_nee(kg, ls.pdf, bsdf_pdf);
+      bsdf_eval_mul(&bsdf_eval, mis_weight);
+    }
+
+    /* Path termination. */
+    const float terminate = path_state_rng_light_termination(kg, rng_state);
+    if (light_sample_terminate(kg, &ls, &bsdf_eval, terminate)) {
+      return;
+    }
+
+    /* Create shadow ray. */
+    light_sample_to_surface_shadow_ray(kg, sd, &ls, &ray);
   }
 
-  /* Path termination. */
-  const float terminate = path_state_rng_light_termination(kg, rng_state);
-  if (light_sample_terminate(kg, &ls, &bsdf_eval, terminate)) {
-    return;
-  }
-
-  /* Create shadow ray. */
-  Ray ray ccl_optional_struct_init;
-  light_sample_to_surface_shadow_ray(kg, sd, &ls, &ray);
   const bool is_light = light_sample_is_light(&ls);
 
   /* Branch off shadow kernel. */
@@ -501,7 +535,7 @@ ccl_device bool integrate_surface(KernelGlobals kg,
 
     /* Direct light. */
     PROFILING_EVENT(PROFILING_SHADE_SURFACE_DIRECT_LIGHT);
-    integrate_surface_direct_light(kg, state, &sd, &rng_state);
+    integrate_surface_direct_light<node_feature_mask>(kg, state, &sd, &rng_state);
 
 #if defined(__AO__)
     /* Ambient occlusion pass. */

intern/cycles/kernel/integrator/shader_eval.h

@@ -157,7 +157,11 @@ ccl_device_inline void shader_prepare_surface_closures(KernelGlobals kg,
    *
    * Blurring of bsdf after bounces, for rays that have a small likelihood
    * of following this particular path (diffuse, rough glossy) */
-  if (kernel_data.integrator.filter_glossy != FLT_MAX) {
+  if (kernel_data.integrator.filter_glossy != FLT_MAX
+#ifdef __MNEE__
+      && !(INTEGRATOR_STATE(state, path, mnee) & PATH_MNEE_VALID)
+#endif
+  ) {
     float blur_pdf = kernel_data.integrator.filter_glossy *
                      INTEGRATOR_STATE(state, path, min_ray_pdf);
 
@@ -605,7 +609,8 @@ ccl_device void shader_eval_surface(KernelGlobals kg,
                                     ConstIntegratorGenericState state,
                                     ccl_private ShaderData *ccl_restrict sd,
                                     ccl_global float *ccl_restrict buffer,
-                                    uint32_t path_flag)
+                                    uint32_t path_flag,
+                                    bool use_caustics_storage = false)
 {
   /* If path is being terminated, we are tracing a shadow ray or evaluating
    * emission, then we don't need to store closures. The emission and shadow
@@ -615,7 +620,7 @@ ccl_device void shader_eval_surface(KernelGlobals kg,
     max_closures = 0;
   }
   else {
-    max_closures = kernel_data.max_closures;
+    max_closures = use_caustics_storage ? CAUSTICS_MAX_CLOSURE : kernel_data.max_closures;
   }
 
   sd->num_closure = 0;

intern/cycles/kernel/integrator/state_template.h

@@ -34,6 +34,8 @@ KERNEL_STRUCT_MEMBER(path, uint32_t, rng_hash, KERNEL_FEATURE_PATH_TRACING)
 KERNEL_STRUCT_MEMBER(path, uint16_t, rng_offset, KERNEL_FEATURE_PATH_TRACING)
 /* enum PathRayFlag */
 KERNEL_STRUCT_MEMBER(path, uint32_t, flag, KERNEL_FEATURE_PATH_TRACING)
+/* enum PathRayMNEE */
+KERNEL_STRUCT_MEMBER(path, uint8_t, mnee, KERNEL_FEATURE_PATH_TRACING)
 /* Multiple importance sampling
  * The PDF of BSDF sampling at the last scatter point, and distance to the
  * last scatter point minus the last ray segment. This distance lets us

intern/cycles/kernel/integrator/subsurface.h

@@ -171,7 +171,12 @@ ccl_device_inline bool subsurface_scatter(KernelGlobals kg, IntegratorState stat
 
   const int shader = intersection_get_shader(kg, &ss_isect.hits[0]);
   const int shader_flags = kernel_tex_fetch(__shaders, shader).flags;
-  if (shader_flags & SD_HAS_RAYTRACE) {
+  const int object_flags = intersection_get_object_flags(kg, &ss_isect.hits[0]);
+  const bool use_caustics = kernel_data.integrator.use_caustics &&
+                            (object_flags & SD_OBJECT_CAUSTICS);
+  const bool use_raytrace_kernel = (shader_flags & SD_HAS_RAYTRACE) || use_caustics;
+
+  if (use_raytrace_kernel) {
     INTEGRATOR_PATH_NEXT_SORTED(DEVICE_KERNEL_INTEGRATOR_INTERSECT_SUBSURFACE,
                                 DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE,
                                 shader);

intern/cycles/kernel/light/light.h

@@ -246,6 +246,15 @@ ccl_device bool lights_intersect(KernelGlobals kg,
       if (!(klight->shader_id & SHADER_USE_MIS)) {
         continue;
       }
+
+#ifdef __MNEE__
+      /* This path should have been resolved with mnee, it will
+       * generate a firefly for small lights since it is improbable. */
+      if ((INTEGRATOR_STATE(state, path, mnee) & PATH_MNEE_CULL_LIGHT_CONNECTION) &&
+          klight->use_caustics) {
+        continue;
+      }
+#endif
     }
 
     if (path_flag & PATH_RAY_SHADOW_CATCHER_PASS) {

intern/cycles/kernel/svm/closure.h

@@ -395,8 +395,10 @@ ccl_device_noinline int svm_node_closure_bsdf(KernelGlobals kg,
             if (kernel_data.integrator.caustics_refractive || (path_flag & PATH_RAY_DIFFUSE) == 0)
 #  endif
             {
+              /* This is to prevent mnee from receiving a null bsdf. */
+              float refraction_fresnel = fmaxf(0.0001f, 1.0f - fresnel);
               ccl_private MicrofacetBsdf *bsdf = (ccl_private MicrofacetBsdf *)bsdf_alloc(
-                  sd, sizeof(MicrofacetBsdf), base_color * glass_weight * (1.0f - fresnel));
+                  sd, sizeof(MicrofacetBsdf), base_color * glass_weight * refraction_fresnel);
               if (bsdf) {
                 bsdf->N = N;
                 bsdf->T = make_float3(0.0f, 0.0f, 0.0f);
@@ -674,8 +676,10 @@ ccl_device_noinline int svm_node_closure_bsdf(KernelGlobals kg,
       if (kernel_data.integrator.caustics_refractive || (path_flag & PATH_RAY_DIFFUSE) == 0)
 #endif
       {
+        /* This is to prevent mnee from receiving a null bsdf. */
+        float refraction_fresnel = fmaxf(0.0001f, 1.0f - fresnel);
         ccl_private MicrofacetBsdf *bsdf = (ccl_private MicrofacetBsdf *)bsdf_alloc(
-            sd, sizeof(MicrofacetBsdf), weight * (1.0f - fresnel));
+            sd, sizeof(MicrofacetBsdf), weight * refraction_fresnel);
 
         if (bsdf) {
           bsdf->N = N;

intern/cycles/kernel/types.h

@@ -102,6 +102,12 @@ CCL_NAMESPACE_BEGIN
 #  undef __BAKING__
 #endif /* __KERNEL_GPU_RAYTRACING__ */
 
+/* MNEE currently causes "Compute function exceeds available temporary registers"
+ * on Metal, disabled for now. */
+#ifndef __KERNEL_METAL__
+#  define __MNEE__
+#endif
+
 /* Scene-based selective features compilation. */
 #ifdef __KERNEL_FEATURES__
 #  if !(__KERNEL_FEATURES & KERNEL_FEATURE_CAMERA_MOTION)
@@ -293,6 +299,13 @@ enum PathRayFlag : uint32_t {
   PATH_RAY_SHADOW_CATCHER_BACKGROUND = (1U << 31U),
 };
 
+// 8bit enum, just in case we need to move more variables in it
+enum PathRayMNEE {
+  PATH_MNEE_VALID = (1U << 0U),
+  PATH_MNEE_RECEIVER_ANCESTOR = (1U << 1U),
+  PATH_MNEE_CULL_LIGHT_CONNECTION = (1U << 2U),
+};
+
 /* Configure ray visibility bits for rays and objects respectively,
  * to make shadow catchers work.
  *
@@ -649,6 +662,17 @@ typedef struct AttributeDescriptor {
 #  define MAX_CLOSURE __MAX_CLOSURE__
 #endif
 
+/* For manifold next event estimation, we need space to store and evaluate
+ * 2 closures (with extra data) on the refractive interfaces, in addition
+ * to keeping the full sd at the current shading point. We need 4 because a
+ * refractive bsdf is instanced with a companion reflection bsdf, even though
+ * we only need the refractive one, and each of them requires 2 slots. */
+#ifndef __CAUSTICS_MAX_CLOSURE__
+#  define CAUSTICS_MAX_CLOSURE 4
+#else
+#  define CAUSTICS_MAX_CLOSURE __CAUSTICS_MAX_CLOSURE__
+#endif
+
 #ifndef __MAX_VOLUME_STACK_SIZE__
 #  define MAX_VOLUME_STACK_SIZE 32
 #else
@@ -779,11 +803,18 @@ enum ShaderDataObjectFlag {
   SD_OBJECT_SHADOW_CATCHER = (1 << 7),
   /* object has volume attributes */
   SD_OBJECT_HAS_VOLUME_ATTRIBUTES = (1 << 8),
+  /* object is caustics caster */
+  SD_OBJECT_CAUSTICS_CASTER = (1 << 9),
+  /* object is caustics receiver */
+  SD_OBJECT_CAUSTICS_RECEIVER = (1 << 10),
+
+  /* object is using caustics */
+  SD_OBJECT_CAUSTICS = (SD_OBJECT_CAUSTICS_CASTER | SD_OBJECT_CAUSTICS_RECEIVER),
 
   SD_OBJECT_FLAGS = (SD_OBJECT_HOLDOUT_MASK | SD_OBJECT_MOTION | SD_OBJECT_TRANSFORM_APPLIED |
                      SD_OBJECT_NEGATIVE_SCALE_APPLIED | SD_OBJECT_HAS_VOLUME |
                      SD_OBJECT_INTERSECTS_VOLUME | SD_OBJECT_SHADOW_CATCHER |
-                     SD_OBJECT_HAS_VOLUME_ATTRIBUTES)
+                     SD_OBJECT_HAS_VOLUME_ATTRIBUTES | SD_OBJECT_CAUSTICS)
 };
 
 typedef struct ccl_align(16) ShaderData
@@ -882,6 +913,15 @@ typedef struct ccl_align(16) ShaderDataTinyStorage
   char pad[sizeof(ShaderData) - sizeof(ShaderClosure) * MAX_CLOSURE];
 }
 ShaderDataTinyStorage;
+
+/* ShaderDataCausticsStorage needs the same alignment as ShaderData, or else
+ * the pointer cast in AS_SHADER_DATA invokes undefined behavior. */
+typedef struct ccl_align(16) ShaderDataCausticsStorage
+{
+  char pad[sizeof(ShaderData) - sizeof(ShaderClosure) * (MAX_CLOSURE - CAUSTICS_MAX_CLOSURE)];
+}
+ShaderDataCausticsStorage;
+
 #define AS_SHADER_DATA(shader_data_tiny_storage) \
   ((ccl_private ShaderData *)shader_data_tiny_storage)
 
@@ -1201,6 +1241,9 @@ typedef struct KernelIntegrator {
   /* mis */
   int use_lamp_mis;
 
+  /* caustics */
+  int use_caustics;
+
   /* sampler */
   int sampling_pattern;
 
@@ -1219,7 +1262,7 @@ typedef struct KernelIntegrator {
   int direct_light_sampling_type;
 
   /* padding */
-  int pad1, pad2;
+  int pad1;
 } KernelIntegrator;
 static_assert_align(KernelIntegrator, 16);
 
@@ -1383,7 +1426,8 @@ typedef struct KernelLight {
   float max_bounces;
   float random;
   float strength[3];
-  float pad1, pad2;
+  int use_caustics;
+  float pad1;
   Transform tfm;
   Transform itfm;
   union {

intern/cycles/scene/light.cpp

@@ -123,6 +123,7 @@ NODE_DEFINE(Light)
   SOCKET_BOOLEAN(use_glossy, "Use Glossy", true);
   SOCKET_BOOLEAN(use_transmission, "Use Transmission", true);
   SOCKET_BOOLEAN(use_scatter, "Use Scatter", true);
+  SOCKET_BOOLEAN(use_caustics, "Shadow Caustics", false);
 
   SOCKET_INT(max_bounces, "Max Bounces", 1024);
   SOCKET_UINT(random_id, "Random ID", 0);
@@ -896,6 +897,7 @@ void LightManager::device_update_points(Device *, DeviceScene *dscene, Scene *sc
 
     klights[light_index].max_bounces = max_bounces;
     klights[light_index].random = random;
+    klights[light_index].use_caustics = light->use_caustics;
 
     klights[light_index].tfm = light->tfm;
     klights[light_index].itfm = transform_inverse(light->tfm);

intern/cycles/scene/light.h

@@ -61,6 +61,7 @@ class Light : public Node {
   NODE_SOCKET_API(bool, use_glossy)
   NODE_SOCKET_API(bool, use_transmission)
   NODE_SOCKET_API(bool, use_scatter)
+  NODE_SOCKET_API(bool, use_caustics)
 
   NODE_SOCKET_API(bool, is_shadow_catcher)
   NODE_SOCKET_API(bool, is_portal)

intern/cycles/scene/object.cpp

@@ -90,6 +90,9 @@ NODE_DEFINE(Object)
 
   SOCKET_BOOLEAN(is_shadow_catcher, "Shadow Catcher", false);
 
+  SOCKET_BOOLEAN(is_caustics_caster, "Cast Shadow Caustics", false);
+  SOCKET_BOOLEAN(is_caustics_receiver, "Receive Shadow Caustics", false);
+
   SOCKET_NODE(particle_system, "Particle System", ParticleSystem::get_node_type());
   SOCKET_INT(particle_index, "Particle Index", 0);
 
@@ -510,6 +513,14 @@ void ObjectManager::device_update_object_transform(UpdateObjectTransformState *s
   kobject.visibility = ob->visibility_for_tracing();
   kobject.primitive_type = geom->primitive_type();
 
+  /* Object shadow caustics flag */
+  if (ob->is_caustics_caster) {
+    flag |= SD_OBJECT_CAUSTICS_CASTER;
+  }
+  if (ob->is_caustics_receiver) {
+    flag |= SD_OBJECT_CAUSTICS_RECEIVER;
+  }
+
   /* Object flag. */
   if (ob->use_holdout) {
     flag |= SD_OBJECT_HOLDOUT_MASK;

intern/cycles/scene/object.h

@@ -55,6 +55,9 @@ class Object : public Node {
   NODE_SOCKET_API(float, shadow_terminator_shading_offset)
   NODE_SOCKET_API(float, shadow_terminator_geometry_offset)
 
+  NODE_SOCKET_API(bool, is_caustics_caster)
+  NODE_SOCKET_API(bool, is_caustics_receiver)
+
   NODE_SOCKET_API(float3, dupli_generated)
   NODE_SOCKET_API(float2, dupli_uv)
 

intern/cycles/scene/scene.cpp

@@ -489,7 +489,21 @@ void Scene::update_kernel_features()
   if (use_motion && camera->use_motion()) {
     kernel_features |= KERNEL_FEATURE_CAMERA_MOTION;
   }
+
+  /* Figure out whether the scene will use shader raytrace we need at least
+   * one caustic light, one caustic caster and one caustic receiver to use
+   * and enable the mnee code path. */
+  bool has_caustics_receiver = false;
+  bool has_caustics_caster = false;
+  bool has_caustics_light = false;
+
   foreach (Object *object, objects) {
+    if (object->get_is_caustics_caster()) {
+      has_caustics_caster = true;
+    }
+    else if (object->get_is_caustics_receiver()) {
+      has_caustics_receiver = true;
+    }
     Geometry *geom = object->get_geometry();
     if (use_motion) {
       if (object->use_motion() || geom->get_use_motion_blur()) {
@@ -518,6 +532,18 @@ void Scene::update_kernel_features()
     }
   }
 
+  foreach (Light *light, lights) {
+    if (light->get_use_caustics()) {
+      has_caustics_light = true;
+    }
+  }
+
+  dscene.data.integrator.use_caustics = false;
+  if (has_caustics_caster && has_caustics_receiver && has_caustics_light) {
+    dscene.data.integrator.use_caustics = true;
+    kernel_features |= KERNEL_FEATURE_NODE_RAYTRACE;
+  }
+
   if (bake_manager->get_baking()) {
     kernel_features |= KERNEL_FEATURE_BAKING;
   }

intern/cycles/util/math_float2.h

@@ -40,7 +40,7 @@ ccl_device_inline float average(const float2 &a);
 ccl_device_inline float distance(const float2 &a, const float2 &b);
 ccl_device_inline float dot(const float2 &a, const float2 &b);
 ccl_device_inline float cross(const float2 &a, const float2 &b);
-ccl_device_inline float len(const float2 &a);
+ccl_device_inline float len(const float2 a);
 ccl_device_inline float2 normalize(const float2 &a);
 ccl_device_inline float2 normalize_len(const float2 &a, float *t);
 ccl_device_inline float2 safe_normalize(const float2 &a);
@@ -187,11 +187,6 @@ ccl_device_inline float cross(const float2 &a, const float2 &b)
   return (a.x * b.y - a.y * b.x);
 }
 
-ccl_device_inline float len(const float2 &a)
-{
-  return sqrtf(dot(a, a));
-}
-
 ccl_device_inline float2 normalize(const float2 &a)
 {
   return a / len(a);
@@ -251,6 +246,11 @@ ccl_device_inline float2 floor(const float2 &a)
 
 #endif /* !__KERNEL_METAL__ */
 
+ccl_device_inline float len(const float2 a)
+{
+  return sqrtf(dot(a, a));
+}
+
 ccl_device_inline float2 safe_divide_float2_float(const float2 a, const float b)
 {
   return (b != 0.0f) ? a / b : zero_float2();

tests/python/cycles_render_tests.py

@@ -32,6 +32,11 @@ BLACKLIST_OPTIX = [
     'T43865.blend',
 ]
 
+BLACKLIST_METAL = [
+    # No MNEE for Metal currently
+    "underwater_caustics.blend",
+]
+
 BLACKLIST_GPU = [
     # Uninvestigated differences with GPU.
     'image_log.blend',
@@ -116,6 +121,8 @@ def main():
         blacklist += BLACKLIST_OSL
     if device == 'OPTIX':
         blacklist += BLACKLIST_OPTIX
+    if device == 'METAL':
+        blacklist += BLACKLIST_METAL
 
     from modules import render_report
     report = render_report.Report('Cycles', output_dir, idiff, device, blacklist)