Cycles: Add optional probabilistic termination of light samples based on their expected contribution

In scenes with many lights, some of them might have a very small contribution to some pixels, but the shadow rays are traced anyways.
To avoid that, this patch adds probabilistic termination to light samples - if the contribution before checking for shadowing is below a user-defined threshold, the sample will be discarded with probability (1 - (contribution / threshold)) and otherwise kept, but weighted more to remain unbiased.
This is the same approach that's also used in path termination based on length.

Note that the rendering remains unbiased with this option, it just adds a bit of noise - but if the setting is used moderately, the speedup gained easily outweighs the additional noise.

Reviewers: #cycles

Subscribers: sergey, brecht

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

15 files changed, 109 insertions, 20 deletions

diff --git a/intern/cycles/blender/addon/properties.py b/intern/cycles/blender/addon/properties.py
index 977d7f75bb7..ede3ece5bce 100644
--- a/intern/cycles/blender/addon/properties.py
+++ b/intern/cycles/blender/addon/properties.py
@@ -266,6 +266,13 @@ class CyclesRenderSettings(bpy.types.PropertyGroup):
                 description="Sample all lights (for indirect samples), rather than randomly picking one",
                 default=True,
                 )
+        cls.light_sampling_threshold = FloatProperty(
+                name="Light Sampling Threshold",
+                description="Probabilistically terminate light samples when the light contribution is below this threshold (more noise but faster rendering). "
+                            "Zero disables the test and never ignores lights.",
+                min=0.0, max=1.0,
+                default=0.05,
+                )
 
         cls.caustics_reflective = BoolProperty(
                 name="Reflective Caustics",
diff --git a/intern/cycles/blender/addon/ui.py b/intern/cycles/blender/addon/ui.py
index 52872d2b83f..4942a71ce4d 100644
--- a/intern/cycles/blender/addon/ui.py
+++ b/intern/cycles/blender/addon/ui.py
@@ -166,6 +166,7 @@ class CyclesRender_PT_sampling(CyclesButtonsPanel, Panel):
 
         sub.prop(cscene, "sample_clamp_direct")
         sub.prop(cscene, "sample_clamp_indirect")
+        sub.prop(cscene, "light_sampling_threshold")
 
         if cscene.progressive == 'PATH' or use_branched_path(context) is False:
             col = split.column()
diff --git a/intern/cycles/blender/blender_sync.cpp b/intern/cycles/blender/blender_sync.cpp
index 4ca202ac40d..bc5c3bb8096 100644
--- a/intern/cycles/blender/blender_sync.cpp
+++ b/intern/cycles/blender/blender_sync.cpp
@@ -284,6 +284,7 @@ void BlenderSync::sync_integrator()
 
 	integrator->sample_all_lights_direct = get_boolean(cscene, "sample_all_lights_direct");
 	integrator->sample_all_lights_indirect = get_boolean(cscene, "sample_all_lights_indirect");
+	integrator->light_sampling_threshold = get_float(cscene, "light_sampling_threshold");
 
 	int diffuse_samples = get_int(cscene, "diffuse_samples");
 	int glossy_samples = get_int(cscene, "glossy_samples");
diff --git a/intern/cycles/kernel/kernel_accumulate.h b/intern/cycles/kernel/kernel_accumulate.h
index 623c1dcaaa1..6c3ee6b8098 100644
--- a/intern/cycles/kernel/kernel_accumulate.h
+++ b/intern/cycles/kernel/kernel_accumulate.h
@@ -96,7 +96,7 @@ ccl_device_inline bool bsdf_eval_is_zero(BsdfEval *eval)
 	}
 }
 
-ccl_device_inline void bsdf_eval_mul(BsdfEval *eval, float3 value)
+ccl_device_inline void bsdf_eval_mul(BsdfEval *eval, float value)
 {
 #ifdef __PASSES__
 	if(eval->use_light_pass) {
@@ -115,6 +115,36 @@ ccl_device_inline void bsdf_eval_mul(BsdfEval *eval, float3 value)
 	}
 }
 
