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.

10 files changed, 806 insertions, 586 deletions

diff --git a/source/blender/blenlib/BLI_voxel.h b/source/blender/blenlib/BLI_voxel.h
new file mode 100644
index 00000000000..091d8e3682d
--- /dev/null
+++ b/source/blender/blenlib/BLI_voxel.h
@@ -0,0 +1,40 @@
+/**
+ *
+ * ***** 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 *****
+ */
+
+#ifndef BLI_VOXEL_H
+#define BLI_VOXEL_H
+
+/* find the index number of a voxel, given x/y/z integer coords and resolution vector */
+#define V_I(x, y, z, res) ( (z)*(res)[1]*(res)[0] + (y)*(res)[0] + (x) )
+
+/* all input coordinates must be in bounding box 0.0 - 1.0 */
+float voxel_sample_nearest(float *data, int *res, float *co);
+float voxel_sample_trilinear(float *data, int *res, float *co);
+float voxel_sample_tricubic(float *data, int *res, float *co);
+
+#endif /* BLI_VOXEL_H */
\ No newline at end of file
diff --git a/source/blender/blenlib/intern/voxel.c b/source/blender/blenlib/intern/voxel.c
new file mode 100644
index 00000000000..8a1ca6e59ff
--- /dev/null
+++ b/source/blender/blenlib/intern/voxel.c
@@ -0,0 +1,311 @@
+/**
+ *
+ * ***** 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 "BLI_voxel.h"
+
+#include "BKE_utildefines.h"
+
+
+#if defined( _MSC_VER ) && !defined( __cplusplus )
+# define inline __inline
+#endif // defined( _MSC_VER ) && !defined( __cplusplus )
+
+static inline 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[ V_I(x, y, z, res) ];
+}
+
+/* *** nearest neighbour *** */
+/* input coordinates must be in bounding box 0.0 - 1.0 */
+float voxel_sample_nearest(float *data, int *res, float *co)
+{
+	int xi, yi, zi;
+	
+	xi = co[0] * res[0];
+	yi = co[1] * res[1];
+	zi = co[2] * res[2];
+	
+	return D(data, res, xi, yi, zi);
+}
+
+
+/* *** trilinear *** */
+/* input coordinates must be in bounding box 0.0 - 1.0 */
+
+static inline float lerp(float t, float v1, float v2) {
+	return (1.f - t) * v1 + t * v2;
+}
+
+/* trilinear interpolation - taken partly from pbrt's implementation: http://www.pbrt.org */
+float voxel_sample_trilinear(float *data, int *res, float *co)
+{
+	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;
+	
+	dx = voxx - vx; dy = voxy - vy; dz = voxz - vz;
+	
+	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);
+	
+	return d_final;
+}
+
+/* *** tricubic *** */
+
+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}};
+
+static int ijk2n(int i, int j, int k) {
+	return(i+4*j+16*k);
+}
+
+static 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];
+		}
+	}
+}
+
+static 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;
+	}
+}
+
+
+static 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];
+	}
+	tricubic_get_coeff_stacked(a,x);
+}
+
+static 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);
+			}
+		}
+	}
+	return(ret);
+}
+
+/* tricubic interpolation
+ * from 'libtricubic': http://www.lekien.com/~francois/software/tricubic/ 
+ * input coordinates must be in bounding box 0.0 - 1.0 */
+float voxel_sample_tricubic(float *data, int *res, float *co)
+{
+	float xx, yy, zz;
+	int xi,yi,zi;
+	int *n = res;
+	float dx,dy,dz;
+	float a[64];
+	
+	xx = co[0] * res[0] - 0.5f;
+	yy = co[1] * res[1] - 0.5f;
+	zz = co[2] * res[2] - 0.5f;
+	
+	xi = (int)xx; yi = (int)yy; zi = (int)zz;
+	
+	{
+		float fval[8]={data[V_I(xi,yi,zi,n)],data[V_I(xi+1,yi,zi,n)],data[V_I(xi,yi+1,zi,n)],data[V_I(xi+1,yi+1,zi,n)],data[V_I(xi,yi,zi+1,n)],data[V_I(xi+1,yi,zi+1,n)],data[V_I(xi,yi+1,zi+1,n)],data[V_I(xi+1,yi+1,zi+1,n)]}; 
+		
+		float dfdxval[8]={0.5f*(data[V_I(xi+1,yi,zi,n)]-data[V_I(xi-1,yi,zi,n)]),0.5f*(data[V_I(xi+2,yi,zi,n)]-data[V_I(xi,yi,zi,n)]),
+			0.5f*(data[V_I(xi+1,yi+1,zi,n)]-data[V_I(xi-1,yi+1,zi,n)]),0.5f*(data[V_I(xi+2,yi+1,zi,n)]-data[V_I(xi,yi+1,zi,n)]),
+			0.5f*(data[V_I(xi+1,yi,zi+1,n)]-data[V_I(xi-1,yi,zi+1,n)]),0.5f*(data[V_I(xi+2,yi,zi+1,n)]-data[V_I(xi,yi,zi+1,n)]),
+			0.5f*(data[V_I(xi+1,yi+1,zi+1,n)]-data[V_I(xi-1,yi+1,zi+1,n)]),
+			0.5f*(data[V_I(xi+2,yi+1,zi+1,n)]-data[V_I(xi,yi+1,zi+1,n)])};						
+		
+		float dfdyval[8]={0.5f*(data[V_I(xi,yi+1,zi,n)]-data[V_I(xi,yi-1,zi,n)]),0.5f*(data[V_I(xi+1,yi+1,zi,n)]-data[V_I(xi+1,yi-1,zi,n)]),
+			0.5f*(data[V_I(xi,yi+2,zi,n)]-data[V_I(xi,yi,zi,n)]),0.5f*(data[V_I(xi+1,yi+2,zi,n)]-data[V_I(xi+1,yi,zi,n)]),
+			0.5f*(data[V_I(xi,yi+1,zi+1,n)]-data[V_I(xi,yi-1,zi+1,n)]),0.5f*(data[V_I(xi+1,yi+1,zi+1,n)]-data[V_I(xi+1,yi-1,zi+1,n)]),
+			0.5f*(data[V_I(xi,yi+2,zi+1,n)]-data[V_I(xi,yi,zi+1,n)]),
+			0.5f*(data[V_I(xi+1,yi+2,zi+1,n)]-data[V_I(xi+1,yi,zi+1,n)])};						 
+		
+		float dfdzval[8]={0.5f*(data[V_I(xi,yi,zi+1,n)]-data[V_I(xi,yi,zi-1,n)]),0.5f*(data[V_I(xi+1,yi,zi+1,n)]-data[V_I(xi+1,yi,zi-1,n)]),
+			0.5f*(data[V_I(xi,yi+1,zi+1,n)]-data[V_I(xi,yi+1,zi-1,n)]),0.5f*(data[V_I(xi+1,yi+1,zi+1,n)]-data[V_I(xi+1,yi+1,zi-1,n)]),
+			0.5f*(data[V_I(xi,yi,zi+2,n)]-data[V_I(xi,yi,zi,n)]),0.5f*(data[V_I(xi+1,yi,zi+2,n)]-data[V_I(xi+1,yi,zi,n)]),
+			0.5f*(data[V_I(xi,yi+1,zi+2,n)]-data[V_I(xi,yi+1,zi,n)]),
+			0.5f*(data[V_I(xi+1,yi+1,zi+2,n)]-data[V_I(xi+1,yi+1,zi,n)])};						 
+		
+		float d2fdxdyval[8]={0.25*(data[V_I(xi+1,yi+1,zi,n)]-data[V_I(xi-1,yi+1,zi,n)]-data[V_I(xi+1,yi-1,zi,n)]+data[V_I(xi-1,yi-1,zi,n)]),
+			0.25*(data[V_I(xi+2,yi+1,zi,n)]-data[V_I(xi,yi+1,zi,n)]-data[V_I(xi+2,yi-1,zi,n)]+data[V_I(xi,yi-1,zi,n)]),
+			0.25*(data[V_I(xi+1,yi+2,zi,n)]-data[V_I(xi-1,yi+2,zi,n)]-data[V_I(xi+1,yi,zi,n)]+data[V_I(xi-1,yi,zi,n)]),
+			0.25*(data[V_I(xi+2,yi+2,zi,n)]-data[V_I(xi,yi+2,zi,n)]-data[V_I(xi+2,yi,zi,n)]+data[V_I(xi,yi,zi,n)]),
+			0.25*(data[V_I(xi+1,yi+1,zi+1,n)]-data[V_I(xi-1,yi+1,zi+1,n)]-data[V_I(xi+1,yi-1,zi+1,n)]+data[V_I(xi-1,yi-1,zi+1,n)]),
+			0.25*(data[V_I(xi+2,yi+1,zi+1,n)]-data[V_I(xi,yi+1,zi+1,n)]-data[V_I(xi+2,yi-1,zi+1,n)]+data[V_I(xi,yi-1,zi+1,n)]),
+			0.25*(data[V_I(xi+1,yi+2,zi+1,n)]-data[V_I(xi-1,yi+2,zi+1,n)]-data[V_I(xi+1,yi,zi+1,n)]+data[V_I(xi-1,yi,zi+1,n)]),
+			0.25*(data[V_I(xi+2,yi+2,zi+1,n)]-data[V_I(xi,yi+2,zi+1,n)]-data[V_I(xi+2,yi,zi+1,n)]+data[V_I(xi,yi,zi+1,n)])};						 
+		
+		float d2fdxdzval[8]={0.25f*(data[V_I(xi+1,yi,zi+1,n)]-data[V_I(xi-1,yi,zi+1,n)]-data[V_I(xi+1,yi,zi-1,n)]+data[V_I(xi-1,yi,zi-1,n)]),
+			0.25f*(data[V_I(xi+2,yi,zi+1,n)]-data[V_I(xi,yi,zi+1,n)]-data[V_I(xi+2,yi,zi-1,n)]+data[V_I(xi,yi,zi-1,n)]),
+			0.25f*(data[V_I(xi+1,yi+1,zi+1,n)]-data[V_I(xi-1,yi+1,zi+1,n)]-data[V_I(xi+1,yi+1,zi-1,n)]+data[V_I(xi-1,yi+1,zi-1,n)]),
+			0.25f*(data[V_I(xi+2,yi+1,zi+1,n)]-data[V_I(xi,yi+1,zi+1,n)]-data[V_I(xi+2,yi+1,zi-1,n)]+data[V_I(xi,yi+1,zi-1,n)]),
+			0.25f*(data[V_I(xi+1,yi,zi+2,n)]-data[V_I(xi-1,yi,zi+2,n)]-data[V_I(xi+1,yi,zi,n)]+data[V_I(xi-1,yi,zi,n)]),
+			0.25f*(data[V_I(xi+2,yi,zi+2,n)]-data[V_I(xi,yi,zi+2,n)]-data[V_I(xi+2,yi,zi,n)]+data[V_I(xi,yi,zi,n)]),
+			0.25f*(data[V_I(xi+1,yi+1,zi+2,n)]-data[V_I(xi-1,yi+1,zi+2,n)]-data[V_I(xi+1,yi+1,zi,n)]+data[V_I(xi-1,yi+1,zi,n)]),
+			0.25f*(data[V_I(xi+2,yi+1,zi+2,n)]-data[V_I(xi,yi+1,zi+2,n)]-data[V_I(xi+2,yi+1,zi,n)]+data[V_I(xi,yi+1,zi,n)])};
+		
+		
+		float d2fdydzval[8]={0.25f*(data[V_I(xi,yi+1,zi+1,n)]-data[V_I(xi,yi-1,zi+1,n)]-data[V_I(xi,yi+1,zi-1,n)]+data[V_I(xi,yi-1,zi-1,n)]),
+			0.25f*(data[V_I(xi+1,yi+1,zi+1,n)]-data[V_I(xi+1,yi-1,zi+1,n)]-data[V_I(xi+1,yi+1,zi-1,n)]+data[V_I(xi+1,yi-1,zi-1,n)]),
+			0.25f*(data[V_I(xi,yi+2,zi+1,n)]-data[V_I(xi,yi,zi+1,n)]-data[V_I(xi,yi+2,zi-1,n)]+data[V_I(xi,yi,zi-1,n)]),
+			0.25f*(data[V_I(xi+1,yi+2,zi+1,n)]-data[V_I(xi+1,yi,zi+1,n)]-data[V_I(xi+1,yi+2,zi-1,n)]+data[V_I(xi+1,yi,zi-1,n)]),
+			0.25f*(data[V_I(xi,yi+1,zi+2,n)]-data[V_I(xi,yi-1,zi+2,n)]-data[V_I(xi,yi+1,zi,n)]+data[V_I(xi,yi-1,zi,n)]),
+			0.25f*(data[V_I(xi+1,yi+1,zi+2,n)]-data[V_I(xi+1,yi-1,zi+2,n)]-data[V_I(xi+1,yi+1,zi,n)]+data[V_I(xi+1,yi-1,zi,n)]),
+			0.25f*(data[V_I(xi,yi+2,zi+2,n)]-data[V_I(xi,yi,zi+2,n)]-data[V_I(xi,yi+2,zi,n)]+data[V_I(xi,yi,zi,n)]),
+			0.25f*(data[V_I(xi+1,yi+2,zi+2,n)]-data[V_I(xi+1,yi,zi+2,n)]-data[V_I(xi+1,yi+2,zi,n)]+data[V_I(xi+1,yi,zi,n)])};
+		
+		
+		float d3fdxdydzval[8]={0.125f*(data[V_I(xi+1,yi+1,zi+1,n)]-data[V_I(xi-1,yi+1,zi+1,n)]-data[V_I(xi+1,yi-1,zi+1,n)]+data[V_I(xi-1,yi-1,zi+1,n)]-data[V_I(xi+1,yi+1,zi-1,n)]+data[V_I(xi-1,yi+1,zi-1,n)]+data[V_I(xi+1,yi-1,zi-1,n)]-data[V_I(xi-1,yi-1,zi-1,n)]),
+			0.125f*(data[V_I(xi+2,yi+1,zi+1,n)]-data[V_I(xi,yi+1,zi+1,n)]-data[V_I(xi+2,yi-1,zi+1,n)]+data[V_I(xi,yi-1,zi+1,n)]-data[V_I(xi+2,yi+1,zi-1,n)]+data[V_I(xi,yi+1,zi-1,n)]+data[V_I(xi+2,yi-1,zi-1,n)]-data[V_I(xi,yi-1,zi-1,n)]),
+			0.125f*(data[V_I(xi+1,yi+2,zi+1,n)]-data[V_I(xi-1,yi+2,zi+1,n)]-data[V_I(xi+1,yi,zi+1,n)]+data[V_I(xi-1,yi,zi+1,n)]-data[V_I(xi+1,yi+2,zi-1,n)]+data[V_I(xi-1,yi+2,zi-1,n)]+data[V_I(xi+1,yi,zi-1,n)]-data[V_I(xi-1,yi,zi-1,n)]),
+			0.125f*(data[V_I(xi+2,yi+2,zi+1,n)]-data[V_I(xi,yi+2,zi+1,n)]-data[V_I(xi+2,yi,zi+1,n)]+data[V_I(xi,yi,zi+1,n)]-data[V_I(xi+2,yi+2,zi-1,n)]+data[V_I(xi,yi+2,zi-1,n)]+data[V_I(xi+2,yi,zi-1,n)]-data[V_I(xi,yi,zi-1,n)]),
+			0.125f*(data[V_I(xi+1,yi+1,zi+2,n)]-data[V_I(xi-1,yi+1,zi+2,n)]-data[V_I(xi+1,yi-1,zi+2,n)]+data[V_I(xi-1,yi-1,zi+2,n)]-data[V_I(xi+1,yi+1,zi,n)]+data[V_I(xi-1,yi+1,zi,n)]+data[V_I(xi+1,yi-1,zi,n)]-data[V_I(xi-1,yi-1,zi,n)]),
+			0.125f*(data[V_I(xi+2,yi+1,zi+2,n)]-data[V_I(xi,yi+1,zi+2,n)]-data[V_I(xi+2,yi-1,zi+2,n)]+data[V_I(xi,yi-1,zi+2,n)]-data[V_I(xi+2,yi+1,zi,n)]+data[V_I(xi,yi+1,zi,n)]+data[V_I(xi+2,yi-1,zi,n)]-data[V_I(xi,yi-1,zi,n)]),
+			0.125f*(data[V_I(xi+1,yi+2,zi+2,n)]-data[V_I(xi-1,yi+2,zi+2,n)]-data[V_I(xi+1,yi,zi+2,n)]+data[V_I(xi-1,yi,zi+2,n)]-data[V_I(xi+1,yi+2,zi,n)]+data[V_I(xi-1,yi+2,zi,n)]+data[V_I(xi+1,yi,zi,n)]-data[V_I(xi-1,yi,zi,n)]),
+			0.125f*(data[V_I(xi+2,yi+2,zi+2,n)]-data[V_I(xi,yi+2,zi+2,n)]-data[V_I(xi+2,yi,zi+2,n)]+data[V_I(xi,yi,zi+2,n)]-data[V_I(xi+2,yi+2,zi,n)]+data[V_I(xi,yi+2,zi,n)]+data[V_I(xi+2,yi,zi,n)]-data[V_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);
+	
+}
+
diff --git a/source/blender/makesdna/DNA_texture_types.h b/source/blender/makesdna/DNA_texture_types.h
index 420850ce396..5cbad6cf424 100644
--- a/source/blender/makesdna/DNA_texture_types.h
+++ b/source/blender/makesdna/DNA_texture_types.h
@@ -165,9 +165,10 @@ typedef struct PointDensity {
 typedef struct VoxelData {
 	int resolX, resolY, resolZ;
 	int interp_type;
-
+	short file_format;
+	short flag;
+	
 	float int_multiplier;
-	float vxpad;
 	
 	int still, still_frame;
 	char source_path[240];
@@ -480,13 +481,20 @@ typedef struct TexMapping {
 #define POINT_DATA_VEL		1
 #define POINT_DATA_LIFE		2
 
-/******************** Voxel Data *****************************/ 
-#define TEX_VD_CUBIC              0
-#define TEX_VD_PARALLELOGRAM      1
- 
+/******************** Voxel Data *****************************/
+/* flag */
+
+
+/* interpolation */
 #define TEX_VD_NEARESTNEIGHBOR		0
-#define TEX_VD_LINEAR			1
+#define TEX_VD_LINEAR				1
 #define TEX_VD_TRICUBIC				2
 
+/* file format */
+#define TEX_VD_BLENDERVOXEL		0
+#define TEX_VD_RAW_8BIT			1
+#define TEX_VD_RAW_16BIT		2
+#define TEX_VD_IMAGE_SEQUENCE	3
+
 #endif
 
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;
diff --git a/source/blender/src/buttons_shading.c b/source/blender/src/buttons_shading.c
index f74afde5e81..54fc639059d 100644
--- a/source/blender/src/buttons_shading.c
+++ b/source/blender/src/buttons_shading.c
@@ -1053,67 +1053,7 @@ static void texture_panel_pointdensity_modify(Tex *tex)
 
 }
 
-static void texture_panel_voxeldata(Tex *tex)
-{
-	uiBlock *block;
-	VoxelData *vd;
-	short yco=PANEL_YMAX;
-
-	block= uiNewBlock(&curarea->uiblocks, "texture_panel_voxeldata", UI_EMBOSS, UI_HELV, curarea->win);
-	if(uiNewPanel(curarea, block, "Voxel Data", "Texture", PANELX, PANELY, PANELW, PANELH+YSPACE)==0) return;
-	uiSetButLock(tex->id.lib!=0, ERROR_LIBDATA_MESSAGE);
-
-	if(tex->vd==NULL) {
-		tex->vd= BKE_add_voxeldata();
-	}
-
-	if(tex->vd) {
-		vd= tex->vd;
-		
-		uiDefBut(block, LABEL, B_NOP, "Data source:",
-			X2CLM1, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
-			
-		uiDefIconTextBut(block, BUT, B_VOXELDATA_LOAD, ICON_FILESEL, "Open",
-			X4CLM1, yco-=BUTH, BUTW4, BUTH, 0, 0, 0, 0, 0, "");
-		uiDefBut(block, TEX, 0, "",
-			X4CLM2+XSPACE, yco, BUTW2+BUTW4+2*XSPACE, BUTH, &vd->source_path, 0.0, 79.0, 0, 0, "File path to the voxel data set");
 
-		yco -= YSPACE;
-		
-		uiDefBut(block, LABEL, B_NOP, "Interpolation:",
-			X2CLM1, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
-		uiDefButI(block, MENU, B_REDR, "None %x0|Linear %x1|Tricubic %x2",
-			X2CLM1, yco-=BUTH, BUTW2, BUTH, &vd->interp_type, 0.0, 0.0, 0, 0, "Interpolation type");		
-		
-		yco -= YSPACE;
-		
-		uiDefButF(block, NUM, B_REDR, "Intensity: ",
-			X2CLM1, yco-=BUTH, BUTW2, BUTH, &(vd->int_multiplier), 0.0001, 10000.0, 0, 0, "Multiplier to scale up or down the texture's intensity");
-		
-		yco = PANEL_YMAX - 2*BUTH - 2*YSPACE;
-		
-		uiDefBut(block, LABEL, B_NOP, "Resolution:",
-			X2CLM2, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
-		uiBlockBeginAlign(block);
-		uiDefButI(block, NUM, B_REDR, "X: ",
-			X2CLM2, yco-=BUTH, BUTW2, BUTH, &(vd->resolX), 1, 10000, 0, 0, "Resolution of the voxel data");
-		uiDefButI(block, NUM, B_REDR, "Y: ",
-			X2CLM2, yco-=BUTH, BUTW2, BUTH, &(vd->resolY), 1, 10000, 0, 0, "Resolution of the voxel data");
-		uiDefButI(block, NUM, B_REDR, "Z: ",
-			X2CLM2, yco-= BUTH, BUTW2, BUTH, &(vd->resolZ), 1, 10000, 0, 0, "Resolution of the voxel data");
-		uiBlockEndAlign(block);
-		
-		yco -= YSPACE;
-		
-		uiBlockBeginAlign(block);
-		uiDefButI(block,TOG, B_REDR, "Still",
-			X2CLM2, yco-=BUTH, BUTW2, BUTH, &(vd->still), 0,1, 0, 0, "Use a still frame from the data sequence for the entire rendered animation");
-		if (vd->still) uiDefButI(block, NUM, B_REDR, "Frame: ",
-			X2CLM2, yco-=BUTH, BUTW2, BUTH, &(vd->still_frame), 1, 10000, 0, 0, "The frame to pause on for the entire rendered animation");	
-		uiBlockEndAlign(block);
-		
-	}
-}
 
 static char *layer_menu(RenderResult *rr, short *curlay)
 {
@@ -1791,6 +1731,144 @@ static void texture_panel_envmap(Tex *tex)
 	}
 }
 
+static void texture_panel_voxeldata(Tex *tex)
+{
+	uiBlock *block;
+	VoxelData *vd;
+	short yco=PANEL_YMAX;
+	
+	block= uiNewBlock(&curarea->uiblocks, "texture_panel_voxeldata", UI_EMBOSS, UI_HELV, curarea->win);
+	if(uiNewPanel(curarea, block, "Voxel Data", "Texture", PANELX, PANELY, PANELW, PANELH+BUTH+2*YSPACE)==0) return;
+	uiSetButLock(tex->id.lib!=0, ERROR_LIBDATA_MESSAGE);
+	
+	if(tex->vd==NULL) {
+		tex->vd= BKE_add_voxeldata();
+	}
+	
+	if(tex->vd) {
+		vd= tex->vd;
+		
+		uiDefBut(block, LABEL, B_NOP, "Data source:",
+				 X2CLM1, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
+		
+		if (vd->file_format != TEX_VD_IMAGE_SEQUENCE) {
+			uiDefIconTextBut(block, BUT, B_VOXELDATA_LOAD, ICON_FILESEL, "Open",
+							 X4CLM1, yco-=BUTH, BUTW4, BUTH, 0, 0, 0, 0, 0, "");
+			uiDefBut(block, TEX, 0, "",
+					 X4CLM2+XSPACE, yco, BUTW2+BUTW4+2*XSPACE, BUTH, &vd->source_path, 0.0, 79.0, 0, 0, "File path to the voxel data set");
+		}
+		else {
+			char *strp;
+			char str[128];
+			ImageUser *iuser = &(tex->iuser);
+			uiBut *but;
+			
+			/* Browse */
+			IMAnames_to_pupstring(&strp, NULL, NULL, &(G.main->image), NULL, &iuser->menunr);
+			
+			uiBlockBeginAlign(block);
+			but= uiDefButS(block, MENU, B_REDR, strp,
+						   X2CLM1, yco-=BUTH, ICONBUTW, BUTH, &iuser->menunr, 0, 0, 0, 0, "Selects an existing Image or Movie");
+			uiButSetFunc(but, image_browse_cb, &(tex->ima), iuser);
+			
+			MEM_freeN(strp);
+			
+			if (tex->ima) {
+				uiSetButLock(tex->ima->id.lib!=NULL, ERROR_LIBDATA_MESSAGE);
+				but= uiDefIconBut(block, BUT, B_REDR, ICON_FILESEL,
+								  X2CLM1+ICONBUTW, yco, X2CLM1+ICONBUTW, BUTH, 0, 0, 0, 0, 0, "Open Fileselect to load new Image");
+				uiButSetFunc(but, image_load_fs_cb, &(tex->ima), iuser);
+				
+				but= uiDefBut(block, TEX, B_IDNAME, "IM:",
+							  X2CLM1+2*ICONBUTW, yco, PANEL_XMAX-4*ICONBUTW-BUTW4, BUTH, tex->ima->id.name+2, 0.0, 21.0, 0, 0, "Current Image Datablock name.");
+				uiButSetFunc(but, test_idbutton_cb, tex->ima->id.name, NULL);
+				
+				but= uiDefBut(block, BUT, B_REDR, "Reload",
+					PANEL_XMAX-2*ICONBUTW-BUTW4, yco, BUTW4, BUTH, NULL, 0, 0, 0, 0, "Reloads Image or Movie");
+				uiButSetFunc(but, image_reload_cb, &(tex->ima), iuser);
+				
+				but= uiDefIconBut(block, BUT, B_REDR, ICON_X,
+					PANEL_XMAX-2*ICONBUTW, yco, ICONBUTW, BUTH, 0, 0, 0, 0, 0, "Unlink Image block");
+				uiButSetFunc(but, image_unlink_cb, &(tex->ima), NULL);
+				
+				sprintf(str, "%d", tex->ima->id.us);
+				uiDefBut(block, BUT, B_NOP, str,
+					PANEL_XMAX-ICONBUTW, yco, ICONBUTW, BUTH, 0, 0, 0, 0, 0, "Only displays number of users of Image block");
+				uiBlockEndAlign(block);
+				
+				yco -= (BUTH);
+				
+				uiBlockBeginAlign(block);
+				uiDefButI(block, NUM, B_NOP, "Slices:",
+					X2CLM1, yco, BUTW2, BUTH, &iuser->frames, 1.0, MAXFRAMEF, 0, 0, "Number of images in the sequence");
+				uiDefButI(block, NUM, B_NOP, "Offset:",
+					X2CLM2, yco, BUTW2, BUTH, &iuser->offset, -MAXFRAMEF, MAXFRAMEF, 0, 0, "Offsets the file number to consider as the start of the sequence");
+				uiBlockEndAlign(block);
+				
+			} else {
+				but= uiDefBut(block, BUT, B_REDR, "Load",
+					X2CLM1+ICONBUTW, yco, BUTW2-ICONBUTW, BUTH, NULL, 0, 0, 0, 0, "Load new Image");
+				uiButSetFunc(but, image_load_fs_cb, &(tex->ima), iuser);
+			}
+		}
+		
+		yco -= YSPACE;
+		
+		uiDefBut(block, LABEL, B_NOP, "Interpolation:",
+				 X2CLM1, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
+		uiDefButI(block, MENU, B_REDR, "None %x0|Linear %x1|Tricubic %x2",
+				  X2CLM1, yco-=BUTH, BUTW2, BUTH, &vd->interp_type, 0.0, 0.0, 0, 0, "Interpolation type");
+		
+		yco -= YSPACE;
+		
+		uiDefBut(block, LABEL, B_NOP, "Extend:",
+				 X2CLM1, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
+		uiDefButS(block, MENU, B_REDR, "Repeat %x3|Clip %x2|Extend %x1",
+				  X2CLM1, yco-=BUTH, BUTW2, BUTH, &(tex->extend), 0.0, 0.0, 0, 0, "Extrapolation type");
+		
+		yco -= YSPACE;
+		
+		uiDefButF(block, NUM, B_REDR, "Intensity: ",
+				  X2CLM1, yco-=BUTH, BUTW2, BUTH, &(vd->int_multiplier), 0.0001, 10000.0, 0, 0, "Multiplier to scale up or down the texture's intensity");
+		
+		yco = PANEL_YMAX - 2*BUTH - YSPACE;
+		if (vd->file_format == TEX_VD_IMAGE_SEQUENCE) yco -= BUTH;
+		
+		uiDefBut(block, LABEL, B_NOP, "File Format:",
+				 X2CLM2, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
+		uiDefButS(block, MENU, B_REDR, "Blender Voxel %x0|Raw (8 bit) %x1|Image Sequence %x3",
+				  X2CLM2, yco-=BUTH, BUTW2, BUTH, &vd->file_format, 0.0, 0.0, 0, 0, "File format of the voxel data file");	
+		
+		if (ELEM(vd->file_format, TEX_VD_RAW_8BIT, TEX_VD_RAW_16BIT)) {
+			uiDefBut(block, LABEL, B_NOP, "Resolution:",
+					 X2CLM2, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
+			uiBlockBeginAlign(block);
+			uiDefButI(block, NUM, B_REDR, "X: ",
+					  X2CLM2, yco-=BUTH, BUTW2, BUTH, &(vd->resolX), 1, 10000, 0, 0, "Resolution of the voxel data");
+			uiDefButI(block, NUM, B_REDR, "Y: ",
+					  X2CLM2, yco-=BUTH, BUTW2, BUTH, &(vd->resolY), 1, 10000, 0, 0, "Resolution of the voxel data");
+			uiDefButI(block, NUM, B_REDR, "Z: ",
+					  X2CLM2, yco-= BUTH, BUTW2, BUTH, &(vd->resolZ), 1, 10000, 0, 0, "Resolution of the voxel data");
+			uiBlockEndAlign(block);
+		
+			yco -= YSPACE;
+		}	
+		
+		if (vd->file_format == TEX_VD_BLENDERVOXEL) {
+			yco -= (BUTH+YSPACE);
+			uiBlockBeginAlign(block);
+			uiDefButI(block,TOG, B_REDR, "Still",
+				  X2CLM2, yco-=BUTH, BUTW2, BUTH, &(vd->still), 0,1, 0, 0, "Use a still frame from the data sequence for the entire rendered animation");
+			if (vd->still)
+				uiDefButI(block, NUM, B_REDR, "Frame: ",
+						  X2CLM2, yco-=BUTH, BUTW2, BUTH, &(vd->still_frame), 1, 10000, 0, 0, "The frame to pause on for the entire rendered animation");	
+			uiBlockEndAlign(block);
+		}
+		
+		
+	}
+}
+
 static void texture_panel_colors(Tex *tex)
 {
 	uiBlock *block;