Cycles: reduce noise using regular path tracing + subsurface scattering with

new cubic and gaussian falloff. Like the branched path tracer, this will now
shade all intersection points instead of using one at random.

5 files changed, 362 insertions, 224 deletions

diff --git a/intern/cycles/blender/blender_sync.cpp b/intern/cycles/blender/blender_sync.cpp
index 38a6c8eac91..4a686487462 100644
--- a/intern/cycles/blender/blender_sync.cpp
+++ b/intern/cycles/blender/blender_sync.cpp
@@ -189,7 +189,7 @@ void BlenderSync::sync_integrator()
 	}
 #endif
 
-	integrator->method = (Integrator::Method)get_int(cscene, "progressive");
+	integrator->method = (Integrator::Method)get_enum(cscene, "progressive");
 
 	int diffuse_samples = get_int(cscene, "diffuse_samples");
 	int glossy_samples = get_int(cscene, "glossy_samples");
diff --git a/intern/cycles/kernel/kernel_path.h b/intern/cycles/kernel/kernel_path.h
index cdf80aed4d1..8cb83e0797e 100644
--- a/intern/cycles/kernel/kernel_path.h
+++ b/intern/cycles/kernel/kernel_path.h
@@ -231,52 +231,24 @@ __device_inline bool shadow_blocked(KernelGlobals *kg, PathState *state, Ray *ra
 	return result;
 }
 
-__device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, Ray ray, __global float *buffer)
-{
-	/* initialize */
-	PathRadiance L;
-	float3 throughput = make_float3(1.0f, 1.0f, 1.0f);
-	float L_transparent = 0.0f;
 
-	path_radiance_init(&L, kernel_data.film.use_light_pass);
+#if defined(__BRANCHED_PATH__) || defined(__BSSRDF__)
 
-	float min_ray_pdf = FLT_MAX;
-	float ray_pdf = 0.0f;
+__device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray ray, __global float *buffer,
+	float3 throughput, int num_samples, int num_total_samples,
+	float min_ray_pdf, float ray_pdf, PathState state, int rng_offset, PathRadiance *L)
+{
 #ifdef __LAMP_MIS__
 	float ray_t = 0.0f;
 #endif
-	PathState state;
-	int rng_offset = PRNG_BASE_NUM;
-#ifdef __CMJ__
-	int num_samples = kernel_data.integrator.aa_samples;
-#else
-	int num_samples = 0;
-#endif
-
-	path_state_init(&state);
 
 	/* path iteration */
 	for(;; rng_offset += PRNG_BOUNCE_NUM) {
 		/* intersect scene */
 		Intersection isect;
 		uint visibility = path_state_ray_visibility(kg, &state);
-
 #ifdef __HAIR__
-		float difl = 0.0f, extmax = 0.0f;
-		uint lcg_state = 0;
-
-		if(kernel_data.bvh.have_curves) {
-			if((kernel_data.cam.resolution == 1) && (state.flag & PATH_RAY_CAMERA)) {	
-				float3 pixdiff = ray.dD.dx + ray.dD.dy;
-				/*pixdiff = pixdiff - dot(pixdiff, ray.D)*ray.D;*/
-				difl = kernel_data.curve.minimum_width * len(pixdiff) * 0.5f;
-			}
-
-			extmax = kernel_data.curve.maximum_width;
-			lcg_state = lcg_init(*rng + rng_offset + sample*0x51633e2d);
-		}
-
-		bool hit = scene_intersect(kg, &ray, visibility, &isect, &lcg_state, difl, extmax);
+		bool hit = scene_intersect(kg, &ray, visibility, &isect, NULL, 0.0f, 0.0f);
 #else
 		bool hit = scene_intersect(kg, &ray, visibility, &isect);
 #endif
@@ -295,29 +267,19 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 			light_ray.dP = ray.dP;
 
 			/* intersect with lamp */
-			float light_t = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT);
+			float light_t = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_LIGHT);
 			float3 emission;
 
 			if(indirect_lamp_emission(kg, &light_ray, state.flag, ray_pdf, light_t, &emission, state.bounce))
-				path_radiance_accum_emission(&L, throughput, emission, state.bounce);
+				path_radiance_accum_emission(L, throughput, emission, state.bounce);
 		}
 #endif
 
 		if(!hit) {
-			/* eval background shader if nothing hit */
-			if(kernel_data.background.transparent && (state.flag & PATH_RAY_CAMERA)) {
-				L_transparent += average(throughput);
-
-#ifdef __PASSES__
-				if(!(kernel_data.film.pass_flag & PASS_BACKGROUND))
-#endif
-					break;
-			}
-
 #ifdef __BACKGROUND__
 			/* sample background shader */
 			float3 L_background = indirect_background(kg, &ray, state.flag, ray_pdf, state.bounce);
-			path_radiance_accum_background(&L, throughput, L_background, state.bounce);
+			path_radiance_accum_background(L, throughput, L_background, state.bounce);
 #endif
 
 			break;
@@ -326,31 +288,9 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 		/* setup shading */
 		ShaderData sd;
 		shader_setup_from_ray(kg, &sd, &isect, &ray, state.bounce);
-		float rbsdf = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF);
-		shader_eval_surface(kg, &sd, rbsdf, state.flag, SHADER_CONTEXT_MAIN);
-
-		/* holdout */
-#ifdef __HOLDOUT__
-		if((sd.flag & (SD_HOLDOUT|SD_HOLDOUT_MASK)) && (state.flag & PATH_RAY_CAMERA)) {
-			if(kernel_data.background.transparent) {
-				float3 holdout_weight;
-				
-				if(sd.flag & SD_HOLDOUT_MASK)
-					holdout_weight = make_float3(1.0f, 1.0f, 1.0f);
-				else
-					holdout_weight = shader_holdout_eval(kg, &sd);
-
-				/* any throughput is ok, should all be identical here */
-				L_transparent += average(holdout_weight*throughput);
-			}
-
-			if(sd.flag & SD_HOLDOUT_MASK)
-				break;
-		}
-#endif
-
-		/* holdout mask objects do not write data passes */
-		kernel_write_data_passes(kg, buffer, &L, &sd, sample, state.flag, throughput);
+		float rbsdf = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_BSDF);
+		shader_eval_surface(kg, &sd, rbsdf, state.flag, SHADER_CONTEXT_INDIRECT);
+		shader_merge_closures(kg, &sd);
 
 		/* blurring of bsdf after bounces, for rays that have a small likelihood
 		 * of following this particular path (diffuse, rough glossy) */
