Move astra (#85)

* removed astra stuff

* removed astra and added dependency on ccpi-reconstruction

* removed astra test

* renamed to ccpi-framework

* little changes to demo

7 files changed, 20 insertions, 608 deletions

diff --git a/Wrappers/Python/ccpi/astra/astra_ops.py b/Wrappers/Python/ccpi/astra/astra_ops.py
deleted file mode 100644
index 45b8238..0000000
--- a/Wrappers/Python/ccpi/astra/astra_ops.py
+++ /dev/null
@@ -1,132 +0,0 @@
-# -*- 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
deleted file mode 100644
index 3de43bf..0000000
--- a/Wrappers/Python/ccpi/astra/astra_processors.py
+++ /dev/null
@@ -1,205 +0,0 @@
-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
deleted file mode 100644
index 9f8fe46..0000000
--- a/Wrappers/Python/ccpi/astra/astra_utils.py
+++ /dev/null
@@ -1,61 +0,0 @@
-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/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml
index fc5f09d..9813f78 100644
--- a/Wrappers/Python/conda-recipe/meta.yaml
+++ b/Wrappers/Python/conda-recipe/meta.yaml
@@ -1,5 +1,5 @@
 package:
-  name: ccpi-common
+  name: ccpi-framework
   version: {{ environ['CIL_VERSION'] }}
 
 
@@ -13,12 +13,11 @@ requirements:
   build:
     - python
     - numpy
-    - setuptools
+    - distutils
 
   run:
     - python
     - numpy
-    - scipy
 	
 about:
   home: http://www.ccpi.ac.uk
diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py
index a06b396..8dd16bd 100644
--- a/Wrappers/Python/setup.py
+++ b/Wrappers/Python/setup.py
@@ -35,7 +35,7 @@ if  cil_version == '':
 setup(
     name="ccpi-common",
     version=cil_version,
-    packages=['ccpi' , 'ccpi.reconstruction' , 'ccpi.astra', 'ccpi.io'],
+    packages=['ccpi' , 'ccpi.reconstruction' , 'ccpi.io'],
 
     # Project uses reStructuredText, so ensure that the docutils get
     # installed or upgraded on the target machine
diff --git a/Wrappers/Python/wip/simple_demo.py b/Wrappers/Python/wip/simple_demo.py
index b183794..d5d6826 100644
--- a/Wrappers/Python/wip/simple_demo.py
+++ b/Wrappers/Python/wip/simple_demo.py
@@ -19,8 +19,8 @@ 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)] = 1.0
-x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0
+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()
@@ -101,7 +101,7 @@ plt.semilogy(criter1)
 plt.show()
 
 # Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
-opt = {'tol': 1e-4, 'iter': 10000}
+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)
@@ -112,24 +112,24 @@ plt.semilogy(criter_fbpd1)
 plt.show()
 
 # Now FBPD for least squares plus TV
-lamtv = 1
-gtv = TV2D(lamtv)
+#lamtv = 1
+#gtv = TV2D(lamtv)
 
-x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt)
+#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt)
 
-plt.imshow(x_fbpdtv.array)
-plt.show()
+#plt.imshow(x_fbpdtv.array)
+#plt.show()
 
-plt.semilogy(criter_fbpdtv)
-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(Aop, b, 1000, x_init)
+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.title('CGLS recon, compare FISTA0')
 plt.show()
 
 plt.semilogy(criter_CGLS)
@@ -139,7 +139,7 @@ plt.show()
 
 #%%
 
-clims = (-0.5,2.5)
+clims = (0,1)
 cols = 3
 rows = 2
 current = 1
@@ -170,10 +170,10 @@ 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])
+#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
diff --git a/Wrappers/Python/wip/simple_demo_astra.py b/Wrappers/Python/wip/simple_demo_astra.py
deleted file mode 100755
index 3365504..0000000
--- a/Wrappers/Python/wip/simple_demo_astra.py
+++ /dev/null
@@ -1,189 +0,0 @@
-#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