try TGV tomo

5 files changed, 82 insertions, 77 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
index ace0ba7..f585c9b 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
@@ -60,7 +60,7 @@ class IndicatorBox(Function):
         
         if out is None:
             return (x.maximum(self.lower)).minimum(self.upper)        
-        else:                   
+        else:               
             x.maximum(self.lower, out=out)
             out.minimum(self.upper, out=out) 
             
diff --git a/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py
index 49d4db6..19cf684 100644
--- a/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py
+++ b/Wrappers/Python/demos/PDHG_TGV_Tomo2D.py
@@ -27,7 +27,7 @@ from ccpi.optimisation.algorithms import PDHG
 
 from ccpi.optimisation.operators import BlockOperator, Gradient, Identity, \
                                     SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
+from ccpi.optimisation.functions import ZeroFunction, IndicatorBox, KullbackLeibler, \
                       MixedL21Norm, BlockFunction
 
 from ccpi.astra.ops import AstraProjectorSimple
@@ -84,13 +84,16 @@ f1 = alpha * MixedL21Norm()
 f2 = beta * MixedL21Norm() 
 f3 = KullbackLeibler(noisy_data)    
 f = BlockFunction(f1, f2, f3)         
-g = ZeroFunction()
+
+g =  BlockFunction(-1 * IndicatorBox(lower=0), ZeroFunction())
+#g = IndicatorBox(lower=0)
+#g = ZeroFunction()
     
 # Compute operator Norm
 normK = operator.norm()
 
 # Primal & dual stepsizes
-sigma = 1
+sigma = 10
 tau = 1/(sigma*normK**2)
 
 
diff --git a/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py
index 63f8955..22972da 100644
--- a/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py
+++ b/Wrappers/Python/demos/PDHG_Tikhonov_Tomo2D.py
@@ -50,6 +50,8 @@ from ccpi.optimisation.operators import BlockOperator, Gradient
 from ccpi.optimisation.functions import IndicatorBox, L2NormSquared, BlockFunction
 from skimage.util import random_noise
 from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.framework import TestData
+import os, sys
 
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
 
@@ -97,9 +99,9 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
 
-#%%
+
 # Regularisation Parameter
-alpha = 100
+alpha = 500
 
 # Create operators
 op1 = Gradient(ig)
diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py
index 57f6fcd..15c0a05 100644
--- a/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py
@@ -28,7 +28,7 @@ Problem:     min_{x} \alpha * ||\nabla x - w||_{2,1} +
                      \frac{1}{2} * || x - g ||_{2}^{2}
 
              \alpha: Regularization parameter
-             \alpha: Regularization parameter
+             \beta: Regularization parameter
              
              \nabla: Gradient operator 
               E: Symmetrized Gradient operator
diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py
index 8f39f86..c44f393 100644
--- a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py
@@ -153,74 +153,74 @@ plt.title('Middle Line Profiles')
 plt.show()
 
 
-##%% Check with CVX solution
+#%% Check with CVX solution
 #
-#from ccpi.optimisation.operators import SparseFiniteDiff
-#import astra
-#import numpy
-#
-#try:
-#    from cvxpy import *
-#    cvx_not_installable = True
-#except ImportError:
-#    cvx_not_installable = False
-#
-#
-#if cvx_not_installable:
-#    
-#
-#    ##Construct problem    
-#    u = Variable(N*N)
-#    q = Variable()
-#    
-#    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-#    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-#
-#    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-#    
-#    # create matrix representation for Astra operator
-#
-#    vol_geom = astra.create_vol_geom(N, N)
-#    proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
-#
-#    proj_id = astra.create_projector('strip', proj_geom, vol_geom)
-#
-#    matrix_id = astra.projector.matrix(proj_id)
-#
-#    ProjMat = astra.matrix.get(matrix_id)
-#    
-#    tmp = noisy_data.as_array().ravel()
-#    
-#    fidelity = sum(kl_div(tmp, ProjMat * u))    
-#    
-#    constraints = [q>=fidelity, u>=0]   
-#    solver = SCS
-#    obj =  Minimize( regulariser +  q)
-#    prob = Problem(obj, constraints)
-#    result = prob.solve(verbose = True, solver = solver)    
-#             
-#    diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N)))
-#           
-#    plt.figure(figsize=(15,15))
-#    plt.subplot(3,1,1)
-#    plt.imshow(pdhg.get_output().as_array())
-#    plt.title('PDHG solution')
-#    plt.colorbar()
-#    plt.subplot(3,1,2)
-#    plt.imshow(np.reshape(u.value, (N, N)))
-#    plt.title('CVX solution')
-#    plt.colorbar()        
-#    plt.subplot(3,1,3)
-#    plt.imshow(diff_cvx)
-#    plt.title('Difference')
-#    plt.colorbar()
-#    plt.show()    
-#    
-#    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-#    plt.plot(np.linspace(0,N,N), np.reshape(u.value, (N, N))[int(N/2),:], label = 'CVX')
-#    plt.legend()
-#    plt.title('Middle Line Profiles')
-#    plt.show()
-#            
-#    print('Primal Objective (CVX) {} '.format(obj.value))
-#    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file
+from ccpi.optimisation.operators import SparseFiniteDiff
+import astra
+import numpy
+
+try:
+    from cvxpy import *
+    cvx_not_installable = True
+except ImportError:
+    cvx_not_installable = False
+
+
+if cvx_not_installable:
+    
+
+    ##Construct problem    
+    u = Variable(N*N)
+    q = Variable()
+    
+    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+
+    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+    
+    # create matrix representation for Astra operator
+
+    vol_geom = astra.create_vol_geom(N, N)
+    proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
+
+    proj_id = astra.create_projector('strip', proj_geom, vol_geom)
+
+    matrix_id = astra.projector.matrix(proj_id)
+
+    ProjMat = astra.matrix.get(matrix_id)
+    
+    tmp = noisy_data.as_array().ravel()
+    
+    fidelity = sum(kl_div(tmp, ProjMat * u))    
+    
+    constraints = [q>=fidelity, u>=0]   
+    solver = SCS
+    obj =  Minimize( regulariser +  q)
+    prob = Problem(obj, constraints)
+    result = prob.solve(verbose = True, solver = solver)    
+             
+    diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N)))
+           
+    plt.figure(figsize=(15,15))
+    plt.subplot(3,1,1)
+    plt.imshow(pdhg.get_output().as_array())
+    plt.title('PDHG solution')
+    plt.colorbar()
+    plt.subplot(3,1,2)
+    plt.imshow(np.reshape(u.value, (N, N)))
+    plt.title('CVX solution')
+    plt.colorbar()        
+    plt.subplot(3,1,3)
+    plt.imshow(diff_cvx)
+    plt.title('Difference')
+    plt.colorbar()
+    plt.show()    
+    
+    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
+    plt.plot(np.linspace(0,N,N), np.reshape(u.value, (N, N))[int(N/2),:], label = 'CVX')
+    plt.legend()
+    plt.title('Middle Line Profiles')
+    plt.show()
+            
+    print('Primal Objective (CVX) {} '.format(obj.value))
+    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file