Merge pull request #93 from wjp/GPULink

GPULink support

8 files changed, 358 insertions, 279 deletions

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