Cycles: equi-angular sampling for homogeneous volumes

This adds an option in the Volume Sampling panel, which helps rendering lamps
inside or near volumes with less noise. It can also increase noise though and
needs improvements to support MIS and heterogeneous volumes, but since it's
useful in some cases already (especially world volumes) it's there now.

Based on the code in the old branch by Stuart, with modifications by Thomas
and Brecht.

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

7 files changed, 87 insertions, 14 deletions

diff --git a/intern/cycles/blender/addon/properties.py b/intern/cycles/blender/addon/properties.py
index c80e8a3250c..3920510e38f 100644
--- a/intern/cycles/blender/addon/properties.py
+++ b/intern/cycles/blender/addon/properties.py
@@ -107,6 +107,11 @@ enum_integrator = (
     ('BRANCHED_PATH', "Branched Path Tracing", "Path tracing integrator that branches on the first bounce, giving more control over the number of light and material samples"),
     ('PATH', "Path Tracing", "Pure path tracing integrator"),
     )
+    
+enum_volume_homogeneous_sampling = (
+    ('DISTANCE', "Distance", "Use Distance Sampling"),
+    ('EQUI_ANGULAR', "Equi-angular", "Use Equi-angular Sampling"),
+    )
 
 
 class CyclesRenderSettings(bpy.types.PropertyGroup):
@@ -140,6 +145,13 @@ class CyclesRenderSettings(bpy.types.PropertyGroup):
                 items=enum_integrator,
                 default='PATH',
                 )
