Merge remote-tracking branch 'origin/demos' into update_demo

26 files changed, 1240 insertions, 1257 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
index 647ae98..c8fd0d4 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
@@ -70,7 +70,8 @@ class FISTA(Algorithm):
         
         self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2)))
         
-        self.x.subtract(self.x_old, out=self.y)
+#        self.x.subtract(self.x_old, out=self.y)
+        self.y = self.x - self.x_old
         self.y.__imul__ ((self.t_old-1)/self.t)
         self.y.__iadd__( self.x )
         
diff --git a/Wrappers/Python/data/shapes.png b/Wrappers/Python/data/shapes.png
new file mode 100644
index 0000000..dd4f680
--- /dev/null
+++ b/Wrappers/Python/data/shapes.png
diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
index d1cbe20..653e191 100644
--- a/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
+++ b/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
@@ -82,18 +82,18 @@ plt.colorbar()
 plt.show()
 
 # 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)
+#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)
 
 #%%
 # Show results
diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
index 0875c20..672d4bc 100644
--- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
+++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
@@ -32,7 +32,7 @@ Problem:     min_x || A x - g ||_{2}^{2}
 """
 
 
-from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry
+from ccpi.framework import ImageData, TestData, AcquisitionGeometry
 
 import numpy as np 
 import numpy                          
@@ -42,17 +42,17 @@ from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA
 
 from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition
 from ccpi.astra.ops import AstraProjectorSimple
-import astra          
+import astra   
+import os, sys
 
-# 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)
+# Load Data 
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))  
+                     
+N = 50
+M = 50
+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)
@@ -69,13 +69,14 @@ 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)
+proj_id = astra.create_projector('linear', proj_geom, vol_geom)
 
 # Create a sinogram from a phantom
-sinogram_id, sinogram = astra.create_sino(x, proj_id)
+sinogram_id, sinogram = astra.create_sino(data.as_array(), proj_id)
 
 # Create a data object for the reconstruction
 rec_id = astra.data2d.create('-vol', vol_geom)
@@ -92,6 +93,7 @@ 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)
@@ -99,12 +101,15 @@ 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)
@@ -112,17 +117,17 @@ 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)
+pdhg.run(1000, verbose=True)
 
+###############################################################################
 # 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 = FISTA(x_init=x_init , f=fidelity, g=regularizer)
 fista.max_iteration = 1000
 fista.update_objective_interval = 200
-fista.run(1000, verbose=False)
+fista.run(1000, verbose=True)
 
 #%% Show results
 
@@ -131,18 +136,22 @@ plt.suptitle('Reconstructions ', fontsize=16)
 
 plt.subplot(2,2,1)
 plt.imshow(cgls.get_output().as_array())
+plt.colorbar()
 plt.title('CGLS reconstruction')
 
 plt.subplot(2,2,2)
 plt.imshow(fista.get_output().as_array())
+plt.colorbar()
 plt.title('FISTA reconstruction')
 
 plt.subplot(2,2,3)
 plt.imshow(pdhg.get_output().as_array())
+plt.colorbar()
 plt.title('PDHG reconstruction')
 
 plt.subplot(2,2,4)
 plt.imshow(recon_cgls_astra)
+plt.colorbar()
 plt.title('CGLS astra')
 
 diff1 = pdhg.get_output() - cgls.get_output()
@@ -167,28 +176,14 @@ plt.title('Diff CLGS astra vs CGLS')
 plt.colorbar()
 
 
+#%%
 
+print('Primal Objective (FISTA) {} '.format(fista.objective[-1]))
+print('Primal Objective (CGLS) {} '.format(cgls.objective[-1]))
+print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
 
 
+true_obj = (Aop.direct(cglsd.get_output())-noisy_data).squared_norm()
+print('True objective {}'.format(true_obj))
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-#
-#
-#
-#
-#
-#
-#
-#
diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
index 68fb649..1a96a2d 100644
--- a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
+++ b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
@@ -21,14 +21,15 @@
 
 from ccpi.framework import AcquisitionGeometry
 from ccpi.optimisation.algorithms import FISTA
-from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, L2NormSquared, FunctionOperatorComposition
-from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, \
+                         L2NormSquared, FunctionOperatorComposition
+from ccpi.astra.operators import AstraProjectorSimple
 
 import numpy as np
 import matplotlib.pyplot as plt
 from ccpi.framework import TestData
 import os, sys
-
+from ccpi.optimisation.funcs import Norm2sq
 
 # Load Data 
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
@@ -117,16 +118,15 @@ plt.show()
 # demonstrated in the rest of this file. In general all methods need an initial 
 # guess and some algorithm options to be set:
 
-f = FunctionOperatorComposition(0.5 * L2NormSquared(b=sin), Aop)
+f = FunctionOperatorComposition(L2NormSquared(b=sin), Aop)
+#f = Norm2sq(Aop, sin, c=0.5, memopt=True)
 g = ZeroFunction()
 
 x_init = ig.allocate()
-opt = {'memopt':True}
-
-fista = FISTA(x_init=x_init, f=f, g=g, opt=opt)
-fista.max_iteration = 100
-fista.update_objective_interval = 1
-fista.run(100, verbose=True)
+fista = FISTA(x_init=x_init, f=f, g=g)
+fista.max_iteration = 1000
+fista.update_objective_interval = 100
+fista.run(1000, verbose=True)
 
 plt.figure()
 plt.imshow(fista.get_output().as_array())
@@ -134,18 +134,12 @@ plt.title('FISTA unconstrained')
 plt.colorbar()
 plt.show()
 
-plt.figure()
-plt.semilogy(fista.objective)
-plt.title('FISTA objective')
-plt.show()
-
-#%%
 
 # Run FISTA for least squares with lower/upper bound 
-fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1), opt=opt)
-fista0.max_iteration = 100
-fista0.update_objective_interval = 1
-fista0.run(100, verbose=True)
+fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1))
+fista0.max_iteration = 1000
+fista0.update_objective_interval = 100
+fista0.run(1000, verbose=True)
 
 plt.figure()
 plt.imshow(fista0.get_output().as_array())
@@ -153,10 +147,10 @@ plt.title('FISTA constrained in [0,1]')
 plt.colorbar()
 plt.show()
 
-plt.figure()
-plt.semilogy(fista0.objective)
-plt.title('FISTA constrained in [0,1]')
-plt.show()
+#plt.figure()
+#plt.semilogy(fista0.objective)
+#plt.title('FISTA constrained in [0,1]')
+#plt.show()
 
 #%% Check with CVX solution
 
@@ -178,10 +172,10 @@ if cvx_not_installable:
     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)
+    proj_id = astra.create_projector('linear', proj_geom, vol_geom)
 
     matrix_id = astra.projector.matrix(proj_id)
-
+    
     ProjMat = astra.matrix.get(matrix_id)
     
     tmp = sin.as_array().ravel()
@@ -216,4 +210,6 @@ if cvx_not_installable:
     plt.legend()
     plt.title('Middle Line Profiles')
     plt.show()
-            
\ No newline at end of file
+            
+    print('Primal Objective (CVX) {} '.format(obj.value))
+    print('Primal Objective (FISTA) {} '.format(fista0.loss[1]))    
\ No newline at end of file
diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
index 41219f4..6007990 100644
--- a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
+++ b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
@@ -21,7 +21,7 @@
 
 """ 
 