@@ -366,22 +306,21 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 #ifdef __EMISSION__
 		/* emission */
 		if(sd.flag & SD_EMISSION) {
-			/* todo: is isect.t wrong here for transparent surfaces? */
 			float3 emission = indirect_primitive_emission(kg, &sd, isect.t, state.flag, ray_pdf);
-			path_radiance_accum_emission(&L, throughput, emission, state.bounce);
+			path_radiance_accum_emission(L, throughput, emission, state.bounce);
 		}
 #endif
 
 		/* path termination. this is a strange place to put the termination, it's
 		 * mainly due to the mixed in MIS that we use. gives too many unneeded
 		 * shader evaluations, only need emission if we are going to terminate */
-		float probability = path_state_terminate_probability(kg, &state, throughput);
+		float probability = path_state_terminate_probability(kg, &state, throughput*num_samples);
 
 		if(probability == 0.0f) {
 			break;
 		}
 		else if(probability != 1.0f) {
-			float terminate = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_TERMINATE);
+			float terminate = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_TERMINATE);
 
 			if(terminate >= probability)
 				break;
@@ -389,38 +328,11 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 			throughput /= probability;
 		}
 
-#ifdef __SUBSURFACE__
-		/* bssrdf scatter to a different location on the same object, replacing
-		 * the closures with a diffuse BSDF */
-		if(sd.flag & SD_BSSRDF) {
-			float bssrdf_probability;
-			ShaderClosure *sc = subsurface_scatter_pick_closure(kg, &sd, &bssrdf_probability);
-
-			/* modify throughput for picking bssrdf or bsdf */
-			throughput *= bssrdf_probability;
-
-			/* do bssrdf scatter step if we picked a bssrdf closure */
-			if(sc) {
-				uint lcg_state = lcg_init(*rng + rng_offset + sample*0x68bc21eb);
-
-				if(old_subsurface_scatter_use(&sd)) {
-					old_subsurface_scatter_step(kg, &sd, state.flag, sc, &lcg_state, false);
-				}
-				else {
-					float bssrdf_u, bssrdf_v;
-					path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF_U, &bssrdf_u, &bssrdf_v);
-					subsurface_scatter_step(kg, &sd, state.flag, sc, &lcg_state, bssrdf_u, bssrdf_v, false);
-				}
-			}
-		}
-#endif
-
 #ifdef __AO__
 		/* ambient occlusion */
 		if(kernel_data.integrator.use_ambient_occlusion || (sd.flag & SD_AO)) {
-			/* todo: solve correlation */
 			float bsdf_u, bsdf_v;
-			path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
+			path_rng_2D(kg, rng, sample, num_total_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
 
 			float ao_factor = kernel_data.background.ao_factor;
 			float3 ao_N;
@@ -444,7 +356,35 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 				light_ray.dD = differential3_zero();
 
 				if(!shadow_blocked(kg, &state, &light_ray, &ao_shadow))
-					path_radiance_accum_ao(&L, throughput, ao_bsdf, ao_shadow, state.bounce);
+					path_radiance_accum_ao(L, throughput, ao_bsdf, ao_shadow, state.bounce);
+			}
+		}
+#endif
+
+#ifdef __SUBSURFACE__
+		/* bssrdf scatter to a different location on the same object, replacing
+		 * the closures with a diffuse BSDF */
+		if(sd.flag & SD_BSSRDF) {
+			float bssrdf_probability;
+			ShaderClosure *sc = subsurface_scatter_pick_closure(kg, &sd, &bssrdf_probability);
+
+			/* modify throughput for picking bssrdf or bsdf */
+			throughput *= bssrdf_probability;
+
+			/* do bssrdf scatter step if we picked a bssrdf closure */
+			if(sc) {
+				uint lcg_state = lcg_init(*rng + rng_offset + sample*0x68bc21eb);
+
+				if(old_subsurface_scatter_use(&sd)) {
+					old_subsurface_scatter_step(kg, &sd, state.flag, sc, &lcg_state, false);
+				}
+				else {
+					float bssrdf_u, bssrdf_v;
+					path_rng_2D(kg, rng, sample, num_total_samples, rng_offset + PRNG_BSDF_U, &bssrdf_u, &bssrdf_v);
+					subsurface_scatter_step(kg, &sd, state.flag, sc, &lcg_state, bssrdf_u, bssrdf_v, false);
+				}
+
+				state.flag |= PATH_RAY_BSSRDF_ANCESTOR;
 			}
 		}
 #endif
