Cycles Volume Render: heterogeneous (textured) volumes support.

Volumes can now have textured colors and density. There is a Volume Sampling
panel in the Render properties with these settings:

* Step size: distance between volume shader samples when rendering the volume.
  Lower values give more accurate and detailed results but also increased render
  time.
* Max steps: maximum number of steps through the volume before giving up, to
  protect from extremely long render times with big objects or small step sizes.

This is much more compute intensive than homogeneous volume, so when you are not
using a texture you should enable the Homogeneous Volume option in the material
or world for faster rendering.

One important missing feature is that Generated texture coordinates are not yet
working in volumes, and they are the default coordinates for nearly all texture
nodes. So until that works you need to plug in object texture coordinates or a
world space position.

This is work by "storm", Stuart Broadfoot, Thomas Dinges and myself.

6 files changed, 242 insertions, 45 deletions

diff --git a/intern/cycles/kernel/kernel_path.h b/intern/cycles/kernel/kernel_path.h
index 204f782e5f4..ee56dbd40d9 100644
--- a/intern/cycles/kernel/kernel_path.h
+++ b/intern/cycles/kernel/kernel_path.h
@@ -94,8 +94,9 @@ ccl_device void kernel_path_indirect(KernelGlobals *kg, RNG *rng, int sample, Ra
 #ifdef __VOLUME__
 		/* volume attenuation, emission, scatter */
 		if(state.volume_stack[0].shader != SHADER_NO_ID) {
-			ray.t = (hit)? isect.t: FLT_MAX;
-			kernel_volume_integrate(kg, &state, &ray, L, &throughput);
+			Ray volume_ray = ray;
+			volume_ray.t = (hit)? isect.t: FLT_MAX;
+			kernel_volume_integrate(kg, &state, &volume_ray, L, &throughput);
 		}
 #endif
 
@@ -462,7 +463,7 @@ ccl_device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample,
 	int num_samples = 0;
 #endif
 
-	path_state_init(kg, &state);
+	path_state_init(kg, &state, rng, sample);
 
 	/* path iteration */
 	for(;; rng_offset += PRNG_BOUNCE_NUM) {
@@ -515,8 +516,9 @@ ccl_device float4 kernel_path_integrate(KernelGlobals *kg, RNG *rng, int sample,
 #ifdef __VOLUME__
 		/* volume attenuation, emission, scatter */
 		if(state.volume_stack[0].shader != SHADER_NO_ID) {
-			ray.t = (hit)? isect.t: FLT_MAX;
-			kernel_volume_integrate(kg, &state, &ray, &L, &throughput);
+			Ray volume_ray = ray;
+			volume_ray.t = (hit)? isect.t: FLT_MAX;
+			kernel_volume_integrate(kg, &state, &volume_ray, &L, &throughput);
 		}
 #endif
 
@@ -988,7 +990,7 @@ ccl_device float4 kernel_branched_path_integrate(KernelGlobals *kg, RNG *rng, in
 	int aa_samples = 0;
 #endif
 
-	path_state_init(kg, &state);
+	path_state_init(kg, &state, rng, sample);
 
 	for(;; rng_offset += PRNG_BOUNCE_NUM) {
 		/* intersect scene */
@@ -1018,8 +1020,9 @@ ccl_device float4 kernel_branched_path_integrate(KernelGlobals *kg, RNG *rng, in
 #ifdef __VOLUME__
 		/* volume attenuation, emission, scatter */
 		if(state.volume_stack[0].shader != SHADER_NO_ID) {
-			ray.t = (hit)? isect.t: FLT_MAX;
-			kernel_volume_integrate(kg, &state, &ray, &L, &throughput);
+			Ray volume_ray = ray;
+			volume_ray.t = (hit)? isect.t: FLT_MAX;
+			kernel_volume_integrate(kg, &state, &volume_ray, &L, &throughput);
 		}
 #endif
 
diff --git a/intern/cycles/kernel/kernel_path_state.h b/intern/cycles/kernel/kernel_path_state.h
index afca28fddcf..3e0ad81c187 100644
--- a/intern/cycles/kernel/kernel_path_state.h
+++ b/intern/cycles/kernel/kernel_path_state.h
@@ -16,7 +16,7 @@
 
 CCL_NAMESPACE_BEGIN
 
-ccl_device_inline void path_state_init(KernelGlobals *kg, PathState *state)
+ccl_device_inline void path_state_init(KernelGlobals *kg, PathState *state, RNG *rng, int sample)
 {
 	state->flag = PATH_RAY_CAMERA|PATH_RAY_SINGULAR|PATH_RAY_MIS_SKIP;
 	state->bounce = 0;
@@ -26,7 +26,15 @@ ccl_device_inline void path_state_init(KernelGlobals *kg, PathState *state)
 	state->transparent_bounce = 0;
 
 #ifdef __VOLUME__
-	kernel_volume_stack_init(kg, state->volume_stack);
+	if(kernel_data.integrator.use_volumes) {
+		/* initialize volume stack with volume we are inside of */
+		kernel_volume_stack_init(kg, state->volume_stack);
+		/* seed RNG for cases where we can't use stratified samples */
+		state->rng_congruential = lcg_init(*rng + sample*0x51633e2d);
+	}
+	else {
+		state->volume_stack[0].shader = SHADER_NO_ID;
+	}
 #endif
 }
 
diff --git a/intern/cycles/kernel/kernel_random.h b/intern/cycles/kernel/kernel_random.h
index 69e7b439e1c..c435c27b7b4 100644
--- a/intern/cycles/kernel/kernel_random.h
+++ b/intern/cycles/kernel/kernel_random.h
@@ -18,8 +18,6 @@
 
 CCL_NAMESPACE_BEGIN
 
-typedef uint RNG;
-
 #ifdef __SOBOL__
 
 /* skip initial numbers that are not as well distributed, especially the
@@ -192,10 +190,6 @@ ccl_device void path_rng_end(KernelGlobals *kg, ccl_global uint *rng_state, RNG
 
 /* Linear Congruential Generator */
 
-ccl_device float path_rng(KernelGlobals *kg, RNG& rng, int sample, int dimension)
-{
-}
-
 ccl_device_inline float path_rng_1D(KernelGlobals *kg, RNG& rng, int sample, int num_samples, int dimension)
 {
 	/* implicit mod 2^32 */
diff --git a/intern/cycles/kernel/kernel_shadow.h b/intern/cycles/kernel/kernel_shadow.h
index 4bf063ee185..9b015c98c40 100644
--- a/intern/cycles/kernel/kernel_shadow.h
+++ b/intern/cycles/kernel/kernel_shadow.h
@@ -74,7 +74,7 @@ ccl_device_inline bool shadow_blocked(KernelGlobals *kg, PathState *state, Ray *
 #ifdef __VOLUME__
 					/* attenuation for last line segment towards light */
 					if(ps.volume_stack[0].shader != SHADER_NO_ID)
-						kernel_volume_get_shadow_attenuation(kg, &ps, ray, &throughput);
+						kernel_volume_shadow(kg, &ps, ray, &throughput);
 #endif
 
 					*shadow *= throughput;
@@ -89,7 +89,7 @@ ccl_device_inline bool shadow_blocked(KernelGlobals *kg, PathState *state, Ray *
 				if(ps.volume_stack[0].shader != SHADER_NO_ID) {
 					Ray segment_ray = *ray;
 					segment_ray.t = isect.t;
-					kernel_volume_get_shadow_attenuation(kg, &ps, &segment_ray, &throughput);
+					kernel_volume_shadow(kg, &ps, &segment_ray, &throughput);
 				}
 #endif
 
@@ -120,7 +120,7 @@ ccl_device_inline bool shadow_blocked(KernelGlobals *kg, PathState *state, Ray *
 #ifdef __VOLUME__
 	else if(!result && state->volume_stack[0].shader != SHADER_NO_ID) {
 		/* apply attenuation from current volume shader */
-		kernel_volume_get_shadow_attenuation(kg, state, ray, shadow);
+		kernel_volume_shadow(kg, state, ray, shadow);
 	}
 #endif
 #endif
diff --git a/intern/cycles/kernel/kernel_types.h b/intern/cycles/kernel/kernel_types.h
index 63cd2b223a3..1463ce98d15 100644
--- a/intern/cycles/kernel/kernel_types.h
+++ b/intern/cycles/kernel/kernel_types.h
@@ -153,6 +153,10 @@ CCL_NAMESPACE_BEGIN
 #error "OpenCL: mismatch between advanced shading flags in device_opencl.cpp and kernel_types.h"
 #endif
 
+/* Random Numbers */
+
+typedef uint RNG;
+
 /* Shader Evaluation */
 
 typedef enum ShaderEvalType {
@@ -511,10 +515,10 @@ enum ShaderDataFlag {
 	SD_HAS_TRANSPARENT_SHADOW = 2048,	/* has transparent shadow */
 	SD_HAS_VOLUME = 4096,				/* has volume shader */
 	SD_HAS_ONLY_VOLUME = 8192,			/* has only volume shader, no surface */
-	SD_HOMOGENEOUS_VOLUME = 16384,		/* has homogeneous volume */
+	SD_HETEROGENEOUS_VOLUME = 16384,	/* has heterogeneous volume */
 	SD_HAS_BSSRDF_BUMP = 32768,			/* bssrdf normal uses bump */
 
-	SD_SHADER_FLAGS = (SD_USE_MIS|SD_HAS_TRANSPARENT_SHADOW|SD_HAS_VOLUME|SD_HAS_ONLY_VOLUME|SD_HOMOGENEOUS_VOLUME|SD_HAS_BSSRDF_BUMP),
+	SD_SHADER_FLAGS = (SD_USE_MIS|SD_HAS_TRANSPARENT_SHADOW|SD_HAS_VOLUME|SD_HAS_ONLY_VOLUME|SD_HETEROGENEOUS_VOLUME|SD_HAS_BSSRDF_BUMP),
 
 	/* object flags */
 	SD_HOLDOUT_MASK = 65536,			/* holdout for camera rays */
@@ -625,6 +629,7 @@ typedef struct PathState {
 	int transparent_bounce;
 
 #ifdef __VOLUME__
+	RNG rng_congruential;
 	VolumeStack volume_stack[VOLUME_STACK_SIZE];
 #endif
 } PathState;
@@ -806,6 +811,9 @@ typedef struct KernelIntegrator {
 
 	/* volume render */
 	int use_volumes;
+	int volume_max_steps;
+	float volume_step_size;
+	int pad1, pad2;
 } KernelIntegrator;
 
 typedef struct KernelBVH {
diff --git a/intern/cycles/kernel/kernel_volume.h b/intern/cycles/kernel/kernel_volume.h
index ee9e28999ed..d0eb6dda3ed 100644
--- a/intern/cycles/kernel/kernel_volume.h
+++ b/intern/cycles/kernel/kernel_volume.h
@@ -88,60 +88,130 @@ ccl_device float3 volume_color_attenuation(float3 sigma, float t)
 	return make_float3(expf(-sigma.x * t), expf(-sigma.y * t), expf(-sigma.z * t));
 }
 
-/* Volumetric Shadows */
+ccl_device bool volume_stack_is_heterogeneous(KernelGlobals *kg, VolumeStack *stack)
+{
+	for(int i = 0; stack[i].shader != SHADER_NO_ID; i++) {
+		int shader_flag = kernel_tex_fetch(__shader_flag, (stack[i].shader & SHADER_MASK)*2);
 
-/* get the volume attenuation over line segment defined by ray, with the
- * assumption that there are no surfaces blocking light between the endpoints */
-ccl_device void kernel_volume_get_shadow_attenuation(KernelGlobals *kg, PathState *state, Ray *ray, float3 *throughput)
+		if(shader_flag & SD_HETEROGENEOUS_VOLUME)
+			return true;
+	}
+
+	return false;
+}
+
+/* Volumetric Shadows
+ *
+ * These functions are used to attenuate shadow rays to lights. Both absorption
+ * and scattering will block light, represented by the extinction coefficient. */
+
+/* homogenous volume: assume shader evaluation at the starts gives
+ * the extinction coefficient for the entire line segment */
+ccl_device void kernel_volume_shadow_homogeneous(KernelGlobals *kg, PathState *state, Ray *ray, ShaderData *sd, float3 *throughput)
 {
-	ShaderData sd;
-	shader_setup_from_volume(kg, &sd, ray, state->bounce);
+	ShaderContext ctx = SHADER_CONTEXT_SHADOW;
+	int path_flag = PATH_RAY_SHADOW;
+	float3 sigma_t;
+
+	if(volume_shader_extinction_sample(kg, sd, state->volume_stack, path_flag, ctx, ray->P, &sigma_t))
+		*throughput *= volume_color_attenuation(sigma_t, ray->t);
+}
 
+/* heterogeneous volume: integrate stepping through the volume until we
+ * reach the end, get absorbed entirely, or run out of iterations */
+ccl_device void kernel_volume_shadow_heterogeneous(KernelGlobals *kg, PathState *state, Ray *ray, ShaderData *sd, float3 *throughput)
+{
 	ShaderContext ctx = SHADER_CONTEXT_SHADOW;
 	int path_flag = PATH_RAY_SHADOW;
+	float3 tp = *throughput;
+	const float tp_eps = 1e-10f; /* todo: this is likely not the right value */
+
+	/* prepare for stepping */
+	int max_steps = kernel_data.integrator.volume_max_steps;
+	float step = kernel_data.integrator.volume_step_size;
+	float random_jitter_offset = lcg_step_float(&state->rng_congruential) * step;
+
+	/* compute extinction at the start */
+	float t = 0.0f;
+	float3 P = ray->P;
 	float3 sigma_t;
 
-	/* homogenous volume: assume shader evaluation at the starts gives
-	 * the extinction coefficient for the entire line segment */
-	if(!volume_shader_extinction_sample(kg, &sd, state->volume_stack, path_flag, ctx, ray->P, &sigma_t))
-		return;
+	if(!volume_shader_extinction_sample(kg, sd, state->volume_stack, path_flag, ctx, P, &sigma_t))
+		sigma_t = make_float3(0.0f, 0.0f, 0.0f);
 
-	*throughput *= volume_color_attenuation(sigma_t, ray->t);
-}
+	for(int i = 0; i < max_steps; i++) {
+		/* advance to new position */
+		float new_t = min(ray->t, t + random_jitter_offset + i * step);
+		float3 new_P = ray->P + ray->D * new_t;
+		float3 new_sigma_t;
 
-/* Volumetric Path */
+		/* compute attenuation over segment */
+		if(volume_shader_extinction_sample(kg, sd, state->volume_stack, path_flag, ctx, new_P, &new_sigma_t)) {
+			/* todo: we could avoid computing expf() for each step by summing,
+			 * because exp(a)*exp(b) = exp(a+b), but we still want a quick
+			 * tp_eps check too */
+			tp *= volume_color_attenuation(0.5f*(sigma_t + new_sigma_t), new_t - t);
 
-/* get the volume attenuation and emission over line segment defined by
- * ray, with the assumption that there are no surfaces blocking light
- * between the endpoints */
-ccl_device VolumeIntegrateResult kernel_volume_integrate(KernelGlobals *kg, PathState *state, Ray *ray, PathRadiance *L, float3 *throughput)
+			/* stop if nearly all light blocked */
+			if(tp.x < tp_eps && tp.y < tp_eps && tp.z < tp_eps)
+				break;
+
+			sigma_t = new_sigma_t;
+		}
+		else {
+			/* skip empty space */
+			sigma_t = make_float3(0.0f, 0.0f, 0.0f);
+		}
+
+		/* stop if at the end of the volume */
+		t = new_t;
+		if(t == ray->t)
+			break;
+	}
+
+	*throughput = tp;
+}
+
+/* get the volume attenuation over line segment defined by ray, with the
+ * assumption that there are no surfaces blocking light between the endpoints */
+ccl_device void kernel_volume_shadow(KernelGlobals *kg, PathState *state, Ray *ray, float3 *throughput)
 {
 	ShaderData sd;
 	shader_setup_from_volume(kg, &sd, ray, state->bounce);
 
+	if(volume_stack_is_heterogeneous(kg, state->volume_stack))
+		kernel_volume_shadow_heterogeneous(kg, state, ray, &sd, throughput);
+	else
+		kernel_volume_shadow_homogeneous(kg, state, ray, &sd, throughput);
+}
+
+/* Volumetric Path */
+
+/* homogenous volume: assume shader evaluation at the starts gives
+ * the volume shading coefficient for the entire line segment */
+ccl_device VolumeIntegrateResult kernel_volume_integrate_homogeneous(KernelGlobals *kg, PathState *state, Ray *ray, ShaderData *sd, PathRadiance *L, float3 *throughput)
+{
 	ShaderContext ctx = SHADER_CONTEXT_VOLUME;
 	int path_flag = PATH_RAY_SHADOW;
 	VolumeShaderCoefficients coeff;
 
-	/* homogenous volume: assume shader evaluation at the starts gives
-	 * the extinction coefficient for the entire line segment */
-	if(!volume_shader_sample(kg, &sd, state->volume_stack, path_flag, ctx, ray->P, &coeff))
+	if(!volume_shader_sample(kg, sd, state->volume_stack, path_flag, ctx, ray->P, &coeff))
 		return VOLUME_PATH_MISSED;
 
 	/* todo: in principle the SD_EMISSION, SD_ABSORPTION and SD_SCATTER flags
 	 * should ensure that one of the components is > 0 and so no division by
 	 * zero occurs, however this needs to be double-checked and tested */
 	
-	int closure_flag = sd.flag;
+	int closure_flag = sd->flag;
 	float t = ray->t;
 
-	/* compute attenuation from absorption (+ scattering for now) */
+	/* compute attenuation from absorption */
 	float3 attenuation;
 
 	if(closure_flag & SD_ABSORPTION)
 		attenuation = volume_color_attenuation(coeff.sigma_a, t);
 
-	/* integrate emission attenuated by absorption 
+	/* integrate emission attenuated by absorption
 	 * integral E * exp(-sigma_a * t) from 0 to t = E * (1 - exp(-sigma_a * t))/sigma_a
 	 * this goes to E * t as sigma_a goes to zero
 	 *
@@ -169,7 +239,121 @@ ccl_device VolumeIntegrateResult kernel_volume_integrate(KernelGlobals *kg, Path
 	return VOLUME_PATH_ATTENUATED;
 }
 
-/* Volume Stack */
+/* heterogeneous volume: integrate stepping through the volume until we
+ * reach the end, get absorbed entirely, or run out of iterations */
+ccl_device VolumeIntegrateResult kernel_volume_integrate_heterogeneous(KernelGlobals *kg, PathState *state, Ray *ray, ShaderData *sd, PathRadiance *L, float3 *throughput)
+{
+	ShaderContext ctx = SHADER_CONTEXT_VOLUME;
+	int path_flag = PATH_RAY_SHADOW;
+	VolumeShaderCoefficients coeff;
+	float3 tp = *throughput;
+	const float tp_eps = 1e-10f; /* todo: this is likely not the right value */
+
+	/* prepare for stepping */
+	int max_steps = kernel_data.integrator.volume_max_steps;
+	float step = kernel_data.integrator.volume_step_size;
+	float random_jitter_offset = lcg_step_float(&state->rng_congruential) * step;
+
+	/* compute coefficients at the start */
+	float t = 0.0f;
+	float3 P = ray->P;
+
+	if(!volume_shader_sample(kg, sd, state->volume_stack, path_flag, ctx, P, &coeff)) {
+		coeff.sigma_a = make_float3(0.0f, 0.0f, 0.0f);
+		coeff.sigma_s = make_float3(0.0f, 0.0f, 0.0f);
+		coeff.emission = make_float3(0.0f, 0.0f, 0.0f);
+	}
+
+	for(int i = 0; i < max_steps; i++) {
+		/* advance to new position */
+		float new_t = min(ray->t, t + random_jitter_offset + i * step);
+		float3 new_P = ray->P + ray->D * new_t;
+		VolumeShaderCoefficients new_coeff;
+
+		/* compute attenuation over segment */
+		if(volume_shader_sample(kg, sd, state->volume_stack, path_flag, ctx, new_P, &new_coeff)) {
+			int closure_flag = sd->flag;
+			float dt = new_t - t;
+
+			/* compute attenuation from absorption */
+			float3 attenuation, sigma_a;
+
+			if(closure_flag & SD_ABSORPTION) {
+				/* todo: we could avoid computing expf() for each step by summing,
+				 * because exp(a)*exp(b) = exp(a+b), but we still want a quick
+				 * tp_eps check too */
+				sigma_a = 0.5f*(coeff.sigma_a + new_coeff.sigma_a);
+				attenuation = volume_color_attenuation(sigma_a, dt);
+			}
+
+			/* integrate emission attenuated by absorption 
+			 * integral E * exp(-sigma_a * t) from 0 to t = E * (1 - exp(-sigma_a * t))/sigma_a
+			 * this goes to E * t as sigma_a goes to zero
+			 *
+			 * todo: we should use an epsilon to avoid precision issues near zero sigma_a */
+			if(closure_flag & SD_EMISSION) {
+				float3 emission = 0.5f*(coeff.emission + new_coeff.emission);
+
+				if(closure_flag & SD_ABSORPTION) {
+					emission.x *= (sigma_a.x > 0.0f)? (1.0f - attenuation.x)/sigma_a.x: dt;
+					emission.y *= (sigma_a.y > 0.0f)? (1.0f - attenuation.y)/sigma_a.y: dt;
+					emission.z *= (sigma_a.z > 0.0f)? (1.0f - attenuation.z)/sigma_a.z: dt;
+				}
+				else
+					emission *= t;
+
+				path_radiance_accum_emission(L, tp, emission, state->bounce);
+			}
+
+			/* modify throughput */
+			if(closure_flag & SD_ABSORPTION) {
+				tp *= attenuation;
+
+				/* stop if nearly all light blocked */
+				if(tp.x < tp_eps && tp.y < tp_eps && tp.z < tp_eps) {
+					tp = make_float3(0.0f, 0.0f, 0.0f);
+					break;
+				}
+			}
+
+			coeff = new_coeff;
+		}
+		else {
+			/* skip empty space */
+			coeff.sigma_a = make_float3(0.0f, 0.0f, 0.0f);
+			coeff.sigma_s = make_float3(0.0f, 0.0f, 0.0f);
+			coeff.emission = make_float3(0.0f, 0.0f, 0.0f);
+		}
+
+		/* stop if at the end of the volume */
+		t = new_t;
+		if(t == ray->t)
+			break;
+	}
+
+	*throughput = tp;
+
+	return VOLUME_PATH_ATTENUATED;
+}
+
+/* get the volume attenuation and emission over line segment defined by
+ * ray, with the assumption that there are no surfaces blocking light
+ * between the endpoints */
+ccl_device VolumeIntegrateResult kernel_volume_integrate(KernelGlobals *kg, PathState *state, Ray *ray, PathRadiance *L, float3 *throughput)
+{
+	ShaderData sd;
+	shader_setup_from_volume(kg, &sd, ray, state->bounce);
+
+	if(volume_stack_is_heterogeneous(kg, state->volume_stack))
+		return kernel_volume_integrate_heterogeneous(kg, state, ray, &sd, L, throughput);
+	else
+		return kernel_volume_integrate_homogeneous(kg, state, ray, &sd, L, throughput);
+}
+
+/* Volume Stack
+ *
+ * This is an array of object/shared ID's that the current segment of the path
+ * is inside of. */
 
 ccl_device void kernel_volume_stack_init(KernelGlobals *kg, VolumeStack *stack)
 {