Merge pull request #93 from wjp/GPULink

GPULink support

46 files changed, 1328 insertions, 1370 deletions

diff --git a/astra_vc14.vcxproj b/astra_vc14.vcxproj
index a28bbfa..36d37ec 100644
--- a/astra_vc14.vcxproj
+++ b/astra_vc14.vcxproj
@@ -512,12 +512,15 @@
     <ClCompile Include="src\Float32Data.cpp" />
     <ClCompile Include="src\Float32Data2D.cpp" />
     <ClCompile Include="src\Float32Data3D.cpp" />
+    <ClCompile Include="src\Float32Data3DGPU.cpp" />
     <ClCompile Include="src\Float32Data3DMemory.cpp" />
     <ClCompile Include="src\Float32ProjectionData2D.cpp" />
     <ClCompile Include="src\Float32ProjectionData3D.cpp" />
+    <ClCompile Include="src\Float32ProjectionData3DGPU.cpp" />
     <ClCompile Include="src\Float32ProjectionData3DMemory.cpp" />
     <ClCompile Include="src\Float32VolumeData2D.cpp" />
     <ClCompile Include="src\Float32VolumeData3D.cpp" />
+    <ClCompile Include="src\Float32VolumeData3DGPU.cpp" />
     <ClCompile Include="src\Float32VolumeData3DMemory.cpp" />
     <ClCompile Include="src\ForwardProjectionAlgorithm.cpp" />
     <ClCompile Include="src\Fourier.cpp" />
@@ -624,12 +627,15 @@
     <ClInclude Include="include\astra\Float32Data.h" />
     <ClInclude Include="include\astra\Float32Data2D.h" />
     <ClInclude Include="include\astra\Float32Data3D.h" />
+    <ClInclude Include="include\astra\Float32Data3DGPU.h" />
     <ClInclude Include="include\astra\Float32Data3DMemory.h" />
     <ClInclude Include="include\astra\Float32ProjectionData2D.h" />
     <ClInclude Include="include\astra\Float32ProjectionData3D.h" />
+    <ClInclude Include="include\astra\Float32ProjectionData3DGPU.h" />
     <ClInclude Include="include\astra\Float32ProjectionData3DMemory.h" />
     <ClInclude Include="include\astra\Float32VolumeData2D.h" />
     <ClInclude Include="include\astra\Float32VolumeData3D.h" />
+    <ClInclude Include="include\astra\Float32VolumeData3DGPU.h" />
     <ClInclude Include="include\astra\Float32VolumeData3DMemory.h" />
     <ClInclude Include="include\astra\ForwardProjectionAlgorithm.h" />
     <ClInclude Include="include\astra\Fourier.h" />
diff --git a/astra_vc14.vcxproj.filters b/astra_vc14.vcxproj.filters
index dd7f574..591a4c7 100644
--- a/astra_vc14.vcxproj.filters
+++ b/astra_vc14.vcxproj.filters
@@ -321,6 +321,15 @@
     <ClCompile Include="src\CudaSirtAlgorithm3D.cpp">
       <Filter>CUDA\astra source</Filter>
     </ClCompile>
+    <ClCompile Include="src\Float32Data3DGPU.cpp">
+      <Filter>CUDA\astra source</Filter>
+    </ClCompile>
+    <ClCompile Include="src\Float32ProjectionData3DGPU.cpp">
+      <Filter>CUDA\astra source</Filter>
+    </ClCompile>
+    <ClCompile Include="src\Float32VolumeData3DGPU.cpp">
+      <Filter>CUDA\astra source</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="include\astra\Algorithm.h">
@@ -581,6 +590,15 @@
     <ClInclude Include="include\astra\CudaSirtAlgorithm3D.h">
       <Filter>CUDA\astra headers</Filter>
     </ClInclude>
+    <ClInclude Include="include\astra\Float32Data3DGPU.h">
+      <Filter>CUDA\astra headers</Filter>
+    </ClInclude>
+    <ClInclude Include="include\astra\Float32ProjectionData3DGPU.h">
+      <Filter>CUDA\astra headers</Filter>
+    </ClInclude>
+    <ClInclude Include="include\astra\Float32VolumeData3DGPU.h">
+      <Filter>CUDA\astra headers</Filter>
+    </ClInclude>
     <ClInclude Include="cuda\2d\algo.h">
       <Filter>CUDA\cuda headers</Filter>
     </ClInclude>
diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in
index 371b656..9066f0a 100644
--- a/build/linux/Makefile.in
+++ b/build/linux/Makefile.in
@@ -199,7 +199,10 @@ CUDA_CXX_OBJECTS=\
 	src/CudaFDKAlgorithm3D.lo \
 	src/CudaSirtAlgorithm3D.lo \
 	src/CudaBackProjectionAlgorithm3D.lo \
-	src/CudaForwardProjectionAlgorithm3D.lo
+	src/CudaForwardProjectionAlgorithm3D.lo \
+	src/Float32Data3DGPU.lo \
+	src/Float32ProjectionData3DGPU.lo \
+	src/Float32VolumeData3DGPU.lo
 
 CUDA_OBJECTS=\
 	cuda/2d/algo.lo \
diff --git a/build/msvc/gen.py b/build/msvc/gen.py
index 8a40c45..9c14ffe 100644
--- a/build/msvc/gen.py
+++ b/build/msvc/gen.py
@@ -274,6 +274,9 @@ P_astra["filters"]["CUDA\\astra source"] = [
 "src\\CudaSartAlgorithm.cpp",
 "src\\CudaSirtAlgorithm.cpp",
 "src\\CudaSirtAlgorithm3D.cpp",
+"src\\Float32Data3DGPU.cpp",
+"src\\Float32ProjectionData3DGPU.cpp",
+"src\\Float32VolumeData3DGPU.cpp",
 ]
 P_astra["filters"]["CUDA\\cuda headers"] = [
 "4e17872e-db7d-41bc-9760-fad1c253b583",
@@ -411,7 +414,9 @@ P_astra["filters"]["CUDA\\astra headers"] = [
 "include\\astra\\CudaSartAlgorithm.h",
 "include\\astra\\CudaSirtAlgorithm.h",
 "include\\astra\\CudaSirtAlgorithm3D.h",
-
+"include\\astra\\Float32Data3DGPU.h",
+"include\\astra\\Float32ProjectionData3DGPU.h",
+"include\\astra\\Float32VolumeData3DGPU.h",
 ]
 P_astra["filters"]["Projectors\\inline"] = [
 "0daffd63-ba49-4a5f-8d7a-5322e0e74f22",
diff --git a/cuda/3d/mem3d.cu b/cuda/3d/mem3d.cu
index 2b26fe1..97be8a4 100644
--- a/cuda/3d/mem3d.cu
+++ b/cuda/3d/mem3d.cu
@@ -118,6 +118,13 @@ MemHandle3D allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, Me
 	return ret;
 }
 
+bool zeroGPUMemory(MemHandle3D handle, unsigned int x, unsigned int y, unsigned int z)
+{
+	SMemHandle3D_internal& hnd = *handle.d.get();
+	cudaError_t err = cudaMemset3D(hnd.ptr, 0, make_cudaExtent(sizeof(float)*x, y, z));
+	return err == cudaSuccess;
+}
+
 bool freeGPUMemory(MemHandle3D handle)
 {
 	size_t free = availableGPUMemory();
@@ -307,6 +314,23 @@ bool FDK(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, co
 
 }
 
+MemHandle3D wrapHandle(float *D_ptr, unsigned int x, unsigned int y, unsigned int z, unsigned int pitch)
+{
+	cudaPitchedPtr ptr;
+	ptr.ptr = D_ptr;
+	ptr.xsize = sizeof(float) * x;
+	ptr.pitch = sizeof(float) * pitch;
+	ptr.ysize = y;
+
+	SMemHandle3D_internal h;
+	h.ptr = ptr;
+
+	MemHandle3D hnd;
+	hnd.d = boost::shared_ptr<SMemHandle3D_internal>(new SMemHandle3D_internal);
+	*hnd.d = h;
+
+	return hnd;
+}
 
 
 
diff --git a/cuda/3d/mem3d.h b/cuda/3d/mem3d.h
index a0829e2..7a87ae6 100644
--- a/cuda/3d/mem3d.h
+++ b/cuda/3d/mem3d.h
@@ -80,6 +80,8 @@ enum Mem3DZeroMode {
 size_t availableGPUMemory();
 int maxBlockDimension();
 
+MemHandle3D wrapHandle(float *D_ptr, unsigned int x, unsigned int y, unsigned int z, unsigned int pitch);
+
 MemHandle3D allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, Mem3DZeroMode zero);
 
 bool copyToGPUMemory(const float *src, MemHandle3D dst, const SSubDimensions3D &pos);
@@ -88,6 +90,8 @@ bool copyFromGPUMemory(float *dst, MemHandle3D src, const SSubDimensions3D &pos)
 
 bool freeGPUMemory(MemHandle3D handle);
 
+bool zeroGPUMemory(MemHandle3D handle, unsigned int x, unsigned int y, unsigned int z);
+
 bool setGPUIndex(int index);
 
 
@@ -97,7 +101,6 @@ bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, con
 
 bool FDK(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, bool bShortScan, const float *pfFilter = 0);
 
-
 }
 
 #endif
diff --git a/include/astra/CompositeGeometryManager.h b/include/astra/CompositeGeometryManager.h
index c0acf4f..08eb7af 100644
--- a/include/astra/CompositeGeometryManager.h
+++ b/include/astra/CompositeGeometryManager.h
@@ -42,9 +42,9 @@ namespace astra {
 
 class CCompositeVolume;
 class CCompositeProjections;
-class CFloat32Data3DMemory;
-class CFloat32ProjectionData3DMemory;
-class CFloat32VolumeData3DMemory;
+class CFloat32Data3D;
+class CFloat32ProjectionData3D;
+class CFloat32VolumeData3D;
 class CVolumeGeometry3D;
 class CProjectionGeometry3D;
 class CProjector3D;
@@ -77,7 +77,7 @@ public:
 			PART_VOL, PART_PROJ
 		} eType;
 
-		CFloat32Data3DMemory* pData;
+		CFloat32Data3D* pData;
 		unsigned int subX;
 		unsigned int subY;
 		unsigned int subZ;
@@ -88,8 +88,11 @@ public:
 		virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div) = 0;
 		virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div) = 0;
 		virtual CPart* reduce(const CPart *other) = 0;
-		virtual void getDims(size_t &x, size_t &y, size_t &z) = 0;
-		size_t getSize();
+		virtual void getDims(size_t &x, size_t &y, size_t &z) const = 0;
+		size_t getSize() const;
+
+		bool canSplitAndReduce() const;
+		bool isFull() const;
 	};
 
 	class CVolumePart : public CPart {
@@ -104,7 +107,7 @@ public:
 		virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div);
 		virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div);
 		virtual CPart* reduce(const CPart *other);
-		virtual void getDims(size_t &x, size_t &y, size_t &z);
+		virtual void getDims(size_t &x, size_t &y, size_t &z) const;
 
 		CVolumePart* clone() const;
 	};
@@ -120,7 +123,7 @@ public:
 		virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div);
 		virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div);
 		virtual CPart* reduce(const CPart *other);
-		virtual void getDims(size_t &x, size_t &y, size_t &z);
+		virtual void getDims(size_t &x, size_t &y, size_t &z) const;
 
 		CProjectionPart* clone() const;
 	};
@@ -150,23 +153,23 @@ public:
 	bool doJobs(TJobList &jobs);
 
 	SJob createJobFP(CProjector3D *pProjector,
-                     CFloat32VolumeData3DMemory *pVolData,
-                     CFloat32ProjectionData3DMemory *pProjData);
+                     CFloat32VolumeData3D *pVolData,
+                     CFloat32ProjectionData3D *pProjData);
 	SJob createJobBP(CProjector3D *pProjector,
-                     CFloat32VolumeData3DMemory *pVolData,
-                     CFloat32ProjectionData3DMemory *pProjData);
+                     CFloat32VolumeData3D *pVolData,
+                     CFloat32ProjectionData3D *pProjData);
 
 	// Convenience functions for creating and running a single FP or BP job
-	bool doFP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
-	          CFloat32ProjectionData3DMemory *pProjData);
-	bool doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
-	          CFloat32ProjectionData3DMemory *pProjData);
-	bool doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
-	          CFloat32ProjectionData3DMemory *pProjData, bool bShortScan,
+	bool doFP(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData,
+	          CFloat32ProjectionData3D *pProjData);
+	bool doBP(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData,
+	          CFloat32ProjectionData3D *pProjData);
+	bool doFDK(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData,
+	          CFloat32ProjectionData3D *pProjData, bool bShortScan,
 	          const float *pfFilter = 0);
 
-	bool doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData);
-	bool doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData);
+	bool doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3D *>& volData, const std::vector<CFloat32ProjectionData3D *>& projData);
+	bool doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3D *>& volData, const std::vector<CFloat32ProjectionData3D *>& projData);
 
 	void setGPUIndices(const std::vector<int>& GPUIndices);
 
diff --git a/include/astra/CudaBackProjectionAlgorithm3D.h b/include/astra/CudaBackProjectionAlgorithm3D.h
index 6738988..114d6f3 100644
--- a/include/astra/CudaBackProjectionAlgorithm3D.h
+++ b/include/astra/CudaBackProjectionAlgorithm3D.h
@@ -69,8 +69,8 @@ public:
 	 * @param _pReconstruction	VolumeData3D object for storing the reconstructed volume.
 	 */
 	CCudaBackProjectionAlgorithm3D(CProjector3D* _pProjector, 
-	                     CFloat32ProjectionData3DMemory* _pProjectionData, 
-	                     CFloat32VolumeData3DMemory* _pReconstruction);
+	                     CFloat32ProjectionData3D* _pProjectionData, 
+	                     CFloat32VolumeData3D* _pReconstruction);
 	
 	/** Copy constructor.
 	 */
@@ -99,8 +99,8 @@ public:
 	 * @return initialization successful?
 	 */
 	bool initialize(CProjector3D* _pProjector, 
-					CFloat32ProjectionData3DMemory* _pSinogram, 
-					CFloat32VolumeData3DMemory* _pReconstruction);
+					CFloat32ProjectionData3D* _pSinogram, 
+					CFloat32VolumeData3D* _pReconstruction);
 
 	/** Get all information parameters
 	 *
diff --git a/include/astra/CudaFDKAlgorithm3D.h b/include/astra/CudaFDKAlgorithm3D.h
index 386129e..1c4c622 100644
--- a/include/astra/CudaFDKAlgorithm3D.h
+++ b/include/astra/CudaFDKAlgorithm3D.h
@@ -81,8 +81,8 @@ public:
 	 * @param _pReconstruction	VolumeData3D object for storing the reconstructed volume.
 	 */
 	CCudaFDKAlgorithm3D(CProjector3D* _pProjector, 
-	                    CFloat32ProjectionData3DMemory* _pProjectionData, 
-	                    CFloat32VolumeData3DMemory* _pReconstruction);
+	                    CFloat32ProjectionData3D* _pProjectionData, 
+	                    CFloat32VolumeData3D* _pReconstruction);
 	
 	/** Copy constructor.
 	 */
@@ -111,8 +111,8 @@ public:
 	 * @return initialization successful?
 	 */
 	bool initialize(CProjector3D* _pProjector, 
-					CFloat32ProjectionData3DMemory* _pSinogram, 
-					CFloat32VolumeData3DMemory* _pReconstruction);
+					CFloat32ProjectionData3D* _pSinogram, 
+					CFloat32VolumeData3D* _pReconstruction);
 
 	/** Get all information parameters
 	 *
diff --git a/include/astra/CudaForwardProjectionAlgorithm3D.h b/include/astra/CudaForwardProjectionAlgorithm3D.h
index 95af73a..9dc889e 100644
--- a/include/astra/CudaForwardProjectionAlgorithm3D.h
+++ b/include/astra/CudaForwardProjectionAlgorithm3D.h
@@ -71,8 +71,8 @@ public:
 	 * @return initialization successful?
 	 */
 	bool initialize(CProjector3D* _pProjector, 
-					CFloat32ProjectionData3DMemory* _pSinogram, 
-					CFloat32VolumeData3DMemory* _pReconstruction,
+					CFloat32ProjectionData3D* _pSinogram, 
+					CFloat32VolumeData3D* _pReconstruction,
 					int _iGPUindex = -1, int _iDetectorSuperSampling = 1);
 
 
@@ -116,8 +116,8 @@ public:
 
 protected:
 	CProjector3D* m_pProjector;
-	CFloat32ProjectionData3DMemory* m_pProjections;
-	CFloat32VolumeData3DMemory* m_pVolume;
+	CFloat32ProjectionData3D* m_pProjections;
+	CFloat32VolumeData3D* m_pVolume;
 	int m_iGPUIndex;
 	int m_iDetectorSuperSampling;
 