+                
+        cls.volume_homogeneous_sampling = EnumProperty(
+                name="Homogeneous Sampling",
+                description="Sampling method to use for homogeneous volumes",
+                items=enum_volume_homogeneous_sampling,
+                default='DISTANCE',
+                )
 
         cls.use_square_samples = BoolProperty(
                 name="Square Samples",
diff --git a/intern/cycles/blender/addon/ui.py b/intern/cycles/blender/addon/ui.py
index c0ce80426c0..cb50db23767 100644
--- a/intern/cycles/blender/addon/ui.py
+++ b/intern/cycles/blender/addon/ui.py
@@ -167,6 +167,10 @@ class CyclesRender_PT_volume_sampling(CyclesButtonsPanel, Panel):
         scene = context.scene
         cscene = scene.cycles
 
+        layout.prop(cscene, "volume_homogeneous_sampling", text="Homogeneous")
+
+        layout.label("Heterogeneous:")
+
         split = layout.split()
         split.prop(cscene, "volume_step_size")
         split.prop(cscene, "volume_max_steps")
diff --git a/intern/cycles/blender/blender_sync.cpp b/intern/cycles/blender/blender_sync.cpp
index 8e2197a2aa6..1d507ed9b1d 100644
--- a/intern/cycles/blender/blender_sync.cpp
+++ b/intern/cycles/blender/blender_sync.cpp
@@ -172,6 +172,7 @@ void BlenderSync::sync_integrator()
 	integrator->transparent_min_bounce = get_int(cscene, "transparent_min_bounces");
 	integrator->transparent_shadows = get_boolean(cscene, "use_transparent_shadows");
 
+	integrator->volume_homogeneous_sampling = RNA_enum_get(&cscene, "volume_homogeneous_sampling");
 	integrator->volume_max_steps = get_int(cscene, "volume_max_steps");
 	integrator->volume_step_size = get_float(cscene, "volume_step_size");
 
diff --git a/intern/cycles/kernel/kernel_types.h b/intern/cycles/kernel/kernel_types.h
index 5ee25a6cb98..8d7adb2b6d7 100644
--- a/intern/cycles/kernel/kernel_types.h
+++ b/intern/cycles/kernel/kernel_types.h
@@ -824,7 +824,6 @@ typedef struct KernelIntegrator {
 	/* clamp */
 	float sample_clamp_direct;
 	float sample_clamp_indirect;
-	float pad1, pad2, pad3;
 
 	/* branched path */
 	int branched;
@@ -843,10 +842,12 @@ typedef struct KernelIntegrator {
 	int sampling_pattern;
 
 	/* volume render */
+	int volume_homogeneous_sampling;
 	int use_volumes;
 	int volume_max_steps;
 	float volume_step_size;
 	int volume_samples;
+	int pad1, pad2;
 } KernelIntegrator;
 
 typedef struct KernelBVH {
diff --git a/intern/cycles/kernel/kernel_volume.h b/intern/cycles/kernel/kernel_volume.h
index dc2ddf1098e..4d058763f22 100644
--- a/intern/cycles/kernel/kernel_volume.h
+++ b/intern/cycles/kernel/kernel_volume.h
@@ -228,21 +228,72 @@ ccl_device VolumeIntegrateResult kernel_volume_integrate_homogeneous(KernelGloba
 		else
 			sample_sigma_t = sigma_t.z;
 
-		/* xi is [0, 1[ so log(0) should never happen, division by zero is
-		 * avoided because sample_sigma_t > 0 when SD_SCATTER is set */
-		float xi = path_state_rng_1D(kg, rng, state, PRNG_SCATTER_DISTANCE);
-		float sample_t = min(t, -logf(1.0f - xi)/sample_sigma_t);
-
-		transmittance = volume_color_attenuation(sigma_t, sample_t);
-
-		if(sample_t < t) {
-			float pdf = dot(sigma_t, transmittance);
-			new_tp = *throughput * coeff.sigma_s * transmittance * (3.0f / pdf);
-			t = sample_t;
+		/* distance sampling */
+		if(kernel_data.integrator.volume_homogeneous_sampling == 0 || !kernel_data.integrator.num_all_lights) { 
+			/* xi is [0, 1[ so log(0) should never happen, division by zero is
+			 * avoided because sample_sigma_t > 0 when SD_SCATTER is set */
+			float xi = path_state_rng_1D(kg, rng, state, PRNG_SCATTER_DISTANCE);
+			float sample_t = min(t, -logf(1.0f - xi)/sample_sigma_t);
+
+			transmittance = volume_color_attenuation(sigma_t, sample_t);
+
+			if(sample_t < t) {
+				float pdf = dot(sigma_t, transmittance);
+				new_tp = *throughput * coeff.sigma_s * transmittance * (3.0f / pdf);
+				t = sample_t;
+			}
+			else {
+				float pdf = (transmittance.x + transmittance.y + transmittance.z);
+				new_tp = *throughput * transmittance * (3.0f / pdf);
+			}
 		}
+		/* equi-angular sampling */
 		else {
-			float pdf = (transmittance.x + transmittance.y + transmittance.z);
-			new_tp = *throughput * transmittance * (3.0f / pdf);
+			/* decide if we are going to scatter or not, based on sigma_t. this
+			 * is not ideal, instead we should perhaps split the path here and
+			 * do both, and at least add multiple importance sampling */
+			float xi = path_state_rng_1D(kg, rng, state, PRNG_SCATTER_DISTANCE);
+			float sample_transmittance = expf(-sample_sigma_t * t);
+
+			if(xi < sample_transmittance) {
+				/* no scattering */
+				float3 transmittance = volume_color_attenuation(sigma_t, t);
+				float pdf = (transmittance.x + transmittance.y + transmittance.z);
+				new_tp = *throughput * transmittance * (3.0f / pdf);
+			}
+			else {
+				/* rescale random number so we can reuse it */
+				xi = (xi - sample_transmittance)/(1.0f - sample_transmittance);
+
+				/* equi-angular scattering somewhere on segment 0..t */
+				/* see "Importance Sampling Techniques for Path Tracing in Participating Media" */
+
+				/* light RNGs */
+				float light_t = path_state_rng_1D(kg, rng, state, PRNG_LIGHT);
+				float light_u, light_v;
+				path_state_rng_2D(kg, rng, state, PRNG_LIGHT_U, &light_u, &light_v);
+
+				/* light sample */
+				LightSample ls;
+				light_sample(kg, light_t, light_u, light_v, ray->time, ray->P, &ls);
+				if(ls.pdf == 0.0f)
+					return VOLUME_PATH_MISSED;
+
+				/* sampling */
+				float delta = dot((ls.P - ray->P) , ray->D);
+				float D = sqrtf(len_squared(ls.P - ray->P) - delta * delta);
+				float theta_a = -atan2f(delta, D);
+				float theta_b = atan2f(t - delta, D);
+				float t_ = D * tan((xi * theta_b) + (1 - xi) * theta_a);
+
+				float pdf = D / ((theta_b - theta_a) * (D * D + t_ * t_));
+				float sample_t = min(t, delta + t_);
+
+				transmittance = volume_color_attenuation(sigma_t, sample_t);
+
+				new_tp = *throughput * coeff.sigma_s * transmittance / ((1.0f - sample_transmittance) * pdf);
+				t = sample_t;
+			}
 		}
 	}
 	else if(closure_flag & SD_ABSORPTION) {
diff --git a/intern/cycles/render/integrator.cpp b/intern/cycles/render/integrator.cpp
index f48e04f31e1..5b26440594a 100644
--- a/intern/cycles/render/integrator.cpp
+++ b/intern/cycles/render/integrator.cpp
@@ -41,6 +41,7 @@ Integrator::Integrator()
 	transparent_probalistic = true;
 	transparent_shadows = false;
 
+	volume_homogeneous_sampling = 0;
 	volume_max_steps = 1024;
 	volume_step_size = 0.1;
 
@@ -104,6 +105,7 @@ void Integrator::device_update(Device *device, DeviceScene *dscene, Scene *scene
 
 	kintegrator->transparent_shadows = transparent_shadows;
 
+	kintegrator->volume_homogeneous_sampling = volume_homogeneous_sampling;
 	kintegrator->volume_max_steps = volume_max_steps;
 	kintegrator->volume_step_size = volume_step_size;
 
@@ -176,6 +178,7 @@ bool Integrator::modified(const Integrator& integrator)
 		transparent_max_bounce == integrator.transparent_max_bounce &&
 		transparent_probalistic == integrator.transparent_probalistic &&
 		transparent_shadows == integrator.transparent_shadows &&
+		volume_homogeneous_sampling == integrator.volume_homogeneous_sampling &&
 		volume_max_steps == integrator.volume_max_steps &&
 		volume_step_size == integrator.volume_step_size &&
 		no_caustics == integrator.no_caustics &&
diff --git a/intern/cycles/render/integrator.h b/intern/cycles/render/integrator.h
index 573b258af60..4a8240c3941 100644
--- a/intern/cycles/render/integrator.h
+++ b/intern/cycles/render/integrator.h
@@ -41,6 +41,7 @@ public:
 	bool transparent_probalistic;
 	bool transparent_shadows;
 
+	int volume_homogeneous_sampling;
 	int volume_max_steps;
 	float volume_step_size;