From 495512324b84a75782f9fbc11921668ad9c170a9 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 5 Mar 2015 15:19:43 +0100 Subject: Added Python support for CUDA projectors --- python/astra/ASTRAProjector.py | 9 +-- python/astra/PyIncludes.pxd | 16 +++++ python/astra/PyProjector3DFactory.pxd | 35 ++++++++++ python/astra/PyProjector3DManager.pxd | 39 +++++++++++ python/astra/__init__.py | 1 + python/astra/creators.py | 66 +++++++------------ python/astra/projector.py | 30 ++++++--- python/astra/projector3d.py | 100 ++++++++++++++++++++++++++++ python/astra/projector3d_c.pyx | 119 ++++++++++++++++++++++++++++++++++ python/astra/projector_c.pyx | 17 +++++ 10 files changed, 375 insertions(+), 57 deletions(-) create mode 100644 python/astra/PyProjector3DFactory.pxd create mode 100644 python/astra/PyProjector3DManager.pxd create mode 100644 python/astra/projector3d.py create mode 100644 python/astra/projector3d_c.pyx diff --git a/python/astra/ASTRAProjector.py b/python/astra/ASTRAProjector.py index 96acb10..f282618 100644 --- a/python/astra/ASTRAProjector.py +++ b/python/astra/ASTRAProjector.py @@ -70,11 +70,9 @@ class ASTRAProjector2D(object): :type vol_geom: :class:`dict` :param proj_type: Projector type, such as ``'line'``, ``'linear'``, ... :type proj_type: :class:`string` - :param useCUDA: If ``True``, use CUDA for calculations, when possible. - :type useCUDA: :class:`bool` """ - def __init__(self, proj_geom, vol_geom, proj_type, useCUDA=False): + def __init__(self, proj_geom, vol_geom, proj_type): self.vol_geom = vol_geom self.recSize = vol_geom['GridColCount'] self.angles = proj_geom['ProjectionAngles'] @@ -84,7 +82,6 @@ class ASTRAProjector2D(object): self.nProj = self.angles.shape[0] self.proj_geom = proj_geom self.proj_id = ac.create_projector(proj_type, proj_geom, vol_geom) - self.useCUDA = useCUDA self.T = ASTRAProjector2DTranspose(self) def backProject(self, data): @@ -96,7 +93,7 @@ class ASTRAProjector2D(object): """ vol_id, vol = ac.create_backprojection( - data, self.proj_id, useCUDA=self.useCUDA, returnData=True) + data, self.proj_id, returnData=True) data2d.delete(vol_id) return vol @@ -108,7 +105,7 @@ class ASTRAProjector2D(object): :returns: :class:`numpy.ndarray` -- The forward projection. """ - sin_id, sino = ac.create_sino(data, self.proj_id, useCUDA=self.useCUDA, returnData=True) + sin_id, sino = ac.create_sino(data, self.proj_id, returnData=True) data2d.delete(sin_id) return sino diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 434546a..7df02c5 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -27,6 +27,8 @@ from libcpp cimport bool from libcpp.string cimport string from .PyXMLDocument cimport XMLNode +include "config.pxi" + cdef extern from "astra/Globals.h" namespace "astra": ctypedef float float32 ctypedef double float64 @@ -150,6 +152,20 @@ cdef extern from "astra/Projector2D.h" namespace "astra": CVolumeGeometry2D* getVolumeGeometry() CSparseMatrix* getMatrix() +cdef extern from "astra/Projector3D.h" namespace "astra": + 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/CudaProjector2D.h" namespace "astra": + cdef cppclass CCudaProjector2D + + cdef extern from "astra/SparseMatrix.h" namespace "astra": cdef cppclass CSparseMatrix: CSparseMatrix(unsigned int,unsigned int,unsigned long) diff --git a/python/astra/PyProjector3DFactory.pxd b/python/astra/PyProjector3DFactory.pxd new file mode 100644 index 0000000..bcbce94 --- /dev/null +++ b/python/astra/PyProjector3DFactory.pxd @@ -0,0 +1,35 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to 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 Python interface to 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 Python interface to the ASTRA Toolbox. If not, see . +# +#----------------------------------------------------------------------- +from libcpp.string cimport string +from libcpp cimport bool +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectFactory.h" namespace "astra": + cdef cppclass CProjector3DFactory: + CProjector3D *create(Config) + +cdef extern from "astra/AstraObjectFactory.h" namespace "astra::CProjector3DFactory": + cdef CProjector3DFactory* getSingletonPtr() diff --git a/python/astra/PyProjector3DManager.pxd b/python/astra/PyProjector3DManager.pxd new file mode 100644 index 0000000..b1eac6b --- /dev/null +++ b/python/astra/PyProjector3DManager.pxd @@ -0,0 +1,39 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to 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 Python interface to 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 Python interface to the ASTRA Toolbox. If not, see . +# +#----------------------------------------------------------------------- +from libcpp.string cimport string + +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectManager.h" namespace "astra": + cdef cppclass CProjector3DManager: + string info() + void clear() + void remove(int i) + int store(CProjector3D *) + CProjector3D * get(int i) + +cdef extern from "astra/AstraObjectManager.h" namespace "astra::CProjector3DManager": + cdef CProjector3DManager* getSingletonPtr() diff --git a/python/astra/__init__.py b/python/astra/__init__.py index a61aafc..c249c01 100644 --- a/python/astra/__init__.py +++ b/python/astra/__init__.py @@ -33,6 +33,7 @@ from . import astra from . import data3d from . import algorithm from . import projector +from . import projector3d from . import matrix import os diff --git a/python/astra/creators.py b/python/astra/creators.py index 9aba464..2e2dc71 100644 --- a/python/astra/creators.py +++ b/python/astra/creators.py @@ -30,15 +30,16 @@ import math from . import data2d from . import data3d from . import projector +from . import projector3d from . import algorithm def astra_dict(intype): """Creates a dict to use with the ASTRA Toolbox. - + :param intype: Type of the ASTRA object. :type intype: :class:`string` :returns: :class:`dict` -- An ASTRA dict of type ``intype``. - + """ if intype == 'SIRT_CUDA2': intype = 'SIRT_CUDA' @@ -255,25 +256,23 @@ This method can be called in a number of ways: raise Exception('not enough variables: astra_create_proj_geom(parallel3d_vec, det_row_count, det_col_count, V)') if not args[2].shape[1] == 12: raise Exception('V should be a Nx12 matrix, with N the number of projections') - return {'type': 'parallel3d_vec','DetectorRowCount':args[0],'DetectorColCount':args[1],'Vectors':args[2]} + return {'type': 'parallel3d_vec','DetectorRowCount':args[0],'DetectorColCount':args[1],'Vectors':args[2]} elif intype == 'sparse_matrix': if len(args) < 4: raise Exception( 'not enough variables: astra_create_proj_geom(sparse_matrix, det_width, det_count, angles, matrix_id)') return {'type': 'sparse_matrix', 'DetectorWidth': args[0], 'DetectorCount': args[1], 'ProjectionAngles': args[2], 'MatrixID': args[3]} else: - raise Exception('Error: unknown type ' + intype) + raise Exception('Error: unknown type ' + intype) -def create_backprojection(data, proj_id, useCUDA=False, returnData=True): +def create_backprojection(data, proj_id, returnData=True): """Create a backprojection of a sinogram (2D). :param data: Sinogram data or ID. :type data: :class:`numpy.ndarray` or :class:`int` :param proj_id: ID of the projector to use. :type proj_id: :class:`int` -:param useCUDA: If ``True``, use CUDA for the calculation. -:type useCUDA: :class:`bool` :param returnData: If False, only return the ID of the backprojection. :type returnData: :class:`bool` :returns: :class:`int` or (:class:`int`, :class:`numpy.ndarray`) -- If ``returnData=False``, returns the ID of the backprojection. Otherwise, returns a tuple containing the ID of the backprojection and the backprojection itself, in that order. @@ -287,13 +286,13 @@ def create_backprojection(data, proj_id, useCUDA=False, returnData=True): sino_id = data vol_id = data2d.create('-vol', vol_geom, 0) - algString = 'BP' - if useCUDA: - algString = algString + '_CUDA' + if projector.is_cuda(proj_id): + algString = 'BP_CUDA' + else: + algString = 'BP' cfg = astra_dict(algString) - if not useCUDA: - cfg['ProjectorId'] = proj_id + cfg['ProjectorId'] = proj_id cfg['ProjectionDataId'] = sino_id cfg['ReconstructionDataId'] = vol_id alg_id = algorithm.create(cfg) @@ -345,20 +344,13 @@ def create_backprojection3d_gpu(data, proj_geom, vol_geom, returnData=True): return vol_id -def create_sino(data, proj_id=None, proj_geom=None, vol_geom=None, - useCUDA=False, returnData=True, gpuIndex=None): +def create_sino(data, proj_id, returnData=True, gpuIndex=None): """Create a forward projection of an image (2D). :param data: Image data or ID. :type data: :class:`numpy.ndarray` or :class:`int` :param proj_id: ID of the projector to use. :type proj_id: :class:`int` - :param proj_geom: Projection geometry. - :type proj_geom: :class:`dict` - :param vol_geom: Volume geometry. - :type vol_geom: :class:`dict` - :param useCUDA: If ``True``, use CUDA for the calculation. - :type useCUDA: :class:`bool` :param returnData: If False, only return the ID of the forward projection. :type returnData: :class:`bool` :param gpuIndex: Optional GPU index. @@ -374,31 +366,20 @@ def create_sino(data, proj_id=None, proj_geom=None, vol_geom=None, ``proj_geom`` and ``vol_geom``. If ``proj_id`` is given, then ``proj_geom`` and ``vol_geom`` must be None and vice versa. """ - if proj_id is not None: - proj_geom = projector.projection_geometry(proj_id) - vol_geom = projector.volume_geometry(proj_id) - elif proj_geom is not None and vol_geom is not None: - if not useCUDA: - # We need more parameters to create projector. - raise ValueError( - """A ``proj_id`` is needed when CUDA is not used.""") - else: - raise Exception("""The geometry setup is not defined. - The geometry of setup is defined by ``proj_id`` or with - ``proj_geom`` and ``vol_geom``. If ``proj_id`` is given, then - ``proj_geom`` and ``vol_geom`` must be None and vice versa.""") + proj_geom = projector.projection_geometry(proj_id) + vol_geom = projector.volume_geometry(proj_id) if isinstance(data, np.ndarray): volume_id = data2d.create('-vol', vol_geom, data) else: volume_id = data sino_id = data2d.create('-sino', proj_geom, 0) - algString = 'FP' - if useCUDA: - algString = algString + '_CUDA' + if projector.is_cuda(proj_id): + algString = 'FP_CUDA' + else: + algString = 'FP' cfg = astra_dict(algString) - if not useCUDA: - cfg['ProjectorId'] = proj_id + cfg['ProjectorId'] = proj_id if gpuIndex is not None: cfg['option'] = {'GPUindex': gpuIndex} cfg['ProjectionDataId'] = sino_id @@ -496,8 +477,7 @@ def create_reconstruction(rec_type, proj_id, sinogram, iterations=1, use_mask='n vol_geom = projector.volume_geometry(proj_id) recon_id = data2d.create('-vol', vol_geom, 0) cfg = astra_dict(rec_type) - if not 'CUDA' in rec_type: - cfg['ProjectorId'] = proj_id + cfg['ProjectorId'] = proj_id cfg['ProjectionDataId'] = sino_id cfg['ReconstructionDataId'] = recon_id cfg['options'] = {} @@ -560,4 +540,8 @@ def create_projector(proj_type, proj_geom, vol_geom): cfg = astra_dict(proj_type) cfg['ProjectionGeometry'] = proj_geom cfg['VolumeGeometry'] = vol_geom - return projector.create(cfg) + types3d = ['linear3d', 'linearcone', 'cuda3d'] + if proj_type in types3d: + return projector3d.create(cfg) + else: + return projector.create(cfg) diff --git a/python/astra/projector.py b/python/astra/projector.py index c916c52..e370e5a 100644 --- a/python/astra/projector.py +++ b/python/astra/projector.py @@ -27,21 +27,21 @@ from . import projector_c as p def create(config): """Create projector object. - + :param config: Projector options. :type config: :class:`dict` :returns: :class:`int` -- the ID of the constructed object. - + """ return p.create(config) def delete(ids): """Delete a projector object. - + :param ids: ID or list of ID's to delete. :type ids: :class:`int` or :class:`list` - + """ return p.delete(ids) @@ -57,22 +57,22 @@ def info(): def projection_geometry(i): """Get projection geometry of a projector. - + :param i: ID of projector. :type i: :class:`int` :returns: :class:`dict` -- projection geometry - + """ return p.projection_geometry(i) def volume_geometry(i): """Get volume geometry of a projector. - + :param i: ID of projector. :type i: :class:`int` :returns: :class:`dict` -- volume geometry - + """ return p.volume_geometry(i) @@ -88,13 +88,23 @@ def weights_projection(i, projection_index): def splat(i, row, col): return p.splat(i, row, col) +def is_cuda(i): + """Check whether a projector is a CUDA projector. + + :param i: ID of projector. + :type i: :class:`int` + :returns: :class:`bool` -- True if the projector is a CUDA projector. + + """ + return p.is_cuda(i) + def matrix(i): """Get sparse matrix of a projector. - + :param i: ID of projector. :type i: :class:`int` :returns: :class:`int` -- ID of sparse matrix. - + """ return p.matrix(i) diff --git a/python/astra/projector3d.py b/python/astra/projector3d.py new file mode 100644 index 0000000..d1086b9 --- /dev/null +++ b/python/astra/projector3d.py @@ -0,0 +1,100 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to 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 Python interface to 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 Python interface to the ASTRA Toolbox. If not, see . +# +#----------------------------------------------------------------------- +from . import projector3d_c as p + +def create(config): + """Create projector object. + + :param config: Projector options. + :type config: :class:`dict` + :returns: :class:`int` -- the ID of the constructed object. + + """ + return p.create(config) + + +def delete(ids): + """Delete a projector object. + + :param ids: ID or list of ID's to delete. + :type ids: :class:`int` or :class:`list` + + """ + return p.delete(ids) + + +def clear(): + """Clear all projector objects.""" + return p.clear() + + +def info(): + """Print info on projector objects in memory.""" + return p.info() + +def projection_geometry(i): + """Get projection geometry of a projector. + + :param i: ID of projector. + :type i: :class:`int` + :returns: :class:`dict` -- projection geometry + + """ + return p.projection_geometry(i) + + +def volume_geometry(i): + """Get volume geometry of a projector. + + :param i: ID of projector. + :type i: :class:`int` + :returns: :class:`dict` -- volume geometry + + """ + return p.volume_geometry(i) + + +def weights_single_ray(i, projection_index, detector_index): + return p.weights_single_ray(i, projection_index, detector_index) + + +def weights_projection(i, projection_index): + return p.weights_projection(i, projection_index) + + +def splat(i, row, col): + return p.splat(i, row, col) + + +def is_cuda(i): + """Check whether a projector is a CUDA projector. + + :param i: ID of projector. + :type i: :class:`int` + :returns: :class:`bool` -- True if the projector is a CUDA projector. + + """ + return p.is_cuda(i) diff --git a/python/astra/projector3d_c.pyx b/python/astra/projector3d_c.pyx new file mode 100644 index 0000000..8b978d7 --- /dev/null +++ b/python/astra/projector3d_c.pyx @@ -0,0 +1,119 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to 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 Python interface to 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 Python interface to the ASTRA Toolbox. If not, see . +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +import six +from .PyIncludes cimport * + +cimport utils +from .utils import wrap_from_bytes + +cimport PyProjector3DFactory +from .PyProjector3DFactory cimport CProjector3DFactory + +cimport PyProjector3DManager +from .PyProjector3DManager cimport CProjector3DManager + +cimport PyXMLDocument +from .PyXMLDocument cimport XMLDocument + +cdef CProjector3DManager * manProj = PyProjector3DManager.getSingletonPtr() + +include "config.pxi" + +IF HAVE_CUDA: + cdef extern from *: + CCudaProjector3D* dynamic_cast_cuda_projector "dynamic_cast" (CProjector3D*) + + +def create(config): + cdef Config * cfg = utils.dictToConfig(six.b('Projector3D'), config) + cdef CProjector3D * proj + proj = PyProjector3DFactory.getSingletonPtr().create(cfg[0]) + if proj == NULL: + del cfg + raise Exception("Error creating Projector3D.") + del cfg + return manProj.store(proj) + + +def delete(ids): + try: + for i in ids: + manProj.remove(i) + except TypeError: + manProj.remove(ids) + + +def clear(): + manProj.clear() + + +def info(): + six.print_(wrap_from_bytes(manProj.info())) + +cdef CProjector3D * getObject(i) except NULL: + cdef CProjector3D * proj = manProj.get(i) + if proj == NULL: + raise Exception("Projector not initialized.") + if not proj.isInitialized(): + raise Exception("Projector not initialized.") + return proj + + +def projection_geometry(i): + cdef CProjector3D * proj = getObject(i) + return utils.configToDict(proj.getProjectionGeometry().getConfiguration()) + + +def volume_geometry(i): + cdef CProjector3D * proj = getObject(i) + return utils.configToDict(proj.getVolumeGeometry().getConfiguration()) + + +def weights_single_ray(i, projection_index, detector_index): + raise Exception("Not yet implemented") + + +def weights_projection(i, projection_index): + raise Exception("Not yet implemented") + + +def splat(i, row, col): + raise Exception("Not yet implemented") + +def is_cuda(i): + cdef CProjector3D * proj = getObject(i) + IF HAVE_CUDA==True: + cdef CCudaProjector3D * cudaproj = NULL + cudaproj = dynamic_cast_cuda_projector(proj) + if cudaproj==NULL: + return False + else: + return True + ELSE: + return False diff --git a/python/astra/projector_c.pyx b/python/astra/projector_c.pyx index f91a8dd..9aa868e 100644 --- a/python/astra/projector_c.pyx +++ b/python/astra/projector_c.pyx @@ -47,6 +47,12 @@ from .PyMatrixManager cimport CMatrixManager cdef CProjector2DManager * manProj = PyProjector2DManager.getSingletonPtr() cdef CMatrixManager * manM = PyMatrixManager.getSingletonPtr() +include "config.pxi" + +IF HAVE_CUDA: + cdef extern from *: + CCudaProjector2D* dynamic_cast_cuda_projector "dynamic_cast" (CProjector2D*) + def create(config): cdef Config * cfg = utils.dictToConfig(six.b('Projector2D'), config) @@ -104,6 +110,17 @@ def weights_projection(i, projection_index): def splat(i, row, col): raise Exception("Not yet implemented") +def is_cuda(i): + cdef CProjector2D * proj = getObject(i) + IF HAVE_CUDA==True: + cdef CCudaProjector2D * cudaproj = NULL + cudaproj = dynamic_cast_cuda_projector(proj) + if cudaproj==NULL: + return False + else: + return True + ELSE: + return False def matrix(i): cdef CProjector2D * proj = getObject(i) -- cgit v1.2.3