Merge pull request #110 from vais-ral/tol

Tolerance-based stopping criterion, fixes for a new structure, new demos

4 files changed, 516 insertions, 14 deletions

diff --git a/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py b/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
index 01491d9..5991989 100644
--- a/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
+++ b/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
@@ -1,15 +1,15 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-This demo scripts support the following publication: 
-"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with 
+This demo scripts support the following publication:
+"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
 proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
  Philip J. Withers; Software X, 2019
 ____________________________________________________________________________
 * Reads real tomographic data (stored at Zenodo)
 --- https://doi.org/10.5281/zenodo.2578893
 * Reconstructs using TomoRec software
-* Saves reconstructed images 
+* Saves reconstructed images
 ____________________________________________________________________________
 >>>>> Dependencies: <<<<<
 1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox
@@ -40,7 +40,7 @@ data_norm = normaliser(dataRaw, flats, darks, log='log')
 del dataRaw, darks, flats
 
 intens_max = 2.3
-plt.figure() 
+plt.figure()
 plt.subplot(131)
 plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max)
 plt.title('2D Projection (analytical)')
@@ -72,7 +72,7 @@ FBPrec = RectoolsDIR.FBP(data_norm[0:100,:,det_y_crop])
 
 sliceSel = 50
 max_val = 0.003
-plt.figure() 
+plt.figure()
 plt.subplot(131)
 plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray")
 plt.title('FBP Reconstruction, axial view')
@@ -108,10 +108,10 @@ RectoolsIR = RecToolsIR(DetectorsDimH =  np.size(det_y_crop),  # DetectorsDimH #
                     DetectorsDimV = 100,  # DetectorsDimV # detector dimension (vertical) for 3D case only
                     AnglesVec = angles_rad, # array of angles in radians
                     ObjSize = N_size, # a scalar to define reconstructed object dimensions
-                    datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
+                    datafidelity='LS',# data fidelity, choose LS, PWLS, GH (wip), Students t (wip)
                     nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
                     OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
-                    tolerance = 1e-08, # tolerance to stop outer iterations earlier
+                    tolerance = 0.0, # tolerance to stop inner (regularisation) iterations earlier
                     device='gpu')
 #%%
 print ("Reconstructing with ADMM method using SB-TV penalty")
@@ -124,7 +124,7 @@ RecADMM_reg_sbtv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
 
 sliceSel = 50
 max_val = 0.003
-plt.figure() 
+plt.figure()
 plt.subplot(131)
 plt.imshow(RecADMM_reg_sbtv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray")
 plt.title('3D ADMM-SB-TV Reconstruction, axial view')
@@ -164,7 +164,7 @@ RecADMM_reg_rofllt = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
 
 sliceSel = 50
 max_val = 0.003
-plt.figure() 
+plt.figure()
 plt.subplot(131)
 plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val)
 plt.title('3D ADMM-ROFLLT Reconstruction, axial view')
@@ -202,7 +202,7 @@ RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
 
 sliceSel = 50
 max_val = 0.003
-plt.figure() 
+plt.figure()
 plt.subplot(131)
 plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val)
 plt.title('3D ADMM-TGV Reconstruction, axial view')
@@ -228,4 +228,4 @@ for i in range(0,np.size(RecADMM_reg_tgv,0)):
 # Saving recpnstructed data with a unique time label
 np.save('Dendr_ADMM_TGV'+str(time_label)+'.npy', RecADMM_reg_tgv)
 del RecADMM_reg_tgv
-#%%
\ No newline at end of file
+#%%
diff --git a/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py b/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
index 59ffc0e..be99afe 100644
--- a/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
+++ b/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
@@ -77,7 +77,7 @@ RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector d
                     datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
                     nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
                     OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
-                    tolerance = 1e-08, # tolerance to stop outer iterations earlier
+                    tolerance = 0.0, # tolerance to stop inner (regularisation) iterations earlier
                     device='gpu')
 #%%
 param_space = 30
diff --git a/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py b/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
index 93b0cef..ae2bfba 100644
--- a/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -78,7 +78,6 @@ plt.title('3D Phantom, coronal (Y-Z) view')
 plt.subplot(133)
 plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1, cmap="PuOr")
 plt.title('3D Phantom, sagittal view')
