Voxel data & volume light cache

* Added support for additional file types in the voxel data texture. I added 
support for 8 bit raw files, but most notably for image sequences.

Image sequences generate the voxel grid by stacking layers of image slices on top 
of each other to generate the voxels in the Z axis - the number of slices in the 
sequence is the resolution of the voxel grid's Z axis. 

i.e. http://mke3.net/blender/devel/rendering/volumetrics/skull_layers.jpg

The image sequence option is particularly useful for loading medical data into 
Blender. 3d medical data such as MRIs or CT scans are often stored as DICOM 
format image sequences. It's not in Blender's scope to support the DICOM format, 
but there are plenty of utilities such as ImageMagick, Photoshop or OsiriX that 
can easily convert DICOM files to formats that Blender supports, such as PNG or JPEG.

Here are some example renderings using these file formats to load medical data:
http://vimeo.com/5242961
http://vimeo.com/5242989
http://vimeo.com/5243228

Currently the 8 bit raw and image sequence formats only support the 'still' 
rendering feature.

* Changed the default texture placement to be centred around 0.5,0.5,0.5, rather 
than within the 0.0,1.0 cube. This is more consistent with image textures, and it 
also means you can easily add a voxel data texture to a default cube without 
having to mess around with mapping.

* Added some more extrapolation modes such as Repeat and Extend rather than just clipping

http://mke3.net/blender/devel/rendering/volumetrics/bradybunch.jpg

* Changed the voxel data storage to use MEM_Mapalloc (memory mapped disk cache) 
rather than storing in ram, to help cut down memory usage.

* Refactored and cleaned up the code a lot. Now the access and interpolation code 
is separated into a separate voxel library inside blenlib. This is now properly 
shared between voxel data texture and light cache (previously there was some 
duplicated code).

* Made volume light cache support non-cubic voxel grids. Now the resolution 
specified in the material properties is used for the longest edge of the volume 
object's bounding box, and the shorter edges are proportional (similar to how 
resolution is calculated for fluid sim domains).

This is *much* more memory efficient for squashed volume regions like clouds 
layer bounding boxes, allowing you to raise the resolution considerably while 
still keeping memory usage acceptable.

6 files changed, 302 insertions, 519 deletions

diff --git a/source/blender/render/intern/include/render_types.h b/source/blender/render/intern/include/render_types.h
index 8dbdde77726..9703ab10821 100644
--- a/source/blender/render/intern/include/render_types.h
+++ b/source/blender/render/intern/include/render_types.h
@@ -290,8 +290,8 @@ typedef struct ObjectInstanceRen {
 	float dupliorco[3], dupliuv[2];
 	float (*duplitexmat)[4];
 	
-	float *volume_precache;
-
+	struct VolumePrecache *volume_precache;
+	
 	float *vectors;
 	int totvector;
 } ObjectInstanceRen;
@@ -426,12 +426,19 @@ typedef struct VolPrecachePart
 	int minx, maxx;
 	int miny, maxy;
 	int minz, maxz;
-	int res;
+	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;
 
 /* ------------------------------------------------------------------------- */
 
diff --git a/source/blender/render/intern/include/voxeldata.h b/source/blender/render/intern/include/voxeldata.h
index 36c7789da80..b291bdc096d 100644
--- a/source/blender/render/intern/include/voxeldata.h
+++ b/source/blender/render/intern/include/voxeldata.h
@@ -29,23 +29,17 @@
 #ifndef VOXELDATA_H
 #define VOXELDATA_H 
 
-/**
- * Load voxel data for all point density textures in the scene
- */
-
 struct Render;
 struct TexResult;
 
 typedef struct VoxelDataHeader
 {
-int resolX, resolY, resolZ;
-int frames;
+	int resolX, resolY, resolZ;
+	int frames;
 } VoxelDataHeader;
 
-__inline int _I(int x, int y, int z, int *n);
 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/shadeinput.c b/source/blender/render/intern/source/shadeinput.c
index 294175460bf..59ab5661c19 100644
--- a/source/blender/render/intern/source/shadeinput.c
+++ b/source/blender/render/intern/source/shadeinput.c
@@ -51,6 +51,7 @@
 #include "shading.h"
 #include "strand.h"
 #include "texture.h"
+#include "volumetric.h"
 #include "zbuf.h"
 
 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
diff --git a/source/blender/render/intern/source/volume_precache.c b/source/blender/render/intern/source/volume_precache.c
index b3cffb45a27..98a1a2ec2c8 100644
--- a/source/blender/render/intern/source/volume_precache.c
+++ b/source/blender/render/intern/source/volume_precache.c
@@ -36,6 +36,7 @@
 #include "BLI_blenlib.h"
 #include "BLI_arithb.h"
 #include "BLI_threads.h"
+#include "BLI_voxel.h"
 
 #include "PIL_time.h"
 
@@ -61,6 +62,7 @@
 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 */
@@ -151,28 +153,28 @@ RayTree *create_raytree_obi(ObjectInstanceRen *obi, float *bbmin, float *bbmax)
 	return tree;
 }
 
-static float get_avg_surrounds(float *cache, int res, int res_2, int res_3, int rgb, int xx, int yy, int zz)
+/* *** 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;
-	int i;
 	
-	for (x=-1; x <= 1; x++) {
-		x_ = xx+x;
-		if (x_ >= 0 && x_ <= res-1) {
+	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) {
+				if (y_ >= 0 && y_ <= res[1]-1) {
 				
-					for (z=-1; z <= 1; z++) {
-						z_ = zz+z;
-						if (z_ >= 0 && z_ <= res-1) {
+					for (x=-1; x <= 1; x++) {
+						x_ = xx+x;
+						if (x_ >= 0 && x_ <= res[0]-1) {
 						
-							i = rgb*res_3 + x_*res_2 + y_*res + z_;
-							if (cache[i] > 0.0f) {
-								tot += cache[i];
+							if (cache[ V_I(x_, y_, z_, res) ] > 0.0f) {
+								tot += cache[ V_I(x_, y_, z_, res) ];
 								added++;
 							}
 							
@@ -192,54 +194,46 @@ static float get_avg_surrounds(float *cache, int res, int res_2, int res_3, int
  * 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(float *cache, int res)
+static void lightcache_filter(VolumePrecache *vp)
 {
-	int x, y, z, rgb;
-	int res_2, res_3;
-	int i;
-	
-	res_2 = res*res;
-	res_3 = res*res*res;
-
-	for (x=0; x < res; x++) {
-		for (y=0; y < res; y++) {
-			for (z=0; z < res; z++) {
-				for (rgb=0; rgb < 3; rgb++) {
-					i = rgb*res_3 + x*res_2 + y*res + z;
+	int x, y, z;
 
-					/* trigger for outside mesh */
-					if (cache[i] < -0.5f) cache[i] = get_avg_surrounds(cache, res, res_2, res_3, rgb, 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 I(int x,int y,int z,int n) //has a pad of 1 voxel surrounding the core for boundary simulation
+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 (x*(n+2)+y)*(n+2)+z;
+	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(float *cache, const int res)
+static float total_ss_energy(VolumePrecache *vp)
 {
-	int x, y, z, rgb;
-	int res_2, res_3;
-	int i;
+	int x, y, z;
+	int *res = vp->res;
 	float energy=0.f;
 	
-	res_2 = res*res;
-	res_3 = res*res*res;
-
-	for (x=0; x < res; x++) {
-		for (y=0; y < res; y++) {
-			for (z=0; z < res; z++) {
-				for (rgb=0; rgb < 3; rgb++) {
-					i = rgb*res_3 + x*res_2 + y*res + z;
-					
-					if (cache[i] > 0.f) energy += cache[i];
-				}
+	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) ];
 			}
 		}
 	}
