Rework of volume shading

After code review and experimentation, this commit makes some changes to the way that volumes are shaded. Previously, there were problems with the 'scattering' component, in that it wasn't physically correct - it didn't conserve energy and was just acting as a brightness multiplier. This has been changed to be more correct, so that as the light is scattered out of the volume, there is less remaining to penetrate through.

Since this behaviour is very similar to absorption but more useful, absorption has been removed and has been replaced by a 'transmission colour' - controlling the colour of light penetrating through the volume after it has been scattered/absorbed. As well as this, there's now 'reflection', a non-physically correct RGB multiplier for out-scattered light. This is handy for tweaking the overall colour of the volume, without having to worry about wavelength dependent absorption, and its effects on transmitted light. Now at least, even though there is the ability to tweak things non-physically, volume shading is physically based by default, and has a better combination of correctness and ease of use.

There's more detailed information and example images here:
http://wiki.blender.org/index.php/User:Broken/VolumeRendering

Also did some tweaks/optimisation:
* Removed shading step size (was a bit annoying, if it comes back, it will be in a different form)
* Removed phase function options, now just one asymmetry slider controls the range between back-scattering, isotropic scattering, and forward scattering. (note, more extreme values gives artifacts with light cache, will fix...)
* Disabled the extra 'bounce lights' from the preview render for volumes, speeds updates significantly
* Enabled voxeldata texture in preview render
* Fixed volume shadows (they were too dark, fixed by avoiding using the shadfac/AddAlphaLight stuff)

More revisions to come later...

8 files changed, 227 insertions, 163 deletions

diff --git a/source/blender/render/intern/include/volume_precache.h b/source/blender/render/intern/include/volume_precache.h
index 9d87a219c82..368d60ee7f0 100644
--- a/source/blender/render/intern/include/volume_precache.h
+++ b/source/blender/render/intern/include/volume_precache.h
@@ -30,4 +30,4 @@ void volume_precache(Render *re);
 void free_volume_precache(Render *re);
 int point_inside_volume_objectinstance(ObjectInstanceRen *obi, float *co);
 
-#define VOL_MS_TIMESTEP	0.1f
\ No newline at end of file
+#define VOL_MS_TIMESTEP	0.1f
diff --git a/source/blender/render/intern/include/volumetric.h b/source/blender/render/intern/include/volumetric.h
index 026b4840ea3..97e7e022fa0 100644
--- a/source/blender/render/intern/include/volumetric.h
+++ b/source/blender/render/intern/include/volumetric.h
@@ -26,9 +26,8 @@
  * ***** END GPL LICENSE BLOCK *****
  */
 
-float vol_get_stepsize(struct ShadeInput *shi, int context);
 float vol_get_density(struct ShadeInput *shi, float *co);
-void vol_get_scattering(ShadeInput *shi, float *scatter_col, float *co, float stepsize, float density);
+void vol_get_scattering(ShadeInput *shi, float *scatter_col, float *co_);
 
 void shade_volume_outside(ShadeInput *shi, ShadeResult *shr);
 void shade_volume_inside(ShadeInput *shi, ShadeResult *shr);
diff --git a/source/blender/render/intern/source/convertblender.c b/source/blender/render/intern/source/convertblender.c
index 48a7c003c9f..b3784f26048 100644
--- a/source/blender/render/intern/source/convertblender.c
+++ b/source/blender/render/intern/source/convertblender.c
@@ -889,6 +889,28 @@ static void free_mesh_orco_hash(Render *re)
 	}
 }
 
+static void check_material_mapto(Material *ma)
+{
+	int a;
+	ma->mapto_textured = 0;
+	
+	/* cache which inputs are actually textured.
+	 * this can avoid a bit of time spent iterating through all the texture slots, map inputs and map tos
+	 * every time a property which may or may not be textured is accessed */
+	
+	for(a=0; a<MAX_MTEX; a++) {
+		if(ma->mtex[a] && ma->mtex[a]->tex) {
+			/* currently used only in volume render, so we'll check for those flags */
+			if(ma->mtex[a]->mapto & MAP_DENSITY) ma->mapto_textured |= MAP_DENSITY;
+			if(ma->mtex[a]->mapto & MAP_EMISSION) ma->mapto_textured |= MAP_EMISSION;
+			if(ma->mtex[a]->mapto & MAP_EMISSION_COL) ma->mapto_textured |= MAP_EMISSION_COL;
+			if(ma->mtex[a]->mapto & MAP_SCATTERING) ma->mapto_textured |= MAP_SCATTERING;
+			if(ma->mtex[a]->mapto & MAP_TRANSMISSION_COL) ma->mapto_textured |= MAP_TRANSMISSION_COL;
+			if(ma->mtex[a]->mapto & MAP_REFLECTION) ma->mapto_textured |= MAP_REFLECTION;
+			if(ma->mtex[a]->mapto & MAP_REFLECTION_COL) ma->mapto_textured |= MAP_REFLECTION_COL;
+		}
+	}
+}
 static void flag_render_node_material(Render *re, bNodeTree *ntree)
 {
 	bNode *node;
@@ -930,6 +952,8 @@ static Material *give_render_material(Render *re, Object *ob, int nr)
 	if(ma->nodetree && ma->use_nodes)
 		flag_render_node_material(re, ma->nodetree);
 	
+	check_material_mapto(ma);
+	
 	return ma;
 }
 
