initial revision

14 files changed, 1103 insertions, 0 deletions

diff --git a/Wrappers/Python/ccpi/astra/__init__.py b/Wrappers/Python/ccpi/astra/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/Wrappers/Python/ccpi/astra/__init__.py
diff --git a/Wrappers/Python/ccpi/astra/astra_ops.py b/Wrappers/Python/ccpi/astra/astra_ops.py
new file mode 100644
index 0000000..45b8238
--- /dev/null
+++ b/Wrappers/Python/ccpi/astra/astra_ops.py
@@ -0,0 +1,132 @@
+# -*- coding: utf-8 -*-
+#    This work is independent part of the Core Imaging Library developed by
+#    Visual Analytics and Imaging System Group of the Science Technology
+#    Facilities Council, STFC
+#    This program 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.
+#
+#    This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from ccpi.reconstruction.ops import Operator
+import numpy
+import astra
+from ccpi.framework import AcquisitionData, ImageData, DataContainer
+from ccpi.reconstruction.ops import PowerMethodNonsquare
+from ccpi.astra.astra_processors import AstraForwardProjector, AstraBackProjector, \
+     AstraForwardProjectorMC, AstraBackProjectorMC
+
+class AstraProjectorSimple(Operator):
+    """ASTRA projector modified to use DataSet and geometry."""
+    def __init__(self, geomv, geomp, device):
+        super(AstraProjectorSimple, self).__init__()
+        
+        # Store volume and sinogram geometries.
+        self.sinogram_geometry = geomp
+        self.volume_geometry = geomv
+        
+        self.fp = AstraForwardProjector(volume_geometry=geomv,
+                                        sinogram_geometry=geomp,
+                                        proj_id=None,
+                                        device=device)
+        
+        self.bp = AstraBackProjector(volume_geometry=geomv,
+                                        sinogram_geometry=geomp,
+                                        proj_id=None,
+                                        device=device)
+                
+        # Initialise empty for singular value.
+        self.s1 = None
+    
+    def direct(self, IM):
+        self.fp.set_input(IM)
+        out = self.fp.get_output()
+        return out
+    
+    def adjoint(self, DATA):
+        self.bp.set_input(DATA)
+        out = self.bp.get_output()
+        return out
+    
+    #def delete(self):
+    #    astra.data2d.delete(self.proj_id)
+    
+    def get_max_sing_val(self):
+        self.s1, sall, svec = PowerMethodNonsquare(self,10)
+        return self.s1
+    
+    def size(self):
+        # Only implemented for 2D
+        return ( (self.sinogram_geometry.angles.size, \
+                  self.sinogram_geometry.pixel_num_h), \
+                 (self.volume_geometry.voxel_num_x, \
+                  self.volume_geometry.voxel_num_y) )
+    
+    def create_image_data(self):
+        inputsize = self.size()[1]
+        return DataContainer(numpy.random.randn(inputsize[0],
+                                                inputsize[1]))
+
+
+class AstraProjectorMC(Operator):
+    """ASTRA Multichannel projector"""
+    def __init__(self, geomv, geomp, device):
+        super(AstraProjectorMC, self).__init__()
+        
+        # Store volume and sinogram geometries.
+        self.sinogram_geometry = geomp
+        self.volume_geometry = geomv
+        
+        self.fp = AstraForwardProjectorMC(volume_geometry=geomv,
+                                        sinogram_geometry=geomp,
+                                        proj_id=None,
+                                        device=device)
+        
+        self.bp = AstraBackProjectorMC(volume_geometry=geomv,
+                                        sinogram_geometry=geomp,
+                                        proj_id=None,
+                                        device=device)
+                
+        # Initialise empty for singular value.
+        self.s1 = None
+    
+    def direct(self, IM):
+        self.fp.set_input(IM)
+        out = self.fp.get_output()
+        return out
+    
+    def adjoint(self, DATA):
+        self.bp.set_input(DATA)
+        out = self.bp.get_output()
+        return out
+    
+    #def delete(self):
+    #    astra.data2d.delete(self.proj_id)
+    
+    def get_max_sing_val(self):
+        if self.s1 is None:
+            self.s1, sall, svec = PowerMethodNonsquare(self,10)
+            return self.s1
+        else:
+            return self.s1
+    
+    def size(self):
+        # Only implemented for 2D
+        return ( (self.sinogram_geometry.angles.size, \
+                  self.sinogram_geometry.pixel_num_h), \
+                 (self.volume_geometry.voxel_num_x, \
+                  self.volume_geometry.voxel_num_y) )
+    
+    def create_image_data(self):
+        inputsize = self.size()[1]
+        return DataContainer(numpy.random.randn(self.volume_geometry.channels,
+                                                inputsize[0],
+                                                inputsize[1]))
+    
diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py
new file mode 100644
index 0000000..3de43bf
--- /dev/null
+++ b/Wrappers/Python/ccpi/astra/astra_processors.py
@@ -0,0 +1,205 @@
+from ccpi.framework import DataSetProcessor, ImageData, AcquisitionData
+from ccpi.astra.astra_utils import convert_geometry_to_astra
+import astra
+
+
+class AstraForwardProjector(DataSetProcessor):
+    '''AstraForwardProjector
+    
+    Forward project ImageData to AcquisitionData using ASTRA proj_id.
+    
+    Input: ImageData
+    Parameter: proj_id
+    Output: AcquisitionData
+    '''
+    
+    def __init__(self,
+                 volume_geometry=None,
+                 sinogram_geometry=None,
+                 proj_id=None,
+                 device='cpu'):
+        kwargs = {
+                  'volume_geometry'  : volume_geometry,
+                  'sinogram_geometry'  : sinogram_geometry,
+                  'proj_id'  : proj_id,
+                  'device'  : device
+                  }
+        
+        #DataSetProcessor.__init__(self, **kwargs)
+        super(AstraForwardProjector, self).__init__(**kwargs)
+        
+        self.set_ImageGeometry(volume_geometry)
+        self.set_AcquisitionGeometry(sinogram_geometry)
+        
+        # Set up ASTRA Volume and projection geometry, not to be stored in self
+        vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+                                                        self.sinogram_geometry)
+        
+        # ASTRA projector, to be stored
+        if device == 'cpu':
+            # Note that 'line' only one option
+            if self.sinogram_geometry.geom_type == 'parallel':
+                self.set_projector(astra.create_projector('line', proj_geom, vol_geom) )
+            elif self.sinogram_geometry.geom_type == 'cone':
+                self.set_projector(astra.create_projector('line_fanflat', proj_geom, vol_geom) )
+            else:
+                NotImplemented    
+        elif device == 'gpu':
+            self.set_projector(astra.create_projector('cuda', proj_geom, vol_geom) )
+        else:
+            NotImplemented
+    
+    def check_input(self, dataset):
+        if dataset.number_of_dimensions == 3 or\
+           dataset.number_of_dimensions == 2:
+               return True
+        else:
+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+                             .format(dataset.number_of_dimensions))
+    
+    def set_projector(self, proj_id):
+        self.proj_id = proj_id
+    
+    def set_ImageGeometry(self, volume_geometry):
+        self.volume_geometry = volume_geometry
+    
+    def set_AcquisitionGeometry(self, sinogram_geometry):
+        self.sinogram_geometry = sinogram_geometry
+    
+    def process(self):
+        IM = self.get_input()
+        DATA = AcquisitionData(geometry=self.sinogram_geometry)
+        #sinogram_id, DATA = astra.create_sino( IM.as_array(), 
+        #                           self.proj_id)
+        sinogram_id, DATA.array = astra.create_sino(IM.as_array(), 
+                                                           self.proj_id)
+        astra.data2d.delete(sinogram_id)
+        #return AcquisitionData(array=DATA, geometry=self.sinogram_geometry)
+        return DATA
+
+class AstraBackProjector(DataSetProcessor):
+    '''AstraBackProjector
+    
+    Back project AcquisitionData to ImageData using ASTRA proj_id.
+    
+    Input: AcquisitionData
+    Parameter: proj_id
+    Output: ImageData
+    '''
+    
+    def __init__(self,
+                 volume_geometry=None,
+                 sinogram_geometry=None,
+                 proj_id=None,
+                 device='cpu'):
+        kwargs = {
+                  'volume_geometry'  : volume_geometry,
+                  'sinogram_geometry'  : sinogram_geometry,
+                  'proj_id'  : proj_id,
+                  'device'  : device
+                  }
+        
+        #DataSetProcessor.__init__(self, **kwargs)
+        super(AstraBackProjector, self).__init__(**kwargs)
+        
+        self.set_ImageGeometry(volume_geometry)
+        self.set_AcquisitionGeometry(sinogram_geometry)
+        
+                # Set up ASTRA Volume and projection geometry, not to be stored in self
+        vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+                                                        self.sinogram_geometry)
+        
+        # ASTRA projector, to be stored
+        if device == 'cpu':
+            # Note that 'line' only one option
+            if self.sinogram_geometry.geom_type == 'parallel':
+                self.set_projector(astra.create_projector('line', proj_geom, vol_geom) )
+            elif self.sinogram_geometry.geom_type == 'cone':
+                self.set_projector(astra.create_projector('line_fanflat', proj_geom, vol_geom) )
+            else:
+                NotImplemented 
+        elif device == 'gpu':
+            self.set_projector(astra.create_projector('cuda', proj_geom, vol_geom) )
+        else:
+            NotImplemented
+    
+    def check_input(self, dataset):
+        if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:
+               return True
+        else:
+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+                             .format(dataset.number_of_dimensions))
+    
+    def set_projector(self, proj_id):
+        self.proj_id = proj_id
+        
+    def set_ImageGeometry(self, volume_geometry):
+        self.volume_geometry = volume_geometry
+        
+    def set_AcquisitionGeometry(self, sinogram_geometry):
+        self.sinogram_geometry = sinogram_geometry
+    
+    def process(self):
+        DATA = self.get_input()
+        IM = ImageData(geometry=self.volume_geometry)
+        rec_id, IM.array = astra.create_backprojection(DATA.as_array(),
+                            self.proj_id)
+        astra.data2d.delete(rec_id)
+        return IM
+
+class AstraForwardProjectorMC(AstraForwardProjector):
+    '''AstraForwardProjector Multi channel
+    
+    Forward project ImageData to AcquisitionDataSet using ASTRA proj_id.
+    
+    Input: ImageDataSet
+    Parameter: proj_id
+    Output: AcquisitionData
+    '''
+    def check_input(self, dataset):
+        if dataset.number_of_dimensions == 2 or \
+           dataset.number_of_dimensions == 3 or \
+           dataset.number_of_dimensions == 4:
+               return True
+        else:
+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+                             .format(dataset.number_of_dimensions))
+    def process(self):
+        IM = self.get_input()
+        #create the output AcquisitionData
+        DATA = AcquisitionData(geometry=self.sinogram_geometry)
+        
+        for k in range(DATA.geometry.channels):
+            sinogram_id, DATA.as_array()[k] = astra.create_sino(IM.as_array()[k], 
+                                                           self.proj_id)
+            astra.data2d.delete(sinogram_id)
+        return DATA
+
+class AstraBackProjectorMC(AstraBackProjector):
+    '''AstraBackProjector Multi channel
+    
+    Back project AcquisitionData to ImageData using ASTRA proj_id.
+    
+    Input: AcquisitionData
+    Parameter: proj_id
+    Output: ImageData
+    '''
+    def check_input(self, dataset):
+        if dataset.number_of_dimensions == 2 or \
+           dataset.number_of_dimensions == 3 or \
+           dataset.number_of_dimensions == 4:
+               return True
+        else:
+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+                             .format(dataset.number_of_dimensions))
+    def process(self):
+        DATA = self.get_input()
+        
+        IM = ImageData(geometry=self.volume_geometry)
+        
+        for k in range(IM.geometry.channels):
+            rec_id, IM.as_array()[k] = astra.create_backprojection(
+                    DATA.as_array()[k], 
+                    self.proj_id)
+            astra.data2d.delete(rec_id)
+        return IM
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/astra/astra_utils.py b/Wrappers/Python/ccpi/astra/astra_utils.py
new file mode 100644
index 0000000..9f8fe46
--- /dev/null
+++ b/Wrappers/Python/ccpi/astra/astra_utils.py
@@ -0,0 +1,61 @@
+import astra
+
+def convert_geometry_to_astra(volume_geometry, sinogram_geometry):
+    # Set up ASTRA Volume and projection geometry, not stored
+    if sinogram_geometry.dimension == '2D':
+        vol_geom = astra.create_vol_geom(volume_geometry.voxel_num_x, 
+                                         volume_geometry.voxel_num_y, 
+                                         volume_geometry.get_min_x(), 
+                                         volume_geometry.get_max_x(), 
+                                         volume_geometry.get_min_y(), 
+                                         volume_geometry.get_max_y())
+        
+        if sinogram_geometry.geom_type == 'parallel':
+            proj_geom = astra.create_proj_geom('parallel',
+                                               sinogram_geometry.pixel_size_h,
+                                               sinogram_geometry.pixel_num_h,
+                                               sinogram_geometry.angles)
+        elif sinogram_geometry.geom_type == 'cone':
+            proj_geom = astra.create_proj_geom('fanflat',
+                                               sinogram_geometry.pixel_size_h,
+                                               sinogram_geometry.pixel_num_h,
+                                               sinogram_geometry.angles,
+                                               sinogram_geometry.dist_source_center,
+                                               sinogram_geometry.dist_center_detector)
+        else:
+            NotImplemented
+            
+    elif sinogram_geometry.dimension == '3D':
+        vol_geom = astra.create_vol_geom(volume_geometry.voxel_num_x, 
+                                         volume_geometry.voxel_num_y, 
+                                         volume_geometry.voxel_num_z, 
+                                         volume_geometry.get_min_x(), 
+                                         volume_geometry.get_max_x(), 
+                                         volume_geometry.get_min_y(), 
+                                         volume_geometry.get_max_y(), 
+                                         volume_geometry.get_min_z(), 
+                                         volume_geometry.get_max_z())
+        
+        if sinogram_geometry.geom_type == 'parallel':
+            proj_geom = astra.create_proj_geom('parallel3d',
+                                               sinogram_geometry.pixel_size_h,
+                                               sinogram_geometry.pixel_size_v,
+                                               sinogram_geometry.pixel_num_v,
+                                               sinogram_geometry.pixel_num_h,
+                                               sinogram_geometry.angles)
+        elif sinogram_geometry.geom_type == 'cone':
+            proj_geom = astra.create_proj_geom('cone',
+                                               sinogram_geometry.pixel_size_h,
+                                               sinogram_geometry.pixel_size_v,
+                                               sinogram_geometry.pixel_num_v,
+                                               sinogram_geometry.pixel_num_h,
+                                               sinogram_geometry.angles,
+                                               sinogram_geometry.dist_source_center,
+                                               sinogram_geometry.dist_center_detector)
+        else:
+            NotImplemented
+            
+    else:
+        NotImplemented
+    
+    return vol_geom, proj_geom
\ No newline at end of file
diff --git a/Wrappers/Python/conda-recipe/bld.bat b/Wrappers/Python/conda-recipe/bld.bat
new file mode 100755
index 0000000..4b4c7f5
--- /dev/null
+++ b/Wrappers/Python/conda-recipe/bld.bat
@@ -0,0 +1,10 @@
+IF NOT DEFINED CIL_VERSION (
+ECHO CIL_VERSION Not Defined.
+exit 1
+)
+
+ROBOCOPY /E "%RECIPE_DIR%\.." "%SRC_DIR%"
+
+%PYTHON% setup.py build_py
+%PYTHON% setup.py install
+if errorlevel 1 exit 1
diff --git a/Wrappers/Python/conda-recipe/build.sh b/Wrappers/Python/conda-recipe/build.sh
new file mode 100755
index 0000000..5dd97b0
--- /dev/null
+++ b/Wrappers/Python/conda-recipe/build.sh
@@ -0,0 +1,9 @@
+if [ -z "$CIL_VERSION" ]; then
+    echo "Need to set CIL_VERSION"
+    exit 1
+fi  
+mkdir ${SRC_DIR}/ccpi
+cp -r "${RECIPE_DIR}/../../../" ${SRC_DIR}/ccpi
+
+cd ${SRC_DIR}/ccpi/Wrappers/Python
+$PYTHON setup.py install
diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml
new file mode 100755
index 0000000..12e83e2
--- /dev/null
+++ b/Wrappers/Python/conda-recipe/meta.yaml
@@ -0,0 +1,28 @@
+package:
+  name: ccpi-astra
+  version: {{ environ['CIL_VERSION'] }}
+
+
+build:
+  preserve_egg_dir: False
+  script_env:
+    - CIL_VERSION   
+#  number: 0
+  
+requirements:
+  build:
+    - python
+    - numpy
+    - setuptools
+
+  run:
+    - python
+    - numpy
+    - scipy
+    - ccpi-common
+    - astra-toolbox
+	
+about:
+  home: http://www.ccpi.ac.uk
+  license:  GPLv3
+  summary: 'CCPi Toolbox'
diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py
new file mode 100755
index 0000000..8985ddd
--- /dev/null
+++ b/Wrappers/Python/setup.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#   This work is part of the Core Imaging Library developed by
+#   Visual Analytics and Imaging System Group of the Science Technology
+#   Facilities Council, STFC
+
+#   Copyright 2018 Edoardo Pasca
+
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+
+#       http://www.apache.org/licenses/LICENSE-2.0
+
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+from distutils.core import setup
+import os
+import sys
+
+cil_version=os.environ['CIL_VERSION']
+if  cil_version == '':
+    print("Please set the environmental variable CIL_VERSION")
+    sys.exit(1)
+
+cil_version=os.environ['CIL_VERSION']
+if  cil_version == '':
+    print("Please set the environmental variable CIL_VERSION")
+    sys.exit(1)
+
+setup(
+    name="ccpi-common",
+    version=cil_version,
+    packages=['ccpi' , 'ccpi.astra'],
+
+    # Project uses reStructuredText, so ensure that the docutils get
+    # installed or upgraded on the target machine
+    #install_requires=['docutils>=0.3'],
+
+#    package_data={
+#        # If any package contains *.txt or *.rst files, include them:
+#        '': ['*.txt', '*.rst'],
+#        # And include any *.msg files found in the 'hello' package, too:
+#        'hello': ['*.msg'],
+#    },
+    # zip_safe = False,
+
+    # metadata for upload to PyPI
+    author="Edoardo Pasca",
+    author_email="edoardo.pasca@stfc.ac.uk",
+    description='CCPi Core Imaging Library - Astra bindings',
+    license="Apache v2.0",
+    keywords="Python Framework",
+    url="http://www.ccpi.ac.uk",   # project home page, if any
+
+    # could also include long_description, download_url, classifiers, etc.
+)
diff --git a/Wrappers/Python/wip/DemoRecIP.py b/Wrappers/Python/wip/DemoRecIP.py
new file mode 100755
index 0000000..442e40e
--- /dev/null
+++ b/Wrappers/Python/wip/DemoRecIP.py
@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Reading multi-channel data and reconstruction using FISTA modular
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+#import sys
+#sys.path.append('../../../data/')
+from read_IPdata import read_IPdata
+
+from ccpi.astra.astra_ops import AstraProjectorSimple, AstraProjectorMC
+from ccpi.reconstruction.funcs import Norm2sq, Norm1, BaseFunction
+from ccpi.reconstruction.algs import FISTA
+#from ccpi.reconstruction.funcs import BaseFunction
+
+from ccpi.framework import ImageData, AcquisitionData, AcquisitionGeometry, ImageGeometry
+
+# read IP paper data into a dictionary
+dataDICT = read_IPdata('..\..\..\data\IP_data70channels.mat')
+
+# Set ASTRA Projection-backprojection class (fan-beam geometry)
+DetWidth = dataDICT.get('im_size')[0] * dataDICT.get('det_width')[0] / \
+            dataDICT.get('detectors_numb')[0]
+SourceOrig = dataDICT.get('im_size')[0] * dataDICT.get('src_to_rotc')[0] / \
+            dataDICT.get('dom_width')[0]
+OrigDetec = dataDICT.get('im_size')[0] * \
+            (dataDICT.get('src_to_det')[0] - dataDICT.get('src_to_rotc')[0]) /\
+            dataDICT.get('dom_width')[0]
+
+N = dataDICT.get('im_size')[0]
+
+vg = ImageGeometry(voxel_num_x=dataDICT.get('im_size')[0],
+                   voxel_num_y=dataDICT.get('im_size')[0],
+                   channels=1)
+
+pg = AcquisitionGeometry('cone',
+                         '2D',
+                         angles=(np.pi/180)*dataDICT.get('theta')[0],
+                         pixel_num_h=dataDICT.get('detectors_numb')[0],
+                         pixel_size_h=DetWidth,
+                         dist_source_center=SourceOrig, 
+                         dist_center_detector=OrigDetec,
+                         channels=1)
+
+
+sino = dataDICT.get('data_norm')[0][:,:,34] # select mid-channel 
+b = AcquisitionData(sino,geometry=pg)
+
+# Initial guess
+x_init = ImageData(np.zeros((N, N)),geometry=vg)
+
+
+
+
+
+Aop = AstraProjectorSimple(vg,pg,'gpu')
+f = Norm2sq(Aop,b,c=0.5)
+
+# Run FISTA for least squares without regularization
+opt = {'tol': 1e-4, 'iter': 10}
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
+
+plt.imshow(x_fista0.array)
+plt.show()
+
+# Now least squares plus 1-norm regularization
+g1 = Norm1(10)
+
+# Run FISTA for least squares plus 1-norm function.
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt)
+
+plt.imshow(x_fista1.array)
+plt.show()
+
+# Multiple channels
+sino_mc = dataDICT.get('data_norm')[0][:,:,32:37] # select mid-channel 
+
+vg_mc = ImageGeometry(voxel_num_x=dataDICT.get('im_size')[0],
+                   voxel_num_y=dataDICT.get('im_size')[0],
+                   channels=5)
+
+pg_mc = AcquisitionGeometry('cone',
+                         '2D',
+                         angles=(np.pi/180)*dataDICT.get('theta')[0],
+                         pixel_num_h=dataDICT.get('detectors_numb')[0],
+                         pixel_size_h=DetWidth,
+                         dist_source_center=SourceOrig, 
+                         dist_center_detector=OrigDetec,
+                         channels=5)
+
+b_mc = AcquisitionData(np.transpose(sino_mc,(2,0,1)),
+                       geometry=pg_mc, 
+                       dimension_labels=("channel","angle","horizontal"))
+
+# ASTRA operator using volume and sinogram geometries
+Aop_mc = AstraProjectorMC(vg_mc, pg_mc, 'gpu')
+
+f_mc = Norm2sq(Aop_mc,b_mc,c=0.5)
+
+# Initial guess
+x_init_mc = ImageData(np.zeros((5, N, N)),geometry=vg_mc)
+
+
+x_fista0_mc, it0_mc, timing0_mc, criter0_mc = FISTA(x_init_mc, f_mc, None, opt)
+
+plt.imshow(x_fista0_mc.as_array()[4,:,:])
+plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/astra_test.py b/Wrappers/Python/wip/astra_test.py
new file mode 100755
index 0000000..c0a7359
--- /dev/null
+++ b/Wrappers/Python/wip/astra_test.py
@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-