diff --git a/include/astra/Float32Data2D.h b/include/astra/Float32Data2D.h
index b33bcd0..1524683 100644
--- a/include/astra/Float32Data2D.h
+++ b/include/astra/Float32Data2D.h
@@ -41,7 +41,7 @@ public:
 };
 
 /** 
- * This class represents a two-dimensional block of float32ing point data.
+ * This class represents a 2-dimensional block of 32-bit floating point data.
  * It contains member functions for accessing this data and for performing 
  * elementary computations on the data.
  * The data block is "owned" by the class, meaning that the class is 
diff --git a/include/astra/Float32Data3D.h b/include/astra/Float32Data3D.h
index 1cd9c20..226689e 100644
--- a/include/astra/Float32Data3D.h
+++ b/include/astra/Float32Data3D.h
@@ -32,10 +32,12 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "Float32Data.h"
 #include "Float32Data2D.h"
 
+#include "../../cuda/3d/mem3d.h"
+
 namespace astra {
 
 /**
- * This class represents a three-dimensional block of float32ing point data.
+ * This class represents a 3-dimensional block of 32-bit floating point data.
  */
 class _AstraExport CFloat32Data3D : public CFloat32Data {
 
@@ -107,22 +109,6 @@ public:
 	 */
 	int getDimensionCount() const;	
 
-	/**
-	 * Clamp data to minimum value
-	 *
-	 * @param _fMin minimum value
-	 * @return l-value
-	 */
-	virtual CFloat32Data3D& clampMin(float32& _fMin) = 0;
-
-	/**
-	 * Clamp data to maximum value
-	 *
-	 * @param _fMax maximum value
-	 * @return l-value
-	 */
-	virtual CFloat32Data3D& clampMax(float32& _fMax) = 0;
-
 	/** get a description of the class
 	 *
 	 * @return description string
diff --git a/include/astra/Float32Data3DGPU.h b/include/astra/Float32Data3DGPU.h
new file mode 100644
index 0000000..1247e65
--- /dev/null
+++ b/include/astra/Float32Data3DGPU.h
@@ -0,0 +1,107 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp
+           2014-2016, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://www.astra-toolbox.com/
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox 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 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+*/
+
+#ifndef _INC_ASTRA_FLOAT32DATA3DGPU
+#define _INC_ASTRA_FLOAT32DATA3DGPU
+
+#ifdef ASTRA_CUDA
+
+#include "Globals.h"
+#include "Float32Data3D.h"
+
+#include "../../cuda/3d/mem3d.h"
+
+namespace astra {
+
+
+astraCUDA3d::MemHandle3D wrapHandle(float *D_ptr, unsigned int x, unsigned int y, unsigned int z, unsigned int pitch);
+
+
+/** 
+ * This class represents a 3-dimensional block of 32-bit floating point data.
+ * The data block is stored on a GPU, and owned by external code.
+ *
+ * TODO: Store/remember which GPU the data is stored on
+ */
+class _AstraExport CFloat32Data3DGPU : public virtual CFloat32Data3D {
+
+protected:
+	/** Handle for the memory block */
+	astraCUDA3d::MemHandle3D m_hnd;
+
+	/** Clear all member variables, setting all numeric variables to 0 and all pointers to NULL. 
+	 */
+	void _clear();
+
+	/** Un-initialize the object, bringing it back in the unitialized state.
+	 */
+	void _unInit();
+
+	/** Initialization. Initializes an instance of the CFloat32Data3DGPU class.
+	 * Can only be called by derived classes.
+	 *
+	 * This function does not set m_bInitialized to true if everything is ok.
+	 *
+	 * @param _iWidth width of the 3D data (x-axis), must be > 0
+	 * @param _iHeight height of the 3D data (y-axis), must be > 0
+	 * @param _iDepth depth of the 3D data (z-axis), must be > 0
+	 * @param _hnd the CUDA memory handle
+	 */
+
+	bool _initialize(int _iWidth, int _iHeight, int _iDepth, astraCUDA3d::MemHandle3D _hnd);
+
+public:
+
+	/** Default constructor. Sets all numeric member variables to 0 and all pointer member variables to NULL.
+	 *
+	 * If an object is constructed using this default constructor, it must always be followed by a call 
+	 * to one of the initialize() methods before the object can be used. Any use before calling init() is not allowed,
+	 * except calling the member function isInitialized().
+	 *
+	 */
+	CFloat32Data3DGPU();
+
+	/** Destructor.
+	 */
+	virtual ~CFloat32Data3DGPU();
+
+	/** which type is this class?
+	 *
+	 * @return DataType: ASTRA_DATATYPE_FLOAT32_PROJECTION or
+	 *					 ASTRA_DATATYPE_FLOAT32_VOLUME
+	 */
+	virtual EDataType getType() const { return BASE; }
+
+	astraCUDA3d::MemHandle3D getHandle() const { return m_hnd; }
+
+};
+
+} // end namespace astra
+
+#endif
+
+#endif // _INC_ASTRA_FLOAT32DATA3DGPU
diff --git a/include/astra/Float32Data3DMemory.h b/include/astra/Float32Data3DMemory.h
index d885101..876aa37 100644
--- a/include/astra/Float32Data3DMemory.h
+++ b/include/astra/Float32Data3DMemory.h
@@ -34,7 +34,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 namespace astra {
 
 /** 
- * This class represents a three-dimensional block of float32ing point data.
+ * This class represents a 3-dimensional block of 32-bit floating point data.
  * It contains member functions for accessing this data and for performing 
  * elementary computations on the data.
  * The data block is "owned" by the class, meaning that the class is 
@@ -106,8 +106,8 @@ protected:
 	 * This function does not set m_bInitialized to true if everything is ok.
 	 *
 	 * @param _iWidth width of the 3D data (x-axis), must be > 0
-	 * @param _iHeight height of the 3D data (y-axis), must be > 0 
-	 * @param _iDepth depth of the 3D data (z-axis), must be > 0 
+	 * @param _iHeight height of the 3D data (y-axis), must be > 0
+	 * @param _iDepth depth of the 3D data (z-axis), must be > 0
 	 * @return initialization of the base class successfull
 	 */
 	bool _initialize(int _iWidth, int _iHeight, int _iDepth);
@@ -121,9 +121,9 @@ protected:
 	 * object is reinitialized and memory is freed and reallocated if necessary.
 	 * This function does not set m_bInitialized to true if everything is ok.
 	 *
-	 * @param _iWidth width of the 2D data (x-axis), must be > 0
-	 * @param _iHeight height of the 2D data (y-axis), must be > 0 
-	 * @param _iDepth depth of the 2D data (z-axis), must be > 0 
+	 * @param _iWidth width of the 3D data (x-axis), must be > 0
+	 * @param _iHeight height of the 3D data (y-axis), must be > 0
+	 * @param _iDepth depth of the 3D data (z-axis), must be > 0
 	 * @param _pfData pointer to a one-dimensional float32 data block
 	 * @return initialization of the base class successfull
 	 */
@@ -138,9 +138,9 @@ protected:
 	 * object is reinitialized and memory is freed and reallocated if necessary.
 	 * This function does not set m_bInitialized to true if everything is ok.
 	 *
-	 * @param _iWidth width of the 2D data (x-axis), must be > 0
-	 * @param _iHeight height of the 2D data (y-axis), must be > 0 
-	 * @param _iDepth depth of the 2D data (z-axis), must be > 0 
+	 * @param _iWidth width of the 3D data (x-axis), must be > 0
+	 * @param _iHeight height of the 3D data (y-axis), must be > 0
+	 * @param _iDepth depth of the 3D data (z-axis), must be > 0
 	 * @param _fScalar scalar value to fill the data
 	 * @return initialization of the base class successfull
 	 */
@@ -157,9 +157,9 @@ protected:
 	 * object is reinitialized and memory is freed and reallocated if necessary.
 	 * This function does not set m_bInitialized to true if everything is ok.
 	 *
-	 * @param _iWidth width of the 2D data (x-axis), must be > 0
-	 * @param _iHeight height of the 2D data (y-axis), must be > 0 
-	 * @param _iDepth depth of the 2D data (z-axis), must be > 0 
+	 * @param _iWidth width of the 3D data (x-axis), must be > 0
+	 * @param _iHeight height of the 3D data (y-axis), must be > 0
+	 * @param _iDepth depth of the 3D data (z-axis), must be > 0
 	 * @param _pCustomMemory the custom memory handle
 	 */
 
@@ -237,35 +237,11 @@ public:
 	 * After the call p = getData3D(), use p[iy][ix] to access element (ix, iy, iz).
 	 * The data memory and pointer array are still "owned" by the CFloat32Data3DMemory 
 	 * instance; this memory may NEVER be freed by the caller of this function. 
-	 * If changes are made to this data, the function updateStatistics() 
-	 * should be called after completion of all changes. 
 	 *
 	 * @return pointer to the 3-dimensional 32-bit floating point data block
  	 */
 	const float32*** getData3DConst() const;
 
-	/** Update data statistics, such as minimum and maximum value, after the data has been modified. 
-	 */
-	virtual void updateStatistics();
-
-	/** Get the minimum value in the data block.
-	 * If the data has been changed after construction, the function
-	 * updateStatistics() must be called at least once before 
-	 * a query can be made on this value.
-	 *
-	 * @return minimum value in the data block
-	 */
-	virtual float32 getGlobalMin() const;
-
-	/** Get the maximum value in the data block
-	 * If the data has been changed after construction, the function
-	 * updateStatistics() must be called at least once before 
-	 * a query can be made on this value.
-	 *
-	 * @return maximum value in the data block
-	 */
-	virtual float32 getGlobalMax() const;
-
 	/** which type is this class?
 	 *
 	 * @return DataType: ASTRA_DATATYPE_FLOAT32_PROJECTION or
@@ -306,26 +282,11 @@ inline CFloat32Data3DMemory::EDataType CFloat32Data3DMemory::getType() const
 }
 
 //----------------------------------------------------------------------------------------
-// Get the minimum value in the data block.
-inline float32 CFloat32Data3DMemory::getGlobalMin() const
-{
-	ASTRA_ASSERT(m_bInitialized);
-	return m_fGlobalMin;
-}
-
-//----------------------------------------------------------------------------------------
-// Get the maximum value in the data block
-inline float32 CFloat32Data3DMemory::getGlobalMax() const
-{
-	ASTRA_ASSERT(m_bInitialized);
-	return m_fGlobalMax;
-}
-
-//----------------------------------------------------------------------------------------
 // Get a pointer to the data block, represented as a 1-dimensional array of float32 values.
 inline float32* CFloat32Data3DMemory::getData()
 {
 	ASTRA_ASSERT(m_bInitialized);
+
 	return m_pfData;
 }
 
@@ -334,6 +295,7 @@ inline float32* CFloat32Data3DMemory::getData()
 inline const float32* CFloat32Data3DMemory::getDataConst() const
 {
 	ASTRA_ASSERT(m_bInitialized);
+
 	return (const float32*)m_pfData;
 }
 
diff --git a/include/astra/Float32ProjectionData3D.h b/include/astra/Float32ProjectionData3D.h
index ae0664b..1634eeb 100644
--- a/include/astra/Float32ProjectionData3D.h
+++ b/include/astra/Float32ProjectionData3D.h
@@ -85,55 +85,6 @@ public:
 	 */
 	virtual CFloat32Data3D::EDataType getType() const;
 
-	/** Fetch a COPY of a projection of the data.  Note that if you update the 2D data slice, the data in the 
-	 * 3d data object will remain unaltered.  To copy the data back in the 3D-volume you must return the data by calling 'returnProjection'.
-	 *
-	 * @param _iProjectionNr projection number
-	 * @return Volume data object
-	 */
-	virtual CFloat32VolumeData2D* fetchProjection(int _iProjectionNr) const = 0;
-	
-	/** Return a projection slice to the 3d data.  The data will be deleted. If the slice was fetched with 
-	 * 'fetchProjection', the data will be stored first. 
-	 *
-	 * @param _iProjectionNr projection number
-	 * @param _pProjection 2D Projection Data
-	 */
-	virtual void returnProjection(int _iProjectionNr, CFloat32VolumeData2D* _pProjection) = 0;
-
-	/** Fetch a COPY of a sinogram slice of the data.  Note that if you update the 2D data slice, the data in the 
-	 * 3d data object will remain unaltered.  To copy the data back in the 3D-volume you must return the data by calling 'returnSlice'.
-	 *
-	 * @param _iSliceNr slice number
-	 * @return Sinogram data object
-	 */
-	virtual CFloat32ProjectionData2D* fetchSinogram(int _iSliceNr) const = 0;
-
-	/** Return a sinogram slice to the 3d data.  The data will be stored in the 3D Data object.
-	 *
-	 * @param _iSliceNr slice number
-	 * @param _pSinogram2D 2D Sinogram Object.
-	 */
-	virtual void returnSinogram(int _iSliceNr, CFloat32ProjectionData2D* _pSinogram2D) = 0;
-
-	/** This SLOW function returns a detector value stored a specific index in the array.
-	 *  Reading values in this way might cause a lot of unnecessar__y memory operations, don't
-	 *  use it in time-critical code.
-	 * 
-	 *  @param _iIndex Index in the array if the data were stored completely in main memory
-	 *  @return The value the location specified by _iIndex
-	 */
-	virtual float32 getDetectorValue(int _iIndex) = 0;
-
-	/** This SLOW function stores a detector value at a specific index in the array.
-	 *  Writing values in this way might cause a lot of unnecessary memory operations, don't
-	 *  use it in time-critical code.
-	 * 
-	 *  @param _iIndex Index in the array if the data were stored completely in main memory
-	 *  @param _fValue The value to be stored at the location specified by _iIndex
-	 */
-	virtual void setDetectorValue(int _iIndex, float32 _fValue) = 0;
-
 	/**
 	 * Overloaded Operator: data += data (pointwise)
 	 *
diff --git a/include/astra/Float32ProjectionData3DGPU.h b/include/astra/Float32ProjectionData3DGPU.h
new file mode 100644
index 0000000..cbf9b6f
--- /dev/null
+++ b/include/astra/Float32ProjectionData3DGPU.h
@@ -0,0 +1,91 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp
+           2014-2016, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://www.astra-toolbox.com/
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox 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 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+*/
+
+#ifndef _INC_ASTRA_FLOAT32PROJECTIONDATA3DGPU
+#define _INC_ASTRA_FLOAT32PROJECTIONDATA3DGPU
+
+#include "Float32Data3DGPU.h"
+#include "ProjectionGeometry3D.h"
+#include "Float32ProjectionData3D.h"
+
+#ifdef ASTRA_CUDA
+
+namespace astra {
+
+/**
+ * This class represents three-dimensional Projection Data where the entire data block is stored in GPU memory.
+ */
+class _AstraExport CFloat32ProjectionData3DGPU : public CFloat32Data3DGPU, public CFloat32ProjectionData3D
+{
+public:
+
+	/** Default constructor. Sets all numeric member variables to 0 and all pointer member variables to NULL.
+	 *
+	 * If an object is constructed using this default constructor, it must always be followed by a call 
+	 * to one of the init() methods before the object can be used. Any use before calling init() is not allowed,
+	 * except calling the member function isInitialized().
+	 *
+	 */
+	CFloat32ProjectionData3DGPU();
+	
+	/** Construction.
+	 *
+	 * @param _pGeometry 3D volume geometry
+	 * @param _hnd the CUDA memory handle
+	 */
+
+	CFloat32ProjectionData3DGPU(CProjectionGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd);
+
+	virtual ~CFloat32ProjectionData3DGPU();
+
+	/** Initialization.
+	 *
+	 * @param _pGeometry 3D volume geometry
+	 * @param _hnd the CUDA memory handle
+	 */
+
+	bool initialize(CProjectionGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd);
+
+	/** Which type is this class?
+	 *
+	 * @return DataType: PROJECTION
+	 */
+	virtual CFloat32Data3D::EDataType getType() const { return PROJECTION; }
+
+	/** Get the volume geometry.
+	 *
+	 * @return pointer to volume geometry.
+	 */
+	CProjectionGeometry3D* getGeometry() const { ASTRA_ASSERT(m_bInitialized); return m_pGeometry; }
+
+};
+
+} // end namesProjection astra
+
+#endif
+
+#endif // _INC_ASTRA_FLOAT32PROJECTIONDATA3DGPU
diff --git a/include/astra/Float32ProjectionData3DMemory.h b/include/astra/Float32ProjectionData3DMemory.h
index 14bae36..732c31a 100644
--- a/include/astra/Float32ProjectionData3DMemory.h
+++ b/include/astra/Float32ProjectionData3DMemory.h
@@ -166,55 +166,6 @@ public:
 
 
 
-	/** Fetch a COPY of a projection of the data.  Note that if you update the 2D data slice, the data in the 
-	 * 3D data object will remain unaltered.  To copy the data back in the 3D-volume you must return the data by calling 'returnProjection'.
-	 *
-	 * @param _iProjectionNr projection number
-	 * @return Volume data object
-	 */
-	virtual CFloat32VolumeData2D* fetchProjection(int _iProjectionNr) const;
-	
-	/** Return a projection slice to the 3D data.  The data will be deleted. If the slice was fetched with 
-	 * 'fetchProjection', the data will be stored first. 
-	 *
-	 * @param _iProjectionNr projection number
-	 * @param _pProjection 2D Projection image
-	 */
-	virtual void returnProjection(int _iProjectionNr, CFloat32VolumeData2D* _pProjection);
-
-	/** Fetch a COPY of a sinogram slice of the data.  Note that if you update the 2D data slice, the data in the 
-	 * 3D data object will remain unaltered.  To copy the data back in the 3D-volume you must return the data by calling 'returnSlice'.
-	 *
-	 * @param _iSliceNr slice number
-	 * @return Sinogram data object
-	 */
-	virtual CFloat32ProjectionData2D* fetchSinogram(int _iSliceNr) const;
-
-	/** This SLOW function returns a detector value stored a specific index in the array.
-	 *  Reading values in this way might cause a lot of unnecessary memory operations, don't
-	 *  use it in time-critical code.
-	 * 
-	 *  @param _iIndex Index in the array if the data were stored completely in main memory
-	 *  @return The value the location specified by _iIndex
-	 */
-	virtual float32 getDetectorValue(int _iIndex);
-
-	/** This SLOW function stores a detector value at a specific index in the array.
-	 *  Writing values in this way might cause a lot of unnecessary memory operations, don't
-	 *  use it in time-critical code.
-	 * 
-	 *  @param _iIndex Index in the array if the data were stored completely in main memory
-	 *  @param _fValue The value to be stored at the location specified by _iIndex
-	 */
-	virtual void setDetectorValue(int _iIndex, float32 _fValue);
-
-	/** Return a sinogram slice to the 3d data.  The data will be stored in the 3D Data object.
-	 *
-	 * @param _iSliceNr slice number
-	 * @param _pSinogram2D 2D Sinogram Object.
-	 */
-	virtual void returnSinogram(int _iSliceNr, CFloat32ProjectionData2D* _pSinogram2D);
-
 	/** Which type is this class?
 	 *
 	 * @return DataType: PROJECTION 
diff --git a/include/astra/Float32VolumeData3D.h b/include/astra/Float32VolumeData3D.h
index 7a32efe..535e960 100644
--- a/include/astra/Float32VolumeData3D.h
+++ b/include/astra/Float32VolumeData3D.h
@@ -85,72 +85,6 @@ public:
 	 */
 	virtual CFloat32Data3D::EDataType getType() const;
 	
-	/** Fetch a slice from the data in the x direction.  Note that if you update the 2D data slice, the data in the 
-	 * 3d data object will remain unaltered.  To copy the data you must return the data by calling 'returnSliceX'.
-	 * You should not delete data fetched with this function yourself, instead call the 'returnSliceX' function.
-	 *
-	 * @param _iColumnIndex slice number
-	 * @return Volume data object
-	 */
-	virtual CFloat32VolumeData2D* fetchSliceX(int _iColumnIndex) const = 0;
-
-	/** Fetch a slice from the data in the y direction.  Note that if you update the 2D data slice, the data in the 
-	 * 3d data object will remain unaltered.  To copy the data you must return the data by calling 'returnSliceY'.
-	 * You should not delete data fetched with this function yourself, instead call the 'returnSliceY' function.
-	 *
-	 * @param _iRowIndex slice number
-	 * @return Volume data object
-	 */
-	virtual CFloat32VolumeData2D* fetchSliceY(int _iRowIndex) const = 0;
-
-	/** Fetch a slice from the data in the z direction.  Note that if you update the 2D data slice, the data in the 
-	 * 3d data object will remain unaltered.  To copy the data you must return the data by calling 'returnSliceZ'.
-	 * You should not delete data fetched with this function yourself, instead call the 'returnSliceZ' function.
-	 *
-	 * @param _iSliceIndex slice number
-	 * @return Volume data object
-	 */
-	virtual CFloat32VolumeData2D* fetchSliceZ(int _iSliceIndex) const = 0;
-
-	/** Return a slice from the data in the x direction to the 3d data.  The data will be deleted. If the slice was 
-	 * fetched with 'fetchSliceX', the data will be stored first. 
-	 *
-	 * @param _iColumnIndex slice number
-	 */
-	virtual void returnSliceX(int _iColumnIndex, CFloat32VolumeData2D * _pSlice) = 0;
-	
-	/** Return a slice from the data in the y direction to the 3d data.  The data will be deleted. If the slice was 
-	 * fetched with 'fetchSliceY', the data will be stored first. 
-	 *
-	 * @param _iRowIndex slice number
-	 */
-	virtual void returnSliceY(int _iRowIndex, CFloat32VolumeData2D * _pSlice) = 0;
-
-	/** Return a slice from the data in the z direction to the 3d data.  The data will be deleted. If the slice was 
-	 * fetched with 'fetchSliceZ', the data will be stored first. 
-	 *
-	 * @param _iSliceIndex slice number
-	 */
-	virtual void returnSliceZ(int _iSliceIndex, CFloat32VolumeData2D * _pSlice) = 0;
-
-	/** This SLOW function returns a voxel value stored at a specific index in the array.
-	 *  Reading values in this way might cause a lot of unnecessary memory operations, don't
-	 *  use it in time-critical code.
-	 * 
-	 *  @param _iIndex Index in the array if the data were stored completely in main memory
-	 *  @return The value stored at the location specified by _iIndex
-	 */
-	virtual float32 getVoxelValue(int _iIndex) = 0;
-
-	/** This SLOW function stores a voxel value at a specific index in the array.
-	 *  Writing values in this way might cause a lot of unnecessary memory operations, don't
-	 *  use it in time-critical code.
-	 * 
-	 *  @param _iIndex Index in the array if the data were stored completely in main memory
-	 *  @param _fValue The value to be stored at the location specified by _iIndex
-	 */
-	virtual void setVoxelValue(int _iIndex, float32 _fValue) = 0;
-
 	/**
 	 * Overloaded Operator: data += data (pointwise)
 	 *
diff --git a/include/astra/Float32VolumeData3DGPU.h b/include/astra/Float32VolumeData3DGPU.h
new file mode 100644
index 0000000..ddeb469
--- /dev/null
+++ b/include/astra/Float32VolumeData3DGPU.h
@@ -0,0 +1,91 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp
+           2014-2016, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://www.astra-toolbox.com/
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox 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 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+*/
+
+#ifndef _INC_ASTRA_FLOAT32VOLUMEDATA3DGPU
+#define _INC_ASTRA_FLOAT32VOLUMEDATA3DGPU
+
+#include "Float32Data3DGPU.h"
+#include "VolumeGeometry3D.h"
+#include "Float32VolumeData3D.h"
+
+#ifdef ASTRA_CUDA
+
+namespace astra {
+
+/**
+ * This class represents three-dimensional Volume Data where the entire data block is stored in GPU memory.
+ */
+class _AstraExport CFloat32VolumeData3DGPU : public CFloat32Data3DGPU, public CFloat32VolumeData3D
+{
+public:
+
+	/** Default constructor. Sets all numeric member variables to 0 and all pointer member variables to NULL.
+	 *
+	 * If an object is constructed using this default constructor, it must always be followed by a call 
+	 * to one of the init() methods before the object can be used. Any use before calling init() is not allowed,
+	 * except calling the member function isInitialized().
+	 *
+	 */
+	CFloat32VolumeData3DGPU();
+	
+	/** Construction.
+	 *
+	 * @param _pGeometry 3D volume geometry
+	 * @param _hnd the CUDA memory handle
+	 */
+
+	CFloat32VolumeData3DGPU(CVolumeGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd);
+
+	virtual ~CFloat32VolumeData3DGPU();
+
+	/** Initialization.
+	 *
+	 * @param _pGeometry 3D volume geometry
+	 * @param _hnd the CUDA memory handle
+	 */
+
+	bool initialize(CVolumeGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd);
+
+	/** Which type is this class?
+	 *
+	 * @return DataType: VOLUME
+	 */
+	virtual CFloat32Data3D::EDataType getType() const { return VOLUME; }
+
+	/** Get the volume geometry.
+	 *
+	 * @return pointer to volume geometry.
+	 */
+	CVolumeGeometry3D* getGeometry() const { ASTRA_ASSERT(m_bInitialized); return m_pGeometry; }
+
+};
+
+} // end namespace astra
+
+#endif
+
+#endif // _INC_ASTRA_FLOAT32VOLUMEDATA3DGPU
diff --git a/include/astra/Float32VolumeData3DMemory.h b/include/astra/Float32VolumeData3DMemory.h
index 70f0939..4a70f65 100644
--- a/include/astra/Float32VolumeData3DMemory.h
+++ b/include/astra/Float32VolumeData3DMemory.h
@@ -159,55 +159,6 @@ public:
 	CVolumeGeometry3D* getGeometry() const;
 
 	/**
-	 * Gets a slice, containing all voxels with a given x (= column) index.
-	 */
-	CFloat32VolumeData2D * fetchSliceX(int _iColumnIndex) const;
-
-	/**
-	 * Gets a slice, containing all voxels with a given y (= row) index.
-	 */
-	CFloat32VolumeData2D * fetchSliceY(int _iRowIndex) const;
-
-	/**
-	 * Gets a slice, containing all voxels with a given z (= slice) index.
-	 */
-	CFloat32VolumeData2D * fetchSliceZ(int _iSliceIndex) const;
-
-	/**
-	 * Gets a slice, containing all voxels with a given x (= column) index.
-	 */
-	void returnSliceX(int _iColumnIndex, CFloat32VolumeData2D * _pSliceData);
-
-	/**
-	 * Gets a slice, containing all voxels with a given y (= row) index.
-	 */
-	void returnSliceY(int _iRowIndex, CFloat32VolumeData2D * _pSliceData);
-
-	/**
-	 * Copies data from a 2D slice containing all voxels with a given z (= slice) index to the
-	 * 3D  memory stored in this class.
-	 */
-	void returnSliceZ(int _iSliceIndex, CFloat32VolumeData2D * _pSliceData);
-
-	/** This SLOW function returns a volume value stored a specific index in the array.
-	 *  Reading values in this way might cause a lot of unnecessary memory operations, don't
-	 *  use it in time-critical code.
-	 * 
-	 *  @param _iIndex Index in the array if the data were stored completely in main memory
-	 *  @return The value the location specified by _iIndex
-	 */
-	virtual float32 getVoxelValue(int _iIndex);
-
-	/** This SLOW function stores a voxel value at a specific index in the array.
-	 *  Writing values in this way might cause a lot of unnecessary memory operations, don't
-	 *  use it in time-critical code.
-	 * 
-	 *  @param _iIndex Index in the array if the data were stored completely in main memory
-	 *  @param _fValue The value to be stored at the location specified by _iIndex
-	 */
-	virtual void setVoxelValue(int _iIndex, float32 _fValue);
-
-	/**
 	 * Overloaded Operator: data = data (pointwise)
 	 *
 	 * @param _dataIn r-value
diff --git a/include/astra/Globals.h b/include/astra/Globals.h
index dec978d..44a77b0 100644
--- a/include/astra/Globals.h
+++ b/include/astra/Globals.h
@@ -53,8 +53,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include <iostream>
 #include <fstream>
 #include <math.h>
-#include <boost/static_assert.hpp>
-#include <boost/throw_exception.hpp>
 
 //----------------------------------------------------------------------------------------
 // macro's
@@ -227,16 +225,6 @@ namespace astra {
 		int m_iSliceIndex;
 	};
 }
-//----------------------------------------------------------------------------------------
-// some toys
-
-// safe reinterpret cast
-template <class To, class From>
-To safe_reinterpret_cast(From from)
-{
-	BOOST_STATIC_ASSERT(sizeof(From) <= sizeof(To));
-	return reinterpret_cast<To>(from);
-}
 
 //----------------------------------------------------------------------------------------
 // functions for testing
diff --git a/matlab/mex/astra_mex_data3d_c.cpp b/matlab/mex/astra_mex_data3d_c.cpp
index 7909a79..f0aa3c8 100644
--- a/matlab/mex/astra_mex_data3d_c.cpp
+++ b/matlab/mex/astra_mex_data3d_c.cpp
@@ -60,7 +60,7 @@ using namespace astra;
 //-----------------------------------------------------------------------------------------
 /**
  * id = astra_mex_io_data('create', datatype, geometry, data);
- *        datatype: ['-vol','-sino','-sinocone'] 
+ *        datatype: ['-vol','-sino']
  */
 void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray* prhs[])
 { 
@@ -103,7 +103,6 @@ void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra
 		copyMexToCFloat32Array(data, pDataObject3D->getData(),
 				pDataObject3D->getSize());
 	}
-	pDataObject3D->updateStatistics();
 
 	// step4: store data object
 	int iIndex = CData3DManager::getSingleton().store(pDataObject3D);
@@ -162,8 +161,6 @@ void astra_mex_data3d_link(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray*
 		return;
 	}
 
-	//pDataObject3D->updateStatistics();
-
 	// step4: store data object
 	int iIndex = CData3DManager::getSingleton().store(pDataObject3D);
 
@@ -234,7 +231,6 @@ void astra_mex_data3d_store(int nlhs, mxArray* plhs[], int nrhs, const mxArray*
 	}
 
 	copyMexToCFloat32Array(prhs[2], pDataObject->getData(), pDataObject->getSize());
-	pDataObject->updateStatistics();
 }
 
 void astra_mex_data3d_dimensions(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
diff --git a/matlab/mex/mexDataManagerHelpFunctions.cpp b/matlab/mex/mexDataManagerHelpFunctions.cpp
index ebf6f8f..c47ee27 100644
--- a/matlab/mex/mexDataManagerHelpFunctions.cpp
+++ b/matlab/mex/mexDataManagerHelpFunctions.cpp
@@ -166,17 +166,6 @@ checkDataSize(const mxArray * const mArray,
 
 //-----------------------------------------------------------------------------------------
 void
-updateStatistics(const std::vector<astra::CFloat32Data3DMemory *> & vecIn)
-{
-	const size_t tot_size = vecIn.size();
-	for (size_t count = 0; count < tot_size; count++)
-	{
-		vecIn[count]->updateStatistics();
-	}
-}
-
-//-----------------------------------------------------------------------------------------
-void
 getDataPointers(const std::vector<astra::CFloat32Data3DMemory *> & vecIn,
 		std::vector<astra::float32 *> & vecOut)
 {
diff --git a/matlab/mex/mexDataManagerHelpFunctions.h b/matlab/mex/mexDataManagerHelpFunctions.h
index fff83a0..b7edb0f 100644
--- a/matlab/mex/mexDataManagerHelpFunctions.h
+++ b/matlab/mex/mexDataManagerHelpFunctions.h
@@ -50,8 +50,6 @@ bool checkDataSize(const mxArray * const, const astra::CProjectionGeometry3D * c
 bool checkDataSize(const mxArray * const, const astra::CVolumeGeometry3D * const,
 		const mwIndex & zOffset);
 
-void updateStatistics(const std::vector<astra::CFloat32Data3DMemory *> &);
-
 void getDataPointers(const std::vector<astra::CFloat32Data3DMemory *> &,
 		std::vector<astra::float32 *> &);
 void getDataSizes(const std::vector<astra::CFloat32Data3DMemory *> &,
diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd
index 39f9039..b40d787 100644
--- a/python/astra/PyIncludes.pxd
+++ b/python/astra/PyIncludes.pxd
@@ -30,242 +30,262 @@ from .PyXMLDocument cimport XMLNode
 include "config.pxi"
 
 cdef extern from "astra/Globals.h" namespace "astra":
-	ctypedef float float32
-	ctypedef double float64
-	ctypedef unsigned short int uint16
-	ctypedef signed short int sint16
-	ctypedef unsigned char uchar8
-	ctypedef signed char schar8
-	ctypedef int int32
-	ctypedef short int int16
+    ctypedef float float32
+    ctypedef double float64
+    ctypedef unsigned short int uint16
+    ctypedef signed short int sint16
+    ctypedef unsigned char uchar8
+    ctypedef signed char schar8
+    ctypedef int int32
+    ctypedef short int int16
 
 cdef extern from "astra/Config.h" namespace "astra":
-	cdef cppclass Config:
-		Config()
-		void initialize(string rootname)
-		XMLNode self
+    cdef cppclass Config:
+        Config()
+        void initialize(string rootname)
+        XMLNode self
 
 cdef extern from "astra/VolumeGeometry2D.h" namespace "astra":
-	cdef cppclass CVolumeGeometry2D:
-		bool initialize(Config)
-		int getGridColCount()
-		int getGridRowCount()
-		int getGridTotCount()
-		float32 getWindowLengthX()
-		float32 getWindowLengthY()
-		float32 getWindowArea()
-		float32 getPixelLengthX()
-		float32 getPixelLengthY()
-		float32 getPixelArea()
-		float32 getWindowMinX()
-		float32 getWindowMinY()
-		float32 getWindowMaxX()
-		float32 getWindowMaxY()
-		Config* getConfiguration()
-		bool isEqual(CVolumeGeometry2D*)
+    cdef cppclass CVolumeGeometry2D:
+        bool initialize(Config)
+        int getGridColCount()
+        int getGridRowCount()
+        int getGridTotCount()
+        float32 getWindowLengthX()
+        float32 getWindowLengthY()
+        float32 getWindowArea()
+        float32 getPixelLengthX()
+        float32 getPixelLengthY()
+        float32 getPixelArea()
+        float32 getWindowMinX()
+        float32 getWindowMinY()
+        float32 getWindowMaxX()
+        float32 getWindowMaxY()
+        Config* getConfiguration()
+        bool isEqual(CVolumeGeometry2D*)
 
 cdef extern from "astra/Float32Data2D.h" namespace "astra":
-	cdef cppclass CFloat32CustomMemory:
-		pass
+    cdef cppclass CFloat32CustomMemory:
+        pass
 
 cdef extern from "astra/Float32VolumeData2D.h" namespace "astra":
-	cdef cppclass CFloat32VolumeData2D:
-		CFloat32VolumeData2D(CVolumeGeometry2D*)
-		CFloat32VolumeData2D(CVolumeGeometry2D*, CFloat32CustomMemory*)
-		CVolumeGeometry2D * getGeometry()
-		int getWidth()
-		int getHeight()
-		void changeGeometry(CVolumeGeometry2D*)
-		Config* getConfiguration()
+    cdef cppclass CFloat32VolumeData2D:
+        CFloat32VolumeData2D(CVolumeGeometry2D*)
+        CFloat32VolumeData2D(CVolumeGeometry2D*, CFloat32CustomMemory*)
+        CVolumeGeometry2D * getGeometry()
+        int getWidth()
+        int getHeight()
+        void changeGeometry(CVolumeGeometry2D*)
+        Config* getConfiguration()
 
 
 
 cdef extern from "astra/ProjectionGeometry2D.h" namespace "astra":
-	cdef cppclass CProjectionGeometry2D:
-		CProjectionGeometry2D()
-		bool initialize(Config)
-		int getDetectorCount()
-		int getProjectionAngleCount()
-		bool isOfType(string)
-		float32 getProjectionAngle(int)
-		float32 getDetectorWidth()
-		Config* getConfiguration()
-		bool isEqual(CProjectionGeometry2D*)
+    cdef cppclass CProjectionGeometry2D:
+        CProjectionGeometry2D()
+        bool initialize(Config)
+        int getDetectorCount()
+        int getProjectionAngleCount()
+        bool isOfType(string)
+        float32 getProjectionAngle(int)
+        float32 getDetectorWidth()
+        Config* getConfiguration()
+        bool isEqual(CProjectionGeometry2D*)
 
 cdef extern from "astra/Float32Data2D.h" namespace "astra::CFloat32Data2D":
-	cdef enum TWOEDataType "astra::CFloat32Data2D::EDataType":
-		TWOPROJECTION "astra::CFloat32Data2D::PROJECTION"
-		TWOVOLUME "astra::CFloat32Data2D::VOLUME"
+    cdef enum TWOEDataType "astra::CFloat32Data2D::EDataType":
+        TWOPROJECTION "astra::CFloat32Data2D::PROJECTION"
+        TWOVOLUME "astra::CFloat32Data2D::VOLUME"
 
 cdef extern from "astra/Float32Data3D.h" namespace "astra::CFloat32Data3D":
-	cdef enum THREEEDataType "astra::CFloat32Data3D::EDataType":
-		THREEPROJECTION "astra::CFloat32Data3D::PROJECTION"
-		THREEVOLUME "astra::CFloat32Data3D::VOLUME"
+    cdef enum THREEEDataType "astra::CFloat32Data3D::EDataType":
+        THREEPROJECTION "astra::CFloat32Data3D::PROJECTION"
+        THREEVOLUME "astra::CFloat32Data3D::VOLUME"
 
 
 cdef extern from "astra/Float32Data2D.h" namespace "astra":
-	cdef cppclass CFloat32Data2D:
-		bool isInitialized()
-		int getSize()
-		float32 *getData()
-		float32 **getData2D()
-		int getWidth()
-		int getHeight()
-		TWOEDataType getType()
+    cdef cppclass CFloat32Data2D:
+        bool isInitialized()
+        int getSize()
+        float32 *getData()
+        float32 **getData2D()
+        int getWidth()
+        int getHeight()
+        TWOEDataType getType()
 
 
 
 
 cdef extern from "astra/SparseMatrixProjectionGeometry2D.h" namespace "astra":
-	cdef cppclass CSparseMatrixProjectionGeometry2D:
-		CSparseMatrixProjectionGeometry2D()
+    cdef cppclass CSparseMatrixProjectionGeometry2D:
+        CSparseMatrixProjectionGeometry2D()
 
 cdef extern from "astra/FanFlatProjectionGeometry2D.h" namespace "astra":
-	cdef cppclass CFanFlatProjectionGeometry2D:
-		CFanFlatProjectionGeometry2D()
+    cdef cppclass CFanFlatProjectionGeometry2D:
+        CFanFlatProjectionGeometry2D()
 
 cdef extern from "astra/FanFlatVecProjectionGeometry2D.h" namespace "astra":
-	cdef cppclass CFanFlatVecProjectionGeometry2D:
-		CFanFlatVecProjectionGeometry2D()
+    cdef cppclass CFanFlatVecProjectionGeometry2D:
+        CFanFlatVecProjectionGeometry2D()
 
 cdef extern from "astra/ParallelProjectionGeometry2D.h" namespace "astra":
-	cdef cppclass CParallelProjectionGeometry2D:
-		CParallelProjectionGeometry2D()
+    cdef cppclass CParallelProjectionGeometry2D:
+        CParallelProjectionGeometry2D()
 
 
 cdef extern from "astra/Float32ProjectionData2D.h" namespace "astra":
-	cdef cppclass CFloat32ProjectionData2D:
-		CFloat32ProjectionData2D(CProjectionGeometry2D*)
-		CFloat32ProjectionData2D(CProjectionGeometry2D*, CFloat32CustomMemory*)
-		CProjectionGeometry2D * getGeometry()
-		void changeGeometry(CProjectionGeometry2D*)
-		int getDetectorCount()
-		int getAngleCount()
+    cdef cppclass CFloat32ProjectionData2D:
+        CFloat32ProjectionData2D(CProjectionGeometry2D*)
+        CFloat32ProjectionData2D(CProjectionGeometry2D*, CFloat32CustomMemory*)
+        CProjectionGeometry2D * getGeometry()
+        void changeGeometry(CProjectionGeometry2D*)
+        int getDetectorCount()
+        int getAngleCount()
 
 cdef extern from "astra/Algorithm.h" namespace "astra":
-	cdef cppclass CAlgorithm:
-		bool initialize(Config)
-		void run(int) nogil
-		bool isInitialized()
+    cdef cppclass CAlgorithm:
+        bool initialize(Config)
+        void run(int) nogil
+        bool isInitialized()
 
 cdef extern from "astra/ReconstructionAlgorithm2D.h" namespace "astra":
-	cdef cppclass CReconstructionAlgorithm2D:
-		bool getResidualNorm(float32&)
+    cdef cppclass CReconstructionAlgorithm2D:
+        bool getResidualNorm(float32&)
 
 cdef extern from "astra/Projector2D.h" namespace "astra":
-	cdef cppclass CProjector2D:
-		bool isInitialized()
-		CProjectionGeometry2D* getProjectionGeometry()
-		CVolumeGeometry2D* getVolumeGeometry()
-		CSparseMatrix* getMatrix()
+    cdef cppclass CProjector2D:
+        bool isInitialized()
+        CProjectionGeometry2D* getProjectionGeometry()
+        CVolumeGeometry2D* getVolumeGeometry()
+        CSparseMatrix* getMatrix()
 
 cdef extern from "astra/Projector3D.h" namespace "astra":
-	cdef cppclass CProjector3D:
-		bool isInitialized()
-		CProjectionGeometry3D* getProjectionGeometry()
-		CVolumeGeometry3D* getVolumeGeometry()
+    cdef cppclass CProjector3D:
+        bool isInitialized()
+        CProjectionGeometry3D* getProjectionGeometry()
+        CVolumeGeometry3D* getVolumeGeometry()
 
 IF HAVE_CUDA==True:
-	cdef extern from "astra/CudaProjector3D.h" namespace "astra":
-		cdef cppclass CCudaProjector3D
+    cdef extern from "astra/CudaProjector3D.h" namespace "astra":
+        cdef cppclass CCudaProjector3D
 
-	cdef extern from "astra/CudaProjector2D.h" namespace "astra":
-		cdef cppclass CCudaProjector2D
+    cdef extern from "astra/CudaProjector2D.h" namespace "astra":
+        cdef cppclass CCudaProjector2D
+
+    cdef extern from "astra/Float32Data3DGPU.h" namespace "astraCUDA3d":
+        cdef cppclass MemHandle3D:
+            pass
+
+    cdef extern from "astra/Float32Data3DGPU.h" namespace "astraCUDA3d":
+        cdef MemHandle3D wrapHandle(float *D_ptr, unsigned int x, unsigned int y, unsigned int z, unsigned int pitch)
 
 
 cdef extern from "astra/SparseMatrix.h" namespace "astra":
-	cdef cppclass CSparseMatrix:
-		CSparseMatrix(unsigned int,unsigned int,unsigned long)
-		unsigned int m_iWidth
-		unsigned int m_iHeight
-		unsigned long m_lSize
-		bool isInitialized()
-		float32* m_pfValues
-		unsigned int* m_piColIndices
-		unsigned long* m_plRowStarts
+    cdef cppclass CSparseMatrix:
+        CSparseMatrix(unsigned int,unsigned int,unsigned long)
+        unsigned int m_iWidth
+        unsigned int m_iHeight
+        unsigned long m_lSize
+        bool isInitialized()
+        float32* m_pfValues
+        unsigned int* m_piColIndices
+        unsigned long* m_plRowStarts
+
+cdef extern from "astra/Float32Data3D.h" namespace "astra":
+    cdef cppclass CFloat32Data3D:
+        bool isInitialized()
+        int getSize()
+        int getWidth()
+        int getHeight()
+        int getDepth()
+
 
 cdef extern from "astra/Float32Data3DMemory.h" namespace "astra":
-	cdef cppclass CFloat32Data3DMemory:
-		CFloat32Data3DMemory()
-		bool isInitialized()
-		int getSize()
-		int getWidth()
-		int getHeight()
-		int getDepth()
-		void updateStatistics()
-		float32 *getData()
-		float32 ***getData3D()
-		THREEEDataType getType()
+    cdef cppclass CFloat32Data3DMemory(CFloat32Data3D):
+        CFloat32Data3DMemory()
+        void updateStatistics()
+        float32 *getData()
+        float32 ***getData3D()
+        THREEEDataType getType()
 
 
 cdef extern from "astra/VolumeGeometry3D.h" namespace "astra":
-	cdef cppclass CVolumeGeometry3D:
-		CVolumeGeometry3D()
-		bool initialize(Config)
-		Config * getConfiguration()
-		int getGridColCount()
-		int getGridRowCount()
-		int getGridSliceCount()
+    cdef cppclass CVolumeGeometry3D:
+        CVolumeGeometry3D()
+        bool initialize(Config)
+        Config * getConfiguration()
+        int getGridColCount()
+        int getGridRowCount()
+        int getGridSliceCount()
 
 cdef extern from "astra/ProjectionGeometry3D.h" namespace "astra":
-	cdef cppclass CProjectionGeometry3D:
-		CProjectionGeometry3D()
-		bool initialize(Config)
-		Config * getConfiguration()
-		int getProjectionCount()
-		int getDetectorColCount()
-		int getDetectorRowCount()
+    cdef cppclass CProjectionGeometry3D:
+        CProjectionGeometry3D()
+        bool initialize(Config)
+        Config * getConfiguration()
+        int getProjectionCount()
+        int getDetectorColCount()
+        int getDetectorRowCount()
 
 
 cdef extern from "astra/Float32VolumeData3DMemory.h" namespace "astra":
-	cdef cppclass CFloat32VolumeData3DMemory:
-		CFloat32VolumeData3DMemory(CVolumeGeometry3D*)
-		CFloat32VolumeData3DMemory(CVolumeGeometry3D*, CFloat32CustomMemory*)
-		CVolumeGeometry3D* getGeometry()
-		void changeGeometry(CVolumeGeometry3D*)
-		int getRowCount()
-		int getColCount()
-		int getSliceCount()
-		bool isInitialized()
-
-
+    cdef cppclass CFloat32VolumeData3DMemory:
+        CFloat32VolumeData3DMemory(CVolumeGeometry3D*)
+        CFloat32VolumeData3DMemory(CVolumeGeometry3D*, CFloat32CustomMemory*)
+        CVolumeGeometry3D* getGeometry()
+        void changeGeometry(CVolumeGeometry3D*)
+        int getRowCount()
+        int getColCount()
+        int getSliceCount()
+        bool isInitialized()
 
 cdef extern from "astra/ParallelProjectionGeometry3D.h" namespace "astra":
-	cdef cppclass CParallelProjectionGeometry3D:
-		CParallelProjectionGeometry3D()
+    cdef cppclass CParallelProjectionGeometry3D:
+        CParallelProjectionGeometry3D()
 
 cdef extern from "astra/ParallelVecProjectionGeometry3D.h" namespace "astra":
-	cdef cppclass CParallelVecProjectionGeometry3D:
-		CParallelVecProjectionGeometry3D()
+    cdef cppclass CParallelVecProjectionGeometry3D:
+        CParallelVecProjectionGeometry3D()
 
 cdef extern from "astra/ConeProjectionGeometry3D.h" namespace "astra":
-	cdef cppclass CConeProjectionGeometry3D:
-		CConeProjectionGeometry3D()
-		bool initialize(Config)
+    cdef cppclass CConeProjectionGeometry3D:
+        CConeProjectionGeometry3D()
+        bool initialize(Config)
 
 cdef extern from "astra/ConeVecProjectionGeometry3D.h" namespace "astra":
-	cdef cppclass CConeVecProjectionGeometry3D:
-		CConeVecProjectionGeometry3D()
+    cdef cppclass CConeVecProjectionGeometry3D:
+        CConeVecProjectionGeometry3D()
 
 cdef extern from "astra/Float32ProjectionData3DMemory.h" namespace "astra":
-	cdef cppclass CFloat32ProjectionData3DMemory:
-		CFloat32ProjectionData3DMemory(CProjectionGeometry3D*)
-		CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*)
-		CFloat32ProjectionData3DMemory(CProjectionGeometry3D*, CFloat32CustomMemory*)
-		CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*, CFloat32CustomMemory*)
-		CProjectionGeometry3D* getGeometry()
-		void changeGeometry(CProjectionGeometry3D*)
-		int getDetectorColCount()
-		int getDetectorRowCount()
-		int getAngleCount()
-		bool isInitialized()
+    cdef cppclass CFloat32ProjectionData3DMemory:
+        CFloat32ProjectionData3DMemory(CProjectionGeometry3D*)
+        CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*)
+        CFloat32ProjectionData3DMemory(CProjectionGeometry3D*, CFloat32CustomMemory*)
+        CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*, CFloat32CustomMemory*)
+        CProjectionGeometry3D* getGeometry()
+        void changeGeometry(CProjectionGeometry3D*)
+        int getDetectorColCount()
+        int getDetectorRowCount()
+        int getAngleCount()
+        bool isInitialized()
 
-cdef extern from "astra/Float32Data3D.h" namespace "astra":
-	cdef cppclass CFloat32Data3D:
-		CFloat32Data3D()
-		bool isInitialized()
-		int getSize()
-		int getWidth()
-		int getHeight()
-		int getDepth()
-		void updateStatistics()
+IF HAVE_CUDA==True:
+    cdef extern from "astra/Float32VolumeData3DGPU.h" namespace "astra":
+        cdef cppclass CFloat32VolumeData3DGPU:
+            CFloat32VolumeData3DGPU(CVolumeGeometry3D*, MemHandle3D)
+            CVolumeGeometry3D* getGeometry()
+            void changeGeometry(CVolumeGeometry3D*)
+            int getRowCount()
+            int getColCount()
+            int getSliceCount()
+            bool isInitialized()
+
+    cdef extern from "astra/Float32ProjectionData3DGPU.h" namespace "astra":
+        cdef cppclass CFloat32ProjectionData3DGPU:
+            CFloat32ProjectionData3DGPU(CProjectionGeometry3D*, MemHandle3D)
+            CProjectionGeometry3D* getGeometry()
+            void changeGeometry(CProjectionGeometry3D*)
+            int getRowCount()
+            int getColCount()
+            int getSliceCount()
+            bool isInitialized()
diff --git a/python/astra/creators.py b/python/astra/creators.py
index f45fd52..aa16221 100644
--- a/python/astra/creators.py
+++ b/python/astra/creators.py
@@ -258,7 +258,7 @@ This method can be called in a number of ways:
     elif intype == 'cone':
         if len(args) < 7:
             raise Exception('not enough variables: astra_create_proj_geom(cone, detector_spacing_x, detector_spacing_y, det_row_count, det_col_count, angles, source_origin, origin_det)')
-        return {'type':	'cone','DetectorSpacingX':args[0], 'DetectorSpacingY':args[1], 'DetectorRowCount':args[2],'DetectorColCount':args[3],'ProjectionAngles':args[4],'DistanceOriginSource':	args[5],'DistanceOriginDetector':args[6]}
+        return {'type': 'cone','DetectorSpacingX':args[0], 'DetectorSpacingY':args[1], 'DetectorRowCount':args[2],'DetectorColCount':args[3],'ProjectionAngles':args[4],'DistanceOriginSource': args[5],'DistanceOriginDetector':args[6]}
     elif intype == 'cone_vec':
         if len(args) < 3:
             raise Exception('not enough variables: astra_create_proj_geom(cone_vec, det_row_count, det_col_count, V)')
diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx
index 203fde2..9c88073 100644
--- a/python/astra/data2d_c.pyx
+++ b/python/astra/data2d_c.pyx
@@ -83,8 +83,10 @@ def create(datatype, geometry, data=None, link=False):
     cdef CFloat32Data2D * pDataObject2D
     cdef CFloat32CustomMemory * pCustom
 
-    if link and data.shape!=geom_size(geometry):
-        raise Exception("The dimensions of the data do not match those specified in the geometry.")
+    if link:
+        geom_shape = geom_size(geometry)
+        if data.shape != geom_shape:
+            raise ValueError("The dimensions of the data do not match those specified in the geometry: {} != {}".format(data.shape, geom_shape))
 
     if datatype == '-vol':
         cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry)
