demos and check 2D_time denosigin

3 files changed, 276 insertions, 100 deletions

diff --git a/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py
index b6d7725..72d0670 100644
--- a/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py
+++ b/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py
@@ -29,7 +29,7 @@ from ccpi.optimisation.algorithms import PDHG
 
 from ccpi.optimisation.operators import BlockOperator, Gradient
 from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
+                      MixedL21Norm, BlockFunction, IndicatorBox
 
 from ccpi.astra.ops import AstraProjectorSimple
 
@@ -67,10 +67,13 @@ Aop = AstraProjectorSimple(ig, ag, 'cpu')
 sin = Aop.direct(data)
 
 # Create noisy data. Apply Poisson noise
-scale = 0.5
-n1 = scale * np.random.poisson(sin.as_array()/scale)
+scale = 0.25
+eta = 0 #np.random.randint(0, sin.as_array().max()/2, ag.shape)
+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)
@@ -83,9 +86,10 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
 
+#%%
 
 # Regularisation Parameter
-alpha = 0.5
+alpha = 2
 
 # Create operators
 op1 = Gradient(ig)
@@ -100,20 +104,23 @@ f1 = alpha * MixedL21Norm()
 f2 = KullbackLeibler(noisy_data)    
 f = BlockFunction(f1, f2)  
                                       
-g = ZeroFunction()
+#g = ZeroFunction()
+g = IndicatorBox(lower=0)
 
 # Compute operator Norm
 normK = operator.norm()
     
 # Primal & dual stepsizes
-sigma = 10
+sigma = 2
 tau = 1/(sigma*normK**2)
+#sigma = 1/normK
+#tau = 1/normK
 
 # Setup and run the PDHG algorithm
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
 pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 200
-pdhg.run(2000)
+pdhg.update_objective_interval = 500
+pdhg.run(2000, verbose = True)
 
 plt.figure(figsize=(15,15))
 plt.subplot(3,1,1)
@@ -137,64 +144,11 @@ 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('F')
-    
-#    fidelity = sum( ProjMat * u - tmp * log(ProjMat * u + 1e-6)) 
-    #constraints = [q>= fidelity, u>=0]
-    constraints = []
-    
-    fidelity = sum(kl_div(tmp, ProjMat * u + 1e-6))
-#    fidelity =  kl_div(cp.multiply(alpha, W),
-#                  cp.multiply(alpha, W + cp.multiply(beta, P))) - \
-#                  cp.multiply(alpha, cp.multiply(beta, P))
-    
-    
-        
-    solver = SCS
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj, constraints)
-    result = prob.solve(verbose = True, solver = solver)    
-     
-
-###%% Check with CVX solution
+##%% Check with CVX solution
 #
 #from ccpi.optimisation.operators import SparseFiniteDiff
+#import astra
+#import numpy
 #
 #try:
 #    from cvxpy import *
@@ -204,48 +158,101 @@ if cvx_not_installable:
 #
 #
 #if cvx_not_installable:
+#    
 #
 #    ##Construct problem    
-#    u = Variable(ig.shape)
+#    u = Variable(N*N)
+#    #q = Variable()
 #    
 #    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
 #    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-#    
-#    # Define Total Variation as a regulariser
+#
 #    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-#    fidelity = pnorm( u - noisy_data.as_array(),1)
 #    
-#    # choose solver
-#    if 'MOSEK' in installed_solvers():
-#        solver = MOSEK
-#    else:
-#        solver = SCS  
+#    # 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('F')
+#    
+##    fidelity = sum( ProjMat * u - tmp * log(ProjMat * u + 1e-6)) 
+#    #constraints = [q>= fidelity, u>=0]
+#    constraints = []
+#    
+#    fidelity = sum(kl_div(tmp, ProjMat * u + 1e-6))
+##    fidelity =  kl_div(cp.multiply(alpha, W),
+##                  cp.multiply(alpha, W + cp.multiply(beta, P))) - \
+##                  cp.multiply(alpha, cp.multiply(beta, P))
+#    
+#    
 #        
+#    solver = SCS
 #    obj =  Minimize( regulariser +  fidelity)
