diff options
3 files changed, 26 insertions, 29 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py index 6cf4f06..1695a73 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -161,12 +161,12 @@ class CGLS(Algorithm): def should_stop(self): + self.update_objective() flag = (self.norms <= self.norms0 * self.tolerance) or (self.normx * self.tolerance >= 1); #if self.gamma<=self.tolerance: if flag == 1 or self.max_iteration_stop_cryterion(): - print('Tolerance is reached: Iter: {}'.format(self.iteration)) - self.update_objective() + print('Tolerance is reached: Iter: {}'.format(self.iteration)) return True diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py index e8f4063..2ac050f 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 = 128 # set dimension of the phantom +N = 64 # 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) @@ -75,7 +75,6 @@ sin = Aop.direct(data) np.random.seed(10) noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,1,ag.shape)) -#noisy_data = AcquisitionData( sin.as_array() ) # Show Ground Truth and Noisy Data plt.figure(figsize=(10,10)) @@ -89,7 +88,6 @@ plt.title('Noisy Data') plt.colorbar() plt.show() - # Setup and run the simple CGLS algorithm x_init = ig.allocate() @@ -98,7 +96,6 @@ cgls1.max_iteration = 20 cgls1.update_objective_interval = 5 cgls1.run(20, verbose = True) - # Setup and run the regularised CGLS algorithm (Tikhonov with Identity) x_init = ig.allocate() diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py index 568df38..09e350b 100644 --- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py +++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py @@ -20,8 +20,8 @@ #========================================================================= """ -Compare solutions of FISTA & PDHG - & CGLS & Astra Built-in algorithms for Least Squares +Compare solutions of FISTA & PDHG & CGLS + & Astra Built-in algorithms for Least Squares Problem: min_x || A x - g ||_{2}^{2} @@ -32,7 +32,7 @@ Problem: min_x || A x - g ||_{2}^{2} """ -from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry import numpy as np import numpy @@ -71,10 +71,9 @@ else: dev = 'cpu' Aop = AstraProjectorSimple(ig, ag, dev) -sin = Aop.direct(data) +sinogram = Aop.direct(data) + -np.random.seed(10) -noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,1,ag.shape)) ############################################################################### # Setup and run Astra CGLS algorithm @@ -82,10 +81,10 @@ vol_geom = astra.create_vol_geom(N, N) proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) proj_id = astra.create_projector('line', proj_geom, vol_geom) -# Create a sinogram from a phantom -sinogram_id = astra.data2d.create('-sino', proj_geom, noisy_data.as_array()) +# Create a sinogram id +sinogram_id = astra.data2d.create('-sino', proj_geom, sinogram.as_array()) -# Create a data object for the reconstruction +# Create a data id rec_id = astra.data2d.create('-vol', vol_geom) cgls_astra = astra.astra_dict('CGLS') @@ -96,7 +95,7 @@ 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, 25) +astra.algorithm.run(alg_id, 500) recon_cgls_astra = ImageData(astra.data2d.get(rec_id)) @@ -105,41 +104,42 @@ recon_cgls_astra = ImageData(astra.data2d.get(rec_id)) # 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) +cgls = CGLS(x_init = x_init, operator = Aop, data = sinogram) +cgls.max_iteration = 500 +cgls.update_objective_interval = 100 +cgls.run(500, verbose = True) #%% ############################################################################### # Setup and run the PDHG algorithm operator = Aop -f = L2NormSquared(b = noisy_data) +f = L2NormSquared(b = sinogram) g = ZeroFunction() ## Compute operator Norm normK = operator.norm() ## Primal & dual stepsizes -sigma = 1 +sigma = 0.02 tau = 1/(sigma*normK**2) + pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 400 -pdhg.update_objective_interval = 20 -pdhg.run(400, verbose=True) +pdhg.max_iteration = 1000 +pdhg.update_objective_interval = 100 +pdhg.run(1000, verbose=True) #%% ############################################################################### # Setup and run the FISTA algorithm -fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop) +fidelity = FunctionOperatorComposition(L2NormSquared(b=sinogram), Aop) regularizer = ZeroFunction() fista = FISTA(x_init=x_init , f=fidelity, g=regularizer) -fista.max_iteration = 20 -fista.update_objective_interval = 5 -fista.run(20, verbose = True) +fista.max_iteration = 500 +fista.update_objective_interval = 100 +fista.run(500, verbose = True) #%% Show results |