@@ -453,14 +393,14 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 		if(kernel_data.integrator.use_direct_light) {
 			/* sample illumination from lights to find path contribution */
 			if(sd.flag & SD_BSDF_HAS_EVAL) {
-				float light_t = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT);
+				float light_t = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_LIGHT);
 #ifdef __MULTI_CLOSURE__
 				float light_o = 0.0f;
 #else
-				float light_o = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT_F);
+				float light_o = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_LIGHT_F);
 #endif
 				float light_u, light_v;
-				path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT_U, &light_u, &light_v);
+				path_rng_2D(kg, rng, sample, num_total_samples, rng_offset + PRNG_LIGHT_U, &light_u, &light_v);
 
 				Ray light_ray;
 				BsdfEval L_light;
@@ -470,13 +410,14 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 				light_ray.time = sd.time;
 #endif
 
+				/* sample random light */
 				if(direct_emission(kg, &sd, -1, light_t, light_o, light_u, light_v, &light_ray, &L_light, &is_lamp, state.bounce)) {
 					/* trace shadow ray */
 					float3 shadow;
 
 					if(!shadow_blocked(kg, &state, &light_ray, &shadow)) {
 						/* accumulate */
-						path_radiance_accum_light(&L, throughput, &L_light, shadow, 1.0f, state.bounce, is_lamp);
+						path_radiance_accum_light(L, throughput, &L_light, shadow, 1.0f, state.bounce, is_lamp);
 					}
 				}
 			}
@@ -493,7 +434,7 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 		float3 bsdf_omega_in;
 		differential3 bsdf_domega_in;
 		float bsdf_u, bsdf_v;
-		path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
+		path_rng_2D(kg, rng, sample, num_total_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
 		int label;
 
 		label = shader_bsdf_sample(kg, &sd, bsdf_u, bsdf_v, &bsdf_eval,
@@ -503,7 +444,7 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 			break;
 
 		/* modify throughput */
-		path_radiance_bsdf_bounce(&L, &throughput, &bsdf_eval, bsdf_pdf, state.bounce, label);
+		path_radiance_bsdf_bounce(L, &throughput, &bsdf_eval, bsdf_pdf, state.bounce, label);
 
 		/* set labels */
 		if(!(label & LABEL_TRANSPARENT)) {
@@ -520,46 +461,155 @@ __device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, R
 		/* setup ray */
 		ray.P = ray_offset(sd.P, (label & LABEL_TRANSMIT)? -sd.Ng: sd.Ng);
 		ray.D = bsdf_omega_in;
-
-		if(state.bounce == 0)
-			ray.t -= sd.ray_length; /* clipping works through transparent */
-		else
-			ray.t = FLT_MAX;
-
+		ray.t = FLT_MAX;
 #ifdef __RAY_DIFFERENTIALS__
 		ray.dP = sd.dP;
 		ray.dD = bsdf_domega_in;
 #endif
 	}
+}
 
-	float3 L_sum = path_radiance_sum(kg, &L);
+#endif
 
-#ifdef __CLAMP_SAMPLE__
-	path_radiance_clamp(&L, &L_sum, kernel_data.integrator.sample_clamp);
+#ifdef __SUBSURFACE__
+
+__device_inline bool kernel_path_integrate_lighting(KernelGlobals *kg, RNG *rng,
+	int sample, int num_samples,
+	ShaderData *sd, float3 *throughput,
+	float *min_ray_pdf, float *ray_pdf, PathState *state,
+	int rng_offset, PathRadiance *L, Ray *ray, float *ray_t)
+{
+#ifdef __EMISSION__
+	if(kernel_data.integrator.use_direct_light) {
+		/* sample illumination from lights to find path contribution */
+		if(sd->flag & SD_BSDF_HAS_EVAL) {
+			float light_t = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT);
+#ifdef __MULTI_CLOSURE__
+			float light_o = 0.0f;
+#else
+			float light_o = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT_F);
 #endif
+			float light_u, light_v;
+			path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT_U, &light_u, &light_v);
 
-	kernel_write_light_passes(kg, buffer, &L, sample);
+			Ray light_ray;
+			BsdfEval L_light;
+			bool is_lamp;
 
-	return make_float4(L_sum.x, L_sum.y, L_sum.z, 1.0f - L_transparent);
+#ifdef __OBJECT_MOTION__
+			light_ray.time = sd->time;
+#endif
+
+			if(direct_emission(kg, sd, -1, light_t, light_o, light_u, light_v, &light_ray, &L_light, &is_lamp, state->bounce)) {
+				/* trace shadow ray */
+				float3 shadow;
+
+				if(!shadow_blocked(kg, state, &light_ray, &shadow)) {
+					/* accumulate */
+					path_radiance_accum_light(L, *throughput, &L_light, shadow, 1.0f, state->bounce, is_lamp);
+				}
+			}
+		}
+	}
+#endif
+
+	/* no BSDF? we can stop here */
+	if(!(sd->flag & SD_BSDF))
+		return false;
+
+	/* sample BSDF */
+	float bsdf_pdf;
+	BsdfEval bsdf_eval;
+	float3 bsdf_omega_in;
+	differential3 bsdf_domega_in;
+	float bsdf_u, bsdf_v;
+	path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
+	int label;
+
+	label = shader_bsdf_sample(kg, sd, bsdf_u, bsdf_v, &bsdf_eval,
+		&bsdf_omega_in, &bsdf_domega_in, &bsdf_pdf);
+
+	if(bsdf_pdf == 0.0f || bsdf_eval_is_zero(&bsdf_eval))
+		return false;
+
+	/* modify throughput */
+	path_radiance_bsdf_bounce(L, throughput, &bsdf_eval, bsdf_pdf, state->bounce, label);
+
+	/* set labels */
+	if(!(label & LABEL_TRANSPARENT)) {
+		*ray_pdf = bsdf_pdf;
+#ifdef __LAMP_MIS__
+		*ray_t = 0.0f;
+#endif
+		*min_ray_pdf = fminf(bsdf_pdf, *min_ray_pdf);
+	}
+
+	/* update path state */
+	path_state_next(kg, state, label);
+
+	/* setup ray */
+	ray->P = ray_offset(sd->P, (label & LABEL_TRANSMIT)? -sd->Ng: sd->Ng);
+	ray->D = bsdf_omega_in;
+
+	if(state->bounce == 0)
+		ray->t -= sd->ray_length; /* clipping works through transparent */
+	else
+		ray->t = FLT_MAX;
+
+#ifdef __RAY_DIFFERENTIALS__
+	ray->dP = sd->dP;
+	ray->dD = bsdf_domega_in;
+#endif
+	
+	return true;
 }
 