-#    prob = Problem(obj)
-#    result = prob.solve(verbose = True, solver = solver)
+#    prob = Problem(obj, constraints)
+#    result = prob.solve(verbose = True, solver = solver)    
+#     
+#
+####%% 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    
+##    u = Variable(ig.shape)
+##    
+##    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+##    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+##    
+##    # Define Total Variation as a regulariser
+##    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+##    fidelity = pnorm( u - noisy_data.as_array(),1)
+##    
+##    # 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 = solver)
+##    
+##           
+#    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.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), u.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 (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file
+#    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
+#    plt.plot(np.linspace(0,N,N), u.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 (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file
diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py
new file mode 100644
index 0000000..045458a
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py
@@ -0,0 +1,169 @@
+# -*- 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
+
+from ccpi.optimisation.operators import BlockOperator, Gradient
+from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
+                      MixedL21Norm, BlockFunction
+
+from ccpi.astra.ops import AstraProjectorMC
+
+import os
+import tomophantom
+from tomophantom import TomoP2D
+
+# Create phantom for TV 2D dynamic tomography 
+
+model = 102  # note that the selected model is temporal (2D + time)
+N = 50 # set dimension of the phantom
+# one can specify an exact path to the parameters file
+# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat'
+path = os.path.dirname(tomophantom.__file__)
+path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
+#This will generate a N_size x N_size x Time frames phantom (2D + time)
+phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D)
+
+plt.close('all')
+plt.figure(1)
+plt.rcParams.update({'font.size': 21})
+plt.title('{}''{}'.format('2D+t phantom using model no.',model))
+for sl in range(0,np.shape(phantom_2Dt)[0]):
+    im = phantom_2Dt[sl,:,:]
+    plt.imshow(im, vmin=0, vmax=1)
+    plt.pause(.1)
+    plt.draw
+
+    
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0])
+data = ImageData(phantom_2Dt, geometry=ig)
+
+detectors = N
+angles = np.linspace(0,np.pi,N)
+
+ag = AcquisitionGeometry('parallel','2D', angles, detectors, channels = np.shape(phantom_2Dt)[0])
+Aop = AstraProjectorMC(ig, ag, 'gpu')
+sin = Aop.direct(data)
+
+scale = 2
+n1 = scale * np.random.poisson(sin.as_array()/scale)
+noisy_data = AcquisitionData(n1, ag)
+
+tindex = [3, 6, 10]
+
+fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10))
+plt.subplot(1,3,1)
+plt.imshow(noisy_data.as_array()[tindex[0],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[0]))
+plt.subplot(1,3,2)
+plt.imshow(noisy_data.as_array()[tindex[1],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[1]))
+plt.subplot(1,3,3)
+plt.imshow(noisy_data.as_array()[tindex[2],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[2]))
+
+fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
+                    wspace=0.02, hspace=0.02)
+
+plt.show()   
+
+#%%
+# Regularisation Parameter
+alpha = 5
+
+# Create operators
+#op1 = Gradient(ig)
+op1 = Gradient(ig, correlation='SpaceChannels')
+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 = ZeroFunction()
+    
+# 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, memopt=True)
+pdhg.max_iteration = 2000
+pdhg.update_objective_interval = 200
+pdhg.run(2000)
+
+
+#%%
+fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
+
+plt.subplot(2,3,1)
+plt.imshow(phantom_2Dt[tindex[0],:,:],vmin=0, vmax=1)
+plt.axis('off')
+plt.title('Time {}'.format(tindex[0]))
+
+plt.subplot(2,3,2)
+plt.imshow(phantom_2Dt[tindex[1],:,:],vmin=0, vmax=1)
+plt.axis('off')
+plt.title('Time {}'.format(tindex[1]))
+
+plt.subplot(2,3,3)
+plt.imshow(phantom_2Dt[tindex[2],:,:],vmin=0, vmax=1)
+plt.axis('off')
+plt.title('Time {}'.format(tindex[2]))
+
+
+plt.subplot(2,3,4)
+plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:])
+plt.axis('off')
+plt.subplot(2,3,5)
+plt.imshow(pdhg.get_output().as_array()[tindex[1],:,:])
+plt.axis('off')
+plt.subplot(2,3,6)
+plt.imshow(pdhg.get_output().as_array()[tindex[2],:,:])
+plt.axis('off')
+im = plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:])
+
+
+fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
+                    wspace=0.02, hspace=0.02)
+
+cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
+cbar = fig.colorbar(im, cax=cb_ax)
+
+
+plt.show()
+
diff --git a/Wrappers/Python/wip/demo_box_constraints_FISTA.py b/Wrappers/Python/wip/demo_box_constraints_FISTA.py
index 2f9e0c6..b15dd45 100644
--- a/Wrappers/Python/wip/demo_box_constraints_FISTA.py
+++ b/Wrappers/Python/wip/demo_box_constraints_FISTA.py
@@ -72,7 +72,7 @@ else:
 
 # Set up Operator object combining the ImageGeometry and AcquisitionGeometry
 # wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
-Aop = AstraProjectorSimple(ig, ag, 'gpu')
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
 
 Aop = Identity(ig,ig)