@@ -92,7 +94,7 @@ def create(datatype, geometry, data=None, link=False):
         if not pGeometry.initialize(cfg[0]):
             del cfg
             del pGeometry
-            raise Exception('Geometry class not initialized.')
+            raise RuntimeError('Geometry class not initialized.')
         if link:
             pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data)
             pDataObject2D = <CFloat32Data2D * > new CFloat32VolumeData2D(pGeometry, pCustom)
@@ -114,7 +116,7 @@ def create(datatype, geometry, data=None, link=False):
         if not ppGeometry.initialize(cfg[0]):
             del cfg
             del ppGeometry
-            raise Exception('Geometry class not initialized.')
+            raise RuntimeError('Geometry class not initialized.')
         if link:
             pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data)
             pDataObject2D = <CFloat32Data2D * > new CFloat32ProjectionData2D(ppGeometry, pCustom)
@@ -123,11 +125,11 @@ def create(datatype, geometry, data=None, link=False):
         del ppGeometry
         del cfg
     else:
-        raise Exception("Invalid datatype.  Please specify '-vol' or '-sino'.")
+        raise ValueError("Invalid datatype.  Please specify '-vol' or '-sino'.")
 
     if not pDataObject2D.isInitialized():
         del pDataObject2D
-        raise Exception("Couldn't initialize data object.")
+        raise RuntimeError("Couldn't initialize data object.")
 
     if not link: fillDataObject(pDataObject2D, data)
 
