add denoising demos

4 files changed, 70 insertions, 72 deletions

diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py
index dffb6d8..7b65c31 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py
@@ -23,7 +23,7 @@ from skimage.util import random_noise
 
 # Create phantom for TGV SaltPepper denoising
 
-N = 300
+N = 100
 
 data = np.zeros((N,N))
 
@@ -47,7 +47,7 @@ noisy_data = ImageData(n1)
 alpha = 0.8
 beta = numpy.sqrt(2)* alpha
 
-method = '0'
+method = '1'
 
 if method == '0':
 
@@ -83,7 +83,7 @@ else:
     f2 = beta * MixedL21Norm()     
     
     f = BlockFunction(f1, f2)         
-    g = BlockFunction(0.5 * L1Norm(b=noisy_data), ZeroFunction())
+    g = BlockFunction(L1Norm(b=noisy_data), ZeroFunction())
      
 ## Compute operator Norm
 normK = operator.norm()
@@ -122,75 +122,73 @@ 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 = pnorm(u - noisy_data.as_array(),1)    
+    solver = MOSEK
+
+        # 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()[0].as_array() - u.value )
+        
+    plt.figure(figsize=(15,15))
+    plt.subplot(3,1,1)
+    plt.imshow(pdhg.get_output()[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), pdhg.get_output()[0].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]))
 
 
 
 
 
-#%% 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
index 10e5621..1a3e0df 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
@@ -48,7 +48,7 @@ noisy_data = ImageData(n1)
 # Regularisation Parameter
 alpha = 0.5
 
-method = '1'
+method = '0'
 
 if method == '0':
 
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
index f2c2023..482b3b4 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
@@ -32,7 +32,7 @@ from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
 from skimage.util import random_noise
 
 # Create phantom for TV Poisson denoising
-N = 300
+N = 100
 
 data = np.zeros((N,N))
 data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
@@ -48,7 +48,7 @@ noisy_data = ImageData(n1)
 # Regularisation Parameter
 alpha = 2
 
-method = '0'
+method = '1'
 
 if method == '0':
 
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
index f96acd5..484fe42 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
@@ -32,7 +32,7 @@ from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
 from skimage.util import random_noise
 
 # Create phantom for TV Salt & Pepper denoising
-N = 300
+N = 100
 
 data = np.zeros((N,N))
 data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
@@ -48,7 +48,7 @@ noisy_data = ImageData(n1)
 # Regularisation Parameter
 alpha = 2
 
-method = '0'
+method = '1'
 
 if method == '0':