-#ifdef __BRANCHED_PATH__
+#endif
 
-__device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray ray, __global float *buffer,
-	float3 throughput, int num_samples, int num_total_samples,
-	float min_ray_pdf, float ray_pdf, PathState state, int rng_offset, PathRadiance *L)
+__device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample, Ray ray, __global float *buffer)
 {
-#ifdef __LAMP_MIS__
+	/* initialize */
+	PathRadiance L;
+	float3 throughput = make_float3(1.0f, 1.0f, 1.0f);
+	float L_transparent = 0.0f;
+
+	path_radiance_init(&L, kernel_data.film.use_light_pass);
+
+	float min_ray_pdf = FLT_MAX;
+	float ray_pdf = 0.0f;
 	float ray_t = 0.0f;
+	PathState state;
+	int rng_offset = PRNG_BASE_NUM;
+#ifdef __CMJ__
+	int num_samples = kernel_data.integrator.aa_samples;
+#else
+	int num_samples = 0;
 #endif
 
+	path_state_init(&state);
+
 	/* path iteration */
 	for(;; rng_offset += PRNG_BOUNCE_NUM) {
 		/* intersect scene */
 		Intersection isect;
 		uint visibility = path_state_ray_visibility(kg, &state);
+
 #ifdef __HAIR__
-		bool hit = scene_intersect(kg, &ray, visibility, &isect, NULL, 0.0f, 0.0f);
+		float difl = 0.0f, extmax = 0.0f;
+		uint lcg_state = 0;
+
+		if(kernel_data.bvh.have_curves) {
+			if((kernel_data.cam.resolution == 1) && (state.flag & PATH_RAY_CAMERA)) {	
+				float3 pixdiff = ray.dD.dx + ray.dD.dy;
+				/*pixdiff = pixdiff - dot(pixdiff, ray.D)*ray.D;*/
+				difl = kernel_data.curve.minimum_width * len(pixdiff) * 0.5f;
+			}
+
+			extmax = kernel_data.curve.maximum_width;
+			lcg_state = lcg_init(*rng + rng_offset + sample*0x51633e2d);
+		}
+
+		bool hit = scene_intersect(kg, &ray, visibility, &isect, &lcg_state, difl, extmax);
 #else
 		bool hit = scene_intersect(kg, &ray, visibility, &isect);
 #endif
@@ -578,19 +628,29 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 			light_ray.dP = ray.dP;
 
 			/* intersect with lamp */
-			float light_t = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_LIGHT);
+			float light_t = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT);
 			float3 emission;
 
 			if(indirect_lamp_emission(kg, &light_ray, state.flag, ray_pdf, light_t, &emission, state.bounce))
