FISTA and PDHG ex

2 files changed, 402 insertions, 0 deletions

diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
new file mode 100644
index 0000000..5b1bb16
--- /dev/null
+++ b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
@@ -0,0 +1,184 @@
+#========================================================================
+# 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.
+#
+#=========================================================================
+
+""" 
+
+Tikhonov for Poisson denoising using FISTA algorithm:
+
+Problem:     min_x, x>0  \alpha * ||\nabla x||_{2} + \int x - g * log(x) 
+
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator  
+             
+             g: Noisy Data with Poisson Noise
+               
+                       
+"""
+
+from ccpi.framework import ImageData
+
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import FISTA
+
+from ccpi.optimisation.operators import Gradient, BlockOperator, Identity
+from ccpi.optimisation.functions import KullbackLeibler, IndicatorBox, BlockFunction, \
+                      L2NormSquared, IndicatorBox, FunctionOperatorComposition
+
+from ccpi.framework import TestData
+import os, sys
+from skimage.util import random_noise
+
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+
+# Load Data                      
+N = 50
+M = 50
+data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
+
+ig = data.geometry
+ag = ig
+
+# Create Noisy data. Add Gaussian noise
+n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10)
+noisy_data = ImageData(n1)
+
+# 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()
+
+#%%
+
+# Regularisation Parameter
+alpha = 20
+
+# Setup and run the FISTA algorithm
+op1 = Gradient(ig)
+op2 = BlockOperator(Identity(ig), Identity(ig), shape=(2,1))
+
+tmp_function = BlockFunction( KullbackLeibler(noisy_data), IndicatorBox(lower=0) )
+
+fid = tmp
+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.max_iteration = 2000
+fista.update_objective_interval = 500
+fista.run(2000, verbose=True)
+
+#%%
+# Show results
+plt.figure(figsize=(15,15))
+plt.subplot(3,1,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(3,1,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(3,1,3)
+plt.imshow(fista.get_output().as_array())
+plt.title('Reconstruction')
+plt.colorbar()
+plt.show()
+
+plt.plot(np.linspace(0,N,M), data.as_array()[int(N/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,N,M), fista.get_output().as_array()[int(N/2),:], label = 'Reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
+
+#%% Check with CVX solution
+
+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    
+    u1 = Variable(ig.shape)
+    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(u1), DY.matrix() * vec(u1)]), 2, axis = 0))    
+    fidelity = sum(kl_div(noisy_data.as_array(), u1))  
+    
+    constraints = [q>=fidelity, u1>=0]    
+    
+    solver = SCS
+    obj =  Minimize( regulariser +  q)
+    prob = Problem(obj, constraints)
+    result = prob.solve(verbose = True, solver = solver)
+    
+    diff_cvx = numpy.abs( fista.get_output().as_array() - u1.value )
+        
+    plt.figure(figsize=(15,15))
+    plt.subplot(3,1,1)
+    plt.imshow(fista.get_output().as_array())
+    plt.title('FISTA solution')
+    plt.colorbar()
+    plt.subplot(3,1,2)
+    plt.imshow(u1.value)
+    plt.title('CVX solution')
+    plt.colorbar()
+    plt.subplot(3,1,3)
+    plt.imshow(diff_cvx)
+    plt.title('Difference')
+    plt.colorbar()
+    plt.show()    
+    
+    plt.plot(np.linspace(0,N,M), fista.get_output().as_array()[int(N/2),:], label = 'FISTA')
+    plt.plot(np.linspace(0,N,M), u1.value[int(N/2),:], label = 'CVX')
+    plt.legend()
+    plt.title('Middle Line Profiles')
+    plt.show()
+            
+    print('Primal Objective (CVX) {} '.format(obj.value))
+    print('Primal Objective (FISTA) {} '.format(fista.objective[-1][0]))
+
+
+
+
+
+
+
diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py
new file mode 100644
index 0000000..830ea00
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py
@@ -0,0 +1,218 @@
+#========================================================================
+# 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.
+#
+#=========================================================================
+
+from ccpi.framework import AcquisitionGeometry, AcquisitionData
+
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
+from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
+                      MixedL21Norm, BlockFunction, IndicatorBox
+
+from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.framework import TestData
+import os, sys
+
+""" 
+
+Total Variation Denoising using PDHG algorithm:
+
+
+Problem:     min_x, x>0  \alpha * ||\nabla x||_{2,1} + int A x -g log(Ax + \eta)
+
+             \nabla: Gradient operator 
+             
+             A: Projection Matrix
+             g: Noisy sinogram corrupted with Poisson Noise
+             
+             \eta: Background Noise
+             \alpha: Regularization parameter
+ 
+"""
+
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+
+# Load Data                      
+N = 50
+M = 50
+data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
+
+ig = data.geometry
+ag = ig
+
+#Create Acquisition Data and apply poisson noise
+
+detectors = N
+angles = np.linspace(0, np.pi, N)
+
+ag = AcquisitionGeometry('parallel','2D',angles, detectors)
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
+sin = Aop.direct(data)
+
+# Create noisy data. Apply Poisson noise
+scale = 0.5
+eta = 0 
+n1 = scale * np.random.poisson(eta + sin.as_array()/scale)
+
+noisy_data = AcquisitionData(n1, ag)
+
+# 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()
+
+# Regularisation Parameter
+alpha = 2
+
+# Create operators
+op1 = Gradient(ig)
+op2 = Aop
+
+# Create BlockOperator
+operator = BlockOperator(op1,op2, shape=(2,1) ) 
+
+# Create functions
+      
+f1 = alpha * MixedL21Norm() 
+f2 = KullbackLeibler(noisy_data)    
+f = BlockFunction(f1, f2)  
+                                      
+g = IndicatorBox(lower=0)
+
+
+# Compute operator Norm
+normK = operator.norm()
+    
+# Primal & dual stepsizes
+sigma = 1
+tau = 1/(sigma*normK**2)
+
+
+# Setup and run the PDHG algorithm
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+pdhg.max_iteration = 3000
+pdhg.update_objective_interval = 500
+pdhg.run(3000, verbose = True)
+
+plt.figure(figsize=(15,15))
+plt.subplot(3,1,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(3,1,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(3,1,3)
+plt.imshow(pdhg.get_output().as_array())
+plt.title('TV Reconstruction')
+plt.colorbar()
+plt.show()
+## 
+plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
+
+
+##%% Check with CVX solution
+#
+#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 = alpha * sum(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)
+#    
+#    tmp = noisy_data.as_array().ravel()
+#    
+#    fidelity = sum(kl_div(tmp, ProjMat * u))    
+#    
+#    constraints = [q>=fidelity, u>=0]   
+#    solver = SCS
+#    obj =  Minimize( regulariser +  q)
+#    prob = Problem(obj, constraints)
+#    result = prob.solve(verbose = True, solver = solver)    
+#             
+#    diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N)))
+#           
+#    plt.figure(figsize=(15,15))
+#    plt.subplot(3,1,1)
+#    plt.imshow(pdhg.get_output().as_array())
+#    plt.title('PDHG solution')
+#    plt.colorbar()
+#    plt.subplot(3,1,2)
+#    plt.imshow(np.reshape(u.value, (N, N)))
+#    plt.title('CVX solution')
+#    plt.colorbar()        
+#    plt.subplot(3,1,3)
+#    plt.imshow(diff_cvx)
+#    plt.title('Difference')
+#    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), np.reshape(u.value, (N, N))[int(N/2),:], label = 'CVX')
+#    plt.legend()
+#    plt.title('Middle Line Profiles')
+#    plt.show()
+#            
+#    print('Primal Objective (CVX) {} '.format(obj.value))
+#    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file