add complete Demos PDHG

6 files changed, 814 insertions, 0 deletions

diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py
new file mode 100644
index 0000000..dffb6d8
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py
@@ -0,0 +1,196 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Feb 22 14:53:03 2019
+
+@author: evangelos
+"""
+
+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, SymmetrizedGradient, ZeroOperator
+from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
+                      MixedL21Norm, BlockFunction
+
+from skimage.util import random_noise
+
+# Create phantom for TGV SaltPepper denoising
+
+N = 300
+
+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)
+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 Parameters
+alpha = 0.8
+beta = numpy.sqrt(2)* alpha
+
+method = '0'
+
+if method == '0':
+
+    # Create operators
+    op11 = Gradient(ig)
+    op12 = Identity(op11.range_geometry())
+    
+    op22 = SymmetrizedGradient(op11.domain_geometry())    
+    op21 = ZeroOperator(ig, op22.range_geometry())
+        
+    op31 = Identity(ig, ag)
+    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 = L1Norm(b=noisy_data)    
+    f = BlockFunction(f1, f2, f3)         
+    g = ZeroFunction()
+        
+else:
+    
+    # Create operators
+    op11 = Gradient(ig)
+    op12 = Identity(op11.range_geometry())
+    op22 = SymmetrizedGradient(op11.domain_geometry())    
+    op21 = ZeroOperator(ig, op22.range_geometry())    
+    
+    operator = BlockOperator(op11, -1*op12, op21, op22, shape=(2,2) )      
+    
+    f1 = alpha * MixedL21Norm()
+    f2 = beta * MixedL21Norm()     
+    
+    f = BlockFunction(f1, f2)         
+    g = BlockFunction(0.5 * L1Norm(b=noisy_data), 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 = '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:    
+#    
+#    u = Variable(ig.shape)
+#    w1 = Variable((N, N))
+#    w2 = Variable((N, N))
+#    
+#    # create TGV regulariser
+#    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+#    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+#    
+#    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \
+#                                           DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \
+#                  beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \
+#                                      0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \
+#                                      0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0  ) )  
+#    
+#    constraints = []
+#    fidelity = 0.5 * sum_squares(u - noisy_data.as_array())    
+#    solver = MOSEK
+#
+#    obj =  Minimize( regulariser +  fidelity)
+#    prob = Problem(obj)
+#    result = prob.solve(verbose = True, solver = solver)
+#    
+#    diff_cvx = numpy.abs( res[0].as_array() - u.value )   
+#    
+#    # Show result
+#    plt.figure(figsize=(15,15))
+#    plt.subplot(3,1,1)
+#    plt.imshow(res[0].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), res[0].as_array()[int(N/2),:], label = 'PDHG')
+#    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
+#    plt.legend()   
+#    
+#    print('Primal Objective (CVX) {} '.format(obj.value))
+#    print('Primal Objective (PDHG) {} '.format(primal[-1])) 
+#    print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max()))
+#    
+    
+   
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
new file mode 100644
index 0000000..10e5621
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
@@ -0,0 +1,177 @@
+# -*- 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, \
+                      MixedL21Norm, BlockFunction
+
+from skimage.util import random_noise
+
+# Create phantom for TV Gaussian denoising
+N = 300
+
+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 = 'gaussian', mean=0, var = 0.05, seed=10)
+noisy_data = ImageData(n1)
+
+# Regularisation Parameter
+alpha = 0.5
+
+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 = 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, 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 = 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/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
new file mode 100644
index 0000000..f2c2023
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
@@ -0,0 +1,177 @@
+# -*- 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, KullbackLeibler, \
+                      MixedL21Norm, BlockFunction
+
+from skimage.util import random_noise
+
+# Create phantom for TV Poisson denoising
+N = 300
+
+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 = '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 = 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
+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
+#%% 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
new file mode 100644
index 0000000..f96acd5
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
@@ -0,0 +1,177 @@
+# -*- 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, L1Norm, \
+                      MixedL21Norm, BlockFunction
+
+from skimage.util import random_noise
+
+# Create phantom for TV Salt & Pepper denoising
+N = 300
+
+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/sinogram_demo_tomography2D.npy b/Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy
new file mode 100644
index 0000000..f37fd4b
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/sinogram_demo_tomography2D.npy
diff --git a/Wrappers/Python/wip/Demos/untitled4.py b/Wrappers/Python/wip/Demos/untitled4.py
new file mode 100644
index 0000000..0cacbd7
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/untitled4.py
@@ -0,0 +1,87 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 24 14:21:08 2019
+
+@author: evangelos
+"""
+
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
+import numpy  
+import numpy as np                        
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG, PDHG_old
+
+from ccpi.optimisation.operators import BlockOperator, Gradient
+from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
+                      MixedL21Norm, BlockFunction
+
+from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple
+from skimage.util import random_noise
+from timeit import default_timer as timer
+
+#N = 75
+#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)
+
+N = 75
+#x = np.zeros((N,N))
+
+vert = 4
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert)
+data = ig.allocate()
+Phantom = data
+# Populate image data by looping over and filling slices
+i = 0
+while i < vert:
+    if vert > 1:
+        x = Phantom.subset(vertical=i).array
+    else:
+        x = Phantom.array
+    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)] = 0.98
+    if vert > 1 :
+        Phantom.fill(x, vertical=i)
+    i += 1
+    
+angles_num = 100
+det_w = 1.0
+det_num = N
+
+angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\
+             180/np.pi
+
+# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, 
+#         horz detector pixel size, vert detector pixel count, 
+#         vert detector pixel size.
+ag = AcquisitionGeometry('parallel',
+                         '3D',
+                         angles,
+                         N, 
+                         det_w,
+                         vert,
+                         det_w)
+
+sino = numpy.load("sinogram_demo_tomography2D.npy", mmap_mode='r')
+plt.imshow(sino)
+plt.title('Sinogram CCPi')
+plt.colorbar()
+plt.show()
+             
+#%%
+Aop = AstraProjector3DSimple(ig, ag)
+sin = Aop.direct(data)
+
+plt.imshow(sin.as_array())
+
+plt.title('Sinogram Astra')
+plt.colorbar()
+plt.show()
+
+
+
+#%%
\ No newline at end of file