Cycles: Subsurface Scattering

New features:

* Bump mapping now works with SSS
* Texture Blur factor for SSS, see the documentation for details:
http://wiki.blender.org/index.php/Doc:2.6/Manual/Render/Cycles/Nodes/Shaders#Subsurface_Scattering

Work in progress for feedback:

Initial implementation of the "BSSRDF Importance Sampling" paper, which uses
a different importance sampling method. It gives better quality results in
many ways, with the availability of both Cubic and Gaussian falloff functions,
but also tends to be more noisy when using the progressive integrator and does
not give great results with some geometry. It works quite well for the
non-progressive integrator and is often less noisy there.

This code may still change a lot, so unless you're testing it may be best to
stick to the Compatible falloff function.

Skin test render and file that takes advantage of the gaussian falloff:
http://www.pasteall.org/pic/show.php?id=57661
http://www.pasteall.org/pic/show.php?id=57662
http://www.pasteall.org/blend/23501

2 files changed, 137 insertions, 85 deletions

diff --git a/intern/cycles/kernel/closure/bsdf_microfacet.h b/intern/cycles/kernel/closure/bsdf_microfacet.h
index 915b9eafbc1..b159f585831 100644
--- a/intern/cycles/kernel/closure/bsdf_microfacet.h
+++ b/intern/cycles/kernel/closure/bsdf_microfacet.h
@@ -37,11 +37,6 @@ CCL_NAMESPACE_BEGIN
 
 /* GGX */
 
-__device_inline float safe_sqrtf(float f)
-{
-	return sqrtf(max(f, 0.0f));
-}
-
 __device int bsdf_microfacet_ggx_setup(ShaderClosure *sc)
 {
 	sc->data0 = clamp(sc->data0, 0.0f, 1.0f); /* m_ag */
diff --git a/intern/cycles/kernel/closure/bssrdf.h b/intern/cycles/kernel/closure/bssrdf.h
index 486de4ca65f..23b932a91c6 100644
--- a/intern/cycles/kernel/closure/bssrdf.h
+++ b/intern/cycles/kernel/closure/bssrdf.h
@@ -21,130 +21,187 @@
 
 CCL_NAMESPACE_BEGIN
 
-__device int bssrdf_setup(ShaderClosure *sc)
+__device int bssrdf_setup(ShaderClosure *sc, ClosureType type)
 {
 	if(sc->data0 < BSSRDF_MIN_RADIUS) {
 		/* revert to diffuse BSDF if radius too small */
 		sc->data0 = 0.0f;
 		sc->data1 = 0.0f;
-		return bsdf_diffuse_setup(sc);
+		int flag = bsdf_diffuse_setup(sc);
+		sc->type = CLOSURE_BSDF_BSSRDF_ID;
+		return flag;
 	}
 	else {
-		/* IOR param */
-		sc->data1 = max(sc->data1, 1.0f);
-		sc->type = CLOSURE_BSSRDF_ID;
+		sc->data1 = clamp(sc->data1, 0.0f, 1.0f); /* texture blur */
+		sc->type = type;
 
 		return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSSRDF;
 	}
 }
 
-/* Simple Cubic BSSRDF falloff */
+/* Planar Truncated Gaussian
+ *
+ * Note how this is different from the typical gaussian, this one integrates
+ * to 1 over the plane (where you get an extra 2*pi*x factor). We are lucky
+ * that integrating x*exp(-x) gives a nice closed form solution. */
+
+/* paper suggests 1/12.46 which is much too small, suspect it's *12.46 */
+#define GAUSS_TRUNCATE 12.46f
 
-__device float bssrdf_cubic(float ld, float r)
+__device float bssrdf_gaussian_eval(ShaderClosure *sc, float r)
 {
-	if(ld == 0.0f)
-		return (r == 0.0f)? 1.0f: 0.0f;
+	/* integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) from 0 to Rm
+	 * = 1 - exp(-Rm*Rm/(2*v)) */
+	const float v = sc->data0;
+	const float Rm = sqrtf(v*GAUSS_TRUNCATE);
+
+	if(r >= Rm)
+		return 0.0f;
 
-	return powf(ld - min(r, ld), 3.0f) * 4.0f/powf(ld, 4.0f);
+	return expf(-r*r/(2.0f*v))/(2.0f*M_PI_F*v);
 }
 
-/* Original BSSRDF fallof function */
-
-typedef struct BSSRDFParams {
-	float eta;		/* index of refraction */
-	float sigma_t_; /* reduced extinction coefficient */
-	float sigma_tr;	/* effective extinction coefficient */
-	float Fdr;		/* diffuse fresnel reflectance */
-	float D;		/* diffusion constant */
-	float A;
-	float alpha_;	/* reduced albedo */
-	float zr;		/* distance of virtual lightsource above surface */
-	float zv;		/* distance of virtual lightsource below surface */
-	float ld;		/* mean free path */
-	float ro;		/* diffuse reflectance */
-} BSSRDFParams;
-
-__device float bssrdf_reduced_albedo_Rd(float alpha_, float A, float ro)
+__device float bssrdf_gaussian_pdf(ShaderClosure *sc, float r)
 {
-	float sq;
+	/* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
+	const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
+
+	return bssrdf_gaussian_eval(sc, r) * (1.0f/(area_truncated));
+}
+
+__device void bssrdf_gaussian_sample(ShaderClosure *sc, float xi, float *r, float *h)
+{
+	/* xi = integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) = -exp(-r^2/(2*v))
+	 * r = sqrt(-2*v*logf(xi)) */
+
+	const float v = sc->data0;
+	const float Rm = sqrtf(v*GAUSS_TRUNCATE);
+
+	/* 1.0 - expf(-Rm*Rm/(2*v)) simplified */
+	const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE);
+
+	/* r(xi) */
+	const float r_squared = -2.0f*v*logf(1.0f - xi*area_truncated);
+	*r = sqrtf(r_squared);
+
+	 /* h^2 + r^2 = Rm^2 */
+	 *h = sqrtf(Rm*Rm - r_squared);
+}
+
+/* Planar Cubic BSSRDF falloff
+ *
+ * This is basically (Rm - x)^3, with some factors to normalize it. For sampling
+ * we integrate 2*pi*x * (Rm - x)^3, which gives us a quintic equation that as
+ * far as I can tell has no closed form solution. So we get an iterative solution
+ * instead with newton-raphson. */
+
+__device float bssrdf_cubic_eval(ShaderClosure *sc, float r)
+{
+	const float Rm = sc->data0;
+
+	if(r >= Rm)
+		return 0.0f;
+	
+	/* integrate (2*pi*r * 10*(R - r)^3)/(pi * R^5) from 0 to R = 1 */
+	const float Rm5 = (Rm*Rm) * (Rm*Rm) * Rm;
+	const float f = Rm - min(r, Rm);
+	const float f3 = f*f*f;
 
-	sq = sqrtf(3.0f*(1.0f - alpha_));
-	return (alpha_/2.0f)*(1.0f + expf((-4.0f/3.0f)*A*sq))*expf(-sq) - ro;
+	return (f3 * 10.0f) / (Rm5 * M_PI_F);
 }
 
