diff options
| -rw-r--r-- | Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py | 117 | ||||
| -rw-r--r-- | Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py | 198 | ||||
| -rw-r--r-- | Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py | 127 | ||||
| -rw-r--r-- | Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py | 152 |
4 files changed, 594 insertions, 0 deletions
diff --git a/Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py b/Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py new file mode 100644 index 0000000..eb62761 --- /dev/null +++ b/Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py @@ -0,0 +1,117 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018-2019 Evangelos Papoutsellis and Edoardo Pasca + +# 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 + +# 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 FISTA & PDHG classes + + +Problem: min_x alpha * ||\grad x ||^{2}_{2} + || x - g ||_{1} + + A: Projection operator + g: Sinogram + +""" + +from ccpi.framework import ImageData, ImageGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import FISTA, PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient, Identity +from ccpi.optimisation.functions import L2NormSquared, L1Norm, \ + FunctionOperatorComposition, BlockFunction, ZeroFunction + +from skimage.util import random_noise + + +# Create Ground truth and noisy data + +N = 100 + +data = np.zeros((N,N)) +data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 5 + + +# Setup and run the FISTA algorithm +operator = Gradient(ig) +fidelity = L1Norm(b=noisy_data) +regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator) + +x_init = ig.allocate() +opt = {'memopt':True} +fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt) +fista.max_iteration = 2000 +fista.update_objective_interval = 50 +fista.run(2000, verbose=False) + +# Setup and run the PDHG algorithm +op1 = Gradient(ig) +op2 = Identity(ig, ag) + +operator = BlockOperator(op1, op2, shape=(2,1) ) +f = BlockFunction(alpha * L2NormSquared(), fidelity) +g = ZeroFunction() + +normK = operator.norm() + +sigma = 1 +tau = 1/(sigma*normK**2) + +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000, verbose=False) + +#%% +# Show results + +plt.figure(figsize=(15,15)) + +plt.subplot(1,2,1) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') + +plt.subplot(1,2,2) +plt.imshow(fista.get_output().as_array()) +plt.title('FISTA reconstruction') + +plt.show() + +diff1 = pdhg.get_output() - fista.get_output() + +plt.imshow(diff1.abs().as_array()) +plt.title('Diff PDHG vs CGLS') +plt.colorbar() +plt.show() + + diff --git a/Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py b/Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py new file mode 100644 index 0000000..c877018 --- /dev/null +++ b/Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py @@ -0,0 +1,198 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018-2019 Evangelos Papoutsellis and Edoardo Pasca + +# 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 + +# 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 Least Squares minimization problem using FISTA, PDHG, CGLS classes +and Astra Built-in CGLS + +Problem: min_x || A x - g ||_{2}^{2} + + A: Projection operator + g: Sinogram + +""" + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA + +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition +from ccpi.astra.ops import AstraProjectorSimple +import astra + +# 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) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D', angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +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) + +# Create a sinogram from a phantom +sinogram_id, sinogram = astra.create_sino(x, proj_id) + +# Create a data object for the reconstruction +rec_id = astra.data2d.create('-vol', vol_geom) + +cgls_astra = astra.astra_dict('CGLS') +cgls_astra['ReconstructionDataId'] = rec_id +cgls_astra['ProjectionDataId'] = sinogram_id +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) + +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) + +# 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) + +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) + +# 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.max_iteration = 1000 +fista.update_objective_interval = 200 +fista.run(1000, verbose=True) + +#%% Show results + +plt.figure(figsize=(15,15)) +plt.suptitle('Reconstructions ') + +plt.subplot(2,2,1) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') + +plt.subplot(2,2,2) +plt.imshow(fista.get_output().as_array()) +plt.title('FISTA reconstruction') + +plt.subplot(2,2,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') + +plt.subplot(2,2,4) +plt.imshow(recon_cgls_astra) +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() + +print( diff1.squared_norm()) +print( diff2.squared_norm()) +print( diff3.squared_norm()) + +plt.figure(figsize=(15,15)) + +plt.subplot(3,1,1) +plt.imshow(diff1.abs().as_array()) +plt.title('Diff PDHG vs CGLS') +plt.colorbar() + +plt.subplot(3,1,2) +plt.imshow(diff2.abs().as_array()) +plt.title('Diff FISTA vs CGLS') +plt.colorbar() + +plt.subplot(3,1,3) +plt.imshow(diff3.abs().as_array()) +plt.title('Diff CLGS astra vs CGLS') +plt.colorbar() + + + + + + + + + + + + + + + + + + + +# +# +# +# +# +# +# +# diff --git a/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py b/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py new file mode 100644 index 0000000..3155654 --- /dev/null +++ b/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py @@ -0,0 +1,127 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018-2019 Evangelos Papoutsellis and Edoardo Pasca + +# 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 + +# 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS + +from ccpi.optimisation.operators import Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition +from skimage.util import random_noise +from ccpi.astra.ops import AstraProjectorSimple + +#%% + +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) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +noisy_data = sin + +x_init = ig.allocate() + +## Setup and run the CGLS algorithm +cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) +cgls.max_iteration = 500 +cgls.update_objective_interval = 50 +cgls.run(500, verbose=True) + +# Create BlockOperator +operator = Aop +f = 0.5 * L2NormSquared(b = noisy_data) +g = ZeroFunction() + +## Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes +sigma = 0.1 +tau = 1/(sigma*normK**2) +# +# +## Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + +#%% + +diff = pdhg.get_output() - cgls.get_output() +print( diff.norm()) +# +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(diff.abs().as_array()) +plt.title('Difference reconstruction') +plt.colorbar() +plt.show() + + + + + + + + + + + + + + + + + + + + + + +# +# +# +# +# +# +# +# diff --git a/Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py b/Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py new file mode 100644 index 0000000..984fca4 --- /dev/null +++ b/Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py @@ -0,0 +1,152 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018-2019 Evangelos Papoutsellis and Edoardo Pasca + +# 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 + +# 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 Tikhonov with PDHG, CGLS classes + + +Problem: min_x alpha * ||\grad x ||^{2}_{2} + || A x - g ||_{2}^{2} + + A: Projection operator + g: Sinogram + +""" + + +from ccpi.framework import ImageData, ImageGeometry, \ + AcquisitionGeometry, BlockDataContainer, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS +from ccpi.optimisation.operators import BlockOperator, Gradient + +from ccpi.optimisation.functions import ZeroFunction, BlockFunction, L2NormSquared +from ccpi.astra.ops import AstraProjectorSimple + +# 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) + +detectors = N +angles = np.linspace(0, np.pi, N, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D', angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,3,ig.shape)) + +# Setup and run the CGLS algorithm + +alpha = 50 +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) + +#%% +#Setup and run the PDHG algorithm + +# Create BlockOperator +op_PDHG = BlockOperator(Grad, Aop, shape=(2,1) ) + +# Create functions + +f1 = 0.5 * alpha**2 * L2NormSquared() +f2 = 0.5 * L2NormSquared(b = noisy_data) +f = BlockFunction(f1, f2) +g = ZeroFunction() + +## Compute operator Norm +normK = op_PDHG.norm() + +## Primal & dual stepsizes +sigma = 10 +tau = 1/(sigma*normK**2) + +pdhg = PDHG(f=f,g=g,operator=op_PDHG, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 1000 +pdhg.update_objective_interval = 200 +pdhg.run(1000) + +#%% +# Show results + +plt.figure(figsize=(15,15)) + +plt.subplot(1,2,1) +plt.imshow(cgls.get_output().as_array()) +plt.title('CGLS reconstruction') + +plt.subplot(1,2,2) +plt.imshow(pdhg.get_output().as_array()) +plt.title('PDHG reconstruction') + +plt.show() + +diff1 = pdhg.get_output() - cgls.get_output() + +plt.imshow(diff1.abs().as_array()) +plt.title('Diff PDHG vs CGLS') +plt.colorbar() +plt.show() + +#%% + +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') +plt.plot(np.linspace(0,N,N), cgls.get_output().as_array()[int(N/2),:], label = 'CGLS') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + + + + + + + + +# +# +# +# +# +# +# +# |