+"""

+Created on Wed Mar  7 15:07:16 2018

+

+@author: ofn77899

+"""

+

+# -----------------------------------------------------------------------

+# Copyright: 2010-2018, imec Vision Lab, University of Antwerp

+#            2013-2018, CWI, Amsterdam

+#

+# Contact: astra@astra-toolbox.com

+# Website: http://www.astra-toolbox.com/

+#

+# This file is part of the ASTRA Toolbox.

+#

+#

+# 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 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 ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.

+#

+# -----------------------------------------------------------------------

+

+import astra

+import numpy as np

+

+vol_geom = astra.create_vol_geom(256, 256)

+proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False))

+

+# For CPU-based algorithms, a "projector" object specifies the projection

+# model used. In this case, we use the "strip" model.

+proj_id = astra.create_projector('strip', proj_geom, vol_geom)

+

+# Create a sinogram from a phantom

+import scipy.io

+P = scipy.io.loadmat('phantom.mat')['phantom256']

+sinogram_id, sinogram = astra.create_sino(P, proj_id)

+

+import pylab

+pylab.gray()

+pylab.figure(1)

+pylab.imshow(P)

+pylab.figure(2)

+pylab.imshow(sinogram)

+

+# Create a data object for the reconstruction

+rec_id = astra.data2d.create('-vol', vol_geom)

