From b2afe5f5d740b5c79556b32ef20a419dbb4366c7 Mon Sep 17 00:00:00 2001 From: Jakob Jorgensen Date: Tue, 24 Apr 2018 11:22:18 +0100 Subject: Rename consistent --- Wrappers/Python/wip/demo_astra_mc.py | 137 ++++++++++++++++ Wrappers/Python/wip/demo_astra_simple.py | 219 ++++++++++++++++++++++++++ Wrappers/Python/wip/demo_astra_sophiabeads.py | 73 +++++++++ Wrappers/Python/wip/demo_sophiabeads.py | 73 --------- Wrappers/Python/wip/simple_demo_astra.py | 219 -------------------------- Wrappers/Python/wip/simple_mc_demo.py | 137 ---------------- 6 files changed, 429 insertions(+), 429 deletions(-) create mode 100755 Wrappers/Python/wip/demo_astra_mc.py create mode 100755 Wrappers/Python/wip/demo_astra_simple.py create mode 100755 Wrappers/Python/wip/demo_astra_sophiabeads.py delete mode 100755 Wrappers/Python/wip/demo_sophiabeads.py delete mode 100755 Wrappers/Python/wip/simple_demo_astra.py delete mode 100755 Wrappers/Python/wip/simple_mc_demo.py (limited to 'Wrappers') diff --git a/Wrappers/Python/wip/demo_astra_mc.py b/Wrappers/Python/wip/demo_astra_mc.py new file mode 100755 index 0000000..3b78eb3 --- /dev/null +++ b/Wrappers/Python/wip/demo_astra_mc.py @@ -0,0 +1,137 @@ + +# 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.algs import FISTA +from ccpi.optimisation.funcs import Norm2sq, Norm1 +from ccpi.astra.ops 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 = ImageData(np.zeros(x.shape),geometry=ig) +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 +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) + +# Display reconstruction and criteration +ff0, axarrf0 = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrf0[k].imshow(x_fista0.as_array()[k],vmin=0,vmax=2.5) +plt.show() + +plt.semilogy(criter0) +plt.title('Criterion vs iterations, least squares') +plt.show() + +# 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 = 0.1 +g0 = Norm1(lam) + +# Run FISTA for least squares plus 1-norm function. +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) + +# Display reconstruction and criteration +ff1, axarrf1 = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrf1[k].imshow(x_fista1.as_array()[k],vmin=0,vmax=2.5) +plt.show() + +plt.semilogy(criter1) +plt.title('Criterion vs iterations, least squares plus 1-norm regu') +plt.show() diff --git a/Wrappers/Python/wip/demo_astra_simple.py b/Wrappers/Python/wip/demo_astra_simple.py new file mode 100755 index 0000000..8529c98 --- /dev/null +++ b/Wrappers/Python/wip/demo_astra_simple.py @@ -0,0 +1,219 @@ + +# 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 +from ccpi.optimisation.algs import FISTA, FBPD, CGLS +from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D +from ccpi.astra.ops 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.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} + +# First a CGLS reconstruction can be done: +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt) + +plt.imshow(x_CGLS.array) +plt.title('CGLS') +plt.show() + +plt.semilogy(criter_CGLS) +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 regularization +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None,opt) + +plt.imshow(x_fista0.array) +plt.title('FISTA Least squares') +plt.show() + +plt.semilogy(criter0) +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) + +# 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': 7000} +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 = 1 +gtv = TV2D(lamtv) + +x_fbpdtv,it_fbpdtv,timing_fbpdtv,criter_fbpdtv=FBPD(x_init,None,f,gtv,opt_FBPD) + +plt.imshow(x_fbpdtv.array) +plt.show() + +plt.semilogy(criter_fbpdtv) +plt.show() + + +# Compare all reconstruction and criteria +clims = (0,1) +cols = 3 +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_fista0.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]) +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') +plt.show() diff --git a/Wrappers/Python/wip/demo_astra_sophiabeads.py b/Wrappers/Python/wip/demo_astra_sophiabeads.py new file mode 100755 index 0000000..4a0f963 --- /dev/null +++ b/Wrappers/Python/wip/demo_astra_sophiabeads.py @@ -0,0 +1,73 @@ + +# 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: +# 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. + +# Do all imports +from ccpi.io.reader import XTEKReader +import numpy as np +import matplotlib.pyplot as plt +from ccpi.framework import ImageGeometry, AcquisitionGeometry, AcquisitionData, ImageData +from ccpi.astra.ops import AstraProjectorSimple +from ccpi.optimisation.algs import CGLS + +# 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() + +# Extract central slice, scale and negative-log transform +sino = -np.log(data.as_array()[:,:,1000]/60000.0) + +# Apply centering correction by zero padding, amount found manually +cor_pad = 30 +sino_pad = np.zeros((sino.shape[0],sino.shape[1]+cor_pad)) +sino_pad[:,cor_pad:] = sino + +# Simple beam hardening correction as done in SophiaBeads coda +sino_pad = sino_pad**2 + +# Extract AcquisitionGeometry for central slice for 2D fanbeam reconstruction +ag2d = AcquisitionGeometry('cone', + '2D', + angles=-np.pi/180*data.geometry.angles, + pixel_num_h=data.geometry.pixel_num_h + cor_pad, + pixel_size_h=data.geometry.pixel_size_h, + dist_source_center=data.geometry.dist_source_center, + dist_center_detector=data.geometry.dist_center_detector) + +# Set up AcquisitionData object for central slice 2D fanbeam +data2d = AcquisitionData(sino_pad,geometry=ag2d) + +# Choose the number of voxels to reconstruct onto as number of detector pixels +N = data.geometry.pixel_num_h + +# Geometric magnification +mag = (np.abs(data.geometry.dist_center_detector) + \ + np.abs(data.geometry.dist_source_center)) / \ + np.abs(data.geometry.dist_source_center) + +# Voxel size is detector pixel size divided by mag +voxel_size_h = data.geometry.pixel_size_h / mag + +# Construct the appropriate ImageGeometry +ig2d = ImageGeometry(voxel_num_x=N, + voxel_num_y=N, + voxel_size_x=voxel_size_h, + voxel_size_y=voxel_size_h) + +# Set up the Projector (AcquisitionModel) using ASTRA on GPU +Aop = AstraProjectorSimple(ig2d, ag2d,"gpu") + +# Set initial guess for CGLS reconstruction +x_init = ImageData(np.zeros((N,N)),geometry=ig2d) + +# Run 50-iteration CGLS reconstruction +num_iter = 50 +x, it, timing, criter = CGLS(Aop,data2d,num_iter,x_init) + +# Display reconstruction +plt.figure() +plt.imshow(x.as_array()) \ No newline at end of file diff --git a/Wrappers/Python/wip/demo_sophiabeads.py b/Wrappers/Python/wip/demo_sophiabeads.py deleted file mode 100755 index 4a0f963..0000000 --- a/Wrappers/Python/wip/demo_sophiabeads.py +++ /dev/null @@ -1,73 +0,0 @@ - -# 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: -# 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. - -# Do all imports -from ccpi.io.reader import XTEKReader -import numpy as np -import matplotlib.pyplot as plt -from ccpi.framework import ImageGeometry, AcquisitionGeometry, AcquisitionData, ImageData -from ccpi.astra.ops import AstraProjectorSimple -from ccpi.optimisation.algs import CGLS - -# 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() - -# Extract central slice, scale and negative-log transform -sino = -np.log(data.as_array()[:,:,1000]/60000.0) - -# Apply centering correction by zero padding, amount found manually -cor_pad = 30 -sino_pad = np.zeros((sino.shape[0],sino.shape[1]+cor_pad)) -sino_pad[:,cor_pad:] = sino - -# Simple beam hardening correction as done in SophiaBeads coda -sino_pad = sino_pad**2 - -# Extract AcquisitionGeometry for central slice for 2D fanbeam reconstruction -ag2d = AcquisitionGeometry('cone', - '2D', - angles=-np.pi/180*data.geometry.angles, - pixel_num_h=data.geometry.pixel_num_h + cor_pad, - pixel_size_h=data.geometry.pixel_size_h, - dist_source_center=data.geometry.dist_source_center, - dist_center_detector=data.geometry.dist_center_detector) - -# Set up AcquisitionData object for central slice 2D fanbeam -data2d = AcquisitionData(sino_pad,geometry=ag2d) - -# Choose the number of voxels to reconstruct onto as number of detector pixels -N = data.geometry.pixel_num_h - -# Geometric magnification -mag = (np.abs(data.geometry.dist_center_detector) + \ - np.abs(data.geometry.dist_source_center)) / \ - np.abs(data.geometry.dist_source_center) - -# Voxel size is detector pixel size divided by mag -voxel_size_h = data.geometry.pixel_size_h / mag - -# Construct the appropriate ImageGeometry -ig2d = ImageGeometry(voxel_num_x=N, - voxel_num_y=N, - voxel_size_x=voxel_size_h, - voxel_size_y=voxel_size_h) - -# Set up the Projector (AcquisitionModel) using ASTRA on GPU -Aop = AstraProjectorSimple(ig2d, ag2d,"gpu") - -# Set initial guess for CGLS reconstruction -x_init = ImageData(np.zeros((N,N)),geometry=ig2d) - -# Run 50-iteration CGLS reconstruction -num_iter = 50 -x, it, timing, criter = CGLS(Aop,data2d,num_iter,x_init) - -# Display reconstruction -plt.figure() -plt.imshow(x.as_array()) \ No newline at end of file diff --git a/Wrappers/Python/wip/simple_demo_astra.py b/Wrappers/Python/wip/simple_demo_astra.py deleted file mode 100755 index 8529c98..0000000 --- a/Wrappers/Python/wip/simple_demo_astra.py +++ /dev/null @@ -1,219 +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 -from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D -from ccpi.astra.ops 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.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} - -# First a CGLS reconstruction can be done: -x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt) - -plt.imshow(x_CGLS.array) -plt.title('CGLS') -plt.show() - -plt.semilogy(criter_CGLS) -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 regularization -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None,opt) - -plt.imshow(x_fista0.array) -plt.title('FISTA Least squares') -plt.show() - -plt.semilogy(criter0) -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) - -# 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': 7000} -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 = 1 -gtv = TV2D(lamtv) - -x_fbpdtv,it_fbpdtv,timing_fbpdtv,criter_fbpdtv=FBPD(x_init,None,f,gtv,opt_FBPD) - -plt.imshow(x_fbpdtv.array) -plt.show() - -plt.semilogy(criter_fbpdtv) -plt.show() - - -# Compare all reconstruction and criteria -clims = (0,1) -cols = 3 -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_fista0.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]) -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') -plt.show() diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py deleted file mode 100755 index 3b78eb3..0000000 --- a/Wrappers/Python/wip/simple_mc_demo.py +++ /dev/null @@ -1,137 +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.algs import FISTA -from ccpi.optimisation.funcs import Norm2sq, Norm1 -from ccpi.astra.ops 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 = ImageData(np.zeros(x.shape),geometry=ig) -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 -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) - -# Display reconstruction and criteration -ff0, axarrf0 = plt.subplots(1,numchannels) -for k in numpy.arange(3): - axarrf0[k].imshow(x_fista0.as_array()[k],vmin=0,vmax=2.5) -plt.show() - -plt.semilogy(criter0) -plt.title('Criterion vs iterations, least squares') -plt.show() - -# 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 = 0.1 -g0 = Norm1(lam) - -# Run FISTA for least squares plus 1-norm function. -x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) - -# Display reconstruction and criteration -ff1, axarrf1 = plt.subplots(1,numchannels) -for k in numpy.arange(3): - axarrf1[k].imshow(x_fista1.as_array()[k],vmin=0,vmax=2.5) -plt.show() - -plt.semilogy(criter1) -plt.title('Criterion vs iterations, least squares plus 1-norm regu') -plt.show() -- cgit v1.2.3