diff options
Diffstat (limited to 'Wrappers/Python')
| -rwxr-xr-x | Wrappers/Python/wip/demo_astra_mc.py | 41 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_astra_nexus.py | 233 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/demo_astra_simple.py | 132 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/demo_astra_sophiabeads.py | 93 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_astra_sophiabeads3D.py | 59 | ||||
| -rw-r--r-- | Wrappers/Python/wip/work_out_adjoint.py | 115 | ||||
| -rw-r--r-- | Wrappers/Python/wip/work_out_adjoint3D.py | 131 | ||||
| -rw-r--r-- | Wrappers/Python/wip/work_out_adjoint_sophiabeads.py | 115 |
8 files changed, 325 insertions, 594 deletions
diff --git a/Wrappers/Python/wip/demo_astra_mc.py b/Wrappers/Python/wip/demo_astra_mc.py index e5999de..fde0120 100755 --- a/Wrappers/Python/wip/demo_astra_mc.py +++ b/Wrappers/Python/wip/demo_astra_mc.py @@ -5,9 +5,10 @@ # regularisation reconstructions. # Do all imports -from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry -from ccpi.optimisation.algs import FISTA -from ccpi.optimisation.funcs import Norm2sq, Norm1 +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 @@ -97,23 +98,28 @@ 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(numpy.zeros(x.shape),geometry=ig) +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) +f = Norm2Sq(Aop,b,c=0.5) # Run FISTA for least squares without regularization -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) +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 criteration +# Display reconstruction and criterion 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) + axarrf0[k].imshow(x_FISTA.as_array()[k],vmin=0,vmax=2.5) plt.show() -plt.semilogy(criter0) +plt.figure() +plt.semilogy(FISTA_alg.objective) plt.title('Criterion vs iterations, least squares') plt.show() @@ -121,17 +127,22 @@ plt.show() # such as 1-norm regularisation with choice of regularisation parameter lam. # Again the regulariser is over all channels: lam = 10 -g0 = Norm1(lam) +g1 = lam * L1Norm() -# Run FISTA for least squares plus 1-norm function. -x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) +# 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 criteration +# 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) + axarrf1[k].imshow(x_FISTA1.as_array()[k],vmin=0,vmax=2.5) plt.show() -plt.semilogy(criter1) +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 index 7892d24..64895dd 100644 --- a/Wrappers/Python/wip/demo_astra_nexus.py +++ b/Wrappers/Python/wip/demo_astra_nexus.py @@ -1,6 +1,6 @@ # This script demonstrates how to load a parallel beam data set in Nexus -# format, apply dark and flat field correction and reconstruct using the +# format and reconstruct using the # modular optimisation framework and the ASTRA Tomography toolbox. # # The data set is available from @@ -8,67 +8,46 @@ # and should be downloaded to a local directory to be specified below. # All own imports -from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry -from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, Norm1 -from ccpi.plugins.ops import CCPiProjectorSimple -from ccpi.processors import Normalizer, CenterOfRotationFinder -from ccpi.plugins.processors import AcquisitionDataPadder -from ccpi.io.reader import NexusReader +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 -import os - -# Define utility function to average over flat and dark images. -def avg_img(image): - shape = list(numpy.shape(image)) - l = shape.pop(0) - avg = numpy.zeros(shape) - for i in range(l): - avg += image[i] / l - return avg - -# Set up a reader object pointing to the Nexus data set. Revise path as needed. -reader = NexusReader(os.path.join(".." ,".." ,".." , "..", "CCPi-ReconstructionFramework","data" , "24737_fd.nxs" )) - -# Read and print the dimensions of the raw projections -dims = reader.get_projection_dimensions() -print (dims) - -# Load and average all flat and dark images in preparation for normalising data. -flat = avg_img(reader.load_flat()) -dark = avg_img(reader.load_dark()) - -# Set up normaliser object for normalising data by flat and dark images. -norm = Normalizer(flat_field=flat, dark_field=dark) - -# Load the raw projections and pass as input to the normaliser. -norm.set_input(reader.get_acquisition_data()) -# Set up CenterOfRotationFinder object to center data. -cor = CenterOfRotationFinder() +## 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="/media/newhd/shared/Data/nexus/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.set_input(norm.get_output()) +cor = CenterOfRotationFinder() +cor.set_input(data) center_of_rotation = cor.get_output() -# Set up AcquisitionDataPadder to pad data for centering using the computed -# center, set the output of the normaliser as input and execute to produce -# padded/centered data. -padder = AcquisitionDataPadder(center_of_rotation=center_of_rotation) -padder.set_input(norm.get_output()) -padded_data = padder.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_data2 = padded_data.subset(dimensions=['vertical','angle','horizontal']) -padded_data2.geometry = padded_data.geometry -padded_data2.geometry.angles = padded_data2.geometry.angles/180*numpy.pi +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_data2.geometry +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) @@ -78,94 +57,150 @@ print ("Define projector") Cop = AstraProjector3DSimple(ig, ag) # Test backprojection and projection -z1 = Cop.adjoint(padded_data2) +z1 = Cop.adjoint(padded_data) z2 = Cop.direct(z1) -plt.imshow(z1.subset(vertical=68).array) +plt.imshow(z1.subset(horizontal_x=68).array) plt.show() -# Set initial guess +# Set initial guess for reconstruction algorithms. print ("Initial guess") x_init = ImageData(geometry=ig) -# Create least squares object instance with projector and data. -print ("Create least squares object instance with projector and data.") -f = Norm2sq(Cop,padded_data2,c=0.5) - -# Run FISTA reconstruction for least squares without regularization -print ("Run FISTA for least squares") +# Set tolerance and number of iterations for reconstruction algorithms. opt = {'tol': 1e-4, 'iter': 100} -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) -plt.imshow(x_fista0.subset(horizontal_x=80).array) -plt.title('FISTA LS') -plt.colorbar() -plt.show() +# 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']) -# Set up 1-norm function for FISTA least squares plus 1-norm regularisation -print ("Run FISTA for least squares plus 1-norm regularisation") -lam = 0.1 -g0 = Norm1(lam) +x_CGLS = CGLS_alg.get_output() -# Run FISTA for least squares plus 1-norm function. -x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0,opt=opt) +# Fix color and slices for display +v1 = -0.01 +v2 = 0.13 +hx=80 +hy=80 +v=68 -plt.imshow(x_fista0.subset(horizontal_x=80).array) -plt.title('FISTA LS+1') -plt.colorbar() -plt.show() +# Display ortho slices of reconstruction +# Display all reconstructions and decay of objective function +cols = 3 +rows = 1 +fig = plt.figure() -# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm -print ("Run FBPD for least squares plus 1-norm regularisation") -x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) +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) -plt.imshow(x_fbpd1.subset(horizontal_x=80).array) -plt.title('FBPD LS+1') +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() -# Run CGLS, which should agree with the FISTA least squares -print ("Run CGLS for least squares") -x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, padded_data2, opt=opt) -plt.imshow(x_CGLS.subset(horizontal_x=80).array) -plt.title('CGLS') -plt.colorbar() +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 -cols = 4 -rows = 1 -current = 1 + fig = plt.figure() -current = current +current = 1 a=fig.add_subplot(rows,cols,current) -a.set_title('FISTA LS') -imgplot = plt.imshow(x_fista0.subset(horizontal_x=80).as_array()) +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('FISTA LS+1') -imgplot = plt.imshow(x_fista1.subset(horizontal_x=80).as_array()) +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('FBPD LS+1') -imgplot = plt.imshow(x_fbpd1.subset(horizontal_x=80).as_array()) +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('CGLS') -imgplot = plt.imshow(x_CGLS.subset(horizontal_x=80).as_array()) +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(criter0 , label='FISTA LS') -imgplot = plt.loglog(criter1 , label='FISTA LS+1') -imgplot = plt.loglog(criter_fbpd1, label='FBPD LS+1') -imgplot = plt.loglog(criter_CGLS, label='CGLS') +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 index 57caaee..925df77 100755 --- a/Wrappers/Python/wip/demo_astra_simple.py +++ b/Wrappers/Python/wip/demo_astra_simple.py @@ -2,12 +2,12 @@ # 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 and FBPD for 1-norm and TV regularisation. +# and 1-norm regularisation. # First make all imports from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry -from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D +from ccpi.optimisation.algorithms import FISTA, CGLS +from ccpi.optimisation.functions import Norm2Sq, L1Norm from ccpi.astra.operators import AstraProjectorSimple import numpy as np @@ -77,22 +77,30 @@ 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(np.zeros(x.shape),geometry=ig) -opt = {'tol': 1e-4, 'iter': 1000} +x_init = ImageData(geometry=ig) +opt = {'tol': 1e-4, 'iter': 100} # First a CGLS reconstruction can be done: -x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt) +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.semilogy(criter_CGLS) +plt.figure() +plt.semilogy(CGLS_alg.objective) plt.title('CGLS criterion') plt.show() @@ -102,72 +110,56 @@ plt.show() # 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 regularization -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None,opt) - -plt.imshow(x_fista0.array) -plt.title('FISTA Least squares') +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.semilogy(criter0) +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 = 0.1 -g0 = Norm1(lam) +lam = 1.0 +g0 = lam * L1Norm() # Run FISTA for least squares plus 1-norm function. -x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) - -plt.imshow(x_fista1.array) -plt.title('FISTA Least squares plus 1-norm regularisation') -plt.show() - -plt.semilogy(criter1) -plt.title('FISTA Least squares plus 1-norm regularisation criterion') -plt.show() - -# The least squares plus 1-norm regularisation problem can also be solved by -# other algorithms such as the Forward Backward Primal Dual algorithm. This -# algorithm minimises the sum of three functions and the least squares and -# 1-norm functions should be given as the second and third function inputs. -# In this test case, this algorithm requires more iterations to converge, so -# new options are specified. -opt_FBPD = {'tol': 1e-4, 'iter': 20000} -x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt_FBPD) - -plt.imshow(x_fbpd1.array) -plt.title('FBPD for least squares plus 1-norm regularisation') -plt.show() - -plt.semilogy(criter_fbpd1) -plt.title('FBPD for least squares plus 1-norm regularisation criterion') -plt.show() - -# The FBPD algorithm can also be used conveniently for TV regularisation: - -# Specify TV function object -lamtv = 10 -gtv = TV2D(lamtv) - -x_fbpdtv,it_fbpdtv,timing_fbpdtv,criter_fbpdtv=FBPD(x_init,None,f,gtv,opt_FBPD) - -plt.imshow(x_fbpdtv.array) +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.semilogy(criter_fbpdtv) +plt.figure() +plt.semilogy(FISTA_alg1.objective) +plt.title('FISTA LS+L1norm criterion') plt.show() # Compare all reconstruction and criteria clims = (0,1) -cols = 3 +cols = 2 rows = 2 current = 1 @@ -186,34 +178,20 @@ plt.axis('off') current = current + 1 a=fig.add_subplot(rows,cols,current) a.set_title('FISTA LS') -imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1]) +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') - -current = current + 1 -a=fig.add_subplot(rows,cols,current) -a.set_title('FBPD LS+1') -imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1]) -plt.axis('off') - -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]) +imgplot = plt.imshow(x_FISTA1.as_array(),vmin=clims[0],vmax=clims[1]) plt.axis('off') fig = plt.figure() -b=fig.add_subplot(1,1,1) -b.set_title('criteria') -imgplot = plt.loglog(criter_CGLS, label='CGLS') -imgplot = plt.loglog(criter0 , label='FISTA LS') -imgplot = plt.loglog(criter1 , label='FISTA LS+1') -imgplot = plt.loglog(criter_fbpd1, label='FBPD LS+1') -imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV') -b.legend(loc='lower left') +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 index b72a5f9..dce3d91 100755 --- a/Wrappers/Python/wip/demo_astra_sophiabeads.py +++ b/Wrappers/Python/wip/demo_astra_sophiabeads.py @@ -1,25 +1,27 @@ # This demo shows how to load a Nikon XTek micro-CT data set and reconstruct -# the central slice using the CGLS method. The SophiaBeads dataset with 256 -# projections is used as test data and can be obtained from here: +# 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 XTEKReader to load in the data. +# input to NikonDataReader to load in the data. # Do all imports -from ccpi.io.reader import XTEKReader 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.algs import CGLS +from ccpi.optimisation.algorithms import FISTA, CGLS +from ccpi.optimisation.functions import Norm2Sq, L1Norm -# Set up reader object and read the data -datareader = XTEKReader("REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_256_averaged.xtekct") -data = datareader.get_acquisition_data() +# 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()[:,:,1000]/60000.0) +sino = -np.log(data.as_array()[:,0,:]/60000.0) # Apply centering correction by zero padding, amount found manually cor_pad = 30 @@ -59,13 +61,74 @@ ig2d = ImageGeometry(voxel_num_x=N, Aop = AstraProjectorSimple(ig2d, ag2d,"gpu") # Set initial guess for CGLS reconstruction -x_init = ImageData(np.zeros((N,N)),geometry=ig2d) +x_init = ImageData(geometry=ig2d) -# Run 50-iteration CGLS reconstruction +# Set tolerance and number of iterations for reconstruction algorithms. opt = {'tol': 1e-4, 'iter': 50} -x, it, timing, criter = CGLS(x_init,Aop,data2d,opt=opt) -# Display reconstruction +# 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.imshow(x.as_array()) -plt.colorbar()
\ No newline at end of file +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 index 8c43657..624f121 100644 --- a/Wrappers/Python/wip/demo_astra_sophiabeads3D.py +++ b/Wrappers/Python/wip/demo_astra_sophiabeads3D.py @@ -1,38 +1,31 @@ # This demo shows how to load a Nikon XTek micro-CT data set and reconstruct -# the central slice using the CGLS method. The SophiaBeads dataset with 256 +# 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 XTEKReader to load in the data. +# input to NikonDataReader to load in the data. # Do all imports -from ccpi.io.reader import XTEKReader import numpy as np import matplotlib.pyplot as plt -from ccpi.framework import ImageGeometry, AcquisitionData, ImageData +from ccpi.io import NikonDataReader +from ccpi.framework import ImageGeometry, ImageData from ccpi.astra.operators import AstraProjector3DSimple -from ccpi.optimisation.algs import CGLS +from ccpi.optimisation.algorithms import CGLS -import numpy - -# Set up reader object and read the data -datareader = XTEKReader("REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_256_averaged.xtekct") -data = datareader.get_acquisition_data() - -# Crop data and fix dimension labels -data = AcquisitionData(data.array[:,:,901:1101], - geometry=data.geometry, - dimension_labels=['angle','horizontal','vertical']) -data.geometry.pixel_num_v = 200 +# 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.array = -np.log(data.as_array()/60000.0) +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]+cor_pad,data.shape[2])) -data_pad[:,cor_pad:,:] = data.array +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 @@ -58,7 +51,7 @@ ig = ImageGeometry(voxel_num_x=N, # 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*numpy.pi +data2.geometry.angles = -data2.geometry.angles/180*np.pi del data # Extract the Acquisition geometry for setting up projector. @@ -67,7 +60,7 @@ ag = data2.geometry # Set up the Projector (AcquisitionModel) using ASTRA on GPU Aop = AstraProjector3DSimple(ig, ag) -# So and show simple backprojection +# Do and show simple backprojection z = Aop.adjoint(data2) plt.figure() plt.imshow(z.subset(horizontal_x=1000).as_array()) @@ -82,17 +75,29 @@ plt.show() # Set initial guess for CGLS reconstruction x_init = ImageData(geometry=ig) -# Run 50-iteration CGLS reconstruction -opt = {'tol': 1e-4, 'iter': 20} -x, it, timing, criter = CGLS(x_init,Aop,data2,opt=opt) +# 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.subset(horizontal_x=1000).as_array()) +plt.imshow(x_CGLS.subset(horizontal_x=1000).as_array()) plt.show() plt.figure() -plt.imshow(x.subset(horizontal_y=1000).as_array()) +plt.imshow(x_CGLS.subset(horizontal_y=1000).as_array()) plt.show() plt.figure() -plt.imshow(x.subset(vertical=100).as_array()) +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 diff --git a/Wrappers/Python/wip/work_out_adjoint.py b/Wrappers/Python/wip/work_out_adjoint.py deleted file mode 100644 index 276fb76..0000000 --- a/Wrappers/Python/wip/work_out_adjoint.py +++ /dev/null @@ -1,115 +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 and FBPD for 1-norm and TV regularisation. - -# First make all imports -from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry, AcquisitionData -from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D -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 = 2 - -# 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 = 300 -x1 = -1 -x2 = 1 -dx = (x2-x1)/N -ig = ImageGeometry(voxel_num_x=N, - voxel_num_y=N, - voxel_size_x=dx, - voxel_size_y=dx) -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(3*N/8):round(5*N/8)] = 1 - -plt.imshow(x) -plt.title('Phantom image') -plt.colorbar() -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 = 40 -det_num = 400 - -SourceOrig = 20 -OrigDetec = 100 - -geo_mag = (SourceOrig+OrigDetec)/SourceOrig - -det_w = geo_mag*dx*1 - -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.colorbar() -plt.show() - -plt.imshow(z.array) -plt.title('Backprojected data') -plt.colorbar() -plt.show() - - -One = ImageData(geometry=ig) -xOne = One.as_array() -xOne[:,:] = 1.0 - -OneD = AcquisitionData(geometry=ag) -y1 = OneD.as_array() -y1[:,:] = 1.0 - -s1 = (OneD*(Aop.direct(One))).sum() -s2 = (One*(Aop.adjoint(OneD))).sum() -print(s1) -print(s2) -print(s2/s1) - -print((b*b).sum()) -print((z*Phantom).sum()) -print((z*Phantom).sum() / (b*b).sum()) - -#print(N/det_num) -#print(0.5*det_w/dx)
\ No newline at end of file diff --git a/Wrappers/Python/wip/work_out_adjoint3D.py b/Wrappers/Python/wip/work_out_adjoint3D.py deleted file mode 100644 index 0e4f847..0000000 --- a/Wrappers/Python/wip/work_out_adjoint3D.py +++ /dev/null @@ -1,131 +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 and FBPD for 1-norm and TV regularisation. - -# First make all imports -from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry, AcquisitionData -from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D -from ccpi.astra.operators import AstraProjector3DSimple - -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 = 2 - -# 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 = 100 -x1 = -1 -x2 = 1 -dx = (x2-x1)/N -ig = ImageGeometry(voxel_num_x=N, - voxel_num_y=N, - voxel_num_z=N, - voxel_size_x=dx, - voxel_size_y=dx, - voxel_size_z=dx) -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(3*N/8):round(5*N/8),:] = 1 - -plt.imshow(x[90,:,:]) -plt.title('Phantom image') -plt.colorbar() -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 = 200 -det_num = 100 - -det_w = dx - -if test_case==1: - angles = np.linspace(0,np.pi,angles_num,endpoint=False) - ag = AcquisitionGeometry('parallel', - '3D', - angles, - det_num, - det_w, - det_num, - det_w) -elif test_case==2: - angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) - - SourceOrig = 20 - OrigDetec = 100 - - geo_mag = (SourceOrig+OrigDetec)/SourceOrig - - det_w = geo_mag*dx - - ag = AcquisitionGeometry('cone', - '3D', - angles, - det_num, - det_w, - 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 = AstraProjector3DSimple(ig, ag) - -# 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) - -print((b*b).sum()) -print((z*Phantom).sum()) -print((z*Phantom).sum() / (b*b).sum()) - -plt.imshow(b.array[:,0,:]) -plt.title('Simulated data') -plt.colorbar() -plt.show() -plt.imshow(b.array[:,1,:]) -plt.title('Simulated data') -plt.colorbar() -plt.show() - -plt.imshow(z.array[50,:,:]) -plt.title('Backprojected data') -plt.colorbar() -plt.show() - - -One = ImageData(geometry=ig) -xOne = One.as_array() -xOne[:,:,:] = 1.0 - -OneD = b.clone() -y1 = OneD.as_array() -y1[:,:,:] = 1.0 - -s1 = (OneD*(Aop.direct(One))).sum() -s2 = (One*(Aop.adjoint(OneD))).sum() -print(s1) -print(s2) -print(s2/s1) - - - -#print(N/det_num) -#print(0.5*det_w/dx)
\ No newline at end of file diff --git a/Wrappers/Python/wip/work_out_adjoint_sophiabeads.py b/Wrappers/Python/wip/work_out_adjoint_sophiabeads.py deleted file mode 100644 index 7f35a16..0000000 --- a/Wrappers/Python/wip/work_out_adjoint_sophiabeads.py +++ /dev/null @@ -1,115 +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 and FBPD for 1-norm and TV regularisation. - -# First make all imports -from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry, AcquisitionData -from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D -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 = 2 - -# 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 = 2000 -x1 = -16.015690364093633 -x2 = 16.015690364093633 -dx = (x2-x1)/N -ig = ImageGeometry(voxel_num_x=N, - voxel_num_y=N, - voxel_size_x=dx, - voxel_size_y=dx) -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(3*N/8):round(5*N/8)] = 1 - -plt.imshow(x) -plt.title('Phantom image') -plt.colorbar() -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 = 4 -det_num = 2500 - -SourceOrig = 80.6392412185669 -OrigDetec = 926.3637587814331 - -geo_mag = (SourceOrig+OrigDetec)/SourceOrig - -det_w = geo_mag*dx*1 - -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.colorbar() -plt.show() - -plt.imshow(z.array) -plt.title('Backprojected data') -plt.colorbar() -plt.show() - - -One = ImageData(geometry=ig) -xOne = One.as_array() -xOne[:,:] = 1.0 - -OneD = AcquisitionData(geometry=ag) -y1 = OneD.as_array() -y1[:,:] = 1.0 - -s1 = (OneD*(Aop.direct(One))).sum() -s2 = (One*(Aop.adjoint(OneD))).sum() -print(s1) -print(s2) -print(s2/s1) - -print((b*b).sum()) -print((z*Phantom).sum()) -print((z*Phantom).sum() / (b*b).sum()) - -#print(N/det_num) -#print(0.5*det_w/dx)
\ No newline at end of file |