-__device float bssrdf_compute_reduced_albedo(float A, float ro)
+__device float bssrdf_cubic_pdf(ShaderClosure *sc, float r)
 {
-	const float tolerance = 1e-8f;
-	const int max_iteration_count = 20;
-	float d, fsub, xn_1 = 0.0f, xn = 1.0f, fxn, fxn_1;
+	return bssrdf_cubic_eval(sc, r);
+}
+
+/* solve 10x^2 - 20x^3 + 15x^4 - 4x^5 - xi == 0 */
+__device float bssrdf_cubic_quintic_root_find(float xi)
+{
+	/* newton-raphson iteration, usually succeeds in 2-4 iterations, except
+	 * outside 0.02 ... 0.98 where it can go up to 10, so overall performance
+	 * should not be too bad */
+	const float tolerance = 1e-6f;
+	const int max_iteration_count = 10;
+	float x = 0.25f;
 	int i;
 
-	/* use secant method to compute reduced albedo using Rd function inverse
-	 * with a given reflectance */
-	fxn = bssrdf_reduced_albedo_Rd(xn, A, ro);
-	fxn_1 = bssrdf_reduced_albedo_Rd(xn_1, A, ro);
+	for (i = 0; i < max_iteration_count; i++) {
+		float x2 = x*x;
+		float x3 = x2*x;
+		float nx = (1.0f - x);
 
-	for (i= 0; i < max_iteration_count; i++) {
-		fsub = (fxn - fxn_1);
-		if (fabsf(fsub) < tolerance)
-			break;
-		d = ((xn - xn_1)/fsub)*fxn;
-		if (fabsf(d) < tolerance)
-			break;
+		float f = 10.0f*x2 - 20.0f*x3 + 15.0f*x2*x2 - 4.0f*x2*x3 - xi;
+		float f_ = 20.0f*(x*nx)*(nx*nx);
 
-		xn_1 = xn;
-		fxn_1 = fxn;
-		xn = xn - d;
+		if(fabsf(f) < tolerance || f_ == 0.0f)
+			break;
 
-		if (xn > 1.0f) xn = 1.0f;
-		if (xn_1 > 1.0f) xn_1 = 1.0f;
-		
-		fxn = bssrdf_reduced_albedo_Rd(xn, A, ro);
+		x = clamp(x - f/f_, 0.0f, 1.0f);
 	}
 
-	/* avoid division by zero later */
-	if (xn <= 0.0f)
-		xn = 0.00001f;
-
-	return xn;
+	return x;
 }
 
-__device void bssrdf_setup_params(BSSRDFParams *ss, float refl, float radius, float ior)
+__device void bssrdf_cubic_sample(ShaderClosure *sc, float xi, float *r, float *h)
 {
-	ss->eta = ior;
-	ss->Fdr = -1.440f/ior*ior + 0.710f/ior + 0.668f + 0.0636f*ior;
-	ss->A = (1.0f + ss->Fdr)/(1.0f - ss->Fdr);
-	ss->ld = radius;
-	ss->ro = min(refl, 0.999f);
+	const float Rm = sc->data0;
+	const float r_ = bssrdf_cubic_quintic_root_find(xi) * Rm;
 
-	ss->alpha_ = bssrdf_compute_reduced_albedo(ss->A, ss->ro);
+	*r = r_;
 
-	ss->sigma_tr = 1.0f/ss->ld;
-	ss->sigma_t_ = ss->sigma_tr/sqrtf(3.0f*(1.0f - ss->alpha_));
+	/* h^2 + r^2 = Rm^2 */
+	*h = sqrtf(Rm*Rm - r_*r_);
+}
 
-	ss->D = 1.0f/(3.0f*ss->sigma_t_);
+/* None BSSRDF falloff
+ * 
+ * Samples distributed over disk with no falloff, for reference. */
 
-	ss->zr = 1.0f/ss->sigma_t_;
-	ss->zv = ss->zr + 4.0f*ss->A*ss->D;
+__device float bssrdf_none_eval(ShaderClosure *sc, float r)
+{
+	const float Rm = sc->data0;
+	return (r < Rm)? 1.0f: 0.0f;
 }
 
-/* exponential falloff function */
+__device float bssrdf_none_pdf(ShaderClosure *sc, float r)
+{
+	/* integrate (2*pi*r)/(pi*Rm*Rm) from 0 to Rm = 1 */
+	const float Rm = sc->data0;
+	const float area = (M_PI_F*Rm*Rm);
+
+	return bssrdf_none_eval(sc, r) / area;
+}
 
-__device float bssrdf_original(const BSSRDFParams *ss, float r)
+__device void bssrdf_none_sample(ShaderClosure *sc, float xi, float *r, float *h)
 {
-	if(ss->ld == 0.0f)
-		return (r == 0.0f)? 1.0f: 0.0f;
+	/* xi = integrate (2*pi*r)/(pi*Rm*Rm) = r^2/Rm^2
+	 * r = sqrt(xi)*Rm */
+	const float Rm = sc->data0;
+	const float r_ = sqrtf(xi)*Rm;
+
+	*r = r_;
 
-	float rr = r*r;
-	float sr, sv, Rdr, Rdv;
+	/* h^2 + r^2 = Rm^2 */
+	*h = sqrtf(Rm*Rm - r_*r_);
+}
 
-	sr = sqrtf(rr + ss->zr*ss->zr);
-	sv = sqrtf(rr + ss->zv*ss->zv);
+/* Generic */
 
-	Rdr = ss->zr*(1.0f + ss->sigma_tr*sr)*expf(-ss->sigma_tr*sr)/(sr*sr*sr);
-	Rdv = ss->zv*(1.0f + ss->sigma_tr*sv)*expf(-ss->sigma_tr*sv)/(sv*sv*sv);
+__device void bssrdf_sample(ShaderClosure *sc, float xi, float *r, float *h)
+{
+	if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
+		bssrdf_cubic_sample(sc, xi, r, h);
+	else
+		bssrdf_gaussian_sample(sc, xi, r, h);
+}
 
-	return ss->alpha_*(1.0f/M_4PI_F)*(Rdr + Rdv);
+__device float bssrdf_pdf(ShaderClosure *sc, float r)
+{
+	if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
+		return bssrdf_cubic_pdf(sc, r);
+	else
+		return bssrdf_gaussian_pdf(sc, r);
 }
 
 CCL_NAMESPACE_END