Cycles: Implement approximate reflectance profiles

Using this paper:

  http://graphics.pixar.com/library/ApproxBSSRDF/paper.pdf

This model gives less blurry results than the Cubic and Gaussian
we had implemented:

- Cubic: https://developer.blender.org/F279670
- Burley: https://developer.blender.org/F279671

The model is called "Christensen-Burley" in the interface, which
actually should be read as "Physically based" or "Realistic".

Reviewers: juicyfruit, dingto, lukasstockner97, brecht

Reviewed By: brecht, dingto

Subscribers: robocyte

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

15 files changed, 174 insertions, 10 deletions

diff --git a/intern/cycles/app/cycles_xml.cpp b/intern/cycles/app/cycles_xml.cpp
index b366e937d94..dd1dae121b1 100644
--- a/intern/cycles/app/cycles_xml.cpp
+++ b/intern/cycles/app/cycles_xml.cpp
@@ -596,8 +596,10 @@ static void xml_read_shader_graph(const XMLReadState& state, Shader *shader, pug
 			xml_read_string(&falloff, node, "falloff");
 			if(falloff == "cubic")
 				sss->closure = CLOSURE_BSSRDF_CUBIC_ID;
-			else
+			else if(falloff == "gaussian")
 				sss->closure = CLOSURE_BSSRDF_GAUSSIAN_ID;
+			else /*if(falloff == "burley")*/
+				sss->closure = CLOSURE_BSSRDF_BURLEY_ID;
 
 			snode = sss;
 		}
diff --git a/intern/cycles/blender/blender_shader.cpp b/intern/cycles/blender/blender_shader.cpp
index 8870507d00a..a40eaf52480 100644
--- a/intern/cycles/blender/blender_shader.cpp
+++ b/intern/cycles/blender/blender_shader.cpp
@@ -405,6 +405,9 @@ static ShaderNode *add_node(Scene *scene,
 			case BL::ShaderNodeSubsurfaceScattering::falloff_GAUSSIAN:
 				subsurface->closure = CLOSURE_BSSRDF_GAUSSIAN_ID;
 				break;
+			case BL::ShaderNodeSubsurfaceScattering::falloff_BURLEY:
+				subsurface->closure = CLOSURE_BSSRDF_BURLEY_ID;
+				break;
 		}
 
 		node = subsurface;
diff --git a/intern/cycles/kernel/closure/bssrdf.h b/intern/cycles/kernel/closure/bssrdf.h
index b986ab6ba39..041832ee5b7 100644
--- a/intern/cycles/kernel/closure/bssrdf.h
+++ b/intern/cycles/kernel/closure/bssrdf.h
@@ -192,6 +192,110 @@ ccl_device void bssrdf_cubic_sample(ShaderClosure *sc, float xi, float *r, float
 	*h = sqrtf(Rm*Rm - r_*r_);
 }
 