-				path_radiance_accum_emission(L, throughput, emission, state.bounce);
+				path_radiance_accum_emission(&L, throughput, emission, state.bounce);
 		}
 #endif
 
 		if(!hit) {
+			/* eval background shader if nothing hit */
+			if(kernel_data.background.transparent && (state.flag & PATH_RAY_CAMERA)) {
+				L_transparent += average(throughput);
+
+#ifdef __PASSES__
+				if(!(kernel_data.film.pass_flag & PASS_BACKGROUND))
+#endif
+					break;
+			}
+
 #ifdef __BACKGROUND__
 			/* sample background shader */
 			float3 L_background = indirect_background(kg, &ray, state.flag, ray_pdf, state.bounce);
-			path_radiance_accum_background(L, throughput, L_background, state.bounce);
+			path_radiance_accum_background(&L, throughput, L_background, state.bounce);
 #endif
 
 			break;
@@ -599,9 +659,31 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 		/* setup shading */
 		ShaderData sd;
 		shader_setup_from_ray(kg, &sd, &isect, &ray, state.bounce);
-		float rbsdf = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_BSDF);
-		shader_eval_surface(kg, &sd, rbsdf, state.flag, SHADER_CONTEXT_INDIRECT);
-		shader_merge_closures(kg, &sd);
+		float rbsdf = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF);
+		shader_eval_surface(kg, &sd, rbsdf, state.flag, SHADER_CONTEXT_MAIN);
+
+		/* holdout */
+#ifdef __HOLDOUT__
+		if((sd.flag & (SD_HOLDOUT|SD_HOLDOUT_MASK)) && (state.flag & PATH_RAY_CAMERA)) {
+			if(kernel_data.background.transparent) {
+				float3 holdout_weight;
+				
+				if(sd.flag & SD_HOLDOUT_MASK)
+					holdout_weight = make_float3(1.0f, 1.0f, 1.0f);
+				else
+					holdout_weight = shader_holdout_eval(kg, &sd);
+
+				/* any throughput is ok, should all be identical here */
+				L_transparent += average(holdout_weight*throughput);
+			}
+
+			if(sd.flag & SD_HOLDOUT_MASK)
+				break;
+		}
+#endif
+
+		/* holdout mask objects do not write data passes */
+		kernel_write_data_passes(kg, buffer, &L, &sd, sample, state.flag, throughput);
 
 		/* blurring of bsdf after bounces, for rays that have a small likelihood
 		 * of following this particular path (diffuse, rough glossy) */
@@ -617,21 +699,22 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 #ifdef __EMISSION__
 		/* emission */
 		if(sd.flag & SD_EMISSION) {
+			/* todo: is isect.t wrong here for transparent surfaces? */
 			float3 emission = indirect_primitive_emission(kg, &sd, isect.t, state.flag, ray_pdf);
-			path_radiance_accum_emission(L, throughput, emission, state.bounce);
+			path_radiance_accum_emission(&L, throughput, emission, state.bounce);
 		}
 #endif
 
 		/* path termination. this is a strange place to put the termination, it's
 		 * mainly due to the mixed in MIS that we use. gives too many unneeded
 		 * shader evaluations, only need emission if we are going to terminate */
-		float probability = path_state_terminate_probability(kg, &state, throughput*num_samples);
+		float probability = path_state_terminate_probability(kg, &state, throughput);
 
 		if(probability == 0.0f) {
 			break;
 		}
 		else if(probability != 1.0f) {
-			float terminate = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_TERMINATE);
+			float terminate = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_TERMINATE);
 
 			if(terminate >= probability)
 				break;
@@ -639,37 +722,12 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 			throughput /= probability;
 		}
 
-#ifdef __SUBSURFACE__
-		/* bssrdf scatter to a different location on the same object, replacing
-		 * the closures with a diffuse BSDF */
-		if(sd.flag & SD_BSSRDF) {
-			float bssrdf_probability;
-			ShaderClosure *sc = subsurface_scatter_pick_closure(kg, &sd, &bssrdf_probability);
-
-			/* modify throughput for picking bssrdf or bsdf */
-			throughput *= bssrdf_probability;
-
-			/* do bssrdf scatter step if we picked a bssrdf closure */
-			if(sc) {
-				uint lcg_state = lcg_init(*rng + rng_offset + sample*0x68bc21eb);
-
-				if(old_subsurface_scatter_use(&sd)) {
-					old_subsurface_scatter_step(kg, &sd, state.flag, sc, &lcg_state, false);
-				}
-				else {
-					float bssrdf_u, bssrdf_v;
-					path_rng_2D(kg, rng, sample, num_total_samples, rng_offset + PRNG_BSDF_U, &bssrdf_u, &bssrdf_v);
-					subsurface_scatter_step(kg, &sd, state.flag, sc, &lcg_state, bssrdf_u, bssrdf_v, false);
-				}
-			}
-		}
-#endif
-
 #ifdef __AO__
 		/* ambient occlusion */
 		if(kernel_data.integrator.use_ambient_occlusion || (sd.flag & SD_AO)) {
+			/* todo: solve correlation */
 			float bsdf_u, bsdf_v;
-			path_rng_2D(kg, rng, sample, num_total_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
+			path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
 
 			float ao_factor = kernel_data.background.ao_factor;
 			float3 ao_N;
@@ -693,23 +751,84 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 				light_ray.dD = differential3_zero();
 
 				if(!shadow_blocked(kg, &state, &light_ray, &ao_shadow))
-					path_radiance_accum_ao(L, throughput, ao_bsdf, ao_shadow, state.bounce);
+					path_radiance_accum_ao(&L, throughput, ao_bsdf, ao_shadow, state.bounce);
 			}
 		}
 #endif
 
