From 1b32573046f33050b9300324e6c74e10abb6caaf Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 9 Apr 2015 15:44:01 +0200 Subject: Add 'link' feature to Python (for 2D and 3D data) --- python/astra/CFloat32CustomPython.h | 17 +++++++++++++++++ python/astra/PyIncludes.pxd | 8 ++++++++ python/astra/data2d.py | 21 +++++++++++++++++++++ python/astra/data2d_c.pyx | 21 +++++++++++++++++---- python/astra/data3d.py | 22 ++++++++++++++++++++++ python/astra/data3d_c.pyx | 27 ++++++++++++++++++++++----- 6 files changed, 107 insertions(+), 9 deletions(-) create mode 100644 python/astra/CFloat32CustomPython.h (limited to 'python/astra') diff --git a/python/astra/CFloat32CustomPython.h b/python/astra/CFloat32CustomPython.h new file mode 100644 index 0000000..d8593fc --- /dev/null +++ b/python/astra/CFloat32CustomPython.h @@ -0,0 +1,17 @@ +class CFloat32CustomPython : public astra::CFloat32CustomMemory { +public: + CFloat32CustomPython(PyObject * arrIn) + { + arr = arrIn; + // Set pointer to numpy data pointer + m_fPtr = (float *)PyArray_DATA(arr); + // Increase reference count since ASTRA has a reference + Py_INCREF(arr); + } + virtual ~CFloat32CustomPython() { + // Decrease reference count since ASTRA object is destroyed + Py_DECREF(arr); + } +private: + PyObject* arr; +}; \ No newline at end of file diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 7df02c5..13329d1 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -63,10 +63,14 @@ cdef extern from "astra/VolumeGeometry2D.h" namespace "astra": float32 getWindowMaxY() Config* getConfiguration() +cdef extern from "astra/Float32Data2D.h" namespace "astra": + cdef cppclass CFloat32CustomMemory: + pass cdef extern from "astra/Float32VolumeData2D.h" namespace "astra": cdef cppclass CFloat32VolumeData2D: CFloat32VolumeData2D(CVolumeGeometry2D*) + CFloat32VolumeData2D(CVolumeGeometry2D*, CFloat32CustomMemory*) CVolumeGeometry2D * getGeometry() int getWidth() int getHeight() @@ -130,6 +134,7 @@ cdef extern from "astra/ParallelProjectionGeometry2D.h" namespace "astra": cdef extern from "astra/Float32ProjectionData2D.h" namespace "astra": cdef cppclass CFloat32ProjectionData2D: CFloat32ProjectionData2D(CProjectionGeometry2D*) + CFloat32ProjectionData2D(CProjectionGeometry2D*, CFloat32CustomMemory*) CProjectionGeometry2D * getGeometry() void changeGeometry(CProjectionGeometry2D*) int getDetectorCount() @@ -207,6 +212,7 @@ cdef extern from "astra/ProjectionGeometry3D.h" namespace "astra": cdef extern from "astra/Float32VolumeData3DMemory.h" namespace "astra": cdef cppclass CFloat32VolumeData3DMemory: CFloat32VolumeData3DMemory(CVolumeGeometry3D*) + CFloat32VolumeData3DMemory(CVolumeGeometry3D*, CFloat32CustomMemory*) CVolumeGeometry3D* getGeometry() @@ -231,6 +237,8 @@ cdef extern from "astra/Float32ProjectionData3DMemory.h" namespace "astra": cdef cppclass CFloat32ProjectionData3DMemory: CFloat32ProjectionData3DMemory(CProjectionGeometry3D*) CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*) + CFloat32ProjectionData3DMemory(CProjectionGeometry3D*, CFloat32CustomMemory*) + CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*, CFloat32CustomMemory*) CProjectionGeometry3D* getGeometry() cdef extern from "astra/Float32Data3D.h" namespace "astra": diff --git a/python/astra/data2d.py b/python/astra/data2d.py index 8c4be03..f119f05 100644 --- a/python/astra/data2d.py +++ b/python/astra/data2d.py @@ -24,6 +24,7 @@ # #----------------------------------------------------------------------- from . import data2d_c as d +import numpy as np def clear(): """Clear all 2D data objects.""" @@ -52,6 +53,26 @@ def create(datatype, geometry, data=None): """ return d.create(datatype,geometry,data) +def link(datatype, geometry, data): + """Link a 2D numpy array with the toolbox. + + :param datatype: Data object type, '-vol' or '-sino'. + :type datatype: :class:`string` + :param geometry: Volume or projection geometry. + :type geometry: :class:`dict` + :param data: Numpy array to link + :type data: :class:`numpy.ndarray` + :returns: :class:`int` -- the ID of the constructed object. + + """ + if not isinstance(data,np.ndarray): + raise ValueError("Input should be a numpy array") + if not data.dtype==np.float32: + raise ValueError("Numpy array should be float32") + if not (data.flags['C_CONTIGUOUS'] and data.flags['ALIGNED']): + raise ValueError("Numpy array should be C_CONTIGUOUS and ALIGNED") + return d.create(datatype,geometry,data,True) + def store(i, data): """Fill existing 2D object with data. diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index b9c105e..ac54898 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -49,6 +49,10 @@ from .utils import wrap_from_bytes cdef CData2DManager * man2d = PyData2DManager.getSingletonPtr() +cdef extern from "CFloat32CustomPython.h": + cdef cppclass CFloat32CustomPython: + CFloat32CustomPython(arrIn) + def clear(): man2d.clear() @@ -61,11 +65,12 @@ def delete(ids): man2d.remove(ids) -def create(datatype, geometry, data=None): +def create(datatype, geometry, data=None, link=False): cdef Config *cfg cdef CVolumeGeometry2D * pGeometry cdef CProjectionGeometry2D * ppGeometry cdef CFloat32Data2D * pDataObject2D + cdef CFloat32CustomMemory * pCustom if datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) pGeometry = new CVolumeGeometry2D() @@ -73,7 +78,11 @@ def create(datatype, geometry, data=None): del cfg del pGeometry raise Exception('Geometry class not initialized.') - pDataObject2D = new CFloat32VolumeData2D(pGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject2D = new CFloat32VolumeData2D(pGeometry, pCustom) + else: + pDataObject2D = new CFloat32VolumeData2D(pGeometry) del cfg del pGeometry elif datatype == '-sino': @@ -91,7 +100,11 @@ def create(datatype, geometry, data=None): del cfg del ppGeometry raise Exception('Geometry class not initialized.') - pDataObject2D = new CFloat32ProjectionData2D(ppGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject2D = new CFloat32ProjectionData2D(ppGeometry, pCustom) + else: + pDataObject2D = new CFloat32ProjectionData2D(ppGeometry) del ppGeometry del cfg else: @@ -101,7 +114,7 @@ def create(datatype, geometry, data=None): del pDataObject2D raise Exception("Couldn't initialize data object.") - fillDataObject(pDataObject2D, data) + if not link: fillDataObject(pDataObject2D, data) return man2d.store(pDataObject2D) diff --git a/python/astra/data3d.py b/python/astra/data3d.py index a2e9201..4fdf9d7 100644 --- a/python/astra/data3d.py +++ b/python/astra/data3d.py @@ -24,6 +24,7 @@ # #----------------------------------------------------------------------- from . import data3d_c as d +import numpy as np def create(datatype,geometry,data=None): """Create a 3D object. @@ -39,6 +40,27 @@ def create(datatype,geometry,data=None): """ return d.create(datatype,geometry,data) +def link(datatype, geometry, data): + """Link a 3D numpy array with the toolbox. + + :param datatype: Data object type, '-vol' or '-sino'. + :type datatype: :class:`string` + :param geometry: Volume or projection geometry. + :type geometry: :class:`dict` + :param data: Numpy array to link + :type data: :class:`numpy.ndarray` + :returns: :class:`int` -- the ID of the constructed object. + + """ + if not isinstance(data,np.ndarray): + raise ValueError("Input should be a numpy array") + if not data.dtype==np.float32: + raise ValueError("Numpy array should be float32") + if not (data.flags['C_CONTIGUOUS'] and data.flags['ALIGNED']): + raise ValueError("Numpy array should be C_CONTIGUOUS and ALIGNED") + return d.create(datatype,geometry,data,True) + + def get(i): """Get a 3D object. diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 4b069f7..f2c6e26 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -50,12 +50,17 @@ cdef CData3DManager * man3d = PyData3DManager.getSingletonPtr cdef extern from *: CFloat32Data3DMemory * dynamic_cast_mem "dynamic_cast" (CFloat32Data3D * ) except NULL -def create(datatype,geometry,data=None): +cdef extern from "CFloat32CustomPython.h": + cdef cppclass CFloat32CustomPython: + CFloat32CustomPython(arrIn) + +def create(datatype,geometry,data=None, link=False): cdef Config *cfg cdef CVolumeGeometry3D * pGeometry cdef CProjectionGeometry3D * ppGeometry cdef CFloat32Data3DMemory * pDataObject3D cdef CConeProjectionGeometry3D* pppGeometry + cdef CFloat32CustomMemory * pCustom if datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) pGeometry = new CVolumeGeometry3D() @@ -63,7 +68,11 @@ def create(datatype,geometry,data=None): del cfg del pGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = new CFloat32VolumeData3DMemory(pGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject3D = new CFloat32VolumeData3DMemory(pGeometry, pCustom) + else: + pDataObject3D = new CFloat32VolumeData3DMemory(pGeometry) del cfg del pGeometry elif datatype == '-sino' or datatype == '-proj3d': @@ -84,7 +93,11 @@ def create(datatype,geometry,data=None): del cfg del ppGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = new CFloat32ProjectionData3DMemory(ppGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject3D = new CFloat32ProjectionData3DMemory(ppGeometry, pCustom) + else: + pDataObject3D = new CFloat32ProjectionData3DMemory(ppGeometry) del ppGeometry del cfg elif datatype == "-sinocone": @@ -94,7 +107,11 @@ def create(datatype,geometry,data=None): del cfg del pppGeometry raise Exception('Geometry class not initialized.') - pDataObject3D = new CFloat32ProjectionData3DMemory(pppGeometry) + if link: + pCustom = new CFloat32CustomPython(data) + pDataObject3D = new CFloat32ProjectionData3DMemory(pppGeometry, pCustom) + else: + pDataObject3D = new CFloat32ProjectionData3DMemory(pppGeometry) else: raise Exception("Invalid datatype. Please specify '-vol' or '-proj3d'.") @@ -102,7 +119,7 @@ def create(datatype,geometry,data=None): del pDataObject3D raise Exception("Couldn't initialize data object.") - fillDataObject(pDataObject3D, data) + if not link: fillDataObject(pDataObject3D, data) pDataObject3D.updateStatistics() -- cgit v1.2.3 From e2c87a5e259c847772c733eefb1b291b2a5b1a6e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 10 Apr 2015 14:29:27 +0200 Subject: Add python data3d.change_geometry --- python/astra/PyIncludes.pxd | 15 +++++++++++++ python/astra/data3d.py | 11 +++++++++ python/astra/data3d_c.pyx | 55 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 81 insertions(+) (limited to 'python/astra') diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 7df02c5..a581f88 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -196,18 +196,29 @@ cdef extern from "astra/VolumeGeometry3D.h" namespace "astra": CVolumeGeometry3D() bool initialize(Config) Config * getConfiguration() + int getGridColCount() + int getGridRowCount() + int getGridSliceCount() cdef extern from "astra/ProjectionGeometry3D.h" namespace "astra": cdef cppclass CProjectionGeometry3D: CProjectionGeometry3D() bool initialize(Config) Config * getConfiguration() + int getProjectionCount() + int getDetectorColCount() + int getDetectorRowCount() cdef extern from "astra/Float32VolumeData3DMemory.h" namespace "astra": cdef cppclass CFloat32VolumeData3DMemory: CFloat32VolumeData3DMemory(CVolumeGeometry3D*) CVolumeGeometry3D* getGeometry() + void changeGeometry(CVolumeGeometry3D*) + int getRowCount() + int getColCount() + int getSliceCount() + cdef extern from "astra/ParallelProjectionGeometry3D.h" namespace "astra": @@ -232,6 +243,10 @@ cdef extern from "astra/Float32ProjectionData3DMemory.h" namespace "astra": CFloat32ProjectionData3DMemory(CProjectionGeometry3D*) CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*) CProjectionGeometry3D* getGeometry() + void changeGeometry(CProjectionGeometry3D*) + int getDetectorColCount() + int getDetectorRowCount() + int getAngleCount() cdef extern from "astra/Float32Data3D.h" namespace "astra": cdef cppclass CFloat32Data3D: diff --git a/python/astra/data3d.py b/python/astra/data3d.py index a2e9201..4679489 100644 --- a/python/astra/data3d.py +++ b/python/astra/data3d.py @@ -90,6 +90,17 @@ def get_geometry(i): """ return d.get_geometry(i) +def change_geometry(i, geometry): + """Change the geometry of a 3D object. + + :param i: ID of object. + :type i: :class:`int` + :param geometry: Volume or projection geometry. + :type geometry: :class:`dict` + + """ + return d.change_geometry(i, geometry) + def dimensions(i): """Get dimensions of a 3D object. diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 4b069f7..48af032 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -122,6 +122,61 @@ def get_geometry(i): raise Exception("Not a known data object") return geom +def change_geometry(i, geom): + cdef CFloat32Data3DMemory * pDataObject = dynamic_cast_mem(getObject(i)) + cdef CFloat32ProjectionData3DMemory * pDataObject2 + cdef CFloat32VolumeData3DMemory * pDataObject3 + if pDataObject.getType() == THREEPROJECTION: + pDataObject2 = pDataObject + # TODO: Reduce code duplication here + cfg = utils.dictToConfig(six.b('ProjectionGeometry'), geom) + tpe = wrap_from_bytes(cfg.self.getAttribute(six.b('type'))) + if (tpe == "parallel3d"): + ppGeometry = new CParallelProjectionGeometry3D(); + elif (tpe == "parallel3d_vec"): + ppGeometry = new CParallelVecProjectionGeometry3D(); + elif (tpe == "cone"): + ppGeometry = new CConeProjectionGeometry3D(); + elif (tpe == "cone_vec"): + ppGeometry = new CConeVecProjectionGeometry3D(); + else: + raise Exception("Invalid geometry type.") + if not ppGeometry.initialize(cfg[0]): + del cfg + del ppGeometry + raise Exception('Geometry class not initialized.') + del cfg + if (ppGeometry.getDetectorColCount() != pDataObject2.getDetectorColCount() or \ + ppGeometry.getProjectionCount() != pDataObject2.getAngleCount() or \ + ppGeometry.getDetectorRowCount() != pDataObject2.getDetectorRowCount()): + del ppGeometry + raise Exception( + "The dimensions of the data do not match those specified in the geometry.") + pDataObject2.changeGeometry(ppGeometry) + del ppGeometry + + elif pDataObject.getType() == THREEVOLUME: + pDataObject3 = pDataObject + cfg = utils.dictToConfig(six.b('VolumeGeometry'), geom) + pGeometry = new CVolumeGeometry3D() + if not pGeometry.initialize(cfg[0]): + del cfg + del pGeometry + raise Exception('Geometry class not initialized.') + del cfg + if (pGeometry.getGridColCount() != pDataObject3.getColCount() or \ + pGeometry.getGridRowCount() != pDataObject3.getRowCount() or \ + pGeometry.getGridSliceCount() != pDataObject3.getSliceCount()): + del pGeometry + raise Exception( + "The dimensions of the data do not match those specified in the geometry.") + pDataObject3.changeGeometry(pGeometry) + del pGeometry + + else: + raise Exception("Not a known data object") + + cdef fillDataObject(CFloat32Data3DMemory * obj, data): if data is None: fillDataObjectScalar(obj, 0) -- cgit v1.2.3 From f69d9f6bc1704560518da3c30c46e495c0228aac Mon Sep 17 00:00:00 2001 From: Daan Pelt Date: Thu, 30 Apr 2015 11:02:50 +0200 Subject: Check data size when using 'link' function in Python --- python/astra/data2d_c.pyx | 16 ++++++++++++ python/astra/data3d_c.pyx | 17 ++++++++++++ python/astra/functions.py | 25 ++---------------- python/astra/pythonutils.py | 63 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 98 insertions(+), 23 deletions(-) create mode 100644 python/astra/pythonutils.py (limited to 'python/astra') diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index ac54898..29548b5 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -47,6 +47,12 @@ from .PyIncludes cimport * cimport utils from .utils import wrap_from_bytes +from .pythonutils import geom_size + +import operator + +from six.moves import reduce + cdef CData2DManager * man2d = PyData2DManager.getSingletonPtr() cdef extern from "CFloat32CustomPython.h": @@ -71,6 +77,16 @@ def create(datatype, geometry, data=None, link=False): cdef CProjectionGeometry2D * ppGeometry cdef CFloat32Data2D * pDataObject2D cdef CFloat32CustomMemory * pCustom + + if link: + geomSize = geom_size(geometry) + if len(data.shape)==1: + if data.size!=reduce(operator.mul,geomSize): + raise Exception("The dimensions of the data do not match those specified in the geometry.") + else: + if data.shape!=geomSize: + raise Exception("The dimensions of the data do not match those specified in the geometry.") + if datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) pGeometry = new CVolumeGeometry2D() diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 84472c1..30745b4 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -45,6 +45,13 @@ from .PyXMLDocument cimport XMLDocument cimport utils from .utils import wrap_from_bytes +from .pythonutils import geom_size + +import operator + +from six.moves import reduce + + cdef CData3DManager * man3d = PyData3DManager.getSingletonPtr() cdef extern from *: @@ -61,6 +68,16 @@ def create(datatype,geometry,data=None, link=False): cdef CFloat32Data3DMemory * pDataObject3D cdef CConeProjectionGeometry3D* pppGeometry cdef CFloat32CustomMemory * pCustom + + if link: + geomSize = geom_size(geometry) + if len(data.shape)==1: + if data.size!=reduce(operator.mul,geomSize): + raise Exception("The dimensions of the data do not match those specified in the geometry.") + else: + if data.shape!=geomSize: + raise Exception("The dimensions of the data do not match those specified in the geometry.") + if datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) pGeometry = new CVolumeGeometry3D() diff --git a/python/astra/functions.py b/python/astra/functions.py index 4025468..b826b86 100644 --- a/python/astra/functions.py +++ b/python/astra/functions.py @@ -38,6 +38,7 @@ from . import data2d from . import data3d from . import projector from . import algorithm +from . import pythonutils @@ -158,29 +159,7 @@ def geom_size(geom, dim=None): :param dim: Optional axis index to return :type dim: :class:`int` """ - - if 'GridSliceCount' in geom: - # 3D Volume geometry? - s = (geom['GridSliceCount'], geom[ - 'GridRowCount'], geom['GridColCount']) - elif 'GridColCount' in geom: - # 2D Volume geometry? - s = (geom['GridRowCount'], geom['GridColCount']) - elif geom['type'] == 'parallel' or geom['type'] == 'fanflat': - s = (len(geom['ProjectionAngles']), geom['DetectorCount']) - elif geom['type'] == 'parallel3d' or geom['type'] == 'cone': - s = (geom['DetectorRowCount'], len( - geom['ProjectionAngles']), geom['DetectorColCount']) - elif geom['type'] == 'fanflat_vec': - s = (geom['Vectors'].shape[0], geom['DetectorCount']) - elif geom['type'] == 'parallel3d_vec' or geom['type'] == 'cone_vec': - s = (geom['DetectorRowCount'], geom[ - 'Vectors'].shape[0], geom['DetectorColCount']) - - if dim != None: - s = s[dim] - - return s + return pythonutils.geom_size(geom,dim) def geom_2vec(proj_geom): diff --git a/python/astra/pythonutils.py b/python/astra/pythonutils.py new file mode 100644 index 0000000..8ea4af5 --- /dev/null +++ b/python/astra/pythonutils.py @@ -0,0 +1,63 @@ +#----------------------------------------------------------------------- +# Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +# Author: Daniel M. Pelt +# Contact: D.M.Pelt@cwi.nl +# Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +# This file is part of the Python interface to the +# All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +# The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +# The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with the Python interface to the ASTRA Toolbox. If not, see . +# +#----------------------------------------------------------------------- +"""Additional purely Python functions for PyAstraToolbox. + +.. moduleauthor:: Daniel M. Pelt + + +""" + +def geom_size(geom, dim=None): + """Returns the size of a volume or sinogram, based on the projection or volume geometry. + + :param geom: Geometry to calculate size from + :type geometry: :class:`dict` + :param dim: Optional axis index to return + :type dim: :class:`int` + """ + + if 'GridSliceCount' in geom: + # 3D Volume geometry? + s = (geom['GridSliceCount'], geom[ + 'GridRowCount'], geom['GridColCount']) + elif 'GridColCount' in geom: + # 2D Volume geometry? + s = (geom['GridRowCount'], geom['GridColCount']) + elif geom['type'] == 'parallel' or geom['type'] == 'fanflat': + s = (len(geom['ProjectionAngles']), geom['DetectorCount']) + elif geom['type'] == 'parallel3d' or geom['type'] == 'cone': + s = (geom['DetectorRowCount'], len( + geom['ProjectionAngles']), geom['DetectorColCount']) + elif geom['type'] == 'fanflat_vec': + s = (geom['Vectors'].shape[0], geom['DetectorCount']) + elif geom['type'] == 'parallel3d_vec' or geom['type'] == 'cone_vec': + s = (geom['DetectorRowCount'], geom[ + 'Vectors'].shape[0], geom['DetectorColCount']) + + if dim != None: + s = s[dim] + + return s -- cgit v1.2.3 From 2bc0d98c413fee4108115f26aa337f65337eec55 Mon Sep 17 00:00:00 2001 From: Daan Pelt Date: Thu, 26 Mar 2015 16:40:38 +0100 Subject: Add SPOT-like object for Python (overrides `__mul__` and works with scipy.sparse.linalg) --- python/astra/ASTRAProjector.py | 135 -------------------------- python/astra/__init__.py | 2 +- python/astra/operator.py | 208 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 209 insertions(+), 136 deletions(-) delete mode 100644 python/astra/ASTRAProjector.py create mode 100644 python/astra/operator.py (limited to 'python/astra') diff --git a/python/astra/ASTRAProjector.py b/python/astra/ASTRAProjector.py deleted file mode 100644 index f282618..0000000 --- a/python/astra/ASTRAProjector.py +++ /dev/null @@ -1,135 +0,0 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam -# -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ -# -# -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). -# -#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify -#it under the terms of the GNU General Public License as published by -#the Free Software Foundation, either version 3 of the License, or -#(at your option) any later version. -# -#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, -#but WITHOUT ANY WARRANTY; without even the implied warranty of -#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -#GNU General Public License for more details. -# -#You should have received a copy of the GNU General Public License -#along with the Python interface to the ASTRA Toolbox. If not, see . -# -#----------------------------------------------------------------------- - -import math -from . import creators as ac -from . import data2d - - -class ASTRAProjector2DTranspose(): - """Implements the ``proj.T`` functionality. - - Do not use directly, since it can be accessed as member ``.T`` of - an :class:`ASTRAProjector2D` object. - - """ - def __init__(self, parentProj): - self.parentProj = parentProj - - def __mul__(self, data): - return self.parentProj.backProject(data) - - -class ASTRAProjector2D(object): - """Helps with various common ASTRA Toolbox 2D operations. - - This class can perform several often used toolbox operations, such as: - - * Forward projecting - * Back projecting - * Reconstructing - - Note that this class has a some computational overhead, because it - copies a lot of data. If you use many repeated operations, directly - using the PyAstraToolbox methods directly is faster. - - You can use this class as an abstracted weight matrix :math:`W`: multiplying an instance - ``proj`` of this class by an image results in a forward projection of the image, and multiplying - ``proj.T`` by a sinogram results in a backprojection of the sinogram:: - - proj = ASTRAProjector2D(...) - fp = proj*image - bp = proj.T*sinogram - - :param proj_geom: The projection geometry. - :type proj_geom: :class:`dict` - :param vol_geom: The volume geometry. - :type vol_geom: :class:`dict` - :param proj_type: Projector type, such as ``'line'``, ``'linear'``, ... - :type proj_type: :class:`string` - """ - - def __init__(self, proj_geom, vol_geom, proj_type): - self.vol_geom = vol_geom - self.recSize = vol_geom['GridColCount'] - self.angles = proj_geom['ProjectionAngles'] - self.nDet = proj_geom['DetectorCount'] - nexpow = int(pow(2, math.ceil(math.log(2 * self.nDet, 2)))) - self.filterSize = nexpow / 2 + 1 - self.nProj = self.angles.shape[0] - self.proj_geom = proj_geom - self.proj_id = ac.create_projector(proj_type, proj_geom, vol_geom) - self.T = ASTRAProjector2DTranspose(self) - - def backProject(self, data): - """Backproject a sinogram. - - :param data: The sinogram data or ID. - :type data: :class:`numpy.ndarray` or :class:`int` - :returns: :class:`numpy.ndarray` -- The backprojection. - - """ - vol_id, vol = ac.create_backprojection( - data, self.proj_id, returnData=True) - data2d.delete(vol_id) - return vol - - def forwardProject(self, data): - """Forward project an image. - - :param data: The image data or ID. - :type data: :class:`numpy.ndarray` or :class:`int` - :returns: :class:`numpy.ndarray` -- The forward projection. - - """ - sin_id, sino = ac.create_sino(data, self.proj_id, returnData=True) - data2d.delete(sin_id) - return sino - - def reconstruct(self, data, method, **kwargs): - """Reconstruct an image from a sinogram. - - :param data: The sinogram data or ID. - :type data: :class:`numpy.ndarray` or :class:`int` - :param method: Name of the reconstruction algorithm. - :type method: :class:`string` - :param kwargs: Additional named parameters to pass to :func:`astra.creators.create_reconstruction`. - :returns: :class:`numpy.ndarray` -- The reconstruction. - - Example of a SIRT reconstruction using CUDA:: - - proj = ASTRAProjector2D(...) - rec = proj.reconstruct(sinogram,'SIRT_CUDA',iterations=1000) - - """ - kwargs['returnData'] = True - rec_id, rec = ac.create_reconstruction( - method, self.proj_id, data, **kwargs) - data2d.delete(rec_id) - return rec - - def __mul__(self, data): - return self.forwardProject(data) diff --git a/python/astra/__init__.py b/python/astra/__init__.py index 063dc16..8c1740c 100644 --- a/python/astra/__init__.py +++ b/python/astra/__init__.py @@ -27,7 +27,6 @@ from . import matlab as m from .creators import astra_dict,create_vol_geom, create_proj_geom, create_backprojection, create_sino, create_reconstruction, create_projector,create_sino3d_gpu, create_backprojection3d_gpu from .functions import data_op, add_noise_to_sino, clear, move_vol_geom from .extrautils import clipCircle -from .ASTRAProjector import ASTRAProjector2D from . import data2d from . import astra from . import data3d @@ -36,6 +35,7 @@ from . import projector from . import projector3d from . import matrix from . import log +from .operator import OpTomo import os try: diff --git a/python/astra/operator.py b/python/astra/operator.py new file mode 100644 index 0000000..a3abd5a --- /dev/null +++ b/python/astra/operator.py @@ -0,0 +1,208 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see . +# +#----------------------------------------------------------------------- + +from . import data2d +from . import data3d +from . import projector +from . import projector3d +from . import creators +from . import algorithm +from . import functions +import numpy as np +from six.moves import range, reduce +import operator +import scipy.sparse.linalg + +class OpTomo(scipy.sparse.linalg.LinearOperator): + """Object that imitates a projection matrix with a given projector. + + This object can do forward projection by using the ``*`` operator:: + + W = astra.OpTomo(proj_id) + fp = W*image + bp = W.T*sinogram + + It can also be used in minimization methods of the :mod:`scipy.sparse.linalg` module:: + + W = astra.OpTomo(proj_id) + output = scipy.sparse.linalg.lsqr(W,sinogram) + + :param proj_id: ID to a projector. + :type proj_id: :class:`int` + """ + + def __init__(self,proj_id): + self.dtype = np.float32 + try: + self.vg = projector.volume_geometry(proj_id) + self.pg = projector.projection_geometry(proj_id) + self.data_mod = data2d + self.appendString = "" + if projector.is_cuda(proj_id): + self.appendString += "_CUDA" + except Exception: + self.vg = projector3d.volume_geometry(proj_id) + self.pg = projector3d.projection_geometry(proj_id) + self.data_mod = data3d + self.appendString = "3D" + if projector3d.is_cuda(proj_id): + self.appendString += "_CUDA" + + self.vshape = functions.geom_size(self.vg) + self.vsize = reduce(operator.mul,self.vshape) + self.sshape = functions.geom_size(self.pg) + self.ssize = reduce(operator.mul,self.sshape) + + self.shape = (self.ssize, self.vsize) + + self.proj_id = proj_id + + self.T = OpTomoTranspose(self) + + def __checkArray(self, arr, shp): + if len(arr.shape)==1: + arr = arr.reshape(shp) + if arr.dtype != np.float32: + arr = arr.astype(np.float32) + if arr.flags['C_CONTIGUOUS']==False: + arr = np.ascontiguousarray(arr) + return arr + + def matvec(self,v): + """Implements the forward operator. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + v = self.__checkArray(v, self.vshape) + vid = self.data_mod.link('-vol',self.vg,v) + s = np.zeros(self.sshape,dtype=np.float32) + sid = self.data_mod.link('-sino',self.pg,s) + + cfg = creators.astra_dict('FP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['VolumeDataId'] = vid + cfg['ProjectorId'] = self.proj_id + fp_id = algorithm.create(cfg) + algorithm.run(fp_id) + + algorithm.delete(fp_id) + self.data_mod.delete([vid,sid]) + return s.flatten() + + def rmatvec(self,s): + """Implements the transpose operator. + + :param s: The projection data. + :type s: :class:`numpy.ndarray` + """ + s = self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + v = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,v) + + cfg = creators.astra_dict('BP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + bp_id = algorithm.create(cfg) + algorithm.run(bp_id) + + algorithm.delete(bp_id) + self.data_mod.delete([vid,sid]) + return v.flatten() + + def matmat(self,m): + """Implements the forward operator with a matrix. + + :param m: Volumes to forward project, arranged in columns. + :type m: :class:`numpy.ndarray` + """ + out = np.zeros((self.ssize,m.shape[1]),dtype=np.float32) + for i in range(m.shape[1]): + out[:,i] = self.matvec(m[:,i].flatten()) + return out + + def __mul__(self,v): + """Provides easy forward operator by *. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + return self.matvec(v) + + def reconstruct(self, method, s, iterations=1, extraOptions = {}): + """Reconstruct an object. + + :param method: Method to use for reconstruction. + :type method: :class:`string` + :param s: The projection data. + :type s: :class:`numpy.ndarray` + :param iterations: Number of iterations to use. + :type iterations: :class:`int` + :param extraOptions: Extra options to use during reconstruction (i.e. for cfg['option']). + :type extraOptions: :class:`dict` + """ + self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + v = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,v) + cfg = creators.astra_dict(method) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + cfg['option'] = extraOptions + alg_id = algorithm.create(cfg) + algorithm.run(alg_id,iterations) + algorithm.delete(alg_id) + self.data_mod.delete([vid,sid]) + return v + +class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): + """This object provides the transpose operation (``.T``) of the OpTomo object. + + Do not use directly, since it can be accessed as member ``.T`` of + an :class:`OpTomo` object. + """ + def __init__(self,parent): + self.parent = parent + self.dtype = np.float32 + self.shape = (parent.shape[1], parent.shape[0]) + + def matvec(self, s): + return self.parent.rmatvec(s) + + def rmatvec(self, v): + return self.parent.matvec(v) + + def matmat(self, m): + out = np.zeros((self.vsize,m.shape[1]),dtype=np.float32) + for i in range(m.shape[1]): + out[:,i] = self.matvec(m[:,i].flatten()) + return out + + def __mul__(self,v): + return self.matvec(v) -- cgit v1.2.3 From 47fe3421585302f2101691a685ab99b0e1ad5cfc Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 1 May 2015 17:48:32 +0200 Subject: Change XMLNode* to XMLNode An XMLNode object is already simply a pointer, so no need to dynamically allocate XMLNodes. --- python/astra/PyIncludes.pxd | 2 +- python/astra/PyXMLDocument.pxd | 8 ++++---- python/astra/utils.pyx | 29 ++++++++++------------------- 3 files changed, 15 insertions(+), 24 deletions(-) (limited to 'python/astra') diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 1d8285b..909f58f 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -43,7 +43,7 @@ cdef extern from "astra/Config.h" namespace "astra": cdef cppclass Config: Config() void initialize(string rootname) - XMLNode *self + XMLNode self cdef extern from "astra/VolumeGeometry2D.h" namespace "astra": cdef cppclass CVolumeGeometry2D: diff --git a/python/astra/PyXMLDocument.pxd b/python/astra/PyXMLDocument.pxd index 69781f1..57c447e 100644 --- a/python/astra/PyXMLDocument.pxd +++ b/python/astra/PyXMLDocument.pxd @@ -44,14 +44,14 @@ cdef extern from "astra/Globals.h" namespace "astra": cdef extern from "astra/XMLNode.h" namespace "astra": cdef cppclass XMLNode: string getName() - XMLNode *addChildNode(string name) - XMLNode *addChildNode(string, string) + XMLNode addChildNode(string name) + XMLNode addChildNode(string, string) void addAttribute(string, string) void addAttribute(string, float32) void addOption(string, string) bool hasOption(string) string getAttribute(string) - list[XMLNode *] getNodes() + list[XMLNode] getNodes() vector[float32] getContentNumericalArray() string getContent() bool hasAttribute(string) @@ -59,7 +59,7 @@ cdef extern from "astra/XMLNode.h" namespace "astra": cdef extern from "astra/XMLDocument.h" namespace "astra": cdef cppclass XMLDocument: void saveToFile(string sFilename) - XMLNode *getRootNode() + XMLNode getRootNode() cdef extern from "astra/XMLDocument.h" namespace "astra::XMLDocument": cdef XMLDocument *createDocument(string rootname) diff --git a/python/astra/utils.pyx b/python/astra/utils.pyx index 0439f1b..8f1e0b7 100644 --- a/python/astra/utils.pyx +++ b/python/astra/utils.pyx @@ -80,9 +80,9 @@ def wrap_from_bytes(value): return s -cdef void readDict(XMLNode * root, _dc): - cdef XMLNode * listbase - cdef XMLNode * itm +cdef void readDict(XMLNode root, _dc): + cdef XMLNode listbase + cdef XMLNode itm cdef int i cdef int j @@ -102,34 +102,29 @@ cdef void readDict(XMLNode * root, _dc): itm.addAttribute(< string > six.b('index'), < float32 > index) itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) index += 1 - del itm elif val.ndim == 1: for i in range(val.shape[0]): itm = listbase.addChildNode(six.b('ListItem')) itm.addAttribute(< string > six.b('index'), < float32 > index) itm.addAttribute(< string > six.b('value'), < float32 > val[i]) index += 1 - del itm else: raise Exception("Only 1 or 2 dimensions are allowed") - del listbase elif isinstance(val, dict): if item == six.b('option') or item == six.b('options') or item == six.b('Option') or item == six.b('Options'): readOptions(root, val) else: itm = root.addChildNode(item) readDict(itm, val) - del itm else: if item == six.b('type'): root.addAttribute(< string > six.b('type'), wrap_to_bytes(val)) else: itm = root.addChildNode(item, wrap_to_bytes(val)) - del itm -cdef void readOptions(XMLNode * node, dc): - cdef XMLNode * listbase - cdef XMLNode * itm +cdef void readOptions(XMLNode node, dc): + cdef XMLNode listbase + cdef XMLNode itm cdef int i cdef int j for item in dc: @@ -150,17 +145,14 @@ cdef void readOptions(XMLNode * node, dc): itm.addAttribute(< string > six.b('index'), < float32 > index) itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) index += 1 - del itm elif val.ndim == 1: for i in range(val.shape[0]): itm = listbase.addChildNode(six.b('ListItem')) itm.addAttribute(< string > six.b('index'), < float32 > index) itm.addAttribute(< string > six.b('value'), < float32 > val[i]) index += 1 - del itm else: raise Exception("Only 1 or 2 dimensions are allowed") - del listbase else: node.addOption(item, wrap_to_bytes(val)) @@ -214,10 +206,10 @@ def stringToPythonValue(inputIn): return str(input) -cdef XMLNode2dict(XMLNode * node): - cdef XMLNode * subnode - cdef list[XMLNode * ] nodes - cdef list[XMLNode * ].iterator it +cdef XMLNode2dict(XMLNode node): + cdef XMLNode subnode + cdef list[XMLNode] nodes + cdef list[XMLNode].iterator it dct = {} opts = {} if node.hasAttribute(six.b('type')): @@ -230,7 +222,6 @@ cdef XMLNode2dict(XMLNode * node): opts[castString(subnode.getAttribute('key'))] = stringToPythonValue(subnode.getAttribute('value')) else: dct[castString(subnode.getName())] = stringToPythonValue(subnode.getContent()) - del subnode inc(it) if len(opts)>0: dct['options'] = opts return dct -- cgit v1.2.3 From fff7470f1d74b0085355130350fa834ea8d37069 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 6 May 2015 13:50:11 +0200 Subject: Make XML array handling consistent setContent and getContent were using different XML formats previously. --- python/astra/PyXMLDocument.pxd | 2 ++ python/astra/utils.pyx | 35 +++++++++-------------------------- 2 files changed, 11 insertions(+), 26 deletions(-) (limited to 'python/astra') diff --git a/python/astra/PyXMLDocument.pxd b/python/astra/PyXMLDocument.pxd index 57c447e..033b8ef 100644 --- a/python/astra/PyXMLDocument.pxd +++ b/python/astra/PyXMLDocument.pxd @@ -53,6 +53,8 @@ cdef extern from "astra/XMLNode.h" namespace "astra": string getAttribute(string) list[XMLNode] getNodes() vector[float32] getContentNumericalArray() + void setContent(double*, int, int, bool) + void setContent(double*, int) string getContent() bool hasAttribute(string) diff --git a/python/astra/utils.pyx b/python/astra/utils.pyx index 8f1e0b7..ddb37aa 100644 --- a/python/astra/utils.pyx +++ b/python/astra/utils.pyx @@ -26,6 +26,7 @@ # distutils: language = c++ # distutils: libraries = astra +cimport numpy as np import numpy as np import six from libcpp.string cimport string @@ -85,6 +86,7 @@ cdef void readDict(XMLNode root, _dc): cdef XMLNode itm cdef int i cdef int j + cdef double* data dc = convert_item(_dc) for item in dc: @@ -93,21 +95,11 @@ cdef void readDict(XMLNode root, _dc): if val.size == 0: break listbase = root.addChildNode(item) - listbase.addAttribute(< string > six.b('listsize'), < float32 > val.size) - index = 0 + data = np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) if val.ndim == 2: - for i in range(val.shape[0]): - for j in range(val.shape[1]): - itm = listbase.addChildNode(six.b('ListItem')) - itm.addAttribute(< string > six.b('index'), < float32 > index) - itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) - index += 1 + listbase.setContent(data, val.shape[1], val.shape[0], False) elif val.ndim == 1: - for i in range(val.shape[0]): - itm = listbase.addChildNode(six.b('ListItem')) - itm.addAttribute(< string > six.b('index'), < float32 > index) - itm.addAttribute(< string > six.b('value'), < float32 > val[i]) - index += 1 + listbase.setContent(data, val.shape[0]) else: raise Exception("Only 1 or 2 dimensions are allowed") elif isinstance(val, dict): @@ -127,6 +119,7 @@ cdef void readOptions(XMLNode node, dc): cdef XMLNode itm cdef int i cdef int j + cdef double* data for item in dc: val = dc[item] if node.hasOption(item): @@ -136,21 +129,11 @@ cdef void readOptions(XMLNode node, dc): break listbase = node.addChildNode(six.b('Option')) listbase.addAttribute(< string > six.b('key'), < string > item) - listbase.addAttribute(< string > six.b('listsize'), < float32 > val.size) - index = 0 + data = np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) if val.ndim == 2: - for i in range(val.shape[0]): - for j in range(val.shape[1]): - itm = listbase.addChildNode(six.b('ListItem')) - itm.addAttribute(< string > six.b('index'), < float32 > index) - itm.addAttribute( < string > six.b('value'), < float32 > val[i, j]) - index += 1 + listbase.setContent(data, val.shape[1], val.shape[0], False) elif val.ndim == 1: - for i in range(val.shape[0]): - itm = listbase.addChildNode(six.b('ListItem')) - itm.addAttribute(< string > six.b('index'), < float32 > index) - itm.addAttribute(< string > six.b('value'), < float32 > val[i]) - index += 1 + listbase.setContent(data, val.shape[0]) else: raise Exception("Only 1 or 2 dimensions are allowed") else: -- cgit v1.2.3 From f730efe78367e8fe8e589c2b43fb0886d384f5c8 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 7 May 2015 11:34:37 +0200 Subject: Do not allow 1D input in Python link method --- python/astra/data2d_c.pyx | 10 ++-------- python/astra/data3d_c.pyx | 10 ++-------- 2 files changed, 4 insertions(+), 16 deletions(-) (limited to 'python/astra') diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index 29548b5..4919bf2 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -78,14 +78,8 @@ def create(datatype, geometry, data=None, link=False): cdef CFloat32Data2D * pDataObject2D cdef CFloat32CustomMemory * pCustom - if link: - geomSize = geom_size(geometry) - if len(data.shape)==1: - if data.size!=reduce(operator.mul,geomSize): - raise Exception("The dimensions of the data do not match those specified in the geometry.") - else: - if data.shape!=geomSize: - raise Exception("The dimensions of the data do not match those specified in the geometry.") + if link and data.shape!=geom_size(geometry): + raise Exception("The dimensions of the data do not match those specified in the geometry.") if datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 30745b4..3b27ab7 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -69,14 +69,8 @@ def create(datatype,geometry,data=None, link=False): cdef CConeProjectionGeometry3D* pppGeometry cdef CFloat32CustomMemory * pCustom - if link: - geomSize = geom_size(geometry) - if len(data.shape)==1: - if data.size!=reduce(operator.mul,geomSize): - raise Exception("The dimensions of the data do not match those specified in the geometry.") - else: - if data.shape!=geomSize: - raise Exception("The dimensions of the data do not match those specified in the geometry.") + if link and data.shape!=geom_size(geometry): + raise Exception("The dimensions of the data do not match those specified in the geometry.") if datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) -- cgit v1.2.3 From b9b9c82f8634f9c77416de5e857d107005cccbdf Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 7 May 2015 15:40:17 +0200 Subject: Use superclass __mul__ in Python OpTomo __mul__ --- python/astra/operator.py | 33 +++++++++++---------------------- 1 file changed, 11 insertions(+), 22 deletions(-) (limited to 'python/astra') diff --git a/python/astra/operator.py b/python/astra/operator.py index a3abd5a..0c37353 100644 --- a/python/astra/operator.py +++ b/python/astra/operator.py @@ -91,7 +91,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): arr = np.ascontiguousarray(arr) return arr - def matvec(self,v): + def _matvec(self,v): """Implements the forward operator. :param v: Volume to forward project. @@ -135,24 +135,16 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): self.data_mod.delete([vid,sid]) return v.flatten() - def matmat(self,m): - """Implements the forward operator with a matrix. - - :param m: Volumes to forward project, arranged in columns. - :type m: :class:`numpy.ndarray` - """ - out = np.zeros((self.ssize,m.shape[1]),dtype=np.float32) - for i in range(m.shape[1]): - out[:,i] = self.matvec(m[:,i].flatten()) - return out - def __mul__(self,v): """Provides easy forward operator by *. :param v: Volume to forward project. :type v: :class:`numpy.ndarray` """ - return self.matvec(v) + # Catch the case of a forward projection of a 2D/3D image + if isinstance(v, np.ndarray) and v.shape==self.vshape: + return self._matvec(v) + return scipy.sparse.linalg.LinearOperator.__mul__(self, v) def reconstruct(self, method, s, iterations=1, extraOptions = {}): """Reconstruct an object. @@ -192,17 +184,14 @@ class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): self.dtype = np.float32 self.shape = (parent.shape[1], parent.shape[0]) - def matvec(self, s): + def _matvec(self, s): return self.parent.rmatvec(s) def rmatvec(self, v): return self.parent.matvec(v) - def matmat(self, m): - out = np.zeros((self.vsize,m.shape[1]),dtype=np.float32) - for i in range(m.shape[1]): - out[:,i] = self.matvec(m[:,i].flatten()) - return out - - def __mul__(self,v): - return self.matvec(v) + def __mul__(self,s): + # Catch the case of a backprojection of 2D/3D data + if isinstance(s, np.ndarray) and s.shape==self.parent.sshape: + return self._matvec(s) + return scipy.sparse.linalg.LinearOperator.__mul__(self, s) -- cgit v1.2.3 From 89da933904262d6b7e80e8adf85ca9d1273881b3 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Fri, 8 May 2015 13:48:39 +0200 Subject: Rename optomo.py for Python 2 support --- python/astra/__init__.py | 2 +- python/astra/operator.py | 197 ----------------------------------------------- python/astra/optomo.py | 197 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 198 insertions(+), 198 deletions(-) delete mode 100644 python/astra/operator.py create mode 100644 python/astra/optomo.py (limited to 'python/astra') diff --git a/python/astra/__init__.py b/python/astra/__init__.py index 8c1740c..6c15d30 100644 --- a/python/astra/__init__.py +++ b/python/astra/__init__.py @@ -35,7 +35,7 @@ from . import projector from . import projector3d from . import matrix from . import log -from .operator import OpTomo +from .optomo import OpTomo import os try: diff --git a/python/astra/operator.py b/python/astra/operator.py deleted file mode 100644 index 0c37353..0000000 --- a/python/astra/operator.py +++ /dev/null @@ -1,197 +0,0 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam -# -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ -# -# -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). -# -#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify -#it under the terms of the GNU General Public License as published by -#the Free Software Foundation, either version 3 of the License, or -#(at your option) any later version. -# -#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, -#but WITHOUT ANY WARRANTY; without even the implied warranty of -#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -#GNU General Public License for more details. -# -#You should have received a copy of the GNU General Public License -#along with the Python interface to the ASTRA Toolbox. If not, see . -# -#----------------------------------------------------------------------- - -from . import data2d -from . import data3d -from . import projector -from . import projector3d -from . import creators -from . import algorithm -from . import functions -import numpy as np -from six.moves import range, reduce -import operator -import scipy.sparse.linalg - -class OpTomo(scipy.sparse.linalg.LinearOperator): - """Object that imitates a projection matrix with a given projector. - - This object can do forward projection by using the ``*`` operator:: - - W = astra.OpTomo(proj_id) - fp = W*image - bp = W.T*sinogram - - It can also be used in minimization methods of the :mod:`scipy.sparse.linalg` module:: - - W = astra.OpTomo(proj_id) - output = scipy.sparse.linalg.lsqr(W,sinogram) - - :param proj_id: ID to a projector. - :type proj_id: :class:`int` - """ - - def __init__(self,proj_id): - self.dtype = np.float32 - try: - self.vg = projector.volume_geometry(proj_id) - self.pg = projector.projection_geometry(proj_id) - self.data_mod = data2d - self.appendString = "" - if projector.is_cuda(proj_id): - self.appendString += "_CUDA" - except Exception: - self.vg = projector3d.volume_geometry(proj_id) - self.pg = projector3d.projection_geometry(proj_id) - self.data_mod = data3d - self.appendString = "3D" - if projector3d.is_cuda(proj_id): - self.appendString += "_CUDA" - - self.vshape = functions.geom_size(self.vg) - self.vsize = reduce(operator.mul,self.vshape) - self.sshape = functions.geom_size(self.pg) - self.ssize = reduce(operator.mul,self.sshape) - - self.shape = (self.ssize, self.vsize) - - self.proj_id = proj_id - - self.T = OpTomoTranspose(self) - - def __checkArray(self, arr, shp): - if len(arr.shape)==1: - arr = arr.reshape(shp) - if arr.dtype != np.float32: - arr = arr.astype(np.float32) - if arr.flags['C_CONTIGUOUS']==False: - arr = np.ascontiguousarray(arr) - return arr - - def _matvec(self,v): - """Implements the forward operator. - - :param v: Volume to forward project. - :type v: :class:`numpy.ndarray` - """ - v = self.__checkArray(v, self.vshape) - vid = self.data_mod.link('-vol',self.vg,v) - s = np.zeros(self.sshape,dtype=np.float32) - sid = self.data_mod.link('-sino',self.pg,s) - - cfg = creators.astra_dict('FP'+self.appendString) - cfg['ProjectionDataId'] = sid - cfg['VolumeDataId'] = vid - cfg['ProjectorId'] = self.proj_id - fp_id = algorithm.create(cfg) - algorithm.run(fp_id) - - algorithm.delete(fp_id) - self.data_mod.delete([vid,sid]) - return s.flatten() - - def rmatvec(self,s): - """Implements the transpose operator. - - :param s: The projection data. - :type s: :class:`numpy.ndarray` - """ - s = self.__checkArray(s, self.sshape) - sid = self.data_mod.link('-sino',self.pg,s) - v = np.zeros(self.vshape,dtype=np.float32) - vid = self.data_mod.link('-vol',self.vg,v) - - cfg = creators.astra_dict('BP'+self.appendString) - cfg['ProjectionDataId'] = sid - cfg['ReconstructionDataId'] = vid - cfg['ProjectorId'] = self.proj_id - bp_id = algorithm.create(cfg) - algorithm.run(bp_id) - - algorithm.delete(bp_id) - self.data_mod.delete([vid,sid]) - return v.flatten() - - def __mul__(self,v): - """Provides easy forward operator by *. - - :param v: Volume to forward project. - :type v: :class:`numpy.ndarray` - """ - # Catch the case of a forward projection of a 2D/3D image - if isinstance(v, np.ndarray) and v.shape==self.vshape: - return self._matvec(v) - return scipy.sparse.linalg.LinearOperator.__mul__(self, v) - - def reconstruct(self, method, s, iterations=1, extraOptions = {}): - """Reconstruct an object. - - :param method: Method to use for reconstruction. - :type method: :class:`string` - :param s: The projection data. - :type s: :class:`numpy.ndarray` - :param iterations: Number of iterations to use. - :type iterations: :class:`int` - :param extraOptions: Extra options to use during reconstruction (i.e. for cfg['option']). - :type extraOptions: :class:`dict` - """ - self.__checkArray(s, self.sshape) - sid = self.data_mod.link('-sino',self.pg,s) - v = np.zeros(self.vshape,dtype=np.float32) - vid = self.data_mod.link('-vol',self.vg,v) - cfg = creators.astra_dict(method) - cfg['ProjectionDataId'] = sid - cfg['ReconstructionDataId'] = vid - cfg['ProjectorId'] = self.proj_id - cfg['option'] = extraOptions - alg_id = algorithm.create(cfg) - algorithm.run(alg_id,iterations) - algorithm.delete(alg_id) - self.data_mod.delete([vid,sid]) - return v - -class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): - """This object provides the transpose operation (``.T``) of the OpTomo object. - - Do not use directly, since it can be accessed as member ``.T`` of - an :class:`OpTomo` object. - """ - def __init__(self,parent): - self.parent = parent - self.dtype = np.float32 - self.shape = (parent.shape[1], parent.shape[0]) - - def _matvec(self, s): - return self.parent.rmatvec(s) - - def rmatvec(self, v): - return self.parent.matvec(v) - - def __mul__(self,s): - # Catch the case of a backprojection of 2D/3D data - if isinstance(s, np.ndarray) and s.shape==self.parent.sshape: - return self._matvec(s) - return scipy.sparse.linalg.LinearOperator.__mul__(self, s) diff --git a/python/astra/optomo.py b/python/astra/optomo.py new file mode 100644 index 0000000..0c37353 --- /dev/null +++ b/python/astra/optomo.py @@ -0,0 +1,197 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see . +# +#----------------------------------------------------------------------- + +from . import data2d +from . import data3d +from . import projector +from . import projector3d +from . import creators +from . import algorithm +from . import functions +import numpy as np +from six.moves import range, reduce +import operator +import scipy.sparse.linalg + +class OpTomo(scipy.sparse.linalg.LinearOperator): + """Object that imitates a projection matrix with a given projector. + + This object can do forward projection by using the ``*`` operator:: + + W = astra.OpTomo(proj_id) + fp = W*image + bp = W.T*sinogram + + It can also be used in minimization methods of the :mod:`scipy.sparse.linalg` module:: + + W = astra.OpTomo(proj_id) + output = scipy.sparse.linalg.lsqr(W,sinogram) + + :param proj_id: ID to a projector. + :type proj_id: :class:`int` + """ + + def __init__(self,proj_id): + self.dtype = np.float32 + try: + self.vg = projector.volume_geometry(proj_id) + self.pg = projector.projection_geometry(proj_id) + self.data_mod = data2d + self.appendString = "" + if projector.is_cuda(proj_id): + self.appendString += "_CUDA" + except Exception: + self.vg = projector3d.volume_geometry(proj_id) + self.pg = projector3d.projection_geometry(proj_id) + self.data_mod = data3d + self.appendString = "3D" + if projector3d.is_cuda(proj_id): + self.appendString += "_CUDA" + + self.vshape = functions.geom_size(self.vg) + self.vsize = reduce(operator.mul,self.vshape) + self.sshape = functions.geom_size(self.pg) + self.ssize = reduce(operator.mul,self.sshape) + + self.shape = (self.ssize, self.vsize) + + self.proj_id = proj_id + + self.T = OpTomoTranspose(self) + + def __checkArray(self, arr, shp): + if len(arr.shape)==1: + arr = arr.reshape(shp) + if arr.dtype != np.float32: + arr = arr.astype(np.float32) + if arr.flags['C_CONTIGUOUS']==False: + arr = np.ascontiguousarray(arr) + return arr + + def _matvec(self,v): + """Implements the forward operator. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + v = self.__checkArray(v, self.vshape) + vid = self.data_mod.link('-vol',self.vg,v) + s = np.zeros(self.sshape,dtype=np.float32) + sid = self.data_mod.link('-sino',self.pg,s) + + cfg = creators.astra_dict('FP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['VolumeDataId'] = vid + cfg['ProjectorId'] = self.proj_id + fp_id = algorithm.create(cfg) + algorithm.run(fp_id) + + algorithm.delete(fp_id) + self.data_mod.delete([vid,sid]) + return s.flatten() + + def rmatvec(self,s): + """Implements the transpose operator. + + :param s: The projection data. + :type s: :class:`numpy.ndarray` + """ + s = self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + v = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,v) + + cfg = creators.astra_dict('BP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + bp_id = algorithm.create(cfg) + algorithm.run(bp_id) + + algorithm.delete(bp_id) + self.data_mod.delete([vid,sid]) + return v.flatten() + + def __mul__(self,v): + """Provides easy forward operator by *. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + # Catch the case of a forward projection of a 2D/3D image + if isinstance(v, np.ndarray) and v.shape==self.vshape: + return self._matvec(v) + return scipy.sparse.linalg.LinearOperator.__mul__(self, v) + + def reconstruct(self, method, s, iterations=1, extraOptions = {}): + """Reconstruct an object. + + :param method: Method to use for reconstruction. + :type method: :class:`string` + :param s: The projection data. + :type s: :class:`numpy.ndarray` + :param iterations: Number of iterations to use. + :type iterations: :class:`int` + :param extraOptions: Extra options to use during reconstruction (i.e. for cfg['option']). + :type extraOptions: :class:`dict` + """ + self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + v = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,v) + cfg = creators.astra_dict(method) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + cfg['option'] = extraOptions + alg_id = algorithm.create(cfg) + algorithm.run(alg_id,iterations) + algorithm.delete(alg_id) + self.data_mod.delete([vid,sid]) + return v + +class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): + """This object provides the transpose operation (``.T``) of the OpTomo object. + + Do not use directly, since it can be accessed as member ``.T`` of + an :class:`OpTomo` object. + """ + def __init__(self,parent): + self.parent = parent + self.dtype = np.float32 + self.shape = (parent.shape[1], parent.shape[0]) + + def _matvec(self, s): + return self.parent.rmatvec(s) + + def rmatvec(self, v): + return self.parent.matvec(v) + + def __mul__(self,s): + # Catch the case of a backprojection of 2D/3D data + if isinstance(s, np.ndarray) and s.shape==self.parent.sshape: + return self._matvec(s) + return scipy.sparse.linalg.LinearOperator.__mul__(self, s) -- cgit v1.2.3