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 (limited to 'python/astra') 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 From f603045f5bb41de6bc1ffa93badd932b891f5f1d Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 6 Mar 2015 10:58:50 +0100 Subject: Adjust docstring and samples to new python create_sino function --- python/astra/creators.py | 4 ---- 1 file changed, 4 deletions(-) (limited to 'python/astra') diff --git a/python/astra/creators.py b/python/astra/creators.py index 2e2dc71..68bc8a2 100644 --- a/python/astra/creators.py +++ b/python/astra/creators.py @@ -361,10 +361,6 @@ def create_sino(data, proj_id, returnData=True, gpuIndex=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. """ proj_geom = projector.projection_geometry(proj_id) vol_geom = projector.volume_geometry(proj_id) -- cgit v1.2.3 From 535564ccd8b9563fb52be0dff247b99495942f51 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Tue, 10 Mar 2015 12:54:08 +0100 Subject: Add logging support to Python --- python/astra/__init__.py | 1 + python/astra/log.py | 153 +++++++++++++++++++++++++++++++++++++++++++++++ python/astra/log_c.pyx | 96 +++++++++++++++++++++++++++++ 3 files changed, 250 insertions(+) create mode 100644 python/astra/log.py create mode 100644 python/astra/log_c.pyx (limited to 'python/astra') diff --git a/python/astra/__init__.py b/python/astra/__init__.py index a61aafc..1d3176c 100644 --- a/python/astra/__init__.py +++ b/python/astra/__init__.py @@ -34,6 +34,7 @@ from . import data3d from . import algorithm from . import projector from . import matrix +from . import log import os try: 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 . +# +#----------------------------------------------------------------------- + +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..969cc06 --- /dev/null +++ b/python/astra/log_c.pyx @@ -0,0 +1,96 @@ +#----------------------------------------------------------------------- +#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 + +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): + debug(six.b(sfile),sline,six.b(message)) + +def log_info(sfile, sline, message): + info(six.b(sfile),sline,six.b(message)) + +def log_warn(sfile, sline, message): + warn(six.b(sfile),sline,six.b(message)) + +def log_error(sfile, sline, message): + error(six.b(sfile),sline,six.b(message)) + +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): + setFormatFile(six.b(fmt)) + +def log_setFormatScreen(fmt): + setFormatScreen(six.b(fmt)) + +enumList = [LOG_DEBUG,LOG_INFO,LOG_WARN,LOG_ERROR] + +def log_setOutputScreen(fd, level): + setOutputScreen(fd, enumList[level]) + +def log_setOutputFile(filename, level): + setOutputFile(six.b(filename), enumList[level]) \ No newline at end of file -- cgit v1.2.3 From 1b32573046f33050b9300324e6c74e10abb6caaf Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 9 Apr 2015 15:44:01 +0200 Subject: Add 'link' feature to Python (for 2D and 3D data) --- python/astra/CFloat32CustomPython.h | 17 +++++++++++++++++ python/astra/PyIncludes.pxd | 8 ++++++++ python/astra/data2d.py | 21 +++++++++++++++++++++ python/astra/data2d_c.pyx | 21 +++++++++++++++++---- python/astra/data3d.py | 22 ++++++++++++++++++++++ python/astra/data3d_c.pyx | 27 ++++++++++++++++++++++----- 6 files changed, 107 insertions(+), 9 deletions(-) create mode 100644 python/astra/CFloat32CustomPython.h (limited to 'python/astra') 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 7df02c5..13329d1 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -63,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() @@ -130,6 +134,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() @@ -207,6 +212,7 @@ cdef extern from "astra/ProjectionGeometry3D.h" namespace "astra": cdef extern from "astra/Float32VolumeData3DMemory.h" namespace "astra": cdef cppclass CFloat32VolumeData3DMemory: CFloat32VolumeData3DMemory(CVolumeGeometry3D*) + CFloat32VolumeData3DMemory(CVolumeGeometry3D*, CFloat32CustomMemory*) CVolumeGeometry3D* getGeometry() @@ -231,6 +237,8 @@ cdef extern from "astra/Float32ProjectionData3DMemory.h" namespace "astra": cdef cppclass CFloat32ProjectionData3DMemory: CFloat32ProjectionData3DMemory(CProjectionGeometry3D*) CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*) + CFloat32ProjectionData3DMemory(CProjectionGeometry3D*, CFloat32CustomMemory*) + CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*, CFloat32CustomMemory*) CProjectionGeometry3D* getGeometry() cdef extern from "astra/Float32Data3D.h" namespace "astra": 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 b9c105e..ac54898 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -49,6 +49,10 @@ from .utils import wrap_from_bytes cdef CData2DManager * man2d = PyData2DManager.getSingletonPtr() +cdef extern from "CFloat32CustomPython.h": + cdef cppclass CFloat32CustomPython: + CFloat32CustomPython(arrIn) + def clear(): man2d.clear() @@ -61,11 +65,12 @@ 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 datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) pGeometry = new CVolumeGeometry2D() @@ -73,7 +78,11 @@ def create(datatype, geometry, data=None): del cfg del pGeometry raise Exception('Geometry class not initialized.') - pDataObject2D = new CFloat32VolumeData2D(pGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject2D = new CFloat32VolumeData2D(pGeometry, pCustom) + else: + pDataObject2D = new CFloat32VolumeData2D(pGeometry) del cfg del pGeometry elif datatype == '-sino': @@ -91,7 +100,11 @@ def create(datatype, geometry, data=None): del cfg del ppGeometry raise Exception('Geometry class not initialized.') - pDataObject2D = new CFloat32ProjectionData2D(ppGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject2D = new CFloat32ProjectionData2D(ppGeometry, pCustom) + else: + pDataObject2D = new CFloat32ProjectionData2D(ppGeometry) del ppGeometry del cfg else: @@ -101,7 +114,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..4fdf9d7 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. diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 4b069f7..f2c6e26 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -50,12 +50,17 @@ cdef CData3DManager * man3d = PyData3DManager.getSingletonPtr cdef extern from *: CFloat32Data3DMemory * dynamic_cast_mem "dynamic_cast" (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 datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) pGeometry = new CVolumeGeometry3D() @@ -63,7 +68,11 @@ def create(datatype,geometry,data=None): del cfg del pGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = new CFloat32VolumeData3DMemory(pGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject3D = new CFloat32VolumeData3DMemory(pGeometry, pCustom) + else: + pDataObject3D = new CFloat32VolumeData3DMemory(pGeometry) del cfg del pGeometry elif datatype == '-sino' or datatype == '-proj3d': @@ -84,7 +93,11 @@ def create(datatype,geometry,data=None): del cfg del ppGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = new CFloat32ProjectionData3DMemory(ppGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject3D = new CFloat32ProjectionData3DMemory(ppGeometry, pCustom) + else: + pDataObject3D = new CFloat32ProjectionData3DMemory(ppGeometry) del ppGeometry del cfg elif datatype == "-sinocone": @@ -94,7 +107,11 @@ def create(datatype,geometry,data=None): del cfg del pppGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = new CFloat32ProjectionData3DMemory(pppGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject3D = new CFloat32ProjectionData3DMemory(pppGeometry, pCustom) + else: + pDataObject3D = new CFloat32ProjectionData3DMemory(pppGeometry) else: raise Exception("Invalid datatype. Please specify '-vol' or '-proj3d'.") @@ -102,7 +119,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() -- cgit v1.2.3 From e2c87a5e259c847772c733eefb1b291b2a5b1a6e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 10 Apr 2015 14:29:27 +0200 Subject: Add python data3d.change_geometry --- python/astra/PyIncludes.pxd | 15 +++++++++++++ python/astra/data3d.py | 11 +++++++++ python/astra/data3d_c.pyx | 55 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 81 insertions(+) (limited to 'python/astra') diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 7df02c5..a581f88 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -196,18 +196,29 @@ 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*) CVolumeGeometry3D* getGeometry() + void changeGeometry(CVolumeGeometry3D*) + int getRowCount() + int getColCount() + int getSliceCount() + cdef extern from "astra/ParallelProjectionGeometry3D.h" namespace "astra": @@ -232,6 +243,10 @@ cdef extern from "astra/Float32ProjectionData3DMemory.h" namespace "astra": CFloat32ProjectionData3DMemory(CProjectionGeometry3D*) CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*) 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/data3d.py b/python/astra/data3d.py index a2e9201..4679489 100644 --- a/python/astra/data3d.py +++ b/python/astra/data3d.py @@ -90,6 +90,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..48af032 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -122,6 +122,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 = 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 = new CParallelProjectionGeometry3D(); + elif (tpe == "parallel3d_vec"): + ppGeometry = new CParallelVecProjectionGeometry3D(); + elif (tpe == "cone"): + ppGeometry = new CConeProjectionGeometry3D(); + elif (tpe == "cone_vec"): + ppGeometry = 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 = 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) -- cgit v1.2.3 From f69d9f6bc1704560518da3c30c46e495c0228aac Mon Sep 17 00:00:00 2001 From: Daan Pelt Date: Thu, 30 Apr 2015 11:02:50 +0200 Subject: Check data size when using 'link' function in Python --- python/astra/data2d_c.pyx | 16 ++++++++++++ python/astra/data3d_c.pyx | 17 ++++++++++++ python/astra/functions.py | 25 ++---------------- python/astra/pythonutils.py | 63 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 98 insertions(+), 23 deletions(-) create mode 100644 python/astra/pythonutils.py (limited to 'python/astra') diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index ac54898..29548b5 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -47,6 +47,12 @@ 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 = PyData2DManager.getSingletonPtr() cdef extern from "CFloat32CustomPython.h": @@ -71,6 +77,16 @@ def create(datatype, geometry, data=None, link=False): cdef CProjectionGeometry2D * ppGeometry cdef CFloat32Data2D * pDataObject2D cdef CFloat32CustomMemory * pCustom + + if link: + geomSize = geom_size(geometry) + if len(data.shape)==1: + if data.size!=reduce(operator.mul,geomSize): + raise Exception("The dimensions of the data do not match those specified in the geometry.") + else: + if data.shape!=geomSize: + 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() diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 84472c1..30745b4 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -45,6 +45,13 @@ 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 = PyData3DManager.getSingletonPtr() cdef extern from *: @@ -61,6 +68,16 @@ def create(datatype,geometry,data=None, link=False): cdef CFloat32Data3DMemory * pDataObject3D cdef CConeProjectionGeometry3D* pppGeometry cdef CFloat32CustomMemory * pCustom + + if link: + geomSize = geom_size(geometry) + if len(data.shape)==1: + if data.size!=reduce(operator.mul,geomSize): + raise Exception("The dimensions of the data do not match those specified in the geometry.") + else: + if data.shape!=geomSize: + 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() diff --git a/python/astra/functions.py b/python/astra/functions.py index 4025468..b826b86 100644 --- a/python/astra/functions.py +++ b/python/astra/functions.py @@ -38,6 +38,7 @@ from . import data2d from . import data3d from . import projector from . import algorithm +from . import pythonutils @@ -158,29 +159,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'] == '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/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 . +# +#----------------------------------------------------------------------- +"""Additional purely Python functions for PyAstraToolbox. + +.. moduleauthor:: Daniel M. Pelt + + +""" + +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 -- cgit v1.2.3 From 2bc0d98c413fee4108115f26aa337f65337eec55 Mon Sep 17 00:00:00 2001 From: Daan Pelt Date: Thu, 26 Mar 2015 16:40:38 +0100 Subject: Add SPOT-like object for Python (overrides `__mul__` and works with scipy.sparse.linalg) --- python/astra/ASTRAProjector.py | 135 -------------------------- python/astra/__init__.py | 2 +- python/astra/operator.py | 208 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 209 insertions(+), 136 deletions(-) delete mode 100644 python/astra/ASTRAProjector.py create mode 100644 python/astra/operator.py (limited to 'python/astra') diff --git a/python/astra/ASTRAProjector.py b/python/astra/ASTRAProjector.py deleted file mode 100644 index f282618..0000000 --- a/python/astra/ASTRAProjector.py +++ /dev/null @@ -1,135 +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 . -# -#----------------------------------------------------------------------- - -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` - """ - - def __init__(self, proj_geom, vol_geom, proj_type): - 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.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, 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, 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/__init__.py b/python/astra/__init__.py index 063dc16..8c1740c 100644 --- a/python/astra/__init__.py +++ b/python/astra/__init__.py @@ -27,7 +27,6 @@ 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 @@ -36,6 +35,7 @@ from . import projector from . import projector3d from . import matrix from . import log +from .operator import OpTomo import os try: diff --git a/python/astra/operator.py b/python/astra/operator.py new file mode 100644 index 0000000..a3abd5a --- /dev/null +++ b/python/astra/operator.py @@ -0,0 +1,208 @@ +#----------------------------------------------------------------------- +#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 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 range, reduce +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 matmat(self,m): + """Implements the forward operator with a matrix. + + :param m: Volumes to forward project, arranged in columns. + :type m: :class:`numpy.ndarray` + """ + out = np.zeros((self.ssize,m.shape[1]),dtype=np.float32) + for i in range(m.shape[1]): + out[:,i] = self.matvec(m[:,i].flatten()) + return out + + def __mul__(self,v): + """Provides easy forward operator by *. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + return self.matvec(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 matmat(self, m): + out = np.zeros((self.vsize,m.shape[1]),dtype=np.float32) + for i in range(m.shape[1]): + out[:,i] = self.matvec(m[:,i].flatten()) + return out + + def __mul__(self,v): + return self.matvec(v) -- cgit v1.2.3 From 47fe3421585302f2101691a685ab99b0e1ad5cfc Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 1 May 2015 17:48:32 +0200 Subject: Change XMLNode* to XMLNode An XMLNode object is already simply a pointer, so no need to dynamically allocate XMLNodes. --- python/astra/PyIncludes.pxd | 2 +- python/astra/PyXMLDocument.pxd | 8 ++++---- python/astra/utils.pyx | 29 ++++++++++------------------- 3 files changed, 15 insertions(+), 24 deletions(-) (limited to 'python/astra') diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 1d8285b..909f58f 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -43,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: diff --git a/python/astra/PyXMLDocument.pxd b/python/astra/PyXMLDocument.pxd index 69781f1..57c447e 100644 --- a/python/astra/PyXMLDocument.pxd +++ b/python/astra/PyXMLDocument.pxd @@ -44,14 +44,14 @@ 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() string getContent() bool hasAttribute(string) @@ -59,7 +59,7 @@ cdef extern from "astra/XMLNode.h" namespace "astra": 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/utils.pyx b/python/astra/utils.pyx index 0439f1b..8f1e0b7 100644 --- a/python/astra/utils.pyx +++ b/python/astra/utils.pyx @@ -80,9 +80,9 @@ 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 @@ -102,34 +102,29 @@ cdef void readDict(XMLNode * root, _dc): itm.addAttribute(< string > six.b('index'), < float32 > index) itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) index += 1 - del itm 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 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'), 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 for item in dc: @@ -150,17 +145,14 @@ cdef void readOptions(XMLNode * node, dc): itm.addAttribute(< string > six.b('index'), < float32 > index) itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) index += 1 - del itm 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 else: raise Exception("Only 1 or 2 dimensions are allowed") - del listbase else: node.addOption(item, wrap_to_bytes(val)) @@ -214,10 +206,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 +222,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 -- cgit v1.2.3 From fff7470f1d74b0085355130350fa834ea8d37069 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 6 May 2015 13:50:11 +0200 Subject: Make XML array handling consistent setContent and getContent were using different XML formats previously. --- python/astra/PyXMLDocument.pxd | 2 ++ python/astra/utils.pyx | 35 +++++++++-------------------------- 2 files changed, 11 insertions(+), 26 deletions(-) (limited to 'python/astra') diff --git a/python/astra/PyXMLDocument.pxd b/python/astra/PyXMLDocument.pxd index 57c447e..033b8ef 100644 --- a/python/astra/PyXMLDocument.pxd +++ b/python/astra/PyXMLDocument.pxd @@ -53,6 +53,8 @@ cdef extern from "astra/XMLNode.h" namespace "astra": string getAttribute(string) list[XMLNode] getNodes() vector[float32] getContentNumericalArray() + void setContent(double*, int, int, bool) + void setContent(double*, int) string getContent() bool hasAttribute(string) diff --git a/python/astra/utils.pyx b/python/astra/utils.pyx index 8f1e0b7..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 @@ -85,6 +86,7 @@ cdef void readDict(XMLNode root, _dc): cdef XMLNode itm cdef int i cdef int j + cdef double* data dc = convert_item(_dc) for item in dc: @@ -93,21 +95,11 @@ 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 = 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 + 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 + listbase.setContent(data, val.shape[0]) else: raise Exception("Only 1 or 2 dimensions are allowed") elif isinstance(val, dict): @@ -127,6 +119,7 @@ cdef void readOptions(XMLNode node, dc): cdef XMLNode itm cdef int i cdef int j + cdef double* data for item in dc: val = dc[item] if node.hasOption(item): @@ -136,21 +129,11 @@ 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 = 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 + 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 + listbase.setContent(data, val.shape[0]) else: raise Exception("Only 1 or 2 dimensions are allowed") else: -- cgit v1.2.3 From f730efe78367e8fe8e589c2b43fb0886d384f5c8 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 7 May 2015 11:34:37 +0200 Subject: Do not allow 1D input in Python link method --- python/astra/data2d_c.pyx | 10 ++-------- python/astra/data3d_c.pyx | 10 ++-------- 2 files changed, 4 insertions(+), 16 deletions(-) (limited to 'python/astra') diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index 29548b5..4919bf2 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -78,14 +78,8 @@ def create(datatype, geometry, data=None, link=False): cdef CFloat32Data2D * pDataObject2D cdef CFloat32CustomMemory * pCustom - if link: - geomSize = geom_size(geometry) - if len(data.shape)==1: - if data.size!=reduce(operator.mul,geomSize): - raise Exception("The dimensions of the data do not match those specified in the geometry.") - else: - if data.shape!=geomSize: - raise Exception("The dimensions of the data do not match those specified in the geometry.") + 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) diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 30745b4..3b27ab7 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -69,14 +69,8 @@ def create(datatype,geometry,data=None, link=False): cdef CConeProjectionGeometry3D* pppGeometry cdef CFloat32CustomMemory * pCustom - if link: - geomSize = geom_size(geometry) - if len(data.shape)==1: - if data.size!=reduce(operator.mul,geomSize): - raise Exception("The dimensions of the data do not match those specified in the geometry.") - else: - if data.shape!=geomSize: - raise Exception("The dimensions of the data do not match those specified in the geometry.") + 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) -- cgit v1.2.3 From b9b9c82f8634f9c77416de5e857d107005cccbdf Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 7 May 2015 15:40:17 +0200 Subject: Use superclass __mul__ in Python OpTomo __mul__ --- python/astra/operator.py | 33 +++++++++++---------------------- 1 file changed, 11 insertions(+), 22 deletions(-) (limited to 'python/astra') diff --git a/python/astra/operator.py b/python/astra/operator.py index a3abd5a..0c37353 100644 --- a/python/astra/operator.py +++ b/python/astra/operator.py @@ -91,7 +91,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): arr = np.ascontiguousarray(arr) return arr - def matvec(self,v): + def _matvec(self,v): """Implements the forward operator. :param v: Volume to forward project. @@ -135,24 +135,16 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): self.data_mod.delete([vid,sid]) return v.flatten() - def matmat(self,m): - """Implements the forward operator with a matrix. - - :param m: Volumes to forward project, arranged in columns. - :type m: :class:`numpy.ndarray` - """ - out = np.zeros((self.ssize,m.shape[1]),dtype=np.float32) - for i in range(m.shape[1]): - out[:,i] = self.matvec(m[:,i].flatten()) - return out - def __mul__(self,v): """Provides easy forward operator by *. :param v: Volume to forward project. :type v: :class:`numpy.ndarray` """ - return self.matvec(v) + # 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. @@ -192,17 +184,14 @@ class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): self.dtype = np.float32 self.shape = (parent.shape[1], parent.shape[0]) - def matvec(self, s): + def _matvec(self, s): return self.parent.rmatvec(s) def rmatvec(self, v): return self.parent.matvec(v) - def matmat(self, m): - out = np.zeros((self.vsize,m.shape[1]),dtype=np.float32) - for i in range(m.shape[1]): - out[:,i] = self.matvec(m[:,i].flatten()) - return out - - def __mul__(self,v): - return self.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) -- cgit v1.2.3 From 89da933904262d6b7e80e8adf85ca9d1273881b3 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Fri, 8 May 2015 13:48:39 +0200 Subject: Rename optomo.py for Python 2 support --- python/astra/__init__.py | 2 +- python/astra/operator.py | 197 ----------------------------------------------- python/astra/optomo.py | 197 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 198 insertions(+), 198 deletions(-) delete mode 100644 python/astra/operator.py create mode 100644 python/astra/optomo.py (limited to 'python/astra') diff --git a/python/astra/__init__.py b/python/astra/__init__.py index 8c1740c..6c15d30 100644 --- a/python/astra/__init__.py +++ b/python/astra/__init__.py @@ -35,7 +35,7 @@ from . import projector from . import projector3d from . import matrix from . import log -from .operator import OpTomo +from .optomo import OpTomo import os try: diff --git a/python/astra/operator.py b/python/astra/operator.py deleted file mode 100644 index 0c37353..0000000 --- a/python/astra/operator.py +++ /dev/null @@ -1,197 +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 . -# -#----------------------------------------------------------------------- - -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 range, reduce -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/optomo.py b/python/astra/optomo.py new file mode 100644 index 0000000..0c37353 --- /dev/null +++ b/python/astra/optomo.py @@ -0,0 +1,197 @@ +#----------------------------------------------------------------------- +#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 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 range, reduce +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) -- cgit v1.2.3 From 0f4969c3011910518639e3dd9ff6c1e4b2539f4c Mon Sep 17 00:00:00 2001 From: Daan Pelt Date: Mon, 18 May 2015 20:32:52 +0200 Subject: Hold reference to results of six.b calls --- python/astra/log_c.pyx | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) (limited to 'python/astra') diff --git a/python/astra/log_c.pyx b/python/astra/log_c.pyx index 969cc06..f16329f 100644 --- a/python/astra/log_c.pyx +++ b/python/astra/log_c.pyx @@ -52,16 +52,20 @@ cdef extern from "astra/Logging.h" namespace "astra::CLogger": void setFormatScreen(const char *fmt) def log_debug(sfile, sline, message): - debug(six.b(sfile),sline,six.b(message)) + cstr = list(map(six.b,(sfile,message))) + debug(cstr[0],sline,cstr[1]) def log_info(sfile, sline, message): - info(six.b(sfile),sline,six.b(message)) + cstr = list(map(six.b,(sfile,message))) + info(cstr[0],sline,cstr[1]) def log_warn(sfile, sline, message): - warn(six.b(sfile),sline,six.b(message)) + cstr = list(map(six.b,(sfile,message))) + warn(cstr[0],sline,cstr[1]) def log_error(sfile, sline, message): - error(six.b(sfile),sline,six.b(message)) + cstr = list(map(six.b,(sfile,message))) + error(cstr[0],sline,cstr[1]) def log_enable(): enable() @@ -82,10 +86,12 @@ def log_disableFile(): disableFile() def log_setFormatFile(fmt): - setFormatFile(six.b(fmt)) + cstr = six.b(fmt) + setFormatFile(cstr) def log_setFormatScreen(fmt): - setFormatScreen(six.b(fmt)) + cstr = six.b(fmt) + setFormatScreen(cstr) enumList = [LOG_DEBUG,LOG_INFO,LOG_WARN,LOG_ERROR] @@ -93,4 +99,5 @@ def log_setOutputScreen(fd, level): setOutputScreen(fd, enumList[level]) def log_setOutputFile(filename, level): - setOutputFile(six.b(filename), enumList[level]) \ No newline at end of file + cstr = six.b(filename) + setOutputFile(cstr, enumList[level]) -- cgit v1.2.3 From aed157812c31df681a7ba7dd4479cb35fbce05bb Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 18 May 2015 20:54:54 +0200 Subject: Fix range import for old six versions --- python/astra/functions.py | 6 +++++- python/astra/matrix_c.pyx | 7 ++++++- python/astra/optomo.py | 8 +++++++- 3 files changed, 18 insertions(+), 3 deletions(-) (limited to 'python/astra') diff --git a/python/astra/functions.py b/python/astra/functions.py index b826b86..e38b5bc 100644 --- a/python/astra/functions.py +++ b/python/astra/functions.py @@ -32,7 +32,11 @@ 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 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 index 0c37353..2937d9c 100644 --- a/python/astra/optomo.py +++ b/python/astra/optomo.py @@ -32,7 +32,13 @@ from . import creators from . import algorithm from . import functions import numpy as np -from six.moves import range, reduce +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 -- cgit v1.2.3