+#ifdef __SUBSURFACE__
+		/* bssrdf scatter to a different location on the same object, replacing
+		 * the closures with a diffuse BSDF */
+		if(sd.flag & SD_BSSRDF) {
+			float bssrdf_probability;
+			ShaderClosure *sc = subsurface_scatter_pick_closure(kg, &sd, &bssrdf_probability);
+
+			/* modify throughput for picking bssrdf or bsdf */
+			throughput *= bssrdf_probability;
+
+			/* do bssrdf scatter step if we picked a bssrdf closure */
+			if(sc) {
+				uint lcg_state = lcg_init(*rng + rng_offset + sample*0x68bc21eb);
+
+				if(old_subsurface_scatter_use(&sd)) {
+					old_subsurface_scatter_step(kg, &sd, state.flag, sc, &lcg_state, false);
+				}
+				else {
+					ShaderData bssrdf_sd[BSSRDF_MAX_HITS];
+					float bssrdf_u, bssrdf_v;
+					path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF_U, &bssrdf_u, &bssrdf_v);
+					int num_hits = subsurface_scatter_multi_step(kg, &sd, bssrdf_sd, state.flag, sc, &lcg_state, bssrdf_u, bssrdf_v, false);
+
+					/* compute lighting with the BSDF closure */
+					for(int hit = 0; hit < num_hits; hit++) {
+						float3 tp = throughput;
+						PathState hit_state = state;
+						Ray hit_ray = ray;
+						float hit_ray_t = ray_t;
+						float hit_ray_pdf = ray_pdf;
+						float hit_min_ray_pdf = min_ray_pdf;
+
+						hit_state.flag |= PATH_RAY_BSSRDF_ANCESTOR;
+						
+						if(kernel_path_integrate_lighting(kg, rng, sample, num_samples, &bssrdf_sd[hit],
+							&tp, &hit_min_ray_pdf, &hit_ray_pdf, &hit_state, rng_offset, &L, &hit_ray, &hit_ray_t)) {
+							kernel_path_indirect(kg, rng, sample, hit_ray, buffer,
+								tp, num_samples, num_samples,
+								hit_min_ray_pdf, hit_ray_pdf, hit_state, rng_offset+PRNG_BOUNCE_NUM, &L);
+
+							/* for render passes, sum and reset indirect light pass variables
+							 * for the next samples */
+							path_radiance_sum_indirect(&L);
+							path_radiance_reset_indirect(&L);
+						}
+					}
+
+					break;
+				}
+			}
+		}
+#endif
+
+#ifdef __SUBSURFACE__
+
+		if(!kernel_path_integrate_lighting(kg, rng, sample, num_samples, &sd,
+			&throughput, &min_ray_pdf, &ray_pdf, &state, rng_offset, &L, &ray, &ray_t))
+			break;
+
+#else
+
 #ifdef __EMISSION__
 		if(kernel_data.integrator.use_direct_light) {
 			/* sample illumination from lights to find path contribution */
 			if(sd.flag & SD_BSDF_HAS_EVAL) {
-				float light_t = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_LIGHT);
+				float light_t = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT);
 #ifdef __MULTI_CLOSURE__
 				float light_o = 0.0f;
 #else
-				float light_o = path_rng_1D(kg, rng, sample, num_total_samples, rng_offset + PRNG_LIGHT_F);
+				float light_o = path_rng_1D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT_F);
 #endif
 				float light_u, light_v;
-				path_rng_2D(kg, rng, sample, num_total_samples, rng_offset + PRNG_LIGHT_U, &light_u, &light_v);
+				path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_LIGHT_U, &light_u, &light_v);
 
 				Ray light_ray;
 				BsdfEval L_light;
@@ -719,14 +838,13 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 				light_ray.time = sd.time;
 #endif
 
-				/* sample random light */
 				if(direct_emission(kg, &sd, -1, light_t, light_o, light_u, light_v, &light_ray, &L_light, &is_lamp, state.bounce)) {
 					/* trace shadow ray */
 					float3 shadow;
 
 					if(!shadow_blocked(kg, &state, &light_ray, &shadow)) {
 						/* accumulate */
-						path_radiance_accum_light(L, throughput, &L_light, shadow, 1.0f, state.bounce, is_lamp);
+						path_radiance_accum_light(&L, throughput, &L_light, shadow, 1.0f, state.bounce, is_lamp);
 					}
 				}
 			}
@@ -743,7 +861,7 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 		float3 bsdf_omega_in;
 		differential3 bsdf_domega_in;
 		float bsdf_u, bsdf_v;
