Merge branch 'master' into blender2.8

28 files changed, 417 insertions, 78 deletions

diff --git a/intern/cycles/blender/addon/version_update.py b/intern/cycles/blender/addon/version_update.py
index 12f1cec4ac9..f1dc5a6e15e 100644
--- a/intern/cycles/blender/addon/version_update.py
+++ b/intern/cycles/blender/addon/version_update.py
@@ -130,10 +130,13 @@ def displacement_nodes_insert():
         if check_is_new_shading_material(material):
             displacement_node_insert(material, material.node_tree, traversed)
 
-def displacement_node_space(node):
+def displacement_principled_nodes(node):
     if node.bl_idname == 'ShaderNodeDisplacement':
         if node.space != 'WORLD':
             node.space = 'OBJECT'
+    if node.bl_idname == 'ShaderNodeBsdfPrincipled':
+        if node.subsurface_method != 'RANDOM_WALK':
+            node.subsurface_method = 'BURLEY'
 
 
 def mapping_node_order_flip(node):
@@ -373,4 +376,4 @@ def do_versions(self):
             if not cmat.is_property_set("displacement_method"):
                 cmat.displacement_method = 'BUMP'
 
-        foreach_cycles_node(displacement_node_space)
+        foreach_cycles_node(displacement_principled_nodes)
diff --git a/intern/cycles/blender/blender_shader.cpp b/intern/cycles/blender/blender_shader.cpp
index d4c34d52011..131e76d106e 100644
--- a/intern/cycles/blender/blender_shader.cpp
+++ b/intern/cycles/blender/blender_shader.cpp
@@ -427,6 +427,9 @@ static ShaderNode *add_node(Scene *scene,
 			case BL::ShaderNodeSubsurfaceScattering::falloff_BURLEY:
 				subsurface->falloff = CLOSURE_BSSRDF_BURLEY_ID;
 				break;
+			case BL::ShaderNodeSubsurfaceScattering::falloff_RANDOM_WALK:
+				subsurface->falloff = CLOSURE_BSSRDF_RANDOM_WALK_ID;
+				break;
 		}
 
 		node = subsurface;
@@ -526,6 +529,14 @@ static ShaderNode *add_node(Scene *scene,
 				principled->distribution = CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID;
 				break;
 		}