+/* Approximate Reflectance Profiles
+ * http://graphics.pixar.com/library/ApproxBSSRDF/paper.pdf
+ */
+
+ccl_device_inline float bssrdf_burley_fitting(float A)
+{
+	/* Diffuse surface transmission, equation (6). */
+	return 1.9f - A + 3.5f * (A - 0.8f) * (A - 0.8f);
+}
+
+/* Scale mean free path length so it gives similar looking result
+ * to Cubic and Gaussian models.
+ */
+ccl_device_inline float bssrdf_burley_compatible_mfp(float r)
+{
+	return 0.5f * M_1_PI_F * r;
+}
+
+ccl_device float bssrdf_burley_eval(ShaderClosure *sc, float r)
+{
+	/* Mean free path length. */
+	const float l = bssrdf_burley_compatible_mfp(sc->data0);
+	/* Surface albedo. */
+	const float A = sc->data2;
+	const float s = bssrdf_burley_fitting(A);
+	/* Burley refletance profile, equation (3).
+	 *
+	 * Note that surface albedo is already included into sc->weight, no need to
+	 * multiply by this term here.
+	 */
+	float exp_r_3_d = expf(-s*r / (3.0f * l));
+	float exp_r_d = exp_r_3_d * exp_r_3_d * exp_r_3_d;
+	return s * (exp_r_d + exp_r_3_d) / (8*M_PI_F*l*r);
+}
+
+ccl_device float bssrdf_burley_pdf(ShaderClosure *sc, float r)
+{
+	return bssrdf_burley_eval(sc, r);
+}
+
+/* Find the radius for desired CDF value.
+ * Returns scaled radius, meaning the result is to be scaled up by d.
+ * Since there's no closed form solution we do Newton-Raphson method to find it.
+ */
+ccl_device float bssrdf_burley_root_find(float xi)
+{
+	const float tolerance = 1e-6f;
+	const int max_iteration_count = 10;
+	/* Do initial guess based on manual curve fitting, this allows us to reduce
+	 * number of iterations to maximum 4 across the [0..1] range. We keep maximum
+	 * number of iteration higher just to be sure we didn't miss root in some
+	 * corner case.
+	 */
+	float r;
+	if (xi <= 0.9f) {
+		r = expf(xi * xi * 2.4f) - 1.0f;
+	}
+	else {
+		float a = expf(xi * xi * 4.0f) - 1.0f;
+		r = a*a;
+	}
+	/* Solve against scaled radius. */
+	for(int i = 0; i < max_iteration_count; i++) {
+		float exp_r_3 = expf(-r / 3.0f);
+		float exp_r = exp_r_3 * exp_r_3 * exp_r_3;
+		float f = 1.0f - 0.25f * exp_r - 0.75f * exp_r_3 - xi;
+		float f_ = 0.25f * exp_r + 0.25f * exp_r_3;
+
+		if(fabsf(f) < tolerance || f_ == 0.0f) {
+			break;
+		}
+
+		r = r - f/f_;
+		if(r < 0.0f) {
+			r = 0.0f;
+		}
+	}
+	return r;
+}
+
+ccl_device void bssrdf_burley_sample(ShaderClosure *sc,
+                                     float xi,
+                                     float *r,
+                                     float *h)
+{
+	/* Mean free path length. */
+	const float l = bssrdf_burley_compatible_mfp(sc->data0);
+	/* Surface albedo. */
+	const float A = sc->data2;
+	const float s = bssrdf_burley_fitting(A);
+	const float d = l / s;
+	/* This is a bit arbitrary, just need big enough radius so it matches
+	 * the mean free length, but still not too big so sampling is still
+	 * effective. Might need some further tweaks.
+	 */
+	const float Rm = 10.0f*d;
+	const float r_ = bssrdf_burley_root_find(xi) * d;
+
+	*r = r_;
+
+	/* h^2 + r^2 = Rm^2 */
+	*h = sqrtf(Rm*Rm - r_*r_);
+}
+
 /* None BSSRDF falloff
  *
  * Samples distributed over disk with no falloff, for reference. */
@@ -230,16 +334,20 @@ ccl_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
+	else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
 		bssrdf_gaussian_sample(sc, xi, r, h);
+	else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID)*/
+		bssrdf_burley_sample(sc, xi, r, h);
 }
 
 ccl_device float bssrdf_pdf(ShaderClosure *sc, float r)
 {
 	if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
 		return bssrdf_cubic_pdf(sc, r);
-	else
+	else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
 		return bssrdf_gaussian_pdf(sc, r);
+	else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID)*/
+		return bssrdf_burley_pdf(sc, r);
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/osl/osl_bssrdf.cpp b/intern/cycles/kernel/osl/osl_bssrdf.cpp
index 7f286c76c46..da4afb138f6 100644
--- a/intern/cycles/kernel/osl/osl_bssrdf.cpp
+++ b/intern/cycles/kernel/osl/osl_bssrdf.cpp
@@ -105,5 +105,34 @@ ClosureParam *closure_bssrdf_gaussian_params()
 
 CCLOSURE_PREPARE(closure_bssrdf_gaussian_prepare, GaussianBSSRDFClosure)
 
+/* Burley */
+
+class BurleyBSSRDFClosure : public CBSSRDFClosure {
+public:
+	BurleyBSSRDFClosure()
+	{}
+
+	void setup()
+	{
+		sc.type = CLOSURE_BSSRDF_BURLEY_ID;
+		sc.data0 = fabsf(average(radius));
+	}
+};
+
+ClosureParam *closure_bssrdf_burley_params()
+{
+	static ClosureParam params[] = {
+		CLOSURE_FLOAT3_PARAM(BurleyBSSRDFClosure, sc.N),
+		CLOSURE_FLOAT3_PARAM(BurleyBSSRDFClosure, radius),
+		CLOSURE_FLOAT_PARAM(BurleyBSSRDFClosure, sc.data1),
+		CLOSURE_FLOAT3_PARAM(BurleyBSSRDFClosure, albedo),
+		CLOSURE_STRING_KEYPARAM(BurleyBSSRDFClosure, label, "label"),
+		CLOSURE_FINISH_PARAM(BurleyBSSRDFClosure)
+	};
+	return params;
+}
+
+CCLOSURE_PREPARE(closure_bssrdf_burley_prepare, BurleyBSSRDFClosure)
+
 CCL_NAMESPACE_END
 