@@ -138,6 +140,10 @@ cdef fillDataObject(CFloat32Data2D * obj, data):
         fillDataObjectScalar(obj, 0)
     else:
         if isinstance(data, np.ndarray):
+            obj_shape = (obj.getHeight(), obj.getWidth())
+            if data.shape != obj_shape:
+                raise ValueError(
+                  "The dimensions of the data do not match those specified in the geometry: {} != {}".format(data.shape, obj_shape))
             fillDataObjectArray(obj, np.ascontiguousarray(data,dtype=np.float32))
         else:
             fillDataObjectScalar(obj, np.float32(data))
@@ -150,18 +156,15 @@ cdef fillDataObjectScalar(CFloat32Data2D * obj, float s):
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef fillDataObjectArray(CFloat32Data2D * obj, float [:,::1] data):
-    if (not data.shape[0] == obj.getHeight()) or (not data.shape[1] == obj.getWidth()):
-        raise Exception(
-            "The dimensions of the data do not match those specified in the geometry.")
     cdef float [:,::1] cView =  <float[:data.shape[0],:data.shape[1]]> obj.getData2D()[0]
     cView[:] = data
 
 cdef CFloat32Data2D * getObject(i) except NULL:
     cdef CFloat32Data2D * pDataObject = man2d.get(i)
     if pDataObject == NULL:
-        raise Exception("Data object not found")
+        raise ValueError("Data object not found")
     if not pDataObject.isInitialized():
-        raise Exception("Data object not initialized properly.")
+        raise RuntimeError("Data object not initialized properly.")
     return pDataObject
 
 
@@ -180,15 +183,15 @@ def get_geometry(i):
         pDataObject3 = <CFloat32VolumeData2D * >pDataObject
         geom = utils.configToDict(pDataObject3.getGeometry().getConfiguration())
     else:
-        raise Exception("Not a known data object")
+        raise RuntimeError("Not a known data object")
     return geom
 
 cdef CProjector2D * getProjector(i) except NULL:
     cdef CProjector2D * proj = manProj.get(i)
     if proj == NULL:
-        raise Exception("Projector not initialized.")
+        raise RuntimeError("Projector not initialized.")
     if not proj.isInitialized():
-        raise Exception("Projector not initialized.")
+        raise RuntimeError("Projector not initialized.")
     return proj
 
 def check_compatible(i, proj_id):
@@ -203,7 +206,7 @@ def check_compatible(i, proj_id):
         pDataObject3 = <CFloat32VolumeData2D * >pDataObject
         return pDataObject3.getGeometry().isEqual(proj.getVolumeGeometry())
     else:
-        raise Exception("Not a known data object")
+        raise RuntimeError("Not a known data object")
 
 def change_geometry(i, geom):
     cdef Config *cfg
@@ -227,12 +230,14 @@ def change_geometry(i, geom):
         if not ppGeometry.initialize(cfg[0]):
             del cfg
             del ppGeometry
-            raise Exception('Geometry class not initialized.')
-        if (ppGeometry.getDetectorCount() != pDataObject2.getDetectorCount() or ppGeometry.getProjectionAngleCount() != pDataObject2.getAngleCount()):
+            raise RuntimeError('Geometry class not initialized.')
+        geom_shape = (ppGeometry.getProjectionAngleCount(), ppGeometry.getDetectorCount())
+        obj_shape = (pDataObject2.getAngleCount(), pDataObject2.getDetectorCount())
+        if geom_shape != obj_shape:
             del ppGeometry
             del cfg
-            raise Exception(
-                "The dimensions of the data do not match those specified in the geometry.")
+            raise ValueError(
+                "The dimensions of the data do not match those specified in the geometry: {} != {}", obj_shape, geom_shape)
         pDataObject2.changeGeometry(ppGeometry)
         del ppGeometry
         del cfg
@@ -243,17 +248,19 @@ def change_geometry(i, geom):
         if not pGeometry.initialize(cfg[0]):
             del cfg
             del pGeometry
-            raise Exception('Geometry class not initialized.')
-        if (pGeometry.getGridColCount() != pDataObject3.getWidth() or pGeometry.getGridRowCount() != pDataObject3.getHeight()):
+            raise RuntimeError('Geometry class not initialized.')
+        geom_shape = (pGeometry.getGridRowCount(), pGeometry.getGridColCount())
+        obj_shape = (pDataObject3.getHeight(), pDataObject3.getWidth())
+        if geom_shape != obj_shape:
             del cfg
             del pGeometry
-            raise Exception(
-                'The dimensions of the data do not match those specified in the geometry.')
+            raise ValueError(
+                "The dimensions of the data do not match those specified in the geometry: {} != {}", obj_shape, geom_shape)
         pDataObject3.changeGeometry(pGeometry)
         del cfg
         del pGeometry
     else:
-        raise Exception("Not a known data object")
+        raise RuntimeError("Not a known data object")
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
@@ -274,7 +281,7 @@ def get_shared(i):
 
 
 def get_single(i):
-    raise Exception("Not yet implemented")
+    raise NotImplementedError("Not yet implemented")
 
 
 def info():
diff --git a/python/astra/data3d.py b/python/astra/data3d.py
index 9c2bea4..ecb99d0 100644
--- a/python/astra/data3d.py
+++ b/python/astra/data3d.py
@@ -26,6 +26,8 @@
 from . import data3d_c as d
 import numpy as np
 
+from .pythonutils import GPULink
+
 def create(datatype,geometry,data=None):
     """Create a 3D object.
 
@@ -52,12 +54,13 @@ def link(datatype, geometry, data):
     :returns: :class:`int` -- the ID of the constructed object.
     
     """
-    if not isinstance(data,np.ndarray):
-        raise ValueError("Input should be a numpy array")
-    if not data.dtype==np.float32:
-        raise ValueError("Numpy array should be float32")
-    if not (data.flags['C_CONTIGUOUS'] and data.flags['ALIGNED']):
-        raise ValueError("Numpy array should be C_CONTIGUOUS and ALIGNED")
+    if not isinstance(data,np.ndarray) and not isinstance(data,GPULink):
+        raise TypeError("Input should be a numpy ndarray or GPULink object")
+    if isinstance(data, np.ndarray):
+        if data.dtype != np.float32:
+            raise ValueError("Numpy array should be float32")
+        if not (data.flags['C_CONTIGUOUS'] and data.flags['ALIGNED']):
+            raise ValueError("Numpy array should be C_CONTIGUOUS and ALIGNED")
     return d.create(datatype,geometry,data,True)
 
 
diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx
index 915c60d..78ed620 100644
--- a/python/astra/data3d_c.pyx
+++ b/python/astra/data3d_c.pyx
@@ -45,12 +45,17 @@ from .PyXMLDocument cimport XMLDocument
 cimport utils
 from .utils import wrap_from_bytes
 
-from .pythonutils import geom_size
+from .pythonutils import geom_size, GPULink
 
 import operator
 
 from six.moves import reduce
 
+include "config.pxi"
+
+cdef extern from "Python.h":
+    void* PyLong_AsVoidPtr(object)
+
 
 cdef CData3DManager * man3d = <CData3DManager * >PyData3DManager.getSingletonPtr()
 
@@ -65,12 +70,22 @@ def create(datatype,geometry,data=None, link=False):
     cdef Config *cfg
     cdef CVolumeGeometry3D * pGeometry
     cdef CProjectionGeometry3D * ppGeometry
