Cycles: random walk subsurface scattering.

It is basically brute force volume scattering within the mesh, but part
of the SSS code for faster performance. The main difference with actual
volume scattering is that we assume the boundaries are diffuse and that
all lighting is coming through this boundary from outside the volume.

This gives much more accurate results for thin features and low density.
Some challenges remain however:

* Significantly more noisy than BSSRDF. Adding Dwivedi sampling may help
  here, but it's unclear still how much it helps in real world cases.
* Due to this being a volumetric method, geometry like eyes or mouth can
  darken the skin on the outside. We may be able to reduce this effect,
  or users can compensate for it by reducing the scattering radius in
  such areas.
* Sharp corners are quite bright. This matches actual volume rendering
  and results in some other renderers, but maybe not so much real world
  objects.

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

3 files changed, 28 insertions, 19 deletions

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..790368ee888 100644
--- a/intern/cycles/kernel/closure/bssrdf.h
+++ b/intern/cycles/kernel/closure/bssrdf.h
@@ -408,7 +408,8 @@ 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)
 		{
 			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__