diff --git a/intern/cycles/kernel/osl/osl_bssrdf.h b/intern/cycles/kernel/osl/osl_bssrdf.h
index 6aee2c28ea8..d81ecade543 100644
--- a/intern/cycles/kernel/osl/osl_bssrdf.h
+++ b/intern/cycles/kernel/osl/osl_bssrdf.h
@@ -49,6 +49,7 @@ class CBSSRDFClosure : public CClosurePrimitive {
 public:
 	ShaderClosure sc;
 	float3 radius;
+	float3 albedo;
 
 	CBSSRDFClosure() : CClosurePrimitive(BSSRDF) { }
 	int scattering() const { return LABEL_DIFFUSE; }
diff --git a/intern/cycles/kernel/osl/osl_closures.cpp b/intern/cycles/kernel/osl/osl_closures.cpp
index 461ce8f7598..95b8cea0922 100644
--- a/intern/cycles/kernel/osl/osl_closures.cpp
+++ b/intern/cycles/kernel/osl/osl_closures.cpp
@@ -236,6 +236,8 @@ void OSLShader::register_closures(OSLShadingSystem *ss_)
 		closure_bssrdf_cubic_params(), closure_bssrdf_cubic_prepare);
 	register_closure(ss, "bssrdf_gaussian", id++,
 		closure_bssrdf_gaussian_params(), closure_bssrdf_gaussian_prepare);
+	register_closure(ss, "bssrdf_burley", id++,
+		closure_bssrdf_burley_params(), closure_bssrdf_burley_prepare);
 
 	register_closure(ss, "hair_reflection", id++,
 		bsdf_hair_reflection_params(), bsdf_hair_reflection_prepare);
diff --git a/intern/cycles/kernel/osl/osl_closures.h b/intern/cycles/kernel/osl/osl_closures.h
index d8a3f77d0dd..526c03557f7 100644
--- a/intern/cycles/kernel/osl/osl_closures.h
+++ b/intern/cycles/kernel/osl/osl_closures.h
@@ -50,6 +50,7 @@ OSL::ClosureParam *closure_bsdf_diffuse_ramp_params();
 OSL::ClosureParam *closure_bsdf_phong_ramp_params();
 OSL::ClosureParam *closure_bssrdf_cubic_params();
 OSL::ClosureParam *closure_bssrdf_gaussian_params();
+OSL::ClosureParam *closure_bssrdf_burley_params();
 OSL::ClosureParam *closure_henyey_greenstein_volume_params();
 
 void closure_emission_prepare(OSL::RendererServices *, int id, void *data);