+

+# Set up the parameters for a reconstruction algorithm using the CPU

+# The main difference with the configuration of a GPU algorithm is the

+# extra ProjectorId setting.

+cfg = astra.astra_dict('SIRT')

+cfg['ReconstructionDataId'] = rec_id

+cfg['ProjectionDataId'] = sinogram_id

+cfg['ProjectorId'] = proj_id

+

+# Available algorithms:

+# ART, SART, SIRT, CGLS, FBP

+

+

+# Create the algorithm object from the configuration structure

+alg_id = astra.algorithm.create(cfg)

+

+# Run 20 iterations of the algorithm

+# This will have a runtime in the order of 10 seconds.

+astra.algorithm.run(alg_id, 20)

+

+# Get the result

+rec = astra.data2d.get(rec_id)

+pylab.figure(3)

+pylab.imshow(rec)

+pylab.show()

+

+# Clean up.

+astra.algorithm.delete(alg_id)

+astra.data2d.delete(rec_id)

+astra.data2d.delete(sinogram_id)

+astra.projector.delete(proj_id)

diff --git a/Wrappers/Python/wip/demo_sophiabeads.py b/Wrappers/Python/wip/demo_sophiabeads.py
new file mode 100755
index 0000000..e3c7f3a
--- /dev/null
+++ b/Wrappers/Python/wip/demo_sophiabeads.py
@@ -0,0 +1,65 @@
+
+from ccpi.io.reader import XTEKReader
+import numpy as np
+import matplotlib.pyplot as plt
+from ccpi.framework import ImageGeometry, AcquisitionGeometry, AcquisitionData, ImageData
+from ccpi.astra.astra_ops import AstraProjectorSimple
+from ccpi.reconstruction.algs import CGLS
+
+# Set up reader object and read the data
+datareader = XTEKReader("C:/Users/mbbssjj2/Documents/SophiaBeads_256_averaged/SophiaBeads_256_averaged.xtekct")
+data = datareader.getAcquisitionData()
+
+# Extract central slice, scale and negative-log transform
+sino = -np.log(data.as_array()[:,:,1000]/60000.0)
+
+# Apply centering correction by zero padding, amount found manually
+cor_pad = 30
+sino_pad = np.zeros((sino.shape[0],sino.shape[1]+cor_pad))
+sino_pad[:,cor_pad:] = sino
+
+# Simple beam hardening correction as done in SophiaBeads coda
+sino_pad = sino_pad**2
+
+# Extract AcquisitionGeometry for central slice for 2D fanbeam reconstruction
+ag2d = AcquisitionGeometry('cone',
+                          '2D',
+                          angles=-np.pi/180*data.geometry.angles,
+                          pixel_num_h=data.geometry.pixel_num_h + cor_pad,
+                          pixel_size_h=data.geometry.pixel_size_h,
+                          dist_source_center=data.geometry.dist_source_center, 
+                          dist_center_detector=data.geometry.dist_center_detector)
+
+# Set up AcquisitionData object for central slice 2D fanbeam
+data2d = AcquisitionData(sino_pad,geometry=ag2d)
+
+# Choose the number of voxels to reconstruct onto as number of detector pixels
+N = data.geometry.pixel_num_h
+
+# Geometric magnification
+mag = (np.abs(data.geometry.dist_center_detector) + \
+      np.abs(data.geometry.dist_source_center)) / \
+      np.abs(data.geometry.dist_source_center)
+
+# Voxel size is detector pixel size divided by mag
+voxel_size_h = data.geometry.pixel_size_h / mag
+
+# Construct the appropriate ImageGeometry
+ig2d = ImageGeometry(voxel_num_x=N,
+                   voxel_num_y=N,
+                   voxel_size_x=voxel_size_h, 
+                   voxel_size_y=voxel_size_h)
+
+# Set up the Projector (AcquisitionModel) using ASTRA on GPU
+Aop = AstraProjectorSimple(ig2d, ag2d,"gpu")
+
+# Set initial guess for CGLS reconstruction
+x_init = ImageData(np.zeros((N,N)),geometry=ig2d)
+
+# Run CGLS reconstruction
+num_iter = 50
+x, it, timing, criter = CGLS(Aop,data2d,num_iter,x_init)
+
+# Display reconstruction
+plt.figure()
+plt.imshow(x.as_array())
\ No newline at end of file
diff --git a/Wrappers/Python/wip/desktop.ini b/Wrappers/Python/wip/desktop.ini
new file mode 100755
index 0000000..bb9f3d6
--- /dev/null
+++ b/Wrappers/Python/wip/desktop.ini
@@ -0,0 +1,4 @@
+[ViewState]

