Remove already moved files (#116)

4 files changed, 0 insertions, 666 deletions

diff --git a/Wrappers/Python/wip/demo_sophiabeads.py b/Wrappers/Python/wip/demo_sophiabeads.py
deleted file mode 100644
index e3c7f3a..0000000
--- a/Wrappers/Python/wip/demo_sophiabeads.py
+++ /dev/null
@@ -1,65 +0,0 @@
-
-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/simple_demo.py b/Wrappers/Python/wip/simple_demo.py
deleted file mode 100644
index d5d6826..0000000
--- a/Wrappers/Python/wip/simple_demo.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
diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py
deleted file mode 100755
index d2cd6f5..0000000
--- a/Wrappers/Python/wip/simple_demo_ccpi.py
+++ /dev/null
@@ -1,270 +0,0 @@
-#import sys
-#sys.path.append("..")
-
-from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry
-from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-from ccpi.optimisation.funcs import Norm2sq, Norm1 , TV2D
-#from ccpi.astra.astra_ops import AstraProjectorSimple
-from ccpi.reconstruction.ops import CCPiProjectorSimple
-from ccpi.reconstruction.parallelbeam import alg as pbalg
-from ccpi.processors import CCPiForwardProjector, CCPiBackwardProjector 
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-test_case = 1   # 1=parallel2D, 2=cone2D, 3=parallel3D
-
-# Set up phantom
-N = 128
-vert = 1
-# Set up measurement geometry
-angles_num = 20; # angles number
-det_w = 1.0
-det_num = N
-SourceOrig = 200
-OrigDetec = 0
-
-if test_case==1:
-    angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*180/np.pi
-    #nangles = angles_num
-    #angles = np.linspace(0,360, nangles, dtype=np.float32)
-
-elif test_case==2:
-    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
-elif test_case == 3:
-    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
-else:
-    NotImplemented
-
-
-def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter):
-    vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z)
-    Phantom_ccpi = ImageData(geometry=vg,
-                        dimension_labels=['horizontal_x','horizontal_y','vertical'])
-    #.subset(['horizontal_x','horizontal_y','vertical'])
-    # ask the ccpi code what dimensions it would like
-        
-    voxel_per_pixel = 1
-    geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(),
-                                                angles,
-                                                voxel_per_pixel )
-    
-    pg = AcquisitionGeometry('parallel',
-                              '3D',
-                              angles,
-                              geoms['n_h'],det_w,
-                              geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick
-                              )
-    
-    center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2
-    ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal'])
-    geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(),
-                                                angles,
-                                                center_of_rotation,
-                                                voxel_per_pixel )
-
-    #print (counter)
-    counter+=1
-    #print (geoms , geoms_i)
-    if counter < 4:
-        if (not ( geoms_i == geoms )):
-            print ("not equal and {0}".format(counter))
-            X = max(geoms['output_volume_x'], geoms_i['output_volume_x'])
-            Y = max(geoms['output_volume_y'], geoms_i['output_volume_y'])
-            Z = max(geoms['output_volume_z'], geoms_i['output_volume_z'])
-            return setupCCPiGeometries(X,Y,Z,angles, counter)
-        else:
-            print ("return geoms {0}".format(geoms))
-            return geoms
-    else:
-        print ("return geoms_i {0}".format(geoms_i))
-        return geoms_i
-
-geoms = setupCCPiGeometries(N,N,vert,angles,0)
-#%%
-#geoms = {'n_v': 12, 'output_volume_y': 128, 'n_angles': 20, 
-# 'output_volume_x': 128, 'output_volume_z': 12, 'n_h': 128}
-vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'],
-                   voxel_num_y=geoms['output_volume_y'], 
-                   voxel_num_z=geoms['output_volume_z'])
-Phantom = ImageData(geometry=vg,dimension_labels=['horizontal_x','horizontal_y','vertical']) + 0.1
-    
-
-#x = Phantom.as_array()
-i = 0
-while i < geoms['n_v']:
-    #x = Phantom.subset(vertical=i, dimensions=['horizontal_x','horizontal_y']).array
-    x = Phantom.subset(vertical=i).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)] = 0.98
-    Phantom.fill(x, vertical=i)
-    i += 1
-
-plt.imshow(Phantom.subset(vertical=0).as_array())
-plt.show()
-
-
-
-# Parallelbeam geometry test
-if test_case==1:
-    #Phantom_ccpi = Phantom.subset(['horizontal_x','horizontal_y','vertical'])
-    #Phantom_ccpi.geometry = vg.clone()
-    center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2
-        
-    pg = AcquisitionGeometry('parallel',
-                          '3D',
-                          angles,
-                          geoms['n_h'],det_w,
-                          geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick
-                          )
-elif test_case==2:
-    raise NotImplemented('cone beam projector not yet available')
-    pg = AcquisitionGeometry('cone',
-                          '2D',
-                          angles,
-                          det_num,
-                          det_w,
-                          vert, det_w, #2D in 3D is a slice 1 pixel thick 
-                          dist_source_center=SourceOrig, 
-                          dist_center_detector=OrigDetec)
-
-# ASTRA operator using volume and sinogram geometries
-#Aop = AstraProjectorSimple(vg, pg, 'cpu')
-Cop = CCPiProjectorSimple(vg, pg)
-
-# Try forward and backprojection
-b = Cop.direct(Phantom)
-out2 = Cop.adjoint(b)
-
-#%%
-for i in range(b.get_dimension_size('vertical')):
-    plt.imshow(b.subset(vertical=i).array)
-    #plt.imshow(Phantom.subset( vertical=i).array)
-    #plt.imshow(b.array[:,i,:])
-    plt.show()
-#%%
-
-plt.imshow(out2.subset( vertical=0).array)
-plt.show()
-
-# Create least squares object instance with projector and data.
-f = Norm2sq(Cop,b,c=0.5)
-
-# Initial guess
-x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y','vertical'])
-#invL = 0.5
-#g = f.grad(x_init)
-#print (g)
-#u = x_init - invL*f.grad(x_init)
-        
-#%%
-# Run FISTA for least squares without regularization
-opt = {'tol': 1e-4, 'iter': 100}
-x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt)
-
-plt.imshow(x_fista0.subset(vertical=0).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,opt=opt)
-
-plt.imshow(x_fista0.subset(vertical=0).array)
-plt.title('FISTA1')
-plt.show()
-
-plt.semilogy(criter1)
-plt.show()
-
-# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
-x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
-
-plt.imshow(x_fbpd1.subset(vertical=0).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.subset(vertical=0).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, Cop, b, opt=opt)
-
-plt.imshow(x_CGLS.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
-
-plt.show()
-#%%
-#current = current + 1
-#a=fig.add_subplot(rows,cols,current)
-#a.set_title('FBPD TV')
-#imgplot = plt.imshow(x_fbpdtv.subset(vertical=0).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
deleted file mode 100755
index f77a678..0000000
--- a/Wrappers/Python/wip/simple_mc_demo.py
+++ /dev/null
@@ -1,142 +0,0 @@
-#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