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') 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 +++++++++++++++++++++++++++++++++++++++ python/docSRC/ASTRAProjector.rst | 8 -- python/docSRC/index.rst | 2 +- python/docSRC/operator.rst | 8 ++ 6 files changed, 218 insertions(+), 145 deletions(-) delete mode 100644 python/astra/ASTRAProjector.py create mode 100644 python/astra/operator.py delete mode 100644 python/docSRC/ASTRAProjector.rst create mode 100644 python/docSRC/operator.rst (limited to 'python') 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) diff --git a/python/docSRC/ASTRAProjector.rst b/python/docSRC/ASTRAProjector.rst deleted file mode 100644 index 1c267e3..0000000 --- a/python/docSRC/ASTRAProjector.rst +++ /dev/null @@ -1,8 +0,0 @@ -Helper class: the :mod:`ASTRAProjector` module -============================================== - -.. automodule:: astra.ASTRAProjector - :members: - :undoc-members: - :show-inheritance: - diff --git a/python/docSRC/index.rst b/python/docSRC/index.rst index 8d17a4a..b7cc6d6 100644 --- a/python/docSRC/index.rst +++ b/python/docSRC/index.rst @@ -18,7 +18,7 @@ Contents: matrix creators functions - ASTRAProjector + operator matlab astra .. astra diff --git a/python/docSRC/operator.rst b/python/docSRC/operator.rst new file mode 100644 index 0000000..f5369fa --- /dev/null +++ b/python/docSRC/operator.rst @@ -0,0 +1,8 @@ +OpTomo class: the :mod:`operator` module +============================================== + +.. automodule:: astra.operator + :members: + :undoc-members: + :show-inheritance: + -- 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') 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') 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') 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') 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') diff --git a/python/astra/__init__.py b/python/astra/__init__.py index 8c1740c..6c15d30 100644 --- a/python/astra/__init__.py +++ b/python/astra/__init__.py @@ -35,7 +35,7 @@ from . import projector from . import projector3d from . import matrix from . import log -from .operator import OpTomo +from .optomo import OpTomo import os try: diff --git a/python/astra/operator.py b/python/astra/operator.py deleted file mode 100644 index 0c37353..0000000 --- a/python/astra/operator.py +++ /dev/null @@ -1,197 +0,0 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam -# -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ -# -# -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). -# -#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify -#it under the terms of the GNU General Public License as published by -#the Free Software Foundation, either version 3 of the License, or -#(at your option) any later version. -# -#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, -#but WITHOUT ANY WARRANTY; without even the implied warranty of -#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -#GNU General Public License for more details. -# -#You should have received a copy of the GNU General Public License -#along with the Python interface to the ASTRA Toolbox. If not, see . -# -#----------------------------------------------------------------------- - -from . import data2d -from . import data3d -from . import projector -from . import projector3d -from . import creators -from . import algorithm -from . import functions -import numpy as np -from six.moves import range, reduce -import operator -import scipy.sparse.linalg - -class OpTomo(scipy.sparse.linalg.LinearOperator): - """Object that imitates a projection matrix with a given projector. - - This object can do forward projection by using the ``*`` operator:: - - W = astra.OpTomo(proj_id) - fp = W*image - bp = W.T*sinogram - - It can also be used in minimization methods of the :mod:`scipy.sparse.linalg` module:: - - W = astra.OpTomo(proj_id) - output = scipy.sparse.linalg.lsqr(W,sinogram) - - :param proj_id: ID to a projector. - :type proj_id: :class:`int` - """ - - def __init__(self,proj_id): - self.dtype = np.float32 - try: - self.vg = projector.volume_geometry(proj_id) - self.pg = projector.projection_geometry(proj_id) - self.data_mod = data2d - self.appendString = "" - if projector.is_cuda(proj_id): - self.appendString += "_CUDA" - except Exception: - self.vg = projector3d.volume_geometry(proj_id) - self.pg = projector3d.projection_geometry(proj_id) - self.data_mod = data3d - self.appendString = "3D" - if projector3d.is_cuda(proj_id): - self.appendString += "_CUDA" - - self.vshape = functions.geom_size(self.vg) - self.vsize = reduce(operator.mul,self.vshape) - self.sshape = functions.geom_size(self.pg) - self.ssize = reduce(operator.mul,self.sshape) - - self.shape = (self.ssize, self.vsize) - - self.proj_id = proj_id - - self.T = OpTomoTranspose(self) - - def __checkArray(self, arr, shp): - if len(arr.shape)==1: - arr = arr.reshape(shp) - if arr.dtype != np.float32: - arr = arr.astype(np.float32) - if arr.flags['C_CONTIGUOUS']==False: - arr = np.ascontiguousarray(arr) - return arr - - def _matvec(self,v): - """Implements the forward operator. - - :param v: Volume to forward project. - :type v: :class:`numpy.ndarray` - """ - v = self.__checkArray(v, self.vshape) - vid = self.data_mod.link('-vol',self.vg,v) - s = np.zeros(self.sshape,dtype=np.float32) - sid = self.data_mod.link('-sino',self.pg,s) - - cfg = creators.astra_dict('FP'+self.appendString) - cfg['ProjectionDataId'] = sid - cfg['VolumeDataId'] = vid - cfg['ProjectorId'] = self.proj_id - fp_id = algorithm.create(cfg) - algorithm.run(fp_id) - - algorithm.delete(fp_id) - self.data_mod.delete([vid,sid]) - return s.flatten() - - def rmatvec(self,s): - """Implements the transpose operator. - - :param s: The projection data. - :type s: :class:`numpy.ndarray` - """ - s = self.__checkArray(s, self.sshape) - sid = self.data_mod.link('-sino',self.pg,s) - v = np.zeros(self.vshape,dtype=np.float32) - vid = self.data_mod.link('-vol',self.vg,v) - - cfg = creators.astra_dict('BP'+self.appendString) - cfg['ProjectionDataId'] = sid - cfg['ReconstructionDataId'] = vid - cfg['ProjectorId'] = self.proj_id - bp_id = algorithm.create(cfg) - algorithm.run(bp_id) - - algorithm.delete(bp_id) - self.data_mod.delete([vid,sid]) - return v.flatten() - - def __mul__(self,v): - """Provides easy forward operator by *. - - :param v: Volume to forward project. - :type v: :class:`numpy.ndarray` - """ - # Catch the case of a forward projection of a 2D/3D image - if isinstance(v, np.ndarray) and v.shape==self.vshape: - return self._matvec(v) - return scipy.sparse.linalg.LinearOperator.__mul__(self, v) - - def reconstruct(self, method, s, iterations=1, extraOptions = {}): - """Reconstruct an object. - - :param method: Method to use for reconstruction. - :type method: :class:`string` - :param s: The projection data. - :type s: :class:`numpy.ndarray` - :param iterations: Number of iterations to use. - :type iterations: :class:`int` - :param extraOptions: Extra options to use during reconstruction (i.e. for cfg['option']). - :type extraOptions: :class:`dict` - """ - self.__checkArray(s, self.sshape) - sid = self.data_mod.link('-sino',self.pg,s) - v = np.zeros(self.vshape,dtype=np.float32) - vid = self.data_mod.link('-vol',self.vg,v) - cfg = creators.astra_dict(method) - cfg['ProjectionDataId'] = sid - cfg['ReconstructionDataId'] = vid - cfg['ProjectorId'] = self.proj_id - cfg['option'] = extraOptions - alg_id = algorithm.create(cfg) - algorithm.run(alg_id,iterations) - algorithm.delete(alg_id) - self.data_mod.delete([vid,sid]) - return v - -class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): - """This object provides the transpose operation (``.T``) of the OpTomo object. - - Do not use directly, since it can be accessed as member ``.T`` of - an :class:`OpTomo` object. - """ - def __init__(self,parent): - self.parent = parent - self.dtype = np.float32 - self.shape = (parent.shape[1], parent.shape[0]) - - def _matvec(self, s): - return self.parent.rmatvec(s) - - def rmatvec(self, v): - return self.parent.matvec(v) - - def __mul__(self,s): - # Catch the case of a backprojection of 2D/3D data - if isinstance(s, np.ndarray) and s.shape==self.parent.sshape: - return self._matvec(s) - return scipy.sparse.linalg.LinearOperator.__mul__(self, s) diff --git a/python/astra/optomo.py b/python/astra/optomo.py new file mode 100644 index 0000000..0c37353 --- /dev/null +++ b/python/astra/optomo.py @@ -0,0 +1,197 @@ +#----------------------------------------------------------------------- +#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# +#Author: Daniel M. Pelt +#Contact: D.M.Pelt@cwi.nl +#Website: http://dmpelt.github.io/pyastratoolbox/ +# +# +#This file is part of the Python interface to the +#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). +# +#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with the Python interface to the ASTRA Toolbox. If not, see . +# +#----------------------------------------------------------------------- + +from . import data2d +from . import data3d +from . import projector +from . import projector3d +from . import creators +from . import algorithm +from . import functions +import numpy as np +from six.moves import range, reduce +import operator +import scipy.sparse.linalg + +class OpTomo(scipy.sparse.linalg.LinearOperator): + """Object that imitates a projection matrix with a given projector. + + This object can do forward projection by using the ``*`` operator:: + + W = astra.OpTomo(proj_id) + fp = W*image + bp = W.T*sinogram + + It can also be used in minimization methods of the :mod:`scipy.sparse.linalg` module:: + + W = astra.OpTomo(proj_id) + output = scipy.sparse.linalg.lsqr(W,sinogram) + + :param proj_id: ID to a projector. + :type proj_id: :class:`int` + """ + + def __init__(self,proj_id): + self.dtype = np.float32 + try: + self.vg = projector.volume_geometry(proj_id) + self.pg = projector.projection_geometry(proj_id) + self.data_mod = data2d + self.appendString = "" + if projector.is_cuda(proj_id): + self.appendString += "_CUDA" + except Exception: + self.vg = projector3d.volume_geometry(proj_id) + self.pg = projector3d.projection_geometry(proj_id) + self.data_mod = data3d + self.appendString = "3D" + if projector3d.is_cuda(proj_id): + self.appendString += "_CUDA" + + self.vshape = functions.geom_size(self.vg) + self.vsize = reduce(operator.mul,self.vshape) + self.sshape = functions.geom_size(self.pg) + self.ssize = reduce(operator.mul,self.sshape) + + self.shape = (self.ssize, self.vsize) + + self.proj_id = proj_id + + self.T = OpTomoTranspose(self) + + def __checkArray(self, arr, shp): + if len(arr.shape)==1: + arr = arr.reshape(shp) + if arr.dtype != np.float32: + arr = arr.astype(np.float32) + if arr.flags['C_CONTIGUOUS']==False: + arr = np.ascontiguousarray(arr) + return arr + + def _matvec(self,v): + """Implements the forward operator. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + v = self.__checkArray(v, self.vshape) + vid = self.data_mod.link('-vol',self.vg,v) + s = np.zeros(self.sshape,dtype=np.float32) + sid = self.data_mod.link('-sino',self.pg,s) + + cfg = creators.astra_dict('FP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['VolumeDataId'] = vid + cfg['ProjectorId'] = self.proj_id + fp_id = algorithm.create(cfg) + algorithm.run(fp_id) + + algorithm.delete(fp_id) + self.data_mod.delete([vid,sid]) + return s.flatten() + + def rmatvec(self,s): + """Implements the transpose operator. + + :param s: The projection data. + :type s: :class:`numpy.ndarray` + """ + s = self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + v = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,v) + + cfg = creators.astra_dict('BP'+self.appendString) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + bp_id = algorithm.create(cfg) + algorithm.run(bp_id) + + algorithm.delete(bp_id) + self.data_mod.delete([vid,sid]) + return v.flatten() + + def __mul__(self,v): + """Provides easy forward operator by *. + + :param v: Volume to forward project. + :type v: :class:`numpy.ndarray` + """ + # Catch the case of a forward projection of a 2D/3D image + if isinstance(v, np.ndarray) and v.shape==self.vshape: + return self._matvec(v) + return scipy.sparse.linalg.LinearOperator.__mul__(self, v) + + def reconstruct(self, method, s, iterations=1, extraOptions = {}): + """Reconstruct an object. + + :param method: Method to use for reconstruction. + :type method: :class:`string` + :param s: The projection data. + :type s: :class:`numpy.ndarray` + :param iterations: Number of iterations to use. + :type iterations: :class:`int` + :param extraOptions: Extra options to use during reconstruction (i.e. for cfg['option']). + :type extraOptions: :class:`dict` + """ + self.__checkArray(s, self.sshape) + sid = self.data_mod.link('-sino',self.pg,s) + v = np.zeros(self.vshape,dtype=np.float32) + vid = self.data_mod.link('-vol',self.vg,v) + cfg = creators.astra_dict(method) + cfg['ProjectionDataId'] = sid + cfg['ReconstructionDataId'] = vid + cfg['ProjectorId'] = self.proj_id + cfg['option'] = extraOptions + alg_id = algorithm.create(cfg) + algorithm.run(alg_id,iterations) + algorithm.delete(alg_id) + self.data_mod.delete([vid,sid]) + return v + +class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): + """This object provides the transpose operation (``.T``) of the OpTomo object. + + Do not use directly, since it can be accessed as member ``.T`` of + an :class:`OpTomo` object. + """ + def __init__(self,parent): + self.parent = parent + self.dtype = np.float32 + self.shape = (parent.shape[1], parent.shape[0]) + + def _matvec(self, s): + return self.parent.rmatvec(s) + + def rmatvec(self, v): + return self.parent.matvec(v) + + def __mul__(self,s): + # Catch the case of a backprojection of 2D/3D data + if isinstance(s, np.ndarray) and s.shape==self.parent.sshape: + return self._matvec(s) + return scipy.sparse.linalg.LinearOperator.__mul__(self, s) -- cgit v1.2.3 From 0f4969c3011910518639e3dd9ff6c1e4b2539f4c Mon Sep 17 00:00:00 2001 From: Daan Pelt Date: Mon, 18 May 2015 20:32:52 +0200 Subject: Hold reference to results of six.b calls --- python/astra/log_c.pyx | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) (limited to 'python') diff --git a/python/astra/log_c.pyx b/python/astra/log_c.pyx index 969cc06..f16329f 100644 --- a/python/astra/log_c.pyx +++ b/python/astra/log_c.pyx @@ -52,16 +52,20 @@ cdef extern from "astra/Logging.h" namespace "astra::CLogger": void setFormatScreen(const char *fmt) def log_debug(sfile, sline, message): - debug(six.b(sfile),sline,six.b(message)) + cstr = list(map(six.b,(sfile,message))) + debug(cstr[0],sline,cstr[1]) def log_info(sfile, sline, message): - info(six.b(sfile),sline,six.b(message)) + cstr = list(map(six.b,(sfile,message))) + info(cstr[0],sline,cstr[1]) def log_warn(sfile, sline, message): - warn(six.b(sfile),sline,six.b(message)) + cstr = list(map(six.b,(sfile,message))) + warn(cstr[0],sline,cstr[1]) def log_error(sfile, sline, message): - error(six.b(sfile),sline,six.b(message)) + cstr = list(map(six.b,(sfile,message))) + error(cstr[0],sline,cstr[1]) def log_enable(): enable() @@ -82,10 +86,12 @@ def log_disableFile(): disableFile() def log_setFormatFile(fmt): - setFormatFile(six.b(fmt)) + cstr = six.b(fmt) + setFormatFile(cstr) def log_setFormatScreen(fmt): - setFormatScreen(six.b(fmt)) + cstr = six.b(fmt) + setFormatScreen(cstr) enumList = [LOG_DEBUG,LOG_INFO,LOG_WARN,LOG_ERROR] @@ -93,4 +99,5 @@ def log_setOutputScreen(fd, level): setOutputScreen(fd, enumList[level]) def log_setOutputFile(filename, level): - setOutputFile(six.b(filename), enumList[level]) \ No newline at end of file + cstr = six.b(filename) + setOutputFile(cstr, enumList[level]) -- cgit v1.2.3 From aed157812c31df681a7ba7dd4479cb35fbce05bb Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 18 May 2015 20:54:54 +0200 Subject: Fix range import for old six versions --- python/astra/functions.py | 6 +++++- python/astra/matrix_c.pyx | 7 ++++++- python/astra/optomo.py | 8 +++++++- 3 files changed, 18 insertions(+), 3 deletions(-) (limited to 'python') diff --git a/python/astra/functions.py b/python/astra/functions.py index b826b86..e38b5bc 100644 --- a/python/astra/functions.py +++ b/python/astra/functions.py @@ -32,7 +32,11 @@ from . import creators as ac import numpy as np -from six.moves import range +try: + from six.moves import range +except ImportError: + # six 1.3.0 + from six.moves import xrange as range from . import data2d from . import data3d diff --git a/python/astra/matrix_c.pyx b/python/astra/matrix_c.pyx index b0d8bc4..d099a75 100644 --- a/python/astra/matrix_c.pyx +++ b/python/astra/matrix_c.pyx @@ -27,7 +27,12 @@ # distutils: libraries = astra import six -from six.moves import range +try: + from six.moves import range +except ImportError: + # six 1.3.0 + from six.moves import xrange as range + import numpy as np import scipy.sparse as ss diff --git a/python/astra/optomo.py b/python/astra/optomo.py index 0c37353..2937d9c 100644 --- a/python/astra/optomo.py +++ b/python/astra/optomo.py @@ -32,7 +32,13 @@ from . import creators from . import algorithm from . import functions import numpy as np -from six.moves import range, reduce +from six.moves import reduce +try: + from six.moves import range +except ImportError: + # six 1.3.0 + from six.moves import xrange as range + import operator import scipy.sparse.linalg -- cgit v1.2.3 From 0f577d1fbc2b0c15d85f18cc38eb14e3cbf6c6a2 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Fri, 22 May 2015 14:36:42 +0200 Subject: Fix optomo reconstruction for 1D input --- python/astra/optomo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'python') diff --git a/python/astra/optomo.py b/python/astra/optomo.py index 2937d9c..0108674 100644 --- a/python/astra/optomo.py +++ b/python/astra/optomo.py @@ -164,7 +164,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): :param extraOptions: Extra options to use during reconstruction (i.e. for cfg['option']). :type extraOptions: :class:`dict` """ - self.__checkArray(s, self.sshape) + 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) -- cgit v1.2.3 From ef47dca509c1c26c5037b9b9024a6a84a3b7fe0b Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Sat, 20 Jun 2015 00:31:15 +0200 Subject: Fix wrong module name in documentation --- python/docSRC/operator.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'python') diff --git a/python/docSRC/operator.rst b/python/docSRC/operator.rst index f5369fa..65a5a99 100644 --- a/python/docSRC/operator.rst +++ b/python/docSRC/operator.rst @@ -1,7 +1,7 @@ OpTomo class: the :mod:`operator` module ============================================== -.. automodule:: astra.operator +.. automodule:: astra.optomo :members: :undoc-members: :show-inheritance: -- cgit v1.2.3 From 26198f7e1fbed9d294fc2856f796f408c85b66f3 Mon Sep 17 00:00:00 2001 From: Valerii Sokolov Date: Tue, 9 Jun 2015 11:26:42 +0200 Subject: Don't take address of temporary. --- python/astra/utils.pyx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) (limited to 'python') diff --git a/python/astra/utils.pyx b/python/astra/utils.pyx index ddb37aa..a8e9e4e 100644 --- a/python/astra/utils.pyx +++ b/python/astra/utils.pyx @@ -95,7 +95,8 @@ cdef void readDict(XMLNode root, _dc): if val.size == 0: break listbase = root.addChildNode(item) - data = np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) + contig_data = np.ascontiguousarray(val,dtype=np.float64) + data = np.PyArray_DATA(contig_data) if val.ndim == 2: listbase.setContent(data, val.shape[1], val.shape[0], False) elif val.ndim == 1: @@ -129,7 +130,8 @@ cdef void readOptions(XMLNode node, dc): break listbase = node.addChildNode(six.b('Option')) listbase.addAttribute(< string > six.b('key'), < string > item) - data = np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) + contig_data = np.ascontiguousarray(val,dtype=np.float64) + data = np.PyArray_DATA(contig_data) if val.ndim == 2: listbase.setContent(data, val.shape[1], val.shape[0], False) elif val.ndim == 1: -- cgit v1.2.3 From a9ea08c4a11592378b320cac45be8eede8addd6f Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 22 Jun 2015 17:54:23 +0200 Subject: Additional module name fix in documentation --- python/docSRC/operator.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'python') diff --git a/python/docSRC/operator.rst b/python/docSRC/operator.rst index 65a5a99..fe500ba 100644 --- a/python/docSRC/operator.rst +++ b/python/docSRC/operator.rst @@ -1,4 +1,4 @@ -OpTomo class: the :mod:`operator` module +OpTomo class: the :mod:`optomo` module ============================================== .. automodule:: astra.optomo -- cgit v1.2.3 From 63d78fbaafa7d247347f9052db86f575d89260b7 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Wed, 24 Jun 2015 20:28:46 +0200 Subject: Fix config to struct/dict translation for array options --- python/astra/utils.pyx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'python') diff --git a/python/astra/utils.pyx b/python/astra/utils.pyx index a8e9e4e..260c308 100644 --- a/python/astra/utils.pyx +++ b/python/astra/utils.pyx @@ -204,7 +204,10 @@ cdef XMLNode2dict(XMLNode node): while it != nodes.end(): subnode = deref(it) if castString(subnode.getName())=="Option": - opts[castString(subnode.getAttribute('key'))] = stringToPythonValue(subnode.getAttribute('value')) + if subnode.hasAttribute('value'): + opts[castString(subnode.getAttribute('key'))] = stringToPythonValue(subnode.getAttribute('value')) + else: + opts[castString(subnode.getAttribute('key'))] = stringToPythonValue(subnode.getContent()) else: dct[castString(subnode.getName())] = stringToPythonValue(subnode.getContent()) inc(it) -- cgit v1.2.3 From 9e3472ea9041b8755050427d8bdb8a4701019c55 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 25 Jun 2015 21:52:07 +0200 Subject: Fix memory leak in configuration code --- python/astra/projector3d_c.pyx | 10 ++++++++-- python/astra/projector_c.pyx | 10 ++++++++-- 2 files changed, 16 insertions(+), 4 deletions(-) (limited to 'python') diff --git a/python/astra/projector3d_c.pyx b/python/astra/projector3d_c.pyx index 8b978d7..aec9cde 100644 --- a/python/astra/projector3d_c.pyx +++ b/python/astra/projector3d_c.pyx @@ -87,12 +87,18 @@ cdef CProjector3D * getObject(i) except NULL: def projection_geometry(i): cdef CProjector3D * proj = getObject(i) - return utils.configToDict(proj.getProjectionGeometry().getConfiguration()) + cdef Config * cfg = proj.getProjectionGeometry().getConfiguration() + dct = utils.configToDict(cfg) + del cfg + return dct def volume_geometry(i): cdef CProjector3D * proj = getObject(i) - return utils.configToDict(proj.getVolumeGeometry().getConfiguration()) + cdef Config * cfg = proj.getVolumeGeometry().getConfiguration() + dct = utils.configToDict(cfg) + del cfg + return dct def weights_single_ray(i, projection_index, detector_index): diff --git a/python/astra/projector_c.pyx b/python/astra/projector_c.pyx index 9aa868e..77c64a4 100644 --- a/python/astra/projector_c.pyx +++ b/python/astra/projector_c.pyx @@ -91,12 +91,18 @@ cdef CProjector2D * getObject(i) except NULL: def projection_geometry(i): cdef CProjector2D * proj = getObject(i) - return utils.configToDict(proj.getProjectionGeometry().getConfiguration()) + cdef Config * cfg = proj.getProjectionGeometry().getConfiguration() + dct = utils.configToDict(cfg) + del cfg + return dct def volume_geometry(i): cdef CProjector2D * proj = getObject(i) - return utils.configToDict(proj.getVolumeGeometry().getConfiguration()) + cdef Config * cfg = proj.getVolumeGeometry().getConfiguration() + dct = utils.configToDict(cfg) + del cfg + return dct def weights_single_ray(i, projection_index, detector_index): -- cgit v1.2.3 From 4d39c35d6c9124c26de64c9d227a25f612903a2a Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Fri, 17 Jul 2015 13:44:59 +0200 Subject: Fix formatting when passing strings to log from high-level code --- python/astra/log_c.pyx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'python') diff --git a/python/astra/log_c.pyx b/python/astra/log_c.pyx index f16329f..55c63e6 100644 --- a/python/astra/log_c.pyx +++ b/python/astra/log_c.pyx @@ -53,19 +53,19 @@ cdef extern from "astra/Logging.h" namespace "astra::CLogger": def log_debug(sfile, sline, message): cstr = list(map(six.b,(sfile,message))) - debug(cstr[0],sline,cstr[1]) + debug(cstr[0],sline,"%s",cstr[1]) def log_info(sfile, sline, message): cstr = list(map(six.b,(sfile,message))) - info(cstr[0],sline,cstr[1]) + info(cstr[0],sline,"%s",cstr[1]) def log_warn(sfile, sline, message): cstr = list(map(six.b,(sfile,message))) - warn(cstr[0],sline,cstr[1]) + warn(cstr[0],sline,"%s",cstr[1]) def log_error(sfile, sline, message): cstr = list(map(six.b,(sfile,message))) - error(cstr[0],sline,cstr[1]) + error(cstr[0],sline,"%s",cstr[1]) def log_enable(): enable() -- cgit v1.2.3 From 43155c03488ca98c24f8f369e5c8699efa20cca3 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 30 Jul 2015 15:28:16 +0200 Subject: Fix Python OpTomo for scipy 0.16 scipy 0.16 also uses .T to define a transpose, which conflicts with the old OpTomo implementation. OpTomo now also defines the _transpose() method, which .T will call in scipy 0.16. --- python/astra/optomo.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) (limited to 'python') diff --git a/python/astra/optomo.py b/python/astra/optomo.py index 0108674..19b07e3 100644 --- a/python/astra/optomo.py +++ b/python/astra/optomo.py @@ -86,7 +86,15 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): self.proj_id = proj_id - self.T = OpTomoTranspose(self) + self.transposeOpTomo = OpTomoTranspose(self) + try: + self.T = self.transposeOpTomo + except AttributeError: + # Scipy >= 0.16 defines self.T using self._transpose() + pass + + def _transpose(self): + return self.transposeOpTomo def __checkArray(self, arr, shp): if len(arr.shape)==1: -- cgit v1.2.3 From b4bd441549ea71dd6c34a9f2158bbebc39ba4e9e Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Thu, 30 Jul 2015 15:34:10 +0200 Subject: Define a transpose for the OpTomo transpose as well Allows for chaining .T calls. --- python/astra/optomo.py | 8 ++++++++ 1 file changed, 8 insertions(+) (limited to 'python') diff --git a/python/astra/optomo.py b/python/astra/optomo.py index 19b07e3..4a64150 100644 --- a/python/astra/optomo.py +++ b/python/astra/optomo.py @@ -197,6 +197,11 @@ class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): self.parent = parent self.dtype = np.float32 self.shape = (parent.shape[1], parent.shape[0]) + try: + self.T = self.parent + except AttributeError: + # Scipy >= 0.16 defines self.T using self._transpose() + pass def _matvec(self, s): return self.parent.rmatvec(s) @@ -204,6 +209,9 @@ class OpTomoTranspose(scipy.sparse.linalg.LinearOperator): def rmatvec(self, v): return self.parent.matvec(v) + def _transpose(self): + return self.parent + 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: -- cgit v1.2.3 From ee790c305942675e94ee66bfd24896d1ef61335a Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 10 Aug 2015 16:22:19 +0200 Subject: Release the gil in algorithm.run --- python/astra/PyIncludes.pxd | 2 +- python/astra/algorithm_c.pyx | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) (limited to 'python') diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 909f58f..35dea5f 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -143,7 +143,7 @@ cdef extern from "astra/Float32ProjectionData2D.h" namespace "astra": cdef extern from "astra/Algorithm.h" namespace "astra": cdef cppclass CAlgorithm: bool initialize(Config) - void run(int) + void run(int) nogil bool isInitialized() cdef extern from "astra/ReconstructionAlgorithm2D.h" namespace "astra": diff --git a/python/astra/algorithm_c.pyx b/python/astra/algorithm_c.pyx index 966d3d7..3231c1f 100644 --- a/python/astra/algorithm_c.pyx +++ b/python/astra/algorithm_c.pyx @@ -73,7 +73,9 @@ cdef CAlgorithm * getAlg(i) except NULL: def run(i, iterations=0): cdef CAlgorithm * alg = getAlg(i) - alg.run(iterations) + cdef int its = iterations + with nogil: + alg.run(its) def get_res_norm(i): -- cgit v1.2.3