Code refactor: store RGB BSSRDF in a single closure.

Previously we stored each color channel in a single closure, which was
convenient for sampling a closure and channel together. But this doesn't
work so well for algorithms where we want to render multiple color
channels together.

6 files changed, 92 insertions, 139 deletions

diff --git a/intern/cycles/kernel/closure/bssrdf.h b/intern/cycles/kernel/closure/bssrdf.h
index ec10e452148..383c168463b 100644
--- a/intern/cycles/kernel/closure/bssrdf.h
+++ b/intern/cycles/kernel/closure/bssrdf.h
@@ -22,11 +22,10 @@ CCL_NAMESPACE_BEGIN
 typedef ccl_addr_space struct Bssrdf {
 	SHADER_CLOSURE_BASE;
 
-	float radius;
+	float3 radius;
+	float3 albedo;
 	float sharpness;
-	float d;
 	float texture_blur;
-	float albedo;
 	float roughness;
 } Bssrdf;
 
@@ -200,7 +199,7 @@ ccl_device_inline float bssrdf_burley_fitting(float A)
 /* 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)
+ccl_device_inline float3 bssrdf_burley_compatible_mfp(float3 r)
 {
 	return 0.25f * M_1_PI_F * r;
 }
@@ -208,13 +207,14 @@ ccl_device_inline float bssrdf_burley_compatible_mfp(float r)
 ccl_device void bssrdf_burley_setup(Bssrdf *bssrdf)
 {
 	/* Mean free path length. */
-	const float l = bssrdf_burley_compatible_mfp(bssrdf->radius);
+	const float3 l = bssrdf_burley_compatible_mfp(bssrdf->radius);
 	/* Surface albedo. */
-	const float A = bssrdf->albedo;
-	const float s = bssrdf_burley_fitting(A);
-	const float d = l / s;
+	const float3 A = bssrdf->albedo;
+	const float3 s = make_float3(bssrdf_burley_fitting(A.x),
+                                 bssrdf_burley_fitting(A.y),
+                                 bssrdf_burley_fitting(A.z));
 
-	bssrdf->d = d;
+	bssrdf->radius = l / s;
 }
 
 ccl_device float bssrdf_burley_eval(const float d, float r)
