From dac9d2ce209487765f66332547ce59f92a6c1d78 Mon Sep 17 00:00:00 2001 From: Jakob Jorgensen Date: Tue, 24 Apr 2018 11:50:53 +0100 Subject: Partial cleanup demo --- Wrappers/Python/wip/simple_demo_ccpi.py | 98 +++++++++++++++------------------ 1 file changed, 45 insertions(+), 53 deletions(-) diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py index 3fdc2d4..6344c06 100755 --- a/Wrappers/Python/wip/simple_demo_ccpi.py +++ b/Wrappers/Python/wip/simple_demo_ccpi.py @@ -1,48 +1,44 @@ -#import sys -#sys.path.append("..") -from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry +# This demo illustrates how CCPi 2D parallel-beam projectors can be used with +# the modular optimisation framework. The demo sets up a small 4-slice 3D test +# case and demonstrates reconstruction using CGLS, as well as FISTA for least +# squares and 1-norm regularisation and FBPD for 1-norm regularisation. + +# First make all imports +from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, \ + AcquisitionGeometry from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, Norm1 , TV2D +from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D from ccpi.plugins.ops import CCPiProjectorSimple from ccpi.plugins.processors import CCPiForwardProjector, CCPiBackwardProjector + from ccpi.reconstruction.parallelbeam import alg as pbalg import numpy as np import matplotlib.pyplot as plt -test_case = 1 # 1=parallel2D, 2=cone2D, 3=parallel3D +# Set up phantom size N x N x vert by creating ImageGeometry, initialising the +# ImageData object with this geometry and empty array and finally put some +# data into its array, and display one slice as image. -# Set up phantom +# Image parameters N = 128 vert = 4 -# Set up measurement geometry -angles_num = 20; # angles number -det_w = 1.0 -det_num = N -SourceOrig = 200 -OrigDetec = 0 - -if test_case==1: - angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*180/np.pi - #nangles = angles_num - #angles = np.linspace(0,360, nangles, dtype=np.float32) - -elif test_case==2: - angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) -elif test_case == 3: - angles = np.linspace(0,np.pi,angles_num,endpoint=False) -else: - NotImplemented -vg = ImageGeometry(voxel_num_x=N, +# Set up image geometry +ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=vert) -Phantom = ImageData(geometry=vg,dimension_labels=['horizontal_x','horizontal_y','vertical']) + 0.1 - +# Set up empty image data +Phantom = ImageData(geometry=ig, + dimension_labels=['horizontal_x', + 'horizontal_y', + 'vertical']) + +# Populate image data by looping over and filling slices i = 0 while i < vert: if vert > 1: @@ -55,6 +51,7 @@ while i < vert: Phantom.fill(x, vertical=i) i += 1 +# Display slice of phantom if vert > 1: plt.imshow(Phantom.subset(vertical=0).as_array()) else: @@ -62,33 +59,28 @@ else: 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, set the width of a detector +# pixel relative to an object pixe and the number of detector pixels. +angles_num = 20; # angles number +det_w = 1.0 +det_num = N -# Parallelbeam geometry test -if test_case==1: - #Phantom_ccpi = Phantom.subset(['horizontal_x','horizontal_y','vertical']) - #Phantom_ccpi.geometry = vg.clone() - center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2 - - pg = AcquisitionGeometry('parallel', - '3D', - angles, - N , det_w, - vert , det_w #2D in 3D is a slice 1 pixel thick - ) -elif test_case==2: - raise NotImplemented('cone beam projector not yet available') - pg = AcquisitionGeometry('cone', - '2D', - angles, - det_num, - det_w, - vert, det_w, #2D in 3D is a slice 1 pixel thick - dist_source_center=SourceOrig, - dist_center_detector=OrigDetec) - -# ASTRA operator using volume and sinogram geometries -#Aop = AstraProjectorSimple(vg, pg, 'cpu') -Cop = CCPiProjectorSimple(vg, pg) +angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*180/np.pi + +#center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2 + +ag = AcquisitionGeometry('parallel', + '3D', + angles, + N, + det_w, + vert, + det_w) + +# CCPi operator using image and acquisition geometries +Cop = CCPiProjectorSimple(ig, ag) # Try forward and backprojection b = Cop.direct(Phantom) -- cgit v1.2.3 From f911ed4a7df32dcad37330016180655721dbbd22 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 24 Apr 2018 13:58:24 +0100 Subject: Tidy up demo, remove prints and old code. --- Wrappers/Python/ccpi/plugins/ops.py | 1 - Wrappers/Python/ccpi/plugins/processors.py | 31 +---- Wrappers/Python/wip/simple_demo_ccpi.py | 176 +++++++++++++++-------------- 3 files changed, 94 insertions(+), 114 deletions(-) diff --git a/Wrappers/Python/ccpi/plugins/ops.py b/Wrappers/Python/ccpi/plugins/ops.py index aeb51af..75c5db9 100755 --- a/Wrappers/Python/ccpi/plugins/ops.py +++ b/Wrappers/Python/ccpi/plugins/ops.py @@ -109,6 +109,5 @@ class CCPiProjectorSimple(Operator): x0 = ImageData(geometry = self.volume_geometry, dimension_labels=self.bp.output_axes_order)#\ #.subset(['horizontal_x','horizontal_y','vertical']) - print (x0) x0.fill(numpy.random.randn(*x0.shape)) return x0 \ No newline at end of file diff --git a/Wrappers/Python/ccpi/plugins/processors.py b/Wrappers/Python/ccpi/plugins/processors.py index df580e0..9938b9e 100755 --- a/Wrappers/Python/ccpi/plugins/processors.py +++ b/Wrappers/Python/ccpi/plugins/processors.py @@ -17,15 +17,10 @@ # See the License for the specific language governing permissions and # limitations under the License -from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\ - AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.framework import DataProcessor, AcquisitionData,\ + AcquisitionGeometry, ImageGeometry, ImageData from ccpi.reconstruction.parallelbeam import alg as pbalg import numpy -import h5py -from scipy import ndimage - -import matplotlib.pyplot as plt - def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter): @@ -53,10 +48,9 @@ def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter): angles, center_of_rotation, voxel_per_pixel ) - - #print (counter) + counter+=1 - #print (geoms , geoms_i) + if counter < 4: if (not ( geoms_i == geoms )): print ("not equal and {0}".format(counter)) @@ -65,10 +59,8 @@ def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter): Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) return setupCCPiGeometries(X,Y,Z,angles, counter) else: - print ("return geoms {0}".format(geoms)) return geoms else: - print ("return geoms_i {0}".format(geoms_i)) return geoms_i @@ -119,18 +111,6 @@ class CCPiForwardProjector(DataProcessor): pixel_per_voxel = 1 # should be estimated from image_geometry and # acquisition_geometry if self.acquisition_geometry.geom_type == 'parallel': - #int msize = ndarray_volume.shape(0) > ndarray_volume.shape(1) ? ndarray_volume.shape(0) : ndarray_volume.shape(1); - #int detector_width = msize; - # detector_width is the max between the shape[0] and shape[1] - - - #double rotation_center = (double)detector_width/2.; - #int detector_height = ndarray_volume.shape(2); - - #int number_of_projections = ndarray_angles.shape(0); - - ##numpy_3d pixels(reinterpret_cast(ndarray_volume.get_data()), - #boost::extents[number_of_projections][detector_height][detector_width]); pixels = pbalg.pb_forward_project(volume.as_array(), self.acquisition_geometry.angles, @@ -179,7 +159,6 @@ class CCPiBackwardProjector(DataProcessor): def check_input(self, dataset): if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: - #number_of_projections][detector_height][detector_width return True else: @@ -198,7 +177,7 @@ class CCPiBackwardProjector(DataProcessor): voxel_num_z = self.acquisition_geometry.pixel_num_v) # input centered/padded acquisitiondata center_of_rotation = projections.get_dimension_size('horizontal') / 2 - #print (center_of_rotation) + if self.acquisition_geometry.geom_type == 'parallel': back = pbalg.pb_backward_project( projections.as_array(), diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py index 6344c06..a8265ce 100755 --- a/Wrappers/Python/wip/simple_demo_ccpi.py +++ b/Wrappers/Python/wip/simple_demo_ccpi.py @@ -5,16 +5,12 @@ # squares and 1-norm regularisation and FBPD for 1-norm regularisation. # First make all imports -from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, \ - AcquisitionGeometry +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.funcs import Norm2sq, Norm1 from ccpi.plugins.ops import CCPiProjectorSimple -from ccpi.plugins.processors import CCPiForwardProjector, CCPiBackwardProjector - -from ccpi.reconstruction.parallelbeam import alg as pbalg import numpy as np import matplotlib.pyplot as plt @@ -63,14 +59,16 @@ plt.show() # setup geometry: # Number of angles, the actual angles from 0 to # pi for parallel beam, set the width of a detector # pixel relative to an object pixe and the number of detector pixels. -angles_num = 20; # angles number +angles_num = 20 det_w = 1.0 det_num = N -angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*180/np.pi - -#center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2 +angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ + 180/np.pi +# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, +# horz detector pixel size, vert detector pixel count, +# vert detector pixel size. ag = AcquisitionGeometry('parallel', '3D', angles, @@ -79,142 +77,146 @@ ag = AcquisitionGeometry('parallel', vert, det_w) -# CCPi operator using image and acquisition geometries +# Set up Operator object combining the ImageGeometry and AcquisitionGeometry +# wrapping calls to CCPi projector. Cop = CCPiProjectorSimple(ig, ag) -# Try forward and backprojection +# Forward and backprojection are available as methods direct and adjoint. Here +# generate test data b and do simple backprojection to obtain z. Display all +# data slices as images, and a single backprojected slice. b = Cop.direct(Phantom) -out2 = Cop.adjoint(b) +z = Cop.adjoint(b) -#%% for i in range(b.get_dimension_size('vertical')): plt.imshow(b.subset(vertical=i).array) - #plt.imshow(Phantom.subset( vertical=i).array) - #plt.imshow(b.array[:,i,:]) plt.show() -#%% -plt.imshow(out2.subset( vertical=0).array) +plt.imshow(z.subset(vertical=0).array) +plt.title('Backprojected data') plt.show() -# Create least squares object instance with projector and data. +# 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. Note that 100 iterations for +# some of the methods is a very low number and 1000 or 10000 iterations may be +# needed if one wants to obtain a converged solution. +x_init = ImageData(geometry=ig, + dimension_labels=['horizontal_x','horizontal_y','vertical']) +opt = {'tol': 1e-4, 'iter': 100} + +# First a CGLS reconstruction can be done: +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, b, opt=opt) + +plt.imshow(x_CGLS.subset(vertical=0).array) +plt.title('CGLS') +plt.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(Cop,b,c=0.5) -# Initial guess -x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y','vertical']) -#invL = 0.5 -#g = f.grad(x_init) -#print (g) -#u = x_init - invL*f.grad(x_init) - -#%% # Run FISTA for least squares without regularization -opt = {'tol': 1e-4, 'iter': 100} x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) plt.imshow(x_fista0.subset(vertical=0).array) -plt.title('FISTA0') +plt.title('FISTA Least squares') plt.show() -# Now least squares plus 1-norm regularization +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=opt) +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) -plt.imshow(x_fista0.subset(vertical=0).array) -plt.title('FISTA1') +plt.imshow(x_fista1.subset(vertical=0).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() -# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm +# 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. x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) plt.imshow(x_fbpd1.subset(vertical=0).array) -plt.title('FBPD1') +plt.title('FBPD for least squares plus 1-norm regularisation') plt.show() plt.semilogy(criter_fbpd1) -plt.show() - -# Now FBPD for least squares plus TV -#lamtv = 1 -#gtv = TV2D(lamtv) - -#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt) - -#plt.imshow(x_fbpdtv.subset(vertical=0).array) -#plt.show() - -#plt.semilogy(criter_fbpdtv) -#plt.show() - - -# Run CGLS, which should agree with the FISTA0 -x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, b, opt=opt) - -plt.imshow(x_CGLS.subset(vertical=0).array) -plt.title('CGLS') -plt.title('CGLS recon, compare FISTA0') -plt.show() - -plt.semilogy(criter_CGLS) -plt.title('CGLS criterion') +plt.title('FBPD for least squares plus 1-norm regularisation criterion') plt.show() -#%% +# Compare all reconstruction and criteria clims = (0,1) cols = 3 rows = 2 current = 1 + fig = plt.figure() -# projections row a=fig.add_subplot(rows,cols,current) a.set_title('phantom {0}'.format(np.shape(Phantom.as_array()))) - -imgplot = plt.imshow(Phantom.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) +imgplot = plt.imshow(Phantom.subset(vertical=0).as_array(), + vmin=clims[0],vmax=clims[1]) +plt.axis('off') current = current + 1 a=fig.add_subplot(rows,cols,current) -a.set_title('FISTA0') -imgplot = plt.imshow(x_fista0.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) +a.set_title('CGLS') +imgplot = plt.imshow(x_CGLS.subset(vertical=0).as_array(), + vmin=clims[0],vmax=clims[1]) +plt.axis('off') current = current + 1 a=fig.add_subplot(rows,cols,current) -a.set_title('FISTA1') -imgplot = plt.imshow(x_fista1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) +a.set_title('FISTA LS') +imgplot = plt.imshow(x_fista0.subset(vertical=0).as_array(), + vmin=clims[0],vmax=clims[1]) +plt.axis('off') current = current + 1 a=fig.add_subplot(rows,cols,current) -a.set_title('FBPD1') -imgplot = plt.imshow(x_fbpd1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) +a.set_title('FISTA LS+1') +imgplot = plt.imshow(x_fista1.subset(vertical=0).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.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) - -plt.show() -#%% -#current = current + 1 -#a=fig.add_subplot(rows,cols,current) -#a.set_title('FBPD TV') -#imgplot = plt.imshow(x_fbpdtv.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) +a.set_title('FBPD LS+1') +imgplot = plt.imshow(x_fbpd1.subset(vertical=0).as_array(), + vmin=clims[0],vmax=clims[1]) +plt.axis('off') fig = plt.figure() -# projections row b=fig.add_subplot(1,1,1) b.set_title('criteria') -imgplot = plt.loglog(criter0 , label='FISTA0') -imgplot = plt.loglog(criter1 , label='FISTA1') -imgplot = plt.loglog(criter_fbpd1, label='FBPD1') imgplot = plt.loglog(criter_CGLS, label='CGLS') -#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV') -b.legend(loc='right') -plt.show() -#%% \ No newline at end of file +imgplot = plt.loglog(criter0 , label='FISTA LS') +imgplot = plt.loglog(criter1 , label='FISTA LS+1') +imgplot = plt.loglog(criter_fbpd1, label='FBPD LS+1') +b.legend(loc='lower left') +plt.show() \ No newline at end of file -- cgit v1.2.3 From e6fb1072bfbc952c01d5f1df826a6371ff1cd2cc Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 24 Apr 2018 13:59:18 +0100 Subject: Renamed demo consistently --- Wrappers/Python/wip/demo_ccpi_simple.py | 222 ++++++++++++++++++++++++++++++++ Wrappers/Python/wip/simple_demo_ccpi.py | 222 -------------------------------- 2 files changed, 222 insertions(+), 222 deletions(-) create mode 100755 Wrappers/Python/wip/demo_ccpi_simple.py delete mode 100755 Wrappers/Python/wip/simple_demo_ccpi.py diff --git a/Wrappers/Python/wip/demo_ccpi_simple.py b/Wrappers/Python/wip/demo_ccpi_simple.py new file mode 100755 index 0000000..a8265ce --- /dev/null +++ b/Wrappers/Python/wip/demo_ccpi_simple.py @@ -0,0 +1,222 @@ + +# This demo illustrates how CCPi 2D parallel-beam projectors can be used with +# the modular optimisation framework. The demo sets up a small 4-slice 3D test +# case and demonstrates reconstruction using CGLS, as well as FISTA for least +# squares and 1-norm regularisation and FBPD for 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 + +from ccpi.plugins.ops import CCPiProjectorSimple + +import numpy as np +import matplotlib.pyplot as plt + +# Set up phantom size N x N x vert by creating ImageGeometry, initialising the +# ImageData object with this geometry and empty array and finally put some +# data into its array, and display one slice as image. + +# Image parameters +N = 128 +vert = 4 + +# Set up image geometry +ig = ImageGeometry(voxel_num_x=N, + voxel_num_y=N, + voxel_num_z=vert) + +# Set up empty image data +Phantom = ImageData(geometry=ig, + dimension_labels=['horizontal_x', + 'horizontal_y', + 'vertical']) + +# Populate image data by looping over and filling slices +i = 0 +while i < vert: + if vert > 1: + x = Phantom.subset(vertical=i).array + else: + x = Phantom.array + x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 + x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 + +# Display slice of phantom +if vert > 1: + plt.imshow(Phantom.subset(vertical=0).as_array()) +else: + plt.imshow(Phantom.as_array()) +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, set the width of a detector +# pixel relative to an object pixe and the number of detector pixels. +angles_num = 20 +det_w = 1.0 +det_num = N + +angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ + 180/np.pi + +# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, +# horz detector pixel size, vert detector pixel count, +# vert detector pixel size. +ag = AcquisitionGeometry('parallel', + '3D', + angles, + N, + det_w, + vert, + det_w) + +# Set up Operator object combining the ImageGeometry and AcquisitionGeometry +# wrapping calls to CCPi projector. +Cop = CCPiProjectorSimple(ig, ag) + +# Forward and backprojection are available as methods direct and adjoint. Here +# generate test data b and do simple backprojection to obtain z. Display all +# data slices as images, and a single backprojected slice. +b = Cop.direct(Phantom) +z = Cop.adjoint(b) + +for i in range(b.get_dimension_size('vertical')): + plt.imshow(b.subset(vertical=i).array) + plt.show() + +plt.imshow(z.subset(vertical=0).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. Note that 100 iterations for +# some of the methods is a very low number and 1000 or 10000 iterations may be +# needed if one wants to obtain a converged solution. +x_init = ImageData(geometry=ig, + dimension_labels=['horizontal_x','horizontal_y','vertical']) +opt = {'tol': 1e-4, 'iter': 100} + +# First a CGLS reconstruction can be done: +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, b, opt=opt) + +plt.imshow(x_CGLS.subset(vertical=0).array) +plt.title('CGLS') +plt.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(Cop,b,c=0.5) + +# Run FISTA for least squares without regularization +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) + +plt.imshow(x_fista0.subset(vertical=0).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.subset(vertical=0).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. +x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) + +plt.imshow(x_fbpd1.subset(vertical=0).array) +plt.title('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() + + +# 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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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') +b.legend(loc='lower left') +plt.show() \ No newline at end of file diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py deleted file mode 100755 index a8265ce..0000000 --- a/Wrappers/Python/wip/simple_demo_ccpi.py +++ /dev/null @@ -1,222 +0,0 @@ - -# This demo illustrates how CCPi 2D parallel-beam projectors can be used with -# the modular optimisation framework. The demo sets up a small 4-slice 3D test -# case and demonstrates reconstruction using CGLS, as well as FISTA for least -# squares and 1-norm regularisation and FBPD for 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 - -from ccpi.plugins.ops import CCPiProjectorSimple - -import numpy as np -import matplotlib.pyplot as plt - -# Set up phantom size N x N x vert by creating ImageGeometry, initialising the -# ImageData object with this geometry and empty array and finally put some -# data into its array, and display one slice as image. - -# Image parameters -N = 128 -vert = 4 - -# Set up image geometry -ig = ImageGeometry(voxel_num_x=N, - voxel_num_y=N, - voxel_num_z=vert) - -# Set up empty image data -Phantom = ImageData(geometry=ig, - dimension_labels=['horizontal_x', - 'horizontal_y', - 'vertical']) - -# Populate image data by looping over and filling slices -i = 0 -while i < vert: - if vert > 1: - x = Phantom.subset(vertical=i).array - else: - x = Phantom.array - x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 - x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98 - if vert > 1 : - Phantom.fill(x, vertical=i) - i += 1 - -# Display slice of phantom -if vert > 1: - plt.imshow(Phantom.subset(vertical=0).as_array()) -else: - plt.imshow(Phantom.as_array()) -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, set the width of a detector -# pixel relative to an object pixe and the number of detector pixels. -angles_num = 20 -det_w = 1.0 -det_num = N - -angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ - 180/np.pi - -# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, -# horz detector pixel size, vert detector pixel count, -# vert detector pixel size. -ag = AcquisitionGeometry('parallel', - '3D', - angles, - N, - det_w, - vert, - det_w) - -# Set up Operator object combining the ImageGeometry and AcquisitionGeometry -# wrapping calls to CCPi projector. -Cop = CCPiProjectorSimple(ig, ag) - -# Forward and backprojection are available as methods direct and adjoint. Here -# generate test data b and do simple backprojection to obtain z. Display all -# data slices as images, and a single backprojected slice. -b = Cop.direct(Phantom) -z = Cop.adjoint(b) - -for i in range(b.get_dimension_size('vertical')): - plt.imshow(b.subset(vertical=i).array) - plt.show() - -plt.imshow(z.subset(vertical=0).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. Note that 100 iterations for -# some of the methods is a very low number and 1000 or 10000 iterations may be -# needed if one wants to obtain a converged solution. -x_init = ImageData(geometry=ig, - dimension_labels=['horizontal_x','horizontal_y','vertical']) -opt = {'tol': 1e-4, 'iter': 100} - -# First a CGLS reconstruction can be done: -x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, b, opt=opt) - -plt.imshow(x_CGLS.subset(vertical=0).array) -plt.title('CGLS') -plt.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(Cop,b,c=0.5) - -# Run FISTA for least squares without regularization -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) - -plt.imshow(x_fista0.subset(vertical=0).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.subset(vertical=0).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. -x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) - -plt.imshow(x_fbpd1.subset(vertical=0).array) -plt.title('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() - - -# 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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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') -b.legend(loc='lower left') -plt.show() \ No newline at end of file -- cgit v1.2.3