@@ -60,6 +61,7 @@ void closure_bsdf_diffuse_ramp_prepare(OSL::RendererServices *, int id, void *da
 void closure_bsdf_phong_ramp_prepare(OSL::RendererServices *, int id, void *data);
 void closure_bssrdf_cubic_prepare(OSL::RendererServices *, int id, void *data);
 void closure_bssrdf_gaussian_prepare(OSL::RendererServices *, int id, void *data);
+void closure_bssrdf_burley_prepare(OSL::RendererServices *, int id, void *data);
 void closure_henyey_greenstein_volume_prepare(OSL::RendererServices *, int id, void *data);
 
 #define CCLOSURE_PREPARE(name, classname)          \
diff --git a/intern/cycles/kernel/osl/osl_shader.cpp b/intern/cycles/kernel/osl/osl_shader.cpp
index 8f459c8b267..8acc042cd2c 100644
--- a/intern/cycles/kernel/osl/osl_shader.cpp
+++ b/intern/cycles/kernel/osl/osl_shader.cpp
@@ -281,11 +281,17 @@ static void flatten_surface_closure_tree(ShaderData *sd, int path_flag,
 							if(path_flag & PATH_RAY_DIFFUSE_ANCESTOR)
 								bssrdf->radius = make_float3(0.0f, 0.0f, 0.0f);
 
+							float3 albedo =
+							        (bssrdf->sc.type == CLOSURE_BSSRDF_BURLEY_ID)
+							                ? bssrdf->albedo
+							                : make_float3(0.0f, 0.0f, 0.0f);
+
 							/* create one closure for each color channel */
 							if(fabsf(weight.x) > 0.0f) {
 								sc.weight = make_float3(weight.x, 0.0f, 0.0f);
 								sc.data0 = bssrdf->radius.x;
 								sc.data1 = 0.0f;
+								sc.data2 = albedo.x;
 								sd->flag |= bssrdf_setup(&sc, sc.type);
 								sd->closure[sd->num_closure++] = sc;
 							}
@@ -294,6 +300,7 @@ static void flatten_surface_closure_tree(ShaderData *sd, int path_flag,
 								sc.weight = make_float3(0.0f, weight.y, 0.0f);
 								sc.data0 = bssrdf->radius.y;
 								sc.data1 = 0.0f;
+								sc.data2 = albedo.y;
 								sd->flag |= bssrdf_setup(&sc, sc.type);
 								sd->closure[sd->num_closure++] = sc;
 							}
@@ -302,6 +309,7 @@ static void flatten_surface_closure_tree(ShaderData *sd, int path_flag,
 								sc.weight = make_float3(0.0f, 0.0f, weight.z);
 								sc.data0 = bssrdf->radius.z;
 								sc.data1 = 0.0f;
+								sc.data2 = albedo.z;
 								sd->flag |= bssrdf_setup(&sc, sc.type);
 								sd->closure[sd->num_closure++] = sc;
 							}
diff --git a/intern/cycles/kernel/shaders/node_subsurface_scattering.osl b/intern/cycles/kernel/shaders/node_subsurface_scattering.osl
index dbbf657776c..0a9a98afd33 100644
--- a/intern/cycles/kernel/shaders/node_subsurface_scattering.osl
+++ b/intern/cycles/kernel/shaders/node_subsurface_scattering.osl
@@ -28,7 +28,9 @@ shader node_subsurface_scattering(
 {
 	if (Falloff == "Gaussian")
 		BSSRDF = Color * bssrdf_gaussian(N, Scale * Radius, TextureBlur);
-	else
+	else if (Falloff == "Cubic")
 		BSSRDF = Color * bssrdf_cubic(N, Scale * Radius, TextureBlur, Sharpness);
+	else
+		BSSRDF = Color * bssrdf_burley(N, Scale * Radius, TextureBlur, Color);
 }
 
diff --git a/intern/cycles/kernel/shaders/stdosl.h b/intern/cycles/kernel/shaders/stdosl.h
index 697a1756119..acf3ae8b1c7 100644
--- a/intern/cycles/kernel/shaders/stdosl.h
+++ b/intern/cycles/kernel/shaders/stdosl.h
@@ -510,6 +510,7 @@ closure color ambient_occlusion() BUILTIN;
 // BSSRDF
 closure color bssrdf_cubic(normal N, vector radius, float texture_blur, float sharpness) BUILTIN;
 closure color bssrdf_gaussian(normal N, vector radius, float texture_blur) BUILTIN;
+closure color bssrdf_burley(normal N, vector radius, float texture_blur, color albedo) BUILTIN;
 
 // Hair
 closure color hair_reflection(normal N, float roughnessu, float roughnessv, vector T, float offset) BUILTIN;
diff --git a/intern/cycles/kernel/svm/svm_closure.h b/intern/cycles/kernel/svm/svm_closure.h
index e3f469b1cb6..b24873d9839 100644
--- a/intern/cycles/kernel/svm/svm_closure.h
+++ b/intern/cycles/kernel/svm/svm_closure.h
@@ -442,8 +442,10 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 #  define sc_next(sc) sc = ccl_fetch_array(sd, closure, ccl_fetch(sd, num_closure))
 #endif
 		case CLOSURE_BSSRDF_CUBIC_ID:
-		case CLOSURE_BSSRDF_GAUSSIAN_ID: {
+		case CLOSURE_BSSRDF_GAUSSIAN_ID:
+		case CLOSURE_BSSRDF_BURLEY_ID: {
 			ShaderClosure *sc = ccl_fetch_array(sd, closure, ccl_fetch(sd, num_closure));
+			float3 albedo = sc->weight;
 			float3 weight = sc->weight * mix_weight;
 			float sample_weight = fabsf(average(weight));
 			
@@ -467,7 +469,7 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 					sc->sample_weight = sample_weight;
 					sc->data0 = radius.x;
 					sc->data1 = texture_blur;
-					sc->data2 = 0.0f;
+					sc->data2 = albedo.x;
 					sc->T.x = sharpness;
 #ifdef __OSL__
 					sc->prim = NULL;
@@ -484,7 +486,7 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 					sc->sample_weight = sample_weight;
 					sc->data0 = radius.y;
 					sc->data1 = texture_blur;
-					sc->data2 = 0.0f;
+					sc->data2 = albedo.y;
 					sc->T.x = sharpness;
 #ifdef __OSL__
 					sc->prim = NULL;
@@ -501,7 +503,7 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 					sc->sample_weight = sample_weight;
 					sc->data0 = radius.z;
 					sc->data1 = texture_blur;
-					sc->data2 = 0.0f;
+					sc->data2 = albedo.z;
 					sc->T.x = sharpness;
 #ifdef __OSL__
 					sc->prim = NULL;
diff --git a/intern/cycles/kernel/svm/svm_types.h b/intern/cycles/kernel/svm/svm_types.h
index 516d8aaa959..21b0cb15a4f 100644
--- a/intern/cycles/kernel/svm/svm_types.h
+++ b/intern/cycles/kernel/svm/svm_types.h
@@ -410,6 +410,7 @@ typedef enum ClosureType {
 	/* BSSRDF */
 	CLOSURE_BSSRDF_CUBIC_ID,
 	CLOSURE_BSSRDF_GAUSSIAN_ID,
+	CLOSURE_BSSRDF_BURLEY_ID,
 
 	/* Other */
 	CLOSURE_EMISSION_ID,
@@ -432,8 +433,8 @@ typedef enum ClosureType {
 #define CLOSURE_IS_BSDF_TRANSMISSION(type) (type >= CLOSURE_BSDF_TRANSMISSION_ID && type <= CLOSURE_BSDF_HAIR_TRANSMISSION_ID)
 #define CLOSURE_IS_BSDF_BSSRDF(type) (type == CLOSURE_BSDF_BSSRDF_ID)
 #define CLOSURE_IS_BSDF_ANISOTROPIC(type) (type >= CLOSURE_BSDF_MICROFACET_GGX_ANISO_ID && type <= CLOSURE_BSDF_ASHIKHMIN_SHIRLEY_ANISO_ID)
-#define CLOSURE_IS_BSDF_OR_BSSRDF(type) (type <= CLOSURE_BSSRDF_GAUSSIAN_ID)
-#define CLOSURE_IS_BSSRDF(type) (type >= CLOSURE_BSSRDF_CUBIC_ID && type <= CLOSURE_BSSRDF_GAUSSIAN_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_VOLUME(type) (type >= CLOSURE_VOLUME_ID && type <= CLOSURE_VOLUME_HENYEY_GREENSTEIN_ID)
 #define CLOSURE_IS_EMISSION(type) (type == CLOSURE_EMISSION_ID)
 #define CLOSURE_IS_HOLDOUT(type) (type == CLOSURE_HOLDOUT_ID)
diff --git a/intern/cycles/render/nodes.cpp b/intern/cycles/render/nodes.cpp
index 08c213ad716..c7b2087cf78 100644
--- a/intern/cycles/render/nodes.cpp
+++ b/intern/cycles/render/nodes.cpp
@@ -2243,6 +2243,7 @@ static ShaderEnum subsurface_falloff_init()
 
 	enm.insert("Cubic", CLOSURE_BSSRDF_CUBIC_ID);
 	enm.insert("Gaussian", CLOSURE_BSSRDF_GAUSSIAN_ID);
+	enm.insert("Burley", CLOSURE_BSSRDF_BURLEY_ID);
 
 	return enm;
 }
diff --git a/source/blender/makesdna/DNA_node_types.h b/source/blender/makesdna/DNA_node_types.h
index 3b35320c803..cf5d2dfc815 100644
--- a/source/blender/makesdna/DNA_node_types.h
+++ b/source/blender/makesdna/DNA_node_types.h
@@ -1057,6 +1057,7 @@ enum {
 #endif
 	SHD_SUBSURFACE_CUBIC			= 1,
 	SHD_SUBSURFACE_GAUSSIAN			= 2,
+	SHD_SUBSURFACE_BURLEY			= 3,
 };
 
 /* blur node */
diff --git a/source/blender/makesrna/intern/rna_nodetree.c b/source/blender/makesrna/intern/rna_nodetree.c
index f06492e4ec5..50ad7a4ceb2 100644
--- a/source/blender/makesrna/intern/rna_nodetree.c
+++ b/source/blender/makesrna/intern/rna_nodetree.c
@@ -4247,6 +4247,7 @@ static void def_sh_subsurface(StructRNA *srna)
 	static EnumPropertyItem prop_subsurface_falloff_items[] = {
 		{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"},
 		{0, NULL, 0, NULL, NULL}
 	};