From f69d9f6bc1704560518da3c30c46e495c0228aac Mon Sep 17 00:00:00 2001
From: Daan Pelt <daan.pelt@gmail.com>
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 = <CData2DManager * >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 = <CData3DManager * >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 <http://www.gnu.org/licenses/>.
+#
+#-----------------------------------------------------------------------
+"""Additional purely Python functions for PyAstraToolbox.
+
+.. moduleauthor:: Daniel M. Pelt <D.M.Pelt@cwi.nl>
+
+
+"""
+
+def geom_size(geom, dim=None):
+    """Returns the size of a volume or sinogram, based on the projection or volume geometry.
+
+    :param geom: Geometry to calculate size from
+    :type geometry: :class:`dict`
+    :param dim: Optional axis index to return
+    :type dim: :class:`int`
+    """
+
+    if 'GridSliceCount' in geom:
+        # 3D Volume geometry?
+        s = (geom['GridSliceCount'], geom[
+             'GridRowCount'], geom['GridColCount'])
+    elif 'GridColCount' in geom:
+        # 2D Volume geometry?
+        s = (geom['GridRowCount'], geom['GridColCount'])
+    elif geom['type'] == 'parallel' or geom['type'] == 'fanflat':
+        s = (len(geom['ProjectionAngles']), geom['DetectorCount'])
+    elif geom['type'] == 'parallel3d' or geom['type'] == 'cone':
+        s = (geom['DetectorRowCount'], len(
+            geom['ProjectionAngles']), geom['DetectorColCount'])
+    elif geom['type'] == 'fanflat_vec':
+        s = (geom['Vectors'].shape[0], geom['DetectorCount'])
+    elif geom['type'] == 'parallel3d_vec' or geom['type'] == 'cone_vec':
+        s = (geom['DetectorRowCount'], geom[
+             'Vectors'].shape[0], geom['DetectorColCount'])
+
+    if dim != None:
+        s = s[dim]
+
+    return s
-- 
cgit v1.2.3


From 2bc0d98c413fee4108115f26aa337f65337eec55 Mon Sep 17 00:00:00 2001
From: Daan Pelt <daan.pelt@gmail.com>
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 <http://www.gnu.org/licenses/>.
-#
-#-----------------------------------------------------------------------
-
-import math
-from . import creators as ac
-from . import data2d
-
-
-class ASTRAProjector2DTranspose():
-    """Implements the ``proj.T`` functionality.
-
-    Do not use directly, since it can be accessed as member ``.T`` of
-    an :class:`ASTRAProjector2D` object.
-
-    """
-    def __init__(self, parentProj):
-        self.parentProj = parentProj
-
-    def __mul__(self, data):
-        return self.parentProj.backProject(data)
-
-
-class ASTRAProjector2D(object):
-    """Helps with various common ASTRA Toolbox 2D operations.
-
-    This class can perform several often used toolbox operations, such as:
-
-    * Forward projecting
-    * Back projecting
-    * Reconstructing
-
-    Note that this class has a some computational overhead, because it
-    copies a lot of data. If you use many repeated operations, directly
-    using the PyAstraToolbox methods directly is faster.
-
-    You can use this class as an abstracted weight matrix :math:`W`: multiplying an instance
-    ``proj`` of this class by an image results in a forward projection of the image, and multiplying
-    ``proj.T`` by a sinogram results in a backprojection of the sinogram::
-
-        proj = ASTRAProjector2D(...)
-        fp = proj*image
-        bp = proj.T*sinogram
-
-    :param proj_geom: The projection geometry.
-    :type proj_geom: :class:`dict`
-    :param vol_geom: The volume geometry.
-    :type vol_geom: :class:`dict`
-    :param proj_type: Projector type, such as ``'line'``, ``'linear'``, ...
-    :type proj_type: :class:`string`
-    """
-
-    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 <http://www.gnu.org/licenses/>.
+#
+#-----------------------------------------------------------------------
+
+from . import data2d
+from . import data3d
+from . import projector
+from . import projector3d
+from . import creators
+from . import algorithm
+from . import functions
+import numpy as np
+from six.moves import 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 <wjp@usecode.org>
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'), <string> wrap_to_bytes(val))
             else:
                 itm = root.addChildNode(item, wrap_to_bytes(val))