@@ -247,16 +241,16 @@ static float total_ss_energy(float *cache, const int res)
 	return energy;
 }
 
-static float total_ms_energy(float *sr, float *sg, float *sb, const int res)
+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;z++) {
-		for (y=1;y<=res;y++) {
-			for (x=1;x<=res;x++) {
+	for (z=1;z<=res[2];z++) {
+		for (y=1;y<=res[1];y++) {
+			for (x=1;x<=res[0];x++) {
 			
-				i = I(x,y,z,res);
+				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];
@@ -267,44 +261,44 @@ static float total_ms_energy(float *sr, float *sg, float *sb, const int res)
 	return energy;
 }
 
-static void ms_diffuse(int b, float* x0, float* x, float diff, int n)
+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*n*n;
+	const float a = dt*diff*n[0]*n[1]*n[2];
 	
 	for (l=0; l<20; l++)
 	{
-		for (k=1; k<=n; k++)
+		for (k=1; k<=n[2]; k++)
 		{
-			for (j=1; j<=n; j++)
+			for (j=1; j<=n[1]; j++)
 			{
-				for (i=1; i<=n; i++)
+				for (i=1; i<=n[0]; i++)
 				{
-					x[I(i,j,k,n)] = (x0[I(i,j,k,n)] + a*(
-														 x[I(i-1,j,k,n)]+x[I(i+1,j,k,n)]+
-														 x[I(i,j-1,k,n)]+x[I(i,j+1,k,n)]+
-														 x[I(i,j,k-1,n)]+x[I(i,j,k+1,n)]))/(1+6*a);
+					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, float *cache, int res, Material *ma)
+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 i, j, k, m;
-	int n = res;
-	const int size = (n+2)*(n+2)*(n+2);
+	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 index;
+	int i;
 	float origf;	/* factor for blending in original light cache */
 	float energy_ss, energy_ms;
 
@@ -315,33 +309,31 @@ void multiple_scattering_diffusion(Render *re, float *cache, int res, Material *
 	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*n*n*simframes);
+	total = (float)(n[0]*n[1]*n[2]*simframes);
 	
-	energy_ss = total_ss_energy(cache, res);
+	energy_ss = total_ss_energy(vp);
 	
 	/* Scattering as diffusion pass */
 	for (m=0; m<simframes; m++)
 	{
 		/* add sources */
-		for (k=1; k<=n; k++)
+		for (z=1; z<=n[2]; z++)
 		{
-			for (j=1; j<=n; j++)
+			for (y=1; y<=n[1]; y++)
 			{
-				for (i=1; i<=n; i++)
+				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];
 					
-					index=(i-1)*n*n + (j-1)*n + k-1;
-					
-					if (cache[index] > 0.0f)
-						sr[I(i,j,k,n)] += cache[index];
-					if (cache[1*n*n*n + index] > 0.0f)
-						sg[I(i,j,k,n)] += cache[1*n*n*n + index];
-					if (cache[2*n*n*n + index] > 0.0f)
-						sb[I(i,j,k,n)] += cache[2*n*n*n + index];
-
-
 					/* Displays progress every second */
 					if(time-lasttime>1.0f) {
 						char str[64];
@@ -368,7 +360,7 @@ void multiple_scattering_diffusion(Render *re, float *cache, int res, Material *
 	}
 	
 	/* normalisation factor to conserve energy */
-	energy_ms = total_ms_energy(sr, sg, sb, res);
+	energy_ms = total_ms_energy(sr, sg, sb, n);
 	fac *= (energy_ss / energy_ms);
 	
 	/* blend multiple scattering back in the light cache */
@@ -380,16 +372,16 @@ void multiple_scattering_diffusion(Render *re, float *cache, int res, Material *
 		origf = 0.0f;
 	}
 
-	for (k=1;k<=n;k++)
+	for (z=1;z<=n[2];z++)
 	{
-		for (j=1;j<=n;j++)
+		for (y=1;y<=n[1];y++)
 		{
-			for (i=1;i<=n;i++)
+			for (x=1;x<=n[0];x++)
 			{
-				index=(i-1)*n*n + (j-1)*n + k-1;
-				cache[index]			= origf * cache[index]  + fac * sr[I(i,j,k,n)];
-				cache[1*n*n*n + index]	= origf * cache[1*n*n*n + index] + fac * sg[I(i,j,k,n)];
-				cache[2*n*n*n + index]	= origf * cache[2*n*n*n + index] + fac * sb[I(i,j,k,n)];
+				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)];
 			}
 		}
 	}
@@ -434,23 +426,23 @@ static void *vol_precache_part(void *data)
 	float density, scatter_col[3] = {0.f, 0.f, 0.f};
 	float co[3];
 	int x, y, z;
-	const int res=pa->res, res_2=pa->res*pa->res, res_3=pa->res*pa->res*pa->res;
+	const int res[3]= {pa->res[0], pa->res[1], pa->res[2]};
 	const float stepsize = vol_get_stepsize(shi, STEPSIZE_VIEW);
 
-	for (x= pa->minx; x < pa->maxx; x++) {
-		co[0] = pa->bbmin[0] + (pa->voxel[0] * x);
+	for (z= pa->minz; z < pa->maxz; z++) {
+		co[2] = pa->bbmin[2] + (pa->voxel[2] * z);
 		
 		for (y= pa->miny; y < pa->maxy; y++) {
 			co[1] = pa->bbmin[1] + (pa->voxel[1] * y);
 			
-			for (z=pa->minz; z < pa->maxz; z++) {
-				co[2] = pa->bbmin[2] + (pa->voxel[2] * z);
+			for (x=pa->minx; x < pa->maxx; x++) {
+				co[0] = pa->bbmin[0] + (pa->voxel[0] * x);
 				
 				// don't bother if the point is not inside the volume mesh
 				if (!point_inside_obi(tree, obi, co)) {
-					obi->volume_precache[0*res_3 + x*res_2 + y*res + z] = -1.0f;
-					obi->volume_precache[1*res_3 + x*res_2 + y*res + z] = -1.0f;
-					obi->volume_precache[2*res_3 + x*res_2 + y*res + z] = -1.0f;
+					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;
 				}
 				
@@ -459,9 +451,9 @@ static void *vol_precache_part(void *data)
 				density = vol_get_density(shi, co);
 				vol_get_scattering(shi, scatter_col, co, stepsize, density);
 			
-				obi->volume_precache[0*res_3 + x*res_2 + y*res + z] = scatter_col[0];
-				obi->volume_precache[1*res_3 + x*res_2 + y*res + z] = scatter_col[1];
-				obi->volume_precache[2*res_3 + x*res_2 + y*res + z] = scatter_col[2];
+				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];
 			}
 		}
 	}
@@ -486,44 +478,52 @@ static void precache_setup_shadeinput(Render *re, ObjectInstanceRen *obi, Materi
 	shi->lay = re->scene->lay;
 }
 
-static void precache_init_parts(Render *re, RayTree *tree, ShadeInput *shi, ObjectInstanceRen *obi, float *bbmin, float *bbmax, int res, int totthread, int *parts)
+static void precache_init_parts(Render *re, RayTree *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);
 	if ((voxel[0] < FLT_EPSILON) || (voxel[1] < FLT_EPSILON) || (voxel[2] < FLT_EPSILON))
 		return;
-	VecMulf(voxel, 1.0f/res);
+	voxel[0] /= res[0];
+	voxel[1] /= res[1];
+	voxel[2] /= res[2];
 
 	for (x=0; x < parts[0]; x++) {
-		sizex = ceil(res / (float)parts[0]);
+		sizex = ceil(res[0] / (float)parts[0]);
 		minx = x * sizex;
 		maxx = minx + sizex;
-		maxx = (maxx>res)?res:maxx;
+		maxx = (maxx>res[0])?res[0]:maxx;
 		
 		for (y=0; y < parts[1]; y++) {
-			sizey = ceil(res / (float)parts[1]);
+			sizey = ceil(res[1] / (float)parts[1]);
 			miny = y * sizey;
 			maxy = miny + sizey;
-			maxy = (maxy>res)?res:maxy;
+			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 / (float)parts[2]);
+				sizez = ceil(res[2] / (float)parts[2]);
 				minz = z * sizez;
 				maxz = minz + sizez;
-				maxz = (maxz>res)?res:maxz;
+				maxz = (maxz>res[2])?res[2]:maxz;
 						
 				pa->done = 0;
 				pa->working = 0;
@@ -534,7 +534,7 @@ static void precache_init_parts(Render *re, RayTree *tree, ShadeInput *shi, Obje
 				pa->obi = obi;
 				VECCOPY(pa->bbmin, bbmin);
 				VECCOPY(pa->voxel, voxel);
-				pa->res = res;
+				VECCOPY(pa->res, res);
 				
 				pa->minx = minx; pa->maxx = maxx;
 				pa->miny = miny; pa->maxy = maxy;
@@ -564,19 +564,37 @@ static VolPrecachePart *precache_get_new_part(Render *re)
 	return nextpa;
 }
 
+static void 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;
+}
+
 /* 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, float *bbmin, float *bbmax)
+void vol_precache_objectinstance_threads(Render *re, ObjectInstanceRen *obi, Material *ma)
 {
+	VolumePrecache *vp;
 	VolPrecachePart *nextpa, *pa;
 	RayTree *tree;
 	ShadeInput shi;
 	ListBase threads;
-	const int res = ma->vol_precache_resolution;
+	float *bbmin=obi->obr->boundbox[0], *bbmax=obi->obr->boundbox[1];
 	int parts[3], totparts;
+	int res[3];
 	
 	int caching=1, counter=0;
 	int totthread = re->r.threads;
@@ -589,13 +607,19 @@ void vol_precache_objectinstance_threads(Render *re, ObjectInstanceRen *obi, Mat
 	 * used for checking if the cached point is inside or outside. */
 	tree = create_raytree_obi(obi, bbmin, bbmax);
 	if (!tree) return;
-	
-	obi->volume_precache = MEM_callocN(sizeof(float)*res*res*res*3, "volume light cache");
+
+	vp = MEM_callocN(sizeof(VolumePrecache), "volume light cache");
+	precache_resolution(vp, bbmin, bbmax, ma->vol_precache_resolution);
+
+	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, bbmin, bbmax, res, totthread, parts);
+	precache_init_parts(re, tree, &shi, obi, totthread, parts);
 	totparts = parts[0] * parts[1] * parts[2];
 	
 	BLI_init_threads(&threads, vol_precache_part, totthread);
@@ -644,11 +668,11 @@ void vol_precache_objectinstance_threads(Render *re, ObjectInstanceRen *obi, Mat
 		tree= NULL;
 	}
 	
-	lightcache_filter(obi->volume_precache, res);
+	lightcache_filter(obi->volume_precache);
 	
 	if (ELEM(ma->vol_shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SINGLEPLUSMULTIPLE))
 	{
-		multiple_scattering_diffusion(re, obi->volume_precache, res, ma);
+		multiple_scattering_diffusion(re, vp, ma);
 	}
 }
 
@@ -672,7 +696,7 @@ void volume_precache(Render *re)
 		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, obi->obr->boundbox[0], obi->obr->boundbox[1]);
+					vol_precache_objectinstance_threads(re, obi, vo->ma);
 				}
 			}
 		}
@@ -687,8 +711,12 @@ void free_volume_precache(Render *re)
 	ObjectInstanceRen *obi;
 	
 	for(obi= re->instancetable.first; obi; obi= obi->next) {
-		if (obi->volume_precache)
+		if (obi->volume_precache != NULL) {
 			MEM_freeN(obi->volume_precache);
+			MEM_freeN(obi->volume_precache->data_r);
+			MEM_freeN(obi->volume_precache->data_g);
+			MEM_freeN(obi->volume_precache->data_b);
+		}
 	}
 	
 	BLI_freelistN(&re->volumes);
diff --git a/source/blender/render/intern/source/volumetric.c b/source/blender/render/intern/source/volumetric.c
index 01c4ac8fa74..65f2916af59 100644
--- a/source/blender/render/intern/source/volumetric.c
+++ b/source/blender/render/intern/source/volumetric.c
@@ -36,6 +36,7 @@
 #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"
@@ -139,79 +140,27 @@ static float vol_get_depth_cutoff(struct ShadeInput *shi)
 	return shi->mat->vol_depth_cutoff;
 }
 
-/* SHADING */
-
-static float D(ShadeInput *shi, int rgb, int x, int y, int z)
-{
-	const int res = shi->mat->vol_precache_resolution;
-	CLAMP(x, 0, res-1);
-	CLAMP(y, 0, res-1);
-	CLAMP(z, 0, res-1);
-	return shi->obi->volume_precache[rgb*res*res*res + x*res*res + y*res + z];
-}
-
-static inline float lerp(float t, float v1, float v2) {
-	return (1.f - t) * v1 + t * v2;
-}
-
 /* trilinear interpolation */
 static void vol_get_precached_scattering(ShadeInput *shi, float *scatter_col, float *co)
 {
-	const int res = shi->mat->vol_precache_resolution;
-	float voxx, voxy, voxz;
-	int vx, vy, vz;
-	float dx, dy, dz;
-	float d00, d10, d01, d11, d0, d1, d_final;
+	VolumePrecache *vp = shi->obi->volume_precache;
 	float bbmin[3], bbmax[3], dim[3];
-	int rgb;
+	float sample_co[3];
 	
-	if (!shi->obi->volume_precache) return;
+	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);
-	
-	voxx = ((co[0] - bbmin[0]) / dim[0]) * res - 0.5f;
-	voxy = ((co[1] - bbmin[1]) / dim[1]) * res - 0.5f;
-	voxz = ((co[2] - bbmin[2]) / dim[2]) * res - 0.5f;
-	
-	vx = (int)voxx; vy = (int)voxy; vz = (int)voxz;
-	
-	dx = voxx - vx; dy = voxy - vy; dz = voxz - vz;
-	
-	for (rgb=0; rgb < 3; rgb++) {
-		d00 = lerp(dx, D(shi, rgb, vx, vy, vz), 		D(shi, rgb, vx+1, vy, vz));
-		d10 = lerp(dx, D(shi, rgb, vx, vy+1, vz), 		D(shi, rgb, vx+1, vy+1, vz));
-		d01 = lerp(dx, D(shi, rgb, vx, vy, vz+1), 		D(shi, rgb, vx+1, vy, vz+1));
-		d11 = lerp(dx, D(shi, rgb, vx, vy+1, vz+1), 	D(shi, rgb, vx+1, vy+1, vz+1));
-		d0 = lerp(dy, d00, d10);
-		d1 = lerp(dy, d01, d11);
-		d_final = lerp(dz, d0, d1);
-		
-		scatter_col[rgb] = d_final;
-	}
-}
 
-/* no interpolation */
-static void vol_get_precached_scattering_nearest(ShadeInput *shi, float *scatter_col, float *co)
-{
-	const int res = shi->mat->vol_precache_resolution;
-	int x,y,z;
-	float bbmin[3], bbmax[3], dim[3];
+	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]);
 