+Mode=

+Vid=

+FolderType=Generic

diff --git a/Wrappers/Python/wip/simple_demo_astra.py b/Wrappers/Python/wip/simple_demo_astra.py
new file mode 100755
index 0000000..3365504
--- /dev/null
+++ b/Wrappers/Python/wip/simple_demo_astra.py
@@ -0,0 +1,189 @@
+#import sys
+#sys.path.append("..")
+
+from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry
+from ccpi.reconstruction.algs import FISTA, FBPD, CGLS
+from ccpi.reconstruction.funcs import Norm2sq, Norm1 , TV2D
+from ccpi.astra.astra_ops import AstraProjectorSimple
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+test_case = 1   # 1=parallel2D, 2=cone2D
+
+# Set up phantom
+N = 128
+
+
+vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
+Phantom = ImageData(geometry=vg)
+
+x = Phantom.as_array()
+x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
+x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
+
+plt.imshow(x)
+plt.show()
+
+# Set up measurement geometry
+angles_num = 20; # angles number
+
+if test_case==1:
+    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
+elif test_case==2:
+    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
+else:
+    NotImplemented
+
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
+
+# Parallelbeam geometry test
+if test_case==1:
+    pg = AcquisitionGeometry('parallel',
+                          '2D',
+                          angles,
+                          det_num,det_w)
+elif test_case==2:
+    pg = AcquisitionGeometry('cone',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          dist_source_center=SourceOrig, 
+                          dist_center_detector=OrigDetec)
+
+# ASTRA operator using volume and sinogram geometries
+Aop = AstraProjectorSimple(vg, pg, 'cpu')
+
+# Unused old astra projector without geometry
+# Aop_old = AstraProjector(det_w, det_num, SourceOrig, 
+#                      OrigDetec, angles, 
+#                      N,'fanbeam','gpu') 
+
+# Try forward and backprojection
+b = Aop.direct(Phantom)
+out2 = Aop.adjoint(b)
+
+plt.imshow(b.array)
+plt.show()
+
+plt.imshow(out2.array)
+plt.show()
+
+# Create least squares object instance with projector and data.
+f = Norm2sq(Aop,b,c=0.5)
+
+# Initial guess
+x_init = ImageData(np.zeros(x.shape),geometry=vg)
+
+# Run FISTA for least squares without regularization
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None)
+
+plt.imshow(x_fista0.array)
+plt.title('FISTA0')
+plt.show()
+
+# Now least squares plus 1-norm regularization
+lam = 0.1
+g0 = Norm1(lam)
+
+# Run FISTA for least squares plus 1-norm function.
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0)
+
+plt.imshow(x_fista1.array)
+plt.title('FISTA1')
+plt.show()
+
+plt.semilogy(criter1)
+plt.show()
+
+# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
+opt = {'tol': 1e-4, 'iter': 100}
+x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
+
+plt.imshow(x_fbpd1.array)
+plt.title('FBPD1')
+plt.show()
+
+plt.semilogy(criter_fbpd1)
+plt.show()
+
+# Now FBPD for least squares plus TV
+#lamtv = 1
+#gtv = TV2D(lamtv)
+
+#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt)
+
+#plt.imshow(x_fbpdtv.array)
+#plt.show()
+
+#plt.semilogy(criter_fbpdtv)
+#plt.show()
+
+
+# Run CGLS, which should agree with the FISTA0
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt )
+
+plt.imshow(x_CGLS.array)
+plt.title('CGLS')
+#plt.title('CGLS recon, compare FISTA0')
+plt.show()
+
+plt.semilogy(criter_CGLS)
+plt.title('CGLS criterion')
+plt.show()
+
+
+#%%
+
+clims = (0,1)
+cols = 3
+rows = 2
+current = 1
+fig = plt.figure()
+# projections row
+a=fig.add_subplot(rows,cols,current)
+a.set_title('phantom {0}'.format(np.shape(Phantom.as_array())))
+
+imgplot = plt.imshow(Phantom.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA0')
+imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA1')
+imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FBPD1')
+imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('CGLS')
+imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1])
+
+#current = current + 1
+#a=fig.add_subplot(rows,cols,current)
+#a.set_title('FBPD TV')
+#imgplot = plt.imshow(x_fbpdtv.as_array(),vmin=clims[0],vmax=clims[1])
+
+fig = plt.figure()
+# projections row
+b=fig.add_subplot(1,1,1)
+b.set_title('criteria')
+imgplot = plt.loglog(criter0 , label='FISTA0')
+imgplot = plt.loglog(criter1 , label='FISTA1')
+imgplot = plt.loglog(criter_fbpd1, label='FBPD1')
+imgplot = plt.loglog(criter_CGLS, label='CGLS')
+#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV')
+b.legend(loc='right')
+plt.show()
+#%%
\ No newline at end of file
diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py
new file mode 100755
index 0000000..f77a678
--- /dev/null
+++ b/Wrappers/Python/wip/simple_mc_demo.py
@@ -0,0 +1,142 @@
+#import sys
+#sys.path.append("..")
+
+from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry
+from ccpi.reconstruction.algs import FISTA
+from ccpi.reconstruction.funcs import Norm2sq, Norm1
+from ccpi.astra.astra_ops import AstraProjectorMC
+
+import numpy
+import matplotlib.pyplot as plt
+
+test_case = 1   # 1=parallel2D, 2=cone2D
+
+# Set up phantom
+N = 128
+M = 100
+numchannels = 3
+
+vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels)
+Phantom = ImageData(geometry=vg)
+
+x = Phantom.as_array()
+x[0 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 1.0
+x[0 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.0
+
+x[1 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 0.7
+x[1 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 1.2
+
+x[2 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 1.5
+x[2 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.2
+
+#x = numpy.zeros((N,N,1,numchannels))
+
+#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0
+#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0
+
+#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7
+#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2
+
+#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5
+#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2
+
+f, axarr = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarr[k].imshow(x[k],vmin=0,vmax=2.5)
+plt.show()
+
+#vg = ImageGeometry(N,N,None, 1,1,None,channels=numchannels)
+
+
+#Phantom = ImageData(x,geometry=vg)
+
+# Set up measurement geometry
+angles_num = 20; # angles number
+
+if test_case==1:
+    angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False)
+elif test_case==2:
+    angles = numpy.linspace(0,2*numpy.pi,angles_num,endpoint=False)
+else:
+    NotImplemented
+
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
+
+# Parallelbeam geometry test
+if test_case==1:
+    pg = AcquisitionGeometry('parallel',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          channels=numchannels)
+elif test_case==2:
+    pg = AcquisitionGeometry('cone',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          dist_source_center=SourceOrig, 
+                          dist_center_detector=OrigDetec,
+                          channels=numchannels)
+
+# ASTRA operator using volume and sinogram geometries
+Aop = AstraProjectorMC(vg, pg, 'cpu')
+
+
+# Try forward and backprojection
+b = Aop.direct(Phantom)
+
+fb, axarrb = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrb[k].imshow(b.as_array()[k],vmin=0,vmax=250)
+plt.show()
+
+out2 = Aop.adjoint(b)
+
+fo, axarro = plt.subplots(1,numchannels)
+for k in range(3):
+    axarro[k].imshow(out2.as_array()[k],vmin=0,vmax=3500)
+plt.show()
+
+# Create least squares object instance with projector and data.
+f = Norm2sq(Aop,b,c=0.5)
+
+# Initial guess
+x_init = ImageData(numpy.zeros(x.shape),geometry=vg)
+
+# FISTA options
+opt = {'tol': 1e-4, 'iter': 200}
+
+# Run FISTA for least squares without regularization
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
+
+
+ff0, axarrf0 = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrf0[k].imshow(x_fista0.as_array()[k],vmin=0,vmax=2.5)
+plt.show()
+
+plt.semilogy(criter0)
+plt.title('Criterion vs iterations, least squares')
+plt.show()
+
+# Now least squares plus 1-norm regularization
+lam = 0.1
+g0 = Norm1(lam)
+
+
+# Run FISTA for least squares plus 1-norm function.
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
+
+ff1, axarrf1 = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrf1[k].imshow(x_fista1.as_array()[k],vmin=0,vmax=2.5)
+plt.show()
+
+plt.semilogy(criter1)
+plt.title('Criterion vs iterations, least squares plus 1-norm regu')
+plt.show()
\ No newline at end of file