-Tikhonov for Poisson denoising using FISTA algorithm:
+"Tikhonov regularization" for Poisson denoising using FISTA algorithm:
 
 Problem:     min_x, x>0  \alpha * ||\nabla x||_{2}^{2} + \int x - g * log(x) 
 
@@ -52,14 +52,14 @@ from skimage.util import random_noise
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
 
 # Load Data                      
-N = 150
-M = 150
+N = 100
+M = 100
 data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
 
 ig = data.geometry
 ag = ig
 
-# Create Noisy data. Add Gaussian noise
+# Create Noisy data with Poisson noise
 n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10)
 noisy_data = ImageData(n1)
 
@@ -75,7 +75,6 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
 
-#%%
 # Regularisation Parameter
 alpha = 10
 
@@ -114,8 +113,7 @@ fid.proximal = KL_Prox_PosCone
 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 = FISTA(x_init=x_init , f=reg, g=fid)
 fista.max_iteration = 2000
 fista.update_objective_interval = 500
 fista.run(2000, verbose=True)
@@ -196,7 +194,6 @@ if cvx_not_installable:
     plt.title('Middle Line Profiles')
     plt.show()
             
-    #TODO what is the output of fista.objective, fista.loss
     print('Primal Objective (CVX) {} '.format(obj.value))
     print('Primal Objective (FISTA) {} '.format(fista.loss[1]))
 
diff --git a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py b/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py
index a735323..e69060f 100644
--- a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py
+++ b/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py
@@ -57,9 +57,8 @@ elif phantom == 'powder':
         arrays[k] = numpy.array(v)
     XX = arrays['S']    
     X = numpy.transpose(XX,(0,2,1,3))
-    X = X[0:20]
+    X = X[100:120]
     
-        
     
 #%% Setup Geometry of Colorbay
 
@@ -125,11 +124,16 @@ plt.show()
 
 #%% CGLS
 
+def callback(iteration, objective, x):
+    plt.imshow(x.as_array()[5])
+    plt.colorbar()
+    plt.show()
+
 x_init = ig2d.allocate()      
 cgls1 = CGLS(x_init=x_init, operator=Aall, data=data2d)
 cgls1.max_iteration = 100
 cgls1.update_objective_interval = 1
-cgls1.run(5,verbose=True)
+cgls1.run(5,verbose=True, callback = callback)
 
 plt.imshow(cgls1.get_output().subset(channel=5).array)
 plt.title('CGLS')
@@ -148,7 +152,7 @@ cgls2 = CGLS(x_init=x_init, operator=op_CGLS, data=block_data)
 cgls2.max_iteration = 100
 cgls2.update_objective_interval = 1
 
-cgls2.run(10,verbose=True)
+cgls2.run(10,verbose=True, callback=callback)
 
 plt.imshow(cgls2.get_output().subset(channel=5).array)
 plt.title('Tikhonov')
@@ -174,8 +178,9 @@ f = BlockFunction(f1, f2)
 g = ZeroFunction()
 
 # Compute operator Norm
-normK = 8.70320267279591 # Run one time no need to compute again takes time
-    
+#normK = 8.70320267279591 # For powder Run one time no need to compute again takes time
+normK = 14.60320657253632 # for carbon
+
 # Primal & dual stepsizes
 sigma = 1
 tau = 1/(sigma*normK**2)
@@ -184,11 +189,11 @@ tau = 1/(sigma*normK**2)
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
 pdhg.max_iteration = 2000
 pdhg.update_objective_interval = 100
-pdhg.run(1000, verbose =True)
+pdhg.run(1000, verbose =True, callback=callback)
 
 
 #%% Show sinograms
-channel_ind = [10,15,15]
+channel_ind = [10,15,19]
 
 plt.figure(figsize=(15,15))
 
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
index 8453d20..9dbcf3e 100755
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
@@ -241,7 +241,7 @@ if cvx_not_installable:
         solver = MOSEK
     else:
         solver = SCS      
-    
+        
     # fidelity
     if noise == 's&p':
         fidelity = pnorm( u - noisy_data.as_array(),1)
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
index 74e7901..d848b9f 100755
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
@@ -239,7 +239,7 @@ if cvx_not_installable:
         solver = MOSEK
     else:
         solver = SCS      
-    
+
     # fidelity
     if noise == 's&p':
         fidelity = pnorm( u - noisy_data.as_array(),1)
@@ -248,8 +248,7 @@ if cvx_not_installable:
         solver = SCS
     elif noise == 'gaussian':
         fidelity = 0.5 * sum_squares(noisy_data.as_array() - u)
