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

6 files changed, 239 insertions, 130 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
index cf1a244..3096191 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
@@ -106,17 +106,27 @@ class KullbackLeibler(Function):
             z = x + tau * self.bnoise
             return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt())
         else:
-
-            z_m = x + tau * self.bnoise -1
-            self.b.multiply(4*tau, out=out)
-            z_m.multiply(z_m, out=z_m)
-            out += z_m
-
+            
+            tmp = x + tau * self.bnoise
+            self.b.multiply(4*tau, out=out)            
+            out.add((tmp-1)**2, out=out)
             out.sqrt(out=out)
-                        
             out *= -1
-            out += tmp2
+            out.add(tmp+1, out=out)
             out *= 0.5
+            
+            
+
+#            z_m = x + tau * self.bnoise -1
+#            self.b.multiply(4*tau, out=out)
+#            z_m.multiply(z_m, out=z_m)
+#            out += z_m
+#
+#            out.sqrt(out=out)
+#                        
+#            out *= -1
+#            out += tmp2
+#            out *= 0.5
 
         
     
diff --git a/Wrappers/Python/demos/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py
index 7b65c31..57f6fcd 100644
--- a/Wrappers/Python/demos/PDHG_TGV_Denoising_SaltPepper.py
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py
@@ -1,9 +1,48 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Feb 22 14:53:03 2019
-
-@author: evangelos
+#========================================================================
+# 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 Generalised Variation (TGV) Denoising using PDHG algorithm:
+
+
+Problem:     min_{x} \alpha * ||\nabla x - w||_{2,1} +
+                     \beta *  || E w ||_{2,1} +
+                     \frac{1}{2} * || x - g ||_{2}^{2}
+
+             \alpha: Regularization parameter
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator 
+              E: Symmetrized Gradient operator
+             
+             g: Noisy Data with Salt & Pepper Noise
+                          
+             Method = 0 ( PDHG - split ) :  K = [ \nabla, - Identity
+                                                  ZeroOperator, E 
+                                                  Identity, ZeroOperator]
+                          
+                                                                    
+             Method = 1 (PDHG - explicit ):  K = [ \nabla, - Identity
+                                                  ZeroOperator, E ] 
+                                                                
 """
 
 from ccpi.framework import ImageData, ImageGeometry
@@ -43,6 +82,18 @@ ag = ig
 n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
 noisy_data = ImageData(n1)
 
+# 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 Parameters
 alpha = 0.8
 beta = numpy.sqrt(2)* alpha
@@ -94,10 +145,10 @@ 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 = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
 pdhg.max_iteration = 2000
 pdhg.update_objective_interval = 50
-pdhg.run(2000)
+pdhg.run(2000, verbose = False)
 
 #%%
 plt.figure(figsize=(15,15))
diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian.py
index c830025..afdb6a2 100644
--- a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian.py
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian.py
@@ -18,14 +18,10 @@
 # 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}
 
@@ -35,12 +31,12 @@ Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}
              
              g: Noisy Data with Gaussian Noise
                           
-             Method = 0:  K = [ \nabla,
-                                 Identity]
+             Method = 0 ( PDHG - split ) :  K = [ \nabla,
+                                                 Identity]
+                          
                                                                     
-             Method = 1:  K = \nabla    
-             
-             
+             Method = 1 (PDHG - explicit ):  K = \nabla  
+                                                                
 """
 
 from ccpi.framework import ImageData, ImageGeometry
@@ -54,20 +50,17 @@ 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.data import camera
-
-
-# Load Data
-data = camera(size=(256,256))                    
-
-N, M = data.shape
-
-# Image and Acquitition Geometries
+      
+# Load Data                      
+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
 np.random.seed(10)
 noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.1, size=ig.shape) )
diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Poisson.py
index 70f6b9b..4d53635 100644
--- a/Wrappers/Python/demos/PDHG_TV_Denoising_Poisson.py
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Poisson.py
@@ -1,41 +1,43 @@
-# -*- 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.
+#========================================================================
+# 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, x>0  \alpha * ||\nabla x||_{2,1} + \int x - g * log(x)
 
-
-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    
+             \nabla: Gradient operator  
+             
+             g: Noisy Data with Poisson Noise
              
              
+             Method = 0 ( PDHG - split ) :  K = [ \nabla,
+                                                 Identity]
+                          
+                                                                    
+             Method = 1 (PDHG - explicit ):  K = \nabla     
+                       
 """
 
 from ccpi.framework import ImageData, ImageGeometry
@@ -66,10 +68,25 @@ ag = ig
 n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10)
 noisy_data = ImageData(n1)
 
+
+# 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 = '1'
+method = '0'
 
 if method == '0':
 
@@ -99,26 +116,15 @@ else:
 # Compute operator Norm
 normK = operator.norm()
 
-# Primal & dual stepsizes
+# 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)
+# 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)
 
 
 plt.figure(figsize=(15,15))
diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_SaltPepper.py
index f5d4ce4..c5709c3 100644
--- a/Wrappers/Python/demos/PDHG_TV_Denoising_SaltPepper.py
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_SaltPepper.py
@@ -1,39 +1,43 @@
-# -*- 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.
+#========================================================================
+# 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, x>0  \alpha * ||\nabla x||_{1} + ||x-g||_{1}
+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
-             \alpha: Regularization parameter
              
-             Method = 0:  K = [ \nabla,
-                                 Identity]
+             
+             Method = 0 ( PDHG - split ) :  K = [ \nabla,
+                                                 Identity]
+                          
                                                                     
-             Method = 1:  K = \nabla    
+             Method = 1 (PDHG - explicit ):  K = \nabla    
              
              
 """
@@ -66,6 +70,18 @@ ag = ig
 n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
 noisy_data = ImageData(n1)
 
+# 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
 
@@ -80,8 +96,7 @@ if method == '0':
     # Create BlockOperator
     operator = BlockOperator(op1, op2, shape=(2,1) ) 
 
-    # Create functions
-      
+    # Create functions      
     f1 = alpha * MixedL21Norm()
     f2 = L1Norm(b = noisy_data)    
     f = BlockFunction(f1, f2)  
diff --git a/Wrappers/Python/demos/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Denoising.py
index 041d4ee..7b73c1a 100644
--- a/Wrappers/Python/demos/PDHG_Tikhonov_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Denoising.py
@@ -1,21 +1,43 @@
-# -*- 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.
+#========================================================================
+# 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.
+#
+#=========================================================================
+""" 
+
+Tikhonov Denoising using PDHG algorithm:
+
+
+Problem:     min_{x} \alpha * ||\nabla x||_{2}^{2} + \frac{1}{2} * || x - g ||_{2}^{2}
+
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator 
+             
+             g: Noisy Data with Gaussian Noise
+                          
+             Method = 0 ( PDHG - split ) :  K = [ \nabla,
+                                                 Identity]
+                          
+                                                                    
+             Method = 1 (PDHG - explicit ):  K = \nabla  
+                                                                
+"""
 
 from ccpi.framework import ImageData, ImageGeometry
 
@@ -41,13 +63,25 @@ 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)
+n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
 noisy_data = ImageData(n1)
 
+# 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 = 4
 
-method = '0'
+method = '1'
 
 if method == '0':