diff options
Diffstat (limited to 'python/astra')
29 files changed, 3409 insertions, 0 deletions
diff --git a/python/astra/ASTRAProjector.py b/python/astra/ASTRAProjector.py new file mode 100644 index 0000000..96acb10 --- /dev/null +++ b/python/astra/ASTRAProjector.py @@ -0,0 +1,138 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- + +import math +from . import creators as ac +from . import data2d + + +class ASTRAProjector2DTranspose(): + """Implements the ``proj.T`` functionality. + + Do not use directly, since it can be accessed as member ``.T`` of + an :class:`ASTRAProjector2D` object. + + """ + def __init__(self, parentProj): + self.parentProj = parentProj + + def __mul__(self, data): + return self.parentProj.backProject(data) + + +class ASTRAProjector2D(object): + """Helps with various common ASTRA Toolbox 2D operations. + + This class can perform several often used toolbox operations, such as: + + * Forward projecting + * Back projecting + * Reconstructing + + Note that this class has a some computational overhead, because it + copies a lot of data. If you use many repeated operations, directly + using the PyAstraToolbox methods directly is faster. + + You can use this class as an abstracted weight matrix :math:`W`: multiplying an instance + ``proj`` of this class by an image results in a forward projection of the image, and multiplying + ``proj.T`` by a sinogram results in a backprojection of the sinogram:: + + proj = ASTRAProjector2D(...) + fp = proj*image + bp = proj.T*sinogram + + :param proj_geom: The projection geometry. + :type proj_geom: :class:`dict` + :param vol_geom: The volume geometry. + :type vol_geom: :class:`dict` + :param proj_type: Projector type, such as ``'line'``, ``'linear'``, ... + :type proj_type: :class:`string` + :param useCUDA: If ``True``, use CUDA for calculations, when possible. + :type useCUDA: :class:`bool` + """ + + def __init__(self, proj_geom, vol_geom, proj_type, useCUDA=False): + self.vol_geom = vol_geom + self.recSize = vol_geom['GridColCount'] + self.angles = proj_geom['ProjectionAngles'] + self.nDet = proj_geom['DetectorCount'] + nexpow = int(pow(2, math.ceil(math.log(2 * self.nDet, 2)))) + self.filterSize = nexpow / 2 + 1 + self.nProj = self.angles.shape[0] + self.proj_geom = proj_geom + self.proj_id = ac.create_projector(proj_type, proj_geom, vol_geom) + self.useCUDA = useCUDA + self.T = ASTRAProjector2DTranspose(self) + + def backProject(self, data): + """Backproject a sinogram. + + :param data: The sinogram data or ID. + :type data: :class:`numpy.ndarray` or :class:`int` + :returns: :class:`numpy.ndarray` -- The backprojection. + + """ + vol_id, vol = ac.create_backprojection( + data, self.proj_id, useCUDA=self.useCUDA, returnData=True) + data2d.delete(vol_id) + return vol + + def forwardProject(self, data): + """Forward project an image. + + :param data: The image data or ID. + :type data: :class:`numpy.ndarray` or :class:`int` + :returns: :class:`numpy.ndarray` -- The forward projection. + + """ + sin_id, sino = ac.create_sino(data, self.proj_id, useCUDA=self.useCUDA, returnData=True) + data2d.delete(sin_id) + return sino + + def reconstruct(self, data, method, **kwargs): + """Reconstruct an image from a sinogram. + + :param data: The sinogram data or ID. + :type data: :class:`numpy.ndarray` or :class:`int` + :param method: Name of the reconstruction algorithm. + :type method: :class:`string` + :param kwargs: Additional named parameters to pass to :func:`astra.creators.create_reconstruction`. + :returns: :class:`numpy.ndarray` -- The reconstruction. + + Example of a SIRT reconstruction using CUDA:: + + proj = ASTRAProjector2D(...) + rec = proj.reconstruct(sinogram,'SIRT_CUDA',iterations=1000) + + """ + kwargs['returnData'] = True + rec_id, rec = ac.create_reconstruction( + method, self.proj_id, data, **kwargs) + data2d.delete(rec_id) + return rec + + def __mul__(self, data): + return self.forwardProject(data) diff --git a/python/astra/PyAlgorithmFactory.pxd b/python/astra/PyAlgorithmFactory.pxd new file mode 100644 index 0000000..256d7b2 --- /dev/null +++ b/python/astra/PyAlgorithmFactory.pxd @@ -0,0 +1,36 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp.string cimport string +from libcpp cimport bool +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectFactory.h" namespace "astra": + cdef cppclass CAlgorithmFactory: + CAlgorithmFactory() + CAlgorithm *create(string) + +cdef extern from "astra/AstraObjectFactory.h" namespace "astra::CAlgorithmFactory": + cdef CAlgorithmFactory* getSingletonPtr() diff --git a/python/astra/PyAlgorithmManager.pxd b/python/astra/PyAlgorithmManager.pxd new file mode 100644 index 0000000..a99a807 --- /dev/null +++ b/python/astra/PyAlgorithmManager.pxd @@ -0,0 +1,40 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp.string cimport string +from libcpp cimport bool +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectManager.h" namespace "astra": + cdef cppclass CAlgorithmManager: + int store(CAlgorithm *) + CAlgorithm * get(int) + void remove(int) + void clear() + string info() + +cdef extern from "astra/AstraObjectManager.h" namespace "astra::CAlgorithmManager": + cdef CAlgorithmManager* getSingletonPtr() + diff --git a/python/astra/PyData2DManager.pxd b/python/astra/PyData2DManager.pxd new file mode 100644 index 0000000..db8ec84 --- /dev/null +++ b/python/astra/PyData2DManager.pxd @@ -0,0 +1,39 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp.string cimport string + +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectManager.h" namespace "astra": + cdef cppclass CData2DManager: + string info() + void clear() + void remove(int i) + int store(CFloat32Data2D *) + CFloat32Data2D * get(int i) + +cdef extern from "astra/AstraObjectManager.h" namespace "astra::CData2DManager": + cdef CData2DManager* getSingletonPtr() diff --git a/python/astra/PyData3DManager.pxd b/python/astra/PyData3DManager.pxd new file mode 100644 index 0000000..9264a82 --- /dev/null +++ b/python/astra/PyData3DManager.pxd @@ -0,0 +1,39 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp.string cimport string + +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectManager.h" namespace "astra": + cdef cppclass CData3DManager: + string info() + void clear() + void remove(int i) + int store(CFloat32Data3D *) + CFloat32Data3D * get(int i) + +cdef extern from "astra/AstraObjectManager.h" namespace "astra::CData3DManager": + cdef CData3DManager* getSingletonPtr() diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd new file mode 100644 index 0000000..fc5980f --- /dev/null +++ b/python/astra/PyIncludes.pxd @@ -0,0 +1,215 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp cimport bool +from libcpp.string cimport string +from .PyXMLDocument cimport XMLNode + +cdef extern from "astra/Globals.h" namespace "astra": + ctypedef float float32 + ctypedef double float64 + ctypedef unsigned short int uint16 + ctypedef signed short int sint16 + ctypedef unsigned char uchar8 + ctypedef signed char schar8 + ctypedef int int32 + ctypedef short int int16 + +cdef extern from "astra/Config.h" namespace "astra": + cdef cppclass Config: + Config() + XMLNode *self + +cdef extern from "astra/VolumeGeometry2D.h" namespace "astra": + cdef cppclass CVolumeGeometry2D: + bool initialize(Config) + int getGridColCount() + int getGridRowCount() + int getGridTotCount() + float32 getWindowLengthX() + float32 getWindowLengthY() + float32 getWindowArea() + float32 getPixelLengthX() + float32 getPixelLengthY() + float32 getPixelArea() + float32 getWindowMinX() + float32 getWindowMinY() + float32 getWindowMaxX() + float32 getWindowMaxY() + + +cdef extern from "astra/Float32VolumeData2D.h" namespace "astra": + cdef cppclass CFloat32VolumeData2D: + CFloat32VolumeData2D(CVolumeGeometry2D*) + CVolumeGeometry2D * getGeometry() + int getWidth() + int getHeight() + void changeGeometry(CVolumeGeometry2D*) + + + +cdef extern from "astra/ProjectionGeometry2D.h" namespace "astra": + cdef cppclass CProjectionGeometry2D: + CProjectionGeometry2D() + bool initialize(Config) + int getDetectorCount() + int getProjectionAngleCount() + bool isOfType(string) + float32 getProjectionAngle(int) + float32 getDetectorWidth() + +cdef extern from "astra/Float32Data2D.h" namespace "astra::CFloat32Data2D": + cdef enum EDataType: + BASE,PROJECTION,VOLUME + +cdef extern from "astra/Float32Data2D.h" namespace "astra": + cdef cppclass CFloat32Data2D: + bool isInitialized() + int getSize() + float32 *getData() + float32 **getData2D() + int getWidth() + int getHeight() + EDataType getType() + + + + +cdef extern from "astra/SparseMatrixProjectionGeometry2D.h" namespace "astra": + cdef cppclass CSparseMatrixProjectionGeometry2D: + CSparseMatrixProjectionGeometry2D() + +cdef extern from "astra/FanFlatProjectionGeometry2D.h" namespace "astra": + cdef cppclass CFanFlatProjectionGeometry2D: + CFanFlatProjectionGeometry2D() + +cdef extern from "astra/FanFlatVecProjectionGeometry2D.h" namespace "astra": + cdef cppclass CFanFlatVecProjectionGeometry2D: + CFanFlatVecProjectionGeometry2D() + +cdef extern from "astra/ParallelProjectionGeometry2D.h" namespace "astra": + cdef cppclass CParallelProjectionGeometry2D: + CParallelProjectionGeometry2D() + + +cdef extern from "astra/Float32ProjectionData2D.h" namespace "astra": + cdef cppclass CFloat32ProjectionData2D: + CFloat32ProjectionData2D(CProjectionGeometry2D*) + CProjectionGeometry2D * getGeometry() + void changeGeometry(CProjectionGeometry2D*) + int getDetectorCount() + int getAngleCount() + +cdef extern from "astra/Algorithm.h" namespace "astra": + cdef cppclass CAlgorithm: + bool initialize(Config) + void run(int) + bool isInitialized() + +cdef extern from "astra/ReconstructionAlgorithm2D.h" namespace "astra": + cdef cppclass CReconstructionAlgorithm2D: + bool getResidualNorm(float32&) + +cdef extern from "astra/Projector2D.h" namespace "astra": + cdef cppclass CProjector2D: + bool isInitialized() + CProjectionGeometry2D* getProjectionGeometry() + CVolumeGeometry2D* getVolumeGeometry() + CSparseMatrix* getMatrix() + +cdef extern from "astra/SparseMatrix.h" namespace "astra": + cdef cppclass CSparseMatrix: + CSparseMatrix(unsigned int,unsigned int,unsigned long) + unsigned int m_iWidth + unsigned int m_iHeight + unsigned long m_lSize + bool isInitialized() + float32* m_pfValues + unsigned int* m_piColIndices + unsigned long* m_plRowStarts + +cdef extern from "astra/Float32Data3DMemory.h" namespace "astra": + cdef cppclass CFloat32Data3DMemory: + CFloat32Data3DMemory() + bool isInitialized() + int getSize() + int getWidth() + int getHeight() + int getDepth() + void updateStatistics() + float32 *getData() + float32 ***getData3D() + + +cdef extern from "astra/VolumeGeometry3D.h" namespace "astra": + cdef cppclass CVolumeGeometry3D: + CVolumeGeometry3D() + bool initialize(Config) + +cdef extern from "astra/ProjectionGeometry3D.h" namespace "astra": + cdef cppclass CProjectionGeometry3D: + CProjectionGeometry3D() + bool initialize(Config) + + +cdef extern from "astra/Float32VolumeData3DMemory.h" namespace "astra": + cdef cppclass CFloat32VolumeData3DMemory: + CFloat32VolumeData3DMemory(CVolumeGeometry3D*) + + +cdef extern from "astra/ParallelProjectionGeometry3D.h" namespace "astra": + cdef cppclass CParallelProjectionGeometry3D: + CParallelProjectionGeometry3D() + +cdef extern from "astra/ParallelVecProjectionGeometry3D.h" namespace "astra": + cdef cppclass CParallelVecProjectionGeometry3D: + CParallelVecProjectionGeometry3D() + +cdef extern from "astra/ConeProjectionGeometry3D.h" namespace "astra": + cdef cppclass CConeProjectionGeometry3D: + CConeProjectionGeometry3D() + bool initialize(Config) + +cdef extern from "astra/ConeVecProjectionGeometry3D.h" namespace "astra": + cdef cppclass CConeVecProjectionGeometry3D: + CConeVecProjectionGeometry3D() + +cdef extern from "astra/Float32ProjectionData3DMemory.h" namespace "astra": + cdef cppclass CFloat32ProjectionData3DMemory: + CFloat32ProjectionData3DMemory(CProjectionGeometry3D*) + CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*) + +cdef extern from "astra/Float32Data3D.h" namespace "astra": + cdef cppclass CFloat32Data3D: + CFloat32Data3D() + bool isInitialized() + int getSize() + int getWidth() + int getHeight() + int getDepth() + void updateStatistics() + + + diff --git a/python/astra/PyMatrixManager.pxd b/python/astra/PyMatrixManager.pxd new file mode 100644 index 0000000..b2b84c4 --- /dev/null +++ b/python/astra/PyMatrixManager.pxd @@ -0,0 +1,39 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp.string cimport string + +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectManager.h" namespace "astra": + cdef cppclass CMatrixManager: + string info() + void clear() + void remove(int i) + int store(CSparseMatrix *) + CSparseMatrix * get(int i) + +cdef extern from "astra/AstraObjectManager.h" namespace "astra::CMatrixManager": + cdef CMatrixManager* getSingletonPtr() diff --git a/python/astra/PyProjector2DFactory.pxd b/python/astra/PyProjector2DFactory.pxd new file mode 100644 index 0000000..3314544 --- /dev/null +++ b/python/astra/PyProjector2DFactory.pxd @@ -0,0 +1,35 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp.string cimport string +from libcpp cimport bool +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectFactory.h" namespace "astra": + cdef cppclass CProjector2DFactory: + CProjector2D *create(Config) + +cdef extern from "astra/AstraObjectFactory.h" namespace "astra::CProjector2DFactory": + cdef CProjector2DFactory* getSingletonPtr() diff --git a/python/astra/PyProjector2DManager.pxd b/python/astra/PyProjector2DManager.pxd new file mode 100644 index 0000000..92176ba --- /dev/null +++ b/python/astra/PyProjector2DManager.pxd @@ -0,0 +1,39 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp.string cimport string + +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectManager.h" namespace "astra": + cdef cppclass CProjector2DManager: + string info() + void clear() + void remove(int i) + int store(CProjector2D *) + CProjector2D * get(int i) + +cdef extern from "astra/AstraObjectManager.h" namespace "astra::CProjector2DManager": + cdef CProjector2DManager* getSingletonPtr() diff --git a/python/astra/PyXMLDocument.pxd b/python/astra/PyXMLDocument.pxd new file mode 100644 index 0000000..69781f1 --- /dev/null +++ b/python/astra/PyXMLDocument.pxd @@ -0,0 +1,65 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +from libcpp.string cimport string +from libcpp cimport bool +from libcpp.list cimport list +from libcpp.vector cimport vector + +cdef extern from "astra/Globals.h" namespace "astra": + ctypedef float float32 + ctypedef double float64 + ctypedef unsigned short int uint16 + ctypedef signed short int sint16 + ctypedef unsigned char uchar8 + ctypedef signed char schar8 + ctypedef int int32 + ctypedef short int int16 + +cdef extern from "astra/XMLNode.h" namespace "astra": + cdef cppclass XMLNode: + string getName() + 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() + vector[float32] getContentNumericalArray() + string getContent() + bool hasAttribute(string) + +cdef extern from "astra/XMLDocument.h" namespace "astra": + cdef cppclass XMLDocument: + void saveToFile(string sFilename) + XMLNode *getRootNode() + +cdef extern from "astra/XMLDocument.h" namespace "astra::XMLDocument": + cdef XMLDocument *createDocument(string rootname) diff --git a/python/astra/__init__.py b/python/astra/__init__.py new file mode 100644 index 0000000..a61aafc --- /dev/null +++ b/python/astra/__init__.py @@ -0,0 +1,42 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from . import matlab as m +from .creators import astra_dict,create_vol_geom, create_proj_geom, create_backprojection, create_sino, create_reconstruction, create_projector,create_sino3d_gpu, create_backprojection3d_gpu +from .functions import data_op, add_noise_to_sino, clear, move_vol_geom +from .extrautils import clipCircle +from .ASTRAProjector import ASTRAProjector2D +from . import data2d +from . import astra +from . import data3d +from . import algorithm +from . import projector +from . import matrix + +import os +try: + astra.set_gpu_index(int(os.environ['ASTRA_GPU_INDEX'])) +except KeyError: + pass diff --git a/python/astra/algorithm.py b/python/astra/algorithm.py new file mode 100644 index 0000000..46cfccc --- /dev/null +++ b/python/astra/algorithm.py @@ -0,0 +1,76 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- + +from . import algorithm_c as a + +def create(config): + """Create algorithm object. + + :param config: Algorithm options. + :type config: :class:`dict` + :returns: :class:`int` -- the ID of the constructed object. + + """ + return a.create(config) + +def run(i, iterations=1): + """Run an algorithm. + + :param i: ID of object. + :type i: :class:`int` + :param iterations: Number of iterations to run. + :type iterations: :class:`int` + + """ + return a.run(i,iterations) + +def get_res_norm(i): + """Get residual norm of algorithm. + + :param i: ID of object. + :type i: :class:`int` + :returns: :class:`float` -- The residual norm. + + """ + + return a.get_res_norm(i) + +def delete(ids): + """Delete a matrix object. + + :param ids: ID or list of ID's to delete. + :type ids: :class:`int` or :class:`list` + + """ + return a.delete(ids) + +def clear(): + """Clear all matrix objects.""" + return a.clear() + +def info(): + """Print info on matrix objects in memory.""" + return a.info() diff --git a/python/astra/algorithm_c.pyx b/python/astra/algorithm_c.pyx new file mode 100644 index 0000000..0c2999c --- /dev/null +++ b/python/astra/algorithm_c.pyx @@ -0,0 +1,106 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- + +# distutils: language = c++ +# distutils: libraries = astra + +import six +from .PyIncludes cimport * + +cimport PyAlgorithmManager +from .PyAlgorithmManager cimport CAlgorithmManager + +cimport PyAlgorithmFactory +from .PyAlgorithmFactory cimport CAlgorithmFactory + +cimport PyXMLDocument +from .PyXMLDocument cimport XMLDocument + +cimport utils +from .utils import wrap_from_bytes + +cdef CAlgorithmManager * manAlg = <CAlgorithmManager * >PyAlgorithmManager.getSingletonPtr() + +cdef extern from *: + CReconstructionAlgorithm2D * dynamic_cast_recAlg "dynamic_cast<CReconstructionAlgorithm2D*>" (CAlgorithm * ) except NULL + + +def create(config): + cdef XMLDocument * xml = utils.dict2XML(six.b('Algorithm'), config) + cdef Config cfg + cdef CAlgorithm * alg + cfg.self = xml.getRootNode() + alg = PyAlgorithmFactory.getSingletonPtr().create(cfg.self.getAttribute(six.b('type'))) + if alg == NULL: + del xml + raise Exception("Unknown algorithm.") + if not alg.initialize(cfg): + del xml + del alg + raise Exception("Algorithm not initialized.") + del xml + return manAlg.store(alg) + +cdef CAlgorithm * getAlg(i) except NULL: + cdef CAlgorithm * alg = manAlg.get(i) + if alg == NULL: + raise Exception("Unknown algorithm.") + if not alg.isInitialized(): + raise Exception("Algorithm not initialized.") + return alg + + +def run(i, iterations=0): + cdef CAlgorithm * alg = getAlg(i) + alg.run(iterations) + + +def get_res_norm(i): + cdef CReconstructionAlgorithm2D * pAlg2D + cdef CAlgorithm * alg = getAlg(i) + cdef float32 res = 0.0 + pAlg2D = dynamic_cast_recAlg(alg) + if pAlg2D == NULL: + raise Exception("Operation not supported.") + if not pAlg2D.getResidualNorm(res): + raise Exception("Operation not supported.") + return res + + +def delete(ids): + try: + for i in ids: + manAlg.remove(i) + except TypeError: + manAlg.remove(ids) + + +def clear(): + manAlg.clear() + + +def info(): + six.print_(wrap_from_bytes(manAlg.info())) diff --git a/python/astra/astra.py b/python/astra/astra.py new file mode 100644 index 0000000..26b1ff0 --- /dev/null +++ b/python/astra/astra.py @@ -0,0 +1,58 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- + +from . import astra_c as a + +def credits(): + """Print credits of the ASTRA Toolbox.""" + return a.credits() + + +def use_cuda(): + """Test if CUDA is enabled. + + :returns: :class:`bool` -- ``True`` if CUDA is enabled. + """ + return a.use_cuda() + + +def version(printToScreen=False): + """Check version of the ASTRA Toolbox. + + :param printToScreen: If ``True``, print version string. If ``False``, return version integer. + :type printToScreen: :class:`bool` + :returns: :class:`string` or :class:`int` -- The version string or integer. + + """ + return a.version(printToScreen) + +def set_gpu_index(idx): + """Set default GPU index to use. + + :param idx: GPU index + :type idx: :class:`int` + """ + a.set_gpu_index(idx) diff --git a/python/astra/astra_c.pyx b/python/astra/astra_c.pyx new file mode 100644 index 0000000..342a214 --- /dev/null +++ b/python/astra/astra_c.pyx @@ -0,0 +1,79 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +include "config.pxi" +import six +from .utils import wrap_from_bytes + +from libcpp.string cimport string +from libcpp cimport bool +cdef extern from "astra/Globals.h" namespace "astra": + int getVersion() + string getVersionString() + bool cudaEnabled() + +IF HAVE_CUDA==True: + cdef extern from "../cuda/2d/darthelper.h" namespace "astraCUDA": + bool setGPUIndex(int) +ELSE: + def setGPUIndex(): + pass + +def credits(): + six.print_(""" +All Scale Tomographic Reconstruction Antwerp Toolbox (ASTRA-Toolbox) +was developed at the University of Antwerp by + * Prof. dr. Joost Batenburg + * Andrei Dabravolski + * Gert Merckx + * Willem Jan Palenstijn + * Tom Roelandts + * Prof. dr. Jan Sijbers + * dr. Wim van Aarle + * Sander van der Maar + * dr. Gert Van Gompel + +Python interface written by + * Daniel M. Pelt (CWI, Amsterdam)""") + + +def use_cuda(): + return cudaEnabled() + + +def version(printToScreen=False): + if printToScreen: + six.print_(wrap_from_bytes(getVersionString())) + else: + return getVersion() + +def set_gpu_index(idx): + if use_cuda()==True: + ret = setGPUIndex(idx) + if not ret: + six.print_("Failed to set GPU " + str(idx)) diff --git a/python/astra/creators.py b/python/astra/creators.py new file mode 100644 index 0000000..9aba464 --- /dev/null +++ b/python/astra/creators.py @@ -0,0 +1,563 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- + +import six +import numpy as np +import math +from . import data2d +from . import data3d +from . import projector +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' + six.print_('SIRT_CUDA2 has been deprecated. Use SIRT_CUDA instead.') + elif intype == 'FP_CUDA2': + intype = 'FP_CUDA' + six.print_('FP_CUDA2 has been deprecated. Use FP_CUDA instead.') + return {'type': intype} + +def create_vol_geom(*varargin): + """Create a volume geometry structure. + +This method can be called in a number of ways: + +``create_vol_geom(N)``: + :returns: A 2D volume geometry of size :math:`N \\times N`. + +``create_vol_geom((M, N))``: + :returns: A 2D volume geometry of size :math:`M \\times N`. + +``create_vol_geom(M, N)``: + :returns: A 2D volume geometry of size :math:`M \\times N`. + +``create_vol_geom(M, N, minx, maxx, miny, maxy)``: + :returns: A 2D volume geometry of size :math:`M \\times N`, windowed as :math:`minx \\leq x \\leq maxx` and :math:`miny \\leq y \\leq maxy`. + +``create_vol_geom((M, N, Z))``: + :returns: A 3D volume geometry of size :math:`M \\times N \\times Z`. + +``create_vol_geom(M, N, Z)``: + :returns: A 3D volume geometry of size :math:`M \\times N \\times Z`. + +""" + vol_geom = {'option': {}} + # astra_create_vol_geom(row_count) + if len(varargin) == 1 and isinstance(varargin[0], int) == 1: + vol_geom['GridRowCount'] = varargin[0] + vol_geom['GridColCount'] = varargin[0] + vol_geom['option']['WindowMinX'] = -varargin[0] / 2. + vol_geom['option']['WindowMaxX'] = varargin[0] / 2. + vol_geom['option']['WindowMinY'] = -varargin[0] / 2. + vol_geom['option']['WindowMaxY'] = varargin[0] / 2. + # astra_create_vol_geom([row_count col_count]) + elif len(varargin) == 1 and len(varargin[0]) == 2: + vol_geom['GridRowCount'] = varargin[0][0] + vol_geom['GridColCount'] = varargin[0][1] + vol_geom['option']['WindowMinX'] = -varargin[0][1] / 2. + vol_geom['option']['WindowMaxX'] = varargin[0][1] / 2. + vol_geom['option']['WindowMinY'] = -varargin[0][0] / 2. + vol_geom['option']['WindowMaxY'] = varargin[0][0] / 2. + # astra_create_vol_geom([row_count col_count slice_count]) + elif len(varargin) == 1 and len(varargin[0]) == 3: + vol_geom['GridRowCount'] = varargin[0][0] + vol_geom['GridColCount'] = varargin[0][1] + vol_geom['GridSliceCount'] = varargin[0][2] + vol_geom['option']['WindowMinX'] = -varargin[0][1] / 2. + vol_geom['option']['WindowMaxX'] = varargin[0][1] / 2. + vol_geom['option']['WindowMinY'] = -varargin[0][0] / 2. + vol_geom['option']['WindowMaxY'] = varargin[0][0] / 2. + vol_geom['option']['WindowMinZ'] = -varargin[0][2] / 2. + vol_geom['option']['WindowMaxZ'] = varargin[0][2] / 2. + # astra_create_vol_geom(row_count, col_count) + elif len(varargin) == 2: + vol_geom['GridRowCount'] = varargin[0] + vol_geom['GridColCount'] = varargin[1] + vol_geom['option']['WindowMinX'] = -varargin[1] / 2. + vol_geom['option']['WindowMaxX'] = varargin[1] / 2. + vol_geom['option']['WindowMinY'] = -varargin[0] / 2. + vol_geom['option']['WindowMaxY'] = varargin[0] / 2. + # astra_create_vol_geom(row_count, col_count, min_x, max_x, min_y, max_y) + elif len(varargin) == 6: + vol_geom['GridRowCount'] = varargin[0] + vol_geom['GridColCount'] = varargin[1] + vol_geom['option']['WindowMinX'] = varargin[2] + vol_geom['option']['WindowMaxX'] = varargin[3] + vol_geom['option']['WindowMinY'] = varargin[4] + vol_geom['option']['WindowMaxY'] = varargin[5] + # astra_create_vol_geom(row_count, col_count, slice_count) + elif len(varargin) == 3: + vol_geom['GridRowCount'] = varargin[0] + vol_geom['GridColCount'] = varargin[1] + vol_geom['GridSliceCount'] = varargin[2] + return vol_geom + + +def create_proj_geom(intype, *args): + """Create a projection geometry. + +This method can be called in a number of ways: + +``create_proj_geom('parallel', detector_spacing, det_count, angles)``: + +:param detector_spacing: Distance between two adjacent detector pixels. +:type detector_spacing: :class:`float` +:param det_count: Number of detector pixels. +:type det_count: :class:`int` +:param angles: Array of angles in radians. +:type angles: :class:`numpy.ndarray` +:returns: A parallel projection geometry. + + +``create_proj_geom('fanflat', det_width, det_count, angles, source_origin, source_det)``: + +:param det_width: Size of a detector pixel. +:type det_width: :class:`float` +:param det_count: Number of detector pixels. +:type det_count: :class:`int` +:param angles: Array of angles in radians. +:type angles: :class:`numpy.ndarray` +:param source_origin: Position of the source. +:param source_det: Position of the detector +:returns: A fan-beam projection geometry. + +``create_proj_geom('fanflat_vec', det_count, V)``: + +:param det_count: Number of detector pixels. +:type det_count: :class:`int` +:param V: Vector array. +:type V: :class:`numpy.ndarray` +:returns: A fan-beam projection geometry. + +``create_proj_geom('parallel3d', detector_spacing_x, detector_spacing_y, det_row_count, det_col_count, angles)``: + +:param detector_spacing_*: Distance between two adjacent detector pixels. +:type detector_spacing_*: :class:`float` +:param det_row_count: Number of detector pixel rows. +:type det_row_count: :class:`int` +:param det_col_count: Number of detector pixel columns. +:type det_col_count: :class:`int` +:param angles: Array of angles in radians. +:type angles: :class:`numpy.ndarray` +:returns: A parallel projection geometry. + +``create_proj_geom('cone', detector_spacing_x, detector_spacing_y, det_row_count, det_col_count, angles, source_origin, source_det)``: + +:param detector_spacing_*: Distance between two adjacent detector pixels. +:type detector_spacing_*: :class:`float` +:param det_row_count: Number of detector pixel rows. +:type det_row_count: :class:`int` +:param det_col_count: Number of detector pixel columns. +:type det_col_count: :class:`int` +:param angles: Array of angles in radians. +:type angles: :class:`numpy.ndarray` +:param source_origin: Distance between point source and origin. +:type source_origin: :class:`float` +:param source_det: Distance between the detector and origin. +:type source_det: :class:`float` +:returns: A cone-beam projection geometry. + +``create_proj_geom('cone_vec', det_row_count, det_col_count, V)``: + +:param det_row_count: Number of detector pixel rows. +:type det_row_count: :class:`int` +:param det_col_count: Number of detector pixel columns. +:type det_col_count: :class:`int` +:param V: Vector array. +:type V: :class:`numpy.ndarray` +:returns: A cone-beam projection geometry. + +``create_proj_geom('parallel3d_vec', det_row_count, det_col_count, V)``: + +:param det_row_count: Number of detector pixel rows. +:type det_row_count: :class:`int` +:param det_col_count: Number of detector pixel columns. +:type det_col_count: :class:`int` +:param V: Vector array. +:type V: :class:`numpy.ndarray` +:returns: A parallel projection geometry. + +``create_proj_geom('sparse_matrix', det_width, det_count, angles, matrix_id)``: + +:param det_width: Size of a detector pixel. +:type det_width: :class:`float` +:param det_count: Number of detector pixels. +:type det_count: :class:`int` +:param angles: Array of angles in radians. +:type angles: :class:`numpy.ndarray` +:param matrix_id: ID of the sparse matrix. +:type matrix_id: :class:`int` +:returns: A projection geometry based on a sparse matrix. + +""" + if intype == 'parallel': + if len(args) < 3: + raise Exception( + 'not enough variables: astra_create_proj_geom(parallel, detector_spacing, det_count, angles)') + return {'type': 'parallel', 'DetectorWidth': args[0], 'DetectorCount': args[1], 'ProjectionAngles': args[2]} + elif intype == 'fanflat': + if len(args) < 5: + raise Exception('not enough variables: astra_create_proj_geom(fanflat, det_width, det_count, angles, source_origin, source_det)') + return {'type': 'fanflat', 'DetectorWidth': args[0], 'DetectorCount': args[1], 'ProjectionAngles': args[2], 'DistanceOriginSource': args[3], 'DistanceOriginDetector': args[4]} + elif intype == 'fanflat_vec': + if len(args) < 2: + raise Exception('not enough variables: astra_create_proj_geom(fanflat_vec, det_count, V)') + if not args[1].shape[1] == 6: + raise Exception('V should be a Nx6 matrix, with N the number of projections') + return {'type':'fanflat_vec', 'DetectorCount':args[0], 'Vectors':args[1]} + elif intype == 'parallel3d': + if len(args) < 5: + raise Exception('not enough variables: astra_create_proj_geom(parallel3d, detector_spacing_x, detector_spacing_y, det_row_count, det_col_count, angles)') + return {'type':'parallel3d', 'DetectorSpacingX':args[0], 'DetectorSpacingY':args[1], 'DetectorRowCount':args[2], 'DetectorColCount':args[3],'ProjectionAngles':args[4]} + elif intype == 'cone': + if len(args) < 7: + raise Exception('not enough variables: astra_create_proj_geom(cone, detector_spacing_x, detector_spacing_y, det_row_count, det_col_count, angles, source_origin, source_det)') + return {'type': 'cone','DetectorSpacingX':args[0], 'DetectorSpacingY':args[1], 'DetectorRowCount':args[2],'DetectorColCount':args[3],'ProjectionAngles':args[4],'DistanceOriginSource': args[5],'DistanceOriginDetector':args[6]} + elif intype == 'cone_vec': + if len(args) < 3: + raise Exception('not enough variables: astra_create_proj_geom(cone_vec, det_row_count, det_col_count, V)') + if not args[2].shape[1] == 12: + raise Exception('V should be a Nx12 matrix, with N the number of projections') + return {'type': 'cone_vec','DetectorRowCount':args[0],'DetectorColCount':args[1],'Vectors':args[2]} + elif intype == 'parallel3d_vec': + if len(args) < 3: + 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]} + 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) + + +def create_backprojection(data, proj_id, useCUDA=False, 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. + +""" + proj_geom = projector.projection_geometry(proj_id) + vol_geom = projector.volume_geometry(proj_id) + if isinstance(data, np.ndarray): + sino_id = data2d.create('-sino', proj_geom, data) + else: + sino_id = data + vol_id = data2d.create('-vol', vol_geom, 0) + + algString = 'BP' + if useCUDA: + algString = algString + '_CUDA' + + cfg = astra_dict(algString) + if not useCUDA: + cfg['ProjectorId'] = proj_id + cfg['ProjectionDataId'] = sino_id + cfg['ReconstructionDataId'] = vol_id + alg_id = algorithm.create(cfg) + algorithm.run(alg_id) + algorithm.delete(alg_id) + + if isinstance(data, np.ndarray): + data2d.delete(sino_id) + + if returnData: + return vol_id, data2d.get(vol_id) + else: + return vol_id + +def create_backprojection3d_gpu(data, proj_geom, vol_geom, returnData=True): + """Create a backprojection of a sinogram (3D) using CUDA. + +:param data: Sinogram data or ID. +:type data: :class:`numpy.ndarray` or :class:`int` +:param proj_geom: Projection geometry. +:type proj_geom: :class:`dict` +:param vol_geom: Volume geometry. +:type vol_geom: :class:`dict` +: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. + +""" + if isinstance(data, np.ndarray): + sino_id = data3d.create('-sino', proj_geom, data) + else: + sino_id = data + + vol_id = data3d.create('-vol', vol_geom, 0) + + cfg = astra_dict('BP3D_CUDA') + cfg['ProjectionDataId'] = sino_id + cfg['ReconstructionDataId'] = vol_id + alg_id = algorithm.create(cfg) + algorithm.run(alg_id) + algorithm.delete(alg_id) + + if isinstance(data, np.ndarray): + data3d.delete(sino_id) + + if returnData: + return vol_id, data3d.get(vol_id) + else: + return vol_id + + +def create_sino(data, proj_id=None, proj_geom=None, vol_geom=None, + useCUDA=False, 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. + :type gpuIndex: :class:`int` + :returns: :class:`int` or (:class:`int`, :class:`numpy.ndarray`) + + If ``returnData=False``, returns the ID of the forward + projection. Otherwise, returns a tuple containing the ID of the + forward projection and the forward projection itself, in that + order. + + The geometry of setup is defined by ``proj_id`` or with + ``proj_geom`` and ``vol_geom``. If ``proj_id`` is given, then + ``proj_geom`` and ``vol_geom`` must be None and vice versa. +""" + if proj_id is not None: + proj_geom = projector.projection_geometry(proj_id) + vol_geom = projector.volume_geometry(proj_id) + elif proj_geom is not None and vol_geom is not None: + if not useCUDA: + # We need more parameters to create projector. + raise ValueError( + """A ``proj_id`` is needed when CUDA is not used.""") + else: + raise Exception("""The geometry setup is not defined. + The geometry of setup is defined by ``proj_id`` or with + ``proj_geom`` and ``vol_geom``. If ``proj_id`` is given, then + ``proj_geom`` and ``vol_geom`` must be None and vice versa.""") + + 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' + cfg = astra_dict(algString) + if not useCUDA: + cfg['ProjectorId'] = proj_id + if gpuIndex is not None: + cfg['option'] = {'GPUindex': gpuIndex} + cfg['ProjectionDataId'] = sino_id + cfg['VolumeDataId'] = volume_id + alg_id = algorithm.create(cfg) + algorithm.run(alg_id) + algorithm.delete(alg_id) + + if isinstance(data, np.ndarray): + data2d.delete(volume_id) + if returnData: + return sino_id, data2d.get(sino_id) + else: + return sino_id + + + +def create_sino3d_gpu(data, proj_geom, vol_geom, returnData=True, gpuIndex=None): + """Create a forward projection of an image (3D). + +:param data: Image data or ID. +:type data: :class:`numpy.ndarray` or :class:`int` +:param proj_geom: Projection geometry. +:type proj_geom: :class:`dict` +:param vol_geom: Volume geometry. +:type vol_geom: :class:`dict` +:param returnData: If False, only return the ID of the forward projection. +:type returnData: :class:`bool` +:param gpuIndex: Optional GPU index. +:type gpuIndex: :class:`int` +:returns: :class:`int` or (:class:`int`, :class:`numpy.ndarray`) -- If ``returnData=False``, returns the ID of the forward projection. Otherwise, returns a tuple containing the ID of the forward projection and the forward projection itself, in that order. + +""" + + if isinstance(data, np.ndarray): + volume_id = data3d.create('-vol', vol_geom, data) + else: + volume_id = data + sino_id = data3d.create('-sino', proj_geom, 0) + algString = 'FP3D_CUDA' + cfg = astra_dict(algString) + if not gpuIndex==None: + cfg['option']={'GPUindex':gpuIndex} + cfg['ProjectionDataId'] = sino_id + cfg['VolumeDataId'] = volume_id + alg_id = algorithm.create(cfg) + algorithm.run(alg_id) + algorithm.delete(alg_id) + + if isinstance(data, np.ndarray): + data3d.delete(volume_id) + if returnData: + return sino_id, data3d.get(sino_id) + else: + return sino_id + + +def create_reconstruction(rec_type, proj_id, sinogram, iterations=1, use_mask='no', mask=np.array([]), use_minc='no', minc=0, use_maxc='no', maxc=255, returnData=True, filterType=None, filterData=None): + """Create a reconstruction of a sinogram (2D). + +:param rec_type: Name of the reconstruction algorithm. +:type rec_type: :class:`string` +:param proj_id: ID of the projector to use. +:type proj_id: :class:`int` +:param sinogram: Sinogram data or ID. +:type sinogram: :class:`numpy.ndarray` or :class:`int` +:param iterations: Number of iterations to run. +:type iterations: :class:`int` +:param use_mask: Whether to use a mask. +:type use_mask: ``'yes'`` or ``'no'`` +:param mask: Mask data or ID +:type mask: :class:`numpy.ndarray` or :class:`int` +:param use_minc: Whether to force a minimum value on the reconstruction pixels. +:type use_minc: ``'yes'`` or ``'no'`` +:param minc: Minimum value to use. +:type minc: :class:`float` +:param use_maxc: Whether to force a maximum value on the reconstruction pixels. +:type use_maxc: ``'yes'`` or ``'no'`` +:param maxc: Maximum value to use. +:type maxc: :class:`float` +:param returnData: If False, only return the ID of the reconstruction. +:type returnData: :class:`bool` +:param filterType: Which type of filter to use for filter-based methods. +:type filterType: :class:`string` +:param filterData: Optional filter data for filter-based methods. +:type filterData: :class:`numpy.ndarray` +:returns: :class:`int` or (:class:`int`, :class:`numpy.ndarray`) -- If ``returnData=False``, returns the ID of the reconstruction. Otherwise, returns a tuple containing the ID of the reconstruction and reconstruction itself, in that order. + +""" + proj_geom = projector.projection_geometry(proj_id) + if isinstance(sinogram, np.ndarray): + sino_id = data2d.create('-sino', proj_geom, sinogram) + else: + sino_id = sinogram + 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['ProjectionDataId'] = sino_id + cfg['ReconstructionDataId'] = recon_id + cfg['options'] = {} + if use_mask == 'yes': + if isinstance(mask, np.ndarray): + mask_id = data2d.create('-vol', vol_geom, mask) + else: + mask_id = mask + cfg['options']['ReconstructionMaskId'] = mask_id + if not filterType == None: + cfg['FilterType'] = filterType + if not filterData == None: + if isinstance(filterData, np.ndarray): + nexpow = int( + pow(2, math.ceil(math.log(2 * proj_geom['DetectorCount'], 2)))) + filtSize = nexpow / 2 + 1 + filt_proj_geom = create_proj_geom( + 'parallel', 1.0, filtSize, proj_geom['ProjectionAngles']) + filt_id = data2d.create('-sino', filt_proj_geom, filterData) + else: + filt_id = filterData + cfg['FilterSinogramId'] = filt_id + cfg['options']['UseMinConstraint'] = use_minc + cfg['options']['MinConstraintValue'] = minc + cfg['options']['UseMaxConstraint'] = use_maxc + cfg['options']['MaxConstraintValue'] = maxc + cfg['options']['ProjectionOrder'] = 'random' + alg_id = algorithm.create(cfg) + algorithm.run(alg_id, iterations) + + algorithm.delete(alg_id) + + if isinstance(sinogram, np.ndarray): + data2d.delete(sino_id) + if use_mask == 'yes' and isinstance(mask, np.ndarray): + data2d.delete(mask_id) + if not filterData == None: + if isinstance(filterData, np.ndarray): + data2d.delete(filt_id) + if returnData: + return recon_id, data2d.get(recon_id) + else: + return recon_id + + +def create_projector(proj_type, proj_geom, vol_geom): + """Create a 2D projector. + +:param proj_type: Projector type, such as ``'line'``, ``'linear'``, ... +:type proj_type: :class:`string` +:param proj_geom: Projection geometry. +:type proj_geom: :class:`dict` +:param vol_geom: Volume geometry. +:type vol_geom: :class:`dict` +:returns: :class:`int` -- The ID of the projector. + +""" + if proj_type == 'blob': + raise Exception('Blob type not yet implemented') + cfg = astra_dict(proj_type) + cfg['ProjectionGeometry'] = proj_geom + cfg['VolumeGeometry'] = vol_geom + return projector.create(cfg) diff --git a/python/astra/data2d.py b/python/astra/data2d.py new file mode 100644 index 0000000..8c4be03 --- /dev/null +++ b/python/astra/data2d.py @@ -0,0 +1,120 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from . import data2d_c as d + +def clear(): + """Clear all 2D data objects.""" + return d.clear() + +def delete(ids): + """Delete a 2D object. + + :param ids: ID or list of ID's to delete. + :type ids: :class:`int` or :class:`list` + + """ + return d.delete(ids) + +def create(datatype, geometry, data=None): + """Create a 2D object. + + :param datatype: Data object type, '-vol' or '-sino'. + :type datatype: :class:`string` + :param geometry: Volume or projection geometry. + :type geometry: :class:`dict` + :param data: Data to fill the constructed object with, either a scalar or array. + :type data: :class:`float` or :class:`numpy.ndarray` + :returns: :class:`int` -- the ID of the constructed object. + + """ + return d.create(datatype,geometry,data) + +def store(i, data): + """Fill existing 2D object with data. + + :param i: ID of object to fill. + :type i: :class:`int` + :param data: Data to fill the object with, either a scalar or array. + :type data: :class:`float` or :class:`numpy.ndarray` + + """ + return d.store(i, data) + +def get_geometry(i): + """Get the geometry of a 2D object. + + :param i: ID of object. + :type i: :class:`int` + :returns: :class:`dict` -- The geometry of object with ID ``i``. + + """ + return d.get_geometry(i) + +def change_geometry(i, geom): + """Change the geometry of a 2D object. + + :param i: ID of object. + :type i: :class:`int` + :param geom: new geometry. + :type geom: :class:`dict` + + """ + return d.change_geometry(i, geom) + +def get(i): + """Get a 2D object. + + :param i: ID of object to get. + :type i: :class:`int` + :returns: :class:`numpy.ndarray` -- The object data. + + """ + return d.get(i) + +def get_shared(i): + """Get a 2D object with memory shared between the ASTRA toolbox and numpy array. + + :param i: ID of object to get. + :type i: :class:`int` + :returns: :class:`numpy.ndarray` -- The object data. + + """ + return d.get_shared(i) + + +def get_single(i): + """Get a 2D object in single precision. + + :param i: ID of object to get. + :type i: :class:`int` + :returns: :class:`numpy.ndarray` -- The object data. + + """ + return d.get_single(i) + +def info(): + """Print info on 2D objects in memory.""" + return d.info() diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx new file mode 100644 index 0000000..ec1b478 --- /dev/null +++ b/python/astra/data2d_c.pyx @@ -0,0 +1,241 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +import six + +cimport cython +from cython cimport view + +cimport PyData2DManager +from .PyData2DManager cimport CData2DManager + +cimport PyXMLDocument +from .PyXMLDocument cimport XMLDocument + +import numpy as np + +cimport numpy as np +np.import_array() + + +from .PyIncludes cimport * +cimport utils +from .utils import wrap_from_bytes + + +cdef CData2DManager * man2d = <CData2DManager * >PyData2DManager.getSingletonPtr() + + +def clear(): + man2d.clear() + + +def delete(ids): + try: + for i in ids: + man2d.remove(i) + except TypeError: + man2d.remove(ids) + + +def create(datatype, geometry, data=None): + cdef XMLDocument * xml + cdef Config cfg + cdef CVolumeGeometry2D * pGeometry + cdef CProjectionGeometry2D * ppGeometry + cdef CFloat32Data2D * pDataObject2D + if datatype == '-vol': + xml = utils.dict2XML(six.b('VolumeGeometry'), geometry) + cfg.self = xml.getRootNode() + pGeometry = new CVolumeGeometry2D() + if not pGeometry.initialize(cfg): + del xml + del pGeometry + raise Exception('Geometry class not initialized.') + pDataObject2D = <CFloat32Data2D * > new CFloat32VolumeData2D(pGeometry) + del xml + del pGeometry + elif datatype == '-sino': + xml = utils.dict2XML(six.b('ProjectionGeometry'), geometry) + cfg.self = xml.getRootNode() + tpe = wrap_from_bytes(cfg.self.getAttribute(six.b('type'))) + if (tpe == 'sparse_matrix'): + ppGeometry = <CProjectionGeometry2D * >new CSparseMatrixProjectionGeometry2D() + elif (tpe == 'fanflat'): + ppGeometry = <CProjectionGeometry2D * >new CFanFlatProjectionGeometry2D() + elif (tpe == 'fanflat_vec'): + ppGeometry = <CProjectionGeometry2D * >new CFanFlatVecProjectionGeometry2D() + else: + ppGeometry = <CProjectionGeometry2D * >new CParallelProjectionGeometry2D() + if not ppGeometry.initialize(cfg): + del xml + del ppGeometry + raise Exception('Geometry class not initialized.') + pDataObject2D = <CFloat32Data2D * > new CFloat32ProjectionData2D(ppGeometry) + del ppGeometry + del xml + else: + raise Exception("Invalid datatype. Please specify '-vol' or '-sino'.") + + if not pDataObject2D.isInitialized(): + del pDataObject2D + raise Exception("Couldn't initialize data object.") + + fillDataObject(pDataObject2D, data) + + return man2d.store(pDataObject2D) + +cdef fillDataObject(CFloat32Data2D * obj, data): + if data is None: + fillDataObjectScalar(obj, 0) + else: + if isinstance(data, np.ndarray): + fillDataObjectArray(obj, np.ascontiguousarray(data,dtype=np.float32)) + else: + fillDataObjectScalar(obj, np.float32(data)) + +cdef fillDataObjectScalar(CFloat32Data2D * obj, float s): + cdef int i + for i in range(obj.getSize()): + obj.getData()[i] = s + +@cython.boundscheck(False) +@cython.wraparound(False) +cdef fillDataObjectArray(CFloat32Data2D * obj, float [:,::1] data): + if (not data.shape[0] == obj.getHeight()) or (not data.shape[1] == obj.getWidth()): + raise Exception( + "The dimensions of the data do not match those specified in the geometry.") + cdef float [:,::1] cView = <float[:data.shape[0],:data.shape[1]]> obj.getData2D()[0] + cView[:] = data + +cdef CFloat32Data2D * getObject(i) except NULL: + cdef CFloat32Data2D * pDataObject = man2d.get(i) + if pDataObject == NULL: + raise Exception("Data object not found") + if not pDataObject.isInitialized(): + raise Exception("Data object not initialized properly.") + return pDataObject + + +def store(i, data): + cdef CFloat32Data2D * pDataObject = getObject(i) + fillDataObject(pDataObject, data) + + +def get_geometry(i): + cdef CFloat32Data2D * pDataObject = getObject(i) + cdef CFloat32ProjectionData2D * pDataObject2 + cdef CFloat32VolumeData2D * pDataObject3 + if pDataObject.getType() == PROJECTION: + pDataObject2 = <CFloat32ProjectionData2D * >pDataObject + geom = utils.createProjectionGeometryStruct(pDataObject2.getGeometry()) + elif pDataObject.getType() == VOLUME: + pDataObject3 = <CFloat32VolumeData2D * >pDataObject + geom = utils.createVolumeGeometryStruct(pDataObject3.getGeometry()) + else: + raise Exception("Not a known data object") + return geom + + +def change_geometry(i, geom): + cdef XMLDocument * xml + cdef Config cfg + cdef CVolumeGeometry2D * pGeometry + cdef CProjectionGeometry2D * ppGeometry + cdef CFloat32Data2D * pDataObject = getObject(i) + cdef CFloat32ProjectionData2D * pDataObject2 + cdef CFloat32VolumeData2D * pDataObject3 + if pDataObject.getType() == PROJECTION: + pDataObject2 = <CFloat32ProjectionData2D * >pDataObject + xml = utils.dict2XML(six.b('ProjectionGeometry'), geom) + cfg.self = xml.getRootNode() + tpe = wrap_from_bytes(cfg.self.getAttribute(six.b('type'))) + if (tpe == 'sparse_matrix'): + ppGeometry = <CProjectionGeometry2D * >new CSparseMatrixProjectionGeometry2D() + elif (tpe == 'fanflat'): + ppGeometry = <CProjectionGeometry2D * >new CFanFlatProjectionGeometry2D() + elif (tpe == 'fanflat_vec'): + ppGeometry = <CProjectionGeometry2D * >new CFanFlatVecProjectionGeometry2D() + else: + ppGeometry = <CProjectionGeometry2D * >new CParallelProjectionGeometry2D() + if not ppGeometry.initialize(cfg): + del xml + del ppGeometry + raise Exception('Geometry class not initialized.') + if (ppGeometry.getDetectorCount() != pDataObject2.getDetectorCount() or ppGeometry.getProjectionAngleCount() != pDataObject2.getAngleCount()): + del ppGeometry + del xml + raise Exception( + "The dimensions of the data do not match those specified in the geometry.") + pDataObject2.changeGeometry(ppGeometry) + del ppGeometry + del xml + elif pDataObject.getType() == VOLUME: + pDataObject3 = <CFloat32VolumeData2D * >pDataObject + xml = utils.dict2XML(six.b('VolumeGeometry'), geom) + cfg.self = xml.getRootNode() + pGeometry = new CVolumeGeometry2D() + if not pGeometry.initialize(cfg): + del xml + del pGeometry + raise Exception('Geometry class not initialized.') + if (pGeometry.getGridColCount() != pDataObject3.getWidth() or pGeometry.getGridRowCount() != pDataObject3.getHeight()): + del xml + del pGeometry + raise Exception( + 'The dimensions of the data do not match those specified in the geometry.') + pDataObject3.changeGeometry(pGeometry) + del xml + del pGeometry + else: + raise Exception("Not a known data object") + +@cython.boundscheck(False) +@cython.wraparound(False) +def get(i): + cdef CFloat32Data2D * pDataObject = getObject(i) + outArr = np.empty((pDataObject.getHeight(), pDataObject.getWidth()),dtype=np.float32,order='C') + cdef float [:,::1] mView = outArr + cdef float [:,::1] cView = <float[:outArr.shape[0],:outArr.shape[1]]> pDataObject.getData2D()[0] + mView[:] = cView + return outArr + +def get_shared(i): + cdef CFloat32Data2D * pDataObject = getObject(i) + cdef np.npy_intp shape[2] + shape[0] = <np.npy_intp> pDataObject.getHeight() + shape[1] = <np.npy_intp> pDataObject.getWidth() + return np.PyArray_SimpleNewFromData(2,shape,np.NPY_FLOAT32,<void *>pDataObject.getData2D()[0]) + + +def get_single(i): + raise Exception("Not yet implemented") + + +def info(): + six.print_(wrap_from_bytes(man2d.info())) diff --git a/python/astra/data3d.py b/python/astra/data3d.py new file mode 100644 index 0000000..33bde51 --- /dev/null +++ b/python/astra/data3d.py @@ -0,0 +1,108 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from . import data3d_c as d + +def create(datatype,geometry,data=None): + """Create a 3D object. + + :param datatype: Data object type, '-vol' or '-sino'. + :type datatype: :class:`string` + :param geometry: Volume or projection geometry. + :type geometry: :class:`dict` + :param data: Data to fill the constructed object with, either a scalar or array. + :type data: :class:`float` or :class:`numpy.ndarray` + :returns: :class:`int` -- the ID of the constructed object. + + """ + return d.create(datatype,geometry,data) + +def get(i): + """Get a 3D object. + + :param i: ID of object to get. + :type i: :class:`int` + :returns: :class:`numpy.ndarray` -- The object data. + + """ + return d.get(i) + +def get_shared(i): + """Get a 3D object with memory shared between the ASTRA toolbox and numpy array. + + :param i: ID of object to get. + :type i: :class:`int` + :returns: :class:`numpy.ndarray` -- The object data. + + """ + return d.get_shared(i) + +def get_single(i): + """Get a 3D object in single precision. + + :param i: ID of object to get. + :type i: :class:`int` + :returns: :class:`numpy.ndarray` -- The object data. + + """ + return g.get_single(i) + +def store(i,data): + """Fill existing 3D object with data. + + :param i: ID of object to fill. + :type i: :class:`int` + :param data: Data to fill the object with, either a scalar or array. + :type data: :class:`float` or :class:`numpy.ndarray` + + """ + return d.store(i,data) + +def dimensions(i): + """Get dimensions of a 3D object. + + :param i: ID of object. + :type i: :class:`int` + :returns: :class:`tuple` -- dimensions of object with ID ``i``. + + """ + return d.dimensions(i) + +def delete(ids): + """Delete a 2D object. + + :param ids: ID or list of ID's to delete. + :type ids: :class:`int` or :class:`list` + + """ + return d.delete(ids) + +def clear(): + """Clear all 3D data objects.""" + return d.clear() + +def info(): + """Print info on 3D objects in memory.""" + return d.info() diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx new file mode 100644 index 0000000..f821aaf --- /dev/null +++ b/python/astra/data3d_c.pyx @@ -0,0 +1,188 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +import six + +cimport cython + +cimport PyData3DManager +from .PyData3DManager cimport CData3DManager + +from .PyIncludes cimport * +import numpy as np + +cimport numpy as np +np.import_array() + +cimport PyXMLDocument +from .PyXMLDocument cimport XMLDocument + +cimport utils +from .utils import wrap_from_bytes + +cdef CData3DManager * man3d = <CData3DManager * >PyData3DManager.getSingletonPtr() + +cdef extern from *: + CFloat32Data3DMemory * dynamic_cast_mem "dynamic_cast<astra::CFloat32Data3DMemory*>" (CFloat32Data3D * ) except NULL + +def create(datatype,geometry,data=None): + cdef XMLDocument * xml + cdef Config cfg + cdef CVolumeGeometry3D * pGeometry + cdef CProjectionGeometry3D * ppGeometry + cdef CFloat32Data3DMemory * pDataObject3D + cdef CConeProjectionGeometry3D* pppGeometry + if datatype == '-vol': + xml = utils.dict2XML(six.b('VolumeGeometry'), geometry) + cfg.self = xml.getRootNode() + pGeometry = new CVolumeGeometry3D() + if not pGeometry.initialize(cfg): + del xml + del pGeometry + raise Exception('Geometry class not initialized.') + pDataObject3D = <CFloat32Data3DMemory * > new CFloat32VolumeData3DMemory(pGeometry) + del xml + del pGeometry + elif datatype == '-sino' or datatype == '-proj3d': + xml = utils.dict2XML(six.b('ProjectionGeometry'), geometry) + cfg.self = xml.getRootNode() + tpe = wrap_from_bytes(cfg.self.getAttribute(six.b('type'))) + if (tpe == "parallel3d"): + ppGeometry = <CProjectionGeometry3D*> new CParallelProjectionGeometry3D(); + elif (tpe == "parallel3d_vec"): + ppGeometry = <CProjectionGeometry3D*> new CParallelVecProjectionGeometry3D(); + elif (tpe == "cone"): + ppGeometry = <CProjectionGeometry3D*> new CConeProjectionGeometry3D(); + elif (tpe == "cone_vec"): + ppGeometry = <CProjectionGeometry3D*> new CConeVecProjectionGeometry3D(); + else: + raise Exception("Invalid geometry type.") + + if not ppGeometry.initialize(cfg): + del xml + del ppGeometry + raise Exception('Geometry class not initialized.') + pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(ppGeometry) + del ppGeometry + del xml + elif datatype == "-sinocone": + xml = utils.dict2XML(six.b('ProjectionGeometry'), geometry) + cfg.self = xml.getRootNode() + pppGeometry = new CConeProjectionGeometry3D() + if not pppGeometry.initialize(cfg): + del xml + del pppGeometry + raise Exception('Geometry class not initialized.') + pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(pppGeometry) + else: + raise Exception("Invalid datatype. Please specify '-vol' or '-proj3d'.") + + if not pDataObject3D.isInitialized(): + del pDataObject3D + raise Exception("Couldn't initialize data object.") + + fillDataObject(pDataObject3D, data) + + pDataObject3D.updateStatistics() + + return man3d.store(<CFloat32Data3D*>pDataObject3D) + + +cdef fillDataObject(CFloat32Data3DMemory * obj, data): + if data is None: + fillDataObjectScalar(obj, 0) + else: + if isinstance(data, np.ndarray): + fillDataObjectArray(obj, np.ascontiguousarray(data,dtype=np.float32)) + else: + fillDataObjectScalar(obj, np.float32(data)) + +cdef fillDataObjectScalar(CFloat32Data3DMemory * obj, float s): + cdef int i + for i in range(obj.getSize()): + obj.getData()[i] = s + +@cython.boundscheck(False) +@cython.wraparound(False) +cdef fillDataObjectArray(CFloat32Data3DMemory * obj, float [:,:,::1] data): + if (not data.shape[0] == obj.getDepth()) or (not data.shape[1] == obj.getHeight()) or (not data.shape[2] == obj.getWidth()): + raise Exception( + "The dimensions of the data do not match those specified in the geometry.") + cdef float [:,:,::1] cView = <float[:data.shape[0],:data.shape[1],:data.shape[2]]> obj.getData3D()[0][0] + cView[:] = data + +cdef CFloat32Data3D * getObject(i) except NULL: + cdef CFloat32Data3D * pDataObject = man3d.get(i) + if pDataObject == NULL: + raise Exception("Data object not found") + if not pDataObject.isInitialized(): + raise Exception("Data object not initialized properly.") + return pDataObject + +@cython.boundscheck(False) +@cython.wraparound(False) +def get(i): + cdef CFloat32Data3DMemory * pDataObject = dynamic_cast_mem(getObject(i)) + outArr = np.empty((pDataObject.getDepth(),pDataObject.getHeight(), pDataObject.getWidth()),dtype=np.float32,order='C') + cdef float [:,:,::1] mView = outArr + cdef float [:,:,::1] cView = <float[:outArr.shape[0],:outArr.shape[1],:outArr.shape[2]]> pDataObject.getData3D()[0][0] + mView[:] = cView + return outArr + +def get_shared(i): + cdef CFloat32Data3DMemory * pDataObject = dynamic_cast_mem(getObject(i)) + outArr = np.empty((pDataObject.getDepth(),pDataObject.getHeight(), pDataObject.getWidth()),dtype=np.float32,order='C') + cdef np.npy_intp shape[3] + shape[0] = <np.npy_intp> pDataObject.getDepth() + shape[1] = <np.npy_intp> pDataObject.getHeight() + shape[2] = <np.npy_intp> pDataObject.getWidth() + return np.PyArray_SimpleNewFromData(3,shape,np.NPY_FLOAT32,<void *>pDataObject.getData3D()[0][0]) + +def get_single(i): + raise Exception("Not yet implemented") + +def store(i,data): + cdef CFloat32Data3D * pDataObject = getObject(i) + fillDataObject(dynamic_cast_mem(pDataObject), data) + +def dimensions(i): + cdef CFloat32Data3D * pDataObject = getObject(i) + return (pDataObject.getWidth(),pDataObject.getHeight(),pDataObject.getDepth()) + +def delete(ids): + try: + for i in ids: + man3d.remove(i) + except TypeError: + man3d.remove(ids) + +def clear(): + man3d.clear() + +def info(): + six.print_(wrap_from_bytes(man3d.info())) diff --git a/python/astra/extrautils.pyx b/python/astra/extrautils.pyx new file mode 100644 index 0000000..998f5cb --- /dev/null +++ b/python/astra/extrautils.pyx @@ -0,0 +1,43 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- + +def clipCircle(img): + cdef int i,j + cdef double x2,y2,mid,bnd + cdef long sz,sz2 + sz = img.shape[0] + sz2 = sz*sz + bnd = sz2/4. + mid = (sz-1.)/2. + nDel=0 + for i in range(sz): + for j in range(sz): + x2 = (i-mid)*(i-mid) + y2 = (j-mid)*(j-mid) + if x2+y2>bnd: + img[i,j]=0 + nDel=nDel+1 + return nDel diff --git a/python/astra/functions.py b/python/astra/functions.py new file mode 100644 index 0000000..4025468 --- /dev/null +++ b/python/astra/functions.py @@ -0,0 +1,270 @@ +#----------------------------------------------------------------------- +# Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +# Author: Daniel M. Pelt +# Contact: D.M.Pelt@cwi.nl +# Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +# This file is part of the Python interface to the +# All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +# The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +# The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +"""Additional functions for PyAstraToolbox. + +.. moduleauthor:: Daniel M. Pelt <D.M.Pelt@cwi.nl> + + +""" + +from . import creators as ac +import numpy as np +from six.moves import range + +from . import data2d +from . import data3d +from . import projector +from . import algorithm + + + +def clear(): + """Clears all used memory of the ASTRA Toolbox. + + .. note:: + This is irreversible. + + """ + data2d.clear() + data3d.clear() + projector.clear() + algorithm.clear() + + +def data_op(op, data, scalar, gpu_core, mask=None): + """Perform data operation on data. + + :param op: Operation to perform. + :param data: Data to perform operation on. + :param scalar: Scalar argument to data operation. + :param gpu_core: GPU core to perform operation on. + :param mask: Optional mask. + + """ + + cfg = ac.astra_dict('DataOperation_CUDA') + cfg['Operation'] = op + cfg['Scalar'] = scalar + cfg['DataId'] = data + if not mask == None: + cfg['MaskId'] = mask + cfg['option']['GPUindex'] = gpu_core + alg_id = algorithm.create(cfg) + algorithm.run(alg_id) + algorithm.delete(alg_id) + + +def add_noise_to_sino(sinogram_in, I0, seed=None): + """Adds Poisson noise to a sinogram. + + :param sinogram_in: Sinogram to add noise to. + :type sinogram_in: :class:`numpy.ndarray` + :param I0: Background intensity. Lower values lead to higher noise. + :type I0: :class:`float` + :returns: :class:`numpy.ndarray` -- the sinogram with added noise. + + """ + + if not seed==None: + curstate = np.random.get_state() + np.random.seed(seed) + + if isinstance(sinogram_in, np.ndarray): + sinogramRaw = sinogram_in + else: + sinogramRaw = data2d.get(sinogram_in) + max_sinogramRaw = sinogramRaw.max() + sinogramRawScaled = sinogramRaw / max_sinogramRaw + # to detector count + sinogramCT = I0 * np.exp(-sinogramRawScaled) + # add poison noise + sinogramCT_C = np.zeros_like(sinogramCT) + for i in range(sinogramCT_C.shape[0]): + for j in range(sinogramCT_C.shape[1]): + sinogramCT_C[i, j] = np.random.poisson(sinogramCT[i, j]) + # to density + sinogramCT_D = sinogramCT_C / I0 + sinogram_out = -max_sinogramRaw * np.log(sinogramCT_D) + + if not isinstance(sinogram_in, np.ndarray): + at.data2d.store(sinogram_in, sinogram_out) + + if not seed==None: + np.random.set_state(curstate) + + return sinogram_out + +def move_vol_geom(geom, pos, is_relative=False): + """Moves center of volume geometry to new position. + + :param geom: Input volume geometry + :type geom: :class:`dict` + :param pos: Tuple (x,y[,z]) for new position, with the center of the image at (0,0[,0]) + :type pos: :class:`tuple` + :param is_relative: Whether new position is relative to the old position + :type is_relative: :class:`bool` + :returns: :class:`dict` -- Volume geometry with the new center + """ + + ret_geom = geom.copy() + ret_geom['option'] = geom['option'].copy() + + if not is_relative: + ret_geom['option']['WindowMinX'] = -geom['GridColCount']/2. + ret_geom['option']['WindowMaxX'] = geom['GridColCount']/2. + ret_geom['option']['WindowMinY'] = -geom['GridRowCount']/2. + ret_geom['option']['WindowMaxY'] = geom['GridRowCount']/2. + if len(pos)>2: + ret_geom['option']['WindowMinZ'] = -geom['GridSliceCount']/2. + ret_geom['option']['WindowMaxZ'] = geom['GridSliceCount']/2. + ret_geom['option']['WindowMinX'] += pos[0] + ret_geom['option']['WindowMaxX'] += pos[0] + ret_geom['option']['WindowMinY'] += pos[1] + ret_geom['option']['WindowMaxY'] += pos[1] + if len(pos)>2: + ret_geom['option']['WindowMinZ'] += pos[2] + ret_geom['option']['WindowMaxZ'] += pos[2] + return ret_geom + + +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 + + +def geom_2vec(proj_geom): + """Returns a vector-based projection geometry from a basic projection geometry. + + :param proj_geom: Projection geometry to convert + :type proj_geom: :class:`dict` + """ + if proj_geom['type'] == 'fanflat': + angles = proj_geom['ProjectionAngles'] + vectors = np.zeros((len(angles), 6)) + for i in range(len(angles)): + + # source + vectors[i, 0] = np.sin(angles[i]) * proj_geom['DistanceOriginSource'] + vectors[i, 1] = -np.cos(angles[i]) * proj_geom['DistanceOriginSource'] + + # center of detector + vectors[i, 2] = -np.sin(angles[i]) * proj_geom['DistanceOriginDetector'] + vectors[i, 3] = np.cos(angles[i]) * proj_geom['DistanceOriginDetector'] + + # vector from detector pixel 0 to 1 + vectors[i, 4] = np.cos(angles[i]) * proj_geom['DetectorWidth'] + vectors[i, 5] = np.sin(angles[i]) * proj_geom['DetectorWidth'] + proj_geom_out = ac.create_proj_geom( + 'fanflat_vec', proj_geom['DetectorCount'], vectors) + + elif proj_geom['type'] == 'cone': + angles = proj_geom['ProjectionAngles'] + vectors = np.zeros((len(angles), 12)) + for i in range(len(angles)): + # source + vectors[i, 0] = np.sin(angles[i]) * proj_geom['DistanceOriginSource'] + vectors[i, 1] = -np.cos(angles[i]) * proj_geom['DistanceOriginSource'] + vectors[i, 2] = 0 + + # center of detector + vectors[i, 3] = -np.sin(angles[i]) * proj_geom['DistanceOriginDetector'] + vectors[i, 4] = np.cos(angles[i]) * proj_geom['DistanceOriginDetector'] + vectors[i, 5] = 0 + + # vector from detector pixel (0,0) to (0,1) + vectors[i, 6] = np.cos(angles[i]) * proj_geom['DetectorSpacingX'] + vectors[i, 7] = np.sin(angles[i]) * proj_geom['DetectorSpacingX'] + vectors[i, 8] = 0 + + # vector from detector pixel (0,0) to (1,0) + vectors[i, 9] = 0 + vectors[i, 10] = 0 + vectors[i, 11] = proj_geom['DetectorSpacingY'] + + proj_geom_out = ac.create_proj_geom( + 'cone_vec', proj_geom['DetectorRowCount'], proj_geom['DetectorColCount'], vectors) + + # PARALLEL + elif proj_geom['type'] == 'parallel3d': + angles = proj_geom['ProjectionAngles'] + vectors = np.zeros((len(angles), 12)) + for i in range(len(angles)): + + # ray direction + vectors[i, 0] = np.sin(angles[i]) + vectors[i, 1] = -np.cos(angles[i]) + vectors[i, 2] = 0 + + # center of detector + vectors[i, 3] = 0 + vectors[i, 4] = 0 + vectors[i, 5] = 0 + + # vector from detector pixel (0,0) to (0,1) + vectors[i, 6] = np.cos(angles[i]) * proj_geom['DetectorSpacingX'] + vectors[i, 7] = np.sin(angles[i]) * proj_geom['DetectorSpacingX'] + vectors[i, 8] = 0 + + # vector from detector pixel (0,0) to (1,0) + vectors[i, 9] = 0 + vectors[i, 10] = 0 + vectors[i, 11] = proj_geom['DetectorSpacingY'] + + proj_geom_out = ac.create_proj_geom( + 'parallel3d_vec', proj_geom['DetectorRowCount'], proj_geom['DetectorColCount'], vectors) + + else: + raise ValueError( + 'No suitable vector geometry found for type: ' + proj_geom['type']) + return proj_geom_out diff --git a/python/astra/matlab.py b/python/astra/matlab.py new file mode 100644 index 0000000..83b345d --- /dev/null +++ b/python/astra/matlab.py @@ -0,0 +1,112 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +"""This module implements a MATLAB-like interface to the ASTRA Toolbox. + +Note that all functions are called with a :class:`string` as the first +argument, specifying the operation to perform. This un-pythonic way +is used to make transitioning from MATLAB code to Python code easier, as +the MATLAB interface uses the same type of method calling. + +After an initial ``import astra``, these functions can be accessed in the +``astra.m`` module. + +""" + +from . import astra_c +from . import data2d_c +from . import data3d_c +from . import projector_c +from . import algorithm_c +from . import matrix_c +import numpy as np + + +def astra(command, *args): + """MATLAB-like interface to the :mod:`astra.astra` module + + For example: + + ``astra.m.astra('use_cuda')`` -- Check if CUDA is enabled. + + """ + return getattr(astra_c, command)(*args) + + +def data2d(command, *args): + """MATLAB-like interface to the :mod:`astra.data2d` module + + For example: + + ``astra.m.data2d('create',type,geometry,data)`` -- Create a 2D object. + + """ + return getattr(data2d_c, command)(*args) + + +def data3d(command, *args): + """MATLAB-like interface to the :mod:`astra.data3d` module + + For example: + + ``astra.m.data3d('get',i)`` -- Get 3D object data. + + """ + return getattr(data3d_c, command)(*args) + + +def projector(command, *args): + """MATLAB-like interface to the :mod:`astra.projector` module + + For example: + + ``astra.m.projector('volume_geometry',i)`` -- Get volume geometry. + + """ + return getattr(projector_c, command)(*args) + + +def matrix(command, *args): + """MATLAB-like interface to the :mod:`astra.matrix` module + + For example: + + ``astra.m.matrix('delete',i)`` -- Delete a matrix. + + """ + return getattr(matrix_c, command)(*args) + + +def algorithm(command, *args): + """MATLAB-like interface to the :mod:`astra.algorithm` module + + For example: + + ``astra.m.algorithm('run',i,1000)`` -- Run an algorithm with 1000 iterations. + + """ + if command == 'iterate': + command = 'run' + return getattr(algorithm_c, command)(*args) diff --git a/python/astra/matrix.py b/python/astra/matrix.py new file mode 100644 index 0000000..27e4823 --- /dev/null +++ b/python/astra/matrix.py @@ -0,0 +1,86 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from . import matrix_c as m + +def delete(ids): + """Delete a matrix object. + + :param ids: ID or list of ID's to delete. + :type ids: :class:`int` or :class:`list` + + """ + return m.delete(ids) + +def clear(): + """Clear all matrix objects.""" + return m.clear() + +def create(data): + """Create matrix object with data. + + :param data: Data to fill the created object with. + :type data: :class:`scipy.sparse.csr_matrix` + :returns: :class:`int` -- the ID of the constructed object. + + """ + return m.create(data) + + +def store(i,data): + """Fill existing matrix object with data. + + :param i: ID of object to fill. + :type i: :class:`int` + :param data: Data to fill the object with. + :type data: :class:`scipy.sparse.csr_matrix` + + """ + return m.store(i,data) + +def get_size(i): + """Get matrix dimensions. + + :param i: ID of object. + :type i: :class:`int` + :returns: :class:`tuple` -- matrix dimensions. + """ + return m.get_size(i) + +def get(i): + """Get a matrix object. + + :param i: ID of object to get. + :type i: :class:`int` + :returns: :class:`scipy.sparse.csr_matrix` -- The object data. + + """ + return m.get(i) + +def info(): + """Print info on matrix objects in memory.""" + return m.info() + + diff --git a/python/astra/matrix_c.pyx b/python/astra/matrix_c.pyx new file mode 100644 index 0000000..b0d8bc4 --- /dev/null +++ b/python/astra/matrix_c.pyx @@ -0,0 +1,116 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +import six +from six.moves import range +import numpy as np +import scipy.sparse as ss + +from libcpp cimport bool + +cimport PyMatrixManager +from .PyMatrixManager cimport CMatrixManager +from .PyIncludes cimport * +from .utils import wrap_from_bytes + +cdef CMatrixManager * manM = <CMatrixManager * >PyMatrixManager.getSingletonPtr() + + +def delete(ids): + try: + for i in ids: + manM.remove(i) + except TypeError: + manM.remove(ids) + +def clear(): + manM.clear() + +cdef int csr_matrix_to_astra(data,CSparseMatrix *mat) except -1: + if isinstance(data,ss.csr_matrix): + csrD = data + else: + csrD = data.tocsr() + if not mat.isInitialized(): + raise Exception("Couldn't initialize data object.") + if csrD.nnz > mat.m_lSize or csrD.shape[0] > mat.m_iHeight: + raise Exception("Matrix too large to store in this object.") + for i in range(len(csrD.indptr)): + mat.m_plRowStarts[i] = csrD.indptr[i] + for i in range(csrD.nnz): + mat.m_piColIndices[i] = csrD.indices[i] + mat.m_pfValues[i] = csrD.data[i] + +cdef astra_to_csr_matrix(CSparseMatrix *mat): + indptr = np.zeros(mat.m_iHeight+1,dtype=np.int) + indices = np.zeros(mat.m_plRowStarts[mat.m_iHeight],dtype=np.int) + data = np.zeros(mat.m_plRowStarts[mat.m_iHeight]) + for i in range(mat.m_iHeight+1): + indptr[i] = mat.m_plRowStarts[i] + for i in range(mat.m_plRowStarts[mat.m_iHeight]): + indices[i] = mat.m_piColIndices[i] + data[i] = mat.m_pfValues[i] + return ss.csr_matrix((data,indices,indptr),shape=(mat.m_iHeight,mat.m_iWidth)) + +def create(data): + cdef CSparseMatrix* pMatrix + pMatrix = new CSparseMatrix(data.shape[0], data.shape[1], data.nnz) + if not pMatrix.isInitialized(): + del pMatrix + raise Exception("Couldn't initialize data object.") + try: + csr_matrix_to_astra(data,pMatrix) + except: + del pMatrix + raise Exception("Failed to create data object.") + + return manM.store(pMatrix) + +cdef CSparseMatrix * getObject(i) except NULL: + cdef CSparseMatrix * pDataObject = manM.get(i) + if pDataObject == NULL: + raise Exception("Data object not found") + if not pDataObject.isInitialized(): + raise Exception("Data object not initialized properly.") + return pDataObject + + +def store(i,data): + cdef CSparseMatrix * pDataObject = getObject(i) + csr_matrix_to_astra(data,pDataObject) + +def get_size(i): + cdef CSparseMatrix * pDataObject = getObject(i) + return (pDataObject.m_iHeight,pDataObject.m_iWidth) + +def get(i): + cdef CSparseMatrix * pDataObject = getObject(i) + return astra_to_csr_matrix(pDataObject) + +def info(): + six.print_(wrap_from_bytes(manM.info())) diff --git a/python/astra/projector.py b/python/astra/projector.py new file mode 100644 index 0000000..c916c52 --- /dev/null +++ b/python/astra/projector.py @@ -0,0 +1,100 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from . import 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) + + +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 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/projector_c.pyx b/python/astra/projector_c.pyx new file mode 100644 index 0000000..978ca09 --- /dev/null +++ b/python/astra/projector_c.pyx @@ -0,0 +1,119 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +import six +from .PyIncludes cimport * + +cimport utils +from .utils import wrap_from_bytes + +cimport PyProjector2DFactory +from .PyProjector2DFactory cimport CProjector2DFactory + +cimport PyProjector2DManager +from .PyProjector2DManager cimport CProjector2DManager + +cimport PyXMLDocument +from .PyXMLDocument cimport XMLDocument + +cimport PyMatrixManager +from .PyMatrixManager cimport CMatrixManager + +cdef CProjector2DManager * manProj = <CProjector2DManager * >PyProjector2DManager.getSingletonPtr() +cdef CMatrixManager * manM = <CMatrixManager * >PyMatrixManager.getSingletonPtr() + + +def create(config): + cdef XMLDocument * xml = utils.dict2XML(six.b('Projector2D'), config) + cdef Config cfg + cdef CProjector2D * proj + cfg.self = xml.getRootNode() + proj = PyProjector2DFactory.getSingletonPtr().create(cfg) + if proj == NULL: + del xml + raise Exception("Error creating projector.") + del xml + 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 CProjector2D * getObject(i) except NULL: + cdef CProjector2D * 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 CProjector2D * proj = getObject(i) + return utils.createProjectionGeometryStruct(proj.getProjectionGeometry()) + + +def volume_geometry(i): + cdef CProjector2D * proj = getObject(i) + return utils.createVolumeGeometryStruct(proj.getVolumeGeometry()) + + +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 matrix(i): + cdef CProjector2D * proj = getObject(i) + cdef CSparseMatrix * mat = proj.getMatrix() + if mat == NULL: + del mat + raise Exception("Data object not found") + if not mat.isInitialized(): + del mat + raise Exception("Data object not initialized properly.") + return manM.store(mat) diff --git a/python/astra/utils.pxd b/python/astra/utils.pxd new file mode 100644 index 0000000..55db9d3 --- /dev/null +++ b/python/astra/utils.pxd @@ -0,0 +1,37 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +from libcpp.string cimport string + +cimport PyXMLDocument +from .PyXMLDocument cimport XMLDocument +from .PyXMLDocument cimport XMLNode + +from .PyIncludes cimport * + +cdef XMLDocument *dict2XML(string rootname, dc) +cdef XML2dict(XMLDocument *) +cdef createVolumeGeometryStruct(CVolumeGeometry2D* geom) +cdef createProjectionGeometryStruct(CProjectionGeometry2D* geom) diff --git a/python/astra/utils.pyx b/python/astra/utils.pyx new file mode 100644 index 0000000..53e84a9 --- /dev/null +++ b/python/astra/utils.pyx @@ -0,0 +1,260 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +#----------------------------------------------------------------------- +# distutils: language = c++ +# distutils: libraries = astra + +import numpy as np +import six +from libcpp.string cimport string +from libcpp.list cimport list +from libcpp.vector cimport vector +from cython.operator cimport dereference as deref +from cpython.version cimport PY_MAJOR_VERSION + +cimport PyXMLDocument +from .PyXMLDocument cimport XMLDocument +from .PyXMLDocument cimport XMLNode +from .PyIncludes cimport * + + +cdef XMLDocument * dict2XML(string rootname, dc): + cdef XMLDocument * doc = PyXMLDocument.createDocument(rootname) + cdef XMLNode * node = doc.getRootNode() + try: + readDict(node, dc) + except: + six.print_('Error reading XML') + del doc + doc = NULL + finally: + del node + return doc + +def convert_item(item): + if isinstance(item, six.string_types): + return item.encode('ascii') + + if type(item) is not dict: + return item + + out_dict = {} + for k in item: + out_dict[convert_item(k)] = convert_item(item[k]) + return out_dict + + +def wrap_to_bytes(value): + if isinstance(value, six.binary_type): + return value + s = str(value) + if PY_MAJOR_VERSION == 3: + s = s.encode('ascii') + return s + + +def wrap_from_bytes(value): + s = value + if PY_MAJOR_VERSION == 3: + s = s.decode('ascii') + return s + + +cdef void readDict(XMLNode * root, _dc): + cdef XMLNode * listbase + cdef XMLNode * itm + cdef int i + cdef int j + + dc = convert_item(_dc) + for item in dc: + val = dc[item] + if isinstance(val, np.ndarray): + if val.size == 0: + break + listbase = root.addChildNode(item) + listbase.addAttribute(< string > six.b('listsize'), < float32 > val.size) + index = 0 + if val.ndim == 2: + for i in range(val.shape[0]): + for j in range(val.shape[1]): + itm = listbase.addChildNode(six.b('ListItem')) + itm.addAttribute(< string > six.b('index'), < float32 > index) + itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) + index += 1 + del itm + 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'), <string> wrap_to_bytes(val)) + else: + itm = root.addChildNode(item, wrap_to_bytes(val)) + del itm + +cdef void readOptions(XMLNode * node, dc): + cdef XMLNode * listbase + cdef XMLNode * itm + cdef int i + cdef int j + for item in dc: + val = dc[item] + if node.hasOption(item): + raise Exception('Duplicate Option: %s' % item) + if isinstance(val, np.ndarray): + if val.size == 0: + 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 + if val.ndim == 2: + for i in range(val.shape[0]): + for j in range(val.shape[1]): + itm = listbase.addChildNode(six.b('ListItem')) + itm.addAttribute(< string > six.b('index'), < float32 > index) + itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) + index += 1 + del itm + 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)) + +cdef vectorToNumpy(vector[float32] inp): + cdef int i + cdef int sz = inp.size() + ret = np.empty(sz) + for i in range(sz): + ret[i] = inp[i] + return ret + +cdef XMLNode2dict(XMLNode * node): + cdef XMLNode * subnode + cdef list[XMLNode * ] nodes + cdef list[XMLNode * ].iterator it + dct = {} + if node.hasAttribute(six.b('type')): + dct['type'] = node.getAttribute(six.b('type')) + nodes = node.getNodes() + it = nodes.begin() + while it != nodes.end(): + subnode = deref(it) + if subnode.hasAttribute(six.b('listsize')): + dct[subnode.getName( + )] = vectorToNumpy(subnode.getContentNumericalArray()) + else: + dct[subnode.getName()] = subnode.getContent() + del subnode + return dct + +cdef XML2dict(XMLDocument * xml): + cdef XMLNode * node = xml.getRootNode() + dct = XMLNode2dict(node) + del node; + return dct; + +cdef createProjectionGeometryStruct(CProjectionGeometry2D * geom): + cdef int i + cdef CFanFlatVecProjectionGeometry2D * fanvecGeom + # cdef SFanProjection* p + dct = {} + dct['DetectorCount'] = geom.getDetectorCount() + if not geom.isOfType(< string > six.b('fanflat_vec')): + dct['DetectorWidth'] = geom.getDetectorWidth() + angles = np.empty(geom.getProjectionAngleCount()) + for i in range(geom.getProjectionAngleCount()): + angles[i] = geom.getProjectionAngle(i) + dct['ProjectionAngles'] = angles + else: + raise Exception("Not yet implemented") + # fanvecGeom = <CFanFlatVecProjectionGeometry2D*> geom + # vecs = np.empty(fanvecGeom.getProjectionAngleCount()*6) + # iDetCount = pVecGeom.getDetectorCount() + # for i in range(fanvecGeom.getProjectionAngleCount()): + # p = &fanvecGeom.getProjectionVectors()[i]; + # out[6*i + 0] = p.fSrcX + # out[6*i + 1] = p.fSrcY + # out[6*i + 2] = p.fDetSX + 0.5f*iDetCount*p.fDetUX + # out[6*i + 3] = p.fDetSY + 0.5f*iDetCount*p.fDetUY + # out[6*i + 4] = p.fDetUX + # out[6*i + 5] = p.fDetUY + # dct['Vectors'] = vecs + if (geom.isOfType(< string > six.b('parallel'))): + dct["type"] = "parallel" + elif (geom.isOfType(< string > six.b('fanflat'))): + raise Exception("Not yet implemented") + # astra::CFanFlatProjectionGeometry2D* pFanFlatGeom = dynamic_cast<astra::CFanFlatProjectionGeometry2D*>(_pProjGeom) + # mGeometryInfo["DistanceOriginSource"] = mxCreateDoubleScalar(pFanFlatGeom->getOriginSourceDistance()) + # mGeometryInfo["DistanceOriginDetector"] = + # mxCreateDoubleScalar(pFanFlatGeom->getOriginDetectorDistance()) + dct["type"] = "fanflat" + elif (geom.isOfType(< string > six.b('sparse_matrix'))): + raise Exception("Not yet implemented") + # astra::CSparseMatrixProjectionGeometry2D* pSparseMatrixGeom = + # dynamic_cast<astra::CSparseMatrixProjectionGeometry2D*>(_pProjGeom); + dct["type"] = "sparse_matrix" + # dct["MatrixID"] = + # mxCreateDoubleScalar(CMatrixManager::getSingleton().getIndex(pSparseMatrixGeom->getMatrix())) + elif(geom.isOfType(< string > six.b('fanflat_vec'))): + dct["type"] = "fanflat_vec" + return dct + +cdef createVolumeGeometryStruct(CVolumeGeometry2D * geom): + mGeometryInfo = {} + mGeometryInfo["GridColCount"] = geom.getGridColCount() + mGeometryInfo["GridRowCount"] = geom.getGridRowCount() + + mGeometryOptions = {} + mGeometryOptions["WindowMinX"] = geom.getWindowMinX() + mGeometryOptions["WindowMaxX"] = geom.getWindowMaxX() + mGeometryOptions["WindowMinY"] = geom.getWindowMinY() + mGeometryOptions["WindowMaxY"] = geom.getWindowMaxY() + + mGeometryInfo["option"] = mGeometryOptions + return mGeometryInfo |