-                del itm
 
-cdef void readOptions(XMLNode * node, dc):
-    cdef XMLNode * listbase
-    cdef XMLNode * itm
+cdef void readOptions(XMLNode node, dc):
+    cdef XMLNode listbase
+    cdef XMLNode itm
     cdef int i
     cdef int j
     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 <Willem.Jan.Palenstijn@cwi.nl>
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 = <double*>np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) 
             if val.ndim == 2:
-                for i in range(val.shape[0]):
-                    for j in range(val.shape[1]):
-                        itm = listbase.addChildNode(six.b('ListItem'))
-                        itm.addAttribute(< string > six.b('index'), < float32 > index)
-                        itm.addAttribute( < string > six.b('value'), < float32 > val[i, j])
-                        index += 1
+                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 = <double*>np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) 
             if val.ndim == 2:
-                for i in range(val.shape[0]):
-                    for j in range(val.shape[1]):
-                        itm = listbase.addChildNode(six.b('ListItem'))
-                        itm.addAttribute(< string > six.b('index'), < float32 > index)
-                        itm.addAttribute( < string > six.b('value'), < float32 > val[i, j])
-                        index += 1
+                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" <D.M.Pelt@cwi.nl>
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" <D.M.Pelt@cwi.nl>
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" <D.M.Pelt@cwi.nl>
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 <http://www.gnu.org/licenses/>.
-#
-#-----------------------------------------------------------------------
-
-from . import data2d
-from . import data3d
-from . import projector
-from . import projector3d
-from . import creators
-from . import algorithm
-from . import functions
-import numpy as np
-from six.moves import 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 <http://www.gnu.org/licenses/>.
+#
+#-----------------------------------------------------------------------
+
+from . import data2d
+from . import data3d
+from . import projector
+from . import projector3d
+from . import creators
+from . import algorithm
+from . import functions
+import numpy as np
+from six.moves import 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 <daan.pelt@gmail.com>
Date: Mon, 18 May 2015 20:32:52 +0200
Subject: Hold reference to results of six.b calls

---
 python/astra/log_c.pyx | 21 ++++++++++++++-------
 1 file changed, 14 insertions(+), 7 deletions(-)

(limited to 'python/astra')

diff --git a/python/astra/log_c.pyx b/python/astra/log_c.pyx
index 969cc06..f16329f 100644
--- a/python/astra/log_c.pyx
+++ b/python/astra/log_c.pyx
@@ -52,16 +52,20 @@ cdef extern from "astra/Logging.h" namespace "astra::CLogger":
     void setFormatScreen(const char *fmt)
 
 def log_debug(sfile, sline, message):
-    debug(six.b(sfile),sline,six.b(message))
+    cstr = list(map(six.b,(sfile,message)))
+    debug(cstr[0],sline,cstr[1])
 
 def log_info(sfile, sline, message):
-    info(six.b(sfile),sline,six.b(message))
+    cstr = list(map(six.b,(sfile,message)))
+    info(cstr[0],sline,cstr[1])
 
 def log_warn(sfile, sline, message):
-    warn(six.b(sfile),sline,six.b(message))
+    cstr = list(map(six.b,(sfile,message)))
+    warn(cstr[0],sline,cstr[1])
 
 def log_error(sfile, sline, message):
-    error(six.b(sfile),sline,six.b(message))
+    cstr = list(map(six.b,(sfile,message)))
+    error(cstr[0],sline,cstr[1])
 
 def log_enable():
     enable()
@@ -82,10 +86,12 @@ def log_disableFile():
     disableFile()
 
 def log_setFormatFile(fmt):
-    setFormatFile(six.b(fmt))
+    cstr = six.b(fmt)
+    setFormatFile(cstr)
 
 def log_setFormatScreen(fmt):
-    setFormatScreen(six.b(fmt))
+    cstr = six.b(fmt)
+    setFormatScreen(cstr)
 
 enumList = [LOG_DEBUG,LOG_INFO,LOG_WARN,LOG_ERROR]
 