-	if (!shi->obi->volume_precache) return;
-	
-	VECCOPY(bbmin, shi->obi->obr->boundbox[0]);
-	VECCOPY(bbmax, shi->obi->obr->boundbox[1]);
-	VecSubf(dim, bbmax, bbmin);
-	
-	x = (int)(((co[0] - bbmin[0]) / dim[0]) * res);
-	y = (int)(((co[1] - bbmin[1]) / dim[1]) * res);
-	z = (int)(((co[2] - bbmin[2]) / dim[2]) * res);
-	
-	scatter_col[0] = shi->obi->volume_precache[0*res*res*res + x*res*res + y*res + z];
-	scatter_col[1] = shi->obi->volume_precache[1*res*res*res + x*res*res + y*res + z];
-	scatter_col[2] = shi->obi->volume_precache[2*res*res*res + x*res*res + y*res + z];
+	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);
 }
 
 float vol_get_density(struct ShadeInput *shi, float *co)
diff --git a/source/blender/render/intern/source/voxeldata.c b/source/blender/render/intern/source/voxeldata.c
index f0816cfddff..241b11e04b9 100644
--- a/source/blender/render/intern/source/voxeldata.c
+++ b/source/blender/render/intern/source/voxeldata.c
@@ -34,8 +34,13 @@
 
 #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 "DNA_texture_types.h"