-    cdef CFloat32Data3DMemory * pDataObject3D
+    cdef CFloat32Data3D * pDataObject3D
     cdef CConeProjectionGeometry3D* pppGeometry
-    cdef CFloat32CustomMemory * pCustom
+    cdef CFloat32CustomMemory * pCustom = NULL
+    IF HAVE_CUDA==True:
+        cdef MemHandle3D hnd
 
-    if link and data.shape!=geom_size(geometry):
-        raise Exception("The dimensions of the data do not match those specified in the geometry.")
+    if link:
+        geom_shape = geom_size(geometry)
+        if isinstance(data, np.ndarray):
+            data_shape = data.shape
+        elif isinstance(data, GPULink):
+            data_shape = ( data.z, data.y, data.x )
+        else:
+            raise TypeError("data should be a numpy.ndarray or a GPULink object")
+        if geom_shape != data_shape:
+            raise ValueError("The dimensions of the data do not match those specified in the geometry: {} != {}".format(data_shape, geom_shape))
 
     if datatype == '-vol':
         cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry)
@@ -78,15 +93,25 @@ def create(datatype,geometry,data=None, link=False):
         if not pGeometry.initialize(cfg[0]):
             del cfg
             del pGeometry
-            raise Exception('Geometry class not initialized.')
+            raise RuntimeError('Geometry class not initialized.')
         if link:
-            pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data)
-            pDataObject3D = <CFloat32Data3DMemory * > new CFloat32VolumeData3DMemory(pGeometry, pCustom)
+            if isinstance(data, np.ndarray):
+                pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data)
+                pDataObject3D = <CFloat32Data3D * > new CFloat32VolumeData3DMemory(pGeometry, pCustom)
+            elif isinstance(data, GPULink):
+                IF HAVE_CUDA==True:
+                    s = geom_size(geometry)
+                    hnd = wrapHandle(<float*>PyLong_AsVoidPtr(data.ptr), data.x, data.y, data.z, data.pitch/4)
+                    pDataObject3D = <CFloat32Data3D * > new CFloat32VolumeData3DGPU(pGeometry, hnd)
+                ELSE:
+                    raise NotImplementedError("CUDA support is not enabled in ASTRA")
+            else:
+                raise TypeError("data should be a numpy.ndarray or a GPULink object")
         else:
-            pDataObject3D = <CFloat32Data3DMemory * > new CFloat32VolumeData3DMemory(pGeometry)
+            pDataObject3D = <CFloat32Data3D * > new CFloat32VolumeData3DMemory(pGeometry)
         del cfg
         del pGeometry
-    elif datatype == '-sino' or datatype == '-proj3d':
+    elif datatype == '-sino' or datatype == '-proj3d' or datatype == '-sinocone':
         cfg = utils.dictToConfig(six.b('ProjectionGeometry'), geometry)
         tpe = wrap_from_bytes(cfg.self.getAttribute(six.b('type')))
         if (tpe == "parallel3d"):
@@ -98,41 +123,38 @@ def create(datatype,geometry,data=None, link=False):
         elif (tpe == "cone_vec"):
             ppGeometry = <CProjectionGeometry3D*> new CConeVecProjectionGeometry3D();
         else:
-            raise Exception("Invalid geometry type.")
+            raise ValueError("Invalid geometry type.")
 
         if not ppGeometry.initialize(cfg[0]):
             del cfg
             del ppGeometry
-            raise Exception('Geometry class not initialized.')
+            raise RuntimeError('Geometry class not initialized.')
         if link:
-            pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data)
-            pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(ppGeometry, pCustom)
+            if isinstance(data, np.ndarray):
+                pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data)
+                pDataObject3D = <CFloat32Data3D * > new CFloat32ProjectionData3DMemory(ppGeometry, pCustom)
+            elif isinstance(data, GPULink):
+                IF HAVE_CUDA==True:
+                    s = geom_size(geometry)
+                    hnd = wrapHandle(<float*>PyLong_AsVoidPtr(data.ptr), data.x, data.y, data.z, data.pitch/4)
+                    pDataObject3D = <CFloat32Data3D * > new CFloat32ProjectionData3DGPU(ppGeometry, hnd)
+                ELSE:
+                    raise NotImplementedError("CUDA support is not enabled in ASTRA")
+            else:
+                raise TypeError("data should be a numpy.ndarray or a GPULink object")
         else:
             pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(ppGeometry)
         del ppGeometry
         del cfg
-    elif datatype == "-sinocone":
-        cfg = utils.dictToConfig(six.b('ProjectionGeometry'), geometry)
-        pppGeometry = new CConeProjectionGeometry3D()
-        if not pppGeometry.initialize(cfg[0]):
-            del cfg
-            del pppGeometry
-            raise Exception('Geometry class not initialized.')
-        if link:
-            pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data)
-            pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(pppGeometry, pCustom)
-        else:
-            pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(pppGeometry)
     else:
-        raise Exception("Invalid datatype.  Please specify '-vol' or '-proj3d'.")
+        raise ValueError("Invalid datatype.  Please specify '-vol' or '-proj3d'.")
 
     if not pDataObject3D.isInitialized():
         del pDataObject3D
-        raise Exception("Couldn't initialize data object.")
-
-    if not link: fillDataObject(pDataObject3D, data)
+        raise RuntimeError("Couldn't initialize data object.")
 
-    pDataObject3D.updateStatistics()
+    if not link:
+        fillDataObject(dynamic_cast_mem(pDataObject3D), data)
 
     return man3d.store(<CFloat32Data3D*>pDataObject3D)
 
@@ -147,7 +169,7 @@ def get_geometry(i):
         pDataObject3 = <CFloat32VolumeData3DMemory * >pDataObject
         geom = utils.configToDict(pDataObject3.getGeometry().getConfiguration())
     else:
-        raise Exception("Not a known data object")
+        raise RuntimeError("Not a known data object")
     return geom
 
 def change_geometry(i, geom):
@@ -168,18 +190,18 @@ def change_geometry(i, geom):
         elif (tpe == "cone_vec"):
             ppGeometry = <CProjectionGeometry3D*> new CConeVecProjectionGeometry3D();
         else:
-            raise Exception("Invalid geometry type.")
+            raise ValueError("Invalid geometry type.")
         if not ppGeometry.initialize(cfg[0]):
             del cfg
             del ppGeometry
-            raise Exception('Geometry class not initialized.')
+            raise RuntimeError('Geometry class not initialized.')
         del cfg
-        if (ppGeometry.getDetectorColCount() != pDataObject2.getDetectorColCount() or \
-            ppGeometry.getProjectionCount() != pDataObject2.getAngleCount() or \
-            ppGeometry.getDetectorRowCount() != pDataObject2.getDetectorRowCount()):
+        geom_shape = (ppGeometry.getDetectorRowCount(), ppGeometry.getProjectionCount(), ppGeometry.getDetectorColCount())
+        obj_shape = (pDataObject2.getDetectorRowCount(), pDataObject2.getAngleCount(), pDataObject2.getDetectorColCount())
+        if geom_shape != obj_shape:
             del ppGeometry
-            raise Exception(
-                "The dimensions of the data do not match those specified in the geometry.")
+            raise ValueError(
+                "The dimensions of the data do not match those specified in the geometry: {} != {}".format(obj_shape, geom_shape))
         pDataObject2.changeGeometry(ppGeometry)
         del ppGeometry
 
@@ -190,19 +212,19 @@ def change_geometry(i, geom):
         if not pGeometry.initialize(cfg[0]):
             del cfg
             del pGeometry
-            raise Exception('Geometry class not initialized.')
+            raise RuntimeError('Geometry class not initialized.')
         del cfg
-        if (pGeometry.getGridColCount() != pDataObject3.getColCount() or \
-            pGeometry.getGridRowCount() != pDataObject3.getRowCount() or \
-            pGeometry.getGridSliceCount() != pDataObject3.getSliceCount()):
+        geom_shape = (pGeometry.getGridSliceCount(), pGeometry.getGridRowCount(), pGeometry.getGridColCount())
+        obj_shape = (pDataObject3.getSliceCount(), pDataObject3.getRowCount(), pDataObject3.getColCount())
+        if geom_shape != obj_shape:
             del pGeometry
-            raise Exception(
-                "The dimensions of the data do not match those specified in the geometry.")
+            raise ValueError(
+                "The dimensions of the data do not match those specified in the geometry.".format(obj_shape, geom_shape))
         pDataObject3.changeGeometry(pGeometry)
         del pGeometry
 
     else:
-        raise Exception("Not a known data object")
+        raise RuntimeError("Not a known data object")
 
 
 cdef fillDataObject(CFloat32Data3DMemory * obj, data):
@@ -210,6 +232,10 @@ cdef fillDataObject(CFloat32Data3DMemory * obj, data):
         fillDataObjectScalar(obj, 0)
     else:
         if isinstance(data, np.ndarray):
+            obj_shape = (obj.getDepth(), obj.getHeight(), obj.getWidth())
+            if data.shape != obj_shape:
+                raise ValueError(
+                  "The dimensions of the data do not match those specified in the geometry: {} != {}".format(data.shape, obj_shape))
             fillDataObjectArray(obj, np.ascontiguousarray(data,dtype=np.float32))
         else:
             fillDataObjectScalar(obj, np.float32(data))
@@ -222,18 +248,15 @@ cdef fillDataObjectScalar(CFloat32Data3DMemory * obj, float s):
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef fillDataObjectArray(CFloat32Data3DMemory * obj, float [:,:,::1] data):
-    if (not data.shape[0] == obj.getDepth()) or (not data.shape[1] == obj.getHeight()) or (not data.shape[2] == obj.getWidth()):
-        raise Exception(
-            "The dimensions of the data do not match those specified in the geometry.")
     cdef float [:,:,::1] cView = <float[:data.shape[0],:data.shape[1],:data.shape[2]]> obj.getData3D()[0][0]
     cView[:] = data
 
 cdef CFloat32Data3D * getObject(i) except NULL:
     cdef CFloat32Data3D * pDataObject = man3d.get(i)
     if pDataObject == NULL:
-        raise Exception("Data object not found")
+        raise ValueError("Data object not found")
     if not pDataObject.isInitialized():
-        raise Exception("Data object not initialized properly.")
+        raise RuntimeError("Data object not initialized properly.")
     return pDataObject
 
 @cython.boundscheck(False)
@@ -256,7 +279,7 @@ def get_shared(i):
     return np.PyArray_SimpleNewFromData(3,shape,np.NPY_FLOAT32,<void *>pDataObject.getData3D()[0][0])
 
 def get_single(i):
-    raise Exception("Not yet implemented")
+    raise NotImplementedError("Not yet implemented")
 
 def store(i,data):
     cdef CFloat32Data3D * pDataObject = getObject(i)
diff --git a/python/astra/experimental.pyx b/python/astra/experimental.pyx
index b6c7881..0af3118 100644
--- a/python/astra/experimental.pyx
+++ b/python/astra/experimental.pyx
@@ -36,12 +36,20 @@ IF HAVE_CUDA==True:
 
     cdef extern from "astra/CompositeGeometryManager.h" namespace "astra":
         cdef cppclass CCompositeGeometryManager:
-            bool doFP(CProjector3D *, vector[CFloat32VolumeData3DMemory *], vector[CFloat32ProjectionData3DMemory *])
-            bool doBP(CProjector3D *, vector[CFloat32VolumeData3DMemory *], vector[CFloat32ProjectionData3DMemory *])
+            bool doFP(CProjector3D *, vector[CFloat32VolumeData3D *], vector[CFloat32ProjectionData3D *])
+            bool doBP(CProjector3D *, vector[CFloat32VolumeData3D *], vector[CFloat32ProjectionData3D *])
 
     cdef extern from *:
-        CFloat32VolumeData3DMemory * dynamic_cast_vol_mem "dynamic_cast<astra::CFloat32VolumeData3DMemory*>" (CFloat32Data3D * ) except NULL
-        CFloat32ProjectionData3DMemory * dynamic_cast_proj_mem "dynamic_cast<astra::CFloat32ProjectionData3DMemory*>" (CFloat32Data3D * ) except NULL
+        CFloat32VolumeData3D * dynamic_cast_vol_mem "dynamic_cast<astra::CFloat32VolumeData3D*>" (CFloat32Data3D * ) except NULL
+        CFloat32ProjectionData3D * dynamic_cast_proj_mem "dynamic_cast<astra::CFloat32ProjectionData3D*>" (CFloat32Data3D * ) except NULL
+
+    cdef extern from "astra/Float32ProjectionData3D.h" namespace "astra":
+        cdef cppclass CFloat32ProjectionData3D:
+            bool isInitialized()
+    cdef extern from "astra/Float32VolumeData3D.h" namespace "astra":
+        cdef cppclass CFloat32VolumeData3D:
+            bool isInitialized()
+
 
     cimport PyProjector3DManager
     from .PyProjector3DManager cimport CProjector3DManager
@@ -52,9 +60,9 @@ IF HAVE_CUDA==True:
     cdef CData3DManager * man3d = <CData3DManager * >PyData3DManager.getSingletonPtr()
 
     def do_composite(projector_id, vol_ids, proj_ids, t):
-        cdef vector[CFloat32VolumeData3DMemory *] vol
-        cdef CFloat32VolumeData3DMemory * pVolObject
-        cdef CFloat32ProjectionData3DMemory * pProjObject
+        cdef vector[CFloat32VolumeData3D *] vol
+        cdef CFloat32VolumeData3D * pVolObject
+        cdef CFloat32ProjectionData3D * pProjObject
         for v in vol_ids:
             pVolObject = dynamic_cast_vol_mem(man3d.get(v))
             if pVolObject == NULL:
@@ -62,7 +70,7 @@ IF HAVE_CUDA==True:
             if not pVolObject.isInitialized():
                 raise Exception("Data object not initialized properly")
             vol.push_back(pVolObject)
-        cdef vector[CFloat32ProjectionData3DMemory *] proj
+        cdef vector[CFloat32ProjectionData3D *] proj
         for v in proj_ids:
             pProjObject = dynamic_cast_proj_mem(man3d.get(v))
             if pProjObject == NULL:
diff --git a/python/astra/extrautils.pyx b/python/astra/extrautils.pyx
index 502cc78..248c6ee 100644
--- a/python/astra/extrautils.pyx
+++ b/python/astra/extrautils.pyx
@@ -26,19 +26,19 @@
 
 
 def clipCircle(img):
-	cdef int i,j
-	cdef double x2,y2,mid,bnd
-	cdef long sz,sz2
-	sz = img.shape[0]
-	sz2 = sz*sz
-	bnd = sz2/4.
-	mid = (sz-1.)/2.
-	nDel=0
-	for i in range(sz):
-		for j in range(sz):
-			x2 = (i-mid)*(i-mid)
-			y2 = (j-mid)*(j-mid)
-			if x2+y2>bnd:
-				img[i,j]=0
-				nDel=nDel+1
-	return nDel
+    cdef int i,j
+    cdef double x2,y2,mid,bnd
+    cdef long sz,sz2
+    sz = img.shape[0]
+    sz2 = sz*sz
+    bnd = sz2/4.
+    mid = (sz-1.)/2.
+    nDel=0
+    for i in range(sz):
+        for j in range(sz):
+            x2 = (i-mid)*(i-mid)
+            y2 = (j-mid)*(j-mid)
+            if x2+y2>bnd:
+                img[i,j]=0
+                nDel=nDel+1
+    return nDel
diff --git a/python/astra/pythonutils.py b/python/astra/pythonutils.py
index 3bd3321..27fa8fd 100644
--- a/python/astra/pythonutils.py
+++ b/python/astra/pythonutils.py
@@ -61,3 +61,21 @@ def geom_size(geom, dim=None):
         s = s[dim]
 
     return s
+
+class GPULink(object):
+    """Utility class for astra.data3d.link with a CUDA pointer
+
+    The CUDA pointer ptr must point to an array of floats.
+
+    x is the fastest-changing coordinate, z the slowest-changing.
+
+    pitch is the width in bytes of the memory block. For a contiguous
+    memory block, pitch is equal to sizeof(float) * x. For a memory block
+    allocated by cudaMalloc3D, pitch is the pitch as returned by cudaMalloc3D.
+    """
+    def __init__(self, ptr, x, y, z, pitch):
+        self.ptr = ptr
+        self.x = x
+        self.y = y
+        self.z = z
+        self.pitch = pitch
diff --git a/samples/python/s021_pygpu.py b/samples/python/s021_pygpu.py
new file mode 100644
index 0000000..db6584b
--- /dev/null
+++ b/samples/python/s021_pygpu.py
@@ -0,0 +1,72 @@
+# -----------------------------------------------------------------------
+# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp
+#            2013-2016, CWI, Amsterdam
+#
+# Contact: astra@uantwerpen.be
+# Website: http://www.astra-toolbox.com/
+#
+# This file is part of the ASTRA Toolbox.
+#
+#
+# The ASTRA Toolbox 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 3 of the License, or
+# (at your option) any later version.
+#
+# The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+#
+# -----------------------------------------------------------------------
+
+import astra
+import numpy as np
+import pygpu
+import pylab
+
+# Initialize pygpu
+ctx = pygpu.init('cuda')
+pygpu.set_default_context(ctx)
+
+vol_geom = astra.create_vol_geom(128, 128, 128)
+angles = np.linspace(0, 2 * np.pi, 180, False)
+proj_geom = astra.create_proj_geom('cone', 1.0, 1.0, 128, 192, angles, 1000, 0)
+
+# Create a simple hollow cube phantom, as a pygpu gpuarray
+vol_gpuarr = pygpu.gpuarray.zeros(astra.functions.geom_size(vol_geom), dtype='float32')
+vol_gpuarr[17:113, 17:113, 17:113] = 1
+vol_gpuarr[33:97, 33:97, 33:97] = 0
+
+# Create a pygpu gpuarray for the output projection data
+proj_gpuarr = pygpu.gpuarray.zeros(astra.functions.geom_size(proj_geom), dtype='float32')
+
+# Create the astra GPULink objects and create astra data3d objects from them
+z, y, x = proj_gpuarr.shape
+proj_data_link = astra.data3d.GPULink(proj_gpuarr.gpudata, x, y, z,
+                                      proj_gpuarr.strides[-2])
+z, y, x = vol_gpuarr.shape
+vol_link = astra.data3d.GPULink(vol_gpuarr.gpudata, x, y, z,
+                                vol_gpuarr.strides[-2])
+
+proj_id = astra.data3d.link('-sino', proj_geom, proj_data_link)
+vol_id = astra.data3d.link('-vol', vol_geom, vol_link)
+
+# Run a 3D FP
+cfg = astra.astra_dict('FP3D_CUDA')
+cfg['VolumeDataId'] = vol_id
+cfg['ProjectionDataId'] = proj_id
+alg_id = astra.algorithm.create(cfg)
+astra.algorithm.run(alg_id)
+
+pylab.figure(1)
+pylab.gray()
+pylab.imshow(proj_gpuarr[:, 20, :])
+pylab.show()
+
+astra.algorithm.delete(alg_id)
+astra.data3d.delete(vol_id)
+astra.data3d.delete(proj_id)
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index c3af228..a4dace2 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -39,6 +39,8 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/CudaProjector3D.h"
 #include "astra/Float32ProjectionData3DMemory.h"
 #include "astra/Float32VolumeData3DMemory.h"
+#include "astra/Float32ProjectionData3DGPU.h"
+#include "astra/Float32VolumeData3DGPU.h"
 #include "astra/Logging.h"
 
 #include "../cuda/3d/mem3d.h"
@@ -97,6 +99,127 @@ CCompositeGeometryManager::CCompositeGeometryManager()
 //   (First approach: 0.5/0.5)
 
 
