diff options
| -rw-r--r-- | Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py | 4 | ||||
| -rw-r--r-- | Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py | 98 |
2 files changed, 56 insertions, 46 deletions
diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py index abcb26c..e8f4063 100644 --- a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py +++ b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py @@ -50,7 +50,7 @@ import os # Load Shepp-Logan phantom model = 1 # select a model number from the library -N = 64 # set dimension of the phantom +N = 128 # set dimension of the phantom path = os.path.dirname(tomophantom.__file__) path_library2D = os.path.join(path, "Phantom2DLibrary.dat") phantom_2D = TomoP2D.Model(model, N, path_library2D) @@ -59,7 +59,7 @@ ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) data = ImageData(phantom_2D) detectors = N -angles = np.linspace(0, np.pi, 90, dtype=np.float32) +angles = np.linspace(0, np.pi, 180, dtype=np.float32) ag = AcquisitionGeometry('parallel','2D', angles, detectors) diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py index 672d4bc..568df38 100644 --- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py +++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py @@ -32,7 +32,7 @@ Problem: min_x || A x - g ||_{2}^{2} """ -from ccpi.framework import ImageData, TestData, AcquisitionGeometry +from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry import numpy as np import numpy @@ -43,40 +43,47 @@ from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition from ccpi.astra.ops import AstraProjectorSimple import astra -import os, sys +import tomophantom +from tomophantom import TomoP2D +import os -# Load Data -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) - -N = 50 -M = 50 -data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) -ig = data.geometry +# Load Shepp-Logan phantom +model = 1 # select a model number from the library +N = 128 # set dimension of the phantom +path = os.path.dirname(tomophantom.__file__) +path_library2D = os.path.join(path, "Phantom2DLibrary.dat") +phantom_2D = TomoP2D.Model(model, N, path_library2D) -detectors = N -angles = np.linspace(0, np.pi, N, dtype=np.float32) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +data = ImageData(phantom_2D) + +detectors = N +angles = np.linspace(0, np.pi, 180, dtype=np.float32) -device = input('Available device: GPU==1 / CPU==0 ') ag = AcquisitionGeometry('parallel','2D', angles, detectors) -if device=='1': + +device = input('Available device: GPU==1 / CPU==0 ') + +if device =='1': dev = 'gpu' else: dev = 'cpu' - -Aop = AstraProjectorSimple(ig, ag, dev) + +Aop = AstraProjectorSimple(ig, ag, dev) sin = Aop.direct(data) -noisy_data = sin +np.random.seed(10) +noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,1,ag.shape)) ############################################################################### # Setup and run Astra CGLS algorithm vol_geom = astra.create_vol_geom(N, N) proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) -proj_id = astra.create_projector('linear', proj_geom, vol_geom) +proj_id = astra.create_projector('line', proj_geom, vol_geom) # Create a sinogram from a phantom -sinogram_id, sinogram = astra.create_sino(data.as_array(), proj_id) +sinogram_id = astra.data2d.create('-sino', proj_geom, noisy_data.as_array()) # Create a data object for the reconstruction rec_id = astra.data2d.create('-vol', vol_geom) @@ -89,17 +96,21 @@ cgls_astra['ProjectorId'] = proj_id # Create the algorithm object from the configuration structure alg_id = astra.algorithm.create(cgls_astra) -astra.algorithm.run(alg_id, 1000) +astra.algorithm.run(alg_id, 25) +recon_cgls_astra = ImageData(astra.data2d.get(rec_id)) -recon_cgls_astra = astra.data2d.get(rec_id) -############################################################################### -# Setup and run the CGLS algorithm -x_init = ig.allocate() -cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) -cgls.max_iteration = 1000 -cgls.update_objective_interval = 200 -cgls.run(1000, verbose=False) +#%% + +# Setup and run the simple CGLS algorithm +x_init = ig.allocate() + +cgls = CGLS(x_init = x_init, operator = Aop, data = noisy_data) +cgls.max_iteration = 25 +cgls.update_objective_interval = 5 +cgls.run(25, verbose = True) + +#%% ############################################################################### # Setup and run the PDHG algorithm @@ -115,19 +126,20 @@ sigma = 1 tau = 1/(sigma*normK**2) pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 1000 -pdhg.update_objective_interval = 200 -pdhg.run(1000, verbose=True) +pdhg.max_iteration = 400 +pdhg.update_objective_interval = 20 +pdhg.run(400, verbose=True) +#%% ############################################################################### # Setup and run the FISTA algorithm fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop) regularizer = ZeroFunction() fista = FISTA(x_init=x_init , f=fidelity, g=regularizer) -fista.max_iteration = 1000 -fista.update_objective_interval = 200 -fista.run(1000, verbose=True) +fista.max_iteration = 20 +fista.update_objective_interval = 5 +fista.run(20, verbose = True) #%% Show results @@ -150,40 +162,38 @@ plt.colorbar() plt.title('PDHG reconstruction') plt.subplot(2,2,4) -plt.imshow(recon_cgls_astra) +plt.imshow(recon_cgls_astra.as_array()) plt.colorbar() plt.title('CGLS astra') -diff1 = pdhg.get_output() - cgls.get_output() -diff2 = fista.get_output() - cgls.get_output() -diff3 = ImageData(recon_cgls_astra) - cgls.get_output() +diff1 = pdhg.get_output() - recon_cgls_astra +diff2 = fista.get_output() - recon_cgls_astra +diff3 = cgls.get_output() - recon_cgls_astra plt.figure(figsize=(15,15)) plt.subplot(3,1,1) plt.imshow(diff1.abs().as_array()) -plt.title('Diff PDHG vs CGLS') +plt.title('Diff PDHG vs CGLS astra') plt.colorbar() plt.subplot(3,1,2) plt.imshow(diff2.abs().as_array()) -plt.title('Diff FISTA vs CGLS') +plt.title('Diff FISTA vs CGLS astra') plt.colorbar() plt.subplot(3,1,3) plt.imshow(diff3.abs().as_array()) -plt.title('Diff CLGS astra vs CGLS') +plt.title('Diff CLGS vs CGLS astra') plt.colorbar() -#%% +cgls_astra_obj = fidelity(ImageData(recon_cgls_astra)) print('Primal Objective (FISTA) {} '.format(fista.objective[-1])) print('Primal Objective (CGLS) {} '.format(cgls.objective[-1])) print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0])) +print('Primal Objective (CGLS_astra) {} '.format(cgls_astra_obj)) -true_obj = (Aop.direct(cglsd.get_output())-noisy_data).squared_norm() -print('True objective {}'.format(true_obj)) - |