diff options
| -rwxr-xr-x | Wrappers/Python/wip/demo_astra_mc.py | 148 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_astra_nexus.py | 207 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/demo_astra_simple.py | 197 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/demo_astra_sophiabeads.py | 134 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_astra_sophiabeads3D.py | 103 |
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 |