diff --git a/source/blender/blenlib/BLI_voxel.h b/source/blender/blenlib/BLI_voxel.h
new file mode 100644
index 0000000000000000000000000000000000000000..091d8e3682dad658434babb5316a60037980f1eb
--- /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 0000000000000000000000000000000000000000..8a1ca6e59ff50c1e219428f3045256dce3b638f1
--- /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 420850ce396ef74f5d6a85b5eda41121a2c450f7..5cbad6cf4247be3382b274aa29b486ac12edd7bb 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 8dbdde777266d7f3d6eb2ca3790b6fcb36185e20..9703ab10821b4aa36ccd702c03e5e519f88f5882 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 36c7789da80d55f80750d26457769fa129caac44..b291bdc096d9d5acb016788a5948258639335c70 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 294175460bff7c1f3ac4124a9999cd4779452605..59ab5661c19e4f9d8552a5ae92845de2389c5b27 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 b3cffb45a27ddb118d059bfba299aecb84521361..98a1a2ec2c8362f3a88a2863f5f03743343ae82d 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 01c4ac8fa74ba74ac6add06af823fd005fa15bd2..65f2916af59993fd2a49570351eb2c3b1524652a 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 f0816cfddffa95c2a9044ec3b02ba0237c72867a..241b11e04b98a56c68ceb804eeb460cb874f44f8 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 f74afde5e81be587fdd97a9e90a77f5464a5b9ef..54fc639059d6ed76a416c26dba392cf7c37777ff 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;