+
+
+
+class _AstraExport CFloat32CustomGPUMemory {
+public:
+    astraCUDA3d::MemHandle3D hnd; // Only required to be valid between allocate/free
+    virtual bool allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, astraCUDA3d::Mem3DZeroMode zero)=0;
+    virtual bool copyToGPUMemory(const astraCUDA3d::SSubDimensions3D &pos)=0;
+    virtual bool copyFromGPUMemory(const astraCUDA3d::SSubDimensions3D &pos)=0;
+    virtual bool freeGPUMemory()=0;
+	virtual ~CFloat32CustomGPUMemory() { }
+};
+
+class CFloat32ExistingGPUMemory : public astra::CFloat32CustomGPUMemory {
+public:
+    CFloat32ExistingGPUMemory(CFloat32Data3DGPU *d);
+    virtual bool allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, astraCUDA3d::Mem3DZeroMode zero);
+    virtual bool copyToGPUMemory(const astraCUDA3d::SSubDimensions3D &pos);
+    virtual bool copyFromGPUMemory(const astraCUDA3d::SSubDimensions3D &pos);
+    virtual bool freeGPUMemory();
+
+protected:
+    unsigned int x, y, z;
+};
+
+class CFloat32DefaultGPUMemory : public astra::CFloat32CustomGPUMemory {
+public:
+	CFloat32DefaultGPUMemory(CFloat32Data3DMemory* d) {
+		ptr = d->getData();
+	}
+	virtual bool allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, astraCUDA3d::Mem3DZeroMode zero) {
+		hnd = astraCUDA3d::allocateGPUMemory(x, y, z, zero);
+		return (bool)hnd;
+	}
+	virtual bool copyToGPUMemory(const astraCUDA3d::SSubDimensions3D &pos) {
+		return astraCUDA3d::copyToGPUMemory(ptr, hnd, pos);
+	}
+	virtual bool copyFromGPUMemory(const astraCUDA3d::SSubDimensions3D &pos) {
+		return astraCUDA3d::copyFromGPUMemory(ptr, hnd, pos);
+	}
+	virtual bool freeGPUMemory() {
+		return astraCUDA3d::freeGPUMemory(hnd);
+	}
+
+protected:
+	float *ptr;
+};
+
+
+
+CFloat32ExistingGPUMemory::CFloat32ExistingGPUMemory(CFloat32Data3DGPU *d)
+{
+	hnd = d->getHandle();
+	x = d->getWidth();
+	y = d->getHeight();
+	z = d->getDepth();
+}
+
+bool CFloat32ExistingGPUMemory::allocateGPUMemory(unsigned int x_, unsigned int y_, unsigned int z_, astraCUDA3d::Mem3DZeroMode zero) {
+    assert(x_ == x);
+    assert(y_ == y);
+    assert(z_ == z);
+
+    if (zero == astraCUDA3d::INIT_ZERO)
+        return astraCUDA3d::zeroGPUMemory(hnd, x, y, z);
+    else
+        return true;
+}
+bool CFloat32ExistingGPUMemory::copyToGPUMemory(const astraCUDA3d::SSubDimensions3D &pos) {
+    assert(pos.nx == x);
+    assert(pos.ny == y);
+    assert(pos.nz == z);
+    assert(pos.pitch == x);
+    assert(pos.subx == 0);
+    assert(pos.suby == 0);
+    assert(pos.subnx == x);
+    assert(pos.subny == y);
+
+    // These are less necessary than x/y, but allowing access to
+    // subvolumes needs an interface change
+    assert(pos.subz == 0);
+    assert(pos.subnz == z);
+
+    return true;
+}
+bool CFloat32ExistingGPUMemory::copyFromGPUMemory(const astraCUDA3d::SSubDimensions3D &pos) {
+    assert(pos.nx == x);
+    assert(pos.ny == y);
+    assert(pos.nz == z);
+    assert(pos.pitch == x);
+    assert(pos.subx == 0);
+    assert(pos.suby == 0);
+    assert(pos.subnx == x);
+    assert(pos.subny == y);
+
+    // These are less necessary than x/y, but allowing access to
+    // subvolumes needs an interface change
+    assert(pos.subz == 0);
+    assert(pos.subnz == z);
+
+    return true;
+}
+bool CFloat32ExistingGPUMemory::freeGPUMemory() {
+    return true;
+}
+
+
+CFloat32CustomGPUMemory * createGPUMemoryHandler(CFloat32Data3D *d) {
+	CFloat32Data3DMemory *dMem = dynamic_cast<CFloat32Data3DMemory*>(d);
+	CFloat32Data3DGPU *dGPU = dynamic_cast<CFloat32Data3DGPU*>(d);
+
+	if (dMem)
+		return new CFloat32DefaultGPUMemory(dMem);
+	else
+		return new CFloat32ExistingGPUMemory(dGPU);
+}
+
+
+
+
+
 bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div, TJobSet &split)
 {
 	int maxBlockDim = astraCUDA3d::maxBlockDimension();
@@ -205,6 +328,7 @@ static std::pair<double, double> reduceProjectionVertical(const CVolumeGeometry3
 	// reduce self to only cover intersection with projection of VolumePart
 	// (Project corners of volume, take bounding box)
 
+	assert(pProjGeom->getProjectionCount() > 0);
 	for (int i = 0; i < pProjGeom->getProjectionCount(); ++i) {
 
 		double vol_u[8];
@@ -280,7 +404,7 @@ CCompositeGeometryManager::CVolumePart::~CVolumePart()
 	delete pGeom;
 }
 
-void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_t &z)
+void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_t &z) const
 {
 	if (!pGeom) {
 		x = y = z = 0;
@@ -292,13 +416,28 @@ void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_
 	z = pGeom->getGridSliceCount();
 }
 
-size_t CCompositeGeometryManager::CPart::getSize()
+size_t CCompositeGeometryManager::CPart::getSize() const
 {
 	size_t x, y, z;
 	getDims(x, y, z);
 	return x * y * z;
 }
 
+bool CCompositeGeometryManager::CPart::isFull() const
+{
+	size_t x, y, z;
+	getDims(x, y, z);
+	return x == (size_t)pData->getWidth() &&
+	       y == (size_t)pData->getHeight() &&
+	       z == (size_t)pData->getDepth();
+}
+
+bool CCompositeGeometryManager::CPart::canSplitAndReduce() const
+{
+	return dynamic_cast<CFloat32Data3DMemory *>(pData) != 0;
+}
+
+
 
 static bool testVolumeRange(const std::pair<double, double>& fullRange,
                             const CVolumeGeometry3D *pVolGeom,
@@ -334,6 +473,9 @@ static bool testVolumeRange(const std::pair<double, double>& fullRange,
 
 CCompositeGeometryManager::CPart* CCompositeGeometryManager::CVolumePart::reduce(const CPart *_other)
 {
+	if (!canSplitAndReduce())
+		return clone();
+
 	const CProjectionPart *other = dynamic_cast<const CProjectionPart *>(_other);
 	assert(other);
 
@@ -654,7 +796,7 @@ static CProjectionGeometry3D* getSubProjectionGeometryV(const CProjectionGeometr
 // - maybe all approximately the same size?
 void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div)
 {
-	if (true) {
+	if (canSplitAndReduce()) {
 		// Split in vertical direction only at first, until we figure out
 		// a model for splitting in other directions
 
@@ -664,7 +806,7 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T
 		size_t blockSize = computeLinearSplit(m, div, sliceCount);
 
 		int rem = blockSize - (sliceCount % blockSize);
-		if (rem == blockSize)
+		if ((size_t)rem == blockSize)
 			rem = 0;
 
 		ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -698,12 +840,14 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T
 
 			out.push_back(boost::shared_ptr<CPart>(sub));
 		}
+	} else {
+		out.push_back(boost::shared_ptr<CPart>(clone()));
 	}
 }
 
 void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div)
 {
-	if (true) {
+	if (canSplitAndReduce()) {
 		// Split in vertical direction only at first, until we figure out
 		// a model for splitting in other directions
 
@@ -713,7 +857,7 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T
 		size_t blockSize = computeLinearSplit(m, div, sliceCount);
 
 		int rem = blockSize - (sliceCount % blockSize);
-		if (rem == blockSize)
+		if ((size_t)rem == blockSize)
 			rem = 0;
 
 		ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -747,12 +891,14 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T
 
 			out.push_back(boost::shared_ptr<CPart>(sub));
 		}
+	} else {
+		out.push_back(boost::shared_ptr<CPart>(clone()));
 	}
 }
 
 void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div)
 {
-	if (true) {
+	if (canSplitAndReduce()) {
 		// Split in vertical direction only at first, until we figure out
 		// a model for splitting in other directions
 
@@ -762,7 +908,7 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T
 		size_t blockSize = computeLinearSplit(m, div, sliceCount);
 
 		int rem = blockSize - (sliceCount % blockSize);
-		if (rem == blockSize)
+		if ((size_t)rem == blockSize)
 			rem = 0;
 
 		ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -796,6 +942,8 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T
 
 			out.push_back(boost::shared_ptr<CPart>(sub));
 		}
+	} else {
+		out.push_back(boost::shared_ptr<CPart>(clone()));
 	}
 }
 
@@ -815,7 +963,7 @@ CCompositeGeometryManager::CProjectionPart::~CProjectionPart()
 	delete pGeom;
 }
 
-void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, size_t &z)
+void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, size_t &z) const
 {
 	if (!pGeom) {
 		x = y = z = 0;
@@ -831,6 +979,9 @@ void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, s
 
 CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::reduce(const CPart *_other)
 {
+	if (!canSplitAndReduce())
+		return clone();
+
 	const CVolumePart *other = dynamic_cast<const CVolumePart *>(_other);
 	assert(other);
 
@@ -868,7 +1019,7 @@ CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::re
 
 void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div)
 {
-	if (true) {
+	if (canSplitAndReduce()) {
 		// Split in vertical direction only at first, until we figure out
 		// a model for splitting in other directions
 
@@ -878,7 +1029,7 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage
 		size_t blockSize = computeLinearSplit(m, div, sliceCount);
 
 		int rem = blockSize - (sliceCount % blockSize);
-		if (rem == blockSize)
+		if ((size_t)rem == blockSize)
 			rem = 0;
 
 		ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -903,6 +1054,8 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage
 
 			out.push_back(boost::shared_ptr<CPart>(sub));
 		}
+	} else {
+		out.push_back(boost::shared_ptr<CPart>(clone()));
 	}
 }
 
@@ -914,7 +1067,7 @@ void CCompositeGeometryManager::CProjectionPart::splitY(CCompositeGeometryManage
 
 void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div)
 {
-	if (true) {
+	if (canSplitAndReduce()) {
 		// Split in vertical direction only at first, until we figure out
 		// a model for splitting in other directions
 
@@ -924,7 +1077,7 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage
 		size_t blockSize = computeLinearSplit(m, div, sliceCount);
 
 		int rem = blockSize - (sliceCount % blockSize);
-		if (rem == blockSize)
+		if ((size_t)rem == blockSize)
 			rem = 0;
 
 		ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -949,6 +1102,8 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage
 
 			out.push_back(boost::shared_ptr<CPart>(sub));
 		}
+	} else {
+		out.push_back(boost::shared_ptr<CPart>(clone()));
 	}
 
 }
@@ -959,8 +1114,8 @@ CCompositeGeometryManager::CProjectionPart* CCompositeGeometryManager::CProjecti
 }
 
 CCompositeGeometryManager::SJob CCompositeGeometryManager::createJobFP(CProjector3D *pProjector,
-                                            CFloat32VolumeData3DMemory *pVolData,
-                                            CFloat32ProjectionData3DMemory *pProjData)
+                                            CFloat32VolumeData3D *pVolData,
+                                            CFloat32ProjectionData3D *pProjData)
 {
 	ASTRA_DEBUG("CCompositeGeometryManager::createJobFP");
 	// Create single job for FP
@@ -992,8 +1147,8 @@ CCompositeGeometryManager::SJob CCompositeGeometryManager::createJobFP(CProjecto
 }
 
 CCompositeGeometryManager::SJob CCompositeGeometryManager::createJobBP(CProjector3D *pProjector,
-                                            CFloat32VolumeData3DMemory *pVolData,
-                                            CFloat32ProjectionData3DMemory *pProjData)
+                                            CFloat32VolumeData3D *pVolData,
+                                            CFloat32ProjectionData3D *pProjData)
 {
 	ASTRA_DEBUG("CCompositeGeometryManager::createJobBP");
 	// Create single job for BP
@@ -1022,8 +1177,8 @@ CCompositeGeometryManager::SJob CCompositeGeometryManager::createJobBP(CProjecto
 	return BP;
 }
 
-bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
-                                     CFloat32ProjectionData3DMemory *pProjData)
+bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData,
+                                     CFloat32ProjectionData3D *pProjData)
 {
 	TJobList L;
 	L.push_back(createJobFP(pProjector, pVolData, pProjData));
@@ -1031,8 +1186,8 @@ bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeDat
 	return doJobs(L);
 }
 
-bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
-                                     CFloat32ProjectionData3DMemory *pProjData)
+bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData,
+                                     CFloat32ProjectionData3D *pProjData)
 {
 	TJobList L;
 	L.push_back(createJobBP(pProjector, pVolData, pProjData));
@@ -1041,8 +1196,8 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeDat
 }
 
 
-bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
-                                     CFloat32ProjectionData3DMemory *pProjData, bool bShortScan,
+bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData,
+                                     CFloat32ProjectionData3D *pProjData, bool bShortScan,
                                      const float *pfFilter)
 {
 	if (!dynamic_cast<CConeProjectionGeometry3D*>(pProjData->getGeometry())) {
@@ -1061,11 +1216,11 @@ bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeDa
 	return doJobs(L);
 }
 
-bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData)
+bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3D *>& volData, const std::vector<CFloat32ProjectionData3D *>& projData)
 {
 	ASTRA_DEBUG("CCompositeGeometryManager::doFP, multi-volume");
 
-	std::vector<CFloat32VolumeData3DMemory *>::const_iterator i;
+	std::vector<CFloat32VolumeData3D *>::const_iterator i;
 	std::vector<boost::shared_ptr<CPart> > inputs;
 
 	for (i = volData.begin(); i != volData.end(); ++i) {
@@ -1079,7 +1234,7 @@ bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector
 		inputs.push_back(boost::shared_ptr<CPart>(input));
 	}
 
-	std::vector<CFloat32ProjectionData3DMemory *>::const_iterator j;
+	std::vector<CFloat32ProjectionData3D *>::const_iterator j;
 	std::vector<boost::shared_ptr<CPart> > outputs;
 
 	for (j = projData.begin(); j != projData.end(); ++j) {
@@ -1115,12 +1270,12 @@ bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector
 	return doJobs(L);
 }
 
-bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData)
+bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3D *>& volData, const std::vector<CFloat32ProjectionData3D *>& projData)
 {
 	ASTRA_DEBUG("CCompositeGeometryManager::doBP, multi-volume");
 
 
-	std::vector<CFloat32VolumeData3DMemory *>::const_iterator i;
+	std::vector<CFloat32VolumeData3D *>::const_iterator i;
 	std::vector<boost::shared_ptr<CPart> > outputs;
 
 	for (i = volData.begin(); i != volData.end(); ++i) {
@@ -1134,7 +1289,7 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector
 		outputs.push_back(boost::shared_ptr<CPart>(output));
 	}
 
-	std::vector<CFloat32ProjectionData3DMemory *>::const_iterator j;
+	std::vector<CFloat32ProjectionData3D *>::const_iterator j;
 	std::vector<boost::shared_ptr<CPart> > inputs;
 
 	for (j = projData.begin(); j != projData.end(); ++j) {
@@ -1188,14 +1343,25 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 	if (L.begin()->eType == CCompositeGeometryManager::SJob::JOB_NOP) {
 		// just zero output?
 		if (zero) {
-			for (size_t z = 0; z < outz; ++z) {
-				for (size_t y = 0; y < outy; ++y) {
-					float* ptr = output->pData->getData();
-					ptr += (z + output->subX) * (size_t)output->pData->getHeight() * (size_t)output->pData->getWidth();
-					ptr += (y + output->subY) * (size_t)output->pData->getWidth();
-					ptr += output->subX;
-					memset(ptr, 0, sizeof(float) * outx);
+			// TODO: This function shouldn't have to know about this difference
+			// between Memory/GPU
+			CFloat32Data3DMemory *hostMem = dynamic_cast<CFloat32Data3DMemory *>(output->pData);
+			if (hostMem) {
+				for (size_t z = 0; z < outz; ++z) {
+					for (size_t y = 0; y < outy; ++y) {
+						float* ptr = hostMem->getData();
+						ptr += (z + output->subX) * (size_t)output->pData->getHeight() * (size_t)output->pData->getWidth();
+						ptr += (y + output->subY) * (size_t)output->pData->getWidth();
+						ptr += output->subX;
+						memset(ptr, 0, sizeof(float) * outx);
+					}
 				}
+			} else {
+				CFloat32Data3DGPU *gpuMem = dynamic_cast<CFloat32Data3DGPU *>(output->pData);
+				assert(gpuMem);
+				assert(output->isFull()); // TODO: zero subset?
+
+				zeroGPUMemory(gpuMem->getHandle(), outx, outy, outz);
 			}
 		}
 		return true;
@@ -1214,10 +1380,11 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 	dstdims.subx = output->subX;
 	dstdims.suby = output->subY;
 	dstdims.subz = output->subZ;
-	float *dst = output->pData->getData();
 
-	astraCUDA3d::MemHandle3D outputMem = astraCUDA3d::allocateGPUMemory(outx, outy, outz, zero ? astraCUDA3d::INIT_ZERO : astraCUDA3d::INIT_NO);
-	bool ok = outputMem;
+	CFloat32CustomGPUMemory *dstMem = createGPUMemoryHandler(output->pData);
+
+	bool ok = dstMem->allocateGPUMemory(outx, outy, outz, zero ? astraCUDA3d::INIT_ZERO : astraCUDA3d::INIT_NO);
+	if (!ok) ASTRA_ERROR("Error allocating GPU memory");
 
 	for (CCompositeGeometryManager::TJobList::const_iterator i = L.begin(); i != L.end(); ++i) {
 		const CCompositeGeometryManager::SJob &j = *i;
@@ -1238,7 +1405,8 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 
 		size_t inx, iny, inz;
 		j.pInput->getDims(inx, iny, inz);
-		astraCUDA3d::MemHandle3D inputMem = astraCUDA3d::allocateGPUMemory(inx, iny, inz, astraCUDA3d::INIT_NO);
+
+		CFloat32CustomGPUMemory *srcMem = createGPUMemoryHandler(j.pInput->pData);
 
 		astraCUDA3d::SSubDimensions3D srcdims;
 		srcdims.nx = j.pInput->pData->getWidth();
@@ -1251,9 +1419,11 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 		srcdims.subx = j.pInput->subX;
 		srcdims.suby = j.pInput->subY;
 		srcdims.subz = j.pInput->subZ;
-		const float *src = j.pInput->pData->getDataConst();
 
-		ok = astraCUDA3d::copyToGPUMemory(src, inputMem, srcdims);
+		ok = srcMem->allocateGPUMemory(inx, iny, inz, astraCUDA3d::INIT_NO);
+		if (!ok) ASTRA_ERROR("Error allocating GPU memory");
+
+		ok = srcMem->copyToGPUMemory(srcdims);
 		if (!ok) ASTRA_ERROR("Error copying input data to GPU");
 
 		switch (j.eType) {
@@ -1264,7 +1434,7 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 
 			ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FP");
 
-			ok = astraCUDA3d::FP(((CCompositeGeometryManager::CProjectionPart*)j.pOutput.get())->pGeom, outputMem, ((CCompositeGeometryManager::CVolumePart*)j.pInput.get())->pGeom, inputMem, detectorSuperSampling, projKernel);
+			ok = astraCUDA3d::FP(((CCompositeGeometryManager::CProjectionPart*)j.pOutput.get())->pGeom, dstMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pInput.get())->pGeom, srcMem->hnd, detectorSuperSampling, projKernel);
 			if (!ok) ASTRA_ERROR("Error performing sub-FP");
 			ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FP done");
 		}
@@ -1276,7 +1446,7 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 
 			ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing BP");
 
-			ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, outputMem, voxelSuperSampling, densityWeighting);
+			ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, srcMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, dstMem->hnd, voxelSuperSampling, densityWeighting);
 			if (!ok) ASTRA_ERROR("Error performing sub-BP");
 			ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done");
 		}
@@ -1292,7 +1462,7 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 			} else {
 				ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FDK");
 
-				ok = astraCUDA3d::FDK(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, outputMem, j.FDKSettings.bShortScan, j.FDKSettings.pfFilter);
+				ok = astraCUDA3d::FDK(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, srcMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, dstMem->hnd, j.FDKSettings.bShortScan, j.FDKSettings.pfFilter);
 				if (!ok) ASTRA_ERROR("Error performing sub-FDK");
 				ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FDK done");
 			}
@@ -1302,17 +1472,20 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 			assert(false);
 		}
 