@@ -345,7 +345,7 @@ ccl_device_inline Bssrdf *bssrdf_alloc(ShaderData *sd, float3 weight)
 
 ccl_device int bssrdf_setup(Bssrdf *bssrdf, ClosureType type)
 {
-	if(bssrdf->radius < BSSRDF_MIN_RADIUS) {
+	if(max3(bssrdf->radius) < BSSRDF_MIN_RADIUS) {
 		/* revert to diffuse BSDF if radius too small */
 		int flag;
 #ifdef __PRINCIPLED__
@@ -393,25 +393,60 @@ ccl_device int bssrdf_setup(Bssrdf *bssrdf, ClosureType type)
 ccl_device void bssrdf_sample(const ShaderClosure *sc, float xi, float *r, float *h)
 {
 	const Bssrdf *bssrdf = (const Bssrdf*)sc;
+	float radius;
 
-	if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
-		bssrdf_cubic_sample(bssrdf->radius, bssrdf->sharpness, xi, r, h);
-	else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
-		bssrdf_gaussian_sample(bssrdf->radius, xi, r, h);
-	else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID || sc->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
-		bssrdf_burley_sample(bssrdf->d, xi, r, h);
+	/* Sample color channel and reuse random number. */
+	if(xi < 1.0f/3.0f) {
+		xi *= 3.0f;
+		radius = bssrdf->radius.x;
+	}
+	else if(xi < 2.0f/3.0f) {
+		xi = (xi - 1.0f/3.0f)*3.0f;
+		radius = bssrdf->radius.y;
+	}
+	else {
+		xi = (xi - 2.0f/3.0f)*3.0f;
+		radius = bssrdf->radius.z;
+	}
+
+	/* Sample BSSRDF. */
+	if(bssrdf->type == CLOSURE_BSSRDF_CUBIC_ID) {
+		bssrdf_cubic_sample(radius, bssrdf->sharpness, xi, r, h);
+	}
+	else if(bssrdf->type == CLOSURE_BSSRDF_GAUSSIAN_ID){
+		bssrdf_gaussian_sample(radius, xi, r, h);
+	}
+	else { /*if(bssrdf->type == CLOSURE_BSSRDF_BURLEY_ID || bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
+		bssrdf_burley_sample(radius, xi, r, h);
+	}
 }
 
-ccl_device_forceinline float bssrdf_pdf(const ShaderClosure *sc, float r)
+ccl_device float bssrdf_channel_pdf(const Bssrdf *bssrdf, float radius, float r)
+{
+	if(bssrdf->type == CLOSURE_BSSRDF_CUBIC_ID) {
+		return bssrdf_cubic_pdf(radius, bssrdf->sharpness, r);
+	}
+	else if(bssrdf->type == CLOSURE_BSSRDF_GAUSSIAN_ID) {
+		return bssrdf_gaussian_pdf(radius, r);
+	}
+	else { /*if(bssrdf->type == CLOSURE_BSSRDF_BURLEY_ID || bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
+		return bssrdf_burley_pdf(radius, r);
+	}
+}
+
+ccl_device_forceinline float3 bssrdf_eval(const ShaderClosure *sc, float r)
 {
 	const Bssrdf *bssrdf = (const Bssrdf*)sc;
 
-	if(sc->type == CLOSURE_BSSRDF_CUBIC_ID)
-		return bssrdf_cubic_pdf(bssrdf->radius, bssrdf->sharpness, r);
-	else if(sc->type == CLOSURE_BSSRDF_GAUSSIAN_ID)
-		return bssrdf_gaussian_pdf(bssrdf->radius, r);
-	else /*if(sc->type == CLOSURE_BSSRDF_BURLEY_ID || sc->type == CLOSURE_BSSRDF_PRINCIPLED_ID)*/
-		return bssrdf_burley_pdf(bssrdf->d, r);
+	return make_float3(
+		bssrdf_channel_pdf(bssrdf, bssrdf->radius.x, r),
+		bssrdf_channel_pdf(bssrdf, bssrdf->radius.y, r),
+		bssrdf_channel_pdf(bssrdf, bssrdf->radius.z, r));
+}
+
+ccl_device_forceinline float bssrdf_pdf(const ShaderClosure *sc, float r)
+{
+	return average(bssrdf_eval(sc, r));
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/kernel_path_branched.h b/intern/cycles/kernel/kernel_path_branched.h
index 63fe7822e2a..fe2a7d179a4 100644
--- a/intern/cycles/kernel/kernel_path_branched.h
+++ b/intern/cycles/kernel/kernel_path_branched.h
@@ -333,7 +333,7 @@ ccl_device void kernel_branched_path_subsurface_scatter(KernelGlobals *kg,
 
 		/* set up random number generator */
 		uint lcg_state = lcg_state_init(state, 0x68bc21eb);
-		int num_samples = kernel_data.integrator.subsurface_samples;
+		int num_samples = kernel_data.integrator.subsurface_samples * 3;
 		float num_samples_inv = 1.0f/num_samples;
 		uint bssrdf_rng_hash = cmj_hash(state->rng_hash, i);
 
diff --git a/intern/cycles/kernel/kernel_subsurface.h b/intern/cycles/kernel/kernel_subsurface.h
index c79992ee99b..582a20704d3 100644
--- a/intern/cycles/kernel/kernel_subsurface.h
+++ b/intern/cycles/kernel/kernel_subsurface.h
@@ -23,11 +23,6 @@ CCL_NAMESPACE_BEGIN
  *
  */
 
-/* TODO:
- * - test using power heuristic for combing bssrdfs
- * - try to reduce one sample model variance
- */
-
 ccl_device_inline float3 subsurface_scatter_eval(ShaderData *sd,
                                                  const ShaderClosure *sc,
                                                  float disk_r,
@@ -63,12 +58,11 @@ ccl_device_inline float3 subsurface_scatter_eval(ShaderData *sd,
 			float sample_weight = (all)? 1.0f: sc->sample_weight * sample_weight_inv;
 
 			/* compute pdf */
-			float pdf = bssrdf_pdf(sc, r);
-			float disk_pdf = bssrdf_pdf(sc, disk_r);
+			float3 eval = bssrdf_eval(sc, r);
+			float pdf = bssrdf_pdf(sc, disk_r);
 
-			/* TODO power heuristic is not working correct here */
-			eval_sum += sc->weight*pdf; //*sample_weight*disk_pdf;
-			pdf_sum += sample_weight*disk_pdf; //*sample_weight*disk_pdf;
+			eval_sum += sc->weight * eval;
+			pdf_sum += sample_weight * pdf;
 		}
 	}
 
@@ -190,20 +184,20 @@ ccl_device_inline int subsurface_scatter_multi_intersect(
 	disk_N = sd->Ng;
 	make_orthonormals(disk_N, &disk_T, &disk_B);
 
-	if(disk_u < 0.5f) {
+	if(disk_v < 0.5f) {
 		pick_pdf_N = 0.5f;
 		pick_pdf_T = 0.25f;
 		pick_pdf_B = 0.25f;
-		disk_u *= 2.0f;
+		disk_v *= 2.0f;
 	}
-	else if(disk_u < 0.75f) {
+	else if(disk_v < 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;
+		disk_v = (disk_v - 0.5f)*4.0f;
 	}
 	else {
 		float3 tmp = disk_N;
@@ -212,15 +206,14 @@ ccl_device_inline int subsurface_scatter_multi_intersect(
 		pick_pdf_N = 0.25f;
 		pick_pdf_T = 0.25f;
 		pick_pdf_B = 0.5f;
-		disk_u = (disk_u - 0.75f)*4.0f;
+		disk_v = (disk_v - 0.75f)*4.0f;
 	}
 
 	/* sample point on disk */
-	float phi = M_2PI_F * disk_u;
-	float disk_r = disk_v;
-	float disk_height;
+	float phi = M_2PI_F * disk_v;
+	float disk_height, disk_r;
 
-	bssrdf_sample(sc, disk_r, &disk_r, &disk_height);
+	bssrdf_sample(sc, disk_u, &disk_r, &disk_height);
 
 	float3 disk_P = (disk_r*cosf(phi)) * disk_T + (disk_r*sinf(phi)) * disk_B;
 
@@ -359,20 +352,20 @@ ccl_device void subsurface_scatter_step(KernelGlobals *kg, ShaderData *sd, ccl_a
 	disk_N = sd->Ng;
 	make_orthonormals(disk_N, &disk_T, &disk_B);
 
-	if(disk_u < 0.5f) {
+	if(disk_v < 0.5f) {
 		pick_pdf_N = 0.5f;
 		pick_pdf_T = 0.25f;
 		pick_pdf_B = 0.25f;
-		disk_u *= 2.0f;
+		disk_v *= 2.0f;
 	}
-	else if(disk_u < 0.75f) {
+	else if(disk_v < 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;
+		disk_v = (disk_v - 0.5f)*4.0f;
 	}
 	else {
 		float3 tmp = disk_N;
@@ -381,15 +374,14 @@ ccl_device void subsurface_scatter_step(KernelGlobals *kg, ShaderData *sd, ccl_a
 		pick_pdf_N = 0.25f;
 		pick_pdf_T = 0.25f;
 		pick_pdf_B = 0.5f;
-		disk_u = (disk_u - 0.75f)*4.0f;
+		disk_v = (disk_v - 0.75f)*4.0f;
 	}
 
 	/* sample point on disk */
-	float phi = M_2PI_F * disk_u;
-	float disk_r = disk_v;
-	float disk_height;
+	float phi = M_2PI_F * disk_v;
+	float disk_height, disk_r;
 
-	bssrdf_sample(sc, disk_r, &disk_r, &disk_height);
+	bssrdf_sample(sc, disk_u, &disk_r, &disk_height);
 
 	float3 disk_P = (disk_r*cosf(phi)) * disk_T + (disk_r*sinf(phi)) * disk_B;
 
diff --git a/intern/cycles/kernel/osl/osl_bssrdf.cpp b/intern/cycles/kernel/osl/osl_bssrdf.cpp
index 27a96720c1e..b79cef244a0 100644
--- a/intern/cycles/kernel/osl/osl_bssrdf.cpp
+++ b/intern/cycles/kernel/osl/osl_bssrdf.cpp
@@ -72,36 +72,12 @@ public:
 			float texture_blur = params.texture_blur;
 
 			/* create one closure per color channel */
-			Bssrdf *bssrdf = bssrdf_alloc(sd, make_float3(weight.x, 0.0f, 0.0f));
+			Bssrdf *bssrdf = bssrdf_alloc(sd, weight);
 			if(bssrdf) {
-				bssrdf->sample_weight = sample_weight;
-				bssrdf->radius = radius.x;
+				bssrdf->sample_weight = sample_weight * 3.0f;
+				bssrdf->radius = radius;
+				bssrdf->albedo = albedo;
 				bssrdf->texture_blur = texture_blur;
-				bssrdf->albedo = albedo.x;
-				bssrdf->sharpness = sharpness;
-				bssrdf->N = params.N;
-				bssrdf->roughness = params.roughness;
-				sd->flag |= bssrdf_setup(bssrdf, (ClosureType)type);
-			}
-
-			bssrdf = bssrdf_alloc(sd, make_float3(0.0f, weight.y, 0.0f));
-			if(bssrdf) {
-				bssrdf->sample_weight = sample_weight;
-				bssrdf->radius = radius.y;
-				bssrdf->texture_blur = texture_blur;
-				bssrdf->albedo = albedo.y;
-				bssrdf->sharpness = sharpness;
-				bssrdf->N = params.N;
-				bssrdf->roughness = params.roughness;
-				sd->flag |= bssrdf_setup(bssrdf, (ClosureType)type);
-			}
-
-			bssrdf = bssrdf_alloc(sd, make_float3(0.0f, 0.0f, weight.z));
-			if(bssrdf) {
-				bssrdf->sample_weight = sample_weight;
-				bssrdf->radius = radius.z;
-				bssrdf->texture_blur = texture_blur;
-				bssrdf->albedo = albedo.z;
 				bssrdf->sharpness = sharpness;
 				bssrdf->N = params.N;
 				bssrdf->roughness = params.roughness;
diff --git a/intern/cycles/kernel/split/kernel_subsurface_scatter.h b/intern/cycles/kernel/split/kernel_subsurface_scatter.h
index 5bf7483e9a2..38dd1dc5654 100644
--- a/intern/cycles/kernel/split/kernel_subsurface_scatter.h
+++ b/intern/cycles/kernel/split/kernel_subsurface_scatter.h
@@ -54,7 +54,7 @@ ccl_device_noinline bool kernel_split_branched_path_subsurface_indirect_light_it
 			branched_state->lcg_state = lcg_state_init_addrspace(&branched_state->path_state,
 			                                                     0x68bc21eb);
 		}
-		int num_samples = kernel_data.integrator.subsurface_samples;
+		int num_samples = kernel_data.integrator.subsurface_samples * 3;
 		float num_samples_inv = 1.0f/num_samples;
 		uint bssrdf_rng_hash = cmj_hash(branched_state->path_state.rng_hash, i);
 
diff --git a/intern/cycles/kernel/svm/svm_closure.h b/intern/cycles/kernel/svm/svm_closure.h
index 47ebe4288e3..b2fc01f617e 100644
--- a/intern/cycles/kernel/svm/svm_closure.h
+++ b/intern/cycles/kernel/svm/svm_closure.h
@@ -191,40 +191,12 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 					float texture_blur = 0.0f;
 
 					/* create one closure per color channel */
-					Bssrdf *bssrdf = bssrdf_alloc(sd, make_float3(subsurf_weight.x, 0.0f, 0.0f));
+					Bssrdf *bssrdf = bssrdf_alloc(sd, subsurf_weight);
 					if(bssrdf) {
-						bssrdf->sample_weight = subsurf_sample_weight;
-						bssrdf->radius = radius.x;
+						bssrdf->sample_weight = subsurf_sample_weight * 3.0f;
+						bssrdf->radius = radius;
+						bssrdf->albedo = subsurface_color;
 						bssrdf->texture_blur = texture_blur;
-						bssrdf->albedo = subsurface_color.x;
-						bssrdf->sharpness = sharpness;
-						bssrdf->N = N;
-						bssrdf->roughness = roughness;
-
-						/* setup bsdf */
-						sd->flag |= bssrdf_setup(bssrdf, (ClosureType)CLOSURE_BSSRDF_PRINCIPLED_ID);
-					}
-
-					bssrdf = bssrdf_alloc(sd, make_float3(0.0f, subsurf_weight.y, 0.0f));
-					if(bssrdf) {
-						bssrdf->sample_weight = subsurf_sample_weight;
-						bssrdf->radius = radius.y;
-						bssrdf->texture_blur = texture_blur;
-						bssrdf->albedo = subsurface_color.y;
-						bssrdf->sharpness = sharpness;
-						bssrdf->N = N;
-						bssrdf->roughness = roughness;
-
-						/* setup bsdf */
-						sd->flag |= bssrdf_setup(bssrdf, (ClosureType)CLOSURE_BSSRDF_PRINCIPLED_ID);
-					}
-
-					bssrdf = bssrdf_alloc(sd, make_float3(0.0f, 0.0f, subsurf_weight.z));
-					if(bssrdf) {
-						bssrdf->sample_weight = subsurf_sample_weight;
-						bssrdf->radius = radius.z;
-						bssrdf->texture_blur = texture_blur;
-						bssrdf->albedo = subsurface_color.z;
 						bssrdf->sharpness = sharpness;
 						bssrdf->N = N;
 						bssrdf->roughness = roughness;
@@ -803,34 +775,12 @@ ccl_device void svm_node_closure_bsdf(KernelGlobals *kg, ShaderData *sd, float *
 				float texture_blur = param2;
 
 				/* create one closure per color channel */
-				Bssrdf *bssrdf = bssrdf_alloc(sd, make_float3(weight.x, 0.0f, 0.0f));
-				if(bssrdf) {
-					bssrdf->sample_weight = sample_weight;
-					bssrdf->radius = radius.x;
-					bssrdf->texture_blur = texture_blur;
-					bssrdf->albedo = albedo.x;
-					bssrdf->sharpness = sharpness;
-					bssrdf->N = N;
-					sd->flag |= bssrdf_setup(bssrdf, (ClosureType)type);
-				}
-
-				bssrdf = bssrdf_alloc(sd, make_float3(0.0f, weight.y, 0.0f));
-				if(bssrdf) {
-					bssrdf->sample_weight = sample_weight;
-					bssrdf->radius = radius.y;
-					bssrdf->texture_blur = texture_blur;
-					bssrdf->albedo = albedo.y;
-					bssrdf->sharpness = sharpness;
-					bssrdf->N = N;
-					sd->flag |= bssrdf_setup(bssrdf, (ClosureType)type);
-				}
-
-				bssrdf = bssrdf_alloc(sd, make_float3(0.0f, 0.0f, weight.z));
+				Bssrdf *bssrdf = bssrdf_alloc(sd, weight);
 				if(bssrdf) {
-					bssrdf->sample_weight = sample_weight;
-					bssrdf->radius = radius.z;
+					bssrdf->sample_weight = sample_weight * 3.0f;
+					bssrdf->radius = radius;
+					bssrdf->albedo = albedo;
 					bssrdf->texture_blur = texture_blur;
-					bssrdf->albedo = albedo.z;
 					bssrdf->sharpness = sharpness;
 					bssrdf->N = N;
 					sd->flag |= bssrdf_setup(bssrdf, (ClosureType)type);