svn merge -r 22571:22800 https://svn.blender.org/svnroot/bf-blender/trunk/blender

svn merge -r 22800:23207 https://svn.blender.org/svnroot/bf-blender/trunk/blender

Merged volumetric with new raytrace code (it compiles and rendered volume-cube.blend withouth problems)

Part1:
	source/blender

28 files changed, 3089 insertions, 186 deletions

diff --git a/source/blender/render/intern/include/pointdensity.h b/source/blender/render/intern/include/pointdensity.h
new file mode 100644
index 00000000000..93cdef3b14e
--- /dev/null
+++ b/source/blender/render/intern/include/pointdensity.h
@@ -0,0 +1,44 @@
+/*
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Matt Ebb
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#ifndef POINTDENSITY_H
+#define POINTDENSITY_H 
+
+/**
+ * Make point density kd-trees for all point density textures in the scene
+ */
+
+struct Render;
+struct TexResult;
+
+void make_pointdensities(struct Render *re);
+void free_pointdensities(struct Render *re);
+int pointdensitytex(struct Tex *tex, float *texvec, struct TexResult *texres);
+
+#endif /* POINTDENSITY_H */
+
diff --git a/source/blender/render/intern/include/rayobject.h b/source/blender/render/intern/include/rayobject.h
index 2a9b42d22fe..19608fba262 100644
--- a/source/blender/render/intern/include/rayobject.h
+++ b/source/blender/render/intern/include/rayobject.h
@@ -152,13 +152,6 @@ typedef struct RayObjectAPI
 } RayObjectAPI;
 
 
-
-
-/*
- * Extend min/max coords so that the rayobject is inside them
- */
-void RE_rayobject_merge_bb(RayObject *ob, float *min, float *max);
-
 /*
  * This function differs from RE_rayobject_raycast
  * RE_rayobject_intersect does NOT perform last-hit optimization
diff --git a/source/blender/render/intern/include/render_types.h b/source/blender/render/intern/include/render_types.h
index c412921eb31..67c69259213 100644
--- a/source/blender/render/intern/include/render_types.h
+++ b/source/blender/render/intern/include/render_types.h
@@ -205,6 +205,9 @@ struct Render
 	ListBase customdata_names;
 
 	struct Object *excludeob;
+	ListBase render_volumes_inside;
+	ListBase volumes;
+	ListBase volume_precache_parts;
 
 	/* arena for allocating data for use during render, for
 		* example dynamic TFaces to go in the VlakRen structure.
@@ -302,7 +305,9 @@ typedef struct ObjectInstanceRen {
 
 	float dupliorco[3], dupliuv[2];
 	float (*duplitexmat)[4];
-
+	
+	struct VolumePrecache *volume_precache;
+	
 	float *vectors;
 	int totvector;
 	
@@ -416,6 +421,46 @@ typedef struct StrandRen {
 	float orco[3];
 } StrandRen;
 
+/* ------------------------------------------------------------------------- */
+
+typedef struct VolumeOb
+{
+	struct VolumeOb *next, *prev;
+	struct Material *ma;
+	struct ObjectRen *obr;
+} VolumeOb;
+
+typedef struct MatInside {
+	struct MatInside *next, *prev;
+	struct Material *ma;
+	struct ObjectInstanceRen *obi;
+} MatInside;
+
+typedef struct VolPrecachePart
+{
+	struct VolPrecachePart *next, *prev;
+	struct RayObject *tree;
+	struct ShadeInput *shi;
+	struct ObjectInstanceRen *obi;
+	int num;
+	int minx, maxx;
+	int miny, maxy;
+	int minz, maxz;
+	int res[3];
+	float bbmin[3];
+	float voxel[3];
+	int working, done;
+} VolPrecachePart;
+
+typedef struct VolumePrecache
+{
+	int res[3];
+	float *data_r;
+	float *data_g;
+	float *data_b;
+} VolumePrecache;
+
+/* ------------------------------------------------------------------------- */
 
 struct LampRen;
 struct MTex;