@@ -44,288 +49,79 @@
 #include "texture.h"
 #include "voxeldata.h"
 
-#if defined( _MSC_VER ) && !defined( __cplusplus )
-# define inline __inline
-#endif // defined( _MSC_VER ) && !defined( __cplusplus )
-
-/*---------------------------Utils----------------------------------------*/
-inline int _I(int x, int y, int z, int *n)
-{
-	return (z*(n[1])+y)*(n[2])+x;
-}
-
-float Linear(float xx, float yy, float zz, float *x0, int *n)
-{
-	float sx1,sx0,sy1,sy0,sz1,sz0,v0,v1;
-	int i0,i1,j0,j1,k0,k1;
-	
-	if (xx<0.5) xx=0.5f; if (xx>n[0]+0.5) xx=n[0]+0.5f; i0=(int)xx; i1=i0+1;
-	if (yy<0.5) yy=0.5f; if (yy>n[1]+0.5) yy=n[1]+0.5f; j0=(int)yy; j1=j0+1;
-	if (zz<0.5) zz=0.5f; if (zz>n[2]+0.5) zz=n[2]+0.5f; k0=(int)zz; k1=k0+1;
-	
-	sx1 = xx-i0; sx0 = 1-sx1;
-	sy1 = yy-j0; sy0 = 1-sy1;
-	sz1 = zz-k0; sz0 = 1-sz1;
-	v0 = sx0*(sy0*x0[_I(i0,j0,k0,n)]+sy1*x0[_I(i0,j1,k0,n)])+sx1*(sy0*x0[_I(i1,j0,k0,n)]+sy1*x0[_I(i1,j1,k0,n)]);
-	v1 = sx0*(sy0*x0[_I(i0,j0,k1,n)]+sy1*x0[_I(i0,j1,k1,n)])+sx1*(sy0*x0[_I(i1,j0,k1,n)]+sy1*x0[_I(i1,j1,k1,n)]);
-	return sz0*v0 + sz1*v1;
-}
-
-
-static float D(float *data, int *res, int x, int y, int z)
-{
-	CLAMP(x, 0, res[0]-1);
-	CLAMP(y, 0, res[1]-1);
-	CLAMP(z, 0, res[2]-1);
-	return data[ _I(x, y, z, res) ];
-}
-
-static inline float lerp(float t, float v1, float v2) {
-	return (1.f - t) * v1 + t * v2;
+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);
 }
 