@@ -93,4 +99,5 @@ def log_setOutputScreen(fd, level):
     setOutputScreen(fd, enumList[level])
 
 def log_setOutputFile(filename, level):
-    setOutputFile(six.b(filename), enumList[level])
\ No newline at end of file
+    cstr = six.b(filename)
+    setOutputFile(cstr, enumList[level])
-- 
cgit v1.2.3


From aed157812c31df681a7ba7dd4479cb35fbce05bb Mon Sep 17 00:00:00 2001
From: "Daniel M. Pelt" <D.M.Pelt@cwi.nl>
Date: Mon, 18 May 2015 20:54:54 +0200
Subject: Fix range import for old six versions

---
 python/astra/functions.py | 6 +++++-
 python/astra/matrix_c.pyx | 7 ++++++-
 python/astra/optomo.py    | 8 +++++++-
 3 files changed, 18 insertions(+), 3 deletions(-)

(limited to 'python/astra')

diff --git a/python/astra/functions.py b/python/astra/functions.py
index b826b86..e38b5bc 100644
--- a/python/astra/functions.py
+++ b/python/astra/functions.py
@@ -32,7 +32,11 @@
 
 from . import creators as ac
 import numpy as np
-from six.moves import range
+try:
+    from six.moves import range
+except ImportError:
+    # six 1.3.0
+    from six.moves import xrange as range
 
 from . import data2d
 from . import data3d
diff --git a/python/astra/matrix_c.pyx b/python/astra/matrix_c.pyx
index b0d8bc4..d099a75 100644
--- a/python/astra/matrix_c.pyx
+++ b/python/astra/matrix_c.pyx
@@ -27,7 +27,12 @@
 # distutils: libraries = astra
 
 import six
-from six.moves import range
+try:
+    from six.moves import range
+except ImportError:
+    # six 1.3.0
+    from six.moves import xrange as range
+
 import numpy as np
 import scipy.sparse as ss
 
diff --git a/python/astra/optomo.py b/python/astra/optomo.py
index 0c37353..2937d9c 100644
--- a/python/astra/optomo.py
+++ b/python/astra/optomo.py
@@ -32,7 +32,13 @@ from . import creators
 from . import algorithm
 from . import functions
 import numpy as np
-from six.moves import range, reduce
+from six.moves import reduce
+try:
+    from six.moves import range
+except ImportError:
+    # six 1.3.0
+    from six.moves import xrange as range
+
 import operator
 import scipy.sparse.linalg
 
-- 
cgit v1.2.3


From 0f577d1fbc2b0c15d85f18cc38eb14e3cbf6c6a2 Mon Sep 17 00:00:00 2001
From: "Daniel M. Pelt" <D.M.Pelt@cwi.nl>
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/astra')

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 26198f7e1fbed9d294fc2856f796f408c85b66f3 Mon Sep 17 00:00:00 2001
From: Valerii Sokolov <valerii.sokolov@uantwerpen.be>
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/astra')

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 = <double*>np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) 
+            contig_data = np.ascontiguousarray(val,dtype=np.float64)
+            data = <double*>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 = <double*>np.PyArray_DATA(np.ascontiguousarray(val,dtype=np.float64)) 
+            contig_data = np.ascontiguousarray(val,dtype=np.float64)
+            data = <double*>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 63d78fbaafa7d247347f9052db86f575d89260b7 Mon Sep 17 00:00:00 2001
From: "Daniel M. Pelt" <D.M.Pelt@cwi.nl>
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/astra')

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" <D.M.Pelt@cwi.nl>
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/astra')

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" <D.M.Pelt@cwi.nl>
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/astra')

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",<char*>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",<char*>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",<char*>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",<char*>cstr[1])
 
 def log_enable():
     enable()
-- 
cgit v1.2.3


From 43155c03488ca98c24f8f369e5c8699efa20cca3 Mon Sep 17 00:00:00 2001
From: "Daniel M. Pelt" <D.M.Pelt@cwi.nl>
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/astra')

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" <D.M.Pelt@cwi.nl>
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/astra')

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" <D.M.Pelt@cwi.nl>
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/astra')

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