-		ok = astraCUDA3d::freeGPUMemory(inputMem);
+		ok = srcMem->freeGPUMemory();
 		if (!ok) ASTRA_ERROR("Error freeing GPU memory");
 
+		delete srcMem;
 	}
 
-	ok = astraCUDA3d::copyFromGPUMemory(dst, outputMem, dstdims);
+	ok = dstMem->copyFromGPUMemory(dstdims);
 	if (!ok) ASTRA_ERROR("Error copying output data from GPU");
 	
-	ok = astraCUDA3d::freeGPUMemory(outputMem);
+	ok = dstMem->freeGPUMemory();
 	if (!ok) ASTRA_ERROR("Error freeing GPU memory");
 
+	delete dstMem;
+
 	return true;
 }
 
@@ -1455,6 +1628,8 @@ void CCompositeGeometryManager::setGPUIndices(const std::vector<int>& GPUIndices
 
 bool CCompositeGeometryManager::doJobs(TJobList &jobs)
 {
+	// TODO: Proper clean up if substeps fail (Or as proper as possible)
+
 	ASTRA_DEBUG("CCompositeGeometryManager::doJobs");
 
 	// Sort job list into job set by output part
diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp
index 223a9a4..27bb968 100644
--- a/src/CudaBackProjectionAlgorithm3D.cpp
+++ b/src/CudaBackProjectionAlgorithm3D.cpp
@@ -60,8 +60,8 @@ CCudaBackProjectionAlgorithm3D::CCudaBackProjectionAlgorithm3D()
 //----------------------------------------------------------------------------------------
 // Constructor with initialization
 CCudaBackProjectionAlgorithm3D::CCudaBackProjectionAlgorithm3D(CProjector3D* _pProjector, 
-								   CFloat32ProjectionData3DMemory* _pProjectionData, 
-								   CFloat32VolumeData3DMemory* _pReconstruction)
+								   CFloat32ProjectionData3D* _pProjectionData, 
+								   CFloat32VolumeData3D* _pReconstruction)
 {
 	_clear();
 	initialize(_pProjector, _pProjectionData, _pReconstruction);
@@ -145,8 +145,8 @@ bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg)
 //----------------------------------------------------------------------------------------
 // Initialize - C++
 bool CCudaBackProjectionAlgorithm3D::initialize(CProjector3D* _pProjector, 
-								  CFloat32ProjectionData3DMemory* _pSinogram, 
-								  CFloat32VolumeData3DMemory* _pReconstruction)
+								  CFloat32ProjectionData3D* _pSinogram, 
+								  CFloat32VolumeData3D* _pReconstruction)
 {
 	// if already initialized, clear first
 	if (m_bIsInitialized) {
@@ -187,17 +187,21 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
 	// check initialized
 	ASTRA_ASSERT(m_bIsInitialized);
 
-	CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram);
+	CFloat32ProjectionData3D* pSinoMem = dynamic_cast<CFloat32ProjectionData3D*>(m_pSinogram);
 	ASTRA_ASSERT(pSinoMem);
-	CFloat32VolumeData3DMemory* pReconMem = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction);
+	CFloat32VolumeData3D* pReconMem = dynamic_cast<CFloat32VolumeData3D*>(m_pReconstruction);
 	ASTRA_ASSERT(pReconMem);
 
 	const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry();
 	const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry();
 
 	if (m_bSIRTWeighting) {
-		astraCudaBP_SIRTWeighted(pReconMem->getData(),
-		                         pSinoMem->getDataConst(),
+		CFloat32ProjectionData3DMemory* pSinoMemory = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram);
+		ASTRA_ASSERT(pSinoMemory);
+		CFloat32VolumeData3DMemory* pReconMemory = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction);
+		ASTRA_ASSERT(pReconMemory);
+		astraCudaBP_SIRTWeighted(pReconMemory->getData(),
+		                         pSinoMemory->getDataConst(),
 		                         &volgeom, projgeom,
 		                         m_iGPUIndex, m_iVoxelSuperSampling);
 	} else {
diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp
index d02db6d..d503351 100644
--- a/src/CudaFDKAlgorithm3D.cpp
+++ b/src/CudaFDKAlgorithm3D.cpp
@@ -59,8 +59,8 @@ CCudaFDKAlgorithm3D::CCudaFDKAlgorithm3D()
 //----------------------------------------------------------------------------------------
 // Constructor with initialization
 CCudaFDKAlgorithm3D::CCudaFDKAlgorithm3D(CProjector3D* _pProjector, 
-								   CFloat32ProjectionData3DMemory* _pProjectionData, 
-								   CFloat32VolumeData3DMemory* _pReconstruction)
+								   CFloat32ProjectionData3D* _pProjectionData, 
+								   CFloat32VolumeData3D* _pReconstruction)
 {
 	_clear();
 	initialize(_pProjector, _pProjectionData, _pReconstruction);
@@ -179,8 +179,8 @@ bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg)
 //----------------------------------------------------------------------------------------
 // Initialize - C++
 bool CCudaFDKAlgorithm3D::initialize(CProjector3D* _pProjector, 
-								  CFloat32ProjectionData3DMemory* _pSinogram, 
-								  CFloat32VolumeData3DMemory* _pReconstruction)
+								  CFloat32ProjectionData3D* _pSinogram, 
+								  CFloat32VolumeData3D* _pReconstruction)
 {
 	// if already initialized, clear first
 	if (m_bIsInitialized) {
@@ -225,9 +225,9 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations)
 
 	ASTRA_ASSERT(conegeom);
 
-	CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram);
+	CFloat32ProjectionData3D* pSinoMem = dynamic_cast<CFloat32ProjectionData3D*>(m_pSinogram);
 	ASTRA_ASSERT(pSinoMem);
-	CFloat32VolumeData3DMemory* pReconMem = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction);
+	CFloat32VolumeData3D* pReconMem = dynamic_cast<CFloat32VolumeData3D*>(m_pReconstruction);
 	ASTRA_ASSERT(pReconMem);
 
 	const float *filter = NULL;
diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp
index 6783093..ce808eb 100644
--- a/src/CudaForwardProjectionAlgorithm3D.cpp
+++ b/src/CudaForwardProjectionAlgorithm3D.cpp
@@ -101,14 +101,14 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg)
 	node = _cfg.self.getSingleNode("ProjectionDataId");
 	ASTRA_CONFIG_CHECK(node, "CudaForwardProjection3D", "No ProjectionDataId tag specified.");
 	id = node.getContentInt();
-	m_pProjections = dynamic_cast<CFloat32ProjectionData3DMemory*>(CData3DManager::getSingleton().get(id));
+	m_pProjections = dynamic_cast<CFloat32ProjectionData3D*>(CData3DManager::getSingleton().get(id));
 	CC.markNodeParsed("ProjectionDataId");
 
 	// reconstruction data
 	node = _cfg.self.getSingleNode("VolumeDataId");
 	ASTRA_CONFIG_CHECK(node, "CudaForwardProjection3D", "No VolumeDataId tag specified.");
 	id = node.getContentInt();
-	m_pVolume = dynamic_cast<CFloat32VolumeData3DMemory*>(CData3DManager::getSingleton().get(id));
+	m_pVolume = dynamic_cast<CFloat32VolumeData3D*>(CData3DManager::getSingleton().get(id));
 	CC.markNodeParsed("VolumeDataId");
 
 	// optional: projector
@@ -140,8 +140,8 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg)
 
 
 bool CCudaForwardProjectionAlgorithm3D::initialize(CProjector3D* _pProjector, 
-                                  CFloat32ProjectionData3DMemory* _pProjections, 
-                                  CFloat32VolumeData3DMemory* _pVolume,
+                                  CFloat32ProjectionData3D* _pProjections, 
+                                  CFloat32VolumeData3D* _pVolume,
                                   int _iGPUindex, int _iDetectorSuperSampling)
 {
 	m_pProjector = _pProjector;
diff --git a/src/Float32Data3D.cpp b/src/Float32Data3D.cpp
index cad1f18..cc824bd 100644
--- a/src/Float32Data3D.cpp
+++ b/src/Float32Data3D.cpp
@@ -28,6 +28,10 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/Float32Data3D.h"
 #include <sstream>
 
+#ifdef ASTRA_CUDA
+#include "../../cuda/3d/mem3d.h"
+#endif
+
 using namespace std;
 
 namespace astra {
@@ -60,7 +64,5 @@ std::string CFloat32Data3D::description() const
 	if (getType() == CFloat32Data3D::VOLUME) res << " volume data \t";
 	return res.str();
 }
-//----------------------------------------------------------------------------------------
-
 
 } // end namespace astra
diff --git a/src/Float32Data3DGPU.cpp b/src/Float32Data3DGPU.cpp
new file mode 100644
index 0000000..cd9c4ad
--- /dev/null
+++ b/src/Float32Data3DGPU.cpp
@@ -0,0 +1,98 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox 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 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#include "astra/Float32Data3DGPU.h"
+
+namespace astra {
+
+//----------------------------------------------------------------------------------------
+// Default constructor.
+CFloat32Data3DGPU::CFloat32Data3DGPU()
+{
+	_clear();
+	m_bInitialized = false;
+}
+
+//----------------------------------------------------------------------------------------
+// Destructor.
+CFloat32Data3DGPU::~CFloat32Data3DGPU() 
+{
+	if (m_bInitialized)
+	{
+		_unInit();
+	}
+}
+
+//----------------------------------------------------------------------------------------
+// Initializes an instance of the CFloat32Data3DGPU class with pre-allocated memory
+bool CFloat32Data3DGPU::_initialize(int _iWidth, int _iHeight, int _iDepth, astraCUDA3d::MemHandle3D _hnd)
+{
+	// basic checks
+	ASTRA_ASSERT(_iWidth > 0);
+	ASTRA_ASSERT(_iHeight > 0);
+	ASTRA_ASSERT(_iDepth > 0);
+	//ASTRA_ASSERT(_pCustomMemory != NULL);
+
+	if (m_bInitialized) {
+		_unInit();
+	}
+
+	// calculate size
+	m_iWidth = _iWidth;
+	m_iHeight = _iHeight;
+	m_iDepth = _iDepth;
+	m_iSize = (size_t)m_iWidth * m_iHeight * m_iDepth;
+
+	m_hnd = _hnd;
+
+	// initialization complete
+	return true;
+}
+//----------------------------------------------------------------------------------------
+// Clear all member variables, setting all numeric variables to 0 and all pointers to NULL. 
+void CFloat32Data3DGPU::_clear()
+{
+	m_iWidth = 0;
+	m_iHeight = 0;
+	m_iDepth = 0;
+	m_iSize = 0;
+
+	m_hnd.d.reset();
+}
+
+//----------------------------------------------------------------------------------------
+// Un-initialize the object, bringing it back in the unitialized state.
+void CFloat32Data3DGPU::_unInit()
+{
+	ASTRA_ASSERT(m_bInitialized);
+
+	_clear();
+	m_bInitialized = false;
+}
+
+} // end namespace astra
diff --git a/src/Float32Data3DMemory.cpp b/src/Float32Data3DMemory.cpp
index 95b1859..5c5c310 100644
--- a/src/Float32Data3DMemory.cpp
+++ b/src/Float32Data3DMemory.cpp
@@ -27,6 +27,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 
 #include "astra/Float32Data3DMemory.h"
 #include <iostream>
+#include <cstdlib>
 
 namespace astra {
 
@@ -75,10 +76,6 @@ bool CFloat32Data3DMemory::_initialize(int _iWidth, int _iHeight, int _iDepth)
 	m_pCustomMemory = 0;
 	_allocateData();
 
-	// set minmax to default values
-	m_fGlobalMin = 0.0;
-	m_fGlobalMax = 0.0;
-
 	// initialization complete
 	return true;
 
@@ -166,7 +163,7 @@ bool CFloat32Data3DMemory::_initialize(int _iWidth, int _iHeight, int _iDepth, C
 	ASTRA_ASSERT(_iWidth > 0);
 	ASTRA_ASSERT(_iHeight > 0);
 	ASTRA_ASSERT(_iDepth > 0);
-	ASTRA_ASSERT(_pCustomMemory != NULL);
+	//ASTRA_ASSERT(_pCustomMemory != NULL);
 
 	if (m_bInitialized) {
 		_unInit();
@@ -272,9 +269,6 @@ void CFloat32Data3DMemory::_clear()
 	m_ppfDataRowInd = NULL;
 	m_pppfDataSliceInd = NULL;
 	m_pCustomMemory = NULL;
-
-	//m_fGlobalMin = 0.0f;
-	//m_fGlobalMax = 0.0f;
 }
 
 //----------------------------------------------------------------------------------------
@@ -289,37 +283,6 @@ void CFloat32Data3DMemory::_unInit()
 }
 
 //----------------------------------------------------------------------------------------
-// Update data statistics, such as minimum and maximum value, after the data has been modified. 
-void CFloat32Data3DMemory::updateStatistics()
-{
-	_computeGlobalMinMax();
-}
-
-//----------------------------------------------------------------------------------------
-// Find the minimum and maximum data value.
-void CFloat32Data3DMemory::_computeGlobalMinMax() 
-{
-	// basic checks
-	ASTRA_ASSERT(m_bInitialized);
-	ASTRA_ASSERT(m_pfData != NULL);
-	ASTRA_ASSERT(m_iSize > 0);
-	
-	// initial values
-	m_fGlobalMin = m_pfData[0];
-	m_fGlobalMax = m_pfData[0];
-
-	// loop
-	size_t i;
-	float32 v;
-	for (i = 0; i < m_iSize; ++i) 
-	{
-		v = m_pfData[i];
-		if (v < m_fGlobalMin) m_fGlobalMin = v;
-		if (v > m_fGlobalMax) m_fGlobalMax = v;
-	}
-}
-
-//----------------------------------------------------------------------------------------
 // Copy the data block pointed to by _pfData to the data block pointed to by m_pfData.
 void CFloat32Data3DMemory::copyData(const float32* _pfData, size_t _iSize)
 {
diff --git a/src/Float32ProjectionData3D.cpp b/src/Float32ProjectionData3D.cpp
index c751c58..ba41e10 100644
--- a/src/Float32ProjectionData3D.cpp
+++ b/src/Float32ProjectionData3D.cpp
@@ -47,232 +47,6 @@ CFloat32ProjectionData3D::~CFloat32ProjectionData3D() {
 	m_pGeometry = 0;
 }
 
-CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const CFloat32ProjectionData3D& _data)
-{
-	CProjectionGeometry3D * pThisGeometry = getGeometry();
-
-	int iProjectionCount = pThisGeometry->getProjectionCount();
-	int iDetectorCount = pThisGeometry->getDetectorTotCount();
-#ifdef _DEBUG
-	CProjectionGeometry3D * pDataGeometry = _data.getGeometry();
-	int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount();
-
-	ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount());
-	ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount);
-#endif
-
-	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
-		CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
-
-			fThisValue += fDataValue;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnProjection(iProjectionIndex, pThisProjection);
-
-		delete pThisProjection;
-		delete pDataProjection;
-	}
-
-	return *this;
-}
-
-CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const CFloat32ProjectionData3D& _data)
-{
-	CProjectionGeometry3D * pThisGeometry = getGeometry();
-
-	int iProjectionCount = pThisGeometry->getProjectionCount();
-	int iDetectorCount = pThisGeometry->getDetectorTotCount();
-#ifdef _DEBUG
-	CProjectionGeometry3D * pDataGeometry = _data.getGeometry();
-	int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount();
-
-	ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount());
-	ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount);
-#endif
-
-	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
-		CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
-
-			fThisValue -= fDataValue;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnProjection(iProjectionIndex, pThisProjection);
-
-		delete pThisProjection;
-		delete pDataProjection;
-	}
-
-	return *this;
-}
-
-CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const CFloat32ProjectionData3D& _data)
-{
-	CProjectionGeometry3D * pThisGeometry = getGeometry();
-
-	int iProjectionCount = pThisGeometry->getProjectionCount();
-	int iDetectorCount = pThisGeometry->getDetectorTotCount();
-#ifdef _DEBUG
-	CProjectionGeometry3D * pDataGeometry = _data.getGeometry();
-	int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount();
-
-	ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount());
-	ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount);
-#endif
-
-	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
-		CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
-
-			fThisValue *= fDataValue;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnProjection(iProjectionIndex, pThisProjection);
-
-		delete pThisProjection;
-		delete pDataProjection;
-	}
-
-	return *this;
-}
-
-CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const float32& _fScalar)
-{
-	CProjectionGeometry3D * pThisGeometry = getGeometry();
-
-	int iProjectionCount = pThisGeometry->getProjectionCount();
-	int iDetectorCount = pThisGeometry->getDetectorTotCount();
-
-	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-
-			fThisValue *= _fScalar;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnProjection(iProjectionIndex, pThisProjection);
-
-		delete pThisProjection;
-	}
-
-	return *this;
-}
-
-CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator/=(const float32& _fScalar)
-{
-	CProjectionGeometry3D * pThisGeometry = getGeometry();
-
-	int iProjectionCount = pThisGeometry->getProjectionCount();
-	int iDetectorCount = pThisGeometry->getDetectorTotCount();
-
-	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-
-			fThisValue /= _fScalar;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnProjection(iProjectionIndex, pThisProjection);
-
-		delete pThisProjection;
-	}
-
-	return *this;
-}
-
-CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const float32& _fScalar)
-{
-	CProjectionGeometry3D * pThisGeometry = getGeometry();
-
-	int iProjectionCount = pThisGeometry->getProjectionCount();
-	int iDetectorCount = pThisGeometry->getDetectorTotCount();
-
-	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-
-			fThisValue += _fScalar;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnProjection(iProjectionIndex, pThisProjection);
-
-		delete pThisProjection;
-	}
-
-	return *this;
-}
-
-CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const float32& _fScalar)
-{
-	CProjectionGeometry3D * pThisGeometry = getGeometry();
-
-	int iProjectionCount = pThisGeometry->getProjectionCount();
-	int iDetectorCount = pThisGeometry->getDetectorTotCount();
-
-	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-
-			fThisValue -= _fScalar;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnProjection(iProjectionIndex, pThisProjection);
-
-		delete pThisProjection;
-	}
-
-	return *this;
-}
-
 void CFloat32ProjectionData3D::changeGeometry(CProjectionGeometry3D* _pGeometry)
 {
 	if (!m_bInitialized) return;
diff --git a/src/Float32ProjectionData3DGPU.cpp b/src/Float32ProjectionData3DGPU.cpp
new file mode 100644
index 0000000..0e063d6
--- /dev/null
+++ b/src/Float32ProjectionData3DGPU.cpp
@@ -0,0 +1,71 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox 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 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#include "astra/Float32ProjectionData3DGPU.h"
+
+using namespace std;
+
+namespace astra
+{
+
+//----------------------------------------------------------------------------------------
+// Default constructor
+CFloat32ProjectionData3DGPU::CFloat32ProjectionData3DGPU() :
+	CFloat32Data3DGPU() 
+{
+	m_pGeometry = NULL;
+	m_bInitialized = false;
+}
+
+//----------------------------------------------------------------------------------------
+// Create an instance of the CFloat32ProjectionData2D class with pre-allocated data
+CFloat32ProjectionData3DGPU::CFloat32ProjectionData3DGPU(CProjectionGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd)
+{
+	m_bInitialized = false;
+	m_bInitialized = initialize(_pGeometry, _hnd);
+}
+
+
+//----------------------------------------------------------------------------------------
+// Destructor
+CFloat32ProjectionData3DGPU::~CFloat32ProjectionData3DGPU() 
+{
+	delete m_pGeometry;
+	m_pGeometry = 0;
+}
+
+//----------------------------------------------------------------------------------------
+// Initialization
+bool CFloat32ProjectionData3DGPU::initialize(CProjectionGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd)
+{
+	m_pGeometry = _pGeometry->clone();
+	m_bInitialized = _initialize(m_pGeometry->getDetectorColCount(), m_pGeometry->getProjectionCount(), m_pGeometry->getDetectorRowCount(), _hnd);
+	return m_bInitialized;
+}
+
+} // end namespace astra
diff --git a/src/Float32ProjectionData3DMemory.cpp b/src/Float32ProjectionData3DMemory.cpp
index f42e216..69033d5 100644
--- a/src/Float32ProjectionData3DMemory.cpp
+++ b/src/Float32ProjectionData3DMemory.cpp
@@ -114,7 +114,7 @@ bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometr
 
 //----------------------------------------------------------------------------------------
 // Initialization
-bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) 
+bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory)
 {
 	m_pGeometry = _pGeometry->clone();
 	m_bInitialized = _initialize(m_pGeometry->getDetectorColCount(), m_pGeometry->getProjectionCount(), m_pGeometry->getDetectorRowCount(), _pCustomMemory);
@@ -131,102 +131,6 @@ CFloat32ProjectionData3DMemory::~CFloat32ProjectionData3DMemory()
 }
 
 //----------------------------------------------------------------------------------------