diff --git a/source/blender/render/intern/include/rendercore.h b/source/blender/render/intern/include/rendercore.h
index 4b28529a147..250fbc000cb 100644
--- a/source/blender/render/intern/include/rendercore.h
+++ b/source/blender/render/intern/include/rendercore.h
@@ -95,6 +95,7 @@ int get_sample_layers(struct RenderPart *pa, struct RenderLayer *rl, struct Rend
 
 extern void freeraytree(Render *re);
 extern void makeraytree(Render *re);
+RayObject* makeraytree_object(Render *re, ObjectInstanceRen *obi);
 
 extern void ray_shadow(ShadeInput *, LampRen *, float *);
 extern void ray_trace(ShadeInput *, ShadeResult *);
diff --git a/source/blender/render/intern/include/shading.h b/source/blender/render/intern/include/shading.h
index d195f32d5ef..95bccd2be1e 100644
--- a/source/blender/render/intern/include/shading.h
+++ b/source/blender/render/intern/include/shading.h
@@ -33,6 +33,7 @@ struct VlakRen;
 struct StrandSegment;
 struct StrandPoint;
 struct ObjectInstanceRen obi;
+struct Isect;
 
 /* shadeinput.c */
 
@@ -52,6 +53,7 @@ typedef struct ShadeSample {
 
 	/* also the node shader callback */
 void shade_material_loop(struct ShadeInput *shi, struct ShadeResult *shr);
+void shade_volume_loop(struct ShadeInput *shi, struct ShadeResult *shr);
 
 void shade_input_set_triangle_i(struct ShadeInput *shi, struct ObjectInstanceRen *obi, struct VlakRen *vlr, short i1, short i2, short i3);
 void shade_input_set_triangle(struct ShadeInput *shi, volatile int obi, volatile int facenr, int normal_flip);
@@ -87,7 +89,11 @@ void shade_color(struct ShadeInput *shi, ShadeResult *shr);
 void ambient_occlusion_to_diffuse(struct ShadeInput *shi, float *diff);
 void ambient_occlusion(struct ShadeInput *shi);
 
+ListBase *get_lights(struct ShadeInput *shi);
 float lamp_get_visibility(struct LampRen *lar, float *co, float *lv, float *dist);
 void lamp_get_shadow(struct LampRen *lar, ShadeInput *shi, float inp, float *shadfac, int do_real);
 
 float	fresnel_fac(float *view, float *vn, float fresnel, float fac);
+
+/* rayshade.c */
+extern void shade_ray(struct Isect *is, struct ShadeInput *shi, struct ShadeResult *shr);
diff --git a/source/blender/render/intern/include/texture.h b/source/blender/render/intern/include/texture.h
index c254b768292..78d6a912af1 100644
--- a/source/blender/render/intern/include/texture.h
+++ b/source/blender/render/intern/include/texture.h
@@ -56,6 +56,7 @@ void do_halo_tex(struct HaloRen *har, float xn, float yn, float *colf);
 void do_sky_tex(float *rco, float *lo, float *dxyview, float *hor, float *zen, float *blend, int skyflag, short thread);
 void do_material_tex(struct ShadeInput *shi);
 void do_lamp_tex(LampRen *la, float *lavec, struct ShadeInput *shi, float *colf, int effect);
+void do_volume_tex(struct ShadeInput *shi, float *xyz, int mapto_flag, float *col, float *val);
 
 void init_render_textures(Render *re);
 void end_render_textures(void);
diff --git a/source/blender/render/intern/include/volume_precache.h b/source/blender/render/intern/include/volume_precache.h
new file mode 100644
index 00000000000..9d87a219c82
--- /dev/null
+++ b/source/blender/render/intern/include/volume_precache.h
@@ -0,0 +1,33 @@
+/**
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Matt Ebb.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+ 
+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
diff --git a/source/blender/render/intern/include/volumetric.h b/source/blender/render/intern/include/volumetric.h
new file mode 100644
index 00000000000..026b4840ea3
--- /dev/null
+++ b/source/blender/render/intern/include/volumetric.h
@@ -0,0 +1,47 @@
+/**
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Matt Ebb.
+ *
+ * ***** 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 shade_volume_outside(ShadeInput *shi, ShadeResult *shr);
+void shade_volume_inside(ShadeInput *shi, ShadeResult *shr);
+void shade_volume_shadow(struct ShadeInput *shi, struct ShadeResult *shr, struct Isect *last_is);
+
+#define STEPSIZE_VIEW	0
+#define STEPSIZE_SHADE	1
+
+#define VOL_IS_BACKFACE			1
+#define VOL_IS_SAMEMATERIAL		2
+
+#define VOL_BOUNDS_DEPTH	0
+#define VOL_BOUNDS_SS		1
+
+#define VOL_SHADE_OUTSIDE	0
+#define VOL_SHADE_INSIDE	1
diff --git a/source/blender/render/intern/include/voxeldata.h b/source/blender/render/intern/include/voxeldata.h
new file mode 100644
index 00000000000..b291bdc096d
--- /dev/null
+++ b/source/blender/render/intern/include/voxeldata.h
@@ -0,0 +1,45 @@
+/**
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Raul Fernandez Hernandez (Farsthary), Matt Ebb.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#ifndef VOXELDATA_H
+#define VOXELDATA_H 
+
+struct Render;
+struct TexResult;
+
+typedef struct VoxelDataHeader
+{
+	int resolX, resolY, resolZ;
+	int frames;
+} VoxelDataHeader;
+
+void make_voxeldata(struct Render *re);
+void free_voxeldata(struct Render *re);
+int voxeldatatex(struct Tex *tex, float *texvec, struct TexResult *texres);
+
+#endif /* VOXELDATA_H */
diff --git a/source/blender/render/intern/source/Makefile b/source/blender/render/intern/source/Makefile
index 3c8d0f637a3..c313549f9b9 100644
--- a/source/blender/render/intern/source/Makefile
+++ b/source/blender/render/intern/source/Makefile
@@ -39,7 +39,6 @@ CFLAGS += $(LEVEL_1_C_WARNINGS)
 # The external modules follow after. There should be a nicer way to say this.
 CPPFLAGS += -I../include 
 CPPFLAGS += -I../../extern/include 
-CPPFLAGS += -I../../../radiosity/extern/include 
 CPPFLAGS += -I../../../blenlib
 CPPFLAGS += -I../../../imbuf
 CPPFLAGS += -I../../../makesdna
@@ -51,6 +50,7 @@ CPPFLAGS += -I$(NAN_GUARDEDALLOC)/include
 # not very neat: the rest of blender..
 CPPFLAGS += -I../../../editors/include 
 CPPFLAGS += $(NAN_SDLCFLAGS)
+CPPFLAGS += -I../../../../../intern/smoke/extern
 
 ifeq ($(WITH_QUICKTIME), true)
 	CPPFLAGS += -DWITH_QUICKTIME
diff --git a/source/blender/render/intern/source/convertblender.c b/source/blender/render/intern/source/convertblender.c
index 2c264ce2337..99825c0c2ff 100644
--- a/source/blender/render/intern/source/convertblender.c
+++ b/source/blender/render/intern/source/convertblender.c
@@ -31,7 +31,7 @@
 #include <string.h>
 #include <limits.h>
 
-#include "MTC_matrixops.h"
+
 
 #include "MEM_guardedalloc.h"
 
@@ -100,6 +100,8 @@
 
 #include "envmap.h"
 #include "occlusion.h"
+#include "pointdensity.h"
+#include "voxeldata.h"
 #include "render_types.h"
 #include "rendercore.h"
 #include "renderdatabase.h"
@@ -108,6 +110,7 @@
 #include "shading.h"
 #include "strand.h"
 #include "texture.h"
+#include "volume_precache.h"
 #include "sss.h"
 #include "strand.h"
 #include "zbuf.h"
@@ -188,8 +191,8 @@ void RE_make_stars(Render *re, Scene *scenev3d, void (*initfunc)(void),
 	if (re) re->flag |= R_HALO;
 	else stargrid *= 1.0;				/* then it draws fewer */
 	
-	if(re) MTC_Mat4Invert(mat, re->viewmat);
-	else MTC_Mat4One(mat);
+	if(re) Mat4Invert(mat, re->viewmat);
+	else Mat4One(mat);
 	
 	/* BOUNDING BOX CALCULATION
 		* bbox goes from z = loc_near_var | loc_far_var,
@@ -237,7 +240,7 @@ void RE_make_stars(Render *re, Scene *scenev3d, void (*initfunc)(void),
 					done++;
 				}
 				else {
-					MTC_Mat4MulVecfl(re->viewmat, vec);
+					Mat4MulVecfl(re->viewmat, vec);
 					
 					/* in vec are global coordinates
 					* calculate distance to camera
@@ -826,7 +829,7 @@ static void autosmooth(Render *re, ObjectRen *obr, float mat[][4], int degr)
 	/* rotate vertices and calculate normal of faces */
 	for(a=0; a<obr->totvert; a++) {
 		ver= RE_findOrAddVert(obr, a);
-		MTC_Mat4MulVecfl(mat, ver->co);
+		Mat4MulVecfl(mat, ver->co);
 	}
 	for(a=0; a<obr->totvlak; a++) {
 		vlr= RE_findOrAddVlak(obr, a);
@@ -917,6 +920,7 @@ static Material *give_render_material(Render *re, Object *ob, int nr)
 	
 	if(re->r.mode & R_SPEED) ma->texco |= NEED_UV;
 	
+	if(ma->material_type == MA_TYPE_VOLUME) ma->mode |= MA_TRANSP;
 	if((ma->mode & MA_TRANSP) && (ma->mode & MA_ZTRANSP))
 		re->flag |= R_ZTRA;
 	
@@ -1276,19 +1280,19 @@ static void particle_billboard(Render *re, ObjectRen *obr, Material *ma, Particl
 
 	VECADD(vlr->v1->co, bb_center, xvec);
 	VECADD(vlr->v1->co, vlr->v1->co, yvec);
-	MTC_Mat4MulVecfl(re->viewmat, vlr->v1->co);
+	Mat4MulVecfl(re->viewmat, vlr->v1->co);
 
 	VECSUB(vlr->v2->co, bb_center, xvec);
 	VECADD(vlr->v2->co, vlr->v2->co, yvec);
-	MTC_Mat4MulVecfl(re->viewmat, vlr->v2->co);
+	Mat4MulVecfl(re->viewmat, vlr->v2->co);
 
 	VECSUB(vlr->v3->co, bb_center, xvec);
 	VECSUB(vlr->v3->co, vlr->v3->co, yvec);
-	MTC_Mat4MulVecfl(re->viewmat, vlr->v3->co);
+	Mat4MulVecfl(re->viewmat, vlr->v3->co);
 
 	VECADD(vlr->v4->co, bb_center, xvec);
 	VECSUB(vlr->v4->co, vlr->v4->co, yvec);
-	MTC_Mat4MulVecfl(re->viewmat, vlr->v4->co);
+	Mat4MulVecfl(re->viewmat, vlr->v4->co);
 
 	CalcNormFloat4(vlr->v4->co, vlr->v3->co, vlr->v2->co, vlr->v1->co, vlr->n);
 	VECCOPY(vlr->v1->n,vlr->n);
@@ -1388,7 +1392,7 @@ static void particle_normal_ren(short ren_as, ParticleSettings *part, Render *re
 	VECCOPY(loc, state->co);
 
 	if(ren_as != PART_DRAW_BB)
-		MTC_Mat4MulVecfl(re->viewmat, loc);
+		Mat4MulVecfl(re->viewmat, loc);
 
 	switch(ren_as) {
 		case PART_DRAW_LINE:
@@ -1397,7 +1401,7 @@ static void particle_normal_ren(short ren_as, ParticleSettings *part, Render *re
 			sd->size = hasize;
 
 			VECCOPY(vel, state->vel);
-			MTC_Mat4Mul3Vecfl(re->viewmat, vel);
+			Mat4Mul3Vecfl(re->viewmat, vel);
 			Normalize(vel);
 
 			if(part->draw & PART_DRAW_VEL_LENGTH)
@@ -1617,8 +1621,8 @@ static int render_new_particle_system(Render *re, ObjectRen *obr, ParticleSystem
 	}
 	
 /* 2.5 setup matrices */
-	MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-	MTC_Mat4Invert(ob->imat, mat);	/* need to be that way, for imat texture */
+	Mat4MulMat4(mat, ob->obmat, re->viewmat);
+	Mat4Invert(ob->imat, mat);	/* need to be that way, for imat texture */
 	Mat3CpyMat4(nmat, ob->imat);
 	Mat3Transp(nmat);
 
@@ -1739,8 +1743,10 @@ static int render_new_particle_system(Render *re, ObjectRen *obr, ParticleSystem
 
 			pa_size = pa->size;
 
-			r_tilt = 1.0f + pa->r_ave[0];
-			r_length = 0.5f * (1.0f + pa->r_ave[1]);
+			BLI_srandom(psys->seed+a);
+
+			r_tilt = 2.0f*(BLI_frand() - 0.5f);
+			r_length = BLI_frand();
 
 			if(path_nbr) {
 				cache = psys->pathcache[a];
@@ -1898,7 +1904,7 @@ static int render_new_particle_system(Render *re, ObjectRen *obr, ParticleSystem
 				time= curlen/strandlen;
 
 				VECCOPY(loc,state.co);
-				MTC_Mat4MulVecfl(re->viewmat,loc);
+				Mat4MulVecfl(re->viewmat,loc);
 
 				if(strandbuf) {
 					VECCOPY(svert->co, loc);
@@ -2040,8 +2046,8 @@ static void make_render_halos(Render *re, ObjectRen *obr, Mesh *me, int totvert,
 	float vec[3], hasize, mat[4][4], imat[3][3];
 	int a, ok, seed= ma->seed1;
 
-	MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-	MTC_Mat3CpyMat4(imat, ob->imat);
+	Mat4MulMat4(mat, ob->obmat, re->viewmat);
+	Mat3CpyMat4(imat, ob->imat);
 
 	re->flag |= R_HALO;
 
@@ -2052,7 +2058,7 @@ static void make_render_halos(Render *re, ObjectRen *obr, Mesh *me, int totvert,
 			hasize= ma->hasize;
 
 			VECCOPY(vec, mvert->co);
-			MTC_Mat4MulVecfl(mat, vec);
+			Mat4MulVecfl(mat, vec);
 
 			if(ma->mode & MA_HALOPUNO) {
 				xn= mvert->no[0];
@@ -2185,7 +2191,7 @@ static void displace_render_vert(Render *re, ObjectRen *obr, ShadeInput *shi, Ve
 	}
 	if (texco & TEXCO_GLOB) {
 		VECCOPY(shi->gl, shi->co);
-		MTC_Mat4MulVecfl(re->viewinv, shi->gl);
+		Mat4MulVecfl(re->viewinv, shi->gl);
 	}
 	if (texco & TEXCO_NORM) {
 		VECCOPY(shi->orn, shi->vn);
@@ -2326,9 +2332,9 @@ static void init_render_mball(Render *re, ObjectRen *obr)
 	if (ob!=find_basis_mball(re->scene, ob))
 		return;
 
-	MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-	MTC_Mat4Invert(ob->imat, mat);
-	MTC_Mat3CpyMat4(imat, ob->imat);
+	Mat4MulMat4(mat, ob->obmat, re->viewmat);
+	Mat4Invert(ob->imat, mat);
+	Mat3CpyMat4(imat, ob->imat);
 
 	ma= give_render_material(re, ob, 1);
 
@@ -2349,7 +2355,7 @@ static void init_render_mball(Render *re, ObjectRen *obr)
 
 		ver= RE_findOrAddVert(obr, obr->totvert++);
 		VECCOPY(ver->co, data);
-		MTC_Mat4MulVecfl(mat, ver->co);
+		Mat4MulVecfl(mat, ver->co);
 
 		/* render normals are inverted */
 		xn= -nors[0];
@@ -2433,7 +2439,7 @@ static int dl_surf_to_renderdata(ObjectRen *obr, DispList *dl, Material **matar,
 		if(orco) {
 			v1->orco= orco; orco+= 3; orcoret++;
 		}	
-		MTC_Mat4MulVecfl(mat, v1->co);
+		Mat4MulVecfl(mat, v1->co);
 		
 		for (v = 1; v < sizev; v++) {
 			ver= RE_findOrAddVert(obr, obr->totvert++);
@@ -2441,7 +2447,7 @@ static int dl_surf_to_renderdata(ObjectRen *obr, DispList *dl, Material **matar,
 			if(orco) {
 				ver->orco= orco; orco+= 3; orcoret++;
 			}	
-			MTC_Mat4MulVecfl(mat, ver->co);
+			Mat4MulVecfl(mat, ver->co);
 		}
 		/* if V-cyclic, add extra vertices at end of the row */
 		if (dl->flag & DL_CYCL_U) {
@@ -2591,8 +2597,8 @@ static void init_render_surf(Render *re, ObjectRen *obr)
 	nu= cu->nurb.first;
 	if(nu==0) return;
 
-	MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-	MTC_Mat4Invert(ob->imat, mat);
+	Mat4MulMat4(mat, ob->obmat, re->viewmat);
+	Mat4Invert(ob->imat, mat);
 
 	/* material array */
 	totmat= ob->totcol+1;
@@ -2652,8 +2658,8 @@ static void init_render_curve(Render *re, ObjectRen *obr, int timeoffset)
 	dl= cu->disp.first;
 	if(cu->disp.first==NULL) return;
 	
-	MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-	MTC_Mat4Invert(ob->imat, mat);
+	Mat4MulMat4(mat, ob->obmat, re->viewmat);
+	Mat4Invert(ob->imat, mat);
 
 	/* material array */
 	totmat= ob->totcol+1;
@@ -2695,7 +2701,7 @@ static void init_render_curve(Render *re, ObjectRen *obr, int timeoffset)
 					ver->flag = 0;
 				}
 
-				MTC_Mat4MulVecfl(mat, ver->co);
+				Mat4MulVecfl(mat, ver->co);
 				
 				if (orco) {
 					ver->orco = orco;
@@ -2747,7 +2753,7 @@ static void init_render_curve(Render *re, ObjectRen *obr, int timeoffset)
 					ver= RE_findOrAddVert(obr, obr->totvert++);
 						
 					VECCOPY(ver->co, fp);
-					MTC_Mat4MulVecfl(mat, ver->co);
+					Mat4MulVecfl(mat, ver->co);
 					fp+= 3;
 
 					if (orco) {
@@ -2976,6 +2982,53 @@ static void use_mesh_edge_lookup(ObjectRen *obr, DerivedMesh *dm, MEdge *medge,
 	}
 }
 
+static void free_camera_inside_volumes(Render *re)
+{
+	BLI_freelistN(&re->render_volumes_inside);
+}
+
+static void init_camera_inside_volumes(Render *re)
+{
+	ObjectInstanceRen *obi;
+	VolumeOb *vo;
+	float co[3] = {0.f, 0.f, 0.f};
+
+	for(vo= re->volumes.first; vo; vo= vo->next) {
+		for(obi= re->instancetable.first; obi; obi= obi->next) {
+			if (obi->obr == vo->obr) {
+				if (point_inside_volume_objectinstance(obi, co)) {
+					MatInside *mi;
+					
+					mi = MEM_mallocN(sizeof(MatInside), "camera inside material");
+					mi->ma = vo->ma;
+					mi->obi = obi;
+					
+					BLI_addtail(&(re->render_volumes_inside), mi);
+				}
+			}
+		}
+	}
+	
+	/* debug {
+	MatInside *m;
+	for (m=re->render_volumes_inside.first; m; m=m->next) {
+		printf("matinside: ma: %s \n", m->ma->id.name+2);
+	}
+	}*/
+}
+
+static void add_volume(Render *re, ObjectRen *obr, Material *ma)
+{
+	struct VolumeOb *vo;
+	
+	vo = MEM_mallocN(sizeof(VolumeOb), "volume object");
+	
+	vo->ma = ma;
+	vo->obr = obr;
+	
+	BLI_addtail(&re->volumes, vo);
+}
+
 static void init_render_mesh(Render *re, ObjectRen *obr, int timeoffset)
 {
 	Object *ob= obr->ob;
@@ -2997,9 +3050,9 @@ static void init_render_mesh(Render *re, ObjectRen *obr, int timeoffset)
 
 	me= ob->data;
 
-	MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-	MTC_Mat4Invert(ob->imat, mat);
-	MTC_Mat3CpyMat4(imat, ob->imat);
+	Mat4MulMat4(mat, ob->obmat, re->viewmat);
+	Mat4Invert(ob->imat, mat);
+	Mat3CpyMat4(imat, ob->imat);
 
 	if(me->totvert==0)
 		return;
@@ -3080,7 +3133,7 @@ static void init_render_mesh(Render *re, ObjectRen *obr, int timeoffset)
 			ver= RE_findOrAddVert(obr, obr->totvert++);
 			VECCOPY(ver->co, mvert->co);
 			if(do_autosmooth==0)	/* autosmooth on original unrotated data to prevent differences between frames */
-				MTC_Mat4MulVecfl(mat, ver->co);
+				Mat4MulVecfl(mat, ver->co);
   
 			if(orco) {
 				ver->orco= orco;
@@ -3301,13 +3354,13 @@ static void initshadowbuf(Render *re, LampRen *lar, float mat[][4])
 	shb->soft= lar->soft;
 	shb->shadhalostep= lar->shadhalostep;
 	
-	MTC_Mat4Ortho(mat);
-	MTC_Mat4Invert(shb->winmat, mat);	/* winmat is temp */
+	Mat4Ortho(mat);
+	Mat4Invert(shb->winmat, mat);	/* winmat is temp */
 	
 	/* matrix: combination of inverse view and lampmat */
 	/* calculate again: the ortho-render has no correct viewinv */
-	MTC_Mat4Invert(viewinv, re->viewmat);
-	MTC_Mat4MulMat4(shb->viewmat, viewinv, shb->winmat);
+	Mat4Invert(viewinv, re->viewmat);
+	Mat4MulMat4(shb->viewmat, viewinv, shb->winmat);
 	
 	/* projection */
 	shb->d= lar->clipsta;
@@ -3385,11 +3438,11 @@ static GroupObject *add_render_lamp(Render *re, Object *ob)
 	BLI_addtail(&re->lampren, lar);
 	go->lampren= lar;
 
-	MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-	MTC_Mat4Invert(ob->imat, mat);
+	Mat4MulMat4(mat, ob->obmat, re->viewmat);
+	Mat4Invert(ob->imat, mat);
 
-	MTC_Mat3CpyMat4(lar->mat, mat);
-	MTC_Mat3CpyMat4(lar->imat, ob->imat);
+	Mat3CpyMat4(lar->mat, mat);
+	Mat3CpyMat4(lar->imat, ob->imat);
 
 	lar->bufsize = la->bufsize;
 	lar->samp = la->samp;
@@ -3545,7 +3598,7 @@ static GroupObject *add_render_lamp(Render *re, Object *ob)
 			lar->sh_invcampos[0]= -lar->co[0];
 			lar->sh_invcampos[1]= -lar->co[1];
 			lar->sh_invcampos[2]= -lar->co[2];
-			MTC_Mat3MulVecfl(lar->imat, lar->sh_invcampos);
+			Mat3MulVecfl(lar->imat, lar->sh_invcampos);
 
 			/* z factor, for a normalized volume */
 			angle= saacos(lar->spotsi);
@@ -3809,6 +3862,7 @@ static void set_fullsample_flag(Render *re, ObjectRen *obr)
 			vlr->flag |= R_FULL_OSA;
 		else if(trace) {
 			if(mode & MA_SHLESS);
+			else if(vlr->mat->material_type == MA_TYPE_VOLUME);
 			else if((mode & MA_RAYMIRROR) || ((mode & MA_TRANSP) && (mode & MA_RAYTRANSP)))
 				/* for blurry reflect/refract, better to take more samples 
 				 * inside the raytrace than as OSA samples */
@@ -4168,7 +4222,7 @@ static void set_dupli_tex_mat(Render *re, ObjectInstanceRen *obi, DupliObject *d
 
 		obi->duplitexmat= BLI_memarena_alloc(re->memArena, sizeof(float)*4*4);
 		Mat4Invert(imat, dob->mat);
-		MTC_Mat4MulSerie(obi->duplitexmat, re->viewmat, dob->omat, imat, re->viewinv, 0, 0, 0, 0);
+		Mat4MulSerie(obi->duplitexmat, re->viewmat, dob->omat, imat, re->viewinv, 0, 0, 0, 0);
 	}
 }
 
@@ -4216,7 +4270,7 @@ static void add_render_object(Render *re, Object *ob, Object *par, DupliObject *
 	ObjectRen *obr;
 	ObjectInstanceRen *obi;
 	ParticleSystem *psys;
-	int show_emitter, allow_render= 1, index, psysindex;
+	int show_emitter, allow_render= 1, index, psysindex, i;
 
 	index= (dob)? dob->index: 0;
 
@@ -4252,6 +4306,12 @@ static void add_render_object(Render *re, Object *ob, Object *par, DupliObject *
 		}
 		else
 			find_dupli_instances(re, obr);
+			
+		for (i=1; i<=ob->totcol; i++) {
+			Material* ma = give_render_material(re, ob, i);
+			if (ma && ma->material_type == MA_TYPE_VOLUME)
+				add_volume(re, obr, ma);
+		}
 	}
 
 	/* and one render object per particle system */
@@ -4292,8 +4352,8 @@ static void init_render_object(Render *re, Object *ob, Object *par, DupliObject
 	else if(render_object_type(ob->type))
 		add_render_object(re, ob, par, dob, timeoffset, vectorlay);
 	else {
-		MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-		MTC_Mat4Invert(ob->imat, mat);
+		Mat4MulMat4(mat, ob->obmat, re->viewmat);
+		Mat4Invert(ob->imat, mat);
 	}
 	
 	time= PIL_check_seconds_timer();
@@ -4334,6 +4394,8 @@ void RE_Database_Free(Render *re)
 		curvemapping_free(lar->curfalloff);
 	}
 	
+	free_volume_precache(re);
+	
 	BLI_freelistN(&re->lampren);
 	BLI_freelistN(&re->lights);
 
@@ -4360,6 +4422,11 @@ void RE_Database_Free(Render *re)
 	end_render_materials();
 	end_render_textures();
 	
+	free_pointdensities(re);
+	free_voxeldata(re);
+	
+	free_camera_inside_volumes(re);
+	
 	if(re->wrld.aosphere) {
 		MEM_freeN(re->wrld.aosphere);
 		re->wrld.aosphere= NULL;
@@ -4536,8 +4603,8 @@ static void database_init_objects(Render *re, unsigned int renderlay, int nolamp
 	for(SETLOOPER(re->scene, base)) {
 		ob= base->object;
 		/* imat objects has to be done here, since displace can have texture using Object map-input */
-		MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-		MTC_Mat4Invert(ob->imat, mat);
+		Mat4MulMat4(mat, ob->obmat, re->viewmat);
+		Mat4Invert(ob->imat, mat);
 		/* each object should only be rendered once */
 		ob->flag &= ~OB_DONE;
 		ob->transflag &= ~OB_RENDER_DUPLI;
@@ -4683,8 +4750,8 @@ static void database_init_objects(Render *re, unsigned int renderlay, int nolamp
 	if(redoimat) {
 		for(SETLOOPER(re->scene, base)) {
 			ob= base->object;
-			MTC_Mat4MulMat4(mat, ob->obmat, re->viewmat);
-			MTC_Mat4Invert(ob->imat, mat);
+			Mat4MulMat4(mat, ob->obmat, re->viewmat);
+			Mat4Invert(ob->imat, mat);
 		}
 	}
 
@@ -4752,6 +4819,8 @@ void RE_Database_FromScene(Render *re, Scene *scene, int use_camera_view)
 
 	/* MAKE RENDER DATA */
 	database_init_objects(re, lay, 0, 0, 0, 0);
+	
+	init_camera_inside_volumes(re);
 
 	if(!re->test_break(re->tbh)) {
 		int tothalo;
@@ -4800,6 +4869,13 @@ void RE_Database_FromScene(Render *re, Scene *scene, int use_camera_view)
 			/* ENVIRONMENT MAPS */
 			if(!re->test_break(re->tbh))
 				make_envmaps(re);
+				
+			/* point density texture */
+			if(!re->test_break(re->tbh))
+				make_pointdensities(re);
+			/* voxel data texture */
+			if(!re->test_break(re->tbh))
+				make_voxeldata(re);
 		}
 		
 		if(!re->test_break(re->tbh))
@@ -4816,6 +4892,11 @@ void RE_Database_FromScene(Render *re, Scene *scene, int use_camera_view)
 		if((re->r.mode & R_SSS) && !re->test_break(re->tbh))
 			if(re->r.renderer==R_INTERN)
 				make_sss_tree(re);
+		
+		if(!re->test_break(re->tbh))
+			if(re->r.mode & R_RAYTRACE)
+				volume_precache(re);
+		
 	}
 	
 	if(re->test_break(re->tbh))
@@ -5096,7 +5177,7 @@ static int load_fluidsimspeedvectors(Render *re, ObjectInstanceRen *obi, float *
 		return 0;
 	
 	Mat4CpyMat4(mat, re->viewmat);
-	MTC_Mat4Invert(imat, mat);
+	Mat4Invert(imat, mat);
 
 	/* set first vertex OK */
 	if(!fss->meshSurfNormals) return 0;
diff --git a/source/blender/render/intern/source/envmap.c b/source/blender/render/intern/source/envmap.c
index a57e38f47c8..b5774d11799 100644
--- a/source/blender/render/intern/source/envmap.c
+++ b/source/blender/render/intern/source/envmap.c
@@ -51,7 +51,7 @@
 #include "BKE_texture.h"
 #include "BKE_utildefines.h"
 
-#include "MTC_matrixops.h"
+
 
 /* this module */
 #include "render_types.h"
@@ -211,9 +211,9 @@ static void envmap_transmatrix(float mat[][4], int part)
 		eul[2]= -M_PI/2.0;
 	}
 	
-	MTC_Mat4CpyMat4(tmat, mat);
+	Mat4CpyMat4(tmat, mat);
 	EulToMat4(eul, rotmat);
-	MTC_Mat4MulSerie(mat, tmat, rotmat,
+	Mat4MulSerie(mat, tmat, rotmat,
 					 0,   0,    0,
 					 0,   0,    0);
 }
@@ -231,12 +231,12 @@ static void env_rotate_scene(Render *re, float mat[][4], int mode)
 	int a;
 	
 	if(mode==0) {
-		MTC_Mat4Invert(tmat, mat);
-		MTC_Mat3CpyMat4(imat, tmat);
+		Mat4Invert(tmat, mat);
+		Mat3CpyMat4(imat, tmat);
 	}
 	else {
-		MTC_Mat4CpyMat4(tmat, mat);
-		MTC_Mat3CpyMat4(imat, mat);
+		Mat4CpyMat4(tmat, mat);
+		Mat3CpyMat4(imat, mat);
 	}
 
 	for(obi=re->instancetable.first; obi; obi=obi->next) {
@@ -267,7 +267,7 @@ static void env_rotate_scene(Render *re, float mat[][4], int mode)
 			if((a & 255)==0) har= obr->bloha[a>>8];
 			else har++;
 		
-			MTC_Mat4MulVecfl(tmat, har->co);
+			Mat4MulVecfl(tmat, har->co);
 		}
 	}
 	
@@ -280,22 +280,22 @@ static void env_rotate_scene(Render *re, float mat[][4], int mode)
 		Mat3CpyMat3(cmat, lar->imat); 
 		Mat3MulMat3(lar->imat, cmat, imat); 
 
-		MTC_Mat3MulVecfl(imat, lar->vec);
-		MTC_Mat4MulVecfl(tmat, lar->co);
+		Mat3MulVecfl(imat, lar->vec);
+		Mat4MulVecfl(tmat, lar->co);
 
 		lar->sh_invcampos[0]= -lar->co[0];
 		lar->sh_invcampos[1]= -lar->co[1];
 		lar->sh_invcampos[2]= -lar->co[2];
-		MTC_Mat3MulVecfl(lar->imat, lar->sh_invcampos);
+		Mat3MulVecfl(lar->imat, lar->sh_invcampos);
 		lar->sh_invcampos[2]*= lar->sh_zfac;
 		
 		if(lar->shb) {
 			if(mode==1) {
-				MTC_Mat4Invert(pmat, mat);
-				MTC_Mat4MulMat4(smat, pmat, lar->shb->viewmat);
-				MTC_Mat4MulMat4(lar->shb->persmat, smat, lar->shb->winmat);
+				Mat4Invert(pmat, mat);
+				Mat4MulMat4(smat, pmat, lar->shb->viewmat);
+				Mat4MulMat4(lar->shb->persmat, smat, lar->shb->winmat);
 			}
-			else MTC_Mat4MulMat4(lar->shb->persmat, lar->shb->viewmat, lar->shb->winmat);
+			else Mat4MulMat4(lar->shb->persmat, lar->shb->viewmat, lar->shb->winmat);
 		}
 	}
 	
@@ -373,8 +373,8 @@ static void env_set_imats(Render *re)
 	
 	base= re->scene->base.first;
 	while(base) {
-		MTC_Mat4MulMat4(mat, base->object->obmat, re->viewmat);
-		MTC_Mat4Invert(base->object->imat, mat);
+		Mat4MulMat4(mat, base->object->obmat, re->viewmat);
+		Mat4Invert(base->object->imat, mat);
 		
 		base= base->next;
 	}
@@ -393,18 +393,18 @@ static void render_envmap(Render *re, EnvMap *env)
 	short part;
 	
 	/* need a recalc: ortho-render has no correct viewinv */
-	MTC_Mat4Invert(oldviewinv, re->viewmat);
+	Mat4Invert(oldviewinv, re->viewmat);
 
 	envre= envmap_render_copy(re, env);
 	
 	/* precalc orthmat for object */
-	MTC_Mat4CpyMat4(orthmat, env->object->obmat);
-	MTC_Mat4Ortho(orthmat);
+	Mat4CpyMat4(orthmat, env->object->obmat);
+	Mat4Ortho(orthmat);
 	
 	/* need imat later for texture imat */
-	MTC_Mat4MulMat4(mat, orthmat, re->viewmat);
-	MTC_Mat4Invert(tmat, mat);
-	MTC_Mat3CpyMat4(env->obimat, tmat);
+	Mat4MulMat4(mat, orthmat, re->viewmat);
+	Mat4Invert(tmat, mat);
+	Mat3CpyMat4(env->obimat, tmat);
 
 	for(part=0; part<6; part++) {
 		if(env->type==ENV_PLANE && part!=1)
@@ -412,17 +412,17 @@ static void render_envmap(Render *re, EnvMap *env)
 		
 		re->display_clear(re->dch, envre->result);
 		
-		MTC_Mat4CpyMat4(tmat, orthmat);
+		Mat4CpyMat4(tmat, orthmat);
 		envmap_transmatrix(tmat, part);
-		MTC_Mat4Invert(mat, tmat);
+		Mat4Invert(mat, tmat);
 		/* mat now is the camera 'viewmat' */
 
-		MTC_Mat4CpyMat4(envre->viewmat, mat);
-		MTC_Mat4CpyMat4(envre->viewinv, tmat);
+		Mat4CpyMat4(envre->viewmat, mat);
+		Mat4CpyMat4(envre->viewinv, tmat);
 		
 		/* we have to correct for the already rotated vertexcoords */
-		MTC_Mat4MulMat4(tmat, oldviewinv, envre->viewmat);
-		MTC_Mat4Invert(env->imat, tmat);
+		Mat4MulMat4(tmat, oldviewinv, envre->viewmat);
+		Mat4Invert(env->imat, tmat);
 		
 		env_rotate_scene(envre, tmat, 1);
 		init_render_world(envre);
@@ -503,13 +503,13 @@ void make_envmaps(Render *re)
 							float orthmat[4][4], mat[4][4], tmat[4][4];
 							
 							/* precalc orthmat for object */
-							MTC_Mat4CpyMat4(orthmat, env->object->obmat);
-							MTC_Mat4Ortho(orthmat);
+							Mat4CpyMat4(orthmat, env->object->obmat);
+							Mat4Ortho(orthmat);
 							
 							/* need imat later for texture imat */
-							MTC_Mat4MulMat4(mat, orthmat, re->viewmat);
-							MTC_Mat4Invert(tmat, mat);
-							MTC_Mat3CpyMat4(env->obimat, tmat);
+							Mat4MulMat4(mat, orthmat, re->viewmat);
+							Mat4Invert(tmat, mat);
+							Mat3CpyMat4(env->obimat, tmat);
 						}
 						else {
 							
@@ -678,20 +678,20 @@ int envmaptex(Tex *tex, float *texvec, float *dxt, float *dyt, int osatex, TexRe
 	
 	/* rotate to envmap space, if object is set */
 	VECCOPY(vec, texvec);
-	if(env->object) MTC_Mat3MulVecfl(env->obimat, vec);
-	else MTC_Mat4Mul3Vecfl(R.viewinv, vec);
+	if(env->object) Mat3MulVecfl(env->obimat, vec);
+	else Mat4Mul3Vecfl(R.viewinv, vec);
 	
 	face= envcube_isect(env, vec, sco);
 	ibuf= env->cube[face];
 	
 	if(osatex) {
 		if(env->object) {
-			MTC_Mat3MulVecfl(env->obimat, dxt);
-			MTC_Mat3MulVecfl(env->obimat, dyt);
+			Mat3MulVecfl(env->obimat, dxt);
+			Mat3MulVecfl(env->obimat, dyt);
 		}
 		else {
-			MTC_Mat4Mul3Vecfl(R.viewinv, dxt);
-			MTC_Mat4Mul3Vecfl(R.viewinv, dyt);
+			Mat4Mul3Vecfl(R.viewinv, dxt);
+			Mat4Mul3Vecfl(R.viewinv, dyt);
 		}
 		set_dxtdyt(dxts, dyts, dxt, dyt, face);
 		imagewraposa(tex, NULL, ibuf, sco, dxts, dyts, texres);
diff --git a/source/blender/render/intern/source/initrender.c b/source/blender/render/intern/source/initrender.c
index 842adcf8520..d388e81a745 100644
--- a/source/blender/render/intern/source/initrender.c
+++ b/source/blender/render/intern/source/initrender.c
@@ -41,7 +41,7 @@
 #include "BLI_blenlib.h"
 #include "BLI_jitter.h"
 
-#include "MTC_matrixops.h"
+
 
 #include "DNA_camera_types.h"
 #include "DNA_group_types.h"
diff --git a/source/blender/render/intern/source/occlusion.c b/source/blender/render/intern/source/occlusion.c
index feef3dd424a..a15377a8c6d 100644
--- a/source/blender/render/intern/source/occlusion.c
+++ b/source/blender/render/intern/source/occlusion.c
@@ -630,7 +630,7 @@ static OcclusionTree *occ_tree_build(Render *re)
 			if((a & 255)==0) vlr= obr->vlaknodes[a>>8].vlak;
 			else vlr++;
 
-			if(vlr->mat->mode & MA_TRACEBLE)
+			if((vlr->mat->mode & MA_TRACEBLE) && (vlr->mat->material_type == MA_TYPE_SURFACE))
 				totface++;
 		}
 	}
@@ -663,7 +663,7 @@ static OcclusionTree *occ_tree_build(Render *re)
 			if((a & 255)==0) vlr= obr->vlaknodes[a>>8].vlak;
 			else vlr++;
 
-			if(vlr->mat->mode & MA_TRACEBLE) {
+			if((vlr->mat->mode & MA_TRACEBLE) && (vlr->mat->material_type == MA_TYPE_SURFACE)) {
 				tree->face[b].obi= c;
 				tree->face[b].facenr= a;
 				tree->occlusion[b]= 1.0f;
diff --git a/source/blender/render/intern/source/pixelshading.c b/source/blender/render/intern/source/pixelshading.c
index 75a2ab257f4..de3a50acddf 100644
--- a/source/blender/render/intern/source/pixelshading.c
+++ b/source/blender/render/intern/source/pixelshading.c
@@ -32,8 +32,8 @@
 /* External modules: */
 #include "IMB_imbuf_types.h"
 #include "IMB_imbuf.h"
-#include "MTC_matrixops.h"
-#include "MTC_vectorops.h"
+
+
 
 #include "DNA_camera_types.h"
 #include "DNA_group_types.h"
@@ -155,7 +155,7 @@ static void render_lighting_halo(HaloRen *har, float *colf)
 					
 					/* rotate view to lampspace */
 					VECCOPY(lvrot, lv);
-					MTC_Mat3MulVecfl(lar->imat, lvrot);
+					Mat3MulVecfl(lar->imat, lvrot);
 					
 					x= MAX2(fabs(lvrot[0]/lvrot[2]) , fabs(lvrot[1]/lvrot[2]));
 					/* 1.0/(sqrt(1+x*x)) is equivalent to cos(atan(x)) */
@@ -553,7 +553,7 @@ void shadeSkyView(float *colf, float *rco, float *view, float *dxyview, short th
 		VECCOPY(lo, view);
 		if(R.wrld.skytype & WO_SKYREAL) {
 			
-			MTC_Mat3MulVecfl(R.imat, lo);
+			Mat3MulVecfl(R.imat, lo);
 			
 			SWAP(float, lo[1],  lo[2]);
 			
@@ -595,7 +595,7 @@ void shadeSunView(float *colf, float *view)
 				
 				VECCOPY(sview, view);
 				Normalize(sview);
-				MTC_Mat3MulVecfl(R.imat, sview);
+				Mat3MulVecfl(R.imat, sview);
 				if (sview[2] < 0.0)
 					sview[2] = 0.0;
 				Normalize(sview);
@@ -678,7 +678,7 @@ void shadeAtmPixel(struct SunSky *sunsky, float *collector, float fx, float fy,
 		
 	calc_view_vector(view, fx, fy);
 	Normalize(view);
-	/*MTC_Mat3MulVecfl(R.imat, view);*/
+	/*Mat3MulVecfl(R.imat, view);*/
 	AtmospherePixleShader(sunsky, view, distance, collector);
 }
 
diff --git a/source/blender/render/intern/source/pointdensity.c b/source/blender/render/intern/source/pointdensity.c
new file mode 100644
index 00000000000..5f8cf5504fa
--- /dev/null
+++ b/source/blender/render/intern/source/pointdensity.c
@@ -0,0 +1,484 @@
+/* 
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * Contributors: Matt Ebb
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_arithb.h"
+#include "BLI_blenlib.h"
+#include "BLI_kdopbvh.h"
+
+#include "BKE_DerivedMesh.h"
+#include "BKE_global.h"
+#include "BKE_lattice.h"
+#include "BKE_main.h"
+#include "BKE_object.h"
+#include "BKE_particle.h"
+#include "BKE_texture.h"
+
+#include "DNA_texture_types.h"
+#include "DNA_particle_types.h"
+
+#include "render_types.h"
+#include "renderdatabase.h"
+#include "texture.h"
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
+/* only to be used here in this file, it's for speed */
+extern struct Render R;
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+
+
+static int point_data_used(PointDensity *pd)
+{
+	int pd_bitflag = 0;
+	
+	if ((pd->noise_influence == TEX_PD_NOISE_VEL) || (pd->color_source == TEX_PD_COLOR_PARTVEL) || (pd->color_source == TEX_PD_COLOR_PARTSPEED))
+		pd_bitflag |= POINT_DATA_VEL;
+	if ((pd->noise_influence == TEX_PD_NOISE_AGE) || (pd->color_source == TEX_PD_COLOR_PARTAGE)) 
+		pd_bitflag |= POINT_DATA_LIFE;
+		
+	return pd_bitflag;
+}
+
+
+/* additional data stored alongside the point density BVH, 
+ * accessible by point index number to retrieve other information 
+ * such as particle velocity or lifetime */
+static void alloc_point_data(PointDensity *pd, int total_particles, int point_data_used)
+{
+	int data_size = 0;
+	
+	if (point_data_used & POINT_DATA_VEL) {
+		/* store 3 channels of velocity data */
+		data_size += 3;
+	}
+	if (point_data_used & POINT_DATA_LIFE) {
+		/* store 1 channel of lifetime data */
+		data_size += 1;
+	}
+
+	if (data_size)
+		pd->point_data = MEM_mallocN(sizeof(float)*data_size*total_particles, "particle point data");
+}
+
+static void pointdensity_cache_psys(Render *re, PointDensity *pd, Object *ob, ParticleSystem *psys)
+{
+	DerivedMesh* dm;
+	ParticleKey state;
+	ParticleData *pa=NULL;
+	float cfra = bsystem_time(re->scene, ob, (float)re->scene->r.cfra, 0.0);
+	int i, childexists;
+	int total_particles, offset=0;
+	int data_used = point_data_used(pd);
+	float partco[3];
+	float obview[4][4];
+	
+	
+	/* init everything */
+	if (!psys || !ob || !pd) return;
+	
+	Mat4MulMat4(obview, re->viewinv, ob->obmat);
+	
+	/* Just to create a valid rendering context for particles */
+	psys_render_set(ob, psys, re->viewmat, re->winmat, re->winx, re->winy, 0);
+	
+	dm = mesh_create_derived_render(re->scene, ob,CD_MASK_BAREMESH|CD_MASK_MTFACE|CD_MASK_MCOL);
+	
+	if ( !psys_check_enabled(ob, psys)) {
+		psys_render_restore(ob, psys);
+		return;
+	}
+	
+	/* in case ob->imat isn't up-to-date */
+	Mat4Invert(ob->imat, ob->obmat);
+	
+	total_particles = psys->totpart+psys->totchild;
+	psys->lattice=psys_get_lattice(re->scene,ob,psys);
+	
+	pd->point_tree = BLI_bvhtree_new(total_particles, 0.0, 4, 6);
+	alloc_point_data(pd, total_particles, data_used);
+	pd->totpoints = total_particles;
+	if (data_used & POINT_DATA_VEL) offset = pd->totpoints*3;
+	
+	if (psys->totchild > 0 && !(psys->part->draw & PART_DRAW_PARENT))
+		childexists = 1;
+	
+	for (i=0, pa=psys->particles; i < total_particles; i++, pa++) {
+
+		state.time = cfra;
+		if(psys_get_particle_state(re->scene, ob, psys, i, &state, 0)) {
+			
+			VECCOPY(partco, state.co);
+			
+			if (pd->psys_cache_space == TEX_PD_OBJECTSPACE)
+				Mat4MulVecfl(ob->imat, partco);
+			else if (pd->psys_cache_space == TEX_PD_OBJECTLOC) {
+				float obloc[3];
+				VECCOPY(obloc, ob->loc);
+				VecSubf(partco, partco, obloc);
+			} else {
+				/* TEX_PD_WORLDSPACE */
+			}
+			
+			BLI_bvhtree_insert(pd->point_tree, i, partco, 1);
+			
+			if (data_used & POINT_DATA_VEL) {
+				pd->point_data[i*3 + 0] = state.vel[0];
+				pd->point_data[i*3 + 1] = state.vel[1];
+				pd->point_data[i*3 + 2] = state.vel[2];
+			} 
+			if (data_used & POINT_DATA_LIFE) {
+				float pa_time;
+				
+				if (i < psys->totpart) {
+					pa_time = (cfra - pa->time)/pa->lifetime;
+				} else {
+					ChildParticle *cpa= (psys->child + i) - psys->totpart;
+					float pa_birthtime, pa_dietime;
+					
+					pa_time = psys_get_child_time(psys, cpa, cfra, &pa_birthtime, &pa_dietime);
+				}
+				
+				pd->point_data[offset + i] = pa_time;
+			}
+		}
+	}
+	
+	BLI_bvhtree_balance(pd->point_tree);
+	dm->release(dm);
+	
+	if(psys->lattice){
+		end_latt_deform(psys->lattice);
+		psys->lattice=0;
+	}
+	
+	psys_render_restore(ob, psys);
+}
+
+
+static void pointdensity_cache_object(Render *re, PointDensity *pd, ObjectRen *obr)
+{
+	int i;
+	
+	if (!obr || !pd) return;
+	if(!obr->vertnodes) return;
+	
+	/* in case ob->imat isn't up-to-date */
+	Mat4Invert(obr->ob->imat, obr->ob->obmat);
+	
+	pd->point_tree = BLI_bvhtree_new(obr->totvert, 0.0, 4, 6);
+	pd->totpoints = obr->totvert;
+	
+	for(i=0; i<obr->totvert; i++) {
+		float ver_co[3];
+		VertRen *ver= RE_findOrAddVert(obr, i);
+		
+		VECCOPY(ver_co, ver->co);
+		Mat4MulVecfl(re->viewinv, ver_co);
+		
+		if (pd->ob_cache_space == TEX_PD_OBJECTSPACE) {
+			Mat4MulVecfl(obr->ob->imat, ver_co);
+		} else if (pd->psys_cache_space == TEX_PD_OBJECTLOC) {
+			VecSubf(ver_co, ver_co, obr->ob->loc);
+		} else {
+			/* TEX_PD_WORLDSPACE */
+		}
+		
+		BLI_bvhtree_insert(pd->point_tree, i, ver_co, 1);
+	}
+	
+	BLI_bvhtree_balance(pd->point_tree);
+
+}
+static void cache_pointdensity(Render *re, Tex *tex)
+{
+	PointDensity *pd = tex->pd;
+
+	if (pd->point_tree) {
+		BLI_bvhtree_free(pd->point_tree);
+		pd->point_tree = NULL;
+	}
+	
+	if (pd->source == TEX_PD_PSYS) {
+		Object *ob = pd->object;
+
+		if (!ob) return;
+		if (!pd->psys) return;
+		
+		pointdensity_cache_psys(re, pd, ob, pd->psys);
+	}
+	else if (pd->source == TEX_PD_OBJECT) {
+		Object *ob = pd->object;
+		ObjectRen *obr;
+		int found=0;
+
+		/* find the obren that corresponds to the object */
+		for (obr=re->objecttable.first; obr; obr=obr->next) {
+			if (obr->ob == ob) {
+				found=1;
+				break;
+			}
+		}
+		if (!found) return;
+		
+		pointdensity_cache_object(re, pd, obr);
+	}
+}
+
+static void free_pointdensity(Render *re, Tex *tex)
+{
+	PointDensity *pd = tex->pd;
+
+	if (!pd) return;
+	
+	if (pd->point_tree) {
+		BLI_bvhtree_free(pd->point_tree);
+		pd->point_tree = NULL;
+	}
+
+	if (pd->point_data) {
+		MEM_freeN(pd->point_data);
+		pd->point_data = NULL;
+	}
+	pd->totpoints = 0;
+}
+
+
+
+void make_pointdensities(Render *re)
+{
+	Tex *tex;
+	
+	if(re->scene->r.scemode & R_PREVIEWBUTS)
+		return;
+	
+	re->i.infostr= "Caching Point Densities";
+	re->stats_draw(re->sdh, &re->i);
+
+	for (tex= G.main->tex.first; tex; tex= tex->id.next) {
+		if(tex->id.us && tex->type==TEX_POINTDENSITY) {
+			cache_pointdensity(re, tex);
+		}
+	}
+	
+	re->i.infostr= NULL;
+	re->stats_draw(re->sdh, &re->i);
+}
+
+void free_pointdensities(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_POINTDENSITY) {
+			free_pointdensity(re, tex);
+		}
+	}
+}
+
+typedef struct PointDensityRangeData
+{
+    float *density;
+    float squared_radius;
+    float *point_data;
+	float *vec;
+	float softness;
+    short falloff_type;
+	short noise_influence;
+	float *age;
+	int point_data_used;
+	int offset;
+} PointDensityRangeData;
+
+void accum_density(void *userdata, int index, float squared_dist)
+{
+	PointDensityRangeData *pdr = (PointDensityRangeData *)userdata;
+	const float dist = (pdr->squared_radius - squared_dist) / pdr->squared_radius * 0.5f;
+	float density = 0.0f;
+	
+	if (pdr->falloff_type == TEX_PD_FALLOFF_STD)
+		density = dist;
+	else if (pdr->falloff_type == TEX_PD_FALLOFF_SMOOTH)
+		density = 3.0f*dist*dist - 2.0f*dist*dist*dist;
+	else if (pdr->falloff_type == TEX_PD_FALLOFF_SOFT)
+		density = pow(dist, pdr->softness);
+	else if (pdr->falloff_type == TEX_PD_FALLOFF_CONSTANT)
+		density = pdr->squared_radius;
+	else if (pdr->falloff_type == TEX_PD_FALLOFF_ROOT)
+		density = sqrt(dist);
+	
+	if (pdr->point_data_used & POINT_DATA_VEL) {
+		pdr->vec[0] += pdr->point_data[index*3 + 0]; //* density;
+		pdr->vec[1] += pdr->point_data[index*3 + 1]; //* density;
+		pdr->vec[2] += pdr->point_data[index*3 + 2]; //* density;
+	}
+	if (pdr->point_data_used & POINT_DATA_LIFE) {
+		*pdr->age += pdr->point_data[pdr->offset + index]; // * density;
+	}
+	
+	*pdr->density += density;
+}
+
+
+static void init_pointdensityrangedata(PointDensity *pd, PointDensityRangeData *pdr, float *density, float *vec, float *age)
+{
+	pdr->squared_radius = pd->radius*pd->radius;
+	pdr->density = density;
+	pdr->point_data = pd->point_data;
+	pdr->falloff_type = pd->falloff_type;
+	pdr->vec = vec;
+	pdr->age = age;
+	pdr->softness = pd->falloff_softness;
+	pdr->noise_influence = pd->noise_influence;
+	pdr->point_data_used = point_data_used(pd);
+	pdr->offset = (pdr->point_data_used & POINT_DATA_VEL)?pd->totpoints*3:0;
+}
+
+
+int pointdensitytex(Tex *tex, float *texvec, TexResult *texres)
+{
+	int retval = TEX_INT;
+	PointDensity *pd = tex->pd;
+	PointDensityRangeData pdr;
+	float density=0.0f, age=0.0f, time=0.0f;
+	float vec[3] = {0.0f, 0.0f, 0.0f}, co[3];
+	float col[4];
+	float turb, noise_fac;
+	int num=0;
+	
+	texres->tin = 0.0f;
+	
+	if ((!pd) || (!pd->point_tree))		
+		return 0;
+		
+	init_pointdensityrangedata(pd, &pdr, &density, vec, &age);
+	noise_fac = pd->noise_fac * 0.5f;	/* better default */
+	
+	VECCOPY(co, texvec);
+	
+	if (point_data_used(pd)) {
+		/* does a BVH lookup to find accumulated density and additional point data *
+		 * stores particle velocity vector in 'vec', and particle lifetime in 'time' */
+		num = BLI_bvhtree_range_query(pd->point_tree, co, pd->radius, accum_density, &pdr);
+		if (num > 0) {
+			age /= num;
+			VecMulf(vec, 1.0f/num);
+		}
+		
+		/* reset */
+		density = vec[0] = vec[1] = vec[2] = 0.0f;
+	}
+	
+	if (pd->flag & TEX_PD_TURBULENCE) {
+	
+		if (pd->noise_influence == TEX_PD_NOISE_AGE) {
+			turb = BLI_gTurbulence(pd->noise_size, texvec[0]+age, texvec[1]+age, texvec[2]+age, pd->noise_depth, 0, pd->noise_basis);
+		}
+		else if (pd->noise_influence == TEX_PD_NOISE_TIME) {
+			time = R.cfra / (float)R.r.efra;
+			turb = BLI_gTurbulence(pd->noise_size, texvec[0]+time, texvec[1]+time, texvec[2]+time, pd->noise_depth, 0, pd->noise_basis);
+			//turb = BLI_turbulence(pd->noise_size, texvec[0]+time, texvec[1]+time, texvec[2]+time, pd->noise_depth);
+		}
+		else {
+			turb = BLI_gTurbulence(pd->noise_size, texvec[0]+vec[0], texvec[1]+vec[1], texvec[2]+vec[2], pd->noise_depth, 0, pd->noise_basis);
+		}
+
+		turb -= 0.5f;	/* re-center 0.0-1.0 range around 0 to prevent offsetting result */
+		
+		/* now we have an offset coordinate to use for the density lookup */
+		co[0] = texvec[0] + noise_fac * turb;
+		co[1] = texvec[1] + noise_fac * turb;
+		co[2] = texvec[2] + noise_fac * turb;
+	}
+
+	/* BVH query with the potentially perturbed coordinates */
+	num = BLI_bvhtree_range_query(pd->point_tree, co, pd->radius, accum_density, &pdr);
+	if (num > 0) {
+		age /= num;
+		VecMulf(vec, 1.0f/num);
+	}
+	
+	texres->tin = density;
+	BRICONT;
+	
+	if (pd->color_source == TEX_PD_COLOR_CONSTANT)
+		return retval;
+	
+	retval |= TEX_RGB;
+	
+	switch (pd->color_source) {
+		case TEX_PD_COLOR_PARTAGE:
+			if (pd->coba) {
+				if (do_colorband(pd->coba, age, col)) {
+					texres->talpha= 1;
+					VECCOPY(&texres->tr, col);
+					texres->tin *= col[3];
+					texres->ta = texres->tin;
+				}
+			}
+			break;
+		case TEX_PD_COLOR_PARTSPEED:
+		{
+			float speed = VecLength(vec) * pd->speed_scale;
+			
+			if (pd->coba) {
+				if (do_colorband(pd->coba, speed, col)) {
+					texres->talpha= 1;	
+					VECCOPY(&texres->tr, col);
+					texres->tin *= col[3];
+					texres->ta = texres->tin;
+				}
+			}
+			break;
+		}
+		case TEX_PD_COLOR_PARTVEL:
+			texres->talpha= 1;
+			VecMulf(vec, pd->speed_scale);
+			VECCOPY(&texres->tr, vec);
+			texres->ta = texres->tin;
+			break;
+		case TEX_PD_COLOR_CONSTANT:
+		default:
+			texres->tr = texres->tg = texres->tb = texres->ta = 1.0f;
+			break;
+	}
+	BRICONTRGB;
+	
+	return retval;
+	
+	/*
+	if (texres->nor!=NULL) {
+		texres->nor[0] = texres->nor[1] = texres->nor[2] = 0.0f;
+	}
+	*/
+}
diff --git a/source/blender/render/intern/source/rayobject.c b/source/blender/render/intern/source/rayobject.c
index 05308361fdb..4bd8a12aa01 100644
--- a/source/blender/render/intern/source/rayobject.c
+++ b/source/blender/render/intern/source/rayobject.c
@@ -156,6 +156,15 @@ static int vlr_check_intersect(Isect *is, ObjectInstanceRen *obi, VlakRen *vlr)
 		return (is->lay & obi->lay);
 }
 
+static int vlr_check_intersect_solid(Isect *is, ObjectInstanceRen* obi, VlakRen *vlr)
+{
+	/* solid material types only */
+	if (vlr->mat->material_type == MA_TYPE_SURFACE)
+		return 1;
+	else
+		return 0;
+}
+
 /* ray - triangle or quad intersection */
 /* this function shall only modify Isect if it detects an hit */
 static int intersect_rayface(RayFace *face, Isect *is)
@@ -174,6 +183,11 @@ static int intersect_rayface(RayFace *face, Isect *is)
 		if(vlr_check_intersect(is, (ObjectInstanceRen*)face->ob, (VlakRen*)face->face ) == 0)
 			return 0;
 	}
+	if(is->skip & RE_SKIP_VLR_NON_SOLID_MATERIAL)
+	{
+		if(vlr_check_intersect_solid(is, (ObjectInstanceRen*)face->ob, (VlakRen*)face->face) == 0)
+			return 0;
+	}
 
 	RE_RC_COUNT(is->raycounter->faces.test);
 
diff --git a/source/blender/render/intern/source/rayshade.c b/source/blender/render/intern/source/rayshade.c
index 7e2cb4c6992..8cb11d11762 100644
--- a/source/blender/render/intern/source/rayshade.c
+++ b/source/blender/render/intern/source/rayshade.c
@@ -55,6 +55,7 @@
 #include "pixelshading.h"
 #include "shading.h"
 #include "texture.h"
+#include "volumetric.h"
 
 #include "RE_raytrace.h"
 #include "rayobject.h"
@@ -117,6 +118,7 @@ RayObject *  RE_rayobject_create(int type, int size)
 RayCounter re_rc_counter[BLENDER_MAX_THREADS] = {};
 #endif
 
+
 void freeraytree(Render *re)
 {
 	ObjectInstanceRen *obi;
@@ -359,9 +361,7 @@ void makeraytree(Render *re)
 	re->i.infostr= "Raytree finished";
 }
 
-
-
-static void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
+void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
 {
 	ObjectInstanceRen *obi= (ObjectInstanceRen*)is->hit.ob;
 	VlakRen *vlr= (VlakRen*)is->hit.face;
@@ -416,8 +416,14 @@ static void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
 		shade_input_flip_normals(shi);
 
 	shade_input_set_shade_texco(shi);
-	
-	if(is->mode==RE_RAY_SHADOW_TRA) {
+	if (shi->mat->material_type == MA_TYPE_VOLUME) {
+		if(ELEM(is->mode, RE_RAY_SHADOW, RE_RAY_SHADOW_TRA)) {
+			shade_volume_shadow(shi, shr, is);
+		} else {
+			shade_volume_outside(shi, shr);
+		}
+	}
+	else if(is->mode==RE_RAY_SHADOW_TRA) {
 		/* temp hack to prevent recursion */
 		if(shi->nodes==0 && shi->mat->nodetree && shi->mat->use_nodes) {
 			ntreeShaderExecTree(shi->mat->nodetree, shi, shr);
@@ -431,9 +437,20 @@ static void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
 			ntreeShaderExecTree(shi->mat->nodetree, shi, shr);
 			shi->mat= vlr->mat;		/* shi->mat is being set in nodetree */
 		}
-		else
-			shade_material_loop(shi, shr);
-		
+		else {
+			int tempdepth;
+			/* XXX dodgy business here, set ray depth to -1
+			 * to ignore raytrace in shade_material_loop()
+			 * this could really use a refactor --Matt */
+			if (shi->volume_depth == 0) {
+				tempdepth = shi->depth;
+				shi->depth = -1;
+				shade_material_loop(shi, shr);
+				shi->depth = tempdepth;
+			} else {
+				shade_material_loop(shi, shr);
+			}
+		}
 		/* raytrace likes to separate the spec color */
 		VECSUB(shr->diff, shr->combined, shr->spec);
 	}	
@@ -1391,15 +1408,20 @@ void ray_trace(ShadeInput *shi, ShadeResult *shr)
 			}
 			
 			if(shi->combinedflag & SCE_PASS_REFLECT) {
+				/* values in shr->spec can be greater then 1.0.
+				 * In this case the mircol uses a zero blending factor, so ignoring it is ok.
+				 * Fixes bug #18837 - when the spec is higher then 1.0,
+				 * diff can become a negative color - Campbell  */
 				
-				f= fr*(1.0f-shr->spec[0]);	f1= 1.0f-i;
-				diff[0]= f*mircol[0] + f1*diff[0];
+				f1= 1.0f-i;
 				
-				f= fg*(1.0f-shr->spec[1]);	f1= 1.0f-i;
-				diff[1]= f*mircol[1] + f1*diff[1];
+				diff[0] *= f1;
+				diff[1] *= f1;
+				diff[2] *= f1;
 				
-				f= fb*(1.0f-shr->spec[2]);	f1= 1.0f-i;
-				diff[2]= f*mircol[2] + f1*diff[2];
+				if(shr->spec[0]<1.0f)	diff[0] += mircol[0] * (fr*(1.0f-shr->spec[0]));
+				if(shr->spec[1]<1.0f)	diff[1] += mircol[1] * (fg*(1.0f-shr->spec[1]));
+				if(shr->spec[2]<1.0f)	diff[2] += mircol[2] * (fb*(1.0f-shr->spec[2]));
 			}
 		}
 	}
@@ -1456,11 +1478,15 @@ static void ray_trace_shadow_tra(Isect *is, ShadeInput *origshi, int depth, int
 		shi.nodes= origshi->nodes;
 		
 		shade_ray(is, &shi, &shr);
-		if (traflag & RAY_TRA)
-			d= shade_by_transmission(is, &shi, &shr);
-		
-		/* mix colors based on shadfac (rgb + amount of light factor) */
-		addAlphaLight(is->col, shr.diff, shr.alpha, d*shi.mat->filter);
+		if (shi.mat->material_type == MA_TYPE_SURFACE) {
+			if (traflag & RAY_TRA)
+				d= shade_by_transmission(is, &shi, &shr);
+			
+			/* 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);
+		}
 		
 		if(depth>0 && is->col[3]>0.0f) {
 			
@@ -1712,7 +1738,7 @@ static void ray_ao_qmc(ShadeInput *shi, float *shadfac)
 	RE_RC_INIT(isec, *shi);
 	isec.orig.ob   = shi->obi;
 	isec.orig.face = shi->vlr;
-	isec.skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK;
+	isec.skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK | RE_SKIP_VLR_NON_SOLID_MATERIAL;
 	isec.hint = 0;
 
 	isec.hit.ob   = 0;
diff --git a/source/blender/render/intern/source/rendercore.c b/source/blender/render/intern/source/rendercore.c
index 94d7651a03c..f3db64295a3 100644
--- a/source/blender/render/intern/source/rendercore.c
+++ b/source/blender/render/intern/source/rendercore.c
@@ -2246,7 +2246,7 @@ static int bake_check_intersect(Isect *is, int ob, RayFace *face)
 	/* no direction checking for now, doesn't always improve the result
 	 * (INPR(shi->facenor, bs->dir) > 0.0f); */
 
-	return (R.objectinstance[ob].obr->ob != bs->actob);
+	return (R.objectinstance[ob & ~RE_RAY_TRANSFORM_OFFS].obr->ob != bs->actob);
 }
 #endif
 
diff --git a/source/blender/render/intern/source/shadbuf.c b/source/blender/render/intern/source/shadbuf.c
index 33085b98095..48305d31e10 100644
--- a/source/blender/render/intern/source/shadbuf.c
+++ b/source/blender/render/intern/source/shadbuf.c
@@ -26,7 +26,7 @@
 #include <math.h>
 #include <string.h>
 
-#include "MTC_matrixops.h"
+
 #include "MEM_guardedalloc.h"
 
 #include "DNA_group_types.h"
@@ -403,7 +403,7 @@ void makeshadowbuf(Render *re, LampRen *lar)
 	wsize= shb->pixsize*(shb->size/2.0);
 	
 	i_window(-wsize, wsize, -wsize, wsize, shb->d, shb->clipend, shb->winmat);
-	MTC_Mat4MulMat4(shb->persmat, shb->viewmat, shb->winmat);
+	Mat4MulMat4(shb->persmat, shb->viewmat, shb->winmat);
 
 	if(ELEM(lar->buftype, LA_SHADBUF_REGULAR, LA_SHADBUF_HALFWAY)) {
 		/* jitter, weights - not threadsafe! */
@@ -673,7 +673,7 @@ float testshadowbuf(Render *re, ShadBuf *shb, float *rco, float *dxco, float *dy
 	VECCOPY(co, rco);
 	co[3]= 1.0f;
 
-	MTC_Mat4MulVec4fl(shb->persmat, co);	/* rational hom co */
+	Mat4MulVec4fl(shb->persmat, co);	/* rational hom co */
 
 	xs1= siz*(1.0f+co[0]/co[3]);
 	ys1= siz*(1.0f+co[1]/co[3]);
@@ -714,7 +714,7 @@ float testshadowbuf(Render *re, ShadBuf *shb, float *rco, float *dxco, float *dy
 	co[1]= rco[1]+dxco[1];
 	co[2]= rco[2]+dxco[2];
 	co[3]= 1.0;
-	MTC_Mat4MulVec4fl(shb->persmat,co);     /* rational hom co */
+	Mat4MulVec4fl(shb->persmat,co);     /* rational hom co */
 	dx[0]= xs1- siz*(1.0+co[0]/co[3]);
 	dx[1]= ys1- siz*(1.0+co[1]/co[3]);
 	
@@ -722,7 +722,7 @@ float testshadowbuf(Render *re, ShadBuf *shb, float *rco, float *dxco, float *dy
 	co[1]= rco[1]+dyco[1];
 	co[2]= rco[2]+dyco[2];
 	co[3]= 1.0;
-	MTC_Mat4MulVec4fl(shb->persmat,co);     /* rational hom co */
+	Mat4MulVec4fl(shb->persmat,co);     /* rational hom co */
 	dy[0]= xs1- siz*(1.0+co[0]/co[3]);
 	dy[1]= ys1- siz*(1.0+co[1]/co[3]);
 	
@@ -858,7 +858,7 @@ float shadow_halo(LampRen *lar, float *p1, float *p2)
 	co[1]= p1[1];
 	co[2]= p1[2]/lar->sh_zfac;
 	co[3]= 1.0;
-	MTC_Mat4MulVec4fl(shb->winmat, co);	/* rational hom co */
+	Mat4MulVec4fl(shb->winmat, co);	/* rational hom co */
 	xf1= siz*(1.0+co[0]/co[3]);
 	yf1= siz*(1.0+co[1]/co[3]);
 	zf1= (co[2]/co[3]);
@@ -868,7 +868,7 @@ float shadow_halo(LampRen *lar, float *p1, float *p2)
 	co[1]= p2[1];
 	co[2]= p2[2]/lar->sh_zfac;
 	co[3]= 1.0;
-	MTC_Mat4MulVec4fl(shb->winmat, co);	/* rational hom co */
+	Mat4MulVec4fl(shb->winmat, co);	/* rational hom co */
 	xf2= siz*(1.0+co[0]/co[3]);
 	yf2= siz*(1.0+co[1]/co[3]);
 	zf2= (co[2]/co[3]);
@@ -1659,7 +1659,7 @@ static int viewpixel_to_lampbuf(ShadBuf *shb, ObjectInstanceRen *obi, VlakRen *v
 	}
 	
 	/* move 3d vector to lampbuf */
-	MTC_Mat4MulVec4fl(shb->persmat, hoco);	/* rational hom co */
+	Mat4MulVec4fl(shb->persmat, hoco);	/* rational hom co */
 	
 	/* clip We can test for -1.0/1.0 because of the properties of the
 	 * coordinate transformations. */
diff --git a/source/blender/render/intern/source/shadeinput.c b/source/blender/render/intern/source/shadeinput.c
index 496d187cf05..0d3c5772d8a 100644
--- a/source/blender/render/intern/source/shadeinput.c
+++ b/source/blender/render/intern/source/shadeinput.c
@@ -29,7 +29,7 @@
 #include <math.h>
 #include <string.h>
 
-#include "MTC_matrixops.h"
+
 #include "BLI_arithb.h"
 #include "BLI_blenlib.h"
 
@@ -53,6 +53,7 @@
 #include "shading.h"
 #include "strand.h"
 #include "texture.h"
+#include "volumetric.h"
 #include "zbuf.h"
 
 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
@@ -167,6 +168,11 @@ void shade_material_loop(ShadeInput *shi, ShadeResult *shr)
 			if((shi->layflag & SCE_LAY_SKY) && (R.r.alphamode==R_ADDSKY))
 				shr->alpha= 1.0f;
 	}	
+	
+	if(R.r.mode & R_RAYTRACE) {
+		if (R.render_volumes_inside.first)
+			shade_volume_inside(shi, shr);
+	}
 }
 
 
@@ -187,7 +193,12 @@ void shade_input_do_shade(ShadeInput *shi, ShadeResult *shr)
 		/* copy all relevant material vars, note, keep this synced with render_types.h */
 		shade_input_init_material(shi);
 		
-		shade_material_loop(shi, shr);
+		if (shi->mat->material_type == MA_TYPE_VOLUME) {
+			if(R.r.mode & R_RAYTRACE)
+				shade_volume_outside(shi, shr);
+		} else { /* MA_TYPE_SURFACE, MA_TYPE_WIRE */
+			shade_material_loop(shi, shr);
+		}
 	}
 	
 	/* copy additional passes */
@@ -214,11 +225,12 @@ void shade_input_do_shade(ShadeInput *shi, ShadeResult *shr)
 	if(shr->alpha!=1.0f || alpha!=1.0f) {
 		float fac= alpha*(shr->alpha);
 		shr->combined[3]= fac;
-		shr->combined[0]*= fac;
-		shr->combined[1]*= fac;
-		shr->combined[2]*= fac;
+		
+		if (shi->mat->material_type!= MA_TYPE_VOLUME)
+			VecMulf(shr->combined, fac);
 	}
-	else shr->combined[3]= 1.0f;
+	else
+		shr->combined[3]= 1.0f;
 	
 	/* add z */
 	shr->z= -shi->co[2];
@@ -462,13 +474,13 @@ void shade_input_set_strand_texco(ShadeInput *shi, StrandRen *strand, StrandVert
 
 		if(texco & TEXCO_GLOB) {
 			VECCOPY(shi->gl, shi->co);
-			MTC_Mat4MulVecfl(R.viewinv, shi->gl);
+			Mat4MulVecfl(R.viewinv, shi->gl);
 			
 			if(shi->osatex) {
 				VECCOPY(shi->dxgl, shi->dxco);
-				MTC_Mat3MulVecfl(R.imat, shi->dxco);
+				Mat3MulVecfl(R.imat, shi->dxco);
 				VECCOPY(shi->dygl, shi->dyco);
-				MTC_Mat3MulVecfl(R.imat, shi->dyco);
+				Mat3MulVecfl(R.imat, shi->dyco);
 			}
 		}
 
@@ -714,6 +726,10 @@ void shade_input_calc_viewco(ShadeInput *shi, float x, float y, float z, float *
 		}
 	}
 	
+	/* set camera coords - for scanline, it's always 0.0,0.0,0.0 (render is in camera space)
+	 * however for raytrace it can be different - the position of the last intersection */
+	shi->camera_co[0] = shi->camera_co[1] = shi->camera_co[2] = 0.0f;
+	
 	/* cannot normalize earlier, code above needs it at viewplane level */
 	Normalize(view);
 }
@@ -1021,15 +1037,15 @@ void shade_input_set_shade_texco(ShadeInput *shi)
 		
 		if(texco & TEXCO_GLOB) {
 			VECCOPY(shi->gl, shi->co);
-			MTC_Mat4MulVecfl(R.viewinv, shi->gl);
+			Mat4MulVecfl(R.viewinv, shi->gl);
 			if(shi->osatex) {
 				VECCOPY(shi->dxgl, shi->dxco);
 				// TXF: bug was here, but probably should be in convertblender.c, R.imat only valid if there is a world
-				//MTC_Mat3MulVecfl(R.imat, shi->dxco);
-				MTC_Mat4Mul3Vecfl(R.viewinv, shi->dxco);
+				//Mat3MulVecfl(R.imat, shi->dxco);
+				Mat4Mul3Vecfl(R.viewinv, shi->dxco);
 				VECCOPY(shi->dygl, shi->dyco);
-				//MTC_Mat3MulVecfl(R.imat, shi->dyco);
-				MTC_Mat4Mul3Vecfl(R.viewinv, shi->dyco);
+				//Mat3MulVecfl(R.imat, shi->dyco);
+				Mat4Mul3Vecfl(R.viewinv, shi->dyco);
 			}
 		}
 		
diff --git a/source/blender/render/intern/source/shadeoutput.c b/source/blender/render/intern/source/shadeoutput.c
index d5c8cf30b30..5e523199755 100644
--- a/source/blender/render/intern/source/shadeoutput.c
+++ b/source/blender/render/intern/source/shadeoutput.c
@@ -30,7 +30,7 @@
 #include <math.h>
 #include <string.h>
 
-#include "MTC_matrixops.h"
+
 #include "BLI_arithb.h"
 
 #include "BKE_colortools.h"
@@ -58,7 +58,7 @@
 extern struct Render R;
 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
 
-static ListBase *get_lights(ShadeInput *shi)
+ListBase *get_lights(ShadeInput *shi)
 {
 	
 	if(R.r.scemode & R_PREVIEWBUTS)
@@ -168,7 +168,7 @@ static void spothalo(struct LampRen *lar, ShadeInput *shi, float *intens)
 		p1[0]= shi->co[0]-lar->co[0];
 		p1[1]= shi->co[1]-lar->co[1];
 		p1[2]= -lar->co[2];
-		MTC_Mat3MulVecfl(lar->imat, p1);
+		Mat3MulVecfl(lar->imat, p1);
 		VECCOPY(npos, p1);	// npos is double!
 		
 		/* pre-scale */
@@ -180,7 +180,7 @@ static void spothalo(struct LampRen *lar, ShadeInput *shi, float *intens)
 	
 	/* rotate view */
 	VECCOPY(nray, shi->view);
-	MTC_Mat3MulVecd(lar->imat, nray);
+	Mat3MulVecd(lar->imat, nray);
 	
 	if(R.wrld.mode & WO_MIST) {
 		/* patchy... */
@@ -1143,7 +1143,7 @@ float lamp_get_visibility(LampRen *lar, float *co, float *lv, float *dist)
 							
 							/* rotate view to lampspace */
 							VECCOPY(lvrot, lv);
-							MTC_Mat3MulVecfl(lar->imat, lvrot);
+							Mat3MulVecfl(lar->imat, lvrot);
 							
 							x= MAX2(fabs(lvrot[0]/lvrot[2]) , fabs(lvrot[1]/lvrot[2]));
 							/* 1.0f/(sqrt(1+x*x)) is equivalent to cos(atan(x)) */
diff --git a/source/blender/render/intern/source/strand.c b/source/blender/render/intern/source/strand.c
index 6e2a0e966dc..d00076a80e8 100644
--- a/source/blender/render/intern/source/strand.c
+++ b/source/blender/render/intern/source/strand.c
@@ -122,14 +122,14 @@ void strand_eval_point(StrandSegment *sseg, StrandPoint *spoint)
 			spoint->dtstrandco *= 0.5f;
 	}
 	else {
-		set_four_ipo(t, data, type);
+		key_curve_position_weights(t, data, type);
 		spoint->co[0]= data[0]*p[0][0] + data[1]*p[1][0] + data[2]*p[2][0] + data[3]*p[3][0];
 		spoint->co[1]= data[0]*p[0][1] + data[1]*p[1][1] + data[2]*p[2][1] + data[3]*p[3][1];
 		spoint->co[2]= data[0]*p[0][2] + data[1]*p[1][2] + data[2]*p[2][2] + data[3]*p[3][2];
 		spoint->strandco= (1.0f-t)*sseg->v[1]->strandco + t*sseg->v[2]->strandco;
 	}
 
-	set_afgeleide_four_ipo(t, data, type);
+	key_curve_tangent_weights(t, data, type);
 	spoint->dtco[0]= data[0]*p[0][0] + data[1]*p[1][0] + data[2]*p[2][0] + data[3]*p[3][0];
 	spoint->dtco[1]= data[0]*p[0][1] + data[1]*p[1][1] + data[2]*p[2][1] + data[3]*p[3][1];
 	spoint->dtco[2]= data[0]*p[0][2] + data[1]*p[1][2] + data[2]*p[2][2] + data[3]*p[3][2];
diff --git a/source/blender/render/intern/source/texture.c b/source/blender/render/intern/source/texture.c
index bb491efdaba..2d2c01e0bf1 100644
--- a/source/blender/render/intern/source/texture.c
+++ b/source/blender/render/intern/source/texture.c
@@ -30,7 +30,7 @@
 #include <string.h>
 #include <math.h>
 
-#include "MTC_matrixops.h"
+
 
 #include "BLI_blenlib.h"
 #include "BLI_arithb.h"
@@ -65,6 +65,8 @@
 #include "BKE_ipo.h"
 
 #include "envmap.h"
+#include "pointdensity.h"
+#include "voxeldata.h"
 #include "renderpipeline.h"
 #include "render_types.h"
 #include "rendercore.h"
@@ -714,12 +716,12 @@ static float voronoiTex(Tex *tex, float *texvec, TexResult *texres)
 
 /* ------------------------------------------------------------------------- */
 
-static int evalnodes(Tex *tex, float *texvec, TexResult *texres, short thread, short which_output)
+static int evalnodes(Tex *tex, float *texvec, float *dxt, float *dyt, TexResult *texres, short thread, short which_output)
 {
 	short rv = TEX_INT;
 	bNodeTree *nodes = tex->nodetree;
 	
-	ntreeTexExecTree(nodes, texres, texvec, 0, thread, tex, which_output, R.r.cfra);
+	ntreeTexExecTree(nodes, texres, texvec, dxt, dyt, thread, tex, which_output, R.r.cfra);
 	
 	if(texres->nor) rv |= TEX_NOR;
 	rv |= TEX_RGB;
@@ -830,7 +832,7 @@ static int cubemap_glob(float *n, float x, float y, float z, float *adr1, float
 	else {
 		VECCOPY(nor, n);
 	}
-	MTC_Mat4Mul3Vecfl(R.viewinv, nor);
+	Mat4Mul3Vecfl(R.viewinv, nor);
 
 	x1= fabs(nor[0]);
 	y1= fabs(nor[1]);
@@ -923,7 +925,7 @@ static int cubemap_ob(Object *ob, float *n, float x, float y, float z, float *ad
 	if(n==NULL) return 0;
 	
 	VECCOPY(nor, n);
-	if(ob) MTC_Mat4Mul3Vecfl(ob->imat, nor);
+	if(ob) Mat4Mul3Vecfl(ob->imat, nor);
 	
 	x1= fabs(nor[0]);
 	y1= fabs(nor[1]);
@@ -1180,7 +1182,7 @@ static int multitex(Tex *tex, float *texvec, float *dxt, float *dyt, int osatex,
 	texres->talpha= 0;	/* is set when image texture returns alpha (considered premul) */
 	
 	if(tex->use_nodes && tex->nodetree) {
-		retval = evalnodes(tex, texvec, texres, thread, which_output);
+		retval = evalnodes(tex, texvec, dxt, dyt, texres, thread, which_output);
 	}
 	else
 	switch(tex->type) {
@@ -1262,6 +1264,13 @@ static int multitex(Tex *tex, float *texvec, float *dxt, float *dyt, int osatex,
 		
 		retval= mg_distNoiseTex(tex, tmpvec, texres);
 		break;
+	case TEX_POINTDENSITY:
+		retval= pointdensitytex(tex, texvec, texres);
+		break;
+	case TEX_VOXELDATA:
+		retval= voxeldatatex(tex, texvec, texres);  
+		break;
+
 	}
 
 	if (tex->flag & TEX_COLORBAND) {
@@ -1272,7 +1281,7 @@ static int multitex(Tex *tex, float *texvec, float *dxt, float *dyt, int osatex,
 			texres->tg= col[1];
 			texres->tb= col[2];
 			texres->ta= col[3];
-			retval |= 1;
+			retval |= TEX_RGB;
 		}
 	}
 	return retval;
@@ -1446,12 +1455,22 @@ void texture_rgb_blend(float *in, float *tex, float *out, float fact, float facg
 		VECCOPY(in, out);
 		ramp_blend(MA_RAMP_COLOR, in, in+1, in+2, fact, tex);
 		break;
+    case MTEX_SOFT_LIGHT: 
+        fact*= facg; 
+        VECCOPY(in, out); 
+        ramp_blend(MA_RAMP_SOFT, in, in+1, in+2, fact, tex); 
+        break; 
+    case MTEX_LIN_LIGHT: 
+        fact*= facg; 
+        VECCOPY(in, out); 
+        ramp_blend(MA_RAMP_LINEAR, in, in+1, in+2, fact, tex); 
+        break; 
 	}
 }
 
 float texture_value_blend(float tex, float out, float fact, float facg, int blendtype, int flip)
 {
-	float in=0.0, facm, col;
+	float in=0.0, facm, col, scf;
 	
 	fact*= facg;
 	facm= 1.0-fact;
@@ -1496,6 +1515,19 @@ float texture_value_blend(float tex, float out, float fact, float facg, int blen
 		col= fact*tex;
 		if(col > out) in= col; else in= out;
 		break;
+
+    case MTEX_SOFT_LIGHT: 
+        col= fact*tex; 
+        scf=1.0 - (1.0 - tex) * (1.0 - out); 
+        in= facm*out + fact * ((1.0 - out) * tex * out) + (out * scf); 
+        break;       
+
+    case MTEX_LIN_LIGHT: 
+        if (tex > 0.5) 
+            in = out + fact*(2*(tex - 0.5)); 
+        else 
+            in = out + fact*(2*tex - 1); 
+        break;
 	}
 	
 	return in;
@@ -1639,13 +1671,13 @@ void do_material_tex(ShadeInput *shi)
 					VECCOPY(tempvec, shi->co);
 					if(mtex->texflag & MTEX_OB_DUPLI_ORIG)
 						if(shi->obi && shi->obi->duplitexmat)
-							MTC_Mat4MulVecfl(shi->obi->duplitexmat, tempvec);
-					MTC_Mat4MulVecfl(ob->imat, tempvec);
+							Mat4MulVecfl(shi->obi->duplitexmat, tempvec);
+					Mat4MulVecfl(ob->imat, tempvec);
 					if(shi->osatex) {
 						VECCOPY(dxt, shi->dxco);
 						VECCOPY(dyt, shi->dyco);
-						MTC_Mat4Mul3Vecfl(ob->imat, dxt);
-						MTC_Mat4Mul3Vecfl(ob->imat, dyt);
+						Mat4Mul3Vecfl(ob->imat, dxt);
+						Mat4Mul3Vecfl(ob->imat, dyt);
 					}
 				}
 				else {
@@ -2248,6 +2280,187 @@ void do_material_tex(ShadeInput *shi)
 	}
 }
 
+
+void do_volume_tex(ShadeInput *shi, float *xyz, int mapto_flag, float *col, float *val)
+{
+	MTex *mtex;
+	Tex *tex;
+	TexResult texres= {0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0, NULL};
+	int tex_nr, rgbnor= 0;
+	float co[3], texvec[3];
+	float fact, stencilTin=1.0;
+	
+	if (R.r.scemode & R_NO_TEX) return;
+	/* here: test flag if there's a tex (todo) */
+	
+	for(tex_nr=0; tex_nr<MAX_MTEX; tex_nr++) {
+		/* separate tex switching */
+		if(shi->mat->septex & (1<<tex_nr)) continue;
+		
+		if(shi->mat->mtex[tex_nr]) {
+			mtex= shi->mat->mtex[tex_nr];
+			tex= mtex->tex;
+			if(tex==0) continue;
+			
+			/* only process if this texture is mapped 
+			 * to one that we're interested in */
+			if (!(mtex->mapto & mapto_flag)) continue;
+			
+			/* which coords */
+			if(mtex->texco==TEXCO_OBJECT) { 
+				Object *ob= mtex->object;
+				ob= mtex->object;
+				if(ob) {						
+					VECCOPY(co, xyz);	
+					if(mtex->texflag & MTEX_OB_DUPLI_ORIG) {
+						if(shi->obi && shi->obi->duplitexmat)
+							Mat4MulVecfl(shi->obi->duplitexmat, co);					
+					} 
+					Mat4MulVecfl(ob->imat, co);
+				}
+			}
+			/* not really orco, but 'local' */
+			else if(mtex->texco==TEXCO_ORCO) {
+				
+				if(mtex->texflag & MTEX_DUPLI_MAPTO) {
+					VECCOPY(co, shi->duplilo);
+				}
+				else {
+					Object *ob= shi->obi->ob;
+					VECCOPY(co, xyz);
+					Mat4MulVecfl(ob->imat, co);
+				}
+			}
+			else if(mtex->texco==TEXCO_GLOB) {							
+			   VECCOPY(co, xyz);
+			   Mat4MulVecfl(R.viewinv, co);
+			}
+			else continue;	// can happen when texco defines disappear and it renders old files
+
+			texres.nor= NULL;
+			
+			if(tex->type==TEX_IMAGE) {
+				continue;	/* not supported yet */				
+				//do_2d_mapping(mtex, texvec, NULL, NULL, dxt, dyt);
+			}
+			else {
+				/* placement */
+				if(mtex->projx) texvec[0]= mtex->size[0]*(co[mtex->projx-1]+mtex->ofs[0]);
+				else texvec[0]= mtex->size[0]*(mtex->ofs[0]);
+
+				if(mtex->projy) texvec[1]= mtex->size[1]*(co[mtex->projy-1]+mtex->ofs[1]);
+				else texvec[1]= mtex->size[1]*(mtex->ofs[1]);
+
+				if(mtex->projz) texvec[2]= mtex->size[2]*(co[mtex->projz-1]+mtex->ofs[2]);
+				else texvec[2]= mtex->size[2]*(mtex->ofs[2]);
+			}
+			
+			rgbnor= multitex(tex, texvec, NULL, NULL, 0, &texres, 0, mtex->which_output);	/* NULL = dxt/dyt, 0 = shi->osatex - not supported */
+			
+			/* texture output */
+
+			if( (rgbnor & TEX_RGB) && (mtex->texflag & MTEX_RGBTOINT)) {
+				texres.tin= (0.35*texres.tr+0.45*texres.tg+0.2*texres.tb);
+				rgbnor-= TEX_RGB;
+			}
+			if(mtex->texflag & MTEX_NEGATIVE) {
+				if(rgbnor & TEX_RGB) {
+					texres.tr= 1.0-texres.tr;
+					texres.tg= 1.0-texres.tg;
+					texres.tb= 1.0-texres.tb;
+				}
+				texres.tin= 1.0-texres.tin;
+			}
+			if(mtex->texflag & MTEX_STENCIL) {
+				if(rgbnor & TEX_RGB) {
+					fact= texres.ta;
+					texres.ta*= stencilTin;
+					stencilTin*= fact;
+				}
+				else {
+					fact= texres.tin;
+					texres.tin*= stencilTin;
+					stencilTin*= fact;
+				}
+			}
+			
+			
+			if((mapto_flag & (MAP_EMISSION_COL+MAP_ABSORPTION_COL)) && (mtex->mapto & (MAP_EMISSION_COL+MAP_ABSORPTION_COL))) {
+				float tcol[3], colfac;
+				
+				/* stencil maps on the texture control slider, not texture intensity value */
+				colfac= mtex->colfac*stencilTin;
+				
+				if((rgbnor & TEX_RGB)==0) {
+					tcol[0]= mtex->r;
+					tcol[1]= mtex->g;
+					tcol[2]= mtex->b;
+				} else {
+					tcol[0]=texres.tr;
+					tcol[1]=texres.tg;
+					tcol[2]=texres.tb;
+					if(texres.talpha)
+						texres.tin= texres.ta;
+				}
+				
+				/* inverse gamma correction */
+				if (R.r.color_mgt_flag & R_COLOR_MANAGEMENT) {
+					color_manage_linearize(tcol, tcol);
+				}
+				
+				/* used for emit */
+				if((mapto_flag & MAP_EMISSION_COL) && (mtex->mapto & MAP_EMISSION_COL)) {
+					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)) {
+					texture_rgb_blend(col, tcol, col, texres.tin, colfac, mtex->blendtype);
+				}
+			}
+			
+			if((mapto_flag & MAP_VARS) && (mtex->mapto & MAP_VARS)) {
+				/* stencil maps on the texture control slider, not texture intensity value */
+				float varfac= mtex->varfac*stencilTin;
+				
+				/* convert RGB to intensity if intensity info isn't provided */
+				if (!(rgbnor & TEX_INT)) {
+					if (rgbnor & TEX_RGB) {
+						if(texres.talpha) texres.tin= texres.ta;
+						else texres.tin= (0.35*texres.tr+0.45*texres.tg+0.2*texres.tb);
+					}
+				}
+				
+				if((mapto_flag & MAP_EMISSION) && (mtex->mapto & MAP_EMISSION)) {
+					int flip= mtex->maptoneg & MAP_EMISSION;
+
+					*val = texture_value_blend(mtex->def_var, *val, texres.tin, varfac, mtex->blendtype, flip);
+					if(*val<0.0) *val= 0.0;
+				}
+				if((mapto_flag & MAP_DENSITY) && (mtex->mapto & MAP_DENSITY)) {
+					int flip= mtex->maptoneg & MAP_DENSITY;
+
+					*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;
+					
+					*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;
+					
+					*val = texture_value_blend(mtex->def_var, *val, texres.tin, varfac, mtex->blendtype, flip);
+					CLAMP(*val, 0.0, 1.0);
+				}
+			}
+		}
+	}
+}
+
+
 /* ------------------------------------------------------------------------- */
 
 void do_halo_tex(HaloRen *har, float xn, float yn, float *colf)
@@ -2449,7 +2662,7 @@ void do_sky_tex(float *rco, float *lo, float *dxyview, float *hor, float *zen, f
 			case TEXCO_OBJECT:
 				if(mtex->object) {
 					VECCOPY(tempvec, lo);
-					MTC_Mat4MulVecfl(mtex->object->imat, tempvec);
+					Mat4MulVecfl(mtex->object->imat, tempvec);
 					co= tempvec;
 				}
 				break;
@@ -2457,16 +2670,16 @@ void do_sky_tex(float *rco, float *lo, float *dxyview, float *hor, float *zen, f
 			case TEXCO_GLOB:
 				if(rco) {
 					VECCOPY(tempvec, rco);
-					MTC_Mat4MulVecfl(R.viewinv, tempvec);
+					Mat4MulVecfl(R.viewinv, tempvec);
 					co= tempvec;
 				}
 				else
 					co= lo;
 				
 //				VECCOPY(shi->dxgl, shi->dxco);
-//				MTC_Mat3MulVecfl(R.imat, shi->dxco);
+//				Mat3MulVecfl(R.imat, shi->dxco);
 //				VECCOPY(shi->dygl, shi->dyco);
-//				MTC_Mat3MulVecfl(R.imat, shi->dyco);
+//				Mat3MulVecfl(R.imat, shi->dyco);
 				break;
 			}
 			
@@ -2593,12 +2806,12 @@ void do_lamp_tex(LampRen *la, float *lavec, ShadeInput *shi, float *colf, int ef
 					dx= dxt;
 					dy= dyt;
 					VECCOPY(tempvec, shi->co);
-					MTC_Mat4MulVecfl(ob->imat, tempvec);
+					Mat4MulVecfl(ob->imat, tempvec);
 					if(shi->osatex) {
 						VECCOPY(dxt, shi->dxco);
 						VECCOPY(dyt, shi->dyco);
-						MTC_Mat4Mul3Vecfl(ob->imat, dxt);
-						MTC_Mat4Mul3Vecfl(ob->imat, dyt);
+						Mat4Mul3Vecfl(ob->imat, dxt);
+						Mat4Mul3Vecfl(ob->imat, dyt);
 					}
 				}
 				else {
@@ -2609,12 +2822,12 @@ void do_lamp_tex(LampRen *la, float *lavec, ShadeInput *shi, float *colf, int ef
 			else if(mtex->texco==TEXCO_GLOB) {
 				co= shi->gl; dx= shi->dxco; dy= shi->dyco;
 				VECCOPY(shi->gl, shi->co);
-				MTC_Mat4MulVecfl(R.viewinv, shi->gl);
+				Mat4MulVecfl(R.viewinv, shi->gl);
 			}
 			else if(mtex->texco==TEXCO_VIEW) {
 				
 				VECCOPY(tempvec, lavec);
-				MTC_Mat3MulVecfl(la->imat, tempvec);
+				Mat3MulVecfl(la->imat, tempvec);
 				
 				if(la->type==LA_SPOT) {
 					tempvec[0]*= la->spottexfac;
@@ -2627,8 +2840,8 @@ void do_lamp_tex(LampRen *la, float *lavec, ShadeInput *shi, float *colf, int ef
 					VECCOPY(dxt, shi->dxlv);
 					VECCOPY(dyt, shi->dylv);
 					/* need some matrix conversion here? la->imat is a [3][3]  matrix!!! **/
-					MTC_Mat3MulVecfl(la->imat, dxt);
-					MTC_Mat3MulVecfl(la->imat, dyt);
+					Mat3MulVecfl(la->imat, dxt);
+					Mat3MulVecfl(la->imat, dyt);
 					
 					VecMulf(dxt, la->spottexfac);
 					VecMulf(dyt, la->spottexfac);
diff --git a/source/blender/render/intern/source/volume_precache.c b/source/blender/render/intern/source/volume_precache.c
new file mode 100644
index 00000000000..ca19f2580aa
--- /dev/null
+++ b/source/blender/render/intern/source/volume_precache.c
@@ -0,0 +1,759 @@
+/**
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Matt Ebb.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_blenlib.h"
+#include "BLI_arithb.h"
+#include "BLI_threads.h"
+#include "BLI_voxel.h"
+
+#include "PIL_time.h"
+
+#include "RE_shader_ext.h"
+#include "RE_raytrace.h"
+
+#include "DNA_material_types.h"
+
+#include "render_types.h"
+#include "rendercore.h"
+#include "renderdatabase.h"
+#include "volumetric.h"
+#include "volume_precache.h"
+
+#if defined( _MSC_VER ) && !defined( __cplusplus )
+# define inline __inline
+#endif // defined( _MSC_VER ) && !defined( __cplusplus )
+
+#include "BKE_global.h"
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
+/* only to be used here in this file, it's for speed */
+extern struct Render R;
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+
+/* *** utility code to set up an individual raytree for objectinstance, for checking inside/outside *** */
+
+/* Recursive test for intersections, from a point inside the mesh, to outside
+ * Number of intersections (depth) determine if a point is inside or outside the mesh */
+int intersect_outside_volume(RayObject *tree, Isect *isect, float *offset, int limit, int depth)
+{
+	if (limit == 0) return depth;
+	
+	if (RE_rayobject_raycast(tree, isect)) {
+		float hitco[3];
+		
+		hitco[0] = isect->start[0] + isect->labda*isect->vec[0];
+		hitco[1] = isect->start[1] + isect->labda*isect->vec[1];
+		hitco[2] = isect->start[2] + isect->labda*isect->vec[2];
+		VecAddf(isect->start, hitco, offset);
+
+		return intersect_outside_volume(tree, isect, offset, limit-1, depth+1);
+	} else {
+		return depth;
+	}
+}
+
+/* Uses ray tracing to check if a point is inside or outside an ObjectInstanceRen */
+int point_inside_obi(RayObject *tree, ObjectInstanceRen *obi, float *co)
+{
+	Isect isect;
+	float vec[3] = {0.0f,0.0f,1.0f};
+	int final_depth=0, depth=0, limit=20;
+	
+	/* set up the isect */
+	memset(&isect, 0, sizeof(isect));
+	VECCOPY(isect.start, co);
+	VECCOPY(isect.vec, vec);
+	isect.labda = FLT_MAX;
+	
+	/*
+	isect.end[0] = co[0] + vec[0] * maxsize;
+	isect.end[1] = co[1] + vec[1] * maxsize;
+	isect.end[2] = co[2] + vec[2] * maxsize;
+	*/
+	
+	/* and give it a little offset to prevent self-intersections */
+	VecMulf(vec, 1e-5);
+	VecAddf(isect.start, isect.start, vec);
+	
+	isect.mode= RE_RAY_MIRROR;
+	isect.last_hit= NULL;
+	isect.lay= -1;
+	
+	final_depth = intersect_outside_volume(tree, &isect, vec, limit, depth);
+	
+	/* even number of intersections: point is outside
+	 * odd number: point is inside */
+	if (final_depth % 2 == 0) return 0;
+	else return 1;
+}
+
+/*
+static int inside_check_func(Isect *is, int ob, RayObject *face)
+{
+	return 1;
+}
+
+static void vlr_face_coords(RayFace *face, float **v1, float **v2, float **v3, float **v4)
+{
+	VlakRen *vlr= (VlakRen*)face;
+
+	*v1 = (vlr->v1)? vlr->v1->co: NULL;
+	*v2 = (vlr->v2)? vlr->v2->co: NULL;
+	*v3 = (vlr->v3)? vlr->v3->co: NULL;
+	*v4 = (vlr->v4)? vlr->v4->co: NULL;
+}
+
+RayObject *create_raytree_obi(ObjectInstanceRen *obi, float *bbmin, float *bbmax)
+{
+	int v;
+	VlakRen *vlr= NULL;
+	
+	/ * create empty raytree * /
+	RayTree *tree = RE_ray_tree_create(64, obi->obr->totvlak, bbmin, bbmax,
+		vlr_face_coords, inside_check_func, NULL, NULL);
+	
+	/ * fill it with faces * /
+	for(v=0; v<obi->obr->totvlak; v++) {
+		if((v & 255)==0)
+			vlr= obi->obr->vlaknodes[v>>8].vlak;
+		else
+			vlr++;
+	
+		RE_ray_tree_add_face(tree, 0, vlr);
+	}
+	
+	RE_ray_tree_done(tree);
+	
+	return tree;
+}
+*/
+
+/* *** light cache filtering *** */
+
+static float get_avg_surrounds(float *cache, int *res, int xx, int yy, int zz)
+{
+	int x, y, z, x_, y_, z_;
+	int added=0;
+	float tot=0.0f;
+	
+	for (z=-1; z <= 1; z++) {
+		z_ = zz+z;
+		if (z_ >= 0 && z_ <= res[2]-1) {
+		
+			for (y=-1; y <= 1; y++) {
+				y_ = yy+y;
+				if (y_ >= 0 && y_ <= res[1]-1) {
+				
+					for (x=-1; x <= 1; x++) {
+						x_ = xx+x;
+						if (x_ >= 0 && x_ <= res[0]-1) {
+						
+							if (cache[ V_I(x_, y_, z_, res) ] > 0.0f) {
+								tot += cache[ V_I(x_, y_, z_, res) ];
+								added++;
+							}
+							
+						}
+					}
+				}
+			}
+		}
+	}
+	
+	tot /= added;
+	
+	return ((added>0)?tot:0.0f);
+}
+
+/* function to filter the edges of the light cache, where there was no volume originally.
+ * For each voxel which was originally external to the mesh, it finds the average values of
+ * the surrounding internal voxels and sets the original external voxel to that average amount.
+ * Works almost a bit like a 'dilate' filter */
+static void lightcache_filter(VolumePrecache *vp)
+{
+	int x, y, z;
+
+	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.5f)
+					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)
+					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)
+					vp->data_b[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_b, vp->res, x, y, z);
+			}
+		}
+	}
+}
+
+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;
+}
+
+
+/* *** multiple scattering approximation *** */
+
+/* get the total amount of light energy in the light cache. used to normalise after multiple scattering */
+static float total_ss_energy(VolumePrecache *vp)
+{
+	int x, y, z;
+	int *res = vp->res;
+	float energy=0.f;
+	
+	for (z=0; z < res[2]; z++) {
+		for (y=0; y < res[1]; y++) {
+			for (x=0; x < res[0]; x++) {
+				if (vp->data_r[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_r[ V_I(x, y, z, res) ];
+				if (vp->data_g[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_g[ V_I(x, y, z, res) ];
+				if (vp->data_b[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_b[ V_I(x, y, z, res) ];
+			}
+		}
+	}
+	
+	return energy;
+}
+
+static float total_ms_energy(float *sr, float *sg, float *sb, int *res)
+{
+	int x, y, z, i;
+	float energy=0.f;
+	
+	for (z=1;z<=res[2];z++) {
+		for (y=1;y<=res[1];y++) {
+			for (x=1;x<=res[0];x++) {
+			
+				i = ms_I(x,y,z,res);
+				if (sr[i] > 0.f) energy += sr[i];
+				if (sg[i] > 0.f) energy += sg[i];
+				if (sb[i] > 0.f) energy += sb[i];
+			}
+		}
+	}
+	
+	return energy;
+}
+
+static void ms_diffuse(int b, float* x0, float* x, float diff, int *n)
+{
+	int i, j, k, l;
+	const float dt = VOL_MS_TIMESTEP;
+	const float a = dt*diff*n[0]*n[1]*n[2];
+	
+	for (l=0; l<20; l++)
+	{
+		for (k=1; k<=n[2]; k++)
+		{
+			for (j=1; j<=n[1]; j++)
+			{
+				for (i=1; i<=n[0]; i++)
+				{
+					x[ms_I(i,j,k,n)] = (x0[ms_I(i,j,k,n)] + a*(
+						 x[ms_I(i-1,j,k,n)]+x[ms_I(i+1,j,k,n)]+
+						 x[ms_I(i,j-1,k,n)]+x[ms_I(i,j+1,k,n)]+
+						 x[ms_I(i,j,k-1,n)]+x[ms_I(i,j,k+1,n)]))/(1+6*a);
+				}
+			}
+		}
+	}
+}
+
+void multiple_scattering_diffusion(Render *re, VolumePrecache *vp, Material *ma)
+{
+	const float diff = ma->vol.ms_diff * 0.001f; 	/* compensate for scaling for a nicer UI range */
+	const float simframes = ma->vol.ms_steps;
+	const int shade_type = ma->vol.shade_type;
+	float fac = ma->vol.ms_intensity;
+	
+	int x, y, z, m;
+	int *n = vp->res;
+	const int size = (n[0]+2)*(n[1]+2)*(n[2]+2);
+	double time, lasttime= PIL_check_seconds_timer();
+	float total;
+	float c=1.0f;
+	int i;
+	float origf;	/* factor for blending in original light cache */
+	float energy_ss, energy_ms;
+
+	float *sr0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+	float *sr=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+	float *sg0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+	float *sg=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+	float *sb0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+	float *sb=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+
+	total = (float)(n[0]*n[1]*n[2]*simframes);
+	
+	energy_ss = total_ss_energy(vp);
+	
+	/* Scattering as diffusion pass */
+	for (m=0; m<simframes; m++)
+	{
+		/* add sources */
+		for (z=1; z<=n[2]; z++)
+		{
+			for (y=1; y<=n[1]; y++)
+			{
+				for (x=1; x<=n[0]; x++)
+				{
+					i = V_I((x-1), (y-1), (z-1), n);
+					time= PIL_check_seconds_timer();
+					c++;
+										
+					if (vp->data_r[i] > 0.f)
+						sr[ms_I(x,y,z,n)] += vp->data_r[i];
+					if (vp->data_g[i] > 0.f)
+						sg[ms_I(x,y,z,n)] += vp->data_g[i];
+					if (vp->data_b[i] > 0.f)
+						sb[ms_I(x,y,z,n)] += vp->data_b[i];
+					
+					/* Displays progress every second */
+					if(time-lasttime>1.0f) {
+						char str[64];
+						sprintf(str, "Simulating multiple scattering: %d%%", (int)
+								(100.0f * (c / total)));
+						re->i.infostr= str;
+						re->stats_draw(re->sdh, &re->i);
+						re->i.infostr= NULL;
+						lasttime= time;
+					}
+				}
+			}
+		}
+		SWAP(float *, sr, sr0);
+		SWAP(float *, sg, sg0);
+		SWAP(float *, sb, sb0);
+
+		/* main diffusion simulation */
+		ms_diffuse(0, sr0, sr, diff, n);
+		ms_diffuse(0, sg0, sg, diff, n);
+		ms_diffuse(0, sb0, sb, diff, n);
+		
+		if (re->test_break(re->tbh)) break;
+	}
+	
+	/* normalisation factor to conserve energy */
+	energy_ms = total_ms_energy(sr, sg, sb, n);
+	fac *= (energy_ss / energy_ms);
+	
+	/* blend multiple scattering back in the light cache */
+	if (shade_type == MA_VOL_SHADE_SINGLEPLUSMULTIPLE) {
+		/* conserve energy - half single, half multiple */
+		origf = 0.5f;
+		fac *= 0.5f;
+	} else {
+		origf = 0.0f;
+	}
+
+	for (z=1;z<=n[2];z++)
+	{
+		for (y=1;y<=n[1];y++)
+		{
+			for (x=1;x<=n[0];x++)
+			{
+				int index=(x-1)*n[1]*n[2] + (y-1)*n[2] + z-1;
+				vp->data_r[index] = origf * vp->data_r[index] + fac * sr[ms_I(x,y,z,n)];
+				vp->data_g[index] = origf * vp->data_g[index] + fac * sg[ms_I(x,y,z,n)];
+				vp->data_b[index] = origf * vp->data_b[index] + fac * sb[ms_I(x,y,z,n)];
+			}
+		}
+	}
+
+	MEM_freeN(sr0);
+	MEM_freeN(sr);
+	MEM_freeN(sg0);
+	MEM_freeN(sg);
+	MEM_freeN(sb0);
+	MEM_freeN(sb);
+}
+
+
+
+#if 0 // debug stuff
+static void *vol_precache_part_test(void *data)
+{
+	VolPrecachePart *pa = data;
+
+	printf("part number: %d \n", pa->num);
+	printf("done: %d \n", pa->done);
+	printf("x min: %d   x max: %d \n", pa->minx, pa->maxx);
+	printf("y min: %d   y max: %d \n", pa->miny, pa->maxy);
+	printf("z min: %d   z max: %d \n", pa->minz, pa->maxz);
+
+	return NULL;
+}
+#endif
+
+/* Iterate over the 3d voxel grid, and fill the voxels with scattering information
+ *
+ * It's stored in memory as 3 big float grids next to each other, one for each RGB channel.
+ * I'm guessing the memory alignment may work out better this way for the purposes
+ * of doing linear interpolation, but I haven't actually tested this theory! :)
+ */
+static void *vol_precache_part(void *data)
+{
+	VolPrecachePart *pa =  (VolPrecachePart *)data;
+	ObjectInstanceRen *obi = pa->obi;
+	RayObject *tree = pa->tree;
+	ShadeInput *shi = pa->shi;
+	float density, 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));
+		
+		for (y= pa->miny; y < pa->maxy; y++) {
+			co[1] = pa->bbmin[1] + (pa->voxel[1] * (y + 0.5f));
+			
+			for (x=pa->minx; x < pa->maxx; x++) {
+				co[0] = pa->bbmin[0] + (pa->voxel[0] * (x + 0.5f));
+				
+				// don't bother if the point is not inside the volume mesh
+				if (!point_inside_obi(tree, obi, co)) {
+					obi->volume_precache->data_r[ V_I(x, y, z, res) ] = -1.0f;
+					obi->volume_precache->data_g[ V_I(x, y, z, res) ] = -1.0f;
+					obi->volume_precache->data_b[ V_I(x, y, z, res) ] = -1.0f;
+					continue;
+				}
+				
+				VecCopyf(shi->view, co);
+				Normalize(shi->view);
+				density = vol_get_density(shi, co);
+				vol_get_scattering(shi, scatter_col, co, stepsize, density);
+			
+				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];
+				obi->volume_precache->data_b[ V_I(x, y, z, res) ] = scatter_col[2];
+			}
+		}
+	}
+	
+	pa->done = 1;
+	
+	return 0;
+}
+
+
+static void precache_setup_shadeinput(Render *re, ObjectInstanceRen *obi, Material *ma, ShadeInput *shi)
+{
+	memset(shi, 0, sizeof(ShadeInput)); 
+	shi->depth= 1;
+	shi->mask= 1;
+	shi->mat = ma;
+	shi->vlr = NULL;
+	memcpy(&shi->r, &shi->mat->r, 23*sizeof(float));	// note, keep this synced with render_types.h
+	shi->har= shi->mat->har;
+	shi->obi= obi;
+	shi->obr= obi->obr;
+	shi->lay = re->scene->lay;
+}
+
+static void precache_init_parts(Render *re, RayObject *tree, ShadeInput *shi, ObjectInstanceRen *obi, int totthread, int *parts)
+{
+	VolumePrecache *vp = obi->volume_precache;
+	int i=0, x, y, z;
+	float voxel[3];
+	int sizex, sizey, sizez;
+	float *bbmin=obi->obr->boundbox[0], *bbmax=obi->obr->boundbox[1];
+	int *res;
+	int minx, maxx;
+	int miny, maxy;
+	int minz, maxz;
+	
+	if (!vp) return;
+
+	BLI_freelistN(&re->volume_precache_parts);
+	
+	/* currently we just subdivide the box, number of threads per side */
+	parts[0] = parts[1] = parts[2] = totthread;
+	res = vp->res;
+	
+	VecSubf(voxel, bbmax, bbmin);
+	
+	voxel[0] /= res[0];
+	voxel[1] /= res[1];
+	voxel[2] /= res[2];
+
+	for (x=0; x < parts[0]; x++) {
+		sizex = ceil(res[0] / (float)parts[0]);
+		minx = x * sizex;
+		maxx = minx + sizex;
+		maxx = (maxx>res[0])?res[0]:maxx;
+		
+		for (y=0; y < parts[1]; y++) {
+			sizey = ceil(res[1] / (float)parts[1]);
+			miny = y * sizey;
+			maxy = miny + sizey;
+			maxy = (maxy>res[1])?res[1]:maxy;
+			
+			for (z=0; z < parts[2]; z++) {
+				VolPrecachePart *pa= MEM_callocN(sizeof(VolPrecachePart), "new precache part");
+				
+				sizez = ceil(res[2] / (float)parts[2]);
+				minz = z * sizez;
+				maxz = minz + sizez;
+				maxz = (maxz>res[2])?res[2]:maxz;
+						
+				pa->done = 0;
+				pa->working = 0;
+				
+				pa->num = i;
+				pa->tree = tree;
+				pa->shi = shi;
+				pa->obi = obi;
+				VECCOPY(pa->bbmin, bbmin);
+				VECCOPY(pa->voxel, voxel);
+				VECCOPY(pa->res, res);
+				
+				pa->minx = minx; pa->maxx = maxx;
+				pa->miny = miny; pa->maxy = maxy;
+				pa->minz = minz; pa->maxz = maxz;
+				
+				
+				BLI_addtail(&re->volume_precache_parts, pa);
+				
+				i++;
+			}
+		}
+	}
+}
+
+static VolPrecachePart *precache_get_new_part(Render *re)
+{
+	VolPrecachePart *pa, *nextpa=NULL;
+	
+	for (pa = re->volume_precache_parts.first; pa; pa=pa->next)
+	{
+		if (pa->done==0 && pa->working==0) {
+			nextpa = pa;
+			break;
+		}
+	}
+
+	return nextpa;
+}
+
+static int precache_resolution(VolumePrecache *vp, float *bbmin, float *bbmax, int res)
+{
+	float dim[3], div;
+	
+	VecSubf(dim, bbmax, bbmin);
+	
+	div = MAX3(dim[0], dim[1], dim[2]);
+	dim[0] /= div;
+	dim[1] /= div;
+	dim[2] /= div;
+	
+	vp->res[0] = dim[0] * (float)res;
+	vp->res[1] = dim[1] * (float)res;
+	vp->res[2] = dim[2] * (float)res;
+	
+	if ((vp->res[0] < 1) || (vp->res[1] < 1) || (vp->res[2] < 1))
+		return 0;
+	
+	return 1;
+}
+
+/* Precache a volume into a 3D voxel grid.
+ * The voxel grid is stored in the ObjectInstanceRen, 
+ * in camera space, aligned with the ObjectRen's bounding box.
+ * Resolution is defined by the user.
+ */
+void vol_precache_objectinstance_threads(Render *re, ObjectInstanceRen *obi, Material *ma)
+{
+	VolumePrecache *vp;
+	VolPrecachePart *nextpa, *pa;
+	RayObject *tree;
+	ShadeInput shi;
+	ListBase threads;
+	float *bbmin=obi->obr->boundbox[0], *bbmax=obi->obr->boundbox[1];
+	int parts[3], totparts;
+	
+	int caching=1, counter=0;
+	int totthread = re->r.threads;
+	
+	double time, lasttime= PIL_check_seconds_timer();
+	
+	R = *re;
+
+	/* create a raytree with just the faces of the instanced ObjectRen, 
+	 * used for checking if the cached point is inside or outside. */
+	//tree = create_raytree_obi(obi, bbmin, bbmax);
+	tree = makeraytree_object(&R, obi);
+	if (!tree) return;
+	INIT_MINMAX(bbmin, bbmax);
+	RE_rayobject_merge_bb( tree, bbmin, bbmax);
+
+	vp = MEM_callocN(sizeof(VolumePrecache), "volume light cache");
+	
+	if (!precache_resolution(vp, bbmin, bbmax, ma->vol.precache_resolution)) {
+		MEM_freeN(vp);
+		vp = NULL;
+		return;
+	}
+
+	vp->data_r = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data red channel");
+	vp->data_g = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data green channel");
+	vp->data_b = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data blue channel");
+	obi->volume_precache = vp;
+
+	/* Need a shadeinput to calculate scattering */
+	precache_setup_shadeinput(re, obi, ma, &shi);
+	
+	precache_init_parts(re, tree, &shi, obi, totthread, parts);
+	totparts = parts[0] * parts[1] * parts[2];
+	
+	BLI_init_threads(&threads, vol_precache_part, totthread);
+	
+	while(caching) {
+
+		if(BLI_available_threads(&threads) && !(re->test_break(re->tbh))) {
+			nextpa = precache_get_new_part(re);
+			if (nextpa) {
+				nextpa->working = 1;
+				BLI_insert_thread(&threads, nextpa);
+			}
+		}
+		else PIL_sleep_ms(50);
+
+		caching=0;
+		counter=0;
+		for(pa= re->volume_precache_parts.first; pa; pa= pa->next) {
+			
+			if(pa->done) {
+				counter++;
+				BLI_remove_thread(&threads, pa);
+			} else
+				caching = 1;
+		}
+		
+		if (re->test_break(re->tbh) && BLI_available_threads(&threads)==totthread)
+			caching=0;
+		
+		time= PIL_check_seconds_timer();
+		if(time-lasttime>1.0f) {
+			char str[64];
+			sprintf(str, "Precaching volume: %d%%", (int)(100.0f * ((float)counter / (float)totparts)));
+			re->i.infostr= str;
+			re->stats_draw(re->sdh, &re->i);
+			re->i.infostr= NULL;
+			lasttime= time;
+		}
+	}
+	
+	BLI_end_threads(&threads);
+	BLI_freelistN(&re->volume_precache_parts);
+	
+	if(tree) {
+		//TODO: makeraytree_object creates a tree and saves it on OBI, if we free this tree we should also clear other pointers to it
+		//RE_rayobject_free(tree);
+		//tree= NULL;
+	}
+	
+	lightcache_filter(obi->volume_precache);
+	
+	if (ELEM(ma->vol.shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SINGLEPLUSMULTIPLE))
+	{
+		multiple_scattering_diffusion(re, vp, ma);
+	}
+}
+
+static int using_lightcache(Material *ma)
+{
+	return (((ma->vol.shadeflag & MA_VOL_PRECACHESHADING) && (ma->vol.shade_type == MA_VOL_SHADE_SINGLE))
+		|| (ELEM(ma->vol.shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SINGLEPLUSMULTIPLE)));
+}
+
+/* loop through all objects (and their associated materials)
+ * marked for pre-caching in convertblender.c, and pre-cache them */
+void volume_precache(Render *re)
+{
+	ObjectInstanceRen *obi;
+	VolumeOb *vo;
+
+	for(vo= re->volumes.first; vo; vo= vo->next) {
+		if (using_lightcache(vo->ma)) {
+			for(obi= re->instancetable.first; obi; obi= obi->next) {
+				if (obi->obr == vo->obr) {
+					vol_precache_objectinstance_threads(re, obi, vo->ma);
+				}
+			}
+		}
+	}
+	
+	re->i.infostr= NULL;
+	re->stats_draw(re->sdh, &re->i);
+}
+
+void free_volume_precache(Render *re)
+{
+	ObjectInstanceRen *obi;
+	
+	for(obi= re->instancetable.first; obi; obi= obi->next) {
+		if (obi->volume_precache != NULL) {
+			MEM_freeN(obi->volume_precache->data_r);
+			MEM_freeN(obi->volume_precache->data_g);
+			MEM_freeN(obi->volume_precache->data_b);
+			MEM_freeN(obi->volume_precache);
+			obi->volume_precache = NULL;
+		}
+	}
+	
+	BLI_freelistN(&re->volumes);
+}
+
+int point_inside_volume_objectinstance(ObjectInstanceRen *obi, float *co)
+{
+	RayObject *tree;
+	int inside=0;
+	
+	tree = makeraytree_object(&R, obi); //create_raytree_obi(obi, obi->obr->boundbox[0], obi->obr->boundbox[1]);
+	if (!tree) return 0;
+	
+	inside = point_inside_obi(tree, obi, co);
+	
+	//TODO: makeraytree_object creates a tree and saves it on OBI, if we free this tree we should also clear other pointers to it
+	//RE_rayobject_free(tree);
+	//tree= NULL;
+	
+	return inside;
+}
+
diff --git a/source/blender/render/intern/source/volumetric.c b/source/blender/render/intern/source/volumetric.c
new file mode 100644
index 00000000000..2808d266d92
--- /dev/null
+++ b/source/blender/render/intern/source/volumetric.c
@@ -0,0 +1,749 @@
+/**
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Matt Ebb, Raul Fernandez Hernandez (Farsthary)
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_blenlib.h"
+#include "BLI_arithb.h"
+#include "BLI_rand.h"
+#include "BLI_voxel.h"
+
+#include "RE_shader_ext.h"
+#include "RE_raytrace.h"
+
+#include "DNA_material_types.h"
+#include "DNA_group_types.h"
+#include "DNA_lamp_types.h"
+#include "DNA_meta_types.h"
+
+#include "BKE_global.h"
+
+#include "render_types.h"
+#include "pixelshading.h"
+#include "shading.h"
+#include "texture.h"
+#include "volumetric.h"
+#include "volume_precache.h"
+
+#if defined( _MSC_VER ) && !defined( __cplusplus )
+# define inline __inline
+#endif // defined( _MSC_VER ) && !defined( __cplusplus )
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
+/* only to be used here in this file, it's for speed */
+extern struct Render R;
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+
+/* luminance rec. 709 */
+inline float luminance(float* col)
+{
+	return (0.212671f*col[0] + 0.71516f*col[1] + 0.072169f*col[2]);
+}
+
+/* tracing */
+
+static int vol_get_bounds(ShadeInput *shi, float *co, float *vec, float *hitco, Isect *isect, int intersect_type)
+{
+	/* XXX TODO - get raytrace max distance from object instance's bounding box */
+	/* need to account for scaling only, but keep coords in camera space...
+	 * below code is WIP and doesn't work!
+	VecSubf(bb_dim, shi->obi->obr->boundbox[1], shi->obi->obr->boundbox[2]);
+	Mat3MulVecfl(shi->obi->nmat, bb_dim);
+	maxsize = VecLength(bb_dim);
+	*/
+	
+	VECCOPY(isect->start, co);
+	VECCOPY(isect->vec, vec );
+	isect->labda = FLT_MAX;
+	/*
+	isect->end[0] = co[0] + vec[0] * maxsize;
+	isect->end[1] = co[1] + vec[1] * maxsize;
+	isect->end[2] = co[2] + vec[2] * maxsize;
+	*/
+	
+	isect->mode= RE_RAY_MIRROR;
+	isect->orig.ob = (void*)shi->obi;
+	isect->last_hit = NULL;
+	isect->lay= -1;
+	
+	if (intersect_type == VOL_BOUNDS_DEPTH) isect->orig.face = (void*)shi->vlr;
+	else if (intersect_type == VOL_BOUNDS_SS) isect->orig.face= NULL;
+	
+	if(RE_rayobject_raycast(R.raytree, isect))
+	{
+		hitco[0] = isect->start[0] + isect->labda*isect->vec[0];
+		hitco[1] = isect->start[1] + isect->labda*isect->vec[1];
+		hitco[2] = isect->start[2] + isect->labda*isect->vec[2];
+		return 1;
+	} else {
+		return 0;
+	}
+}
+
+static void shade_intersection(ShadeInput *shi, float *col, Isect *is)
+{
+	ShadeInput shi_new;
+	ShadeResult shr_new;
+	
+	memset(&shi_new, 0, sizeof(ShadeInput)); 
+	
+	shi_new.mask= shi->mask;
+	shi_new.osatex= shi->osatex;
+	shi_new.thread= shi->thread;
+	shi_new.depth = shi->depth + 1;
+	shi_new.volume_depth= shi->volume_depth + 1;
+	shi_new.xs= shi->xs;
+	shi_new.ys= shi->ys;
+	shi_new.lay= shi->lay;
+	shi_new.passflag= SCE_PASS_COMBINED; /* result of tracing needs no pass info */
+	shi_new.combinedflag= 0xFFFFFF;		 /* ray trace does all options */
+	shi_new.light_override= shi->light_override;
+	shi_new.mat_override= shi->mat_override;
+	
+	VECCOPY(shi_new.camera_co, is->start);
+	
+	memset(&shr_new, 0, sizeof(ShadeResult));
+	
+	/* hardcoded limit of 100 for now - prevents problems in weird geometry */
+	if (shi->volume_depth < 100) {
+		shade_ray(is, &shi_new, &shr_new);
+	}
+	
+	VecCopyf(col, shr_new.combined);
+	col[3] = shr_new.alpha;
+}
+
+static void vol_trace_behind(ShadeInput *shi, VlakRen *vlr, float *co, float *col)
+{
+	Isect isect;
+	
+	VECCOPY(isect.start, co);
+	VECCOPY(isect.vec, shi->view);
+	isect.labda = FLT_MAX;
+	
+	isect.mode= RE_RAY_MIRROR;
+	isect.orig.ob = (void*) shi->obi;
+	isect.orig.face = (void*)vlr;
+	isect.last_hit = NULL;
+	isect.lay= -1;
+	
+	/* check to see if there's anything behind the volume, otherwise shade the sky */
+	if(RE_rayobject_raycast(R.raytree, &isect)) {
+		shade_intersection(shi, col, &isect);
+	} else {
+		shadeSkyView(col, co, shi->view, NULL, shi->thread);
+		shadeSunView(col, shi->view);
+	}
+}
+
+/* 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)
+{
+	VolumePrecache *vp = shi->obi->volume_precache;
+	float bbmin[3], bbmax[3], dim[3];
+	float sample_co[3];
+	
+	if (!vp) return;
+	
+	/* convert input coords to 0.0, 1.0 */
+	VECCOPY(bbmin, shi->obi->obr->boundbox[0]);
+	VECCOPY(bbmax, shi->obi->obr->boundbox[1]);
+	VecSubf(dim, bbmax, bbmin);
+
+	sample_co[0] = ((co[0] - bbmin[0]) / dim[0]);
+	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);
+}
+
+/* Meta object density, brute force for now 
+ * (might be good enough anyway, don't need huge number of metaobs to model volumetric objects */
+static float metadensity(Object* ob, float* co)
+{
+	float mat[4][4], imat[4][4], dens = 0.f;
+	MetaBall* mb = (MetaBall*)ob->data;
+	MetaElem* ml;
+	
+	/* transform co to meta-element */
+	float tco[3] = {co[0], co[1], co[2]};
+	Mat4MulMat4(mat, ob->obmat, R.viewmat);
+	Mat4Invert(imat, mat);
+	Mat4MulVecfl(imat, tco);
+	
+	for (ml = mb->elems.first; ml; ml=ml->next) {
+		float bmat[3][3], dist2;
+		
+		/* element rotation transform */
+		float tp[3] = {ml->x - tco[0], ml->y - tco[1], ml->z - tco[2]};
+		QuatToMat3(ml->quat, bmat);
+		Mat3Transp(bmat);	// rot.only, so inverse == transpose
+		Mat3MulVecfl(bmat, tp);
+		
+		/* MB_BALL default */
+		switch (ml->type) {
+			case MB_ELIPSOID:
+				tp[0] /= ml->expx, tp[1] /= ml->expy, tp[2] /= ml->expz;
+				break;
+			case MB_CUBE:
+				tp[2] = (tp[2] > ml->expz) ? (tp[2] - ml->expz) : ((tp[2] < -ml->expz) ? (tp[2] + ml->expz) : 0.f);
+				// no break, xy as plane
+			case MB_PLANE:
+				tp[1] = (tp[1] > ml->expy) ? (tp[1] - ml->expy) : ((tp[1] < -ml->expy) ? (tp[1] + ml->expy) : 0.f);
+				// no break, x as tube
+			case MB_TUBE:
+				tp[0] = (tp[0] > ml->expx) ? (tp[0] - ml->expx) : ((tp[0] < -ml->expx) ? (tp[0] + ml->expx) : 0.f);
+		}
+		
+		/* ml->rad2 is not set */
+		dist2 = 1.f - ((tp[0]*tp[0] + tp[1]*tp[1] + tp[2]*tp[2]) / (ml->rad*ml->rad));
+		if (dist2 > 0.f)
+			dens += (ml->flag & MB_NEGATIVE) ? -ml->s*dist2*dist2*dist2 : ml->s*dist2*dist2*dist2;
+	}
+	
+	dens -= mb->thresh;
+	return (dens < 0.f) ? 0.f : dens;
+}
+
+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 meta-object, modulate by metadensity without increasing it
+	if (shi->obi->obr->ob->type == OB_MBALL) {
+		const float md = metadensity(shi->obi->obr->ob, co);
+		if (md < 1.f) density *= md;
+ 	}
+	
+	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)
+{
+	float scatter = shi->mat->vol.scattering;
+	float col[3] = {0.0, 0.0, 0.0};
+	
+	do_volume_tex(shi, co, MAP_SCATTERING, col, &scatter);
+	
+	return 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)
+{
+	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);
+	
+	emission_col[0] = emission_col[0] * emission * density;
+	emission_col[1] = emission_col[1] * emission * density;
+	emission_col[2] = emission_col[2] * emission * density;
+}
+
+void vol_get_absorption(ShadeInput *shi, float *absorb_col, float *co)
+{
+	float absorption = shi->mat->vol.absorption;
+	VECCOPY(absorb_col, shi->mat->vol.absorption_col);
+	
+	do_volume_tex(shi, co, MAP_ABSORPTION+MAP_ABSORPTION_COL, absorb_col, &absorption);
+	
+	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;
+}
+
+/* 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)
+{
+	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
+	 * until Blender's shading system supports this better. --matt
+	 */
+	
+	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);
+		case MA_VOL_PH_MIEMURKY:
+			return scale * (0.5f + 16.5f * powf(0.5 * (1.f + costheta), 32.f)) / (4.f*M_PI);
+		case MA_VOL_PH_RAYLEIGH:
+			return scale * 3.f/(16.f*M_PI) * (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));
+		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)));
+		}
+		case MA_VOL_PH_ISOTROPIC:
+		default:
+			return scale * (1.f / (4.f * M_PI));
+	}
+}
+
+/* Compute transmittance = e^(-attenuation) */
+void vol_get_transmittance_seg(ShadeInput *shi, float *tr, float stepsize, float *co, float density)
+{
+	/* 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);
+	
+	/* homogenous volume within the sampled distance */
+	tau[0] += scatter_dens * absorb[0];
+	tau[1] += scatter_dens * absorb[1];
+	tau[2] += scatter_dens * absorb[2];
+	
+	tr[0] *= exp(-tau[0]);
+	tr[1] *= exp(-tau[1]);
+	tr[2] *= exp(-tau[2]);
+}
+
+/* Compute transmittance = e^(-attenuation) */
+static void vol_get_transmittance(ShadeInput *shi, float *tr, float *co, float *endco)
+{
+	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));
+	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);
+
+	for (; t0 < t1; pt0 = t0, t0 += shi->mat->vol.shade_stepsize) {
+		float absorb[3];
+		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);
+		
+		tau[0] += scatter_dens * absorb[0];
+		tau[1] += scatter_dens * absorb[1];
+		tau[2] += scatter_dens * absorb[2];
+		
+		//if (luminance(tau) >= ambtau) break;
+		VecAddf(p, p, step_vec);
+	}
+	
+	/* return transmittance */
+	tr[0] = expf(-tau[0]);
+	tr[1] = expf(-tau[1]);
+	tr[2] = expf(-tau[2]);
+}
+
+void vol_shade_one_lamp(struct ShadeInput *shi, float *co, LampRen *lar, float *lacol)
+{
+	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 shade_stepsize = vol_get_stepsize(shi, STEPSIZE_SHADE);
+	
+	if (lar->mode & LA_LAYER) if((lar->lay & shi->obi->lay)==0) return;
+	if ((lar->lay & shi->lay)==0) return;
+	if (lar->energy == 0.0) return;
+	
+	if ((visifac= lamp_get_visibility(lar, co, lv, &lampdist)) == 0.f) return;
+	
+	VecCopyf(lacol, &lar->r);
+	
+	if(lar->mode & LA_TEXTURE) {
+		shi->osatex= 0;
+		do_lamp_tex(lar, lv, shi, lacol, LA_TEXTURE);
+	}
+
+	VecMulf(lacol, visifac);
+
+	if (ELEM(lar->type, LA_SUN, LA_HEMI))
+		VECCOPY(lv, lar->vec);
+	VecMulf(lv, -1.0f);
+	
+	if (shi->mat->vol.shade_type != MA_VOL_SHADE_NONE) {
+		Isect is;
+		
+		/* find minimum of volume bounds, or lamp coord */
+		if (vol_get_bounds(shi, co, lv, hitco, &is, VOL_BOUNDS_SS)) {
+			float dist = VecLenf(co, hitco);
+			VlakRen *vlr = (VlakRen *)is.hit.face;
+			
+			/* simple internal shadowing */
+			if (vlr->mat->material_type == MA_TYPE_SURFACE) {
+				lacol[0] = lacol[1] = lacol[2] = 0.0f;
+				return;
+			}
+
+			if (ELEM(lar->type, LA_SUN, LA_HEMI))
+				/* infinite lights, can never be inside volume */
+				atten_co = hitco;
+			else if ( lampdist < dist ) {
+				atten_co = lar->co;
+			} else
+				atten_co = hitco;
+			
+			vol_get_transmittance(shi, tr, co, atten_co);
+			
+			VecMulVecf(lacol, lacol, tr);
+		}
+		else {
+			/* Point is on the outside edge of the volume,
+			 * therefore no attenuation, full transmission.
+			 * Radiance from lamp remains unchanged */
+		}
+	}
+	
+	p = vol_get_phasefunc(shi, shi->mat->vol.phasefunc_type, shi->mat->vol.phasefunc_g, shi->view, lv);
+	VecMulf(lacol, p);
+	
+	scatter_fac = vol_get_scattering_fac(shi, co);
+	VecMulf(lacol, scatter_fac);
+}
+
+/* single scattering only for now */
+void vol_get_scattering(ShadeInput *shi, float *scatter_col, float *co, float stepsize, float density)
+{
+	ListBase *lights;
+	GroupObject *go;
+	LampRen *lar;
+	
+	scatter_col[0] = scatter_col[1] = scatter_col[2] = 0.f;
+	
+	lights= get_lights(shi);
+	for(go=lights->first; go; go= go->next)
+	{
+		float lacol[3] = {0.f, 0.f, 0.f};
+		lar= go->lampren;
+		
+		if (lar) {
+			vol_shade_one_lamp(shi, co, lar, lacol);
+			VecAddf(scatter_col, scatter_col, lacol);
+		}
+	}
+}
+
+	
+/*
+The main volumetric integrator, using an emission/absorption/scattering model.
+
+Incoming radiance = 
+
+outgoing radiance from behind surface * beam transmittance/attenuation
++ added radiance from all points along the ray due to participating media
+	--> radiance for each segment = 
+		(radiance added by scattering + radiance added by emission) * beam transmittance/attenuation
+*/
+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);
+		
+		/* there's only any use in shading here if there's actually some density to shade! */
+		if (density > 0.01f) {
+		
+			/* 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);
+			
+			VecMulf(scatter_col, density);
+			VecAddf(d_radiance, emit_col, scatter_col);
+			
+			/*   Lv += Tr * (Lve() + Ld) */
+			VecMulVecf(d_radiance, tr, d_radiance);
+			VecMulf(d_radiance, stepsize);
+			
+			VecAddf(radiance, radiance, d_radiance);	
+		}
+
+		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;	
+	}
+	
+	/* multiply original color (behind volume) with beam 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]);
+}
+
+/* the main entry point for volume shading */
+static void volume_trace(struct ShadeInput *shi, struct ShadeResult *shr, int inside_volume)
+{
+	float hitco[3], col[4] = {0.f,0.f,0.f,0.f};
+	float *startco, *endco;
+	int trace_behind = 1;
+	const int ztransp= ((shi->depth==0) && (shi->mat->mode & MA_TRANSP) && (shi->mat->mode & MA_ZTRANSP));
+	Isect is;
+
+	/* check for shading an internal face a volume object directly */
+	if (inside_volume == VOL_SHADE_INSIDE)
+		trace_behind = 0;
+	else if (inside_volume == VOL_SHADE_OUTSIDE) {
+		if (shi->flippednor)
+			inside_volume = VOL_SHADE_INSIDE;
+	}
+	
+	if (ztransp && inside_volume == VOL_SHADE_INSIDE) {
+		MatInside *mi;
+		int render_this=0;
+		
+		/* 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
+		 * intersection on the other side, as part of the shading of
+		 * the front face.
+		 
+		 * Because ztransp renders both front and back faces independently
+		 * this will double up, so here we prevent rendering the backface as well, 
+		 * which would otherwise render the volume in between the camera and the backface
+		 * --matt */
+		
+		for (mi=R.render_volumes_inside.first; mi; mi=mi->next) {
+			/* weak... */
+			if (mi->ma == shi->mat) render_this=1;
+		}
+		if (!render_this) return;
+	}
+	
+
+	if (inside_volume == VOL_SHADE_INSIDE)
+	{
+		startco = shi->camera_co;
+		endco = shi->co;
+		
+		if (trace_behind) {
+			if (!ztransp)
+				/* trace behind the volume object */
+				vol_trace_behind(shi, shi->vlr, endco, col);
+		} else {
+			/* we're tracing through the volume between the camera 
+			 * and a solid surface, so use that pre-shaded radiance */
+			QUATCOPY(col, shr->combined);
+		}
+		
+		/* shade volume from 'camera' to 1st hit point */
+		volumeintegrate(shi, col, startco, endco);
+	}
+	/* trace to find a backface, the other side bounds of the volume */
+	/* (ray intersect ignores front faces here) */
+	else if (vol_get_bounds(shi, shi->co, shi->view, hitco, &is, VOL_BOUNDS_DEPTH))
+	{
+		VlakRen *vlr = (VlakRen *)is.hit.face;
+		
+		startco = shi->co;
+		endco = hitco;
+		
+		if (!ztransp) {
+			/* if it's another face in the same material */
+			if (vlr->mat == shi->mat) {
+				/* trace behind the 2nd (raytrace) hit point */
+				vol_trace_behind(shi, (VlakRen *)is.hit.face, endco, col);
+			} else {
+				shade_intersection(shi, col, &is);
+			}
+		}
+		
+		/* shade volume from 1st hit point to 2nd hit point */
+		volumeintegrate(shi, col, startco, endco);
+	}
+	
+	if (ztransp)
+		col[3] = col[3]>1.f?1.f:col[3];
+	else
+		col[3] = 1.f;
+	
+	VecCopyf(shr->combined, col);
+	shr->alpha = col[3];
+	
+	VECCOPY(shr->diff, shr->combined);
+}
+
+/* Traces a shadow through the object, 
+ * pretty much gets the transmission over a ray path */
+void shade_volume_shadow(struct ShadeInput *shi, struct ShadeResult *shr, struct Isect *last_is)
+{
+	float hitco[3];
+	float tr[3] = {1.0,1.0,1.0};
+	Isect is;
+	float shade_stepsize = vol_get_stepsize(shi, STEPSIZE_SHADE);
+	float *startco, *endco;
+	float density=0.f;
+
+	memset(shr, 0, sizeof(ShadeResult));
+	
+	/* if 1st hit normal is facing away from the camera, 
+	 * then we're inside the volume already. */
+	if (shi->flippednor) {
+		startco = last_is->start;
+		endco = shi->co;
+	}
+	/* trace to find a backface, the other side bounds of the volume */
+	/* (ray intersect ignores front faces here) */
+	else if (vol_get_bounds(shi, shi->co, shi->view, hitco, &is, VOL_BOUNDS_DEPTH)) {
+		startco = shi->co;
+		endco = hitco;
+	}
+	else {
+		shr->combined[0] = shr->combined[1] = shr->combined[2] = 0.f;
+		shr->alpha = shr->combined[3] = 1.f;
+		return;
+	}
+	
+	density = vol_get_density(shi, startco);
+	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];
+}
+
+
+/* delivers a fully filled in ShadeResult, for all passes */
+void shade_volume_outside(ShadeInput *shi, ShadeResult *shr)
+{
+	memset(shr, 0, sizeof(ShadeResult));
+	volume_trace(shi, shr, VOL_SHADE_OUTSIDE);
+}
+
+
+void shade_volume_inside(ShadeInput *shi, ShadeResult *shr)
+{
+	MatInside *m;
+	Material *mat_backup;
+	ObjectInstanceRen *obi_backup;
+	float prev_alpha = shr->alpha;
+	
+	//if (BLI_countlist(&R.render_volumes_inside) == 0) return;
+	
+	/* XXX: extend to multiple volumes perhaps later */
+	mat_backup = shi->mat;
+	obi_backup = shi->obi;
+	
+	m = R.render_volumes_inside.first;
+	shi->mat = m->ma;
+	shi->obi = m->obi;
+	shi->obr = m->obi->obr;
+	
+	volume_trace(shi, shr, VOL_SHADE_INSIDE);
+	shr->alpha += prev_alpha;
+	CLAMP(shr->alpha, 0.f, 1.f);
+	
+	shi->mat = mat_backup;
+	shi->obi = obi_backup;
+	shi->obr = obi_backup->obr;
+}
\ No newline at end of file
diff --git a/source/blender/render/intern/source/voxeldata.c b/source/blender/render/intern/source/voxeldata.c
new file mode 100644
index 00000000000..17858e55e3d
--- /dev/null
+++ b/source/blender/render/intern/source/voxeldata.c
@@ -0,0 +1,346 @@
+/**
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Raul Fernandez Hernandez (Farsthary), Matt Ebb.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_arithb.h"
+#include "BLI_blenlib.h"
+#include "BLI_voxel.h"
+
+#include "IMB_imbuf.h"
+#include "IMB_imbuf_types.h"
+
+#include "BKE_global.h"
+#include "BKE_image.h"
+#include "BKE_main.h"
+#include "BKE_modifier.h"
+
+#include "smoke_API.h"
+
+#include "DNA_texture_types.h"
+#include "DNA_object_types.h"
+#include "DNA_modifier_types.h"
+#include "DNA_smoke_types.h"
+
+
+#include "render_types.h"
+#include "renderdatabase.h"
+#include "texture.h"
+#include "voxeldata.h"
+
+void load_frame_blendervoxel(FILE *fp, float *F, int size, int frame, int offset)
+{	
+	fseek(fp,frame*size*sizeof(float)+offset,0);
+	fread(F,sizeof(float),size,fp);
+}
+
+void load_frame_raw8(FILE *fp, float *F, int size, int frame)
+{
+	char *tmp;
+	int i;
+	
+	tmp = (char *)MEM_mallocN(sizeof(char)*size, "temporary voxel file reading storage");
+	
+	fseek(fp,(frame-1)*size*sizeof(char),0);
+	fread(tmp, sizeof(char), size, fp);
+	
+	for (i=0; i<size; i++) {
+		F[i] = (float)tmp[i] / 256.f;
+	}
+	MEM_freeN(tmp);
+}
+
+void load_frame_image_sequence(Render *re, VoxelData *vd, Tex *tex)
+{
+	ImBuf *ibuf;
+	Image *ima = tex->ima;
+	ImageUser *iuser = &tex->iuser;
+	int x=0, y=0, z=0;
+	float *rf;
+
+	if (!ima || !iuser) return;
+	
+	ima->source = IMA_SRC_SEQUENCE;
+	iuser->framenr = 1 + iuser->offset;
+
+	/* find the first valid ibuf and use it to initialise the resolution of the data set */
+	/* need to do this in advance so we know how much memory to allocate */
+	ibuf= BKE_image_get_ibuf(ima, iuser);
+	while (!ibuf && (iuser->framenr < iuser->frames)) {
+		iuser->framenr++;
+		ibuf= BKE_image_get_ibuf(ima, iuser);
+	}
+	if (!ibuf) return;
+	if (!ibuf->rect_float) IMB_float_from_rect(ibuf);
+	
+	vd->flag |= TEX_VD_STILL;
+	vd->resol[0] = ibuf->x;
+	vd->resol[1] = ibuf->y;
+	vd->resol[2] = iuser->frames;
+	vd->dataset = MEM_mapallocN(sizeof(float)*(vd->resol[0])*(vd->resol[1])*(vd->resol[2]), "voxel dataset");
+	
+	for (z=0; z < iuser->frames; z++)
+	{	
+		/* get a new ibuf for each frame */
+		if (z > 0) {
+			iuser->framenr++;
+			ibuf= BKE_image_get_ibuf(ima, iuser);
+			if (!ibuf) break;
+			if (!ibuf->rect_float) IMB_float_from_rect(ibuf);
+		}
+		rf = ibuf->rect_float;
+		
+		for (y=0; y < ibuf->y; y++)
+		{
+			for (x=0; x < ibuf->x; x++)
+			{
+				/* currently converted to monchrome */
+				vd->dataset[ V_I(x, y, z, vd->resol) ] = (rf[0] + rf[1] + rf[2])*0.333f;
+				rf +=4;
+			}
+		}
+		
+		BKE_image_free_anim_ibufs(ima, iuser->framenr);
+	}
+}
+
+void write_voxeldata_header(struct VoxelDataHeader *h, FILE *fp)
+{
+	fwrite(h,sizeof(struct VoxelDataHeader),1,fp);
+}
+
+void read_voxeldata_header(FILE *fp, struct VoxelData *vd)
+{
+	VoxelDataHeader *h=(VoxelDataHeader *)MEM_mallocN(sizeof(VoxelDataHeader), "voxel data header");
+	
+	rewind(fp);
+	fread(h,sizeof(VoxelDataHeader),1,fp);
+	
+	vd->resol[0]=h->resolX;
+	vd->resol[1]=h->resolY;
+	vd->resol[2]=h->resolZ;
+
+	MEM_freeN(h);
+}
+
+void init_frame_smoke(Render *re, VoxelData *vd, Tex *tex)
+{
+	Object *ob;
+	ModifierData *md;
+	
+	vd->dataset = NULL;
+	if (vd->object == NULL)	return;	
+	ob= vd->object;
+	
+	/* draw code for smoke */
+	if( (md = (ModifierData *)modifiers_findByType(ob, eModifierType_Smoke)) )
+	{
+		SmokeModifierData *smd = (SmokeModifierData *)md;
+		
+		if(smd->domain && smd->domain->fluid) {
+			//int big = (smd->domain->flags & MOD_SMOKE_HIGHRES);
+			int big=0;
+			
+			if (big) {
+				//smoke_turbulence_get_res(smd->domain->wt, vd->resol);
+				//vd->dataset = smoke_turbulence_get_density(smd->domain->wt);
+			} else {
+				VECCOPY(vd->resol, smd->domain->res);
+				vd->dataset = smoke_get_density(smd->domain->fluid);
+			}
+		}
+	}
+}
+
+void cache_voxeldata(struct Render *re,Tex *tex)
+{	
+	VoxelData *vd = tex->vd;
+	FILE *fp;
+	int size;
+	int curframe;
+	
+	if (!vd) return;
+	
+	/* image sequence gets special treatment */
+	if (vd->file_format == TEX_VD_IMAGE_SEQUENCE) {
+		load_frame_image_sequence(re, vd, tex);
+		return;
+	} else if (vd->file_format == TEX_VD_SMOKE) {
+		init_frame_smoke(re, vd, tex);
+		return;
+	}
+
+	if (!BLI_exists(vd->source_path)) return;
+	fp = fopen(vd->source_path,"rb");
+	if (!fp) return;
+
+	if (vd->file_format == TEX_VD_BLENDERVOXEL)
+		read_voxeldata_header(fp, vd);
+	
+	size = (vd->resol[0])*(vd->resol[1])*(vd->resol[2]);
+	vd->dataset = MEM_mapallocN(sizeof(float)*size, "voxel dataset");
+		
+	if (vd->flag & TEX_VD_STILL) curframe = vd->still_frame;
+	else curframe = re->r.cfra;
+	
+	switch(vd->file_format) {
+		case TEX_VD_BLENDERVOXEL:
+			load_frame_blendervoxel(fp, vd->dataset, size, curframe-1, sizeof(VoxelDataHeader));
+			break;
+		case TEX_VD_RAW_8BIT:
+			load_frame_raw8(fp, vd->dataset, size, curframe);
+			break;
+	}
+	
+	fclose(fp);
+}
+
+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);
+	
+	/* XXX: should be doing only textures used in this render */
+	for (tex= G.main->tex.first; tex; tex= tex->id.next) {
+		if(tex->id.us && tex->type==TEX_VOXELDATA) {
+			cache_voxeldata(re, tex);
+		}
+	}
+	
+	re->i.infostr= NULL;
+	re->stats_draw(re->sdh, &re->i);
+	
+}
+
+static void free_voxeldata_one(Render *re, Tex *tex)
+{
+	VoxelData *vd = tex->vd;
+	
+	if (vd->dataset) {
+		MEM_freeN(vd->dataset);
+		vd->dataset = NULL;
+	}
+}
+
+
+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);
+		}
+	}
+}
+
+int voxeldatatex(struct Tex *tex, float *texvec, struct TexResult *texres)
+{	 
+    int retval = TEX_INT;
+	VoxelData *vd = tex->vd;	
+	float co[3], offset[3] = {0.5, 0.5, 0.5};
+
+	if ((!vd) || (vd->dataset==NULL)) {
+		texres->tin = 0.0f;
+		return 0;
+	}
+	
+	/* scale lookup from 0.0-1.0 (original location) to -1.0, 1.0, consistent with image texture tex coords */
+	/* in implementation this works backwards, bringing sample locations from -1.0, 1.0
+	 * to the range 0.0, 1.0, before looking up in the voxel structure. */
+	VecCopyf(co, texvec);
+	VecMulf(co, 0.5f);
+	VecAddf(co, co, offset);
+
+	/* co is now in the range 0.0, 1.0 */
+	switch (tex->extend) {
+		case TEX_CLIP:
+		{
+			if ((co[0] < 0.f || co[0] > 1.f) || (co[1] < 0.f || co[1] > 1.f) || (co[2] < 0.f || co[2] > 1.f)) {
+				texres->tin = 0.f;
+				return retval;
+			}
+			break;
+		}
+		case TEX_REPEAT:
+		{
+			co[0] = co[0] - floor(co[0]);
+			co[1] = co[1] - floor(co[1]);
+			co[2] = co[2] - floor(co[2]);
+			break;
+		}
+		case TEX_EXTEND:
+		{
+			CLAMP(co[0], 0.f, 1.f);
+			CLAMP(co[1], 0.f, 1.f);
+			CLAMP(co[2], 0.f, 1.f);
+			break;
+		}
+	}
+	
+	switch (vd->interp_type) {
+		case TEX_VD_NEARESTNEIGHBOR:
+			texres->tin = voxel_sample_nearest(vd->dataset, vd->resol, co);
+			break;  
+		case TEX_VD_LINEAR:
+			texres->tin = voxel_sample_trilinear(vd->dataset, vd->resol, co);
+			break;					
+		case TEX_VD_QUADRATIC:
+			texres->tin = voxel_sample_triquadratic(vd->dataset, vd->resol, co);
+			break;
+		case TEX_VD_TRICUBIC_CATROM:
+		case TEX_VD_TRICUBIC_BSPLINE:
+			texres->tin = voxel_sample_tricubic(vd->dataset, vd->resol, co, (vd->interp_type == TEX_VD_TRICUBIC_BSPLINE));
+			break;
+	}
+	
+	texres->tin *= vd->int_multiplier;
+	BRICONT;
+	
+	texres->tr = texres->tin;
+	texres->tg = texres->tin;
+	texres->tb = texres->tin;
+	texres->ta = texres->tin;
+	BRICONTRGB;
+	
+	return retval;	
+}
+
+
diff --git a/source/blender/render/intern/source/zbuf.c b/source/blender/render/intern/source/zbuf.c
index 53a05dd0d67..3b3a8568933 100644
--- a/source/blender/render/intern/source/zbuf.c
+++ b/source/blender/render/intern/source/zbuf.c
@@ -42,7 +42,7 @@
 #include "BLI_jitter.h"
 #include "BLI_threads.h"
 
-#include "MTC_matrixops.h"
+
 #include "MEM_guardedalloc.h"
 
 #include "DNA_lamp_types.h"