-/* trilinear interpolation */
-static float trilinear(float *data, int *res, float *co)
+void load_frame_raw8(FILE *fp, float *F, int size, int frame)
 {
-	float voxx, voxy, voxz;
-	int vx, vy, vz;
-	float dx, dy, dz;
-	float d00, d10, d01, d11, d0, d1, d_final;
-
-	if (!data) return 0.f;
-	
-	voxx = co[0] * res[0] - 0.5f;
-	voxy = co[1] * res[1] - 0.5f;
-	voxz = co[2] * res[2] - 0.5f;
-	
-	vx = (int)voxx; vy = (int)voxy; vz = (int)voxz;
+	char *tmp;
+	int i;
 	
-	dx = voxx - vx; dy = voxy - vy; dz = voxz - vz;
+	tmp = (char *)MEM_mallocN(sizeof(char)*size, "temporary voxel file reading storage");
 	
-	d00 = lerp(dx, D(data, res, vx, vy, vz), 		D(data, res, vx+1, vy, vz));
-	d10 = lerp(dx, D(data, res, vx, vy+1, vz), 		D(data, res, vx+1, vy+1, vz));
-	d01 = lerp(dx, D(data, res, vx, vy, vz+1), 		D(data, res, vx+1, vy, vz+1));
-	d11 = lerp(dx, D(data, res, vx, vy+1, vz+1), 	D(data, res, vx+1, vy+1, vz+1));
-	d0 = lerp(dy, d00, d10);
-	d1 = lerp(dy, d01, d11);
-	d_final = lerp(dz, d0, d1);
+	fseek(fp,(frame-1)*size*sizeof(char),0);
+	fread(tmp, sizeof(char), size, fp);
 	
-	return d_final;
-}
-
-
-int C[64][64] = {
-	{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 9,-9,-9, 9, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{-6, 6, 6,-6, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{-6, 6, 6,-6, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 4,-4,-4, 4, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0},
-	{-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 9,-9, 0, 0,-9, 9, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{-6, 6, 0, 0, 6,-6, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9, 0, 0,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0},
-	{ 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0},
-	{-27,27,27,-27,27,-27,-27,27,-18,-9,18, 9,18, 9,-18,-9,-18,18,-9, 9,18,-18, 9,-9,-18,18,18,-18,-9, 9, 9,-9,-12,-6,-6,-3,12, 6, 6, 3,-12,-6,12, 6,-6,-3, 6, 3,-12,12,-6, 6,-6, 6,-3, 3,-8,-4,-4,-2,-4,-2,-2,-1},
-	{18,-18,-18,18,-18,18,18,-18, 9, 9,-9,-9,-9,-9, 9, 9,12,-12, 6,-6,-12,12,-6, 6,12,-12,-12,12, 6,-6,-6, 6, 6, 6, 3, 3,-6,-6,-3,-3, 6, 6,-6,-6, 3, 3,-3,-3, 8,-8, 4,-4, 4,-4, 2,-2, 4, 4, 2, 2, 2, 2, 1, 1},
-	{-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0},
-	{18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6, 9,-9, 9,-9,-9, 9,-9, 9,12,-12,-12,12, 6,-6,-6, 6, 6, 3, 6, 3,-6,-3,-6,-3, 8, 4,-8,-4, 4, 2,-4,-2, 6,-6, 6,-6, 3,-3, 3,-3, 4, 2, 4, 2, 2, 1, 2, 1},
-	{-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-6, 6,-6, 6, 6,-6, 6,-6,-8, 8, 8,-8,-4, 4, 4,-4,-3,-3,-3,-3, 3, 3, 3, 3,-4,-4, 4, 4,-2,-2, 2, 2,-4, 4,-4, 4,-2, 2,-2, 2,-2,-2,-2,-2,-1,-1,-1,-1},
-	{ 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{-6, 6, 0, 0, 6,-6, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 4,-4, 0, 0,-4, 4, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0},
-	{-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0},
-	{18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6,12,-12, 6,-6,-12,12,-6, 6, 9,-9,-9, 9, 9,-9,-9, 9, 8, 4, 4, 2,-8,-4,-4,-2, 6, 3,-6,-3, 6, 3,-6,-3, 6,-6, 3,-3, 6,-6, 3,-3, 4, 2, 2, 1, 4, 2, 2, 1},
-	{-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-8, 8,-4, 4, 8,-8, 4,-4,-6, 6, 6,-6,-6, 6, 6,-6,-4,-4,-2,-2, 4, 4, 2, 2,-3,-3, 3, 3,-3,-3, 3, 3,-4, 4,-2, 2,-4, 4,-2, 2,-2,-2,-1,-1,-2,-2,-1,-1},
-	{ 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
-	{ 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0},
-	{-12,12,12,-12,12,-12,-12,12,-8,-4, 8, 4, 8, 4,-8,-4,-6, 6,-6, 6, 6,-6, 6,-6,-6, 6, 6,-6,-6, 6, 6,-6,-4,-2,-4,-2, 4, 2, 4, 2,-4,-2, 4, 2,-4,-2, 4, 2,-3, 3,-3, 3,-3, 3,-3, 3,-2,-1,-2,-1,-2,-1,-2,-1},
-	{ 8,-8,-8, 8,-8, 8, 8,-8, 4, 4,-4,-4,-4,-4, 4, 4, 4,-4, 4,-4,-4, 4,-4, 4, 4,-4,-4, 4, 4,-4,-4, 4, 2, 2, 2, 2,-2,-2,-2,-2, 2, 2,-2,-2, 2, 2,-2,-2, 2,-2, 2,-2, 2,-2, 2,-2, 1, 1, 1, 1, 1, 1, 1, 1}};
-
-int ijk2n(int i, int j, int k) {
-	return(i+4*j+16*k);
-}
-
-void tricubic_get_coeff_stacked(float a[64], float x[64]) {
-	int i,j;
-	for (i=0;i<64;i++) {
-		a[i]=(float)(0.0);
-		for (j=0;j<64;j++) {
-			a[i]+=C[i][j]*x[j];
-		}
-	}
-}
-
-void point2xyz(int p, int *x, int *y, int *z) {
-	switch (p) {
-		case 0: *x=0; *y=0; *z=0; break;
-		case 1: *x=1; *y=0; *z=0; break;
-		case 2: *x=0; *y=1; *z=0; break;
-		case 3: *x=1; *y=1; *z=0; break;
-		case 4: *x=0; *y=0; *z=1; break;
-		case 5: *x=1; *y=0; *z=1; break;
-		case 6: *x=0; *y=1; *z=1; break;
-		case 7: *x=1; *y=1; *z=1; break;
-		default:*x=0; *y=0; *z=0;
+	for (i=0; i<size; i++) {
+		F[i] = (float)tmp[i] / 256.f;
 	}
+	MEM_freeN(tmp);
 }
 
-
-void tricubic_get_coeff(float a[64], float f[8], float dfdx[8], float dfdy[8], float dfdz[8], float d2fdxdy[8], float d2fdxdz[8], float d2fdydz[8], float d3fdxdydz[8]) {
-	int i;
-	float x[64];
-	for (i=0;i<8;i++) {
-		x[0+i]=f[i];
-		x[8+i]=dfdx[i];
-		x[16+i]=dfdy[i];
-		x[24+i]=dfdz[i];
-		x[32+i]=d2fdxdy[i];
-		x[40+i]=d2fdxdz[i];
-		x[48+i]=d2fdydz[i];
-		x[56+i]=d3fdxdydz[i];
+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 r, g, b;
+	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);
 	}
-	tricubic_get_coeff_stacked(a,x);
-}
-
-float tricubic_eval(float a[64], float x, float y, float z) {
-	int i,j,k;
-	float ret=(float)(0.0);
-	
-	for (i=0;i<4;i++) {
-		for (j=0;j<4;j++) {
-			for (k=0;k<4;k++) {
-				ret+=a[ijk2n(i,j,k)]*pow(x,i)*pow(y,j)*pow(z,k);
+	if (!ibuf) return;
+	if (!ibuf->rect_float) IMB_float_from_rect(ibuf);
+	
+	vd->still = 1;
+	vd->resolX = ibuf->x;
+	vd->resolY = ibuf->y;
+	vd->resolZ = iuser->frames;
+	vd->dataset = MEM_mapallocN(sizeof(float)*(vd->resolX)*(vd->resolY)*(vd->resolZ), "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->resolX) ] = (rf[0] + rf[1] + rf[2])*0.333f;
+				rf +=4;
 			}
 		}
 	}
-	return(ret);
-}
-
-
-float tricubic(float xx, float yy, float zz, float *heap, int *n)
-{
-	
-	int xi,yi,zi;
-	float dx,dy,dz;
-	float a[64];
-
-	if (xx<0.5) xx=0.5f; if (xx>n[0]+0.5) xx=n[0]+0.5f; xi=(int)xx;
-	if (yy<0.5) yy=0.5f; if (yy>n[1]+0.5) yy=n[1]+0.5f; yi=(int)yy;
-	if (zz<0.5) zz=0.5f; if (zz>n[2]+0.5) zz=n[2]+0.5f; zi=(int)zz;
-	
-	{
-	float fval[8]={heap[_I(xi,yi,zi,n)],heap[_I(xi+1,yi,zi,n)],heap[_I(xi,yi+1,zi,n)],heap[_I(xi+1,yi+1,zi,n)],heap[_I(xi,yi,zi+1,n)],heap[_I(xi+1,yi,zi+1,n)],heap[_I(xi,yi+1,zi+1,n)],heap[_I(xi+1,yi+1,zi+1,n)]}; 
-	
-	float dfdxval[8]={0.5f*(heap[_I(xi+1,yi,zi,n)]-heap[_I(xi-1,yi,zi,n)]),0.5f*(heap[_I(xi+2,yi,zi,n)]-heap[_I(xi,yi,zi,n)]),
-		0.5f*(heap[_I(xi+1,yi+1,zi,n)]-heap[_I(xi-1,yi+1,zi,n)]),0.5f*(heap[_I(xi+2,yi+1,zi,n)]-heap[_I(xi,yi+1,zi,n)]),
-		0.5f*(heap[_I(xi+1,yi,zi+1,n)]-heap[_I(xi-1,yi,zi+1,n)]),0.5f*(heap[_I(xi+2,yi,zi+1,n)]-heap[_I(xi,yi,zi+1,n)]),
-		0.5f*(heap[_I(xi+1,yi+1,zi+1,n)]-heap[_I(xi-1,yi+1,zi+1,n)]),
-	0.5f*(heap[_I(xi+2,yi+1,zi+1,n)]-heap[_I(xi,yi+1,zi+1,n)])};						
-	
-	float dfdyval[8]={0.5f*(heap[_I(xi,yi+1,zi,n)]-heap[_I(xi,yi-1,zi,n)]),0.5f*(heap[_I(xi+1,yi+1,zi,n)]-heap[_I(xi+1,yi-1,zi,n)]),
-		0.5f*(heap[_I(xi,yi+2,zi,n)]-heap[_I(xi,yi,zi,n)]),0.5f*(heap[_I(xi+1,yi+2,zi,n)]-heap[_I(xi+1,yi,zi,n)]),
-		0.5f*(heap[_I(xi,yi+1,zi+1,n)]-heap[_I(xi,yi-1,zi+1,n)]),0.5f*(heap[_I(xi+1,yi+1,zi+1,n)]-heap[_I(xi+1,yi-1,zi+1,n)]),
-		0.5f*(heap[_I(xi,yi+2,zi+1,n)]-heap[_I(xi,yi,zi+1,n)]),
-	0.5f*(heap[_I(xi+1,yi+2,zi+1,n)]-heap[_I(xi+1,yi,zi+1,n)])};						 
-	
-	float dfdzval[8]={0.5f*(heap[_I(xi,yi,zi+1,n)]-heap[_I(xi,yi,zi-1,n)]),0.5f*(heap[_I(xi+1,yi,zi+1,n)]-heap[_I(xi+1,yi,zi-1,n)]),
-		0.5f*(heap[_I(xi,yi+1,zi+1,n)]-heap[_I(xi,yi+1,zi-1,n)]),0.5f*(heap[_I(xi+1,yi+1,zi+1,n)]-heap[_I(xi+1,yi+1,zi-1,n)]),
-		0.5f*(heap[_I(xi,yi,zi+2,n)]-heap[_I(xi,yi,zi,n)]),0.5f*(heap[_I(xi+1,yi,zi+2,n)]-heap[_I(xi+1,yi,zi,n)]),
-		0.5f*(heap[_I(xi,yi+1,zi+2,n)]-heap[_I(xi,yi+1,zi,n)]),
-	0.5f*(heap[_I(xi+1,yi+1,zi+2,n)]-heap[_I(xi+1,yi+1,zi,n)])};						 
-	
-	float d2fdxdyval[8]={0.25*(heap[_I(xi+1,yi+1,zi,n)]-heap[_I(xi-1,yi+1,zi,n)]-heap[_I(xi+1,yi-1,zi,n)]+heap[_I(xi-1,yi-1,zi,n)]),
-		0.25*(heap[_I(xi+2,yi+1,zi,n)]-heap[_I(xi,yi+1,zi,n)]-heap[_I(xi+2,yi-1,zi,n)]+heap[_I(xi,yi-1,zi,n)]),
-		0.25*(heap[_I(xi+1,yi+2,zi,n)]-heap[_I(xi-1,yi+2,zi,n)]-heap[_I(xi+1,yi,zi,n)]+heap[_I(xi-1,yi,zi,n)]),
-		0.25*(heap[_I(xi+2,yi+2,zi,n)]-heap[_I(xi,yi+2,zi,n)]-heap[_I(xi+2,yi,zi,n)]+heap[_I(xi,yi,zi,n)]),
-		0.25*(heap[_I(xi+1,yi+1,zi+1,n)]-heap[_I(xi-1,yi+1,zi+1,n)]-heap[_I(xi+1,yi-1,zi+1,n)]+heap[_I(xi-1,yi-1,zi+1,n)]),
-		0.25*(heap[_I(xi+2,yi+1,zi+1,n)]-heap[_I(xi,yi+1,zi+1,n)]-heap[_I(xi+2,yi-1,zi+1,n)]+heap[_I(xi,yi-1,zi+1,n)]),
-		0.25*(heap[_I(xi+1,yi+2,zi+1,n)]-heap[_I(xi-1,yi+2,zi+1,n)]-heap[_I(xi+1,yi,zi+1,n)]+heap[_I(xi-1,yi,zi+1,n)]),
-	0.25*(heap[_I(xi+2,yi+2,zi+1,n)]-heap[_I(xi,yi+2,zi+1,n)]-heap[_I(xi+2,yi,zi+1,n)]+heap[_I(xi,yi,zi+1,n)])};						 
-	
-	float d2fdxdzval[8]={0.25f*(heap[_I(xi+1,yi,zi+1,n)]-heap[_I(xi-1,yi,zi+1,n)]-heap[_I(xi+1,yi,zi-1,n)]+heap[_I(xi-1,yi,zi-1,n)]),
-		0.25f*(heap[_I(xi+2,yi,zi+1,n)]-heap[_I(xi,yi,zi+1,n)]-heap[_I(xi+2,yi,zi-1,n)]+heap[_I(xi,yi,zi-1,n)]),
-		0.25f*(heap[_I(xi+1,yi+1,zi+1,n)]-heap[_I(xi-1,yi+1,zi+1,n)]-heap[_I(xi+1,yi+1,zi-1,n)]+heap[_I(xi-1,yi+1,zi-1,n)]),
-		0.25f*(heap[_I(xi+2,yi+1,zi+1,n)]-heap[_I(xi,yi+1,zi+1,n)]-heap[_I(xi+2,yi+1,zi-1,n)]+heap[_I(xi,yi+1,zi-1,n)]),
-		0.25f*(heap[_I(xi+1,yi,zi+2,n)]-heap[_I(xi-1,yi,zi+2,n)]-heap[_I(xi+1,yi,zi,n)]+heap[_I(xi-1,yi,zi,n)]),
-		0.25f*(heap[_I(xi+2,yi,zi+2,n)]-heap[_I(xi,yi,zi+2,n)]-heap[_I(xi+2,yi,zi,n)]+heap[_I(xi,yi,zi,n)]),
-		0.25f*(heap[_I(xi+1,yi+1,zi+2,n)]-heap[_I(xi-1,yi+1,zi+2,n)]-heap[_I(xi+1,yi+1,zi,n)]+heap[_I(xi-1,yi+1,zi,n)]),
-	0.25f*(heap[_I(xi+2,yi+1,zi+2,n)]-heap[_I(xi,yi+1,zi+2,n)]-heap[_I(xi+2,yi+1,zi,n)]+heap[_I(xi,yi+1,zi,n)])};
-	
-	
-	float d2fdydzval[8]={0.25f*(heap[_I(xi,yi+1,zi+1,n)]-heap[_I(xi,yi-1,zi+1,n)]-heap[_I(xi,yi+1,zi-1,n)]+heap[_I(xi,yi-1,zi-1,n)]),
-		0.25f*(heap[_I(xi+1,yi+1,zi+1,n)]-heap[_I(xi+1,yi-1,zi+1,n)]-heap[_I(xi+1,yi+1,zi-1,n)]+heap[_I(xi+1,yi-1,zi-1,n)]),
-		0.25f*(heap[_I(xi,yi+2,zi+1,n)]-heap[_I(xi,yi,zi+1,n)]-heap[_I(xi,yi+2,zi-1,n)]+heap[_I(xi,yi,zi-1,n)]),
-		0.25f*(heap[_I(xi+1,yi+2,zi+1,n)]-heap[_I(xi+1,yi,zi+1,n)]-heap[_I(xi+1,yi+2,zi-1,n)]+heap[_I(xi+1,yi,zi-1,n)]),
-		0.25f*(heap[_I(xi,yi+1,zi+2,n)]-heap[_I(xi,yi-1,zi+2,n)]-heap[_I(xi,yi+1,zi,n)]+heap[_I(xi,yi-1,zi,n)]),
-		0.25f*(heap[_I(xi+1,yi+1,zi+2,n)]-heap[_I(xi+1,yi-1,zi+2,n)]-heap[_I(xi+1,yi+1,zi,n)]+heap[_I(xi+1,yi-1,zi,n)]),
-		0.25f*(heap[_I(xi,yi+2,zi+2,n)]-heap[_I(xi,yi,zi+2,n)]-heap[_I(xi,yi+2,zi,n)]+heap[_I(xi,yi,zi,n)]),
-	0.25f*(heap[_I(xi+1,yi+2,zi+2,n)]-heap[_I(xi+1,yi,zi+2,n)]-heap[_I(xi+1,yi+2,zi,n)]+heap[_I(xi+1,yi,zi,n)])};
-	
-	
-	float d3fdxdydzval[8]={0.125f*(heap[_I(xi+1,yi+1,zi+1,n)]-heap[_I(xi-1,yi+1,zi+1,n)]-heap[_I(xi+1,yi-1,zi+1,n)]+heap[_I(xi-1,yi-1,zi+1,n)]-heap[_I(xi+1,yi+1,zi-1,n)]+heap[_I(xi-1,yi+1,zi-1,n)]+heap[_I(xi+1,yi-1,zi-1,n)]-heap[_I(xi-1,yi-1,zi-1,n)]),
-		0.125f*(heap[_I(xi+2,yi+1,zi+1,n)]-heap[_I(xi,yi+1,zi+1,n)]-heap[_I(xi+2,yi-1,zi+1,n)]+heap[_I(xi,yi-1,zi+1,n)]-heap[_I(xi+2,yi+1,zi-1,n)]+heap[_I(xi,yi+1,zi-1,n)]+heap[_I(xi+2,yi-1,zi-1,n)]-heap[_I(xi,yi-1,zi-1,n)]),
-		0.125f*(heap[_I(xi+1,yi+2,zi+1,n)]-heap[_I(xi-1,yi+2,zi+1,n)]-heap[_I(xi+1,yi,zi+1,n)]+heap[_I(xi-1,yi,zi+1,n)]-heap[_I(xi+1,yi+2,zi-1,n)]+heap[_I(xi-1,yi+2,zi-1,n)]+heap[_I(xi+1,yi,zi-1,n)]-heap[_I(xi-1,yi,zi-1,n)]),
-		0.125f*(heap[_I(xi+2,yi+2,zi+1,n)]-heap[_I(xi,yi+2,zi+1,n)]-heap[_I(xi+2,yi,zi+1,n)]+heap[_I(xi,yi,zi+1,n)]-heap[_I(xi+2,yi+2,zi-1,n)]+heap[_I(xi,yi+2,zi-1,n)]+heap[_I(xi+2,yi,zi-1,n)]-heap[_I(xi,yi,zi-1,n)]),
-		0.125f*(heap[_I(xi+1,yi+1,zi+2,n)]-heap[_I(xi-1,yi+1,zi+2,n)]-heap[_I(xi+1,yi-1,zi+2,n)]+heap[_I(xi-1,yi-1,zi+2,n)]-heap[_I(xi+1,yi+1,zi,n)]+heap[_I(xi-1,yi+1,zi,n)]+heap[_I(xi+1,yi-1,zi,n)]-heap[_I(xi-1,yi-1,zi,n)]),
-		0.125f*(heap[_I(xi+2,yi+1,zi+2,n)]-heap[_I(xi,yi+1,zi+2,n)]-heap[_I(xi+2,yi-1,zi+2,n)]+heap[_I(xi,yi-1,zi+2,n)]-heap[_I(xi+2,yi+1,zi,n)]+heap[_I(xi,yi+1,zi,n)]+heap[_I(xi+2,yi-1,zi,n)]-heap[_I(xi,yi-1,zi,n)]),
-		0.125f*(heap[_I(xi+1,yi+2,zi+2,n)]-heap[_I(xi-1,yi+2,zi+2,n)]-heap[_I(xi+1,yi,zi+2,n)]+heap[_I(xi-1,yi,zi+2,n)]-heap[_I(xi+1,yi+2,zi,n)]+heap[_I(xi-1,yi+2,zi,n)]+heap[_I(xi+1,yi,zi,n)]-heap[_I(xi-1,yi,zi,n)]),
-	0.125f*(heap[_I(xi+2,yi+2,zi+2,n)]-heap[_I(xi,yi+2,zi+2,n)]-heap[_I(xi+2,yi,zi+2,n)]+heap[_I(xi,yi,zi+2,n)]-heap[_I(xi+2,yi+2,zi,n)]+heap[_I(xi,yi+2,zi,n)]+heap[_I(xi+2,yi,zi,n)]-heap[_I(xi,yi,zi,n)])};
-
-	
-	tricubic_get_coeff(a,fval,dfdxval,dfdyval,dfdzval,d2fdxdyval,d2fdxdzval,d2fdydzval,d3fdxdydzval);
-	}
-		
-	dx = xx-xi;
-	dy = yy-yi;
-	dz = zz-zi;
-	
-	return tricubic_eval(a,dx,dy,dz);
-	
-}
-
-void load_frame (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 write_voxeldata_header(struct VoxelDataHeader *h, FILE *fp)
@@ -352,20 +148,37 @@ 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;
+	}
+
 	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);
 	
-	read_voxeldata_header(fp, vd);
 	size = (vd->resolX)*(vd->resolY)*(vd->resolZ);
-	vd->dataset = MEM_mallocN(sizeof(float)*size, "voxel dataset");
-	
-	//here improve the dataset loading function for more dataset types
-	if (vd->still) load_frame(fp, vd->dataset, size, vd->still_frame, sizeof(VoxelDataHeader));
-	else load_frame(fp, vd->dataset, size, re->r.cfra, sizeof(VoxelDataHeader));
+	vd->dataset = MEM_mapallocN(sizeof(float)*size, "voxel dataset");
+		
+	if (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);
 }
@@ -380,6 +193,7 @@ void make_voxeldata(struct Render *re)
 	re->i.infostr= "Loading voxel datasets";
 	re->stats_draw(&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);
@@ -420,69 +234,59 @@ int voxeldatatex(struct Tex *tex, float *texvec, struct TexResult *texres)
 {	 
     int retval = TEX_INT;
 	VoxelData *vd = tex->vd;	
-	float vec[3] = {0.0, 0.0, 0.0};	
-	float co[3];
-	float dx, dy, dz;
-	int xi, yi, zi;
-	float xf, yf, zf;
-	int i=0, fail=0;
-	int resol[3];
-	
+	float co[3], offset[3] = {0.5, 0.5, 0.5};
+
 	if ((!vd) || (vd->dataset==NULL)) {
 		texres->tin = 0.0f;
 		return 0;
 	}
 	
-	resol[0] = vd->resolX;
-	resol[1] = vd->resolY;
-	resol[2] = vd->resolZ;
-	
-	VECCOPY(co, texvec);	
-	
-	dx=1.0f/(resol[0]);
-	dy=1.0f/(resol[1]);
-	dz=1.0f/(resol[2]);
-	
-	xi=co[0]/dx;
-	yi=co[1]/dy;
-	zi=co[2]/dz;
-		
-	xf=co[0]/dx;
-	yf=co[1]/dy;
-	zf=co[2]/dz;
-	
-	if (xi>1 && xi<resol[0])
-	{
-		if (yi>1 && yi<resol[1])
+	/* 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 (zi>1 && zi<resol[2])
-			{
-				switch (vd->interp_type)
-				{
-					case TEX_VD_NEARESTNEIGHBOR:
-					{
-						texres->tin = vd->dataset[_I(xi,yi,zi,resol)];
-						break;  
-					}
-					case TEX_VD_LINEAR:
-					{
-						texres->tin = trilinear(vd->dataset, resol, co);
-						break;					
-					}
-					case TEX_VD_TRICUBIC:
-					{
-						texres->tin = tricubic(xf, yf, zf, vd->dataset, resol);
-						break;
-					}
-				}				  				   
-			} else fail++;
-		} else fail++;
-	} else fail++;
+			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;
+		}
+	}
 	
-	if (fail) texres->tin=0.0f;
-
-	texres->tin *= vd->int_multiplier;
+	switch (vd->interp_type) {
+		case TEX_VD_NEARESTNEIGHBOR:
+			texres->tin = voxel_sample_nearest(vd->dataset, &vd->resolX, co);
+			break;  
+		case TEX_VD_LINEAR:
+			texres->tin = voxel_sample_trilinear(vd->dataset, &vd->resolX, co);
+			break;					
+		case TEX_VD_TRICUBIC:
+			texres->tin = voxel_sample_tricubic(vd->dataset, &vd->resolX, co);
+			break;
+	}
 	
+	texres->tin *= vd->int_multiplier;
 	BRICONT;
 	
 	texres->tr = texres->tin;