-// Fetch a projection
-CFloat32VolumeData2D* CFloat32ProjectionData3DMemory::fetchProjection(int _iProjectionNr) const
-{
-	// fetch slice of the geometry
-	CVolumeGeometry2D volGeom(m_pGeometry->getDetectorColCount(), m_pGeometry->getDetectorRowCount());
-	// create new volume data
-	CFloat32VolumeData2D* res = new CFloat32VolumeData2D(&volGeom);
-	// copy data
-	int row, col;
-	for (row = 0; row < m_pGeometry->getDetectorRowCount(); ++row) {
-		for (col = 0; col < m_pGeometry->getDetectorColCount(); ++col) {
-			res->getData()[row*m_pGeometry->getDetectorColCount() + col] = 
-				m_pfData[_iProjectionNr * m_pGeometry->getDetectorColCount() + m_pGeometry->getDetectorColCount()* m_pGeometry->getProjectionCount() * row + col];
-		}
-	}
-	// return
-	return res;
-}
-
-//----------------------------------------------------------------------------------------
-// Return a projection
-void CFloat32ProjectionData3DMemory::returnProjection(int _iProjectionNr, CFloat32VolumeData2D* _pProjection) 
-{
-	/// TODO: check geometry
-	// copy data
-	int row, col;
-	for (row = 0; row < m_pGeometry->getDetectorRowCount(); ++row) {
-		for (col = 0; col < m_pGeometry->getDetectorColCount(); ++col) {
-			m_pfData[_iProjectionNr * m_pGeometry->getDetectorColCount() + m_pGeometry->getDetectorColCount()* m_pGeometry->getProjectionCount() * row + col] = 
-				_pProjection->getData()[row*m_pGeometry->getDetectorColCount() + col];
-		}
-	}
-}
-
-//----------------------------------------------------------------------------------------
-// Fetch a sinogram
-CFloat32ProjectionData2D* CFloat32ProjectionData3DMemory::fetchSinogram(int _iSliceNr) const
-{
-	CParallelProjectionGeometry3D * pParallelProjGeo = (CParallelProjectionGeometry3D *)m_pGeometry;
-	CParallelProjectionGeometry2D * pProjGeo2D = pParallelProjGeo->createProjectionGeometry2D();
-
-	// create new projection data
-	CFloat32ProjectionData2D* res = new CFloat32ProjectionData2D(pProjGeo2D);
-	// copy data
-	int row, col;
-
-	int iDetectorColumnCount = m_pGeometry->getDetectorColCount();
-	int iProjectionAngleCount = m_pGeometry->getProjectionCount();
-
-	for (row = 0; row < m_pGeometry->getProjectionCount(); ++row) {
-		for (col = 0; col < m_pGeometry->getDetectorColCount(); ++col)
-		{
-			int iTargetIndex = row * iDetectorColumnCount + col;
-			int iSourceIndex = _iSliceNr * iDetectorColumnCount * iProjectionAngleCount + row * iDetectorColumnCount + col;
-
-			float32 fStoredValue = m_pfData[iSourceIndex];
-
-			res->getData()[iTargetIndex] = fStoredValue;
-		}
-	}
-
-	delete pProjGeo2D;
-
-	// return
-	return res;
-}
-
-//----------------------------------------------------------------------------------------
-// Return a sinogram
-void CFloat32ProjectionData3DMemory::returnSinogram(int _iSliceNr, CFloat32ProjectionData2D* _pSinogram2D) 
-{
-	/// TODO: check geometry
-	// copy data
-	int row, col;
-	for (row = 0; row < m_pGeometry->getProjectionCount(); ++row) {
-		for (col = 0; col < m_pGeometry->getDetectorColCount(); ++col) {
-			m_pfData[_iSliceNr*m_pGeometry->getDetectorColCount()*m_pGeometry->getProjectionCount() + row*m_pGeometry->getDetectorColCount() + col] =
-				_pSinogram2D->getData()[row*m_pGeometry->getDetectorColCount() + col];
-		}
-	}
-}
-
-//----------------------------------------------------------------------------------------
-// Returns a specific value
-float32 CFloat32ProjectionData3DMemory::getDetectorValue(int _iIndex)
-{
-	return m_pfData[_iIndex];
-}
-
-//----------------------------------------------------------------------------------------
-// Sets a specific value
-void CFloat32ProjectionData3DMemory::setDetectorValue(int _iIndex, float32 _fValue)
-{
-	m_pfData[_iIndex] = _fValue;
-}
-//----------------------------------------------------------------------------------------
 
 CFloat32ProjectionData3DMemory& CFloat32ProjectionData3DMemory::operator=(const CFloat32ProjectionData3DMemory& _dataIn)
 {
diff --git a/src/Float32VolumeData3D.cpp b/src/Float32VolumeData3D.cpp
index bf7eced..96e6940 100644
--- a/src/Float32VolumeData3D.cpp
+++ b/src/Float32VolumeData3D.cpp
@@ -43,228 +43,6 @@ CFloat32VolumeData3D::~CFloat32VolumeData3D() {
 
 }
 
-CFloat32VolumeData3D& CFloat32VolumeData3D::operator+=(const CFloat32VolumeData3D& _data)
-{
-	CVolumeGeometry3D * pThisGeometry = getGeometry();
-
-	int iSliceCount = pThisGeometry->getGridSliceCount();
-#ifdef _DEBUG
-	CVolumeGeometry3D * pDataGeometry = _data.getGeometry();
-	int iThisSlicePixelCount = pThisGeometry->getGridRowCount() * pThisGeometry->getGridColCount();
-	int iDataSlicePixelCount = pDataGeometry->getGridRowCount() * pDataGeometry->getGridColCount();
-
-	ASTRA_ASSERT(iSliceCount == pDataGeometry->getGridSliceCount());
-	ASTRA_ASSERT(iThisSlicePixelCount == iDataSlicePixelCount);
-#endif
-
-	for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex);
-		CFloat32VolumeData2D * pDataProjection = _data.fetchSliceZ(iSliceIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
-
-			fThisValue += fDataValue;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnSliceZ(iSliceIndex, pThisProjection);
-
-		delete pThisProjection;
-		delete pDataProjection;
-	}
-
-	return *this;
-}
-
-CFloat32VolumeData3D& CFloat32VolumeData3D::operator-=(const CFloat32VolumeData3D& _data)
-{
-	CVolumeGeometry3D * pThisGeometry = getGeometry();
-
-	int iSliceCount = pThisGeometry->getGridSliceCount();
-#ifdef _DEBUG
-	CVolumeGeometry3D * pDataGeometry = _data.getGeometry();
-	int iThisSlicePixelCount = pThisGeometry->getGridRowCount() * pThisGeometry->getGridColCount();
-	int iDataSlicePixelCount = pDataGeometry->getGridRowCount() * pDataGeometry->getGridColCount();
-
-	ASTRA_ASSERT(iSliceCount == pDataGeometry->getGridSliceCount());
-	ASTRA_ASSERT(iThisSlicePixelCount == iDataSlicePixelCount);
-#endif
-
-	for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex);
-		CFloat32VolumeData2D * pDataProjection = _data.fetchSliceZ(iSliceIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
-
-			fThisValue -= fDataValue;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnSliceZ(iSliceIndex, pThisProjection);
-
-		delete pThisProjection;
-		delete pDataProjection;
-	}
-
-	return *this;
-}
-
-CFloat32VolumeData3D& CFloat32VolumeData3D::operator*=(const CFloat32VolumeData3D& _data)
-{
-	CVolumeGeometry3D * pThisGeometry = getGeometry();
-
-	int iSliceCount = pThisGeometry->getGridSliceCount();
-#ifdef _DEBUG
-	CVolumeGeometry3D * pDataGeometry = _data.getGeometry();
-	int iThisSlicePixelCount = pThisGeometry->getGridRowCount() * pThisGeometry->getGridColCount();
-	int iDataSlicePixelCount = pDataGeometry->getGridRowCount() * pDataGeometry->getGridColCount();
-
-	ASTRA_ASSERT(iSliceCount == pDataGeometry->getGridSliceCount());
-	ASTRA_ASSERT(iThisSlicePixelCount == iDataSlicePixelCount);
-#endif
-
-	for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex);
-		CFloat32VolumeData2D * pDataProjection = _data.fetchSliceZ(iSliceIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
-
-			fThisValue *= fDataValue;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnSliceZ(iSliceIndex, pThisProjection);
-
-		delete pThisProjection;
-		delete pDataProjection;
-	}
-
-	return *this;
-}
-
-CFloat32VolumeData3D& CFloat32VolumeData3D::operator*=(const float32& _fScalar)
-{
-	CVolumeGeometry3D * pThisGeometry = getGeometry();
-
-	int iSliceCount = pThisGeometry->getGridSliceCount();
-
-	for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-
-			fThisValue *= _fScalar;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnSliceZ(iSliceIndex, pThisProjection);
-
-		delete pThisProjection;
-	}
-
-	return *this;
-}
-
-CFloat32VolumeData3D& CFloat32VolumeData3D::operator/=(const float32& _fScalar)
-{
-	CVolumeGeometry3D * pThisGeometry = getGeometry();
-
-	int iSliceCount = pThisGeometry->getGridSliceCount();
-
-	for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-
-			fThisValue /= _fScalar;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnSliceZ(iSliceIndex, pThisProjection);
-
-		delete pThisProjection;
-	}
-
-	return *this;
-}
-
-CFloat32VolumeData3D& CFloat32VolumeData3D::operator+=(const float32& _fScalar)
-{
-	CVolumeGeometry3D * pThisGeometry = getGeometry();
-
-	int iSliceCount = pThisGeometry->getGridSliceCount();
-
-	for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-
-			fThisValue += _fScalar;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnSliceZ(iSliceIndex, pThisProjection);
-
-		delete pThisProjection;
-	}
-
-	return *this;
-}
-
-CFloat32VolumeData3D& CFloat32VolumeData3D::operator-=(const float32& _fScalar)
-{
-	CVolumeGeometry3D * pThisGeometry = getGeometry();
-
-	int iSliceCount = pThisGeometry->getGridSliceCount();
-
-	for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++)
-	{
-		CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex);
-
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
-		{
-			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
-
-			fThisValue -= _fScalar;
-
-			pThisProjection->getData()[iDetectorIndex] = fThisValue;
-		}
-
-		returnSliceZ(iSliceIndex, pThisProjection);
-
-		delete pThisProjection;
-	}
-
-	return *this;
-}
-
 void CFloat32VolumeData3D::changeGeometry(CVolumeGeometry3D* _pGeometry)
 {
 	if (!m_bInitialized) return;
diff --git a/src/Float32VolumeData3DGPU.cpp b/src/Float32VolumeData3DGPU.cpp
new file mode 100644
index 0000000..82a222f
--- /dev/null
+++ b/src/Float32VolumeData3DGPU.cpp
@@ -0,0 +1,71 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox 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 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#include "astra/Float32VolumeData3DGPU.h"
+
+using namespace std;
+
+namespace astra
+{
+
+//----------------------------------------------------------------------------------------
+// Default constructor
+CFloat32VolumeData3DGPU::CFloat32VolumeData3DGPU() :
+	CFloat32Data3DGPU() 
+{
+	m_pGeometry = NULL;
+	m_bInitialized = false;
+}
+
+//----------------------------------------------------------------------------------------
+// Create an instance of the CFloat32VolumeData2D class with pre-allocated data
+CFloat32VolumeData3DGPU::CFloat32VolumeData3DGPU(CVolumeGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd)
+{
+	m_bInitialized = false;
+	m_bInitialized = initialize(_pGeometry, _hnd);
+}
+
+
+//----------------------------------------------------------------------------------------
+// Destructor
+CFloat32VolumeData3DGPU::~CFloat32VolumeData3DGPU() 
+{
+	delete m_pGeometry;
+	m_pGeometry = 0;
+}
+
+//----------------------------------------------------------------------------------------
+// Initialization
+bool CFloat32VolumeData3DGPU::initialize(CVolumeGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd)
+{
+	m_pGeometry = _pGeometry->clone();
+	m_bInitialized = _initialize(m_pGeometry->getGridColCount(), m_pGeometry->getGridRowCount(), m_pGeometry->getGridSliceCount(), _hnd);
+	return m_bInitialized;
+}
+
+} // end namespace astra
diff --git a/src/Float32VolumeData3DMemory.cpp b/src/Float32VolumeData3DMemory.cpp
index ba39eaa..27ae87b 100644
--- a/src/Float32VolumeData3DMemory.cpp
+++ b/src/Float32VolumeData3DMemory.cpp
@@ -114,103 +114,13 @@ bool CFloat32VolumeData3DMemory::initialize(CVolumeGeometry3D* _pGeometry, float
 }
 //----------------------------------------------------------------------------------------
 // Initialization
-bool CFloat32VolumeData3DMemory::initialize(CVolumeGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) 
+bool CFloat32VolumeData3DMemory::initialize(CVolumeGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory)
 {
 	m_pGeometry = _pGeometry->clone();
 	m_bInitialized = _initialize(m_pGeometry->getGridColCount(), m_pGeometry->getGridRowCount(), m_pGeometry->getGridSliceCount(), _pCustomMemory);
 	return m_bInitialized;
 }
 
-
-//----------------------------------------------------------------------------------------
-// Fetch a slice
-CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceZ(int _iSliceIndex) const
-{
-	// fetch slice of the geometry
-	int iRowCount = m_pGeometry->getGridRowCount();
-	int iColumnCount = m_pGeometry->getGridColCount();
-	CVolumeGeometry2D volGeom(iColumnCount, iRowCount);
-
-	// create new volume data
-	CFloat32VolumeData2D* res = new CFloat32VolumeData2D(&volGeom);
-
-	// copy data
-	float * pfTargetData = res->getData();
-	for(int iRowIndex = 0; iRowIndex < iRowCount; iRowIndex++)
-	{
-		for(int iColumnIndex = 0; iColumnIndex < iColumnCount; iColumnIndex++)
-		{
-			int iSourceIndex = _iSliceIndex * iColumnCount * iRowCount + iRowIndex * iColumnCount + iColumnIndex;
-			int iTargetIndex = iRowIndex * iColumnCount + iColumnIndex;
-			float fStoredValue = m_pfData[iSourceIndex];
-			pfTargetData[iTargetIndex] = fStoredValue;
-		}
-	}
-	// return
-	return res;
-}
-
-//----------------------------------------------------------------------------------------
-// Return a slice
-void CFloat32VolumeData3DMemory::returnSliceZ(int _iSliceIndex, CFloat32VolumeData2D * _pSlice)
-{
-	int iRowCount = _pSlice->getGeometry()->getGridRowCount();
-	int iColumnCount = _pSlice->getGeometry()->getGridColCount();
-
-	assert(iRowCount == m_pGeometry->getGridRowCount());
-	assert(iColumnCount == m_pGeometry->getGridColCount());
-
-	for(int iRowIndex = 0; iRowIndex < iRowCount; iRowIndex++)
-	{
-		for(int iColumnIndex = 0; iColumnIndex < iColumnCount; iColumnIndex++)
-		{
-			int iSourceIndex = iRowIndex * iColumnCount + iColumnIndex;
-			int iTargetIndex = _iSliceIndex * iColumnCount * iRowCount + iRowIndex * iColumnCount + iColumnIndex;
-			float fStoredValue = _pSlice->getDataConst()[iSourceIndex];
-			m_pfData[iTargetIndex] = fStoredValue;
-		}
-	}
-}
-
-CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceX(int _iColumnIndex) const
-{
-	// TODO:
-	assert(false);
-	return NULL;
-}
-
-CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceY(int _iRowIndex) const
-{
-	// TODO:
-	assert(false);
-	return NULL;
-}
-
-void CFloat32VolumeData3DMemory::returnSliceX(int _iColumnIndex, CFloat32VolumeData2D * _pSliceData)
-{
-	// TODO:
-	assert(false);
-}
-
-void CFloat32VolumeData3DMemory::returnSliceY(int _iRowIndex, CFloat32VolumeData2D * _pSliceData)
-{
-	// TODO:
-	assert(false);
-}
-
-//----------------------------------------------------------------------------------------
-// Returns a specific value
-float32 CFloat32VolumeData3DMemory::getVoxelValue(int _iIndex)
-{
-	return m_pfData[_iIndex];
-}
-
-//----------------------------------------------------------------------------------------
-// Sets a specific value
-void CFloat32VolumeData3DMemory::setVoxelValue(int _iIndex, float32 _fValue)
-{
-	m_pfData[_iIndex] = _fValue;
-}
 //----------------------------------------------------------------------------------------
 
 CFloat32VolumeData3DMemory& CFloat32VolumeData3DMemory::operator=(const CFloat32VolumeData3DMemory& _dataIn)