diff options
Diffstat (limited to 'python')
26 files changed, 1163 insertions, 284 deletions
diff --git a/python/astra/ASTRAProjector.py b/python/astra/ASTRAProjector.py deleted file mode 100644 index 96acb10..0000000 --- a/python/astra/ASTRAProjector.py +++ /dev/null @@ -1,138 +0,0 @@ -#----------------------------------------------------------------------- -#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 <http://www.gnu.org/licenses/>. -# -#----------------------------------------------------------------------- - -import math -from . import creators as ac -from . import data2d - - -class ASTRAProjector2DTranspose(): - """Implements the ``proj.T`` functionality. - - Do not use directly, since it can be accessed as member ``.T`` of - an :class:`ASTRAProjector2D` object. - - """ - def __init__(self, parentProj): - self.parentProj = parentProj - - def __mul__(self, data): - return self.parentProj.backProject(data) - - -class ASTRAProjector2D(object): - """Helps with various common ASTRA Toolbox 2D operations. - - This class can perform several often used toolbox operations, such as: - - * Forward projecting - * Back projecting - * Reconstructing - - Note that this class has a some computational overhead, because it - copies a lot of data. If you use many repeated operations, directly - using the PyAstraToolbox methods directly is faster. - - You can use this class as an abstracted weight matrix :math:`W`: multiplying an instance - ``proj`` of this class by an image results in a forward projection of the image, and multiplying - ``proj.T`` by a sinogram results in a backprojection of the sinogram:: - - proj = ASTRAProjector2D(...) - fp = proj*image - bp = proj.T*sinogram - - :param proj_geom: The projection geometry. - :type proj_geom: :class:`dict` - :param vol_geom: The volume geometry. - :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): - self.vol_geom = vol_geom - self.recSize = vol_geom['GridColCount'] - self.angles = proj_geom['ProjectionAngles'] - self.nDet = proj_geom['DetectorCount'] - nexpow = int(pow(2, math.ceil(math.log(2 * self.nDet, 2)))) - self.filterSize = nexpow / 2 + 1 - 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): - """Backproject a sinogram. - - :param data: The sinogram data or ID. - :type data: :class:`numpy.ndarray` or :class:`int` - :returns: :class:`numpy.ndarray` -- The backprojection. - - """ - vol_id, vol = ac.create_backprojection( - data, self.proj_id, useCUDA=self.useCUDA, returnData=True) - data2d.delete(vol_id) - return vol - - def forwardProject(self, data): - """Forward project an image. - - :param data: The image data or ID. - :type data: :class:`numpy.ndarray` or :class:`int` - :returns: :class:`numpy.ndarray` -- The forward projection. - - """ - sin_id, sino = ac.create_sino(data, self.proj_id, useCUDA=self.useCUDA, returnData=True) - data2d.delete(sin_id) - return sino - - def reconstruct(self, data, method, **kwargs): - """Reconstruct an image from a sinogram. - - :param data: The sinogram data or ID. - :type data: :class:`numpy.ndarray` or :class:`int` - :param method: Name of the reconstruction algorithm. - :type method: :class:`string` - :param kwargs: Additional named parameters to pass to :func:`astra.creators.create_reconstruction`. - :returns: :class:`numpy.ndarray` -- The reconstruction. - - Example of a SIRT reconstruction using CUDA:: - - proj = ASTRAProjector2D(...) - rec = proj.reconstruct(sinogram,'SIRT_CUDA',iterations=1000) - - """ - kwargs['returnData'] = True - rec_id, rec = ac.create_reconstruction( - method, self.proj_id, data, **kwargs) - data2d.delete(rec_id) - return rec - - def __mul__(self, data): - return self.forwardProject(data) diff --git a/python/astra/CFloat32CustomPython.h b/python/astra/CFloat32CustomPython.h new file mode 100644 index 0000000..d8593fc --- /dev/null +++ b/python/astra/CFloat32CustomPython.h @@ -0,0 +1,17 @@ +class CFloat32CustomPython : public astra::CFloat32CustomMemory { +public: + CFloat32CustomPython(PyObject * arrIn) + { + arr = arrIn; + // Set pointer to numpy data pointer + m_fPtr = (float *)PyArray_DATA(arr); + // Increase reference count since ASTRA has a reference + Py_INCREF(arr); + } + virtual ~CFloat32CustomPython() { + // Decrease reference count since ASTRA object is destroyed + Py_DECREF(arr); + } +private: + PyObject* arr; +};
\ No newline at end of file diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index ea624fc..540aad4 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 @@ -41,7 +43,7 @@ cdef extern from "astra/Config.h" namespace "astra": cdef cppclass Config: Config() void initialize(string rootname) - XMLNode *self + XMLNode self cdef extern from "astra/VolumeGeometry2D.h" namespace "astra": cdef cppclass CVolumeGeometry2D: @@ -61,10 +63,14 @@ cdef extern from "astra/VolumeGeometry2D.h" namespace "astra": float32 getWindowMaxY() Config* getConfiguration() +cdef extern from "astra/Float32Data2D.h" namespace "astra": + 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() @@ -132,6 +138,7 @@ cdef extern from "astra/ParallelProjectionGeometry2D.h" namespace "astra": cdef extern from "astra/Float32ProjectionData2D.h" namespace "astra": cdef cppclass CFloat32ProjectionData2D: CFloat32ProjectionData2D(CProjectionGeometry2D*) + CFloat32ProjectionData2D(CProjectionGeometry2D*, CFloat32CustomMemory*) CProjectionGeometry2D * getGeometry() void changeGeometry(CProjectionGeometry2D*) int getDetectorCount() @@ -154,6 +161,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) @@ -184,18 +205,30 @@ cdef extern from "astra/VolumeGeometry3D.h" namespace "astra": 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 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() + cdef extern from "astra/ParallelProjectionGeometry3D.h" namespace "astra": @@ -219,7 +252,13 @@ 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() cdef extern from "astra/Float32Data3D.h" namespace "astra": cdef cppclass CFloat32Data3D: 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 <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +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/PyXMLDocument.pxd b/python/astra/PyXMLDocument.pxd index 69781f1..033b8ef 100644 --- a/python/astra/PyXMLDocument.pxd +++ b/python/astra/PyXMLDocument.pxd @@ -44,22 +44,24 @@ cdef extern from "astra/Globals.h" namespace "astra": cdef extern from "astra/XMLNode.h" namespace "astra": cdef cppclass XMLNode: string getName() - XMLNode *addChildNode(string name) - XMLNode *addChildNode(string, string) + XMLNode addChildNode(string name) + XMLNode addChildNode(string, string) void addAttribute(string, string) void addAttribute(string, float32) void addOption(string, string) bool hasOption(string) string getAttribute(string) - list[XMLNode *] getNodes() + list[XMLNode] getNodes() vector[float32] getContentNumericalArray() + void setContent(double*, int, int, bool) + void setContent(double*, int) string getContent() bool hasAttribute(string) cdef extern from "astra/XMLDocument.h" namespace "astra": cdef cppclass XMLDocument: void saveToFile(string sFilename) - XMLNode *getRootNode() + XMLNode getRootNode() cdef extern from "astra/XMLDocument.h" namespace "astra::XMLDocument": cdef XMLDocument *createDocument(string rootname) diff --git a/python/astra/__init__.py b/python/astra/__init__.py index a61aafc..6c15d30 100644 --- a/python/astra/__init__.py +++ b/python/astra/__init__.py @@ -27,13 +27,15 @@ from . import matlab as m from .creators import astra_dict,create_vol_geom, create_proj_geom, create_backprojection, create_sino, create_reconstruction, create_projector,create_sino3d_gpu, create_backprojection3d_gpu from .functions import data_op, add_noise_to_sino, clear, move_vol_geom from .extrautils import clipCircle -from .ASTRAProjector import ASTRAProjector2D from . import data2d from . import astra from . import data3d from . import algorithm from . import projector +from . import projector3d from . import matrix +from . import log +from .optomo import OpTomo import os try: diff --git a/python/astra/creators.py b/python/astra/creators.py index b5ecd5a..4cd7f2d 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' @@ -268,25 +269,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. @@ -300,13 +299,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) @@ -358,20 +357,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. @@ -382,36 +374,21 @@ def create_sino(data, proj_id=None, proj_geom=None, vol_geom=None, projection. Otherwise, returns a tuple containing the ID of the forward projection and the forward projection itself, in that order. - - 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. """ - 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 @@ -509,8 +486,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'] = {} @@ -573,4 +549,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/data2d.py b/python/astra/data2d.py index 8c4be03..f119f05 100644 --- a/python/astra/data2d.py +++ b/python/astra/data2d.py @@ -24,6 +24,7 @@ # #----------------------------------------------------------------------- from . import data2d_c as d +import numpy as np def clear(): """Clear all 2D data objects.""" @@ -52,6 +53,26 @@ def create(datatype, geometry, data=None): """ return d.create(datatype,geometry,data) +def link(datatype, geometry, data): + """Link a 2D numpy array with the toolbox. + + :param datatype: Data object type, '-vol' or '-sino'. + :type datatype: :class:`string` + :param geometry: Volume or projection geometry. + :type geometry: :class:`dict` + :param data: Numpy array to link + :type data: :class:`numpy.ndarray` + :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") + return d.create(datatype,geometry,data,True) + def store(i, data): """Fill existing 2D object with data. diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index 882051c..42854b3 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -47,8 +47,18 @@ from .PyIncludes cimport * cimport utils from .utils import wrap_from_bytes +from .pythonutils import geom_size + +import operator + +from six.moves import reduce + cdef CData2DManager * man2d = <CData2DManager * >PyData2DManager.getSingletonPtr() +cdef extern from "CFloat32CustomPython.h": + cdef cppclass CFloat32CustomPython: + CFloat32CustomPython(arrIn) + def clear(): man2d.clear() @@ -61,11 +71,16 @@ def delete(ids): man2d.remove(ids) -def create(datatype, geometry, data=None): +def create(datatype, geometry, data=None, link=False): cdef Config *cfg cdef CVolumeGeometry2D * pGeometry cdef CProjectionGeometry2D * ppGeometry 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 datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) pGeometry = new CVolumeGeometry2D() @@ -73,7 +88,11 @@ def create(datatype, geometry, data=None): del cfg del pGeometry raise Exception('Geometry class not initialized.') - pDataObject2D = <CFloat32Data2D * > new CFloat32VolumeData2D(pGeometry) + if link: + pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) + pDataObject2D = <CFloat32Data2D * > new CFloat32VolumeData2D(pGeometry, pCustom) + else: + pDataObject2D = <CFloat32Data2D * > new CFloat32VolumeData2D(pGeometry) del cfg del pGeometry elif datatype == '-sino': @@ -93,7 +112,11 @@ def create(datatype, geometry, data=None): del cfg del ppGeometry raise Exception('Geometry class not initialized.') - pDataObject2D = <CFloat32Data2D * > new CFloat32ProjectionData2D(ppGeometry) + if link: + pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) + pDataObject2D = <CFloat32Data2D * > new CFloat32ProjectionData2D(ppGeometry, pCustom) + else: + pDataObject2D = <CFloat32Data2D * > new CFloat32ProjectionData2D(ppGeometry) del ppGeometry del cfg else: @@ -103,7 +126,7 @@ def create(datatype, geometry, data=None): del pDataObject2D raise Exception("Couldn't initialize data object.") - fillDataObject(pDataObject2D, data) + if not link: fillDataObject(pDataObject2D, data) return man2d.store(pDataObject2D) diff --git a/python/astra/data3d.py b/python/astra/data3d.py index a2e9201..e5ef6b0 100644 --- a/python/astra/data3d.py +++ b/python/astra/data3d.py @@ -24,6 +24,7 @@ # #----------------------------------------------------------------------- from . import data3d_c as d +import numpy as np def create(datatype,geometry,data=None): """Create a 3D object. @@ -39,6 +40,27 @@ def create(datatype,geometry,data=None): """ return d.create(datatype,geometry,data) +def link(datatype, geometry, data): + """Link a 3D numpy array with the toolbox. + + :param datatype: Data object type, '-vol' or '-sino'. + :type datatype: :class:`string` + :param geometry: Volume or projection geometry. + :type geometry: :class:`dict` + :param data: Numpy array to link + :type data: :class:`numpy.ndarray` + :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") + return d.create(datatype,geometry,data,True) + + def get(i): """Get a 3D object. @@ -90,6 +112,17 @@ def get_geometry(i): """ return d.get_geometry(i) +def change_geometry(i, geometry): + """Change the geometry of a 3D object. + + :param i: ID of object. + :type i: :class:`int` + :param geometry: Volume or projection geometry. + :type geometry: :class:`dict` + + """ + return d.change_geometry(i, geometry) + def dimensions(i): """Get dimensions of a 3D object. diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 4b069f7..3b27ab7 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -45,17 +45,33 @@ from .PyXMLDocument cimport XMLDocument cimport utils from .utils import wrap_from_bytes +from .pythonutils import geom_size + +import operator + +from six.moves import reduce + + cdef CData3DManager * man3d = <CData3DManager * >PyData3DManager.getSingletonPtr() cdef extern from *: CFloat32Data3DMemory * dynamic_cast_mem "dynamic_cast<astra::CFloat32Data3DMemory*>" (CFloat32Data3D * ) except NULL -def create(datatype,geometry,data=None): +cdef extern from "CFloat32CustomPython.h": + cdef cppclass CFloat32CustomPython: + CFloat32CustomPython(arrIn) + +def create(datatype,geometry,data=None, link=False): cdef Config *cfg cdef CVolumeGeometry3D * pGeometry cdef CProjectionGeometry3D * ppGeometry cdef CFloat32Data3DMemory * pDataObject3D cdef CConeProjectionGeometry3D* pppGeometry + 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 datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) pGeometry = new CVolumeGeometry3D() @@ -63,7 +79,11 @@ def create(datatype,geometry,data=None): del cfg del pGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = <CFloat32Data3DMemory * > new CFloat32VolumeData3DMemory(pGeometry) + if link: + pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) + pDataObject3D = <CFloat32Data3DMemory * > new CFloat32VolumeData3DMemory(pGeometry, pCustom) + else: + pDataObject3D = <CFloat32Data3DMemory * > new CFloat32VolumeData3DMemory(pGeometry) del cfg del pGeometry elif datatype == '-sino' or datatype == '-proj3d': @@ -84,7 +104,11 @@ def create(datatype,geometry,data=None): del cfg del ppGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(ppGeometry) + if link: + pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) + pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(ppGeometry, pCustom) + else: + pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(ppGeometry) del ppGeometry del cfg elif datatype == "-sinocone": @@ -94,7 +118,11 @@ def create(datatype,geometry,data=None): del cfg del pppGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(pppGeometry) + 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'.") @@ -102,7 +130,7 @@ def create(datatype,geometry,data=None): del pDataObject3D raise Exception("Couldn't initialize data object.") - fillDataObject(pDataObject3D, data) + if not link: fillDataObject(pDataObject3D, data) pDataObject3D.updateStatistics() @@ -122,6 +150,61 @@ def get_geometry(i): raise Exception("Not a known data object") return geom +def change_geometry(i, geom): + cdef CFloat32Data3DMemory * pDataObject = dynamic_cast_mem(getObject(i)) + cdef CFloat32ProjectionData3DMemory * pDataObject2 + cdef CFloat32VolumeData3DMemory * pDataObject3 + if pDataObject.getType() == THREEPROJECTION: + pDataObject2 = <CFloat32ProjectionData3DMemory * >pDataObject + # TODO: Reduce code duplication here + cfg = utils.dictToConfig(six.b('ProjectionGeometry'), geom) + tpe = wrap_from_bytes(cfg.self.getAttribute(six.b('type'))) + if (tpe == "parallel3d"): + ppGeometry = <CProjectionGeometry3D*> new CParallelProjectionGeometry3D(); + elif (tpe == "parallel3d_vec"): + ppGeometry = <CProjectionGeometry3D*> new CParallelVecProjectionGeometry3D(); + elif (tpe == "cone"): + ppGeometry = <CProjectionGeometry3D*> new CConeProjectionGeometry3D(); + elif (tpe == "cone_vec"): + ppGeometry = <CProjectionGeometry3D*> new CConeVecProjectionGeometry3D(); + else: + raise Exception("Invalid geometry type.") + if not ppGeometry.initialize(cfg[0]): + del cfg + del ppGeometry + raise Exception('Geometry class not initialized.') + del cfg + if (ppGeometry.getDetectorColCount() != pDataObject2.getDetectorColCount() or \ + ppGeometry.getProjectionCount() != pDataObject2.getAngleCount() or \ + ppGeometry.getDetectorRowCount() != pDataObject2.getDetectorRowCount()): + del ppGeometry + raise Exception( + "The dimensions of the data do not match those specified in the geometry.") + pDataObject2.changeGeometry(ppGeometry) + del ppGeometry + + elif pDataObject.getType() == THREEVOLUME: + pDataObject3 = <CFloat32VolumeData3DMemory * >pDataObject + cfg = utils.dictToConfig(six.b('VolumeGeometry'), geom) + pGeometry = new CVolumeGeometry3D() + if not pGeometry.initialize(cfg[0]): + del cfg + del pGeometry + raise Exception('Geometry class not initialized.') + del cfg + if (pGeometry.getGridColCount() != pDataObject3.getColCount() or \ + pGeometry.getGridRowCount() != pDataObject3.getRowCount() or \ + pGeometry.getGridSliceCount() != pDataObject3.getSliceCount()): + del pGeometry + raise Exception( + "The dimensions of the data do not match those specified in the geometry.") + pDataObject3.changeGeometry(pGeometry) + del pGeometry + + else: + raise Exception("Not a known data object") + + cdef fillDataObject(CFloat32Data3DMemory * obj, data): if data is None: fillDataObjectScalar(obj, 0) diff --git a/python/astra/functions.py b/python/astra/functions.py index bbbb355..13ca911 100644 --- a/python/astra/functions.py +++ b/python/astra/functions.py @@ -32,12 +32,17 @@ from . import creators as ac import numpy as np -from six.moves import range +try: + from six.moves import range +except ImportError: + # six 1.3.0 + from six.moves import xrange as range from . import data2d from . import data3d from . import projector from . import algorithm +from . import pythonutils @@ -158,29 +163,7 @@ def geom_size(geom, dim=None): :param dim: Optional axis index to return :type dim: :class:`int` """ - - if 'GridSliceCount' in geom: - # 3D Volume geometry? - s = (geom['GridSliceCount'], geom[ - 'GridRowCount'], geom['GridColCount']) - elif 'GridColCount' in geom: - # 2D Volume geometry? - s = (geom['GridRowCount'], geom['GridColCount']) - elif geom['type'] == 'parallel' or geom['type'] == 'fanflat': - s = (len(geom['ProjectionAngles']), geom['DetectorCount']) - elif geom['type'] == 'parallel3d' or geom['type'] == 'cone': - s = (geom['DetectorRowCount'], len( - geom['ProjectionAngles']), geom['DetectorColCount']) - elif geom['type'] == 'parallel_vec' or geom['type'] == 'fanflat_vec': - s = (geom['Vectors'].shape[0], geom['DetectorCount']) - elif geom['type'] == 'parallel3d_vec' or geom['type'] == 'cone_vec': - s = (geom['DetectorRowCount'], geom[ - 'Vectors'].shape[0], geom['DetectorColCount']) - - if dim != None: - s = s[dim] - - return s + return pythonutils.geom_size(geom,dim) def geom_2vec(proj_geom): diff --git a/python/astra/log.py b/python/astra/log.py new file mode 100644 index 0000000..3ec0df5 --- /dev/null +++ b/python/astra/log.py @@ -0,0 +1,153 @@ +#----------------------------------------------------------------------- +#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 <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- + +from . import log_c as l + +import inspect + +def debug(message): + """Log a debug message. + + :param message: Message to log. + :type message: :class:`string` + """ + prev_f = inspect.getframeinfo(inspect.currentframe().f_back) + l.log_debug(prev_f.filename,prev_f.lineno,message) + +def info(message): + """Log an info message. + + :param message: Message to log. + :type message: :class:`string` + """ + prev_f = inspect.getframeinfo(inspect.currentframe().f_back) + l.log_info(prev_f.filename,prev_f.lineno,message) + +def warn(message): + """Log a warning message. + + :param message: Message to log. + :type message: :class:`string` + """ + prev_f = inspect.getframeinfo(inspect.currentframe().f_back) + l.log_warn(prev_f.filename,prev_f.lineno,message) + +def error(message): + """Log an error message. + + :param message: Message to log. + :type message: :class:`string` + """ + prev_f = inspect.getframeinfo(inspect.currentframe().f_back) + l.log_error(prev_f.filename,prev_f.lineno,message) + +def enable(): + """Enable logging to screen and file.""" + l.log_enable() + +def enableScreen(): + """Enable logging to screen.""" + l.log_enableScreen() + +def enableFile(): + """Enable logging to file (note that a file has to be set).""" + l.log_enableFile() + +def disable(): + """Disable all logging.""" + l.log_disable() + +def disableScreen(): + """Disable logging to screen.""" + l.log_disableScreen() + +def disableFile(): + """Disable logging to file.""" + l.log_disableFile() + +def setFormatFile(fmt): + """Set the format string for log messages. Here are the substitutions you may use: + + %f: Source file name generating the log call. + %n: Source line number where the log call was made. + %m: The message text sent to the logger (after printf formatting). + %d: The current date, formatted using the logger's date format. + %t: The current time, formatted using the logger's time format. + %l: The log level (one of "DEBUG", "INFO", "WARN", or "ERROR"). + %%: A literal percent sign. + + The default format string is "%d %t %f(%n): %l: %m\n". + + :param fmt: Format to use, end with "\n". + :type fmt: :class:`string` + """ + l.log_setFormatFile(fmt) + +def setFormatScreen(fmt): + """Set the format string for log messages. Here are the substitutions you may use: + + %f: Source file name generating the log call. + %n: Source line number where the log call was made. + %m: The message text sent to the logger (after printf formatting). + %d: The current date, formatted using the logger's date format. + %t: The current time, formatted using the logger's time format. + %l: The log level (one of "DEBUG", "INFO", "WARN", or "ERROR"). + %%: A literal percent sign. + + The default format string is "%d %t %f(%n): %l: %m\n". + + :param fmt: Format to use, end with "\n". + :type fmt: :class:`string` + """ + l.log_setFormatScreen(fmt) + +STDOUT=1 +STDERR=2 + +DEBUG=0 +INFO=1 +WARN=2 +ERROR=3 + +def setOutputScreen(fd, level): + """Set which screen to output to, and which level to use. + + :param fd: File descriptor of output screen (STDOUT or STDERR). + :type fd: :class:`int` + :param level: Logging level to use (DEBUG, INFO, WARN, or ERROR). + :type level: :class:`int` + """ + l.log_setOutputScreen(fd, level) + +def setOutputFile(filename, level): + """Set which file to output to, and which level to use. + + :param filename: File name of output file. + :type filename: :class:`string` + :param level: Logging level to use (DEBUG, INFO, WARN, or ERROR). + :type level: :class:`int` + """ + l.log_setOutputFile(filename, level)
\ No newline at end of file diff --git a/python/astra/log_c.pyx b/python/astra/log_c.pyx new file mode 100644 index 0000000..f16329f --- /dev/null +++ b/python/astra/log_c.pyx @@ -0,0 +1,103 @@ +#----------------------------------------------------------------------- +#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 <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +import six + +cdef extern from "astra/Logging.h" namespace "astra": + cdef enum log_level: + LOG_DEBUG + LOG_INFO + LOG_WARN + LOG_ERROR + +cdef extern from "astra/Logging.h" namespace "astra::CLogger": + void debug(const char *sfile, int sline, const char *fmt, ...) + void info(const char *sfile, int sline, const char *fmt, ...) + void warn(const char *sfile, int sline, const char *fmt, ...) + void error(const char *sfile, int sline, const char *fmt, ...) + void setOutputScreen(int fd, log_level m_eLevel) + void setOutputFile(const char *filename, log_level m_eLevel) + void enable() + void enableScreen() + void enableFile() + void disable() + void disableScreen() + void disableFile() + void setFormatFile(const char *fmt) + void setFormatScreen(const char *fmt) + +def log_debug(sfile, sline, message): + cstr = list(map(six.b,(sfile,message))) + debug(cstr[0],sline,cstr[1]) + +def log_info(sfile, sline, message): + cstr = list(map(six.b,(sfile,message))) + info(cstr[0],sline,cstr[1]) + +def log_warn(sfile, sline, message): + cstr = list(map(six.b,(sfile,message))) + warn(cstr[0],sline,cstr[1]) + +def log_error(sfile, sline, message): + cstr = list(map(six.b,(sfile,message))) + error(cstr[0],sline,cstr[1]) + +def log_enable(): + enable() + +def log_enableScreen(): + enableScreen() + +def log_enableFile(): + enableFile() + +def log_disable(): + disable() + +def log_disableScreen(): + disableScreen() + +def log_disableFile(): + disableFile() + +def log_setFormatFile(fmt): + cstr = six.b(fmt) + setFormatFile(cstr) + +def log_setFormatScreen(fmt): + cstr = six.b(fmt) + setFormatScreen(cstr) + +enumList = [LOG_DEBUG,LOG_INFO,LOG_WARN,LOG_ERROR] + +def log_setOutputScreen(fd, level): + setOutputScreen(fd, enumList[level]) + +def log_setOutputFile(filename, level): + cstr = six.b(filename) + setOutputFile(cstr, enumList[level]) diff --git a/python/astra/matrix_c.pyx b/python/astra/matrix_c.pyx index b0d8bc4..d099a75 100644 --- a/python/astra/matrix_c.pyx +++ b/python/astra/matrix_c.pyx @@ -27,7 +27,12 @@ # distutils: libraries = astra import six -from six.moves import range +try: + from six.moves import range +except ImportError: + # six 1.3.0 + from six.moves import xrange as range + import numpy as np import scipy.sparse as ss diff --git a/python/astra/optomo.py b/python/astra/optomo.py new file mode 100644 index 0000000..2937d9c --- /dev/null +++ b/python/astra/optomo.py @@ -0,0 +1,203 @@ +#----------------------------------------------------------------------- +#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 <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- + +from . import data2d +from . import data3d +from . import projector +from . import projector3d +from . import creators +from . import algorithm +from . import functions +import numpy as np +from six.moves import reduce +try: + from six.moves import range +except ImportError: + # six 1.3.0 + from six.moves import xrange as range + +import operator +import scipy.sparse.linalg + +class OpTomo(scipy.sparse.linalg.LinearOperator): + """Object that imitates a projection matrix with a given projector. + + This object can do forward projection by using the ``*`` operator:: + + W = astra.OpTomo(proj_id) + fp = W*image + bp = W.T*sinogram + + It can also be used in minimization methods of the :mod:`scipy.sparse.linalg` module:: + + W = astra.OpTomo(proj_id) + output = scipy.sparse.linalg.lsqr(W,sinogram) + + :param proj_id: ID to a projector. + :type proj_id: :class:`int` + """ + + def __init__(self,proj_id): + self.dtype = np.float32 + try: + self.vg = projector.volume_geometry(proj_id) + self.pg = projector.projection_geometry(proj_id) + self.data_mod = data2d + self.appendString = "" + if projector.is_cuda(proj_id): + self.appendString += "_CUDA" + except Exception: + self.vg = projector3d.volume_geometry(proj_id) + self.pg = projector3d.projection_geometry(proj_id) + self.data_mod = data3d + self.appendString = "3D" + if projector3d.is_cuda(proj_id): + self.appendString += "_CUDA" + + self.vshape = functions.geom_size(self.vg) + self.vsize = reduce(operator.mul,self.vshape) + self.sshape = functions.geom_size(self.pg) + self.ssize = reduce(operator.mul,self.sshape) + + self.shape = (self.ssize, self.vsize) + + self.proj_id = proj_id + + self.T = OpTomoTranspose(self) + + def __checkArray(self, arr, shp): + if len(arr.shape)==1: + arr = arr.reshape(shp) + if arr.dtype != np.float32: + arr = arr.astype(np.float32) + if arr.flags['C_CONTIGUOUS']==False: + arr = np.ascontiguousarray(arr) + return arr + + def _matvec(self,v): + """Implements the forward operator. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + v = self.__checkArray(v, self.vshape) + vid = self.data_mod.link('-vol',self.vg,v) + s = np.zeros(self.sshape,dtype=np.float32) + sid = self.data_mod.link('-sino',self.pg,s) + + cfg = creators.astra_dict('FP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['VolumeDataId'] = vid + cfg['ProjectorId'] = self.proj_id + fp_id = algorithm.create(cfg) + algorithm.run(fp_id) + + algorithm.delete(fp_id) + self.data_mod.delete([vid,sid]) + return s.flatten() + + def rmatvec(self,s): + """Implements the transpose operator. + + :param s: The projection data. + :type s: :class:`numpy.ndarray` + """ + s = self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + v = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,v) + + cfg = creators.astra_dict('BP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + bp_id = algorithm.create(cfg) + algorithm.run(bp_id) + + algorithm.delete(bp_id) + self.data_mod.delete([vid,sid]) + return v.flatten() + + def __mul__(self,v): + """Provides easy forward operator by *. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + # Catch the case of a forward projection of a 2D/3D image + if isinstance(v, np.ndarray) and v.shape==self.vshape: + return self._matvec(v) + return scipy.sparse.linalg.LinearOperator.__mul__(self, v) + + def reconstruct(self, method, s, iterations=1, extraOptions = {}): + """Reconstruct an object. + + :param method: Method to use for reconstruction. + :type method: :class:`string` + :param s: The projection data. + :type s: :class:`numpy.ndarray` + :param iterations: Number of iterations to use. + :type iterations: :class:`int` + :param extraOptions: Extra options to use during reconstruction (i.e. for cfg['option']). + :type extraOptions: :class:`dict` + """ + self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + v = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,v) + cfg = creators.astra_dict(method) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + cfg['option'] = extraOptions + alg_id = algorithm.create(cfg) + algorithm.run(alg_id,iterations) + algorithm.delete(alg_id) + self.data_mod.delete([vid,sid]) + return v + +class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): + """This object provides the transpose operation (``.T``) of the OpTomo object. + + Do not use directly, since it can be accessed as member ``.T`` of + an :class:`OpTomo` object. + """ + def __init__(self,parent): + self.parent = parent + self.dtype = np.float32 + self.shape = (parent.shape[1], parent.shape[0]) + + def _matvec(self, s): + return self.parent.rmatvec(s) + + def rmatvec(self, v): + return self.parent.matvec(v) + + def __mul__(self,s): + # Catch the case of a backprojection of 2D/3D data + if isinstance(s, np.ndarray) and s.shape==self.parent.sshape: + return self._matvec(s) + return scipy.sparse.linalg.LinearOperator.__mul__(self, s) 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 <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# 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 = <CProjector3DManager * >PyProjector3DManager.getSingletonPtr() + +include "config.pxi" + +IF HAVE_CUDA: + cdef extern from *: + CCudaProjector3D* dynamic_cast_cuda_projector "dynamic_cast<astra::CCudaProjector3D*>" (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 = <CProjector2DManager * >PyProjector2DManager.getSingletonPtr() cdef CMatrixManager * manM = <CMatrixManager * >PyMatrixManager.getSingletonPtr() +include "config.pxi" + +IF HAVE_CUDA: + cdef extern from *: + CCudaProjector2D* dynamic_cast_cuda_projector "dynamic_cast<astra::CCudaProjector2D*>" (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) diff --git a/python/astra/pythonutils.py b/python/astra/pythonutils.py new file mode 100644 index 0000000..8ea4af5 --- /dev/null +++ b/python/astra/pythonutils.py @@ -0,0 +1,63 @@ +#----------------------------------------------------------------------- +# 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 <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +"""Additional purely Python functions for PyAstraToolbox. + +.. moduleauthor:: Daniel M. Pelt <D.M.Pelt@cwi.nl> + + +""" + +def geom_size(geom, dim=None): + """Returns the size of a volume or sinogram, based on the projection or volume geometry. + + :param geom: Geometry to calculate size from + :type geometry: :class:`dict` + :param dim: Optional axis index to return + :type dim: :class:`int` + """ + + if 'GridSliceCount' in geom: + # 3D Volume geometry? + s = (geom['GridSliceCount'], geom[ + 'GridRowCount'], geom['GridColCount']) + elif 'GridColCount' in geom: + # 2D Volume geometry? + s = (geom['GridRowCount'], geom['GridColCount']) + elif geom['type'] == 'parallel' or geom['type'] == 'fanflat': + s = (len(geom['ProjectionAngles']), geom['DetectorCount']) + elif geom['type'] == 'parallel3d' or geom['type'] == 'cone': + s = (geom['DetectorRowCount'], len( + geom['ProjectionAngles']), geom['DetectorColCount']) + elif geom['type'] == 'fanflat_vec': + s = (geom['Vectors'].shape[0], geom['DetectorCount']) + elif geom['type'] == 'parallel3d_vec' or geom['type'] == 'cone_vec': + s = (geom['DetectorRowCount'], geom[ + 'Vectors'].shape[0], geom['DetectorColCount']) + + if dim != None: + s = s[dim] + + return s diff --git a/python/astra/utils.pyx b/python/astra/utils.pyx index 0439f1b..ddb37aa 100644 --- a/python/astra/utils.pyx +++ b/python/astra/utils.pyx @@ -26,6 +26,7 @@ # distutils: language = c++ # distutils: libraries = astra +cimport numpy as np import numpy as np import six from libcpp.string cimport string @@ -80,11 +81,12 @@ def wrap_from_bytes(value): return s -cdef void readDict(XMLNode * root, _dc): - cdef XMLNode * listbase - cdef XMLNode * itm +cdef void readDict(XMLNode root, _dc): + cdef XMLNode listbase + cdef XMLNode itm cdef int i cdef int j + cdef double* data dc = convert_item(_dc) for item in dc: @@ -93,45 +95,31 @@ cdef void readDict(XMLNode * root, _dc): if val.size == 0: break listbase = root.addChildNode(item) - listbase.addAttribute(< string > six.b('listsize'), < float32 > val.size) - index = 0 + data = <double*>np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) if val.ndim == 2: - for i in range(val.shape[0]): - for j in range(val.shape[1]): - itm = listbase.addChildNode(six.b('ListItem')) - itm.addAttribute(< string > six.b('index'), < float32 > index) - itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) - index += 1 - del itm + listbase.setContent(data, val.shape[1], val.shape[0], False) elif val.ndim == 1: - for i in range(val.shape[0]): - itm = listbase.addChildNode(six.b('ListItem')) - itm.addAttribute(< string > six.b('index'), < float32 > index) - itm.addAttribute(< string > six.b('value'), < float32 > val[i]) - index += 1 - del itm + listbase.setContent(data, val.shape[0]) else: raise Exception("Only 1 or 2 dimensions are allowed") - del listbase elif isinstance(val, dict): if item == six.b('option') or item == six.b('options') or item == six.b('Option') or item == six.b('Options'): readOptions(root, val) else: itm = root.addChildNode(item) readDict(itm, val) - del itm else: if item == six.b('type'): root.addAttribute(< string > six.b('type'), <string> wrap_to_bytes(val)) else: itm = root.addChildNode(item, wrap_to_bytes(val)) - del itm -cdef void readOptions(XMLNode * node, dc): - cdef XMLNode * listbase - cdef XMLNode * itm +cdef void readOptions(XMLNode node, dc): + cdef XMLNode listbase + cdef XMLNode itm cdef int i cdef int j + cdef double* data for item in dc: val = dc[item] if node.hasOption(item): @@ -141,26 +129,13 @@ cdef void readOptions(XMLNode * node, dc): break listbase = node.addChildNode(six.b('Option')) listbase.addAttribute(< string > six.b('key'), < string > item) - listbase.addAttribute(< string > six.b('listsize'), < float32 > val.size) - index = 0 + data = <double*>np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) if val.ndim == 2: - for i in range(val.shape[0]): - for j in range(val.shape[1]): - itm = listbase.addChildNode(six.b('ListItem')) - itm.addAttribute(< string > six.b('index'), < float32 > index) - itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) - index += 1 - del itm + listbase.setContent(data, val.shape[1], val.shape[0], False) elif val.ndim == 1: - for i in range(val.shape[0]): - itm = listbase.addChildNode(six.b('ListItem')) - itm.addAttribute(< string > six.b('index'), < float32 > index) - itm.addAttribute(< string > six.b('value'), < float32 > val[i]) - index += 1 - del itm + listbase.setContent(data, val.shape[0]) else: raise Exception("Only 1 or 2 dimensions are allowed") - del listbase else: node.addOption(item, wrap_to_bytes(val)) @@ -214,10 +189,10 @@ def stringToPythonValue(inputIn): return str(input) -cdef XMLNode2dict(XMLNode * node): - cdef XMLNode * subnode - cdef list[XMLNode * ] nodes - cdef list[XMLNode * ].iterator it +cdef XMLNode2dict(XMLNode node): + cdef XMLNode subnode + cdef list[XMLNode] nodes + cdef list[XMLNode].iterator it dct = {} opts = {} if node.hasAttribute(six.b('type')): @@ -230,7 +205,6 @@ cdef XMLNode2dict(XMLNode * node): opts[castString(subnode.getAttribute('key'))] = stringToPythonValue(subnode.getAttribute('value')) else: dct[castString(subnode.getName())] = stringToPythonValue(subnode.getContent()) - del subnode inc(it) if len(opts)>0: dct['options'] = opts return dct diff --git a/python/builder.py b/python/builder.py index ddca795..cfdb7d1 100644 --- a/python/builder.py +++ b/python/builder.py @@ -41,9 +41,22 @@ try: usecuda=True except KeyError: pass -cfg = open('astra/config.pxi','w') -cfg.write('DEF HAVE_CUDA=' + str(usecuda) + "\n") -cfg.close() + +cfgToWrite = 'DEF HAVE_CUDA=' + str(usecuda) + "\n" +cfgHasToBeUpdated = True +try: + cfg = open('astra/config.pxi','r') + cfgIn = cfg.read() + cfg.close() + if cfgIn==cfgToWrite: + cfgHasToBeUpdated = False +except IOError: + pass + +if cfgHasToBeUpdated: + cfg = open('astra/config.pxi','w') + cfg.write(cfgToWrite) + cfg.close() cmdclass = { } ext_modules = [ ] diff --git a/python/docSRC/index.rst b/python/docSRC/index.rst index 8d17a4a..b7cc6d6 100644 --- a/python/docSRC/index.rst +++ b/python/docSRC/index.rst @@ -18,7 +18,7 @@ Contents: matrix creators functions - ASTRAProjector + operator matlab astra .. astra diff --git a/python/docSRC/ASTRAProjector.rst b/python/docSRC/operator.rst index 1c267e3..f5369fa 100644 --- a/python/docSRC/ASTRAProjector.rst +++ b/python/docSRC/operator.rst @@ -1,7 +1,7 @@ -Helper class: the :mod:`ASTRAProjector` module +OpTomo class: the :mod:`operator` module ============================================== -.. automodule:: astra.ASTRAProjector +.. automodule:: astra.operator :members: :undoc-members: :show-inheritance: |