-    
-            
+                
     obj =  Minimize( regulariser +  fidelity)
     prob = Problem(obj)
     result = prob.solve(verbose = True, solver = solver)
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py
new file mode 100644
index 0000000..febe76d
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py
@@ -0,0 +1,243 @@
+#========================================================================
+# 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.
+#
+#=========================================================================
+""" 
+
+Total Variation (Dynamic) Denoising using PDHG algorithm and Tomophantom:
+
+
+Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
+
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator 
+             
+             g: 2D Dynamic noisy data with Gaussian Noise
+                          
+             K = \nabla  
+                                                                
+"""
+
+from ccpi.framework import ImageData, ImageGeometry
+
+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, L2NormSquared,  \
+                      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 = 128 # set dimension of the phantom
+
+path = os.path.dirname(tomophantom.__file__)
+path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
+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
+
+# Setup geometries    
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0])
+data = ImageData(phantom_2Dt, geometry=ig)
+ag = ig
+                      
+# Create noisy data. Apply Gaussian noise
+np.random.seed(10)
+noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.25, size=ig.shape) )
+
+# time-frames index
+tindex = [8, 16, 24]
+
+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 = 0.3
+
+# Create operators
+op1 = Gradient(ig, correlation='Space')
+op2 = Gradient(ig, correlation='SpaceChannels')
+op3 = Identity(ig, ag)
+
+# Create BlockOperator
+operator1 = BlockOperator(op1, op3, shape=(2,1) ) 
+operator2 = BlockOperator(op2, op3, shape=(2,1) )
+
+# Create functions
+      
+f1 = alpha * MixedL21Norm()
+f2 = 0.5 * L2NormSquared(b = noisy_data)    
+f = BlockFunction(f1, f2)  
+                                                                        
+g = ZeroFunction()
+    
+# Compute operator Norm
+normK1 = operator1.norm()
+normK2 = operator2.norm()
+
+#%%
+# Primal & dual stepsizes
+sigma1 = 1
+tau1 = 1/(sigma1*normK1**2)
+
+sigma2 = 1
+tau2 = 1/(sigma2*normK2**2)
+
+# Setup and run the PDHG algorithm
+pdhg1 = PDHG(f=f,g=g,operator=operator1, tau=tau1, sigma=sigma1)
+pdhg1.max_iteration = 2000
+pdhg1.update_objective_interval = 200
+pdhg1.run(2000)
+
+# Setup and run the PDHG algorithm
+pdhg2 = PDHG(f=f,g=g,operator=operator2, tau=tau2, sigma=sigma2)
+pdhg2.max_iteration = 2000
+pdhg2.update_objective_interval = 200
+pdhg2.run(2000)
+
+
+#%%
+
+tindex = [8, 16, 24]
+fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
+
+plt.subplot(3,3,1)
+plt.imshow(phantom_2Dt[tindex[0],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[0]))
+
+plt.subplot(3,3,2)
+plt.imshow(phantom_2Dt[tindex[1],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[1]))
+
+plt.subplot(3,3,3)
+plt.imshow(phantom_2Dt[tindex[2],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[2]))
+
+plt.subplot(3,3,4)
+plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:])
+plt.axis('off')
+plt.subplot(3,3,5)
+plt.imshow(pdhg1.get_output().as_array()[tindex[1],:,:])
+plt.axis('off')
+plt.subplot(3,3,6)
+plt.imshow(pdhg1.get_output().as_array()[tindex[2],:,:])
+plt.axis('off')
+
+
+plt.subplot(3,3,7)
+plt.imshow(pdhg2.get_output().as_array()[tindex[0],:,:])
+plt.axis('off')
+plt.subplot(3,3,8)
+plt.imshow(pdhg2.get_output().as_array()[tindex[1],:,:])
+plt.axis('off')
+plt.subplot(3,3,9)
+plt.imshow(pdhg2.get_output().as_array()[tindex[2],:,:])
+plt.axis('off')
+
+
+im = plt.imshow(pdhg1.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()
+
+#%%
+import matplotlib.animation as animation
+fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 30))
+ims1 = []
+ims2 = []
+ims3 = []
+for sl in range(0,np.shape(phantom_2Dt)[0]):
+    
+    plt.subplot(1,3,1)
+    im1 = plt.imshow(phantom_2Dt[sl,:,:], animated=True)
+    
+    plt.subplot(1,3,2)
+    im2 = plt.imshow(pdhg1.get_output().as_array()[sl,:,:])
+    
+    plt.subplot(1,3,3)
+    im3 = plt.imshow(pdhg2.get_output().as_array()[sl,:,:])
+    
+    ims1.append([im1])
+    ims2.append([im2])
+    ims3.append([im3])
+    
+    
+ani1 = animation.ArtistAnimation(fig, ims1, interval=500,
+                                repeat_delay=10)
+
+ani2 = animation.ArtistAnimation(fig, ims2, interval=500,
+                                repeat_delay=10)
+
+ani3 = animation.ArtistAnimation(fig, ims3, interval=500,
+                                repeat_delay=10)
+plt.show()        
+#    plt.pause(0.25)
+#    plt.show()
+
+
+
+
+
+
+
+
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py
index c86ddc9..15709cd 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py
@@ -1,42 +1,58 @@
-# -*- 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
-
-import numpy as np 
-import numpy                          
+#========================================================================
+# 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.
+#
+#=========================================================================
+""" 
+
+Total Variation (3D) Denoising using PDHG algorithm and Tomophantom:
+
+
+Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
+
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator 
+             
+             g: 3D Noisy Data with Gaussian Noise
+                          
+             Method = 0 ( PDHG - split ) :  K = [ \nabla,
+                                                 Identity]
+                          
+                                                                    
+             Method = 1 (PDHG - explicit ):  K = \nabla  
+                                                                
+"""
+
+from ccpi.framework import ImageData, ImageGeometry                       
 import matplotlib.pyplot as plt
-
 from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
+from ccpi.optimisation.operators import Gradient
+from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm
 
 from skimage.util import random_noise
 
-# Create phantom for TV Gaussian denoising
 import timeit
 import os
 from tomophantom import TomoP3D
 import tomophantom
 
+# Create a phantom from Tomophantom
 print ("Building 3D phantom using TomoPhantom software")
 tic=timeit.default_timer()
 model = 13 # select a model number from the library
@@ -47,12 +63,13 @@ path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
 #This will generate a N x N x N phantom (3D)
 phantom_tm = TomoP3D.Model(model, N, path_library3D)
 
-# Create noisy data. Add Gaussian noise
+# Create noisy data. Apply Gaussian noise
 ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N)
 ag = ig
 n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10)
 noisy_data = ImageData(n1)
 
+# Show results
 sliceSel = int(0.5*N)
 plt.figure(figsize=(15,15))
 plt.subplot(3,1,1)
@@ -69,52 +86,27 @@ plt.title('Sagittal View')
 plt.colorbar()
 plt.show()   
 
