Merge pull request #29 from vais-ral/Remove_demos

Deleted demos folder (wip)

5 files changed, 0 insertions, 789 deletions

diff --git a/Wrappers/Python/wip/demo_astra_mc.py b/Wrappers/Python/wip/demo_astra_mc.py
deleted file mode 100755
index fde0120..0000000
--- a/Wrappers/Python/wip/demo_astra_mc.py
+++ /dev/null
@@ -1,148 +0,0 @@
-
-# This demo demonstrates a simple multichannel reconstruction case. A 
-# synthetic 3-channel phantom image is set up, data is simulated and the FISTA 
-# algorithm is used to compute least squares and least squares with 1-norm 
-# regularisation reconstructions.
-
-# Do all imports
-from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, \
-                           AcquisitionGeometry
-from ccpi.optimisation.algorithms import FISTA
-from ccpi.optimisation.functions import Norm2Sq, L1Norm
-from ccpi.astra.operators import AstraProjectorMC
-
-import numpy
-import matplotlib.pyplot as plt
-
-# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
-test_case = 1
-
-# Set up phantom NxN pixels and 3 channels. Set up the ImageGeometry and fill 
-# some test data in each of the channels. Display each channel as image.
-N = 128
-numchannels = 3
-
-ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels)
-Phantom = ImageData(geometry=ig)
-
-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
-
-f, axarr = plt.subplots(1,numchannels)
-for k in numpy.arange(3):
-    axarr[k].imshow(x[k],vmin=0,vmax=2.5)
-plt.show()
-
-# Set up AcquisitionGeometry object to hold the parameters of the measurement
-# setup geometry: # Number of angles, the actual angles from 0 to 
-# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector 
-# pixel relative to an object pixel, the number of detector pixels, and the 
-# source-origin and origin-detector distance (here the origin-detector distance 
-# set to 0 to simulate a "virtual detector" with same detector pixel size as
-# object pixel size).
-angles_num = 20
-det_w = 1.0
-det_num = N
-SourceOrig = 200
-OrigDetec = 0
-
-if test_case==1:
-    angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False)
-    ag = AcquisitionGeometry('parallel',
-                             '2D',
-                             angles,
-                             det_num,
-                             det_w,
-                             channels=numchannels)
-elif test_case==2:
-    angles = numpy.linspace(0,2*numpy.pi,angles_num,endpoint=False)
-    ag = AcquisitionGeometry('cone',
-                            '2D',
-                            angles,
-                            det_num,
-                            det_w,
-                            dist_source_center=SourceOrig, 
-                            dist_center_detector=OrigDetec,
-                            channels=numchannels)
-else:
-    NotImplemented
-
-# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
-# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
-Aop = AstraProjectorMC(ig, ag, 'gpu')
-
-# Forward and backprojection are available as methods direct and adjoint. Here 
-# generate test data b and do simple backprojection to obtain z. Applies 
-# channel by channel
-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()
-
-z = Aop.adjoint(b)
-
-fo, axarro = plt.subplots(1,numchannels)
-for k in range(3):
-    axarro[k].imshow(z.as_array()[k],vmin=0,vmax=3500)
-plt.show()
-
-# Using the test data b, different reconstruction methods can now be set up as
-# demonstrated in the rest of this file. In general all methods need an initial 
-# guess and some algorithm options to be set:
-x_init = ig.allocate(0.0)
-opt = {'tol': 1e-4, 'iter': 200}
-
-# Create least squares object instance with projector, test data and a constant 
-# coefficient of 0.5. Note it is least squares over all channels:
-f = Norm2Sq(Aop,b,c=0.5)
-
-# Run FISTA for least squares without regularization
-FISTA_alg = FISTA()
-FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
-FISTA_alg.max_iteration = 2000
-FISTA_alg.run(opt['iter'])
-x_FISTA = FISTA_alg.get_output()
-
-# Display reconstruction and criterion
-ff0, axarrf0 = plt.subplots(1,numchannels)
-for k in numpy.arange(3):
-    axarrf0[k].imshow(x_FISTA.as_array()[k],vmin=0,vmax=2.5)
-plt.show()
-
-plt.figure()
-plt.semilogy(FISTA_alg.objective)
-plt.title('Criterion vs iterations, least squares')
-plt.show()
-
-# FISTA can also solve regularised forms by specifying a second function object
-# such as 1-norm regularisation with choice of regularisation parameter lam. 
-# Again the regulariser is over all channels:
-lam = 10
-g1 = lam * L1Norm()
-
-# Run FISTA for least squares plus 1-norm regularisation.
-FISTA_alg1 = FISTA()
-FISTA_alg1.set_up(x_init=x_init, f=f, g=g1, opt=opt)
-FISTA_alg1.max_iteration = 2000
-FISTA_alg1.run(opt['iter'])
-x_FISTA1 = FISTA_alg1.get_output()
-
-# Display reconstruction and criterion
-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.figure()
-plt.semilogy(FISTA_alg1.objective)
-plt.title('Criterion vs iterations, least squares plus 1-norm regu')
-plt.show()
diff --git a/Wrappers/Python/wip/demo_astra_nexus.py b/Wrappers/Python/wip/demo_astra_nexus.py
deleted file mode 100644
index bb23e81..0000000
--- a/Wrappers/Python/wip/demo_astra_nexus.py
+++ /dev/null
@@ -1,207 +0,0 @@
-
-# This script demonstrates how to load a parallel beam data set in Nexus 
-# format and reconstruct using the
-# modular optimisation framework and the ASTRA Tomography toolbox.
-# 
-# The data set is available from
-# https://github.com/DiamondLightSource/Savu/blob/master/test_data/data/24737_fd.nxs
-# and should be downloaded to a local directory to be specified below.
-# The data should be flat-field and dark-field corrected.
-
-# All own imports
-from ccpi.framework import ImageData, ImageGeometry
-from ccpi.optimisation.algorithms import CGLS, FISTA
-from ccpi.optimisation.functions import Norm2Sq, L1Norm
-from ccpi.processors import CenterOfRotationFinder 
-from ccpi.io import NEXUSDataReader
-from ccpi.astra.operators import AstraProjector3DSimple
-
-# All external imports
-import numpy
-import matplotlib.pyplot as plt
-
-## Set up a reader object pointing to the Nexus data set. Revise path as needed.
-# The data is already  corrected for by flat and dark field.
-myreader = NEXUSDataReader(nexus_file="REPLACE_THIS_BY_PATH_TO_DATASET/24737_fd_normalised.nxs" )
-data= myreader.load_data()
-
-# Negative logarithm transoform.
-data.fill( -numpy.log(data.as_array() ))
-
-# Set up CenterOfRotationFinder object to center data.
-# Set the output of the normaliser as the input and execute to determine center.
-cor = CenterOfRotationFinder()
-cor.set_input(data)
-center_of_rotation = cor.get_output()
-
-# From computed center, determine amount of zero-padding to apply, apply
-# and update geometry to wider detector.
-cor_pad = int(2*(center_of_rotation - data.shape[2]/2))
-data_pad = numpy.zeros((data.shape[0],data.shape[1],data.shape[2]+cor_pad))
-data_pad[:,:,:-cor_pad] = data.as_array()
-data.geometry.pixel_num_h = data.geometry.pixel_num_h + cor_pad
-data.array = data_pad
-
-# Permute array and convert angles to radions for ASTRA
-padded_data = data.subset(dimensions=['vertical','angle','horizontal'])
-padded_data.geometry = data.geometry
-padded_data.geometry.angles = padded_data.geometry.angles/180*numpy.pi
-
-# Create Acquisition and Image Geometries for setting up projector.
-ag = padded_data.geometry
-ig = ImageGeometry(voxel_num_x=ag.pixel_num_h,
-                   voxel_num_y=ag.pixel_num_h, 
-                   voxel_num_z=ag.pixel_num_v)
-
-# Define the projector object
-print ("Define projector")
-Cop = AstraProjector3DSimple(ig, ag)
-
-# Test backprojection and projection
-z1 = Cop.adjoint(padded_data)
-z2 = Cop.direct(z1)
-
-plt.imshow(z1.subset(horizontal_x=68).array)
-plt.show()
-
-# Set initial guess for reconstruction algorithms.
-print ("Initial guess")
-x_init = ImageData(geometry=ig)
-
-# Set tolerance and number of iterations for reconstruction algorithms.
-opt = {'tol': 1e-4, 'iter': 100}
-
-# First a CGLS reconstruction can be done:
-CGLS_alg = CGLS()
-CGLS_alg.set_up(x_init, Cop, padded_data)
-CGLS_alg.max_iteration = 2000
-CGLS_alg.run(opt['iter'])
-
-x_CGLS = CGLS_alg.get_output()
-
-# Fix color and slices for display
-v1 = -0.01
-v2 = 0.13
-hx=80
-hy=80
-v=68
-
-# Display ortho slices of reconstruction
-# Display all reconstructions and decay of objective function
-cols = 3
-rows = 1
-fig = plt.figure()
-
-current = 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('horizontal_x')
-imgplot = plt.imshow(x_CGLS.subset(horizontal_x=hx).as_array(),vmin=v1,vmax=v2)
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('horizontal_y')
-imgplot = plt.imshow(x_CGLS.subset(horizontal_y=hy).as_array(),vmin=v1,vmax=v2)
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('vertical')
-imgplot = plt.imshow(x_CGLS.subset(vertical=v).as_array(),vmin=v1,vmax=v2)
-plt.colorbar()
-
-plt.suptitle('CGLS reconstruction slices')
-plt.show()
-
-plt.figure()
-plt.semilogy(CGLS_alg.objective)
-plt.title('CGLS criterion')
-plt.show()
-
-# Create least squares object instance with projector and data.
-print ("Create least squares object instance with projector and data.")
-f = Norm2Sq(Cop,padded_data,c=0.5)
-
-
-# Run FISTA for least squares without constraints
-FISTA_alg = FISTA()
-FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
-FISTA_alg.max_iteration = 2000
-FISTA_alg.run(opt['iter'])
-x_FISTA = FISTA_alg.get_output()
-
-# Display ortho slices of reconstruction
-# Display all reconstructions and decay of objective function
-
-fig = plt.figure()
-
-current = 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('horizontal_x')
-imgplot = plt.imshow(x_FISTA.subset(horizontal_x=hx).as_array(),vmin=v1,vmax=v2)
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('horizontal_y')
-imgplot = plt.imshow(x_FISTA.subset(horizontal_y=hy).as_array(),vmin=v1,vmax=v2)
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('vertical')
-imgplot = plt.imshow(x_FISTA.subset(vertical=v).as_array(),vmin=v1,vmax=v2)
-plt.colorbar()
-
-plt.suptitle('FISTA Least squares reconstruction slices')
-plt.show()
-
-plt.figure()
-plt.semilogy(FISTA_alg.objective)
-plt.title('FISTA Least squares criterion')
-plt.show()
-
-# Create 1-norm function object
-lam = 30.0
-g0 = lam * L1Norm()
-
-# Run FISTA for least squares plus 1-norm function.
-FISTA_alg1 = FISTA()
-FISTA_alg1.set_up(x_init=x_init, f=f, g=g0, opt=opt)
-FISTA_alg1.max_iteration = 2000
-FISTA_alg1.run(opt['iter'])
-x_FISTA1 = FISTA_alg1.get_output()
-
-# Display all reconstructions and decay of objective function
-fig = plt.figure()
-
-current = 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('horizontal_x')
-imgplot = plt.imshow(x_FISTA1.subset(horizontal_x=hx).as_array(),vmin=v1,vmax=v2)
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('horizontal_y')
-imgplot = plt.imshow(x_FISTA1.subset(horizontal_y=hy).as_array(),vmin=v1,vmax=v2)
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('vertical')
-imgplot = plt.imshow(x_FISTA1.subset(vertical=v).as_array(),vmin=v1,vmax=v2)
-plt.colorbar()
-
-plt.suptitle('FISTA LS+1 reconstruction slices')
-plt.show()
-
-
-plt.figure()
-plt.semilogy(FISTA_alg1.objective)
-plt.title('FISTA LS+1 criterion')
-plt.show()
-
-
-fig = plt.figure()
-b=fig.add_subplot(1,1,1)
-b.set_title('criteria')
-imgplot = plt.loglog(CGLS_alg.objective , label='CGLS')
-imgplot = plt.loglog(FISTA_alg.objective , label='FISTA LS')
-imgplot = plt.loglog(FISTA_alg1.objective, label='FISTA LS+1')
-b.legend(loc='right')
-plt.show()
diff --git a/Wrappers/Python/wip/demo_astra_simple.py b/Wrappers/Python/wip/demo_astra_simple.py
deleted file mode 100755
index 925df77..0000000
--- a/Wrappers/Python/wip/demo_astra_simple.py
+++ /dev/null
@@ -1,197 +0,0 @@
-
-# This demo illustrates how ASTRA 2D projectors can be used with
-# the modular optimisation framework. The demo sets up a 2D test case and 
-# demonstrates reconstruction using CGLS, as well as FISTA for least squares 
-# and 1-norm regularisation.
-
-# First make all imports
-from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry
-from ccpi.optimisation.algorithms import FISTA, CGLS
-from ccpi.optimisation.functions import Norm2Sq, L1Norm
-from ccpi.astra.operators import AstraProjectorSimple
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
-test_case = 1
-
-# Set up phantom size NxN by creating ImageGeometry, initialising the 
-# ImageData object with this geometry and empty array and finally put some
-# data into its array, and display as image.
-N = 128
-ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
-Phantom = ImageData(geometry=ig)
-
-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.title('Phantom image')
-plt.show()
-
-# Set up AcquisitionGeometry object to hold the parameters of the measurement
-# setup geometry: # Number of angles, the actual angles from 0 to 
-# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector 
-# pixel relative to an object pixel, the number of detector pixels, and the 
-# source-origin and origin-detector distance (here the origin-detector distance 
-# set to 0 to simulate a "virtual detector" with same detector pixel size as
-# object pixel size).
-angles_num = 20
-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)
-    ag = AcquisitionGeometry('parallel',
-                             '2D',
-                             angles,
-                             det_num,det_w)
-elif test_case==2:
-    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
-    ag = AcquisitionGeometry('cone',
-                             '2D',
-                             angles,
-                             det_num,
-                             det_w,
-                             dist_source_center=SourceOrig, 
-                             dist_center_detector=OrigDetec)
-else:
-    NotImplemented
-
-# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
-# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
-Aop = AstraProjectorSimple(ig, ag, 'gpu')
-
-# Forward and backprojection are available as methods direct and adjoint. Here 
-# generate test data b and do simple backprojection to obtain z.
-b = Aop.direct(Phantom)
-z = Aop.adjoint(b)
-
-plt.imshow(b.array)
-plt.title('Simulated data')
-plt.show()
-
-plt.imshow(z.array)
-plt.title('Backprojected data')
-plt.colorbar()
-plt.show()
-
-# Using the test data b, different reconstruction methods can now be set up as
-# demonstrated in the rest of this file. In general all methods need an initial 
-# guess and some algorithm options to be set:
-x_init = ImageData(geometry=ig)
-opt = {'tol': 1e-4, 'iter': 100}
-
-# First a CGLS reconstruction can be done:
-CGLS_alg = CGLS()
-CGLS_alg.set_up(x_init, Aop, b )
-CGLS_alg.max_iteration = 2000
-CGLS_alg.run(opt['iter'])
-
-x_CGLS = CGLS_alg.get_output()
-
-plt.figure()
-plt.imshow(x_CGLS.array)
-plt.title('CGLS')
-plt.show()
-
-plt.figure()
-plt.semilogy(CGLS_alg.objective)
-plt.title('CGLS criterion')
-plt.show()
-
-# CGLS solves the simple least-squares problem. The same problem can be solved 
-# by FISTA by setting up explicitly a least squares function object and using 
-# no regularisation:
-
-# Create least squares object instance with projector, test data and a constant 
-# coefficient of 0.5:
-f = Norm2Sq(Aop,b,c=0.5)
-
-# Run FISTA for least squares without constraints
-FISTA_alg = FISTA()
-FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
-FISTA_alg.max_iteration = 2000
-FISTA_alg.run(opt['iter'])
-x_FISTA = FISTA_alg.get_output()
-
-plt.figure()
-plt.imshow(x_FISTA.array)
-plt.title('FISTA Least squares reconstruction')
-plt.colorbar()
-plt.show()
-
-plt.figure()
-plt.semilogy(FISTA_alg.objective)
-plt.title('FISTA Least squares criterion')
-plt.show()
-
-
-# FISTA can also solve regularised forms by specifying a second function object
-# such as 1-norm regularisation with choice of regularisation parameter lam:
-
-# Create 1-norm function object
-lam = 1.0
-g0 = lam * L1Norm()
-
-# Run FISTA for least squares plus 1-norm function.
-FISTA_alg1 = FISTA()
-FISTA_alg1.set_up(x_init=x_init, f=f, g=g0, opt=opt)
-FISTA_alg1.max_iteration = 2000
-FISTA_alg1.run(opt['iter'])
-x_FISTA1 = FISTA_alg1.get_output()
-
-plt.figure()
-plt.imshow(x_FISTA1.array)
-plt.title('FISTA LS+L1Norm reconstruction')
-plt.colorbar()
-plt.show()
-
-plt.figure()
-plt.semilogy(FISTA_alg1.objective)
-plt.title('FISTA LS+L1norm criterion')
-plt.show()
-
-
-# Compare all reconstruction and criteria
-clims = (0,1)
-cols = 2
-rows = 2
-current = 1
-
-fig = plt.figure()
-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])
-plt.axis('off')
-
-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])
-plt.axis('off')
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA LS')
-imgplot = plt.imshow(x_FISTA.as_array(),vmin=clims[0],vmax=clims[1])
-plt.axis('off')
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA LS+1')
-imgplot = plt.imshow(x_FISTA1.as_array(),vmin=clims[0],vmax=clims[1])
-plt.axis('off')
-
-fig = plt.figure()
-a=fig.add_subplot(1,1,1)
-a.set_title('criteria')
-imgplot = plt.loglog(CGLS_alg.objective, label='CGLS')
-imgplot = plt.loglog(FISTA_alg.objective , label='FISTA LS')
-imgplot = plt.loglog(FISTA_alg1.objective , label='FISTA LS+1')
-a.legend(loc='lower left')
-plt.show()
diff --git a/Wrappers/Python/wip/demo_astra_sophiabeads.py b/Wrappers/Python/wip/demo_astra_sophiabeads.py
deleted file mode 100755
index dce3d91..0000000
--- a/Wrappers/Python/wip/demo_astra_sophiabeads.py
+++ /dev/null
@@ -1,134 +0,0 @@
-
-# This demo shows how to load a Nikon XTek micro-CT data set and reconstruct
-# the central slice using the CGLS and FISTA methods. The SophiaBeads dataset 
-# with 64 projections is used as test data and can be obtained from here:
-# https://zenodo.org/record/16474
-# The filename with full path to the .xtekct file should be given as string 
-# input to NikonDataReader to  load in the data.
-
-# Do all imports
-import numpy as np
-import matplotlib.pyplot as plt
-from ccpi.io import NikonDataReader
-from ccpi.framework import ImageGeometry, AcquisitionGeometry, AcquisitionData, ImageData
-from ccpi.astra.operators import AstraProjectorSimple
-from ccpi.optimisation.algorithms import FISTA, CGLS
-from ccpi.optimisation.functions import Norm2Sq, L1Norm
-
-# Set up reader object and read in central slice the data
-datareader = NikonDataReader(xtek_file="REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_64_averaged.xtekct",
-                             roi=[(1000,1001),(0,2000)])
-data = datareader.load_projections()
-
-# Extract central slice, scale and negative-log transform
-sino = -np.log(data.as_array()[:,0,:]/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
-
-# 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(geometry=ig2d)
-
-# Set tolerance and number of iterations for reconstruction algorithms.
-opt = {'tol': 1e-4, 'iter': 50}
-
-# First a CGLS reconstruction can be done:
-CGLS_alg = CGLS()
-CGLS_alg.set_up(x_init, Aop, data2d)
-CGLS_alg.max_iteration = 2000
-CGLS_alg.run(opt['iter'])
-
-x_CGLS = CGLS_alg.get_output()
-
-plt.figure()
-plt.imshow(x_CGLS.array)
-plt.title('CGLS')
-plt.colorbar()
-plt.show()
-
-plt.figure()
-plt.semilogy(CGLS_alg.objective)
-plt.title('CGLS criterion')
-plt.show()
-
-# CGLS solves the simple least-squares problem. The same problem can be solved 
-# by FISTA by setting up explicitly a least squares function object and using 
-# no regularisation:
-
-# Create least squares object instance with projector, test data and a constant 
-# coefficient of 0.5:
-f = Norm2Sq(Aop,data2d)
-
-# Run FISTA for least squares without constraints
-FISTA_alg = FISTA()
-FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
-FISTA_alg.max_iteration = 2000
-FISTA_alg.run(opt['iter'])
-x_FISTA = FISTA_alg.get_output()
-
-plt.figure()
-plt.imshow(x_FISTA.array)
-plt.title('FISTA Least squares reconstruction')
-plt.colorbar()
-plt.show()
-
-plt.figure()
-plt.semilogy(FISTA_alg.objective)
-plt.title('FISTA Least squares criterion')
-plt.show()
-
-# Create 1-norm function object
-lam = 1.0
-g0 = lam * L1Norm()
-
-# Run FISTA for least squares plus 1-norm function.
-FISTA_alg1 = FISTA()
-FISTA_alg1.set_up(x_init=x_init, f=f, g=g0, opt=opt)
-FISTA_alg1.max_iteration = 2000
-FISTA_alg1.run(opt['iter'])
-x_FISTA1 = FISTA_alg1.get_output()
-
-plt.figure()
-plt.imshow(x_FISTA1.array)
-plt.title('FISTA LS+L1Norm reconstruction')
-plt.colorbar()
-plt.show()
-
-plt.figure()
-plt.semilogy(FISTA_alg1.objective)
-plt.title('FISTA LS+L1norm criterion')
-plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/demo_astra_sophiabeads3D.py b/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
deleted file mode 100644
index 624f121..0000000
--- a/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
+++ /dev/null
@@ -1,103 +0,0 @@
-
-# This demo shows how to load a Nikon XTek micro-CT data set and reconstruct
-# a volume of 200 slices using the CGLS method. The SophiaBeads dataset with 64 
-# projections is used as test data and can be obtained from here:
-# https://zenodo.org/record/16474
-# The filename with full path to the .xtekct file should be given as string 
-# input to NikonDataReader to  load in the data.
-
-# Do all imports
-import numpy as np
-import matplotlib.pyplot as plt
-from ccpi.io import NikonDataReader
-from ccpi.framework import ImageGeometry, ImageData
-from ccpi.astra.operators import AstraProjector3DSimple
-from ccpi.optimisation.algorithms import CGLS
-
-# Set up reader object and read in 200 central slices of the data
-datareader= NikonDataReader(xtek_file="REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_64_averaged.xtekct",
-                            roi=[(901,1101),(0,2000)])
-data = datareader.load_projections()
-
-# Scale and negative-log transform
-data.fill(-np.log(data.as_array()/60000.0))
-
-# Apply centering correction by zero padding, amount found manually
-cor_pad = 30
-data_pad = np.zeros((data.shape[0],data.shape[1],data.shape[2]+cor_pad))
-data_pad[:,:,cor_pad:] = data.array
-data.geometry.pixel_num_h = data.geometry.pixel_num_h + cor_pad
-data.array = data_pad
-
-# 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
-ig = ImageGeometry(voxel_num_x=N,
-                   voxel_num_y=N,
-                   voxel_num_z=data.geometry.pixel_num_v,
-                   voxel_size_x=voxel_size_h, 
-                   voxel_size_y=voxel_size_h,
-                   voxel_size_z=voxel_size_h)
-
-# Permute array and convert angles to radions for ASTRA; delete old data.
-data2 = data.subset(dimensions=['vertical','angle','horizontal'])
-data2.geometry = data.geometry
-data2.geometry.angles = -data2.geometry.angles/180*np.pi
-del data
-
-# Extract the Acquisition geometry for setting up projector.
-ag = data2.geometry
-
-# Set up the Projector (AcquisitionModel) using ASTRA on GPU
-Aop = AstraProjector3DSimple(ig, ag)
-
-# Do and show simple backprojection
-z = Aop.adjoint(data2)
-plt.figure()
-plt.imshow(z.subset(horizontal_x=1000).as_array())
-plt.show()
-plt.figure()
-plt.imshow(z.subset(horizontal_y=1000).as_array())
-plt.show()
-plt.figure()
-plt.imshow(z.subset(vertical=100).array)
-plt.show()
-
-# Set initial guess for CGLS reconstruction
-x_init = ImageData(geometry=ig)
-
-# Set tolerance and number of iterations for reconstruction algorithms.
-opt = {'tol': 1e-4, 'iter': 30}
-
-# Run a CGLS reconstruction can be done:
-CGLS_alg = CGLS()
-CGLS_alg.set_up(x_init, Aop, data2)
-CGLS_alg.max_iteration = 2000
-CGLS_alg.run(opt['iter'])
-
-x_CGLS = CGLS_alg.get_output()
-
-# Display ortho slices of reconstruction
-plt.figure()
-plt.imshow(x_CGLS.subset(horizontal_x=1000).as_array())
-plt.show()
-plt.figure()
-plt.imshow(x_CGLS.subset(horizontal_y=1000).as_array())
-plt.show()
-plt.figure()
-plt.imshow(x_CGLS.subset(vertical=100).as_array())
-plt.show()
-
-plt.figure()
-plt.semilogy(CGLS_alg.objective)
-plt.title('CGLS criterion')
-plt.show()
\ No newline at end of file