diff --git a/source/blender/render/intern/source/rayshade.c b/source/blender/render/intern/source/rayshade.c
index c1d0c943ca9..cce99d64b39 100644
--- a/source/blender/render/intern/source/rayshade.c
+++ b/source/blender/render/intern/source/rayshade.c
@@ -1343,7 +1343,8 @@ static void ray_trace_shadow_tra(Isect *is, ShadeInput *origshi, int depth, int
 			/* mix colors based on shadfac (rgb + amount of light factor) */
 			addAlphaLight(is->col, shr.diff, shr.alpha, d*shi.mat->filter);
 		} else if (shi.mat->material_type == MA_TYPE_VOLUME) {
-			addAlphaLight(is->col, shr.combined, shr.alpha, 1.0f);
+			QUATCOPY(is->col, shr.combined);
+			is->col[3] = 1.f;
 		}
 		
 		if(depth>0 && is->col[3]>0.0f) {
diff --git a/source/blender/render/intern/source/texture.c b/source/blender/render/intern/source/texture.c
index 4830dab7d37..02b7f6aad86 100644
--- a/source/blender/render/intern/source/texture.c
+++ b/source/blender/render/intern/source/texture.c
@@ -2385,7 +2385,7 @@ void do_volume_tex(ShadeInput *shi, float *xyz, int mapto_flag, float *col, floa
 			}
 			
 			
-			if((mapto_flag & (MAP_EMISSION_COL+MAP_ABSORPTION_COL)) && (mtex->mapto & (MAP_EMISSION_COL+MAP_ABSORPTION_COL))) {
+			if((mapto_flag & (MAP_EMISSION_COL+MAP_TRANSMISSION_COL+MAP_REFLECTION_COL)) && (mtex->mapto & (MAP_EMISSION_COL+MAP_TRANSMISSION_COL+MAP_REFLECTION_COL))) {
 				float tcol[3], colfac;
 				
 				/* stencil maps on the texture control slider, not texture intensity value */
@@ -2413,8 +2413,11 @@ void do_volume_tex(ShadeInput *shi, float *xyz, int mapto_flag, float *col, floa
 					texture_rgb_blend(col, tcol, col, texres.tin, colfac, mtex->blendtype);
 				}
 				
-				/* MAP_COLMIR is abused for absorption colour at the moment */
-				if((mapto_flag & MAP_ABSORPTION_COL) && (mtex->mapto & MAP_ABSORPTION_COL)) {
+				if((mapto_flag & MAP_REFLECTION_COL) && (mtex->mapto & MAP_REFLECTION_COL)) {
+					texture_rgb_blend(col, tcol, col, texres.tin, colfac, mtex->blendtype);
+				}
+				
+				if((mapto_flag & MAP_TRANSMISSION_COL) && (mtex->mapto & MAP_TRANSMISSION_COL)) {
 					texture_rgb_blend(col, tcol, col, texres.tin, colfac, mtex->blendtype);
 				}
 			}
@@ -2443,14 +2446,14 @@ void do_volume_tex(ShadeInput *shi, float *xyz, int mapto_flag, float *col, floa
 					*val = texture_value_blend(mtex->def_var, *val, texres.tin, varfac, mtex->blendtype, flip);
 					CLAMP(*val, 0.0, 1.0);
 				}
-				if((mapto_flag & MAP_ABSORPTION) && (mtex->mapto & MAP_ABSORPTION)) {
-					int flip= mtex->maptoneg & MAP_ABSORPTION;
+				if((mapto_flag & MAP_SCATTERING) && (mtex->mapto & MAP_SCATTERING)) {
+					int flip= mtex->maptoneg & MAP_SCATTERING;
 					
 					*val = texture_value_blend(mtex->def_var, *val, texres.tin, varfac, mtex->blendtype, flip);
 					CLAMP(*val, 0.0, 1.0);
 				}
-				if((mapto_flag & MAP_SCATTERING) && (mtex->mapto & MAP_SCATTERING)) {
-					int flip= mtex->maptoneg & MAP_SCATTERING;
+				if((mapto_flag & MAP_REFLECTION) && (mtex->mapto & MAP_REFLECTION)) {
+					int flip= mtex->maptoneg & MAP_REFLECTION;
 					
 					*val = texture_value_blend(mtex->def_var, *val, texres.tin, varfac, mtex->blendtype, flip);
 					CLAMP(*val, 0.0, 1.0);
diff --git a/source/blender/render/intern/source/volume_precache.c b/source/blender/render/intern/source/volume_precache.c
index 15d8643fea4..7ecaf83ae27 100644
--- a/source/blender/render/intern/source/volume_precache.c
+++ b/source/blender/render/intern/source/volume_precache.c
@@ -185,9 +185,9 @@ static float get_avg_surrounds(float *cache, int *res, int xx, int yy, int zz)
 		}
 	}
 	
-	tot /= added;
+	if (added > 0) tot /= added;
 	
-	return ((added>0)?tot:0.0f);
+	return tot;
 }
 
 /* function to filter the edges of the light cache, where there was no volume originally.
@@ -202,17 +202,54 @@ static void lightcache_filter(VolumePrecache *vp)
 		for (y=0; y < vp->res[1]; y++) {
 			for (x=0; x < vp->res[0]; x++) {
 				/* trigger for outside mesh */
-				if (vp->data_r[ V_I(x, y, z, vp->res) ] < -0.5f)
+				if (vp->data_r[ V_I(x, y, z, vp->res) ] < -0.f)
 					vp->data_r[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_r, vp->res, x, y, z);
-				if (vp->data_g[ V_I(x, y, z, vp->res) ] < -0.5f)
+				if (vp->data_g[ V_I(x, y, z, vp->res) ] < -0.f)
 					vp->data_g[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_g, vp->res, x, y, z);
-				if (vp->data_b[ V_I(x, y, z, vp->res) ] < -0.5f)
+				if (vp->data_b[ V_I(x, y, z, vp->res) ] < -0.f)
 					vp->data_b[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_b, vp->res, x, y, z);
 			}
 		}
 	}
 }
 
