add Algs comparisons

3 files changed, 471 insertions, 0 deletions

diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
new file mode 100644
index 0000000..0875c20
--- /dev/null
+++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
@@ -0,0 +1,194 @@
+#========================================================================
+# 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 FISTA & PDHG 
+                   & CGLS  & Astra Built-in algorithms for Least Squares
+
+
+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)
+
+device = input('Available device: GPU==1 / CPU==0 ')
+ag = AcquisitionGeometry('parallel','2D', angles, detectors)
+if device=='1':
+    dev = 'gpu'
+else:
+    dev = 'cpu'
+    
+Aop = AstraProjectorSimple(ig, ag, dev)
+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, 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)
+
+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)
+
+# 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=False)
+
+#%% Show results
+
+plt.figure(figsize=(10,10))
+plt.suptitle('Reconstructions ', fontsize=16)
+
+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()
+
+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/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py
new file mode 100644
index 0000000..942d328
--- /dev/null
+++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py
@@ -0,0 +1,153 @@
+#========================================================================
+# 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 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       
+
+# 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)
+device = input('Available device: GPU==1 / CPU==0 ')
+ag = AcquisitionGeometry('parallel','2D', angles, detectors)
+if device=='1':
+    dev = 'gpu'
+else:
+    dev = '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,verbose=False)
+
+
+#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, verbose=False)
+
+
+#%%
+# Show results
+plt.figure(figsize=(10,10))
+
+plt.subplot(2,1,1)
+plt.imshow(cgls.get_output().as_array())
+plt.title('CGLS reconstruction')
+
+plt.subplot(2,1,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()
+            
+
+
+
+
+
+
+
+
+#
+#
+#
+#
+#
+#
+#
+#
diff --git a/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py b/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py
new file mode 100644
index 0000000..c24ebac
--- /dev/null
+++ b/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py
@@ -0,0 +1,124 @@
+#========================================================================
+# 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 FISTA & PDHG algorithms
+
+
+Problem:     min_x alpha * ||\grad x ||^{2}_{2} + || x - g ||_{1}
+
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator
+             
+             g: Noisy Data with Salt & Pepper Noise
+
+"""
+
+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 = 200
+pdhg.run(2000, verbose=False)
+###############################################################################
+
+# Show results
+
+plt.figure(figsize=(10,10))
+
+plt.subplot(2,1,1)
+plt.imshow(pdhg.get_output().as_array())
+plt.title('PDHG reconstruction')
+
+plt.subplot(2,1,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 FISTA')
+plt.colorbar()
+plt.show()
+
+