diff options
Diffstat (limited to 'Wrappers/Python')
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py | 2 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py | 53 | ||||
| -rw-r--r-- | Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py | 153 | ||||
| -rw-r--r-- | Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py | 214 | ||||
| -rw-r--r-- | Wrappers/Python/demos/CGLS_examples/LinearSystem.py | 82 | ||||
| -rw-r--r-- | Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py | 146 | ||||
| -rwxr-xr-x | Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat | bin | 5583 -> 0 bytes |
7 files changed, 486 insertions, 164 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index c62d0ea..bf851d7 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -84,6 +84,8 @@ class Algorithm(object): calling this method triggers update and update_objective ''' if self.should_stop(): + self.update_objective() + print(self.verbose_output()) raise StopIteration() else: time0 = time.time() diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py index 28b19f6..6cf4f06 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -34,19 +34,31 @@ class CGLS(Algorithm): x_init: initial guess operator: operator for forward/backward projections data: data to operate on + tolerance: tolerance to stop algorithm + + Reference: + https://web.stanford.edu/group/SOL/software/cgls/ + ''' + + + def __init__(self, **kwargs): + super(CGLS, self).__init__() self.x = kwargs.get('x_init', None) self.operator = kwargs.get('operator', None) self.data = kwargs.get('data', None) + self.tolerance = kwargs.get('tolerance', None) if self.x is not None and self.operator is not None and \ self.data is not None: - print ("Calling from creator") self.set_up(x_init =kwargs['x_init'], operator=kwargs['operator'], data =kwargs['data']) + if self.tolerance is None: + self.tolerance = 1e-6 + def set_up(self, x_init, operator , data ): self.x = x_init @@ -55,10 +67,19 @@ class CGLS(Algorithm): self.p = self.s self.norms0 = self.s.norm() + + ## + self.norms = self.s.norm() + ## + + self.gamma = self.norms0**2 self.normx = self.x.norm() self.xmax = self.normx - self.configured = True + + n = Norm2Sq(self.operator, self.data) + self.loss.append(n(self.x)) + self.configured = True # def set_up(self, x_init, operator , data ): # @@ -80,15 +101,15 @@ class CGLS(Algorithm): # self.loss.append(n(x_init)) # self.configured = True - def update(self): - self.update_new() + #def update(self): + # self.update_new() - def update_new(self): + def update(self): self.q = self.operator.direct(self.p) delta = self.q.squared_norm() alpha = self.gamma/delta - + self.x += alpha * self.p self.r -= alpha * self.q @@ -102,10 +123,7 @@ class CGLS(Algorithm): self.normx = self.x.norm() self.xmax = numpy.maximum(self.xmax, self.normx) - - - if self.gamma<=1e-6: - raise StopIteration() + # def update_new(self): # # Ad = self.operator.direct(self.d) @@ -141,7 +159,20 @@ class CGLS(Algorithm): raise StopIteration() self.loss.append(a) -# def should_stop(self): + def should_stop(self): + + 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() + return True + + + #raise StopIteration() + + # if self.iteration > 0: # x = self.get_last_objective() # a = x > 0 diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py deleted file mode 100644 index 63a2254..0000000 --- a/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py +++ /dev/null @@ -1,153 +0,0 @@ -#======================================================================== -# Copyright 2019 Science Technology Facilities Council -# Copyright 2019 University of Manchester -# -# This work is part of the Core Imaging Library developed by Science Technology -# Facilities Council and University of Manchester -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0.txt -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -#========================================================================= - -""" -Compare solutions of PDHG & "Block CGLS" algorithms for - - -Problem: min_x alpha * ||\grad x ||^{2}_{2} + || A x - g ||_{2}^{2} - - - A: Projection operator - g: Sinogram - -""" - - -from ccpi.framework import AcquisitionGeometry, BlockDataContainer, AcquisitionData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import CGLS -from ccpi.optimisation.operators import BlockOperator, Gradient - -from ccpi.framework import TestData -import os, sys -from ccpi.astra.ops import AstraProjectorSimple - -# Load Data -loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) -N = 64 -M = 64 -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) - -ag = AcquisitionGeometry('parallel','2D', angles, detectors) - -device = input('Available device: GPU==1 / CPU==0 ') -if device=='1': - dev = 'gpu' -else: - dev = 'cpu' - -Aop = AstraProjectorSimple(ig, ag, dev) -sin = Aop.direct(data) - -noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,3,ig.shape)) - -# Show Ground Truth and Noisy Data -plt.figure(figsize=(10,10)) -plt.subplot(2,1,1) -plt.imshow(data.as_array()) -plt.title('Ground Truth') -plt.colorbar() -plt.subplot(2,1,2) -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Data') -plt.colorbar() -plt.show() - -# Setup and run the CGLS algorithm -alpha = 5 -Grad = Gradient(ig) -# -## Form Tikhonov as a Block CGLS structure -op_CGLS = BlockOperator( Aop, alpha * Grad, shape=(2,1)) -block_data = BlockDataContainer(noisy_data, Grad.range_geometry().allocate()) -# -x_init = ig.allocate() -cgls = CGLS(x_init=x_init, operator=op_CGLS, data=block_data) -cgls.max_iteration = 1000 -cgls.update_objective_interval = 200 -cgls.run(1000,verbose=True) - -# Show results -plt.figure(figsize=(5,5)) -plt.imshow(cgls.get_output().as_array()) -plt.title('CGLS reconstruction') -plt.colorbar() -plt.show() - -#%% -from ccpi.optimisation.operators import SparseFiniteDiff -import astra -import numpy - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - - ##Construct problem - u = Variable(N*M) - #q = Variable() - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - - # create matrix representation for Astra operator - - 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) - - matrix_id = astra.projector.matrix(proj_id) - - ProjMat = astra.matrix.get(matrix_id) - - fidelity = sum_squares( ProjMat * u - noisy_data.as_array().ravel()) - - solver = SCS - obj = Minimize( regulariser + fidelity) - prob = Problem(obj, constraints) - result = prob.solve(verbose = True, solver = solver) - - - - - - - - diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py new file mode 100644 index 0000000..abcb26c --- /dev/null +++ b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py @@ -0,0 +1,214 @@ +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= + +""" +Conjugate Gradient for (Regularized) Least Squares for Tomography + + +Problem: min_u alpha * || L x ||^{2}_{2} + || A u - g ||_{2}^{2} + + A: Projection operator + g: Sinogram + L: Identity or Gradient Operator + +""" + + +from ccpi.framework import ImageGeometry, ImageData, \ + AcquisitionGeometry, BlockDataContainer, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import CGLS +from ccpi.optimisation.operators import BlockOperator, Gradient, Identity + +import tomophantom +from tomophantom import TomoP2D +from ccpi.astra.operators import AstraProjectorSimple +import os + + +# Load Shepp-Logan phantom +model = 1 # select a model number from the library +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) + +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) + +ag = AcquisitionGeometry('parallel','2D', angles, detectors) + +device = input('Available device: GPU==1 / CPU==0 ') + +if device =='1': + dev = 'gpu' +else: + dev = 'cpu' + +Aop = AstraProjectorSimple(ig, ag, dev) +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)) +plt.subplot(2,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(2,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.show() + + +# Setup and run the simple CGLS algorithm +x_init = ig.allocate() + +cgls1 = CGLS(x_init = x_init, operator = Aop, data = noisy_data) +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() +alpha1 = 50 +op1 = Identity(ig) + +block_op1 = BlockOperator( Aop, alpha1 * op1, shape=(2,1)) +block_data1 = BlockDataContainer(noisy_data, op1.range_geometry().allocate()) + +cgls2 = CGLS(x_init = x_init, operator = block_op1, data = block_data1) +cgls2.max_iteration = 200 +cgls2.update_objective_interval = 10 +cgls2.run(200, verbose=True) + +# Setup and run the regularised CGLS algorithm (Tikhonov with Gradient) + +x_init = ig.allocate() +alpha2 = 25 +op2 = Gradient(ig) + +block_op2 = BlockOperator( Aop, alpha2 * op2, shape=(2,1)) +block_data2 = BlockDataContainer(noisy_data, op2.range_geometry().allocate()) + +cgls3 = CGLS(x_init=x_init, operator = block_op2, data = block_data2) +cgls3.max_iteration = 200 +cgls3.update_objective_interval = 5 +cgls3.run(200, verbose = True) + +# Show results +plt.figure(figsize=(8,8)) + +plt.subplot(2,2,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') + +plt.subplot(2,2,2) +plt.imshow(cgls1.get_output().as_array()) +plt.title('GGLS') + +plt.subplot(2,2,3) +plt.imshow(cgls2.get_output().as_array()) +plt.title('Regularised GGLS L = {} * Identity'.format(alpha1)) + +plt.subplot(2,2,4) +plt.imshow(cgls3.get_output().as_array()) +plt.title('Regularised GGLS L = {} * Gradient'.format(alpha2)) + +plt.show() + + + +print('Compare CVX vs Regularised CG with L = Gradient') + +from ccpi.optimisation.operators import SparseFiniteDiff +import astra +import numpy + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u = Variable(N*N) + #q = Variable() + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + regulariser = alpha2**2 * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + + # create matrix representation for Astra operator + + 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) + + matrix_id = astra.projector.matrix(proj_id) + + ProjMat = astra.matrix.get(matrix_id) + + fidelity = sum_squares( ProjMat * u - noisy_data.as_array().ravel()) + + solver = MOSEK + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = solver) + + +plt.figure(figsize=(10,20)) + +plt.subplot(1,3,1) +plt.imshow(np.reshape(u.value, (N, N))) +plt.colorbar() + +plt.subplot(1,3,2) +plt.imshow(cgls3.get_output().as_array()) +plt.colorbar() + +plt.subplot(1,3,3) +plt.imshow(np.abs(cgls3.get_output().as_array() - np.reshape(u.value, (N, N)) )) +plt.colorbar() + +plt.show() + +print('Primal Objective (CVX) {} '.format(obj.value)) +print('Primal Objective (CGLS) {} '.format(cgls3.objective[-1])) + diff --git a/Wrappers/Python/demos/CGLS_examples/LinearSystem.py b/Wrappers/Python/demos/CGLS_examples/LinearSystem.py new file mode 100644 index 0000000..cc398cb --- /dev/null +++ b/Wrappers/Python/demos/CGLS_examples/LinearSystem.py @@ -0,0 +1,82 @@ +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= + +import numpy +from ccpi.optimisation.operators import * +from ccpi.optimisation.algorithms import * +from ccpi.optimisation.functions import * +from ccpi.framework import * + +# Problem dimension. +m = 200 +n = 200 + +numpy.random.seed(10) + +# Create matrix A and data b +Amat = numpy.asarray( numpy.random.randn(m, n), dtype = numpy.float32) +bmat = numpy.asarray( numpy.random.randn(m), dtype=numpy.float32) + +# Geometries for data and solution +vgb = VectorGeometry(m) +vgx = VectorGeometry(n) + +b = vgb.allocate(0, dtype=numpy.float32) +b.fill(bmat) +A = LinearOperatorMatrix(Amat) + +# Setup and run CGLS +x_init = vgx.allocate() + +cgls = CGLS(x_init = x_init, operator = A, data = b) +cgls.max_iteration = 2000 +cgls.update_objective_interval = 200 +cgls.run(2000, verbose = True) + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + +# Compare with CVX +x = Variable(n) +obj = Minimize(sum_squares(A.A*x - b.as_array())) +prob = Problem(obj) + +# choose solver +if 'MOSEK' in installed_solvers(): + solver = MOSEK +else: + solver = SCS + +result = prob.solve(solver = MOSEK) + + +diff_sol = x.value - cgls.get_output().as_array() + +print('Error |CVX - CGLS| = {}'.format(numpy.sum(numpy.abs(diff_sol)))) +print('CVX objective = {}'.format(obj.value)) +print('CGLS objective = {}'.format(cgls.objective[-1])) + + + + diff --git a/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py b/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py new file mode 100644 index 0000000..c124fa1 --- /dev/null +++ b/Wrappers/Python/demos/CGLS_examples/Regularised_CGLS_Denoising.py @@ -0,0 +1,146 @@ +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0.txt +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +#========================================================================= + +""" +Conjugate Gradient for (Regularized) Least Squares for Tomography + + +Problem: min_u alpha * || L x ||^{2}_{2} + || A u - g ||_{2}^{2} + + A: Identity operator + g: Sinogram + L: Identity or Gradient Operator + +""" + + +from ccpi.framework import ImageGeometry, ImageData, \ + AcquisitionGeometry, BlockDataContainer, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import CGLS +from ccpi.optimisation.operators import BlockOperator, Gradient, Identity +from ccpi.framework import TestData + +import os, sys + +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) +data = loader.load(TestData.SHAPES) +ig = data.geometry +ag = ig + +noisy_data = ImageData(TestData.random_noise(data.as_array(), mode = 'gaussian', seed = 1)) +#noisy_data = ImageData(data.as_array()) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(10,10)) +plt.subplot(2,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(2,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.show() + +# Setup and run the regularised CGLS algorithm (Tikhonov with Gradient) +x_init = ig.allocate() +alpha = 2 +op = Gradient(ig) + +block_op = BlockOperator( Identity(ig), alpha * op, shape=(2,1)) +block_data = BlockDataContainer(noisy_data, op.range_geometry().allocate()) + +cgls = CGLS(x_init=x_init, operator = block_op, data = block_data) +cgls.max_iteration = 200 +cgls.update_objective_interval = 5 +cgls.run(200, verbose = True) + +# Show results +plt.figure(figsize=(20,10)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy') +plt.subplot(3,1,3) +plt.imshow(cgls.get_output().as_array()) +plt.title('Regularised GGLS with Gradient') +plt.show() + +#%% + +print('Compare CVX vs Regularised CG with L = Gradient') + +from ccpi.optimisation.operators import SparseFiniteDiff + + +try: + from cvxpy import * + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + + +if cvx_not_installable: + + ##Construct problem + u = Variable(ig.shape) + #q = Variable() + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + regulariser = alpha**2 * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + + fidelity = sum_squares( u - noisy_data.as_array()) + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + + obj = Minimize( regulariser + fidelity) + prob = Problem(obj) + result = prob.solve(verbose = True, solver = MOSEK) + + +plt.figure(figsize=(20,20)) +plt.subplot(3,1,1) +plt.imshow(np.reshape(u.value, ig.shape)) +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(cgls.get_output().as_array()) +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(np.abs(cgls.get_output().as_array() - np.reshape(u.value, ig.shape) )) +plt.colorbar() +plt.show() + +print('Primal Objective (CVX) {} '.format(obj.value)) +print('Primal Objective (CGLS) {} '.format(cgls.objective[-1])) + diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat b/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat Binary files differdeleted file mode 100755 index c465bbe..0000000 --- a/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat +++ /dev/null |