-
 """
 plt.show()
 #%%
@@ -164,7 +163,7 @@ RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector d
                     datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
                     nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
                     OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
-                    tolerance = 1e-08, # tolerance to stop outer iterations earlier
+                    tolerance = 0.0, # tolerance to stop inner (regularisation) iterations earlier
                     device='gpu')
 #%%
 print ("Reconstructing with ADMM method using SB-TV penalty")
diff --git a/demos/SoftwareX_supp/Demo_VolumeDenoise.py b/demos/SoftwareX_supp/Demo_VolumeDenoise.py
new file mode 100644
index 0000000..e128127
--- /dev/null
+++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py
@@ -0,0 +1,503 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This demo scripts support the following publication: 
+"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with 
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
+ Philip J. Withers; Software X, 2019
+____________________________________________________________________________
+* Generates phantom using TomoPhantom software 
+* Denoise using closely to optimal parameters
+____________________________________________________________________________
+>>>>> Dependencies: <<<<<
+1. TomoPhantom software for phantom and data generation
+
+@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
+Apache 2.0.
+"""
+import timeit
+import matplotlib.pyplot as plt
+# import matplotlib.gridspec as gridspec
+import numpy as np
+import os
+import tomophantom
+from tomophantom import TomoP3D
+from tomophantom.supp.artifacts import ArtifactsClass
+from ccpi.supp.qualitymetrics import QualityTools
+from scipy.signal import gaussian
+from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, LLT_ROF, TGV, NDF, Diff4th
+#%%
+print ("Building 3D phantom using TomoPhantom software")
+tic=timeit.default_timer()
+model = 16 # select a model number from the library
+N_size = 128 # Define phantom dimensions using a scalar value (cubic phantom)
+path = os.path.dirname(tomophantom.__file__)
+path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
+#This will generate a N_size x N_size x N_size phantom (3D)
+phantom_tm = TomoP3D.Model(model, N_size, path_library3D)
+toc=timeit.default_timer()
+Run_time = toc - tic
+print("Phantom has been built in {} seconds".format(Run_time))
+
+# adding normally distributed noise
+artifacts_add = ArtifactsClass(phantom_tm)
+phantom_noise = artifacts_add.noise(sigma=0.1,noisetype='Gaussian')
+
+sliceSel = int(0.5*N_size)
+#plt.gray()
+plt.figure() 
+plt.subplot(131)
+plt.imshow(phantom_noise[sliceSel,:,:],vmin=0, vmax=1.4)
+plt.title('3D Phantom, axial view')
+
+plt.subplot(132)
+plt.imshow(phantom_noise[:,sliceSel,:],vmin=0, vmax=1.4)
+plt.title('3D Phantom, coronal view')
+
+plt.subplot(133)
+plt.imshow(phantom_noise[:,:,sliceSel],vmin=0, vmax=1.4)
+plt.title('3D Phantom, sagittal view')
+plt.show()
+#%%
+print ("____________________Applying regularisers_______________________")
+print ("________________________________________________________________")
+
+print ("#############ROF TV CPU####################")
+# set parameters
+pars = {'algorithm': ROF_TV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06,\
+        'number_of_iterations': 1000,\
+        'time_marching_parameter': 0.00025,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(rof_cpu3D, infcpu) = ROF_TV(pars['input'],
+             pars['regularisation_parameter'],
+             pars['number_of_iterations'],
+             pars['time_marching_parameter'],
+             pars['tolerance_constant'],'cpu')
+
+toc=timeit.default_timer()
+
+Run_time_rof = toc - tic
+Qtools = QualityTools(phantom_tm, rof_cpu3D)
+RMSE_rof = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, rof_cpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_rof = Qtools.ssim(win2d)
+
+print("ROF-TV (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_rof[0],Run_time_rof))
+#%%
+print ("#############ROF TV GPU####################")
+# set parameters
+pars = {'algorithm': ROF_TV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06,\
+        'number_of_iterations': 8330,\
+        'time_marching_parameter': 0.00025,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(rof_gpu3D, infogpu) = ROF_TV(pars['input'],
+             pars['regularisation_parameter'],
+             pars['number_of_iterations'],
+             pars['time_marching_parameter'],
+             pars['tolerance_constant'],'gpu')
+
+toc=timeit.default_timer()
+
+Run_time_rof = toc - tic
+Qtools = QualityTools(phantom_tm, rof_gpu3D)
+RMSE_rof = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, rof_gpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_rof = Qtools.ssim(win2d)
+
+print("ROF-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_rof[0],Run_time_rof))
+#%%
+print ("#############FGP TV CPU####################")
+# set parameters
+pars =  {'algorithm' : FGP_TV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06, \
+        'number_of_iterations' : 930 ,\
+        'tolerance_constant':0.0,\
+        'methodTV': 0 ,\
+        'nonneg': 0}
+
+tic=timeit.default_timer()
+(fgp_cpu3D, infoFGP) = FGP_TV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'],
+              pars['nonneg'],'cpu')
+toc=timeit.default_timer()
+
+Run_time_fgp = toc - tic
+Qtools = QualityTools(phantom_tm, fgp_cpu3D)
+RMSE_rof = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, fgp_cpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_fgp = Qtools.ssim(win2d)
+
+print("FGP-TV (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_fgp[0],Run_time_fgp))
+#%%
+print ("#############FGP TV GPU####################")
+# set parameters
+pars = {'algorithm' : FGP_TV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06, \
+        'number_of_iterations' :930 ,\
+        'tolerance_constant':0.0,\
+        'methodTV': 0 ,\
+        'nonneg': 0}
+
+tic=timeit.default_timer()
+(fgp_gpu3D,infogpu)  = FGP_TV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'],
+              pars['nonneg'],'gpu')
+toc=timeit.default_timer()
+
+Run_time_fgp = toc - tic
+Qtools = QualityTools(phantom_tm, fgp_gpu3D)
+RMSE_rof = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, fgp_gpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_fgp = Qtools.ssim(win2d)
+
+print("FGP-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_fgp[0],Run_time_fgp))
+#%%
+print ("#############SB TV CPU####################")
+# set parameters
+pars =  {'algorithm' : SB_TV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06, \
+        'number_of_iterations' :225 ,\
+        'tolerance_constant':0.0,\
+        'methodTV': 0}
+
+tic=timeit.default_timer()
+(sb_cpu3D, info_vec_cpu) = SB_TV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'], 'cpu')
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, sb_cpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, sb_cpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("SB-TV (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############SB TV GPU####################")
+# set parameters
+pars = {'algorithm' : SB_TV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06, \
+        'number_of_iterations' :225 ,\
+        'tolerance_constant':0.0,\
+        'methodTV': 0}
+
+tic=timeit.default_timer()
+(sb_gpu3D,info_vec_gpu) = SB_TV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'], 'gpu')
+
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, sb_gpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, sb_gpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("SB-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############NDF CPU####################")
+# set parameters
+pars = {'algorithm' : NDF, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06, \
+        'edge_parameter':0.017,\
+        'number_of_iterations' :530 ,\
+        'time_marching_parameter':0.01,\
+        'penalty_type':1,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(ndf_cpu3D, info_vec_cpu) = NDF(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['edge_parameter'], 
+              pars['number_of_iterations'],
+              pars['time_marching_parameter'], 
+              pars['penalty_type'],
+              pars['tolerance_constant'],'cpu')
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, ndf_cpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, ndf_cpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("NDF (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############NDF GPU####################")
+# set parameters
+pars = {'algorithm' : NDF, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06, \
+        'edge_parameter':0.017,\
+        'number_of_iterations' :530 ,\
+        'time_marching_parameter':0.01,\
+        'penalty_type':1,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(ndf_gpu3D,info_vec_gpu) = NDF(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['edge_parameter'], 
+              pars['number_of_iterations'],
+              pars['time_marching_parameter'], 
+              pars['penalty_type'],
+              pars['tolerance_constant'],'gpu')
+
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, ndf_gpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, ndf_gpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("NDF (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############Diff4th CPU####################")
+# set parameters
+pars = {'algorithm' : Diff4th, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':4.5, \
+        'edge_parameter':0.035,\
+        'number_of_iterations' :2425 ,\
+        'time_marching_parameter':0.001,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(diff4th_cpu3D, info_vec_cpu) = Diff4th(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['edge_parameter'], 
+              pars['number_of_iterations'],
+              pars['time_marching_parameter'],
+              pars['tolerance_constant'],'cpu')
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, diff4th_cpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, diff4th_cpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("Diff4th (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############Diff4th GPU####################")
+# set parameters
+pars = {'algorithm' : Diff4th, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':4.5, \
+        'edge_parameter':0.035,\
+        'number_of_iterations' :2425 ,\
+        'time_marching_parameter':0.001,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(diff4th_gpu3D,info_vec_gpu) =  Diff4th(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['edge_parameter'], 
+              pars['number_of_iterations'],
+              pars['time_marching_parameter'],
+              pars['tolerance_constant'],'gpu')
+
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, diff4th_gpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, diff4th_gpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("Diff4th (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############TGV CPU####################")
+# set parameters
+pars = {'algorithm' : TGV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06,\
+        'alpha1':1.0,\
+        'alpha0':2.0,\
+        'number_of_iterations' :1000,\
+        'LipshitzConstant' :12,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(tgv_cpu3D, info_vec_cpu) = TGV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['alpha1'],
+              pars['alpha0'],
+              pars['number_of_iterations'],
+              pars['LipshitzConstant'],
+              pars['tolerance_constant'],'cpu')
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, tgv_cpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, tgv_cpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("TGV (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############TGV GPU####################")
+# set parameters
+pars = {'algorithm' : TGV, \
+        'input' : phantom_noise,\
+        'regularisation_parameter':0.06,\
+        'alpha1':1.0,\
+        'alpha0':2.0,\
+        'number_of_iterations' :7845,\
+        'LipshitzConstant' :12,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(tgv_gpu3D,info_vec_gpu) =  TGV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['alpha1'],
+              pars['alpha0'],
+              pars['number_of_iterations'],
+              pars['LipshitzConstant'],
+              pars['tolerance_constant'],'gpu')
+
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, tgv_gpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, tgv_gpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("TGV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############ROF-LLT CPU####################")
+# set parameters
+pars = {'algorithm' : LLT_ROF, \
+        'input' : phantom_noise,\
+        'regularisation_parameterROF':0.03, \
+        'regularisation_parameterLLT':0.015, \
+        'number_of_iterations' : 1000 ,\
+        'time_marching_parameter' :0.00025 ,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(rofllt_cpu3D, info_vec_cpu) = LLT_ROF(pars['input'], 
+              pars['regularisation_parameterROF'],
+              pars['regularisation_parameterLLT'],
+              pars['number_of_iterations'],
+              pars['time_marching_parameter'],
+              pars['tolerance_constant'], 'cpu')
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, rofllt_cpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, rofllt_cpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("ROF-LLT  (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))
+#%%
+print ("#############ROF-LLT  GPU####################")
+# set parameters
+pars = {'algorithm' : LLT_ROF, \
+        'input' : phantom_noise,\
+        'regularisation_parameterROF':0.03, \
+        'regularisation_parameterLLT':0.015, \
+        'number_of_iterations' : 8000 ,\
+        'time_marching_parameter' :0.00025 ,\
+        'tolerance_constant':0.0}
+
+tic=timeit.default_timer()
+(rofllt_gpu3D,info_vec_gpu) = LLT_ROF(pars['input'], 
+              pars['regularisation_parameterROF'],
+              pars['regularisation_parameterLLT'],
+              pars['number_of_iterations'],
+              pars['time_marching_parameter'],
+              pars['tolerance_constant'], 'gpu')
+toc=timeit.default_timer()
+
+Run_time = toc - tic
+Qtools = QualityTools(phantom_tm, rofllt_gpu3D)
+RMSE = Qtools.rmse()
+
+# SSIM measure
+Qtools = QualityTools(phantom_tm[sliceSel,:,:]*255, rofllt_gpu3D[sliceSel,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim = Qtools.ssim(win2d)
+
+print("ROF-LLT  (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE,ssim[0],Run_time))