diff options
5 files changed, 61 insertions, 79 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py index 04e7c38..419958a 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -67,7 +67,8 @@ class FISTA(Algorithm): self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) - self.x.subtract(self.x_old, out=self.y) +# self.x.subtract(self.x_old, out=self.y) + self.y = self.x - self.x_old self.y.__imul__ ((self.t_old-1)/self.t) self.y.__iadd__( self.x ) diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py index 0875c20..672d4bc 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, ImageGeometry, AcquisitionGeometry +from ccpi.framework import ImageData, TestData, AcquisitionGeometry import numpy as np import numpy @@ -42,17 +42,17 @@ 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 astra +import os, sys -# Create Ground truth phantom and Sinogram -N = 128 -x = np.zeros((N,N)) -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 - -data = ImageData(x) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +# 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 detectors = N angles = np.linspace(0, np.pi, N, dtype=np.float32) @@ -69,13 +69,14 @@ sin = Aop.direct(data) noisy_data = sin +############################################################################### # 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('strip', proj_geom, vol_geom) +proj_id = astra.create_projector('linear', proj_geom, vol_geom) # Create a sinogram from a phantom -sinogram_id, sinogram = astra.create_sino(x, proj_id) +sinogram_id, sinogram = astra.create_sino(data.as_array(), proj_id) # Create a data object for the reconstruction rec_id = astra.data2d.create('-vol', vol_geom) @@ -92,6 +93,7 @@ astra.algorithm.run(alg_id, 1000) 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) @@ -99,12 +101,15 @@ cgls.max_iteration = 1000 cgls.update_objective_interval = 200 cgls.run(1000, verbose=False) +############################################################################### # Setup and run the PDHG algorithm operator = Aop f = L2NormSquared(b = noisy_data) g = ZeroFunction() + ## Compute operator Norm normK = operator.norm() + ## Primal & dual stepsizes sigma = 1 tau = 1/(sigma*normK**2) @@ -112,17 +117,17 @@ 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=False) +pdhg.run(1000, verbose=True) +############################################################################### # Setup and run the FISTA algorithm fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop) regularizer = ZeroFunction() -opt = {'memopt':True} -fista = FISTA(x_init=x_init , f=fidelity, g=regularizer, opt=opt) +fista = FISTA(x_init=x_init , f=fidelity, g=regularizer) fista.max_iteration = 1000 fista.update_objective_interval = 200 -fista.run(1000, verbose=False) +fista.run(1000, verbose=True) #%% Show results @@ -131,18 +136,22 @@ plt.suptitle('Reconstructions ', fontsize=16) plt.subplot(2,2,1) plt.imshow(cgls.get_output().as_array()) +plt.colorbar() plt.title('CGLS reconstruction') plt.subplot(2,2,2) plt.imshow(fista.get_output().as_array()) +plt.colorbar() plt.title('FISTA reconstruction') plt.subplot(2,2,3) plt.imshow(pdhg.get_output().as_array()) +plt.colorbar() plt.title('PDHG reconstruction') plt.subplot(2,2,4) plt.imshow(recon_cgls_astra) +plt.colorbar() plt.title('CGLS astra') diff1 = pdhg.get_output() - cgls.get_output() @@ -167,28 +176,14 @@ plt.title('Diff CLGS astra vs CGLS') plt.colorbar() +#%% +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])) +true_obj = (Aop.direct(cglsd.get_output())-noisy_data).squared_norm() +print('True objective {}'.format(true_obj)) - - - - - - - - - - - - -# -# -# -# -# -# -# -# diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py index 68fb649..1a96a2d 100644 --- a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py +++ b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py @@ -21,14 +21,15 @@ from ccpi.framework import AcquisitionGeometry from ccpi.optimisation.algorithms import FISTA -from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, L2NormSquared, FunctionOperatorComposition -from ccpi.astra.ops import AstraProjectorSimple +from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, \ + L2NormSquared, FunctionOperatorComposition +from ccpi.astra.operators import AstraProjectorSimple import numpy as np import matplotlib.pyplot as plt from ccpi.framework import TestData import os, sys - +from ccpi.optimisation.funcs import Norm2sq # Load Data loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) @@ -117,16 +118,15 @@ plt.show() # demonstrated in the rest of this file. In general all methods need an initial # guess and some algorithm options to be set: -f = FunctionOperatorComposition(0.5 * L2NormSquared(b=sin), Aop) +f = FunctionOperatorComposition(L2NormSquared(b=sin), Aop) +#f = Norm2sq(Aop, sin, c=0.5, memopt=True) g = ZeroFunction() x_init = ig.allocate() -opt = {'memopt':True} - -fista = FISTA(x_init=x_init, f=f, g=g, opt=opt) -fista.max_iteration = 100 -fista.update_objective_interval = 1 -fista.run(100, verbose=True) +fista = FISTA(x_init=x_init, f=f, g=g) +fista.max_iteration = 1000 +fista.update_objective_interval = 100 +fista.run(1000, verbose=True) plt.figure() plt.imshow(fista.get_output().as_array()) @@ -134,18 +134,12 @@ plt.title('FISTA unconstrained') plt.colorbar() plt.show() -plt.figure() -plt.semilogy(fista.objective) -plt.title('FISTA objective') -plt.show() - -#%% # Run FISTA for least squares with lower/upper bound -fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1), opt=opt) -fista0.max_iteration = 100 -fista0.update_objective_interval = 1 -fista0.run(100, verbose=True) +fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1)) +fista0.max_iteration = 1000 +fista0.update_objective_interval = 100 +fista0.run(1000, verbose=True) plt.figure() plt.imshow(fista0.get_output().as_array()) @@ -153,10 +147,10 @@ plt.title('FISTA constrained in [0,1]') plt.colorbar() plt.show() -plt.figure() -plt.semilogy(fista0.objective) -plt.title('FISTA constrained in [0,1]') -plt.show() +#plt.figure() +#plt.semilogy(fista0.objective) +#plt.title('FISTA constrained in [0,1]') +#plt.show() #%% Check with CVX solution @@ -178,10 +172,10 @@ if cvx_not_installable: vol_geom = astra.create_vol_geom(N, N) proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) - proj_id = astra.create_projector('strip', proj_geom, vol_geom) + proj_id = astra.create_projector('linear', proj_geom, vol_geom) matrix_id = astra.projector.matrix(proj_id) - + ProjMat = astra.matrix.get(matrix_id) tmp = sin.as_array().ravel() @@ -216,4 +210,6 @@ if cvx_not_installable: plt.legend() plt.title('Middle Line Profiles') plt.show() -
\ No newline at end of file + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (FISTA) {} '.format(fista0.loss[1]))
\ No newline at end of file diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py index 41219f4..6007990 100644 --- a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py +++ b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py @@ -21,7 +21,7 @@ """ -Tikhonov for Poisson denoising using FISTA algorithm: +"Tikhonov regularization" for Poisson denoising using FISTA algorithm: Problem: min_x, x>0 \alpha * ||\nabla x||_{2}^{2} + \int x - g * log(x) @@ -52,14 +52,14 @@ from skimage.util import random_noise loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) # Load Data -N = 150 -M = 150 +N = 100 +M = 100 data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) ig = data.geometry ag = ig -# Create Noisy data. Add Gaussian noise +# Create Noisy data with Poisson noise n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10) noisy_data = ImageData(n1) @@ -75,7 +75,6 @@ plt.title('Noisy Data') plt.colorbar() plt.show() -#%% # Regularisation Parameter alpha = 10 @@ -114,8 +113,7 @@ fid.proximal = KL_Prox_PosCone reg = FunctionOperatorComposition(alpha * L2NormSquared(), operator) x_init = ig.allocate() -opt = {'memopt':True} -fista = FISTA(x_init=x_init , f=reg, g=fid, opt=opt) +fista = FISTA(x_init=x_init , f=reg, g=fid) fista.max_iteration = 2000 fista.update_objective_interval = 500 fista.run(2000, verbose=True) @@ -196,7 +194,6 @@ if cvx_not_installable: plt.title('Middle Line Profiles') plt.show() - #TODO what is the output of fista.objective, fista.loss print('Primal Objective (CVX) {} '.format(obj.value)) print('Primal Objective (FISTA) {} '.format(fista.loss[1])) diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py index bc0081f..422deb9 100644 --- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py +++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py @@ -85,15 +85,8 @@ data = ImageData(data/data.max()) ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) ag = ig -# Create noisy data. Add Gaussian noise -scale = 0.5 -eta = 0 -n1 = scale * np.random.poisson(eta + sin.as_array()/scale) - -noisy_data = AcquisitionData(n1, ag) #Create Acquisition Data and apply poisson noise - detectors = N angles = np.linspace(0, np.pi, N) |