+static void lightcache_filter2(VolumePrecache *vp)
+{
+	int x, y, z;
+	float *new_r, *new_g, *new_b;
+	int field_size = vp->res[0]*vp->res[1]*vp->res[2]*sizeof(float);
+	
+	new_r = MEM_mallocN(field_size, "temp buffer for light cache filter r channel");
+	new_g = MEM_mallocN(field_size, "temp buffer for light cache filter g channel");
+	new_b = MEM_mallocN(field_size, "temp buffer for light cache filter b channel");
+	
+	memcpy(new_r, vp->data_r, field_size);
+	memcpy(new_g, vp->data_g, field_size);
+	memcpy(new_b, vp->data_b, field_size);
+	
+	for (z=0; z < vp->res[2]; z++) {
+		for (y=0; y < vp->res[1]; y++) {
+			for (x=0; x < vp->res[0]; x++) {
+				/* trigger for outside mesh */
+				if (vp->data_r[ V_I(x, y, z, vp->res) ] < -0.f)
+					new_r[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_r, vp->res, x, y, z);
+				if (vp->data_g[ V_I(x, y, z, vp->res) ] < -0.f)
+					new_g[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_g, vp->res, x, y, z);
+				if (vp->data_b[ V_I(x, y, z, vp->res) ] < -0.f)
+					new_b[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_b, vp->res, x, y, z);
+			}
+		}
+	}
+	
+	SWAP(float *, vp->data_r, new_r);
+	SWAP(float *, vp->data_g, new_g);
+	SWAP(float *, vp->data_b, new_b);
+	
+	if (new_r) { MEM_freeN(new_r); new_r=NULL; }
+	if (new_g) { MEM_freeN(new_g); new_g=NULL; }
+	if (new_b) { MEM_freeN(new_b); new_b=NULL; }
+}
+
 static inline int ms_I(int x, int y, int z, int *n) //has a pad of 1 voxel surrounding the core for boundary simulation
 { 
 	return z*(n[1]+2)*(n[0]+2) + y*(n[0]+2) + x;
@@ -423,11 +460,10 @@ static void *vol_precache_part(void *data)
 	ObjectInstanceRen *obi = pa->obi;
 	RayTree *tree = pa->tree;
 	ShadeInput *shi = pa->shi;
-	float density, scatter_col[3] = {0.f, 0.f, 0.f};
+	float scatter_col[3] = {0.f, 0.f, 0.f};
 	float co[3];
 	int x, y, z;
 	const int res[3]= {pa->res[0], pa->res[1], pa->res[2]};
-	const float stepsize = vol_get_stepsize(shi, STEPSIZE_VIEW);
 
 	for (z= pa->minz; z < pa->maxz; z++) {
 		co[2] = pa->bbmin[2] + (pa->voxel[2] * (z + 0.5f));
@@ -448,8 +484,7 @@ static void *vol_precache_part(void *data)
 				
 				VecCopyf(shi->view, co);
 				Normalize(shi->view);
-				density = vol_get_density(shi, co);
-				vol_get_scattering(shi, scatter_col, co, stepsize, density);
+				vol_get_scattering(shi, scatter_col, co);
 			
 				obi->volume_precache->data_r[ V_I(x, y, z, res) ] = scatter_col[0];
 				obi->volume_precache->data_g[ V_I(x, y, z, res) ] = scatter_col[1];
diff --git a/source/blender/render/intern/source/volumetric.c b/source/blender/render/intern/source/volumetric.c
index b2692c25b99..381a32de027 100644
--- a/source/blender/render/intern/source/volumetric.c
+++ b/source/blender/render/intern/source/volumetric.c
@@ -170,29 +170,6 @@ static void vol_trace_behind(ShadeInput *shi, VlakRen *vlr, float *co, float *co
 	}
 }
 
-/* input shader data */
-
-float vol_get_stepsize(struct ShadeInput *shi, int context)
-{
-	if (shi->mat->vol.stepsize_type == MA_VOL_STEP_RANDOMIZED) {
-		/* range between 0.75 and 1.25 */
-		const float rnd = 0.5f * BLI_thread_frand(shi->thread) + 0.75f;
-	
-		if (context == STEPSIZE_VIEW)
-			return shi->mat->vol.stepsize * rnd;
-		else if (context == STEPSIZE_SHADE)
-			return shi->mat->vol.shade_stepsize * rnd;
-	}
-	else {	// MA_VOL_STEP_CONSTANT
-		
-		if (context == STEPSIZE_VIEW)
-			return shi->mat->vol.stepsize;
-		else if (context == STEPSIZE_SHADE)
-			return shi->mat->vol.shade_stepsize;
-	}
-	
-	return shi->mat->vol.stepsize;
-}
 
 /* trilinear interpolation */
 static void vol_get_precached_scattering(ShadeInput *shi, float *scatter_col, float *co)
@@ -212,9 +189,9 @@ static void vol_get_precached_scattering(ShadeInput *shi, float *scatter_col, fl
 	sample_co[1] = ((co[1] - bbmin[1]) / dim[1]);
 	sample_co[2] = ((co[2] - bbmin[2]) / dim[2]);
 
-	scatter_col[0] = voxel_sample_trilinear(vp->data_r, vp->res, sample_co);
-	scatter_col[1] = voxel_sample_trilinear(vp->data_g, vp->res, sample_co);
-	scatter_col[2] = voxel_sample_trilinear(vp->data_b, vp->res, sample_co);
+	scatter_col[0] = voxel_sample_triquadratic(vp->data_r, vp->res, sample_co);
+	scatter_col[1] = voxel_sample_triquadratic(vp->data_g, vp->res, sample_co);
+	scatter_col[2] = voxel_sample_triquadratic(vp->data_b, vp->res, sample_co);
 }
 
 /* Meta object density, brute force for now 
@@ -270,7 +247,8 @@ float vol_get_density(struct ShadeInput *shi, float *co)
 	float density = shi->mat->vol.density;
 	float density_scale = shi->mat->vol.density_scale;
 		
-	do_volume_tex(shi, co, MAP_DENSITY, NULL, &density);
+	if (shi->mat->mapto_textured & MAP_DENSITY)
+		do_volume_tex(shi, co, MAP_DENSITY, NULL, &density);
 	
 	// if meta-object, modulate by metadensity without increasing it
 	if (shi->obi->obr->ob->type == OB_MBALL) {
@@ -281,79 +259,110 @@ float vol_get_density(struct ShadeInput *shi, float *co)
 	return density * density_scale;
 }
 
-/* scattering multiplier, values above 1.0 are non-physical, 
- * but can be useful to tweak lighting */
-float vol_get_scattering_fac(ShadeInput *shi, float *co)
+
+/* Color of light that gets scattered out by the volume */
+/* Uses same physically based scattering parameter as in transmission calculations, 
+ * along with artificial reflection scale/reflection color tint */
+void vol_get_reflection_color(ShadeInput *shi, float *ref_col, float *co)
 {
 	float scatter = shi->mat->vol.scattering;
-	float col[3] = {0.0, 0.0, 0.0};
+	float reflection= shi->mat->vol.reflection;
+	VECCOPY(ref_col, shi->mat->vol.reflection_col);
+	
+	if (shi->mat->mapto_textured & (MAP_SCATTERING+MAP_REFLECTION_COL))
+		do_volume_tex(shi, co, MAP_SCATTERING+MAP_REFLECTION_COL, ref_col, &scatter);
 	
-	do_volume_tex(shi, co, MAP_SCATTERING, col, &scatter);
+	/* only one single float parameter at a time... :s */
+	if (shi->mat->mapto_textured & (MAP_REFLECTION))
+		do_volume_tex(shi, co, MAP_REFLECTION, NULL, &reflection);
 	
-	return scatter;
+	ref_col[0] = reflection * ref_col[0] * scatter;
+	ref_col[1] = reflection * ref_col[1] * scatter;
+	ref_col[2] = reflection * ref_col[2] * scatter;
 }
 
 /* compute emission component, amount of radiance to add per segment
  * can be textured with 'emit' */
-void vol_get_emission(ShadeInput *shi, float *emission_col, float *co, float density)
+void vol_get_emission(ShadeInput *shi, float *emission_col, float *co)
 {
 	float emission = shi->mat->vol.emission;
 	VECCOPY(emission_col, shi->mat->vol.emission_col);
 	
-	do_volume_tex(shi, co, MAP_EMISSION+MAP_EMISSION_COL, emission_col, &emission);
+	if (shi->mat->mapto_textured & (MAP_EMISSION+MAP_EMISSION_COL))
+		do_volume_tex(shi, co, MAP_EMISSION+MAP_EMISSION_COL, emission_col, &emission);
 	
-	emission_col[0] = emission_col[0] * emission * density;
-	emission_col[1] = emission_col[1] * emission * density;
-	emission_col[2] = emission_col[2] * emission * density;
+	emission_col[0] = emission_col[0] * emission;
+	emission_col[1] = emission_col[1] * emission;
+	emission_col[2] = emission_col[2] * emission;
 }
 
-void vol_get_absorption(ShadeInput *shi, float *absorb_col, float *co)
+
+/* A combination of scattering and absorption -> known as sigma T.
+ * This can possibly use a specific scattering colour, 
+ * and absorption multiplier factor too, but these parameters are left out for simplicity.
+ * It's easy enough to get a good wide range of results with just these two parameters. */
+void vol_get_sigma_t(ShadeInput *shi, float *sigma_t, float *co)
 {
-	float absorption = shi->mat->vol.absorption;
-	VECCOPY(absorb_col, shi->mat->vol.absorption_col);
+	/* technically absorption, but named transmission color 
+	 * since it describes the effect of the coloring *after* absorption */
+	float transmission_col[3] = {shi->mat->vol.transmission_col[0], shi->mat->vol.transmission_col[1], shi->mat->vol.transmission_col[2]};
+	float scattering = shi->mat->vol.scattering;
 	
-	do_volume_tex(shi, co, MAP_ABSORPTION+MAP_ABSORPTION_COL, absorb_col, &absorption);
+	if (shi->mat->mapto_textured & (MAP_SCATTERING+MAP_TRANSMISSION_COL))
+		do_volume_tex(shi, co, MAP_SCATTERING+MAP_TRANSMISSION_COL, transmission_col, &scattering);
 	
-	absorb_col[0] = (1.0f - absorb_col[0]) * absorption;
-	absorb_col[1] = (1.0f - absorb_col[1]) * absorption;
-	absorb_col[2] = (1.0f - absorb_col[2]) * absorption;
+	sigma_t[0] = (1.0f - transmission_col[0]) + scattering;
+	sigma_t[1] = (1.0f - transmission_col[1]) + scattering;
+	sigma_t[2] = (1.0f - transmission_col[2]) + scattering;
 }
 
 /* phase function - determines in which directions the light 
  * is scattered in the volume relative to incoming direction 
  * and view direction */
-float vol_get_phasefunc(ShadeInput *shi, short phasefunc_type, float g, float *w, float *wp)
+float vol_get_phasefunc(ShadeInput *shi, float g, float *w, float *wp)
 {
-	const float costheta = Inpf(w, wp);
-	const float scale = M_PI;
-	
-	/*
-	 * Scale constant is required, since Blender's shading system doesn't normalise for
-	 * energy conservation - eg. scaling by 1/pi for a lambert shader.
-	 * This makes volumes darker than other solid objects, for the same lighting intensity.
-	 * To correct this, scale up the phase function values
+	const float normalize = 0.25f; // = 1.f/4.f = M_PI/(4.f*M_PI)
+	
+	/* normalization constant is 1/4 rather than 1/4pi, since
+	 * Blender's shading system doesn't normalise for
+	 * energy conservation - eg. multiplying by pdf ( 1/pi for a lambert brdf ).
+	 * This means that lambert surfaces in Blender are pi times brighter than they 'should be'
+	 * and therefore, with correct energy conservation, volumes will darker than other solid objects,
+	 * for the same lighting intensity.
+	 * To correct this, scale up the phase function values by pi
 	 * until Blender's shading system supports this better. --matt
 	 */
 	
+	if (g == 0.f) {	/* isotropic */
+		return normalize * 1.f;
+	} else {		/* schlick */
+		const float k = 1.55f * g - .55f * g * g * g;
+		const float kcostheta = k * Inpf(w, wp);
+		return normalize * (1.f - k*k) / ((1.f - kcostheta) * (1.f - kcostheta));
+	}
+	
+	/*
+	 * not used, but here for reference:
 	switch (phasefunc_type) {
 		case MA_VOL_PH_MIEHAZY:
-			return scale * (0.5f + 4.5f * powf(0.5 * (1.f + costheta), 8.f)) / (4.f*M_PI);
+			return normalize * (0.5f + 4.5f * powf(0.5 * (1.f + costheta), 8.f));
 		case MA_VOL_PH_MIEMURKY:
-			return scale * (0.5f + 16.5f * powf(0.5 * (1.f + costheta), 32.f)) / (4.f*M_PI);
+			return normalize * (0.5f + 16.5f * powf(0.5 * (1.f + costheta), 32.f));
 		case MA_VOL_PH_RAYLEIGH:
-			return scale * 3.f/(16.f*M_PI) * (1 + costheta * costheta);
+			return normalize * 3.f/4.f * (1 + costheta * costheta);
 		case MA_VOL_PH_HG:
-			return scale * (1.f / (4.f * M_PI) * (1.f - g*g) / powf(1.f + g*g - 2.f * g * costheta, 1.5f));
+			return normalize * (1.f - g*g) / powf(1.f + g*g - 2.f * g * costheta, 1.5f));
 		case MA_VOL_PH_SCHLICK:
 		{
 			const float k = 1.55f * g - .55f * g * g * g;
 			const float kcostheta = k * costheta;
-			return scale * (1.f / (4.f * M_PI) * (1.f - k*k) / ((1.f - kcostheta) * (1.f - kcostheta)));
+			return normalize * (1.f - k*k) / ((1.f - kcostheta) * (1.f - kcostheta));
 		}
 		case MA_VOL_PH_ISOTROPIC:
 		default:
-			return scale * (1.f / (4.f * M_PI));
+			return normalize * 1.f;
 	}
+	*/
 }
 
 /* Compute transmittance = e^(-attenuation) */
@@ -361,15 +370,15 @@ void vol_get_transmittance_seg(ShadeInput *shi, float *tr, float stepsize, float
 {
 	/* input density = density at co */
 	float tau[3] = {0.f, 0.f, 0.f};
-	float absorb[3];
-	const float scatter_dens = vol_get_scattering_fac(shi, co) * density * stepsize;
-
-	vol_get_absorption(shi, absorb, co);
+	const float stepd = density * stepsize;
+	float sigma_t[3];
+	
+	vol_get_sigma_t(shi, sigma_t, co);
 	
 	/* homogenous volume within the sampled distance */
-	tau[0] += scatter_dens * absorb[0];
-	tau[1] += scatter_dens * absorb[1];
-	tau[2] += scatter_dens * absorb[2];
+	tau[0] += stepd * sigma_t[0];
+	tau[1] += stepd * sigma_t[1];
+	tau[2] += stepd * sigma_t[2];
 	
 	tr[0] *= exp(-tau[0]);
 	tr[1] *= exp(-tau[1]);
@@ -381,31 +390,29 @@ static void vol_get_transmittance(ShadeInput *shi, float *tr, float *co, float *
 {
 	float p[3] = {co[0], co[1], co[2]};
 	float step_vec[3] = {endco[0] - co[0], endco[1] - co[1], endco[2] - co[2]};
-	//const float ambtau = -logf(shi->mat->vol.depth_cutoff);	// never zero
 	float tau[3] = {0.f, 0.f, 0.f};
 
 	float t0 = 0.f;
 	float t1 = Normalize(step_vec);
 	float pt0 = t0;
 	
-	t0 += shi->mat->vol.shade_stepsize * ((shi->mat->vol.stepsize_type == MA_VOL_STEP_CONSTANT) ? 0.5f : BLI_thread_frand(shi->thread));
+	t0 += shi->mat->vol.stepsize * ((shi->mat->vol.stepsize_type == MA_VOL_STEP_CONSTANT) ? 0.5f : BLI_thread_frand(shi->thread));
 	p[0] += t0 * step_vec[0];
 	p[1] += t0 * step_vec[1];
 	p[2] += t0 * step_vec[2];
-	VecMulf(step_vec, shi->mat->vol.shade_stepsize);
+	VecMulf(step_vec, shi->mat->vol.stepsize);
 
-	for (; t0 < t1; pt0 = t0, t0 += shi->mat->vol.shade_stepsize) {
-		float absorb[3];
+	for (; t0 < t1; pt0 = t0, t0 += shi->mat->vol.stepsize) {
 		const float d = vol_get_density(shi, p);
 		const float stepd = (t0 - pt0) * d;
-		const float scatter_dens = vol_get_scattering_fac(shi, p) * stepd;
-		vol_get_absorption(shi, absorb, p);
+		float sigma_t[3];
+		
+		vol_get_sigma_t(shi, sigma_t, co);
 		
-		tau[0] += scatter_dens * absorb[0];
-		tau[1] += scatter_dens * absorb[1];
-		tau[2] += scatter_dens * absorb[2];
+		tau[0] += stepd * sigma_t[0];
+		tau[1] += stepd * sigma_t[1];
+		tau[2] += stepd * sigma_t[2];
 		
-		//if (luminance(tau) >= ambtau) break;
 		VecAddf(p, p, step_vec);
 	}
 	
@@ -420,8 +427,7 @@ void vol_shade_one_lamp(struct ShadeInput *shi, float *co, LampRen *lar, float *
 	float visifac, lv[3], lampdist;
 	float tr[3]={1.0,1.0,1.0};
 	float hitco[3], *atten_co;
-	float p;
-	float scatter_fac;
+	float p, ref_col[3];
 	
 	if (lar->mode & LA_LAYER) if((lar->lay & shi->obi->lay)==0) return;
 	if ((lar->lay & shi->lay)==0) return;
@@ -475,15 +481,20 @@ void vol_shade_one_lamp(struct ShadeInput *shi, float *co, LampRen *lar, float *
 		}
 	}
 	
-	p = vol_get_phasefunc(shi, shi->mat->vol.phasefunc_type, shi->mat->vol.phasefunc_g, shi->view, lv);
-	VecMulf(lacol, p);
+	if (luminance(lacol) < 0.001f) return;
+	
+	p = vol_get_phasefunc(shi, shi->mat->vol.asymmetry, shi->view, lv);
+	
+	/* physically based scattering with non-physically based RGB gain */
+	vol_get_reflection_color(shi, ref_col, co);
 	
-	scatter_fac = vol_get_scattering_fac(shi, co);
-	VecMulf(lacol, scatter_fac);
+	lacol[0] *= p * ref_col[0];
+	lacol[1] *= p * ref_col[1];
+	lacol[2] *= p * ref_col[2];
 }
 
 /* single scattering only for now */
-void vol_get_scattering(ShadeInput *shi, float *scatter_col, float *co, float stepsize, float density)
+void vol_get_scattering(ShadeInput *shi, float *scatter_col, float *co)
 {
 	ListBase *lights;
 	GroupObject *go;
@@ -515,71 +526,69 @@ outgoing radiance from behind surface * beam transmittance/attenuation
 	--> radiance for each segment = 
 		(radiance added by scattering + radiance added by emission) * beam transmittance/attenuation
 */
+
+/* For ease of use, I've also introduced a 'reflection' and 'reflection color' parameter, which isn't 
+ * physically correct. This works as an RGB tint/gain on out-scattered light, but doesn't affect the light 
+ * that is transmitted through the volume. While having wavelength dependent absorption/scattering is more correct,
+ * it also makes it harder to control the overall look of the volume since colouring the outscattered light results
+ * in the inverse colour being transmitted through the rest of the volume.
+ */
 static void volumeintegrate(struct ShadeInput *shi, float *col, float *co, float *endco)
 {
-	float tr[3] = {1.0f, 1.0f, 1.0f};
-	float radiance[3] = {0.f, 0.f, 0.f}, d_radiance[3] = {0.f, 0.f, 0.f};
-	float stepsize = vol_get_stepsize(shi, STEPSIZE_VIEW);
-	int nsteps, s;
-	float emit_col[3], scatter_col[3] = {0.0, 0.0, 0.0};
-	float stepvec[3], step_sta[3], step_end[3], step_mid[3];
-	float density;
-	const float depth_cutoff = shi->mat->vol.depth_cutoff;
-
-	/* ray marching */
-	nsteps = (int)((VecLenf(co, endco) / stepsize) + 0.5);
-	
-	VecSubf(stepvec, endco, co);
-	VecMulf(stepvec, 1.0f / nsteps);
-	VecCopyf(step_sta, co);
-	VecAddf(step_end, step_sta, stepvec);
-	
-	/* get radiance from all points along the ray due to participating media */
-	for (s = 0; s < nsteps; s++) {
-
-		density = vol_get_density(shi, step_sta);
+	float radiance[3] = {0.f, 0.f, 0.f};
+	float tr[3] = {1.f, 1.f, 1.f};
+	float p[3] = {co[0], co[1], co[2]};
+	float step_vec[3] = {endco[0] - co[0], endco[1] - co[1], endco[2] - co[2]};
+	const float stepsize = shi->mat->vol.stepsize;
+	
+	float t0 = 0.f;
+	float pt0 = t0;
+	float t1 = Normalize(step_vec);	/* returns vector length */
+	
+	t0 += stepsize * ((shi->mat->vol.stepsize_type == MA_VOL_STEP_CONSTANT) ? 0.5f : BLI_thread_frand(shi->thread));
+	p[0] += t0 * step_vec[0];
+	p[1] += t0 * step_vec[1];
+	p[2] += t0 * step_vec[2];
+	VecMulf(step_vec, stepsize);
+	
+	for (; t0 < t1; pt0 = t0, t0 += stepsize) {
+		const float density = vol_get_density(shi, p);
 		
-		/* there's only any use in shading here if there's actually some density to shade! */
 		if (density > 0.01f) {
-		
+			float scatter_col[3], emit_col[3];
+			const float stepd = (t0 - pt0) * density;
+			
 			/* transmittance component (alpha) */
 			vol_get_transmittance_seg(shi, tr, stepsize, co, density);
-
-			step_mid[0] = step_sta[0] + (stepvec[0] * 0.5);
-			step_mid[1] = step_sta[1] + (stepvec[1] * 0.5);
-			step_mid[2] = step_sta[2] + (stepvec[2] * 0.5);
-		
-			/* incoming light via emission or scattering (additive) */
-			vol_get_emission(shi, emit_col, step_mid, density);
 			
-			if (shi->obi->volume_precache)
-				vol_get_precached_scattering(shi, scatter_col, step_mid);
-			else
-				vol_get_scattering(shi, scatter_col, step_mid, stepsize, density);
+			if (luminance(tr) < shi->mat->vol.depth_cutoff) break;
 			
-			VecMulf(scatter_col, density);
-			VecAddf(d_radiance, emit_col, scatter_col);
+			vol_get_emission(shi, emit_col, p);
 			
-			/*   Lv += Tr * (Lve() + Ld) */
-			VecMulVecf(d_radiance, tr, d_radiance);
-			VecMulf(d_radiance, stepsize);
+			if (shi->obi->volume_precache) {
+				float p2[3];
+				
+				p2[0] = p[0] + (step_vec[0] * 0.5);
+				p2[1] = p[1] + (step_vec[1] * 0.5);
+				p2[2] = p[2] + (step_vec[2] * 0.5);
+				
+				vol_get_precached_scattering(shi, scatter_col, p2);
+			} else
+				vol_get_scattering(shi, scatter_col, p);
 			
-			VecAddf(radiance, radiance, d_radiance);	
+			radiance[0] += stepd * tr[0] * (emit_col[0] + scatter_col[0]);
+			radiance[1] += stepd * tr[1] * (emit_col[1] + scatter_col[1]);
+			radiance[2] += stepd * tr[2] * (emit_col[2] + scatter_col[2]);
 		}
-
-		VecCopyf(step_sta, step_end);
-		VecAddf(step_end, step_end, stepvec);
-		
-		/* luminance rec. 709 */
-		if ((0.2126*tr[0] + 0.7152*tr[1] + 0.0722*tr[2]) < depth_cutoff) break;	
+		VecAddf(p, p, step_vec);
 	}
 	
-	/* multiply original color (behind volume) with beam transmittance over entire distance */
-	VecMulVecf(col, tr, col);	
+	/* multiply original color (from behind volume) with transmittance over entire distance */
+	VecMulVecf(col, tr, col);
 	VecAddf(col, col, radiance);
 	
 	/* alpha <-- transmission luminance */
-	col[3] = 1.0f -(0.2126*tr[0] + 0.7152*tr[1] + 0.0722*tr[2]);
+	col[3] = 1.0f - luminance(tr);
 }
 
 /* the main entry point for volume shading */
@@ -606,7 +615,7 @@ static void volume_trace(struct ShadeInput *shi, struct ShadeResult *shr, int in
 		/* don't render the backfaces of ztransp volume materials.
 		 
 		 * volume shading renders the internal volume from between the
-		 * near view intersection of the solid volume to the
+		 * ' view intersection of the solid volume to the
 		 * intersection on the other side, as part of the shading of
 		 * the front face.
 		 
@@ -709,8 +718,7 @@ void shade_volume_shadow(struct ShadeInput *shi, struct ShadeResult *shr, struct
 	vol_get_transmittance(shi, tr, startco, endco);
 	
 	VecCopyf(shr->combined, tr);
-	shr->combined[3] = 1.0f -(0.2126*tr[0] + 0.7152*tr[1] + 0.0722*tr[2]);
-	shr->alpha = shr->combined[3];
+	shr->combined[3] = 1.0f - luminance(tr);
 }
 
 
diff --git a/source/blender/render/intern/source/voxeldata.c b/source/blender/render/intern/source/voxeldata.c
index 9318d37620c..479f33c9ff2 100644
--- a/source/blender/render/intern/source/voxeldata.c
+++ b/source/blender/render/intern/source/voxeldata.c
@@ -226,9 +226,6 @@ void make_voxeldata(struct Render *re)
 {
     Tex *tex;
 	
-	if(re->scene->r.scemode & R_PREVIEWBUTS)
-		return;
-	
 	re->i.infostr= "Loading voxel datasets";
 	re->stats_draw(re->sdh, &re->i);
 	
@@ -259,9 +256,6 @@ void free_voxeldata(Render *re)
 {
 	Tex *tex;
 	
-	if(re->scene->r.scemode & R_PREVIEWBUTS)
-		return;
-	
 	for (tex= G.main->tex.first; tex; tex= tex->id.next) {
 		if(tex->id.us && tex->type==TEX_VOXELDATA) {
 			free_voxeldata_one(re, tex);