-		path_rng_2D(kg, rng, sample, num_total_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
+		path_rng_2D(kg, rng, sample, num_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
 		int label;
 
 		label = shader_bsdf_sample(kg, &sd, bsdf_u, bsdf_v, &bsdf_eval,
@@ -753,7 +871,7 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 			break;
 
 		/* modify throughput */
-		path_radiance_bsdf_bounce(L, &throughput, &bsdf_eval, bsdf_pdf, state.bounce, label);
+		path_radiance_bsdf_bounce(&L, &throughput, &bsdf_eval, bsdf_pdf, state.bounce, label);
 
 		/* set labels */
 		if(!(label & LABEL_TRANSPARENT)) {
@@ -770,59 +888,39 @@ __device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ray
 		/* setup ray */
 		ray.P = ray_offset(sd.P, (label & LABEL_TRANSMIT)? -sd.Ng: sd.Ng);
 		ray.D = bsdf_omega_in;
-		ray.t = FLT_MAX;
+
+		if(state.bounce == 0)
+			ray.t -= sd.ray_length; /* clipping works through transparent */
+		else
+			ray.t = FLT_MAX;
+
 #ifdef __RAY_DIFFERENTIALS__
 		ray.dP = sd.dP;
 		ray.dD = bsdf_domega_in;
 #endif
+
+#endif
 	}
+
+	float3 L_sum = path_radiance_sum(kg, &L);
+
+#ifdef __CLAMP_SAMPLE__
+	path_radiance_clamp(&L, &L_sum, kernel_data.integrator.sample_clamp);
+#endif
+
+	kernel_write_light_passes(kg, buffer, &L, sample);
+
+	return make_float4(L_sum.x, L_sum.y, L_sum.z, 1.0f - L_transparent);
 }
 
+#ifdef __BRANCHED_PATH__
+
 __device_noinline void kernel_branched_path_integrate_lighting(KernelGlobals *kg, RNG *rng,
 	int sample, int aa_samples,
 	ShaderData *sd, float3 throughput, float num_samples_adjust,
 	float min_ray_pdf, float ray_pdf, PathState state,
 	int rng_offset, PathRadiance *L, __global float *buffer)
 {
-#ifdef __AO__
-	/* ambient occlusion */
-	if(kernel_data.integrator.use_ambient_occlusion || (sd->flag & SD_AO)) {
-		int num_samples = ceil_to_int(kernel_data.integrator.ao_samples*num_samples_adjust);
-		float num_samples_inv = num_samples_adjust/num_samples;
-		float ao_factor = kernel_data.background.ao_factor;
-		float3 ao_N;
-		float3 ao_bsdf = shader_bsdf_ao(kg, sd, ao_factor, &ao_N);
-
-		for(int j = 0; j < num_samples; j++) {
-			float bsdf_u, bsdf_v;
-			path_rng_2D(kg, rng, sample*num_samples + j, aa_samples*num_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
-
-			float3 ao_D;
-			float ao_pdf;
-
-			sample_cos_hemisphere(ao_N, bsdf_u, bsdf_v, &ao_D, &ao_pdf);
-
-			if(dot(sd->Ng, ao_D) > 0.0f && ao_pdf != 0.0f) {
-				Ray light_ray;
-				float3 ao_shadow;
-
-				light_ray.P = ray_offset(sd->P, sd->Ng);
-				light_ray.D = ao_D;
-				light_ray.t = kernel_data.background.ao_distance;
-#ifdef __OBJECT_MOTION__
-				light_ray.time = sd->time;
-#endif
-				light_ray.dP = sd->dP;
-				light_ray.dD = differential3_zero();
-
-				if(!shadow_blocked(kg, &state, &light_ray, &ao_shadow))
-					path_radiance_accum_ao(L, throughput*num_samples_inv, ao_bsdf, ao_shadow, state.bounce);
-			}
-		}
-	}
-#endif
-
-
 #ifdef __EMISSION__
 	/* sample illumination from lights to find path contribution */
 	if(sd->flag & SD_BSDF_HAS_EVAL) {
@@ -901,8 +999,10 @@ __device_noinline void kernel_branched_path_integrate_lighting(KernelGlobals *kg
 
 		int num_samples;
 
-		if(CLOSURE_IS_BSDF_DIFFUSE(sc->type) || CLOSURE_IS_BSDF_BSSRDF(sc->type))
+		if(CLOSURE_IS_BSDF_DIFFUSE(sc->type))
 			num_samples = kernel_data.integrator.diffuse_samples;
+		else if(CLOSURE_IS_BSDF_BSSRDF(sc->type))
+			num_samples = 1;
 		else if(CLOSURE_IS_BSDF_GLOSSY(sc->type))
 			num_samples = kernel_data.integrator.glossy_samples;
 		else
@@ -1086,6 +1186,44 @@ __device float4 kernel_branched_path_integrate(KernelGlobals *kg, RNG *rng, int
 			}
 		}
 
+#ifdef __AO__
+		/* ambient occlusion */
+		if(kernel_data.integrator.use_ambient_occlusion || (sd.flag & SD_AO)) {
+			int num_samples = kernel_data.integrator.ao_samples;
+			float num_samples_inv = 1.0f/num_samples;
+			float ao_factor = kernel_data.background.ao_factor;
+			float3 ao_N;
+			float3 ao_bsdf = shader_bsdf_ao(kg, &sd, ao_factor, &ao_N);
+
+			for(int j = 0; j < num_samples; j++) {
+				float bsdf_u, bsdf_v;
+				path_rng_2D(kg, rng, sample*num_samples + j, aa_samples*num_samples, rng_offset + PRNG_BSDF_U, &bsdf_u, &bsdf_v);
+
+				float3 ao_D;
+				float ao_pdf;
+
+				sample_cos_hemisphere(ao_N, bsdf_u, bsdf_v, &ao_D, &ao_pdf);
+
+				if(dot(sd.Ng, ao_D) > 0.0f && ao_pdf != 0.0f) {
+					Ray light_ray;
+					float3 ao_shadow;
+
+					light_ray.P = ray_offset(sd.P, sd.Ng);
+					light_ray.D = ao_D;
+					light_ray.t = kernel_data.background.ao_distance;
+#ifdef __OBJECT_MOTION__
+					light_ray.time = sd.time;
+#endif
+					light_ray.dP = sd.dP;
+					light_ray.dD = differential3_zero();
+
+					if(!shadow_blocked(kg, &state, &light_ray, &ao_shadow))
+						path_radiance_accum_ao(&L, throughput*num_samples_inv, ao_bsdf, ao_shadow, state.bounce);
+				}
+			}
+		}
+#endif
+
 #ifdef __SUBSURFACE__
 		/* bssrdf scatter to a different location on the same object */
 		if(sd.flag & SD_BSSRDF) {
@@ -1101,6 +1239,8 @@ __device float4 kernel_branched_path_integrate(KernelGlobals *kg, RNG *rng, int
 				float num_samples_inv = 1.0f/num_samples;
 				RNG bssrdf_rng = cmj_hash(*rng, i);
 
+				state.flag |= PATH_RAY_BSSRDF_ANCESTOR;
+
 				/* do subsurface scatter step with copy of shader data, this will
 				 * replace the BSSRDF with a diffuse BSDF closure */
 				for(int j = 0; j < num_samples; j++) {
@@ -1128,6 +1268,8 @@ __device float4 kernel_branched_path_integrate(KernelGlobals *kg, RNG *rng, int
 								ray_pdf, ray_pdf, state, rng_offset, &L, buffer);
 					}
 				}
+
+				state.flag &= ~PATH_RAY_BSSRDF_ANCESTOR;
 			}
 		}
 #endif
diff --git a/intern/cycles/kernel/kernel_subsurface.h b/intern/cycles/kernel/kernel_subsurface.h
index 7127775669b..5eb641f5af4 100644
--- a/intern/cycles/kernel/kernel_subsurface.h
+++ b/intern/cycles/kernel/kernel_subsurface.h
@@ -221,20 +221,18 @@ __device int subsurface_scatter_multi_step(KernelGlobals *kg, ShaderData *sd, Sh
 	disk_N = sd->Ng;
 	make_orthonormals(disk_N, &disk_T, &disk_B);
 
-	if(disk_u < 0.5f) {
+	if(sd->randb_closure < 0.5f) {
 		pick_pdf_N = 0.5f;
 		pick_pdf_T = 0.25f;
 		pick_pdf_B = 0.25f;
-		disk_u *= 2.0f;
 	}
-	else if(disk_u < 0.75f) {
+	else if(sd->randb_closure < 0.75f) {
 		float3 tmp = disk_N;
 		disk_N = disk_T;
 		disk_T = tmp;
 		pick_pdf_N = 0.25f;
 		pick_pdf_T = 0.5f;
 		pick_pdf_B = 0.25f;
-		disk_u = (disk_u - 0.5f)*4.0f;
 	}
 	else {
 		float3 tmp = disk_N;
@@ -243,7 +241,6 @@ __device int subsurface_scatter_multi_step(KernelGlobals *kg, ShaderData *sd, Sh
 		pick_pdf_N = 0.25f;
 		pick_pdf_T = 0.25f;
 		pick_pdf_B = 0.5f;
-		disk_u = (disk_u - 0.75f)*4.0f;
 	}
 
 	/* sample point on disk */
@@ -323,20 +320,18 @@ __device void subsurface_scatter_step(KernelGlobals *kg, ShaderData *sd,
 	disk_N = sd->Ng;
 	make_orthonormals(disk_N, &disk_T, &disk_B);
 
-	if(disk_u < 0.5f) {
+	if(sd->randb_closure < 0.5f) {
 		pick_pdf_N = 0.5f;
 		pick_pdf_T = 0.25f;
 		pick_pdf_B = 0.25f;
-		disk_u *= 2.0f;
 	}
-	else if(disk_u < 0.75f) {
+	else if(sd->randb_closure < 0.75f) {
 		float3 tmp = disk_N;
 		disk_N = disk_T;
 		disk_T = tmp;
 		pick_pdf_N = 0.25f;
 		pick_pdf_T = 0.5f;
 		pick_pdf_B = 0.25f;
-		disk_u = (disk_u - 0.5f)*4.0f;
 	}
 	else {
 		float3 tmp = disk_N;
@@ -345,7 +340,6 @@ __device void subsurface_scatter_step(KernelGlobals *kg, ShaderData *sd,
 		pick_pdf_N = 0.25f;
 		pick_pdf_T = 0.25f;
 		pick_pdf_B = 0.5f;
-		disk_u = (disk_u - 0.75f)*4.0f;
 	}
 
 	/* sample point on disk */
diff --git a/intern/cycles/kernel/kernel_types.h b/intern/cycles/kernel/kernel_types.h
index caf68e913e8..f64277564c4 100644
--- a/intern/cycles/kernel/kernel_types.h
+++ b/intern/cycles/kernel/kernel_types.h
@@ -220,6 +220,7 @@ enum PathRayFlag {
 	PATH_RAY_MIS_SKIP = 1024,
 	PATH_RAY_DIFFUSE_ANCESTOR = 2048,
 	PATH_RAY_GLOSSY_ANCESTOR = 4096,
+	PATH_RAY_BSSRDF_ANCESTOR = 8192,
 
 	/* this gives collisions with localview bits
 	 * see: blender_util.h, grr - Campbell */
diff --git a/intern/cycles/render/osl.cpp b/intern/cycles/render/osl.cpp
index 038da269d91..d0ed95cca34 100644
--- a/intern/cycles/render/osl.cpp
+++ b/intern/cycles/render/osl.cpp
@@ -209,6 +209,7 @@ void OSLShaderManager::shading_system_init()
 			"__unused__",
 			"diffuse_ancestor", /* PATH_RAY_DIFFUSE_ANCESTOR */
 			"glossy_ancestor",  /* PATH_RAY_GLOSSY_ANCESTOR */
+			"bssrdf_ancestor",  /* PATH_RAY_BSSRDF_ANCESTOR */
 		};
 
 		const int nraytypes = sizeof(raytypes)/sizeof(raytypes[0]);