+ccl_device_inline void bsdf_eval_mul3(BsdfEval *eval, float3 value)
+{
+#ifdef __PASSES__
+	if(eval->use_light_pass) {
+		eval->diffuse *= value;
+		eval->glossy *= value;
+		eval->transmission *= value;
+		eval->subsurface *= value;
+		eval->scatter *= value;
+
+		/* skipping transparent, this function is used by for eval(), will be zero then */
+	}
+	else
+		eval->diffuse *= value;
+#else
+	eval->diffuse *= value;
+#endif
+}
+
+ccl_device_inline float3 bsdf_eval_sum(BsdfEval *eval)
+{
+#ifdef __PASSES__
+	if(eval->use_light_pass) {
+		return eval->diffuse + eval->glossy + eval->transmission + eval->subsurface + eval->scatter;
+	}
+	else
+#endif
+	return eval->diffuse;
+}
+
 /* Path Radiance
  *
  * We accumulate different render passes separately. After summing at the end
@@ -193,8 +223,7 @@ ccl_device_inline void path_radiance_bsdf_bounce(PathRadiance *L, ccl_addr_space
 		}
 		else {
 			/* transparent bounce before first hit, or indirectly visible through BSDF */
-			float3 sum = (bsdf_eval->diffuse + bsdf_eval->glossy + bsdf_eval->transmission + bsdf_eval->transparent +
-						  bsdf_eval->subsurface + bsdf_eval->scatter) * inverse_pdf;
+			float3 sum = (bsdf_eval_sum(bsdf_eval) + bsdf_eval->transparent) * inverse_pdf;
 			*throughput *= sum;
 		}
 	}
@@ -264,8 +293,7 @@ ccl_device_inline void path_radiance_accum_light(PathRadiance *L, float3 through
 		}
 		else {
 			/* indirectly visible lighting after BSDF bounce */
-			float3 sum = bsdf_eval->diffuse + bsdf_eval->glossy + bsdf_eval->transmission + bsdf_eval->subsurface + bsdf_eval->scatter;
-			L->indirect += throughput*sum*shadow;
+			L->indirect += throughput*bsdf_eval_sum(bsdf_eval)*shadow;
 		}
 	}
 	else
