fix ROF_TV (#28)

* fix ROF_TV

* fix ROF_TV

*  add TGV proximal

* add new regs

* added FGP_dTV

* fix out

* remove test to wip

3 files changed, 499 insertions, 0 deletions

diff --git a/Wrappers/Python/wip/test_LLT_ROF_framework_ccpi_regulariser.py b/Wrappers/Python/wip/test_LLT_ROF_framework_ccpi_regulariser.py
new file mode 100644
index 0000000..3368732
--- /dev/null
+++ b/Wrappers/Python/wip/test_LLT_ROF_framework_ccpi_regulariser.py
@@ -0,0 +1,179 @@
+from ccpi.framework import ImageData, TestData
+
+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, KullbackLeibler, L2NormSquared
+import os
+import sys
+from ccpi.plugins.regularisers import TGV, LLT_ROF
+
+# user supplied input
+if len(sys.argv) > 1:
+    which_noise = int(sys.argv[1])
+else:
+    which_noise = 0
+print ("Applying {} noise")
+
+if len(sys.argv) > 2:
+    method = sys.argv[2]
+else:
+    method = '0'
+print ("method ", method)
+
+
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+data = loader.load(TestData.SHAPES)
+ig = data.geometry
+ag = ig
+
+# Create noisy data. 
+noises = ['gaussian', 'poisson', 's&p']
+noise = noises[which_noise]
+if noise == 's&p':
+    n1 = TestData.random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2, seed=10)
+elif noise == 'poisson':
+    scale = 5
+    n1 = TestData.random_noise(data.as_array()/scale, mode = noise, seed = 10)*scale
+elif noise == 'gaussian':
+    n1 = TestData.random_noise(data.as_array(), mode = noise, seed = 10)
+else:
+    raise ValueError('Unsupported Noise ', noise)
+noisy_data = ImageData(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 depending on the noise distribution
+if noise == 's&p':
+    alpha = 0.8
+elif noise == 'poisson':
+    alpha = .3
+elif noise == 'gaussian':
+    alpha = .2
+
+beta = 2 * alpha
+
+# Fidelity
+if noise == 's&p':
+    f3 = L1Norm(b=noisy_data)
+elif noise == 'poisson':
+    f3 = KullbackLeibler(noisy_data)
+elif noise == 'gaussian':
+    f3 = 0.5 * L2NormSquared(b=noisy_data)
+
+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() 
+    
+    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(f3, ZeroFunction())
+     
+# Compute operator Norm
+normK = operator.norm()
+
+# Primal & dual stepsizes
+#sigma = 1/normK
+#tau = 1/normK
+
+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(5000)
+
+# Show results
+plt.figure(figsize=(20,5))
+plt.subplot(1,4,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(1,4,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(1,4,3)
+plt.imshow(pdhg.get_output()[0].as_array())
+plt.title('TGV Reconstruction')
+plt.colorbar()
+plt.subplot(1,4,4)
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'TGV reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
+
+#%%  Run CCPi-regulariser
+# The TGV implementation is using PDHG algorithm with fixed 
+# sigma = tau = 1/sqrt(12)
+
+# There is an early stopping criteria 
+# https://github.com/vais-ral/CCPi-Regularisation-Toolkit/blob/master/src/Core/regularisers_CPU/TGV_core.c#L168
+
+g = TGV(1, alpha, beta, 2000, normK**2, 1e-6, 'gpu')
+#alphaROF = 0.1
+#alphaLLT = 0.05
+#g = LLT_ROF(alphaROF, alphaLLT, 500, 0.001, 1e-6, 'gpu')
+sol = g.proximal(noisy_data, 1)
+
+plt.imshow(sol.as_array())
+plt.show()
+
+plt.imshow(pdhg.get_output()[0].as_array())
+plt.show()
+
+plt.imshow(np.abs(sol.as_array() - pdhg.get_output()[0].as_array()))
+plt.colorbar()
+plt.show()
+
+#%%
+
+plt.plot(np.linspace(0,299,300),sol.as_array()[100,:], np.linspace(0,299,300),pdhg.get_output()[0].as_array()[100,:])
+
diff --git a/Wrappers/Python/wip/test_TGV_framework_ccpi_regulariser.py b/Wrappers/Python/wip/test_TGV_framework_ccpi_regulariser.py
new file mode 100644
index 0000000..a6f6094
--- /dev/null
+++ b/Wrappers/Python/wip/test_TGV_framework_ccpi_regulariser.py
@@ -0,0 +1,173 @@
+from ccpi.framework import ImageData, TestData
+
+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, KullbackLeibler, L2NormSquared
+import os
+import sys
+from ccpi.plugins.regularisers import TGV
+
+# user supplied input
+if len(sys.argv) > 1:
+    which_noise = int(sys.argv[1])
+else:
+    which_noise = 0
+print ("Applying {} noise")
+
+if len(sys.argv) > 2:
+    method = sys.argv[2]
+else:
+    method = '1'
+print ("method ", method)
+
+
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+data = loader.load(TestData.SHAPES)
+ig = data.geometry
+ag = ig
+
+# Create noisy data. 
+noises = ['gaussian', 'poisson', 's&p']
+noise = noises[which_noise]
+if noise == 's&p':
+    n1 = TestData.random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2, seed=10)
+elif noise == 'poisson':
+    scale = 5
+    n1 = TestData.random_noise(data.as_array()/scale, mode = noise, seed = 10)*scale
+elif noise == 'gaussian':
+    n1 = TestData.random_noise(data.as_array(), mode = noise, seed = 10)
+else:
+    raise ValueError('Unsupported Noise ', noise)
+noisy_data = ImageData(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 depending on the noise distribution
+if noise == 's&p':
+    alpha = 0.8
+elif noise == 'poisson':
+    alpha = .3
+elif noise == 'gaussian':
+    alpha = .2
+
+beta = 2 * alpha
+
+# Fidelity
+if noise == 's&p':
+    f3 = L1Norm(b=noisy_data)
+elif noise == 'poisson':
+    f3 = KullbackLeibler(noisy_data)
+elif noise == 'gaussian':
+    f3 = 0.5 * L2NormSquared(b=noisy_data)
+
+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() 
+    
+    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(f3, ZeroFunction())
+     
+# Compute operator Norm
+normK = operator.norm()
+
+# Primal & dual stepsizes
+sigma = 1/normK
+tau = 1/normK
+
+# 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(5000)
+
+# Show results
+plt.figure(figsize=(20,5))
+plt.subplot(1,4,1)
+plt.imshow(data.subset(channel=0).as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(1,4,2)
+plt.imshow(noisy_data.subset(channel=0).as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(1,4,3)
+plt.imshow(pdhg.get_output()[0].as_array())
+plt.title('TGV Reconstruction')
+plt.colorbar()
+plt.subplot(1,4,4)
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'TGV reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
+
+#  Run CCPi-regulariser
+# The TGV implementation is using PDHG algorithm with fixed 
+# sigma = tau = 1/sqrt(12)
+
+# There is an early stopping criteria 
+# https://github.com/vais-ral/CCPi-Regularisation-Toolkit/blob/master/src/Core/regularisers_CPU/TGV_core.c#L168
+
+g = TGV(1, alpha, beta, 5000, 12, 1e-6, 'gpu')
+sol = g.proximal(noisy_data, 1)
+
+plt.imshow(sol.as_array())
+plt.show()
+
+plt.imshow(pdhg.get_output()[0].as_array())
+plt.show()
+
+plt.imshow(np.abs(sol.as_array() - pdhg.get_output()[0].as_array()))
+plt.colorbar()
+plt.show()
+
+#%%
+
+plt.plot(np.linspace(0,299,300),sol.as_array()[100,:], np.linspace(0,299,300),pdhg.get_output()[0].as_array()[100,:])
+
diff --git a/Wrappers/Python/wip/test_TVN_framework_ccpi_regulariser.py b/Wrappers/Python/wip/test_TVN_framework_ccpi_regulariser.py
new file mode 100644
index 0000000..a0c893d
--- /dev/null
+++ b/Wrappers/Python/wip/test_TVN_framework_ccpi_regulariser.py
@@ -0,0 +1,147 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb 22 11:39:43 2018
+
+Demonstration of CPU regularisers 
+
+@authors: Daniil Kazantsev, Edoardo Pasca
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import timeit
+from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, Diff4th
+from ccpi.filters.regularisers import PatchSelect, NLTV
+from ccpi.supp.qualitymetrics import QualityTools
+
+from ccpi.plugins.regularisers import TNV as TNV_new
+###############################################################################
+def printParametersToString(pars):
+        txt = r''
+        for key, value in pars.items():
+            if key== 'algorithm' :
+                txt += "{0} = {1}".format(key, value.__name__)
+            elif key == 'input':
+                txt += "{0} = {1}".format(key, np.shape(value))
+            elif key == 'refdata':
+                txt += "{0} = {1}".format(key, np.shape(value))
+            else:
+                txt += "{0} = {1}".format(key, value)
+            txt += '\n'
+        return txt
+###############################################################################
+
+filename = 'lena_gray_512.tif'
+
+# read image
+Im = plt.imread(filename)
+Im = np.asarray(Im, dtype='float32')
+
+Im = Im/255.0
+perc = 0.05
+u0 = Im + np.random.normal(loc = 0 ,
+                                  scale = perc * Im , 
+                                  size = np.shape(Im))
+u_ref = Im + np.random.normal(loc = 0 ,
+                                  scale = 0.01 * Im , 
+                                  size = np.shape(Im))
+(N,M) = np.shape(u0)
+# map the u0 u0->u0>0
+# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = u0.astype('float32')
+u_ref = u_ref.astype('float32')
+
+# change dims to check that modules work with non-squared images
+"""
+M = M-100
+u_ref2 = np.zeros([N,M],dtype='float32')
+u_ref2[:,0:M] = u_ref[:,0:M]
+u_ref = u_ref2
+del u_ref2
+
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
+u0 = u02
+del u02
+
+Im2 = np.zeros([N,M],dtype='float32')
+Im2[:,0:M] = Im[:,0:M]
+Im = Im2
+del Im2
+"""
+
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("__________Total nuclear Variation__________")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot 
+fig = plt.figure()
+plt.suptitle('Performance of TNV regulariser using the CPU')
+a=fig.add_subplot(1,2,1)
+a.set_title('Noisy Image')
+imgplot = plt.imshow(u0,cmap="gray")
+
+channelsNo = 5
+noisyVol = np.zeros((channelsNo,N,M),dtype='float32')
+idealVol = np.zeros((channelsNo,N,M),dtype='float32')
+
+for i in range (channelsNo):
+    noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
+    idealVol[i,:,:] = Im
+
+# set parameters
+pars = {'algorithm' : TNV, \
+        'input' : noisyVol,\
+        'regularisation_parameter': 0.04, \
+        'number_of_iterations' : 200 ,\
+        'tolerance_constant':1e-05
+        }
+        
+print ("#############TNV CPU#################")
+start_time = timeit.default_timer()
+tnv_cpu = TNV(pars['input'],           
+              pars['regularisation_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'])
+
+Qtools = QualityTools(idealVol, tnv_cpu)
+pars['rmse'] = Qtools.rmse()
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,2,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(tnv_cpu[3,:,:], cmap="gray")
+plt.title('{}'.format('CPU results'))
+
+#%%
+from ccpi.framework import ImageData
+g = TNV_new(0.04, 200, 1e-5)
+sol = g.proximal(ImageData(noisyVol), 1)
+
+#%%
+plt.figure(figsize=(10,10))
+plt.subplot(3,1,1)
+plt.imshow(sol.as_array()[2])
+plt.colorbar()
+
+plt.subplot(3,1,2)
+plt.imshow(tnv_cpu[2])
+plt.colorbar()
+
+plt.subplot(3,1,3)
+plt.imshow(np.abs(tnv_cpu[2] - sol.as_array()[2]))
+plt.colorbar()
+
+plt.show()
+
+#plt.imshow(sol.as_array())
\ No newline at end of file