-
 # Regularisation Parameter
 alpha = 0.05
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * MixedL21Norm()
-    f2 = 0.5 * L2NormSquared(b = noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
     
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  0.5 * L2NormSquared(b = noisy_data)
-        
-    
-# Compute operator Norm
+# Setup and run the PDHG algorithm
+operator = Gradient(ig)
+f =  alpha * MixedL21Norm()
+g =  0.5 * L2NormSquared(b = noisy_data)
+            
 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.max_iteration = 1000
 pdhg.update_objective_interval = 200
-pdhg.run(2000)
+pdhg.run(1000, verbose = True)
 
+# Show results
 fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
+fig.suptitle('TV Reconstruction',fontsize=20)
 
 plt.subplot(2,3,1)
 plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1)
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
new file mode 100644
index 0000000..f179e95
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
@@ -0,0 +1,173 @@
+#========================================================================
+# 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.
+#
+#=========================================================================
+
+""" 
+
+Total Variation 2D Tomography Reconstruction using PDHG algorithm:
+
+
+Problem:     min_u  \alpha * ||\nabla u||_{2,1} + \frac{1}{2}||Au - g||^{2}
+             min_u, u>0  \alpha * ||\nabla u||_{2,1} + \int A u  - g log (Au + \eta)
+
+             \nabla: Gradient operator              
+             A: System Matrix
+             g: Noisy sinogram 
+             \eta: Background noise
+             
+             \alpha: Regularization parameter
+ 
+"""
+
+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, L2NormSquared, \
+                      MixedL21Norm, BlockFunction, KullbackLeibler, IndicatorBox
+                      
+from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.framework import TestData
+from PIL import Image
+import os, sys
+if int(numpy.version.version.split('.')[1]) > 12:
+    from skimage.util import random_noise
+else:
+    from demoutil import random_noise
+
+import scipy.io
+
+# user supplied input
+if len(sys.argv) > 1:
+    which_noise = int(sys.argv[1])
+else:
+    which_noise = 1
+
+# Load 256 shepp-logan
+data256 = scipy.io.loadmat('phantom.mat')['phantom256']
+data = ImageData(numpy.array(Image.fromarray(data256).resize((256,256))))
+N, M = data.shape
+ig = ImageGeometry(voxel_num_x=N, voxel_num_y=M)
+
+# Add it to testdata or use tomophantom
+#loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+#data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(50, 50))
+#ig = data.geometry
+
+# Create acquisition data and geometry
+detectors = N
+angles = np.linspace(0, np.pi, 180)
+ag = AcquisitionGeometry('parallel','2D',angles, detectors)
+
+# Select device
+device = '0'
+#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)
+
+# Create noisy data. Apply Gaussian noise
+noises = ['gaussian', 'poisson']
+noise = noises[which_noise]
+
+if noise == 'poisson':
+    scale = 5
+    eta = 0
+    noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag)
+elif noise == 'gaussian':
+    n1 = np.random.normal(0, 1, size = ag.shape)
+    noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
+else:
+    raise ValueError('Unsupported Noise ', noise)
+
+# Show Ground Truth and Noisy Data
+plt.figure(figsize=(10,10))
+plt.subplot(1,2,2)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(1,2,1)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.show()
+
+# Create operators
+op1 = Gradient(ig)
+op2 = Aop
+
+# Create BlockOperator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+# Compute operator Norm
+normK = operator.norm()
+
+# Create functions
+if noise == 'poisson':
+    alpha = 3
+    f2 = KullbackLeibler(noisy_data)  
+    g =  IndicatorBox(lower=0) 
+    sigma = 1
+    tau = 1/(sigma*normK**2)    
+    
+elif noise == 'gaussian':   
+    alpha = 20
+    f2 = 0.5 * L2NormSquared(b=noisy_data)                                         
+    g = ZeroFunction()
+    sigma = 10
+    tau = 1/(sigma*normK**2)     
+    
+f1 = alpha * MixedL21Norm() 
+f = BlockFunction(f1, f2)    
+
+# 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)
+
+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,ig.shape[1],ig.shape[1]), data.as_array()[int(N/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
index e16c5a6..78d4980 100644
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
@@ -195,7 +195,6 @@ try:
 except ImportError:
     cvx_not_installable = False
 
-
 if cvx_not_installable:
 
     ##Construct problem    
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat b/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat
new file mode 100755
index 0000000..c465bbe
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat
diff --git a/Wrappers/Python/demos/PDHG_examples/IMATDemo.py b/Wrappers/Python/demos/PDHG_examples/IMATDemo.py
new file mode 100644
index 0000000..2051860
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/IMATDemo.py
@@ -0,0 +1,339 @@
+
+
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 25 12:50:27 2019
+
+@author: vaggelis
+"""
+
+from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData
+from astropy.io import fits
+import numpy
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import CGLS, PDHG
+from ccpi.optimisation.functions import MixedL21Norm, L2NormSquared, BlockFunction, ZeroFunction, KullbackLeibler, IndicatorBox
+from ccpi.optimisation.operators import Gradient, BlockOperator
+
+from ccpi.astra.operators import  AstraProjectorMC, AstraProjectorSimple
+
+import pickle
+
+
+# load file
+
+#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_282.fits' 
+#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_564.fits'
+#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_141.fits'
+filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_80_channels.fits'
+
+sino_handler = fits.open(filename_sino)
+sino = numpy.array(sino_handler[0].data, dtype=float)
+
+# change axis order: channels, angles, detectors
+sino_new = numpy.rollaxis(sino, 2)
+sino_handler.close()
+
+
+sino_shape = sino_new.shape
+
+num_channels = sino_shape[0] # channelss
+num_pixels_h = sino_shape[2] # detectors
+num_pixels_v = sino_shape[2] # detectors
+num_angles = sino_shape[1] # angles
+
+
+ig = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v, channels = num_channels)
+
+with open("/media/newhd/vaggelis/CCPi/IMAT_reconstruction/CCPi-Framework/Wrappers/Python/ccpi/optimisation/IMAT_data/golden_angles_new.txt") as f:
+    angles_string = [line.rstrip() for line in f]
+    angles = numpy.array(angles_string).astype(float)
+                                                                                 
+
+ag = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h, channels = num_channels)
+op_MC = AstraProjectorMC(ig, ag, 'gpu')
+
+sino_aqdata = AcquisitionData(sino_new, ag)
+result_bp = op_MC.adjoint(sino_aqdata)
+
+#%%
+
+channel = [40, 60]
+for j in range(2):        
+    z4 = sino_aqdata.as_array()[channel[j]] 
+    plt.figure(figsize=(10,6))
+    plt.imshow(z4, cmap='viridis')
+    plt.axis('off')
+    plt.savefig('Sino_141/Sinogram_ch_{}_.png'.format(channel[j]), bbox_inches='tight', transparent=True)
+    plt.show() 
+
+#%%
+
+def callback(iteration, objective, x):
+    plt.imshow(x.as_array()[40])
+    plt.colorbar()
+    plt.show()
+
+#%%
+# CGLS 
+
+x_init = ig.allocate()      
+cgls1 = CGLS(x_init=x_init, operator=op_MC, data=sino_aqdata)
+cgls1.max_iteration = 100
+cgls1.update_objective_interval = 2
+cgls1.run(20,verbose=True, callback=callback)
+
+plt.imshow(cgls1.get_output().subset(channel=20).array)
+plt.title('CGLS')
+plt.colorbar()
+plt.show()
+
+#%%
+with open('Sino_141/CGLS/CGLS_{}_iter.pkl'.format(20), 'wb') as f:
+    z = cgls1.get_output()
+    pickle.dump(z, f)  
+    
+#%% 
+#% Tikhonov Space
+    
+x_init = ig.allocate() 
+alpha = [1,3,5,10,20,50]
+
+for a in alpha:
+    
+    Grad = Gradient(ig, correlation = Gradient.CORRELATION_SPACE)
+    operator = BlockOperator(op_MC, a * Grad, shape=(2,1))
+    blockData = BlockDataContainer(sino_aqdata, \
+                                   Grad.range_geometry().allocate())
+    cgls2 = CGLS()
+    cgls2.max_iteration = 500
+    cgls2.set_up(x_init, operator, blockData)
+    cgls2.update_objective_interval = 50
+    cgls2.run(100,verbose=True)
+        
+    with open('Sino_141/CGLS_Space/CGLS_a_{}.pkl'.format(a), 'wb') as f:
+        z = cgls2.get_output()
+        pickle.dump(z, f)
+   
+#% Tikhonov SpaceChannels
+
+for a1 in alpha:
+    
+    Grad1 = Gradient(ig, correlation = Gradient.CORRELATION_SPACECHANNEL)
+    operator1 = BlockOperator(op_MC, a1 * Grad1, shape=(2,1))
+    blockData1 = BlockDataContainer(sino_aqdata, \
+                                   Grad1.range_geometry().allocate())
+    cgls3 = CGLS()
+    cgls3.max_iteration = 500
+    cgls3.set_up(x_init, operator1, blockData1)
+    cgls3.update_objective_interval = 10
+    cgls3.run(100, verbose=True)
+    
+    with open('Sino_141/CGLS_SpaceChannels/CGLS_a_{}.pkl'.format(a1), 'wb') as f1:
+        z1 = cgls3.get_output()
+        pickle.dump(z1, f1)
+    
+  
+    
+#%%
+#
+ig_tmp = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v)
+ag_tmp = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h)
+op_tmp = AstraProjectorSimple(ig_tmp, ag_tmp, 'gpu')
+normK1 = op_tmp.norm()
+
+alpha_TV = [2, 5, 10] # for powder
+
+# Create operators
+op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL)
+op2 = op_MC
+
+# Create BlockOperator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+
+for alpha in alpha_TV:
+# Create functions      
+    f1 = alpha * MixedL21Norm()
+    
+    f2 = KullbackLeibler(sino_aqdata)    
+    f = BlockFunction(f1, f2)  
+    g = IndicatorBox(lower=0)
+    
+    # Compute operator Norm
+    normK = numpy.sqrt(8 + normK1**2) 
+    
+    # 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 = 5000
+    pdhg.update_objective_interval = 500
+#    pdhg.run(2000, verbose=True, callback=callback)   
+    pdhg.run(5000, verbose=True, callback=callback)  
+#    
+    with open('Sino_141/TV_SpaceChannels/TV_a = {}.pkl'.format(alpha), 'wb') as f3:
+        z3 = pdhg.get_output()
+        pickle.dump(z3, f3)    
+#    
+#
+#
+#
+##%%
+#        
+#ig_tmp = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v)
+#ag_tmp = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h)
+#op_tmp = AstraProjectorSimple(ig_tmp, ag_tmp, 'gpu')
+#normK1 = op_tmp.norm()
+#
+#alpha_TV = 10 # for powder
+#
+## Create operators
+#op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL)
+#op2 = op_MC
+#
+## Create BlockOperator
+#operator = BlockOperator(op1, op2, shape=(2,1) ) 
+#
+#
+## Create functions      
+#f1 = alpha_TV * MixedL21Norm()
+#f2 = 0.5 * L2NormSquared(b=sino_aqdata)    
+#f = BlockFunction(f1, f2)  
+#g = ZeroFunction()
+#
+## Compute operator Norm
+##normK = 8.70320267279591 # For powder Run one time no need to compute again takes time
+#normK = numpy.sqrt(8 + normK1**2) # for carbon
+#
+## Primal & dual stepsizes
+#sigma = 0.1
+#tau = 1/(sigma*normK**2)
+#
+#def callback(iteration, objective, x):
+#    plt.imshow(x.as_array()[100])
+#    plt.colorbar()
+#    plt.show()
+#
+## 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 = 100
+#pdhg.run(2000, verbose=True)       
+#        
+#        
+#        
+#        
+#        
+#        
+        
+        
+        
+        
+        
+        
+        
+        
+        
+        
+#%%
+    
+#with open('/media/newhd/vaggelis/CCPi/IMAT_reconstruction/CCPi-Framework/Wrappers/Python/ccpi/optimisation/CGLS_Tikhonov/CGLS_Space/CGLS_Space_a = 50.pkl', 'wb') as f:
+#    z = cgls2.get_output()
+#    pickle.dump(z, f)
+#      
+        
+        #%%
+with open('Sino_141/CGLS_Space/CGLS_Space_a_20.pkl', 'rb') as f1:
+    x = pickle.load(f1)
+    
+with open('Sino_141/CGLS_SpaceChannels/CGLS_SpaceChannels_a_20.pkl', 'rb') as f1:
+    x1 = pickle.load(f1)    
+    
+   
+    
+#    
+plt.imshow(x.as_array()[40]*mask)
+plt.colorbar()
+plt.show()
+
+plt.imshow(x1.as_array()[40]*mask)
+plt.colorbar()
+plt.show()
+
+plt.plot(x.as_array()[40,100,:])
+plt.plot(x1.as_array()[40,100,:])
+plt.show()
+
+#%%
+
+# Show results
+
+def circ_mask(h, w, center=None, radius=None):
+
+    if center is None: # use the middle of the image
+        center = [int(w/2), int(h/2)]
+    if radius is None: # use the smallest distance between the center and image walls
+        radius = min(center[0], center[1], w-center[0], h-center[1])
+
+    Y, X = numpy.ogrid[:h, :w]
+    dist_from_center = numpy.sqrt((X - center[0])**2 + (Y-center[1])**2)
+
+    mask = dist_from_center <= radius
+    return mask
+
+mask = circ_mask(141, 141, center=None, radius = 55)
+plt.imshow(numpy.multiply(x.as_array()[40],mask))
+plt.show()
+#%%
+#channel = [100, 200, 300]
+#
+#for i in range(3):     
+#    tmp = cgls1.get_output().as_array()[channel[i]] 
+#    
+#    z =  tmp * mask
+#    plt.figure(figsize=(10,6))
+#    plt.imshow(z, vmin=0, cmap='viridis')
+#    plt.axis('off')
+##    plt.clim(0, 0.02)
+##    plt.colorbar()
+##    del z
+#    plt.savefig('CGLS_282/CGLS_Chan_{}.png'.format(channel[i]), bbox_inches='tight', transparent=True)
+#    plt.show()
+#    
+#    
+##%% Line Profiles
+#    
+#n1, n2, n3 =  cgs.get_output().as_array().shape
+#mask = circ_mask(564, 564, center=None, radius = 220)
+#material = ['Cu', 'Fe', 'Ni']
+#ycoords = [200, 300, 380]
+#    
+#for i in range(3):
+#    z = cgs.get_output().as_array()[channel[i]]  * mask
+#    
+#    for k1 in range(len(ycoords)):
+#        plt.plot(numpy.arange(0,n2), z[ycoords[k1],:])
+#        plt.title('Channel {}: {}'.format(channel[i], material[k1]))
+#        plt.savefig('CGLS/line_profile_chan_{}_material_{}.png'.\
+#                        format(channel[i], material[k1]), bbox_inches='tight')
+#        plt.show()
+#    
+#    
+#
+#
+#
+##%%
+#
+#%%
+        
+
+
+#%%
+
+#plt.imshow(pdhg.get_output().subset(channel=100).as_array())
+#plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py
new file mode 100644
index 0000000..ddf5ace
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py
@@ -0,0 +1,115 @@
+#========================================================================
+# 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.
+#
+#=========================================================================
+
+""" 
+Total Variation Denoising using PDHG algorithm:
+Problem:     min_x, x>0  \alpha * ||\nabla x||_{2,1} + ||x-g||_{1}
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator 
+             
+             g: Noisy Data with Salt & Pepper Noise
+             
+             
+             Method = 0 ( PDHG - split ) :  K = [ \nabla,
+                                                 Identity]
+                          
+                                                                    
+             Method = 1 (PDHG - explicit ):  K = \nabla    
+             
+             
+"""
+
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff
+from ccpi.optimisation.functions import MixedL21Norm, MixedL11Norm, L2NormSquared, BlockFunction, L1Norm                      
+from ccpi.framework import TestData, ImageGeometry
+import os, sys
+if int(numpy.version.version.split('.')[1]) > 12:
+    from skimage.util import random_noise
+else:
+    from demoutil import random_noise
+
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+data = loader.load(TestData.PEPPERS, size=(256,256))
+ig = data.geometry
+ag = ig
+
+# Create noisy data. 
+n1 = random_noise(data.as_array(), mode = 'gaussian', var = 0.15, seed = 50)
+noisy_data = ig.allocate()
+noisy_data.fill(n1)
+
+# Show Ground Truth and Noisy Data
+plt.figure(figsize=(10,5))
+plt.subplot(1,2,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(1,2,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.show()
+
+# Regularisation Parameter
+operator = Gradient(ig, correlation=Gradient.CORRELATION_SPACE)
+f1 =  5 * MixedL21Norm()
+g = 0.5 * L2NormSquared(b=noisy_data)
+            
+# Compute operator Norm
+normK = operator.norm()
+
+# Primal & dual stepsizes
+sigma = 1
+tau = 1/(sigma*normK**2)
+
+# Setup and run the PDHG algorithm
+pdhg1 = PDHG(f=f1,g=g,operator=operator, tau=tau, sigma=sigma)
+pdhg1.max_iteration = 2000
+pdhg1.update_objective_interval = 200
+pdhg1.run(1000)
+
+
+# Show results
+plt.figure(figsize=(10,10))
+plt.subplot(2,2,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(2,2,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(2,2,3)
+plt.imshow(pdhg1.get_output().as_array())
+plt.title('TV Reconstruction')
+plt.colorbar()
+plt.subplot(2,2,4)
+plt.imshow(pdhg2.get_output().as_array())
+plt.title('TV Reconstruction')
+plt.colorbar()
+plt.show()
+
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
index bc0081f..b2af753 100644
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
+++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
@@ -23,18 +23,20 @@
 Total Generalised Variation (TGV) Tomography 2D using PDHG algorithm:
 
 
-Problem:     min_{x>0} \alpha * ||\nabla x - w||_{2,1} +
-                     \beta *  || E w ||_{2,1} +
-                     int A x - g log(Ax + \eta)
+Problem:     min_{x>0} \alpha * ||\nabla x - w||_{2,1} + \beta *  || E w ||_{2,1} +
+                       \frac{1}{2}||Au - g||^{2}
+                     
+            min_{u>0} \alpha * ||\nabla u - w||_{2,1} + \beta *  || E w ||_{2,1} +
+                      int A u - g log(Au + \eta)                     
 
              \alpha: Regularization parameter
              \beta: Regularization parameter
              
              \nabla: Gradient operator 
               E: Symmetrized Gradient operator
-              A: Projection Matrix
+              A: System Matrix
              
-             g: Noisy Data with Poisson Noise
+             g: Noisy Sinogram 
                           
               K = [ \nabla, - Identity
                    ZeroOperator, E 
@@ -58,9 +60,12 @@ from ccpi.optimisation.functions import IndicatorBox, KullbackLeibler, ZeroFunct
 from ccpi.astra.ops import AstraProjectorSimple
 from ccpi.framework import TestData
 import os, sys
-from skimage.util import random_noise
+if int(numpy.version.version.split('.')[1]) > 12:
+    from skimage.util import random_noise
+else:
+    from demoutil import random_noise
 
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+#loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
 
 # Load Data                      
 #N = 50
@@ -85,15 +90,8 @@ data = ImageData(data/data.max())
 ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
 ag = ig
 
-# Create noisy data. Add Gaussian noise
-scale = 0.5
-eta = 0 
-n1 = scale * np.random.poisson(eta + sin.as_array()/scale)
-
-noisy_data = AcquisitionData(n1, ag)
 
 #Create Acquisition Data and apply poisson noise
-
 detectors = N
 angles = np.linspace(0, np.pi, N)
 
@@ -109,28 +107,31 @@ else:
 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)
+# Create noisy data.
+noises = ['gaussian', 'poisson']
+noise = noises[which_noise]
+
+if noise == 'poisson':
+    scale = 5
+    eta = 0
+    noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag)
+elif noise == 'gaussian':
+    n1 = np.random.normal(0, 1, size = ag.shape)
+    noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
+else:
+    raise ValueError('Unsupported Noise ', noise)
 
 # Show Ground Truth and Noisy Data
 plt.figure(figsize=(10,10))
-plt.subplot(2,1,1)
+plt.subplot(1,2,2)
 plt.imshow(data.as_array())
 plt.title('Ground Truth')
 plt.colorbar()
-plt.subplot(2,1,2)
+plt.subplot(1,2,1)
 plt.imshow(noisy_data.as_array())
 plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
-#%%
-# Regularisation Parameters
-alpha = 1
-beta =  5
 
 # Create Operators
 op11 = Gradient(ig)
@@ -143,28 +144,41 @@ op31 = Aop
 op32 = ZeroOperator(op22.domain_geometry(), ag)
 
 operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) 
+
+# Create functions
+if noise == 'poisson':
+    alpha = 1
+    beta = 5
+    f3 = KullbackLeibler(noisy_data)    
+    g =  BlockFunction(IndicatorBox(lower=0), ZeroFunction()) 
+    
+    # Primal & dual stepsizes
+    sigma = 1
+    tau = 1/(sigma*normK**2)    
+    
+elif noise == 'gaussian':   
+    alpha = 20
+    f3 = 0.5 * L2NormSquared(b=noisy_data)                                         
+    g = BlockFunction(ZeroFunction(), ZeroFunction())
+    
+    # Primal & dual stepsizes
+    sigma = 10
+    tau = 1/(sigma*normK**2)     
     
 f1 = alpha * MixedL21Norm()
-f2 = beta * MixedL21Norm() 
-f3 = KullbackLeibler(noisy_data)    
+f2 = beta * MixedL21Norm()  
 f = BlockFunction(f1, f2, f3)         
-
-g =  BlockFunction(IndicatorBox(lower=0), 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)
 pdhg.max_iteration = 3000
 pdhg.update_objective_interval = 500
 pdhg.run(3000)
 
+#%%
 plt.figure(figsize=(15,15))
 plt.subplot(3,1,1)
 plt.imshow(data.as_array())
@@ -179,9 +193,8 @@ plt.imshow(pdhg.get_output()[0].as_array())
 plt.title('TGV 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()[0].as_array()[int(N/2),:], label = 'TGV reconstruction')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(N/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TGV reconstruction')
 plt.legend()
 plt.title('Middle Line Profiles')
 plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py
new file mode 100644
index 0000000..f179e95
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py
@@ -0,0 +1,173 @@
+#========================================================================
+# 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.
+#
+#=========================================================================
+
+""" 
+
+Total Variation 2D Tomography Reconstruction using PDHG algorithm:
+
+
+Problem:     min_u  \alpha * ||\nabla u||_{2,1} + \frac{1}{2}||Au - g||^{2}
+             min_u, u>0  \alpha * ||\nabla u||_{2,1} + \int A u  - g log (Au + \eta)
+
+             \nabla: Gradient operator              
+             A: System Matrix
+             g: Noisy sinogram 
+             \eta: Background noise
+             
+             \alpha: Regularization parameter
+ 
+"""
+
+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, L2NormSquared, \
+                      MixedL21Norm, BlockFunction, KullbackLeibler, IndicatorBox
+                      
+from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.framework import TestData
+from PIL import Image
+import os, sys
+if int(numpy.version.version.split('.')[1]) > 12:
+    from skimage.util import random_noise
+else:
+    from demoutil import random_noise
+
+import scipy.io
+
+# user supplied input
+if len(sys.argv) > 1:
+    which_noise = int(sys.argv[1])
+else:
+    which_noise = 1
+
+# Load 256 shepp-logan
+data256 = scipy.io.loadmat('phantom.mat')['phantom256']
+data = ImageData(numpy.array(Image.fromarray(data256).resize((256,256))))
+N, M = data.shape
+ig = ImageGeometry(voxel_num_x=N, voxel_num_y=M)
+
+# Add it to testdata or use tomophantom
+#loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+#data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(50, 50))
+#ig = data.geometry
+
+# Create acquisition data and geometry
+detectors = N
+angles = np.linspace(0, np.pi, 180)
+ag = AcquisitionGeometry('parallel','2D',angles, detectors)
+
+# Select device
+device = '0'
+#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)
+
+# Create noisy data. Apply Gaussian noise
+noises = ['gaussian', 'poisson']
+noise = noises[which_noise]
+
+if noise == 'poisson':
+    scale = 5
+    eta = 0
+    noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag)
+elif noise == 'gaussian':
+    n1 = np.random.normal(0, 1, size = ag.shape)
+    noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
+else:
+    raise ValueError('Unsupported Noise ', noise)
+
+# Show Ground Truth and Noisy Data
+plt.figure(figsize=(10,10))
+plt.subplot(1,2,2)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(1,2,1)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.show()
+
+# Create operators
+op1 = Gradient(ig)
+op2 = Aop
+
+# Create BlockOperator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+# Compute operator Norm
+normK = operator.norm()
+
+# Create functions
+if noise == 'poisson':
+    alpha = 3
+    f2 = KullbackLeibler(noisy_data)  
+    g =  IndicatorBox(lower=0) 
+    sigma = 1
+    tau = 1/(sigma*normK**2)    
+    
+elif noise == 'gaussian':   
+    alpha = 20
+    f2 = 0.5 * L2NormSquared(b=noisy_data)                                         
+    g = ZeroFunction()
+    sigma = 10
+    tau = 1/(sigma*normK**2)     
+    
+f1 = alpha * MixedL21Norm() 
+f = BlockFunction(f1, f2)    
+
+# 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)
+
+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,ig.shape[1],ig.shape[1]), data.as_array()[int(N/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat b/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat
new file mode 100755
index 0000000..c465bbe
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat
diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py
deleted file mode 100644
index 49d4db6..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py
+++ /dev/null
@@ -1,124 +0,0 @@
-# -*- 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, Identity, \
-                                    SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.ops import AstraProjectorSimple
-
-# Create phantom for TV 2D tomography 
-N = 75
-
-data = np.zeros((N,N))
-
-x1 = np.linspace(0, int(N/2), N)
-x2 = np.linspace(int(N/2), 0., N)
-xv, yv = np.meshgrid(x1, x2)
-
-xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T
-data = xv
-data = ImageData(data/data.max())
-
-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, 'gpu')
-sin = Aop.direct(data)
-
-# Create noisy data. Apply Poisson noise
-scale = 0.1
-np.random.seed(5)
-n1 = scale * np.random.poisson(sin.as_array()/scale)
-noisy_data = AcquisitionData(n1, ag)
-
-
-plt.imshow(noisy_data.as_array())
-plt.show()
-#%%
-# Regularisation Parameters
-alpha = 0.7
-beta =  2
-
-# Create Operators
-op11 = Gradient(ig)
-op12 = Identity(op11.range_geometry())
-
-op22 = SymmetrizedGradient(op11.domain_geometry())    
-op21 = ZeroOperator(ig, op22.range_geometry())
-    
-op31 = Aop
-op32 = ZeroOperator(op22.domain_geometry(), ag)
-
-operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) 
-    
-f1 = alpha * MixedL21Norm()
-f2 = beta * MixedL21Norm() 
-f3 = KullbackLeibler(noisy_data)    
-f = BlockFunction(f1, f2, f3)         
-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 = 50
-pdhg.run(2000)
-
-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()[0].as_array())
-plt.title('TGV 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()[0].as_array()[int(N/2),:], label = 'TGV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
deleted file mode 100644
index 5df02b1..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
+++ /dev/null
@@ -1,225 +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.
-#
-#=========================================================================
-
-
-""" 
-
-Total Variation Denoising using PDHG algorithm:
-
-             min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) 
-
-
-Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator 
-             
-             g: Noisy Data with Gaussian Noise
-                          
-             Method = 0:  K = [ \nabla,
-                                 Identity]
-                                                                    
-             Method = 1:  K = \nabla    
-             
-             
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-                      
-  
-
-
-from Data import *
-
-#%%
-
-
-data = ImageData(plt.imread('camera.png'))
-
-#
-## Create phantom for TV Gaussian denoising
-#N = 200
-#
-#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
-#
-#
-#
-## Replace with http://sipi.usc.edu/database/database.php?volume=misc&image=36#top
-
-
-
-# Create noisy data. Add Gaussian noise
-np.random.seed(10)
-noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.05, size=ig.shape) )
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(15,15))
-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
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions      
-    f1 = alpha * MixedL21Norm()
-    f2 = 0.5 * L2NormSquared(b = noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator =  Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  0.5 * L2NormSquared(b = noisy_data)
-        
-# 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 = 200
-pdhg.run(3000, verbose=False)
-
-# 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(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
-
-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 = 0.5 * 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)
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-        
-    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(u.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,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.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'Truth')
-    
-    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]))
-
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
deleted file mode 100644
index 70f6b9b..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
+++ /dev/null
@@ -1,207 +0,0 @@
-# -*- 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 STFC,  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
-
-#   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.
-
-""" 
-
-Total Variation Denoising using PDHG algorithm:
-
-             min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) 
-
-
-Problem:     min_x, x>0  \alpha * ||\nabla x||_{1} + \int x - g * log(x)
-
-             \nabla: Gradient operator              
-             g: Noisy Data with Poisson Noise
-             \alpha: Regularization parameter
-             
-             Method = 0:  K = [ \nabla,
-                                 Identity]
-                                                                    
-             Method = 1:  K = \nabla    
-             
-             
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Poisson denoising
-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
-
-# Create noisy data. Apply Poisson noise
-n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 2
-
-method = '1'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * MixedL21Norm()
-    f2 = KullbackLeibler(noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  KullbackLeibler(noisy_data)   
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-opt = {'niter':2000, 'memopt': True}
-
-# 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
-
-def pdgap_objectives(niter, objective, solution):
-    
-
-     print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\
-                      format(niter, pdhg.max_iteration,'', \
-                             objective[0],'',\
-                             objective[1],'',\
-                             objective[2]))
-
-pdhg.run(2000, callback = pdgap_objectives)
-
-
-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
-
-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')
-    
-    # Define Total Variation as a regulariser
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0))
-    
-    fidelity = sum( u1 - multiply(noisy_data.as_array(), log(u1)) )
-    constraints = [q>= fidelity, u1>=0]
-        
-    solver = ECOS
-    obj =  Minimize( regulariser +  q)
-    prob = Problem(obj, constraints)
-    result = prob.solve(verbose = True, solver = solver)
-
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u1.value )
-        
-    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(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,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), 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 (PDHG) {} '.format(pdhg.objective[-1][0]))
-
-
-
-
-
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
deleted file mode 100644
index f5d4ce4..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
+++ /dev/null
@@ -1,198 +0,0 @@
-# -*- 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.
-
-""" 
-
-Total Variation Denoising using PDHG algorithm:
-
-             min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) 
-
-
-Problem:     min_x, x>0  \alpha * ||\nabla x||_{1} + ||x-g||_{1}
-
-             \nabla: Gradient operator 
-             g: Noisy Data with Salt & Pepper Noise
-             \alpha: Regularization parameter
-             
-             Method = 0:  K = [ \nabla,
-                                 Identity]
-                                                                    
-             Method = 1:  K = \nabla    
-             
-             
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Salt & Pepper denoising
-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
-
-# Create noisy data. Apply Salt & Pepper noise
-n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 2
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * MixedL21Norm()
-    f2 = L1Norm(b = noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  L1Norm(b = noisy_data)
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-opt = {'niter':2000, 'memopt': True}
-
-# 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)
-
-
-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
-
-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)
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-        
-    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(u.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,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]))
-
-
-
-
-
diff --git a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py
deleted file mode 100644
index 041d4ee..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py
+++ /dev/null
@@ -1,176 +0,0 @@
-# -*- 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
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared,  BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Salt & Pepper denoising
-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
-
-# Create noisy data. Apply Salt & Pepper noise
-n1 = random_noise(data.as_array(), mode = 'gaussian', mean=0, var = 0.05, seed=10)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 4
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * L2NormSquared()
-    f2 = 0.5 * L2NormSquared(b = noisy_data)    
-    f = BlockFunction(f1, f2)                                       
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * L2NormSquared()
-    g =  0.5 * L2NormSquared(b = noisy_data)
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-opt = {'niter':2000, 'memopt': True}
-
-# 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)
-
-
-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('Tikhonov 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 = 'Tikhonov 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    
-    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_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    fidelity = 0.5 * 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 = solver)
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-        
-    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(u.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,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]))
-
-
-
-
-
diff --git a/Wrappers/Python/wip/Demos/fista_test.py b/Wrappers/Python/wip/Demos/fista_test.py
deleted file mode 100644
index dd1f6fa..0000000
--- a/Wrappers/Python/wip/Demos/fista_test.py
+++ /dev/null
@@ -1,127 +0,0 @@
-# -*- 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
-
-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, \
-                                        MixedL21Norm, FunctionOperatorComposition, BlockFunction, ZeroFunction
-                                                
-from skimage.util import random_noise
-
-# Create phantom for TV Gaussian denoising
-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
-
-# Create noisy data. Add Gaussian noise
-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
-
-operator = Gradient(ig)
-
-#fidelity = L1Norm(b=noisy_data)
-#regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator)
-
-fidelity = FunctionOperatorComposition(alpha * MixedL21Norm(), operator)
-regulariser = 0.5 * L2NormSquared(b = noisy_data)
-
-x_init = ig.allocate()
-
-## Setup and run the PDHG algorithm
-opt = {'tol': 1e-4, '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=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(fista.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.show()
-
-# Compare with PDHG
-method = '0'
-#
-if method == '0':
-#
-#    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-#
-#    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) )   
-    f = BlockFunction(alpha * L2NormSquared(), fidelity)                                        
-    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 = 50
-pdhg.run(2000)
-#
-#%%
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(fista.get_output().as_array())
-plt.title('FISTA')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('PDHG')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(np.abs(pdhg.get_output().as_array()-fista.get_output().as_array()))
-plt.title('Diff FISTA-PDHG')
-plt.colorbar()
-plt.show()
-
-
diff --git a/Wrappers/Python/wip/demo_colourbay.py b/Wrappers/Python/wip/demo_colourbay.py
index 5dbf2e1..0536b07 100644
--- a/Wrappers/Python/wip/demo_colourbay.py
+++ b/Wrappers/Python/wip/demo_colourbay.py
@@ -18,7 +18,7 @@ from ccpi.optimisation.funcs import Norm2sq, Norm1
 # Permute (numpy.transpose) puts into our default ordering which is 
 # (channel, angle, vertical, horizontal).
 
-pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/ColourBay/spectral_data_sets/CarbonPd/'
+pathname = '/media/newhd/shared/Data/ColourBay/spectral_data_sets/CarbonPd/'
 filename = 'carbonPd_full_sinogram_stripes_removed.mat'
 
 X = loadmat(pathname + filename)