+		switch (b_principled_node.subsurface_method()) {
+			case BL::ShaderNodeBsdfPrincipled::subsurface_method_BURLEY:
+				principled->subsurface_method = CLOSURE_BSSRDF_PRINCIPLED_ID;
+				break;
+			case BL::ShaderNodeBsdfPrincipled::subsurface_method_RANDOM_WALK:
+				principled->subsurface_method = CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID;
+				break;
+		}
 		node = principled;
 	}
 	else if(b_node.is_a(&RNA_ShaderNodeBsdfTranslucent)) {
diff --git a/intern/cycles/kernel/closure/bsdf.h b/intern/cycles/kernel/closure/bsdf.h
index 6f0bdb3fa38..e3beff6675a 100644
--- a/intern/cycles/kernel/closure/bsdf.h
+++ b/intern/cycles/kernel/closure/bsdf.h
@@ -30,9 +30,7 @@
 #include "kernel/closure/bsdf_principled_diffuse.h"
 #include "kernel/closure/bsdf_principled_sheen.h"
 #include "kernel/closure/bssrdf.h"
-#ifdef __VOLUME__
-#  include "kernel/closure/volume.h"
-#endif
+#include "kernel/closure/volume.h"
 
 CCL_NAMESPACE_BEGIN
 
diff --git a/intern/cycles/kernel/closure/bssrdf.h b/intern/cycles/kernel/closure/bssrdf.h
index c8f505e8418..8578767b07e 100644
--- a/intern/cycles/kernel/closure/bssrdf.h
+++ b/intern/cycles/kernel/closure/bssrdf.h
@@ -373,7 +373,9 @@ ccl_device int bssrdf_setup(ShaderData *sd, Bssrdf *bssrdf, ClosureType type)
 	if(bssrdf_channels < 3) {
 		/* Add diffuse BSDF if any radius too small. */
 #ifdef __PRINCIPLED__
-		if(type == CLOSURE_BSSRDF_PRINCIPLED_ID) {
+		if(type == CLOSURE_BSSRDF_PRINCIPLED_ID ||
+		   type == CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID)
+		{
 			float roughness = bssrdf->roughness;
 			float3 N = bssrdf->N;
 
@@ -408,7 +410,9 @@ ccl_device int bssrdf_setup(ShaderData *sd, Bssrdf *bssrdf, ClosureType type)
 		bssrdf->sharpness = saturate(bssrdf->sharpness);
 
 		if(type == CLOSURE_BSSRDF_BURLEY_ID ||
-		   type == CLOSURE_BSSRDF_PRINCIPLED_ID)
+		   type == CLOSURE_BSSRDF_PRINCIPLED_ID ||
+		   type == CLOSURE_BSSRDF_RANDOM_WALK_ID ||
+		   type == CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID)
 		{
 			bssrdf_burley_setup(bssrdf);
 		}
diff --git a/intern/cycles/kernel/closure/volume.h b/intern/cycles/kernel/closure/volume.h
index 4bb5e680723..da791e9aa73 100644
--- a/intern/cycles/kernel/closure/volume.h
+++ b/intern/cycles/kernel/closure/volume.h
@@ -83,35 +83,45 @@ ccl_device float3 volume_henyey_greenstein_eval_phase(const ShaderClosure *sc, c
 	return make_float3(*pdf, *pdf, *pdf);
 }
 
-ccl_device int volume_henyey_greenstein_sample(const ShaderClosure *sc, float3 I, float3 dIdx, float3 dIdy, float randu, float randv,
-	float3 *eval, float3 *omega_in, float3 *domega_in_dx, float3 *domega_in_dy, float *pdf)
+ccl_device float3 henyey_greenstrein_sample(float3 D, float g, float randu, float randv, float *pdf)
 {
-	const HenyeyGreensteinVolume *volume = (const HenyeyGreensteinVolume*)sc;
-	float g = volume->g;
-	float cos_phi, sin_phi, cos_theta;
-
 	/* match pdf for small g */
-	if(fabsf(g) < 1e-3f) {
+	float cos_theta;
+	bool isotropic = fabsf(g) < 1e-3f;
+
+	if(isotropic) {
 		cos_theta = (1.0f - 2.0f * randu);
-		*pdf = M_1_PI_F * 0.25f;
+		if(pdf) {
+			*pdf = M_1_PI_F * 0.25f;
+		}
 	}
 	else {
 		float k = (1.0f - g * g) / (1.0f - g + 2.0f * g * randu);
 		cos_theta = (1.0f + g * g - k * k) / (2.0f * g);
-		*pdf = single_peaked_henyey_greenstein(cos_theta, g);
+		if(pdf) {
+			*pdf = single_peaked_henyey_greenstein(cos_theta, g);
+		}
 	}
 
 	float sin_theta = safe_sqrtf(1.0f - cos_theta * cos_theta);
-
 	float phi = M_2PI_F * randv;
-	cos_phi = cosf(phi);
-	sin_phi = sinf(phi);
+	float3 dir = make_float3(sin_theta * cosf(phi), sin_theta * sinf(phi), cos_theta);
 
-	/* note that I points towards the viewer and so is used negated */
 	float3 T, B;
-	make_orthonormals(-I, &T, &B);
-	*omega_in = sin_theta * cos_phi * T + sin_theta * sin_phi * B + cos_theta * (-I);
+	make_orthonormals(D, &T, &B);
+	dir = dir.x * T + dir.y * B + dir.z * D;
+
+	return dir;
+}
+
+ccl_device int volume_henyey_greenstein_sample(const ShaderClosure *sc, float3 I, float3 dIdx, float3 dIdy, float randu, float randv,
+	float3 *eval, float3 *omega_in, float3 *domega_in_dx, float3 *domega_in_dy, float *pdf)
+{
+	const HenyeyGreensteinVolume *volume = (const HenyeyGreensteinVolume*)sc;
+	float g = volume->g;
 
+	/* note that I points towards the viewer and so is used negated */
+	*omega_in = henyey_greenstrein_sample(-I, g, randu, randv, pdf);
 	*eval = make_float3(*pdf, *pdf, *pdf); /* perfect importance sampling */
 
 #ifdef __RAY_DIFFERENTIALS__
diff --git a/intern/cycles/kernel/geom/geom_motion_triangle_intersect.h b/intern/cycles/kernel/geom/geom_motion_triangle_intersect.h
index 4a19f2ba031..542843edc84 100644
--- a/intern/cycles/kernel/geom/geom_motion_triangle_intersect.h
+++ b/intern/cycles/kernel/geom/geom_motion_triangle_intersect.h
@@ -248,23 +248,30 @@ ccl_device_inline void motion_triangle_intersect_local(
 	motion_triangle_vertices(kg, local_object, prim, time, verts);
 	/* Ray-triangle intersection, unoptimized. */
 	float t, u, v;
-	if(ray_triangle_intersect(P,
-	                          dir,
-	                          tmax,
+	if(!ray_triangle_intersect(P,
+	                           dir,
+	                           tmax,
 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
-	                          (ssef*)verts,
+	                           (ssef*)verts,
 #else
-	                          verts[0], verts[1], verts[2],
+	                           verts[0], verts[1], verts[2],
 #endif
-	                          &u, &v, &t))
+	                           &u, &v, &t))
 	{
+		return;
+	}
+
+	int hit;
+	if(lcg_state) {
+		/* Record up to max_hits intersections. */
 		for(int i = min(max_hits, local_isect->num_hits) - 1; i >= 0; --i) {
 			if(local_isect->hits[i].t == t) {
 				return;
 			}
 		}
+
 		local_isect->num_hits++;
-		int hit;
+
 		if(local_isect->num_hits <= max_hits) {
 			hit = local_isect->num_hits - 1;
 		}
@@ -277,18 +284,29 @@ ccl_device_inline void motion_triangle_intersect_local(
 			if(hit >= max_hits)
 				return;
 		}
-		/* Record intersection. */
-		Intersection *isect = &local_isect->hits[hit];
-		isect->t = t;
-		isect->u = u;
-		isect->v = v;
-		isect->prim = prim_addr;
-		isect->object = object;
-		isect->type = PRIMITIVE_MOTION_TRIANGLE;
-		/* Record geometric normal. */
-		local_isect->Ng[hit] = normalize(cross(verts[1] - verts[0],
-		                                    verts[2] - verts[0]));
 	}
+	else {
+		/* Record closest intersection only. */
+		if(local_isect->num_hits && t > local_isect->hits[0].t) {
+			return;
+		}
+
+		hit = 0;
+		local_isect->num_hits = 1;
+	}
+
+	/* Record intersection. */
+	Intersection *isect = &local_isect->hits[hit];
+	isect->t = t;
+	isect->u = u;
+	isect->v = v;
+	isect->prim = prim_addr;
+	isect->object = object;
+	isect->type = PRIMITIVE_MOTION_TRIANGLE;
+
+	/* Record geometric normal. */
+	local_isect->Ng[hit] = normalize(cross(verts[1] - verts[0],
+	                                       verts[2] - verts[0]));
 }
 #endif  /* __BVH_LOCAL__ */
 
diff --git a/intern/cycles/kernel/geom/geom_triangle_intersect.h b/intern/cycles/kernel/geom/geom_triangle_intersect.h
index 7daa378f81e..a3b23115ae4 100644
--- a/intern/cycles/kernel/geom/geom_triangle_intersect.h
+++ b/intern/cycles/kernel/geom/geom_triangle_intersect.h
@@ -118,28 +118,40 @@ ccl_device_inline void triangle_intersect_local(
 		return;
 	}
 
-	for(int i = min(max_hits, local_isect->num_hits) - 1; i >= 0; --i) {
-		if(local_isect->hits[i].t == t) {
-			return;
+	int hit;
+	if(lcg_state) {
+		/* Record up to max_hits intersections. */
+		for(int i = min(max_hits, local_isect->num_hits) - 1; i >= 0; --i) {
+			if(local_isect->hits[i].t == t) {
+				return;
+			}
 		}
-	}
 
-	local_isect->num_hits++;
-	int hit;
+		local_isect->num_hits++;
+
+		if(local_isect->num_hits <= max_hits) {
+			hit = local_isect->num_hits - 1;
+		}
+		else {
+			/* reservoir sampling: if we are at the maximum number of
+			 * hits, randomly replace element or skip it */
+			hit = lcg_step_uint(lcg_state) % local_isect->num_hits;
 
-	if(local_isect->num_hits <= max_hits) {
-		hit = local_isect->num_hits - 1;
+			if(hit >= max_hits)
+				return;
+		}
 	}
 	else {
-		/* reservoir sampling: if we are at the maximum number of
-		 * hits, randomly replace element or skip it */
-		hit = lcg_step_uint(lcg_state) % local_isect->num_hits;
-
-		if(hit >= max_hits)
+		/* Record closest intersection only. */
+		if(local_isect->num_hits && t > local_isect->hits[0].t) {
 			return;
+		}
+
+		hit = 0;
+		local_isect->num_hits = 1;
 	}
 
-	/* record intersection */
+	/* Record intersection. */
 	Intersection *isect = &local_isect->hits[hit];
 	isect->prim = prim_addr;
 	isect->object = object;
diff --git a/intern/cycles/kernel/kernel_path.h b/intern/cycles/kernel/kernel_path.h
index afca4575331..dbbb80ca37f 100644
--- a/intern/cycles/kernel/kernel_path.h
+++ b/intern/cycles/kernel/kernel_path.h
@@ -32,12 +32,12 @@
 #include "kernel/kernel_light.h"
 #include "kernel/kernel_passes.h"
 
-#ifdef __SUBSURFACE__
-#  include "kernel/kernel_subsurface.h"
+#if defined(__VOLUME__) || defined(__SUBSURFACE__)
+#  include "kernel/kernel_volume.h"
 #endif
 
-#ifdef __VOLUME__
-#  include "kernel/kernel_volume.h"
+#ifdef __SUBSURFACE__
+#  include "kernel/kernel_subsurface.h"
 #endif
 
 #include "kernel/kernel_path_state.h"
diff --git a/intern/cycles/kernel/kernel_path_branched.h b/intern/cycles/kernel/kernel_path_branched.h
index 5f917d509ec..6fb55bda08d 100644
--- a/intern/cycles/kernel/kernel_path_branched.h
+++ b/intern/cycles/kernel/kernel_path_branched.h
@@ -350,6 +350,7 @@ ccl_device void kernel_branched_path_subsurface_scatter(KernelGlobals *kg,
 			int num_hits = subsurface_scatter_multi_intersect(kg,
 			                                                  &ss_isect,
 			                                                  sd,
+			                                                  &hit_state,
 			                                                  sc,
 			                                                  &lcg_state,
 			                                                  bssrdf_u, bssrdf_v,
diff --git a/intern/cycles/kernel/kernel_path_subsurface.h b/intern/cycles/kernel/kernel_path_subsurface.h
index a48bde6443e..71aea9e3b27 100644
--- a/intern/cycles/kernel/kernel_path_subsurface.h
+++ b/intern/cycles/kernel/kernel_path_subsurface.h
@@ -51,6 +51,7 @@ bool kernel_path_subsurface_scatter(
 		int num_hits = subsurface_scatter_multi_intersect(kg,
 		                                                  &ss_isect,
 		                                                  sd,
+		                                                  state,
 		                                                  sc,
 		                                                  &lcg_state,
 		                                                  bssrdf_u, bssrdf_v,
diff --git a/intern/cycles/kernel/kernel_subsurface.h b/intern/cycles/kernel/kernel_subsurface.h
index f4759b26191..80dda31c61e 100644
--- a/intern/cycles/kernel/kernel_subsurface.h
+++ b/intern/cycles/kernel/kernel_subsurface.h
@@ -20,7 +20,6 @@ CCL_NAMESPACE_BEGIN
  *
  * BSSRDF Importance Sampling, SIGGRAPH 2013
  * http://library.imageworks.com/pdfs/imageworks-library-BSSRDF-sampling.pdf
- *
  */
 
 ccl_device_inline float3 subsurface_scatter_eval(ShaderData *sd,
@@ -41,7 +40,7 @@ ccl_device_inline float3 subsurface_scatter_eval(ShaderData *sd,
 		for(int i = 0; i < sd->num_closure; i++) {
 			sc = &sd->closure[i];
 
-			if(CLOSURE_IS_BSSRDF(sc->type)) {
+			if(CLOSURE_IS_DISK_BSSRDF(sc->type)) {
 				sample_weight_sum += sc->sample_weight;
 			}
 		}
@@ -52,7 +51,7 @@ ccl_device_inline float3 subsurface_scatter_eval(ShaderData *sd,
 	for(int i = 0; i < sd->num_closure; i++) {
 		sc = &sd->closure[i];
 		
-		if(CLOSURE_IS_BSSRDF(sc->type)) {
+		if(CLOSURE_IS_DISK_BSSRDF(sc->type)) {
 			/* in case of branched path integrate we sample all bssrdf's once,
 			 * for path trace we pick one, so adjust pdf for that */
 			float sample_weight = (all)? 1.0f: sc->sample_weight * sample_weight_inv;
@@ -79,7 +78,9 @@ ccl_device void subsurface_scatter_setup_diffuse_bsdf(KernelGlobals *kg, ShaderD
 	if(hit) {
 		Bssrdf *bssrdf = (Bssrdf *)sc;
 #ifdef __PRINCIPLED__
-		if(bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID) {
+		if(bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID ||
+		   bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID)
+		{
 			PrincipledDiffuseBsdf *bsdf = (PrincipledDiffuseBsdf*)bsdf_alloc(sd, sizeof(PrincipledDiffuseBsdf), weight);
 
 			if(bsdf) {
@@ -166,7 +167,7 @@ ccl_device void subsurface_color_bump_blur(KernelGlobals *kg,
 /* Subsurface scattering step, from a point on the surface to other
  * nearby points on the same object.
  */
-ccl_device_inline int subsurface_scatter_multi_intersect(
+ccl_device_inline int subsurface_scatter_disk(
         KernelGlobals *kg,
         LocalIntersection *ss_isect,
         ShaderData *sd,
@@ -433,5 +434,202 @@ ccl_device void subsurface_scatter_step(KernelGlobals *kg, ShaderData *sd, ccl_a
 	subsurface_scatter_setup_diffuse_bsdf(kg, sd, sc, eval, (ss_isect.num_hits > 0), N);
 }
 
+/* Random walk subsurface scattering.
+ *
+ * "Practical and Controllable Subsurface Scattering for Production Path
+ *  Tracing". Matt Jen-Yuan Chiang, Peter Kutz, Brent Burley. SIGGRAPH 2016. */
+
+ccl_device void subsurface_random_walk_remap(
+        const float A,
+        const float d,
+        float *sigma_t,
+        float *sigma_s)
+{
+	/* Compute attenuation and scattering coefficients from albedo. */
+	const float a = 1.0f - expf(A * (-5.09406f + A * (2.61188f - A * 4.31805f)));
+	const float s = 1.9f - A + 3.5f * sqr(A - 0.8f);
+
+	*sigma_t = 1.0f / fmaxf(d * s, 1e-16f);
+	*sigma_s = *sigma_t * a;
+}
+
+ccl_device void subsurface_random_walk_coefficients(
+        const ShaderClosure *sc,
+        float3 *sigma_t,
+        float3 *sigma_s,
+        float3 *weight)
+{
+	const Bssrdf *bssrdf = (const Bssrdf*)sc;
+	const float3 A = bssrdf->albedo;
+	const float3 d = bssrdf->radius;
+	float sigma_t_x, sigma_t_y, sigma_t_z;
+	float sigma_s_x, sigma_s_y, sigma_s_z;
+
+	subsurface_random_walk_remap(A.x, d.x, &sigma_t_x, &sigma_s_x);
+	subsurface_random_walk_remap(A.y, d.y, &sigma_t_y, &sigma_s_y);
+	subsurface_random_walk_remap(A.z, d.z, &sigma_t_z, &sigma_s_z);
+
+	*sigma_t = make_float3(sigma_t_x, sigma_t_y, sigma_t_z);
+	*sigma_s = make_float3(sigma_s_x, sigma_s_y, sigma_s_z);
+
+	/* Closure mixing and Fresnel weights separate from albedo. */
+	*weight = safe_divide_color(bssrdf->weight, A);
+}
+
+ccl_device_noinline bool subsurface_random_walk(
+        KernelGlobals *kg,
+        LocalIntersection *ss_isect,
+        ShaderData *sd,
+        ccl_addr_space PathState *state,
+        const ShaderClosure *sc,
+        const float bssrdf_u,
+        const float bssrdf_v)
+{
+	/* Sample diffuse surface scatter into the object. */
+	float3 D;
+	float pdf;
+	sample_cos_hemisphere(-sd->N, bssrdf_u, bssrdf_v, &D, &pdf);
+	if(dot(-sd->Ng, D) <= 0.0f) {
+		return 0;
+	}
+
+	/* Convert subsurface to volume coefficients. */
+	float3 sigma_t, sigma_s;
+	float3 throughput = make_float3(1.0f, 1.0f, 1.0f);
+	subsurface_random_walk_coefficients(sc, &sigma_t, &sigma_s, &throughput);
+
+	/* Setup ray. */
+#ifdef __SPLIT_KERNEL__
+	Ray ray_object = ss_isect->ray;
+	Ray *ray = &ray_object;
+#else
+	Ray *ray = &ss_isect->ray;
+#endif
+	ray->P = ray_offset(sd->P, -sd->Ng);
+	ray->D = D;
+	ray->t = FLT_MAX;
+	ray->time = sd->time;
+
+	/* Modify state for RNGs, decorrelated from other paths. */
+	uint prev_rng_offset = state->rng_offset;
+	uint prev_rng_hash = state->rng_hash;
+	state->rng_hash = cmj_hash(state->rng_hash + state->rng_offset, 0xdeadbeef);
+
+	/* Random walk until we hit the surface again. */
+	bool hit = false;
+
+	for(int bounce = 0; bounce < BSSRDF_MAX_BOUNCES; bounce++) {
+		/* Advance random number offset. */
+		state->rng_offset += PRNG_BOUNCE_NUM;
+
+		if(bounce > 0) {
+			/* Sample scattering direction. */
+			const float anisotropy = 0.0f;
+			float scatter_u, scatter_v;
+			path_state_rng_2D(kg, state, PRNG_BSDF_U, &scatter_u, &scatter_v);
+			ray->D = henyey_greenstrein_sample(ray->D, anisotropy, scatter_u, scatter_v, NULL);
+		}
+
+		/* Sample color channel, use MIS with balance heuristic. */
+		float rphase = path_state_rng_1D(kg, state, PRNG_PHASE_CHANNEL);
+		float3 albedo = safe_divide_color(sigma_s, sigma_t);
+		float3 channel_pdf;
+		int channel = kernel_volume_sample_channel(albedo, throughput, rphase, &channel_pdf);
+
+		/* Distance sampling. */
+		float rdist = path_state_rng_1D(kg, state, PRNG_SCATTER_DISTANCE);
+		float sample_sigma_t = kernel_volume_channel_get(sigma_t, channel);
+		float t = -logf(1.0f - rdist)/sample_sigma_t;
+
+		ray->t = t;
+		scene_intersect_local(kg, *ray, ss_isect, sd->object, NULL, 1);
+		hit = (ss_isect->num_hits > 0);
+
+		if(hit) {
+			/* Compute world space distance to surface hit. */
+			float3 D = ray->D;
+			object_inverse_dir_transform(kg, sd, &D);
+			D = normalize(D) * ss_isect->hits[0].t;
+			object_dir_transform(kg, sd, &D);
+			t = len(D);
+		}
+
+		/* Advance to new scatter location. */
+		ray->P += t * ray->D;
+
+		/* Update throughput. */
+		float3 transmittance = volume_color_transmittance(sigma_t, t);
+		float pdf = dot(channel_pdf, (hit)? transmittance: sigma_t * transmittance);
+		throughput *= ((hit)? transmittance: sigma_s * transmittance) / pdf;
+
+		if(hit) {
+			/* If we hit the surface, we are done. */
+			break;
+		}
+
+		/* Russian roulette. */
+		float terminate = path_state_rng_1D(kg, state, PRNG_TERMINATE);
+		float probability = min(max3(fabs(throughput)), 1.0f);
+		if(terminate >= probability) {
+			break;
+		}
+		throughput /= probability;
+	}
+
+	kernel_assert(isfinite_safe(throughput.x) &&
+	              isfinite_safe(throughput.y) &&
+	              isfinite_safe(throughput.z));
+
+	state->rng_offset = prev_rng_offset;
+	state->rng_hash = prev_rng_hash;
+
+	/* Return number of hits in ss_isect. */
+	if(!hit) {
+		return 0;
+	}
+
+	/* TODO: gain back performance lost from merging with disk BSSRDF. We
+	 * only need to return on hit so this indirect ray push/pop overhead
+	 * is not actually needed, but it does keep the code simpler. */
+	ss_isect->weight[0] = throughput;
+#ifdef __SPLIT_KERNEL__
+	ss_isect->ray = *ray;
+#endif
+
+	return 1;
+}
+
+ccl_device_inline int subsurface_scatter_multi_intersect(
+        KernelGlobals *kg,
+        LocalIntersection *ss_isect,
+        ShaderData *sd,
+        ccl_addr_space PathState *state,
+        const ShaderClosure *sc,
+        uint *lcg_state,
+        float bssrdf_u,
+        float bssrdf_v,
+        bool all)
+{
+	if(CLOSURE_IS_DISK_BSSRDF(sc->type)) {
+		return subsurface_scatter_disk(kg,
+		                               ss_isect,
+		                               sd,
+		                               sc,
+		                               lcg_state,
+		                               bssrdf_u,
+		                               bssrdf_v,
+		                               all);
+	}
+	else {
+		return subsurface_random_walk(kg,
+		                              ss_isect,
+		                              sd,
+		                              state,
+		                              sc,
+		                              bssrdf_u,
+		                              bssrdf_v);
+	}
+}
+
 CCL_NAMESPACE_END
 
diff --git a/intern/cycles/kernel/kernel_types.h b/intern/cycles/kernel/kernel_types.h
index 1e1a4f34650..cd3b450932f 100644
--- a/intern/cycles/kernel/kernel_types.h
+++ b/intern/cycles/kernel/kernel_types.h
@@ -46,6 +46,7 @@ CCL_NAMESPACE_BEGIN
 
 #define BSSRDF_MIN_RADIUS			1e-8f
 #define BSSRDF_MAX_HITS				4
+#define BSSRDF_MAX_BOUNCES			256
 #define LOCAL_MAX_HITS				4
 
 #define BECKMANN_TABLE_SIZE		256
diff --git a/intern/cycles/kernel/kernel_volume.h b/intern/cycles/kernel/kernel_volume.h
index 2af4c9a5e7a..7b67a37adc5 100644
--- a/intern/cycles/kernel/kernel_volume.h
+++ b/intern/cycles/kernel/kernel_volume.h
@@ -35,6 +35,8 @@ typedef struct VolumeShaderCoefficients {
 	float3 emission;
 } VolumeShaderCoefficients;
 
+#ifdef __VOLUME__
+
 /* evaluate shader to get extinction coefficient at P */
 ccl_device_inline bool volume_shader_extinction_sample(KernelGlobals *kg,
                                                        ShaderData *sd,
@@ -92,6 +94,8 @@ ccl_device_inline bool volume_shader_sample(KernelGlobals *kg,
 	return true;
 }
 
+#endif /* __VOLUME__ */
+
 ccl_device float3 volume_color_transmittance(float3 sigma, float t)
 {
 	return make_float3(expf(-sigma.x * t), expf(-sigma.y * t), expf(-sigma.z * t));
@@ -102,6 +106,8 @@ ccl_device float kernel_volume_channel_get(float3 value, int channel)
 	return (channel == 0)? value.x: ((channel == 1)? value.y: value.z);
 }
 
+#ifdef __VOLUME__
+
 ccl_device bool volume_stack_is_heterogeneous(KernelGlobals *kg, ccl_addr_space VolumeStack *stack)
 {
 	for(int i = 0; stack[i].shader != SHADER_NONE; i++) {
@@ -239,6 +245,8 @@ ccl_device_noinline void kernel_volume_shadow(KernelGlobals *kg,
 		kernel_volume_shadow_homogeneous(kg, state, ray, shadow_sd, throughput);
 }
 
+#endif /* __VOLUME__ */
+
 /* Equi-angular sampling as in:
  * "Importance Sampling Techniques for Path Tracing in Participating Media" */
 
@@ -369,6 +377,8 @@ ccl_device int kernel_volume_sample_channel(float3 albedo, float3 throughput, fl
 	}
 }
 
+#ifdef __VOLUME__
+
 /* homogeneous volume: assume shader evaluation at the start gives
  * the volume shading coefficient for the entire line segment */
 ccl_device VolumeIntegrateResult kernel_volume_integrate_homogeneous(
@@ -1346,4 +1356,6 @@ ccl_device_inline void kernel_volume_clean_stack(KernelGlobals *kg,
 	}
 }
 
+#endif /* __VOLUME__ */
+
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/osl/osl_bssrdf.cpp b/intern/cycles/kernel/osl/osl_bssrdf.cpp
index db6426f60e5..da7368bbc61 100644
--- a/intern/cycles/kernel/osl/osl_bssrdf.cpp
+++ b/intern/cycles/kernel/osl/osl_bssrdf.cpp
@@ -52,6 +52,8 @@ static ustring u_cubic("cubic");
 static ustring u_gaussian("gaussian");
 static ustring u_burley("burley");
 static ustring u_principled("principled");
+static ustring u_random_walk("random_walk");
+static ustring u_principled_random_walk("principled_random_walk");
 
 class CBSSRDFClosure : public CClosurePrimitive {
 public:
@@ -79,6 +81,12 @@ public:
 		else if (method == u_principled) {
 			alloc(sd, path_flag, weight, CLOSURE_BSSRDF_PRINCIPLED_ID);
 		}
+		else if (method == u_random_walk) {
+			alloc(sd, path_flag, weight, CLOSURE_BSSRDF_RANDOM_WALK_ID);
+		}
+		else if (method == u_principled_random_walk) {
+			alloc(sd, path_flag, weight, CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID);
+		}
 	}
 
 	void alloc(ShaderData *sd, int path_flag, float3 weight, ClosureType type)
diff --git a/intern/cycles/kernel/shaders/node_principled_bsdf.osl b/intern/cycles/kernel/shaders/node_principled_bsdf.osl
index 0e31dcedee4..fc0a1c894da 100644
--- a/intern/cycles/kernel/shaders/node_principled_bsdf.osl
+++ b/intern/cycles/kernel/shaders/node_principled_bsdf.osl
@@ -19,6 +19,7 @@
 
 shader node_principled_bsdf(
 	string distribution = "Multiscatter GGX",
+	string subsurface_method = "burley",
 	color BaseColor = color(0.8, 0.8, 0.8),
 	float Subsurface = 0.0,
 	vector SubsurfaceRadius = vector(1.0, 1.0, 1.0),
@@ -58,8 +59,14 @@ shader node_principled_bsdf(
 	if (diffuse_weight > 1e-5) {
 		if (Subsurface > 1e-5) {
 			color mixed_ss_base_color = SubsurfaceColor * Subsurface + BaseColor * (1.0 - Subsurface);
-			BSDF = mixed_ss_base_color * bssrdf("principled", Normal, Subsurface * SubsurfaceRadius, SubsurfaceColor, "roughness", Roughness);
-		} else {
+			if (subsurface_method == "burley") {
+				BSDF = mixed_ss_base_color * bssrdf("principled", Normal, Subsurface * SubsurfaceRadius, SubsurfaceColor, "roughness", Roughness);
+			}
+			else {
+				BSDF = mixed_ss_base_color * bssrdf("principled_random_walk", Normal, Subsurface * SubsurfaceRadius, SubsurfaceColor, "roughness", Roughness);
+			}
+		}
+		else {
 			BSDF = BaseColor * principled_diffuse(Normal, Roughness);
 		}
 
diff --git a/intern/cycles/kernel/shaders/node_subsurface_scattering.osl b/intern/cycles/kernel/shaders/node_subsurface_scattering.osl
index c9983fcd5dd..0df3256e1fd 100644
--- a/intern/cycles/kernel/shaders/node_subsurface_scattering.osl
+++ b/intern/cycles/kernel/shaders/node_subsurface_scattering.osl
@@ -30,7 +30,9 @@ shader node_subsurface_scattering(
 		BSSRDF = Color * bssrdf("gaussian", Normal, Scale * Radius, Color, "texture_blur", TextureBlur);
 	else if (falloff == "cubic")
 		BSSRDF = Color * bssrdf("cubic", Normal, Scale * Radius, Color, "texture_blur", TextureBlur, "sharpness", Sharpness);
-	else
+	else if (falloff == "burley")
 		BSSRDF = Color * bssrdf("burley", Normal, Scale * Radius, Color, "texture_blur", TextureBlur);
+	else
+		BSSRDF = Color * bssrdf("random_walk", Normal, Scale * Radius, Color, "texture_blur", TextureBlur);
 }
 
diff --git a/intern/cycles/kernel/split/kernel_subsurface_scatter.h b/intern/cycles/kernel/split/kernel_subsurface_scatter.h
index f902d000918..e50d63ea3bc 100644
--- a/intern/cycles/kernel/split/kernel_subsurface_scatter.h
+++ b/intern/cycles/kernel/split/kernel_subsurface_scatter.h
@@ -85,6 +85,7 @@ ccl_device_noinline bool kernel_split_branched_path_subsurface_indirect_light_it
 				branched_state->num_hits = subsurface_scatter_multi_intersect(kg,
 				                                                              &ss_isect_private,
 				                                                              sd,
+				                                                              hit_state,
 				                                                              sc,
 				                                                              &lcg_state,
 				                                                              bssrdf_u, bssrdf_v,
diff --git a/intern/cycles/kernel/svm/svm_closure.h b/intern/cycles/kernel/svm/svm_closure.h
index 578434792e4..fa43e1b60d0 100644
--- a/intern/cycles/kernel/svm/svm_closure.h
+++ b/intern/cycles/kernel/svm/svm_closure.h
@@ -114,7 +114,8 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 			float transmission_roughness = stack_load_float(stack, transmission_roughness_offset);
 			float eta = fmaxf(stack_load_float(stack, eta_offset), 1e-5f);
 
-			ClosureType distribution = stack_valid(data_node2.y) ? (ClosureType) data_node2.y : CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID;
+			ClosureType distribution = (ClosureType) data_node2.y;
+			ClosureType subsurface_method = (ClosureType) data_node2.z;
 
 			/* rotate tangent */
 			if(anisotropic_rotation != 0.0f)
@@ -193,7 +194,7 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 						bssrdf->roughness = roughness;
 
 						/* setup bsdf */
-						sd->flag |= bssrdf_setup(sd, bssrdf, (ClosureType)CLOSURE_BSSRDF_PRINCIPLED_ID);
+						sd->flag |= bssrdf_setup(sd, bssrdf, subsurface_method);
 					}
 				}
 			}
@@ -764,7 +765,8 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 #ifdef __SUBSURFACE__
 		case CLOSURE_BSSRDF_CUBIC_ID:
 		case CLOSURE_BSSRDF_GAUSSIAN_ID:
-		case CLOSURE_BSSRDF_BURLEY_ID: {
+		case CLOSURE_BSSRDF_BURLEY_ID:
+		case CLOSURE_BSSRDF_RANDOM_WALK_ID: {
 			float3 weight = sd->svm_closure_weight * mix_weight;
 			Bssrdf *bssrdf = bssrdf_alloc(sd, weight);
 
@@ -780,6 +782,7 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 				bssrdf->texture_blur = param2;
 				bssrdf->sharpness = stack_load_float(stack, data_node.w);
 				bssrdf->N = N;
+				bssrdf->roughness = 0.0f;
 				sd->flag |= bssrdf_setup(sd, bssrdf, (ClosureType)type);
 			}
 
diff --git a/intern/cycles/kernel/svm/svm_types.h b/intern/cycles/kernel/svm/svm_types.h
index 6ff04a65462..c0ce0f52cd0 100644
--- a/intern/cycles/kernel/svm/svm_types.h
+++ b/intern/cycles/kernel/svm/svm_types.h
@@ -446,6 +446,8 @@ typedef enum ClosureType {
 	CLOSURE_BSSRDF_GAUSSIAN_ID,
 	CLOSURE_BSSRDF_PRINCIPLED_ID,
 	CLOSURE_BSSRDF_BURLEY_ID,
+	CLOSURE_BSSRDF_RANDOM_WALK_ID,
+	CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID,
 
 	/* Other */
 	CLOSURE_HOLDOUT_ID,
@@ -477,8 +479,9 @@ typedef enum ClosureType {
 #define CLOSURE_IS_BSDF_MICROFACET(type) ((type >= CLOSURE_BSDF_MICROFACET_GGX_ID && type <= CLOSURE_BSDF_ASHIKHMIN_SHIRLEY_ANISO_ID) ||\
                                           (type >= CLOSURE_BSDF_MICROFACET_BECKMANN_REFRACTION_ID && type <= CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID) ||\
                                           (type == CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_FRESNEL_ID))
-#define CLOSURE_IS_BSDF_OR_BSSRDF(type) (type <= CLOSURE_BSSRDF_BURLEY_ID)
-#define CLOSURE_IS_BSSRDF(type) (type >= CLOSURE_BSSRDF_CUBIC_ID && type <= CLOSURE_BSSRDF_BURLEY_ID)
+#define CLOSURE_IS_BSDF_OR_BSSRDF(type) (type <= CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID)
+#define CLOSURE_IS_BSSRDF(type) (type >= CLOSURE_BSSRDF_CUBIC_ID && type <= CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID)
+#define CLOSURE_IS_DISK_BSSRDF(type) (type >= CLOSURE_BSSRDF_CUBIC_ID && type <= CLOSURE_BSSRDF_BURLEY_ID)
 #define CLOSURE_IS_VOLUME(type) (type >= CLOSURE_VOLUME_ID && type <= CLOSURE_VOLUME_HENYEY_GREENSTEIN_ID)
 #define CLOSURE_IS_VOLUME_SCATTER(type) (type == CLOSURE_VOLUME_HENYEY_GREENSTEIN_ID)
 #define CLOSURE_IS_VOLUME_ABSORPTION(type) (type == CLOSURE_VOLUME_ABSORPTION_ID)
diff --git a/intern/cycles/render/integrator.cpp b/intern/cycles/render/integrator.cpp
index 0dc1a9aa053..c337c19ced1 100644
--- a/intern/cycles/render/integrator.cpp
+++ b/intern/cycles/render/integrator.cpp
@@ -187,7 +187,10 @@ void Integrator::device_update(Device *device, DeviceScene *dscene, Scene *scene
 		max_samples = max(max_samples, volume_samples);
 	}
 
-	max_samples *= (max_bounce + transparent_max_bounce + 3 + BSSRDF_MAX_HITS);
+	uint total_bounces = max_bounce + transparent_max_bounce + 3 +
+	                     max(BSSRDF_MAX_HITS, BSSRDF_MAX_BOUNCES);
+
+	max_samples *= total_bounces;
 
 	int dimensions = PRNG_BASE_NUM + max_samples*PRNG_BOUNCE_NUM;
 	dimensions = min(dimensions, SOBOL_MAX_DIMENSIONS);
diff --git a/intern/cycles/render/nodes.cpp b/intern/cycles/render/nodes.cpp
index acfe07bf112..cb884ba9231 100644
--- a/intern/cycles/render/nodes.cpp
+++ b/intern/cycles/render/nodes.cpp
@@ -2312,6 +2312,12 @@ NODE_DEFINE(PrincipledBsdfNode)
 	distribution_enum.insert("GGX", CLOSURE_BSDF_MICROFACET_GGX_GLASS_ID);
 	distribution_enum.insert("Multiscatter GGX", CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID);
 	SOCKET_ENUM(distribution, "Distribution", distribution_enum, CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID);
+
+	static NodeEnum subsurface_method_enum;
+	subsurface_method_enum.insert("burley", CLOSURE_BSSRDF_PRINCIPLED_ID);
+	subsurface_method_enum.insert("random_walk", CLOSURE_BSSRDF_PRINCIPLED_RANDOM_WALK_ID);
+	SOCKET_ENUM(subsurface_method, "Subsurface Method", subsurface_method_enum, CLOSURE_BSSRDF_PRINCIPLED_ID);
+
 	SOCKET_IN_COLOR(base_color, "Base Color", make_float3(0.8f, 0.8f, 0.8f));
 	SOCKET_IN_COLOR(subsurface_color, "Subsurface Color", make_float3(0.8f, 0.8f, 0.8f));
 	SOCKET_IN_FLOAT(metallic, "Metallic", 0.0f);
@@ -2410,7 +2416,7 @@ void PrincipledBsdfNode::compile(SVMCompiler& compiler, ShaderInput *p_metallic,
 		compiler.encode_uchar4(sheen_offset, sheen_tint_offset, clearcoat_offset, clearcoat_roughness_offset));
 
 	compiler.add_node(compiler.encode_uchar4(ior_offset, transmission_offset, anisotropic_rotation_offset, transmission_roughness_offset),
-		distribution, SVM_STACK_INVALID, SVM_STACK_INVALID);
+		distribution, subsurface_method, SVM_STACK_INVALID);
 
 	float3 bc_default = get_float3(base_color_in->socket_type);
 
@@ -2442,6 +2448,7 @@ void PrincipledBsdfNode::compile(SVMCompiler& compiler)
 void PrincipledBsdfNode::compile(OSLCompiler& compiler)
 {
 	compiler.parameter(this, "distribution");
+	compiler.parameter(this, "subsurface_method");
 	compiler.add(this, "node_principled_bsdf");
 }
 
@@ -2525,6 +2532,7 @@ NODE_DEFINE(SubsurfaceScatteringNode)
 	falloff_enum.insert("cubic", CLOSURE_BSSRDF_CUBIC_ID);
 	falloff_enum.insert("gaussian", CLOSURE_BSSRDF_GAUSSIAN_ID);
 	falloff_enum.insert("burley", CLOSURE_BSSRDF_BURLEY_ID);
+	falloff_enum.insert("random_walk", CLOSURE_BSSRDF_RANDOM_WALK_ID);
 	SOCKET_ENUM(falloff, "Falloff", falloff_enum, CLOSURE_BSSRDF_BURLEY_ID);
 	SOCKET_IN_FLOAT(scale, "Scale", 0.01f);
 	SOCKET_IN_VECTOR(radius, "Radius", make_float3(0.1f, 0.1f, 0.1f));
diff --git a/intern/cycles/render/nodes.h b/intern/cycles/render/nodes.h
index a00b48ca5bc..f664ebf545d 100644
--- a/intern/cycles/render/nodes.h
+++ b/intern/cycles/render/nodes.h
@@ -390,6 +390,7 @@ public:
 	float3 normal, clearcoat_normal, tangent;
 	float surface_mix_weight;
 	ClosureType distribution, distribution_orig;
+	ClosureType subsurface_method;
 
 	bool has_integrator_dependency();
 	void attributes(Shader *shader, AttributeRequestSet *attributes);
diff --git a/intern/cycles/util/util_math_float3.h b/intern/cycles/util/util_math_float3.h
index e73e5bc17a2..e0a89b539b0 100644
--- a/intern/cycles/util/util_math_float3.h
+++ b/intern/cycles/util/util_math_float3.h
@@ -59,6 +59,7 @@ ccl_device_inline float3 mix(const float3& a, const float3& b, float t);
 ccl_device_inline float3 rcp(const float3& a);
 #endif  /* !__KERNEL_OPENCL__ */
 
+ccl_device_inline float min3(float3 a);
 ccl_device_inline float max3(float3 a);
 ccl_device_inline float len(const float3 a);
 ccl_device_inline float len_squared(const float3 a);
@@ -285,6 +286,11 @@ ccl_device_inline float3 rcp(const float3& a)
 }
 #endif  /* !__KERNEL_OPENCL__ */
 
+ccl_device_inline float min3(float3 a)
+{
+	return min(min(a.x, a.y), a.z);
+}
+
 ccl_device_inline float max3(float3 a)
 {
 	return max(max(a.x, a.y), a.z);
diff --git a/source/blender/editors/space_node/drawnode.c b/source/blender/editors/space_node/drawnode.c
index ead0e848716..27c4eb2fa3c 100644
--- a/source/blender/editors/space_node/drawnode.c
+++ b/source/blender/editors/space_node/drawnode.c
@@ -1032,6 +1032,12 @@ static void node_shader_buts_glossy(uiLayout *layout, bContext *UNUSED(C), Point
 	uiItemR(layout, ptr, "distribution", 0, "", ICON_NONE);
 }
 
+static void node_shader_buts_principled(uiLayout *layout, bContext *UNUSED(C), PointerRNA *ptr)
+{
+	uiItemR(layout, ptr, "distribution", 0, "", ICON_NONE);
+	uiItemR(layout, ptr, "subsurface_method", 0, "", ICON_NONE);
+}
+
 static void node_shader_buts_anisotropic(uiLayout *layout, bContext *UNUSED(C), PointerRNA *ptr)
 {
 	uiItemR(layout, ptr, "distribution", 0, "", ICON_NONE);
@@ -1204,9 +1210,11 @@ static void node_shader_set_butfunc(bNodeType *ntype)
 		case SH_NODE_BSDF_GLOSSY:
 		case SH_NODE_BSDF_GLASS:
 		case SH_NODE_BSDF_REFRACTION:
-		case SH_NODE_BSDF_PRINCIPLED:
 			ntype->draw_buttons = node_shader_buts_glossy;
 			break;
+		case SH_NODE_BSDF_PRINCIPLED:
+			ntype->draw_buttons = node_shader_buts_principled;
+			break;
 		case SH_NODE_BSDF_ANISOTROPIC:
 			ntype->draw_buttons = node_shader_buts_anisotropic;
 			break;
diff --git a/source/blender/gpu/intern/gpu_material.c b/source/blender/gpu/intern/gpu_material.c
index 2e6c1cbf9df..4efb8ca07b2 100644
--- a/source/blender/gpu/intern/gpu_material.c
+++ b/source/blender/gpu/intern/gpu_material.c
@@ -557,7 +557,9 @@ static float eval_profile(float r, short falloff_type, float sharpness, float pa
 {
 	r = fabsf(r);
 
-	if (falloff_type == SHD_SUBSURFACE_BURLEY) {
+	if (falloff_type == SHD_SUBSURFACE_BURLEY ||
+	    falloff_type == SHD_SUBSURFACE_RANDOM_WALK)
+	{
 		return burley_profile(r, param) / BURLEY_TRUNCATE_CDF;
 	}
 	else if (falloff_type == SHD_SUBSURFACE_CUBIC) {
@@ -598,7 +600,9 @@ static void compute_sss_kernel(
 	/* Christensen-Burley fitting */
 	float l[3], d[3];
 
-	if (falloff_type == SHD_SUBSURFACE_BURLEY) {
+	if (falloff_type == SHD_SUBSURFACE_BURLEY ||
+	    falloff_type == SHD_SUBSURFACE_RANDOM_WALK)
+	{
 		mul_v3_v3fl(l, rad, 0.25f * M_1_PI);
 		const float A = 1.0f;
 		const float s = 1.9f - A + 3.5f * (A - 0.8f) * (A - 0.8f);
diff --git a/source/blender/makesdna/DNA_node_types.h b/source/blender/makesdna/DNA_node_types.h
index bddbf63d8b1..63eecebf44a 100644
--- a/source/blender/makesdna/DNA_node_types.h
+++ b/source/blender/makesdna/DNA_node_types.h
@@ -1077,6 +1077,7 @@ enum {
 	SHD_SUBSURFACE_CUBIC			= 1,
 	SHD_SUBSURFACE_GAUSSIAN			= 2,
 	SHD_SUBSURFACE_BURLEY			= 3,
+	SHD_SUBSURFACE_RANDOM_WALK		= 4,
 };
 
 /* blur node */
diff --git a/source/blender/makesrna/intern/rna_nodetree.c b/source/blender/makesrna/intern/rna_nodetree.c
index 59570c10ebd..a148efb859d 100644
--- a/source/blender/makesrna/intern/rna_nodetree.c
+++ b/source/blender/makesrna/intern/rna_nodetree.c
@@ -3312,11 +3312,18 @@ static const EnumPropertyItem node_script_mode_items[] = {
 };
 
 static const EnumPropertyItem node_principled_distribution_items[] = {
-	{ SHD_GLOSSY_GGX, "GGX", 0, "GGX", "" },
-	{ SHD_GLOSSY_MULTI_GGX, "MULTI_GGX", 0, "Multiscatter GGX", "" },
+	{SHD_GLOSSY_GGX, "GGX", 0, "GGX", ""},
+	{SHD_GLOSSY_MULTI_GGX, "MULTI_GGX", 0, "Multiscatter GGX", ""},
+	{0, NULL, 0, NULL, NULL}
+};
+
+static const EnumPropertyItem node_subsurface_method_items[] = {
+	{SHD_SUBSURFACE_BURLEY, "BURLEY", 0, "Christensen-Burley", "Approximation to physically based volume scattering"},
+	{SHD_SUBSURFACE_RANDOM_WALK, "RANDOM_WALK", 0, "Random Walk", "Volumetric approximation to physically based volume scattering"},
 	{ 0, NULL, 0, NULL, NULL }
 };
 
+
 /* -- Common nodes ---------------------------------------------------------- */
 
 static void def_group_input(StructRNA *srna)
@@ -4263,6 +4270,12 @@ static void def_principled(StructRNA *srna)
 	RNA_def_property_enum_items(prop, node_principled_distribution_items);
 	RNA_def_property_ui_text(prop, "Distribution", "");
 	RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_ShaderNodePrincipled_update");
+
+	prop = RNA_def_property(srna, "subsurface_method", PROP_ENUM, PROP_NONE);
+	RNA_def_property_enum_sdna(prop, NULL, "custom2");
+	RNA_def_property_enum_items(prop, node_subsurface_method_items);
+	RNA_def_property_ui_text(prop, "Subsurface Method", "Method for rendering subsurface scattering");
+	RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_ShaderNodePrincipled_update");
 }
 
 static void def_refraction(StructRNA *srna)
@@ -4466,6 +4479,7 @@ static void def_sh_subsurface(StructRNA *srna)
 		{SHD_SUBSURFACE_CUBIC, "CUBIC", 0, "Cubic", "Simple cubic falloff function"},
 		{SHD_SUBSURFACE_GAUSSIAN, "GAUSSIAN", 0, "Gaussian", "Normal distribution, multiple can be combined to fit more complex profiles"},
 		{SHD_SUBSURFACE_BURLEY, "BURLEY", 0, "Christensen-Burley", "Approximation to physically based volume scattering"},
+		{SHD_SUBSURFACE_RANDOM_WALK, "RANDOM_WALK", 0, "Random Walk", "Volumetric approximation to physically based volume scattering"},
 		{0, NULL, 0, NULL, NULL}
 	};
 
diff --git a/source/blender/nodes/shader/nodes/node_shader_bsdf_principled.c b/source/blender/nodes/shader/nodes/node_shader_bsdf_principled.c
index b0ca4128d19..1c291cb5b83 100644
--- a/source/blender/nodes/shader/nodes/node_shader_bsdf_principled.c
+++ b/source/blender/nodes/shader/nodes/node_shader_bsdf_principled.c
@@ -61,6 +61,7 @@ static bNodeSocketTemplate sh_node_bsdf_principled_out[] = {
 static void node_shader_init_principled(bNodeTree *UNUSED(ntree), bNode *node)
 {
 	node->custom1 = SHD_GLOSSY_MULTI_GGX;
+	node->custom2 = SHD_SUBSURFACE_BURLEY;
 }
 
 static int node_shader_gpu_bsdf_principled(GPUMaterial *mat, bNode *node, bNodeExecData *UNUSED(execdata), GPUNodeStack *in, GPUNodeStack *out)