diff --git a/intern/cycles/kernel/kernel_emission.h b/intern/cycles/kernel/kernel_emission.h
index 9e4a631b998..8c7c651a053 100644
--- a/intern/cycles/kernel/kernel_emission.h
+++ b/intern/cycles/kernel/kernel_emission.h
@@ -94,7 +94,8 @@ ccl_device_noinline bool direct_emission(KernelGlobals *kg,
                                          ccl_addr_space PathState *state,
                                          Ray *ray,
                                          BsdfEval *eval,
-                                         bool *is_lamp)
+                                         bool *is_lamp,
+                                         float rand_terminate)
 {
 	if(ls->pdf == 0.0f)
 		return false;
@@ -134,7 +135,7 @@ ccl_device_noinline bool direct_emission(KernelGlobals *kg,
 	shader_bsdf_eval(kg, sd, ls->D, eval, ls->pdf, ls->shader & SHADER_USE_MIS);
 #endif
 
-	bsdf_eval_mul(eval, light_eval/ls->pdf);
+	bsdf_eval_mul3(eval, light_eval/ls->pdf);
 
 #ifdef __PASSES__
 	/* use visibility flag to skip lights */
@@ -155,6 +156,16 @@ ccl_device_noinline bool direct_emission(KernelGlobals *kg,
 	if(bsdf_eval_is_zero(eval))
 		return false;
 
+	if(kernel_data.integrator.light_inv_rr_threshold > 0.0f) {
+		float probability = max3(bsdf_eval_sum(eval)) * kernel_data.integrator.light_inv_rr_threshold;
+		if(probability < 1.0f) {
+			if(rand_terminate >= probability) {
+				return false;
+			}
+			bsdf_eval_mul(eval, 1.0f / probability);
+		}
+	}
+
 	if(ls->shader & SHADER_CAST_SHADOW) {
 		/* setup ray */
 		bool transmit = (dot(ccl_fetch(sd, Ng), ls->D) < 0.0f);
diff --git a/intern/cycles/kernel/kernel_path.h b/intern/cycles/kernel/kernel_path.h
index 7558fb94478..7ef79815ad5 100644
--- a/intern/cycles/kernel/kernel_path.h
+++ b/intern/cycles/kernel/kernel_path.h
@@ -851,7 +851,6 @@ ccl_device_inline float4 kernel_path_integrate(KernelGlobals *kg,
 		}
 		else if(probability != 1.0f) {
 			float terminate = path_state_rng_1D_for_decision(kg, rng, &state, PRNG_TERMINATE);
-
 			if(terminate >= probability)
 				break;
 
diff --git a/intern/cycles/kernel/kernel_path_surface.h b/intern/cycles/kernel/kernel_path_surface.h
index 45f0f1cbfaa..fea503d06e5 100644
--- a/intern/cycles/kernel/kernel_path_surface.h
+++ b/intern/cycles/kernel/kernel_path_surface.h
@@ -49,6 +49,7 @@ ccl_device_noinline void kernel_branched_path_surface_connect_light(KernelGlobal
 			for(int j = 0; j < num_samples; j++) {
 				float light_u, light_v;
 				path_branched_rng_2D(kg, &lamp_rng, state, j, num_samples, PRNG_LIGHT_U, &light_u, &light_v);
+				float terminate = path_branched_rng_light_termination(kg, &lamp_rng, state, j, num_samples);
 
 				LightSample ls;
 				if(lamp_light_sample(kg, i, light_u, light_v, ccl_fetch(sd, P), &ls)) {
@@ -57,7 +58,7 @@ ccl_device_noinline void kernel_branched_path_surface_connect_light(KernelGlobal
 					if(kernel_data.integrator.pdf_triangles != 0.0f)
 						ls.pdf *= 2.0f;
 
-					if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp)) {
+					if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 						/* trace shadow ray */
 						float3 shadow;
 
@@ -79,6 +80,7 @@ ccl_device_noinline void kernel_branched_path_surface_connect_light(KernelGlobal
 				float light_t = path_branched_rng_1D(kg, rng, state, j, num_samples, PRNG_LIGHT);
 				float light_u, light_v;
 				path_branched_rng_2D(kg, rng, state, j, num_samples, PRNG_LIGHT_U, &light_u, &light_v);
+				float terminate = path_branched_rng_light_termination(kg, rng, state, j, num_samples);
 
 				/* only sample triangle lights */
 				if(kernel_data.integrator.num_all_lights)
@@ -90,7 +92,7 @@ ccl_device_noinline void kernel_branched_path_surface_connect_light(KernelGlobal
 					if(kernel_data.integrator.num_all_lights)
 						ls.pdf *= 2.0f;
 
-					if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp)) {
+					if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 						/* trace shadow ray */
 						float3 shadow;
 
@@ -108,11 +110,12 @@ ccl_device_noinline void kernel_branched_path_surface_connect_light(KernelGlobal
 		float light_t = path_state_rng_1D(kg, rng, state, PRNG_LIGHT);
 		float light_u, light_v;
 		path_state_rng_2D(kg, rng, state, PRNG_LIGHT_U, &light_u, &light_v);
+		float terminate = path_state_rng_light_termination(kg, rng, state);
 
 		LightSample ls;
 		if(light_sample(kg, light_t, light_u, light_v, ccl_fetch(sd, time), ccl_fetch(sd, P), state->bounce, &ls)) {
 			/* sample random light */
-			if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp)) {
+			if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 				/* trace shadow ray */
 				float3 shadow;
 
@@ -210,7 +213,8 @@ ccl_device_inline void kernel_path_surface_connect_light(KernelGlobals *kg, ccl_
 
 	LightSample ls;
 	if(light_sample(kg, light_t, light_u, light_v, ccl_fetch(sd, time), ccl_fetch(sd, P), state->bounce, &ls)) {
-		if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp)) {
+		float terminate = path_state_rng_light_termination(kg, rng, state);
+		if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 			/* trace shadow ray */
 			float3 shadow;
 
diff --git a/intern/cycles/kernel/kernel_path_volume.h b/intern/cycles/kernel/kernel_path_volume.h
index 5ee1912c913..3d3b7385d8b 100644
--- a/intern/cycles/kernel/kernel_path_volume.h
+++ b/intern/cycles/kernel/kernel_path_volume.h
@@ -48,7 +48,8 @@ ccl_device_inline void kernel_path_volume_connect_light(
 
 	if(light_sample(kg, light_t, light_u, light_v, sd->time, sd->P, state->bounce, &ls))
 	{
-		if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp)) {
+		float terminate = path_state_rng_light_termination(kg, rng, state);
+		if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 			/* trace shadow ray */
 			float3 shadow;
 
@@ -161,7 +162,8 @@ ccl_device void kernel_branched_path_volume_connect_light(KernelGlobals *kg, RNG
 					if(kernel_data.integrator.pdf_triangles != 0.0f)
 						ls.pdf *= 2.0f;
 
-					if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp)) {
+					float terminate = path_branched_rng_light_termination(kg, rng, state, j, num_samples);
+					if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 						/* trace shadow ray */
 						float3 shadow;
 
@@ -209,7 +211,8 @@ ccl_device void kernel_branched_path_volume_connect_light(KernelGlobals *kg, RNG
 					if(kernel_data.integrator.num_all_lights)
 						ls.pdf *= 2.0f;
 
-					if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp)) {
+					float terminate = path_branched_rng_light_termination(kg, rng, state, j, num_samples);
+					if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 						/* trace shadow ray */
 						float3 shadow;
 
@@ -246,7 +249,8 @@ ccl_device void kernel_branched_path_volume_connect_light(KernelGlobals *kg, RNG
 		/* todo: split up light_sample so we don't have to call it again with new position */
 		if(light_sample(kg, light_t, light_u, light_v, sd->time, sd->P, state->bounce, &ls)) {
 			/* sample random light */
-			if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp)) {
+			float terminate = path_state_rng_light_termination(kg, rng, state);
+			if(direct_emission(kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 				/* trace shadow ray */
 				float3 shadow;
 
diff --git a/intern/cycles/kernel/kernel_random.h b/intern/cycles/kernel/kernel_random.h
index 4a76ffddbe7..2372b07d974 100644
--- a/intern/cycles/kernel/kernel_random.h
+++ b/intern/cycles/kernel/kernel_random.h
@@ -300,6 +300,23 @@ ccl_device_inline void path_branched_rng_2D(KernelGlobals *kg, ccl_addr_space RN
 	path_rng_2D(kg, rng, state->sample*num_branches + branch, state->num_samples*num_branches, state->rng_offset + dimension, fx, fy);
 }
 
+/* Utitility functions to get light termination value, since it might not be needed in many cases. */
+ccl_device_inline float path_state_rng_light_termination(KernelGlobals *kg, ccl_addr_space RNG *rng, const PathState *state)
+{
+	if(kernel_data.integrator.light_inv_rr_threshold > 0.0f) {
+		return path_state_rng_1D_for_decision(kg, rng, state, PRNG_LIGHT_TERMINATE);
+	}
+	return 0.0f;
+}
+
+ccl_device_inline float path_branched_rng_light_termination(KernelGlobals *kg, ccl_addr_space RNG *rng, const PathState *state, int branch, int num_branches)
+{
+	if(kernel_data.integrator.light_inv_rr_threshold > 0.0f) {
+		return path_branched_rng_1D_for_decision(kg, rng, state, branch, num_branches, PRNG_LIGHT_TERMINATE);
+	}
+	return 0.0f;
+}
+
 ccl_device_inline void path_state_branch(PathState *state, int branch, int num_branches)
 {
 	/* path is splitting into a branch, adjust so that each branch
diff --git a/intern/cycles/kernel/kernel_shader.h b/intern/cycles/kernel/kernel_shader.h
index a8070a133de..b4b980c4e90 100644
--- a/intern/cycles/kernel/kernel_shader.h
+++ b/intern/cycles/kernel/kernel_shader.h
@@ -572,7 +572,7 @@ void shader_bsdf_eval(KernelGlobals *kg,
 		_shader_bsdf_multi_eval(kg, sd, omega_in, &pdf, -1, eval, 0.0f, 0.0f);
 		if(use_mis) {
 			float weight = power_heuristic(light_pdf, pdf);
-			bsdf_eval_mul(eval, make_float3(weight, weight, weight));
+			bsdf_eval_mul(eval, weight);
 		}
 	}
 }
diff --git a/intern/cycles/kernel/kernel_types.h b/intern/cycles/kernel/kernel_types.h
index 91b86618745..a6c31d4a518 100644
--- a/intern/cycles/kernel/kernel_types.h
+++ b/intern/cycles/kernel/kernel_types.h
@@ -250,7 +250,7 @@ enum PathTraceDimension {
 	PRNG_LIGHT = 3,
 	PRNG_LIGHT_U = 4,
 	PRNG_LIGHT_V = 5,
-	PRNG_UNUSED_3 = 6,
+	PRNG_LIGHT_TERMINATE = 6,
 	PRNG_TERMINATE = 7,
 
 #ifdef __VOLUME__
@@ -1123,8 +1123,9 @@ typedef struct KernelIntegrator {
 	float volume_step_size;
 	int volume_samples;
 
+	float light_inv_rr_threshold;
+
 	int pad1;
-	int pad2;
 } KernelIntegrator;
 static_assert_align(KernelIntegrator, 16);
 
diff --git a/intern/cycles/kernel/split/kernel_direct_lighting.h b/intern/cycles/kernel/split/kernel_direct_lighting.h
index 6ad736fc2c1..82ca18829d3 100644
--- a/intern/cycles/kernel/split/kernel_direct_lighting.h
+++ b/intern/cycles/kernel/split/kernel_direct_lighting.h
@@ -72,6 +72,7 @@ ccl_device char kernel_direct_lighting(
 			float light_t = path_state_rng_1D(kg, rng, state, PRNG_LIGHT);
 			float light_u, light_v;
 			path_state_rng_2D(kg, rng, state, PRNG_LIGHT_U, &light_u, &light_v);
+			float terminate = path_state_rng_light_termination(kg, rng, state);
 
 			LightSample ls;
 			if(light_sample(kg,
@@ -88,7 +89,7 @@ ccl_device char kernel_direct_lighting(
 
 				BsdfEval L_light;
 				bool is_lamp;
-				if(direct_emission(kg, sd, kg->sd_input, &ls, state, &light_ray, &L_light, &is_lamp)) {
+				if(direct_emission(kg, sd, kg->sd_input, &ls, state, &light_ray, &L_light, &is_lamp, terminate)) {
 					/* Write intermediate data to global memory to access from
 					 * the next kernel.
 					 */
diff --git a/intern/cycles/render/integrator.cpp b/intern/cycles/render/integrator.cpp
index 63914e57319..a9a33d2e789 100644
--- a/intern/cycles/render/integrator.cpp
+++ b/intern/cycles/render/integrator.cpp
@@ -65,6 +65,7 @@ NODE_DEFINE(Integrator)
 
 	SOCKET_BOOLEAN(sample_all_lights_direct, "Sample All Lights Direct", true);
 	SOCKET_BOOLEAN(sample_all_lights_indirect, "Sample All Lights Indirect", true);
+	SOCKET_FLOAT(light_sampling_threshold, "Light Sampling Threshold", 0.05f);
 
 	static NodeEnum method_enum;
 	method_enum.insert("path", PATH);
@@ -164,6 +165,13 @@ void Integrator::device_update(Device *device, DeviceScene *dscene, Scene *scene
 	kintegrator->sampling_pattern = sampling_pattern;
 	kintegrator->aa_samples = aa_samples;
 
+	if(light_sampling_threshold > 0.0f) {
+		kintegrator->light_inv_rr_threshold = 1.0f / light_sampling_threshold;
+	}
+	else {
+		kintegrator->light_inv_rr_threshold = 0.0f;
+	}
+
 	/* sobol directions table */
 	int max_samples = 1;
 
diff --git a/intern/cycles/render/integrator.h b/intern/cycles/render/integrator.h
index 39eaaf246d4..17fdd0ef1db 100644
--- a/intern/cycles/render/integrator.h
+++ b/intern/cycles/render/integrator.h
@@ -64,8 +64,10 @@ public:
 	int mesh_light_samples;
 	int subsurface_samples;
 	int volume_samples;
+
 	bool sample_all_lights_direct;
 	bool sample_all_lights_indirect;
+	float light_sampling_threshold;
 
 	enum Method {
 		BRANCHED_PATH = 0,
diff --git a/intern/cycles/util/util_math.h b/intern/cycles/util/util_math.h
index e2abfcde702..3f4d3e06c0b 100644
--- a/intern/cycles/util/util_math.h
+++ b/intern/cycles/util/util_math.h
@@ -162,6 +162,11 @@ ccl_device_inline float max4(float a, float b, float c, float d)
 	return max(max(a, b), max(c, d));
 }
 
+ccl_device_inline float max3(float3 a)
+{
+	return max(max(a.x, a.y), a.z);
+}
+
 #ifndef __KERNEL_OPENCL__
 
 ccl_device_inline int clamp(int a, int mn, int mx)