volume denoising complete

2 files changed, 331 insertions, 34 deletions

diff --git a/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py b/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
index 93b0cef..99b9fe8 100644
--- a/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -164,7 +164,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 = 1e-08, # tolerance to stop inner 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
index 07e3133..e128127 100644
--- a/demos/SoftwareX_supp/Demo_VolumeDenoise.py
+++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py
@@ -25,12 +25,12 @@ 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, NDF, Diff4th
+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 = 256 # Define phantom dimensions using a scalar value (cubic phantom)
+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)
@@ -66,9 +66,9 @@ print ("#############ROF TV CPU####################")
 # set parameters
 pars = {'algorithm': ROF_TV, \
         'input' : phantom_noise,\
-        'regularisation_parameter':0.02,\
+        'regularisation_parameter':0.06,\
         'number_of_iterations': 1000,\
-        'time_marching_parameter': 0.001,\
+        'time_marching_parameter': 0.00025,\
         'tolerance_constant':0.0}
 
 tic=timeit.default_timer()
@@ -85,7 +85,7 @@ Qtools = QualityTools(phantom_tm, rof_cpu3D)
 RMSE_rof = Qtools.rmse()
 
 # SSIM measure
-Qtools = QualityTools(phantom_tm[128,:,:]*255, rof_cpu3D[128,:,:]*235)
+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)
@@ -97,9 +97,9 @@ print ("#############ROF TV GPU####################")
 pars = {'algorithm': ROF_TV, \
         'input' : phantom_noise,\
         'regularisation_parameter':0.06,\
-        'number_of_iterations': 10000,\
+        'number_of_iterations': 8330,\
         'time_marching_parameter': 0.00025,\
-        'tolerance_constant':1e-06}
+        'tolerance_constant':0.0}
 
 tic=timeit.default_timer()
 (rof_gpu3D, infogpu) = ROF_TV(pars['input'],
@@ -114,38 +114,21 @@ Run_time_rof = toc - tic
 Qtools = QualityTools(phantom_tm, rof_gpu3D)
 RMSE_rof = Qtools.rmse()
 
-sliceNo = 128
 # SSIM measure
-Qtools = QualityTools(phantom_tm[sliceNo,:,:]*255, rof_gpu3D[sliceNo,:,:]*235)
+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)
 
-sliceSel = int(0.5*N_size)
-#plt.gray()
-plt.figure() 
-plt.subplot(131)
-plt.imshow(rof_gpu3D[sliceSel,:,:],vmin=0, vmax=1.4)
-plt.title('3D ROF-TV, axial view')
-
-plt.subplot(132)
-plt.imshow(rof_gpu3D[:,sliceSel,:],vmin=0, vmax=1.4)
-plt.title('3D ROF-TV, coronal view')
-
-plt.subplot(133)
-plt.imshow(rof_gpu3D[:,:,sliceSel],vmin=0, vmax=1.4)
-plt.title('3D ROF-TV, sagittal view')
-plt.show()
-
 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.05, \
-        'number_of_iterations' :100 ,\
-        'tolerance_constant':1e-04,\
+        'regularisation_parameter':0.06, \
+        'number_of_iterations' : 930 ,\
+        'tolerance_constant':0.0,\
         'methodTV': 0 ,\
         'nonneg': 0}
 
@@ -163,7 +146,7 @@ Qtools = QualityTools(phantom_tm, fgp_cpu3D)
 RMSE_rof = Qtools.rmse()
 
 # SSIM measure
-Qtools = QualityTools(phantom_tm[128,:,:]*255, fgp_cpu3D[128,:,:]*235)
+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)
@@ -174,9 +157,9 @@ print ("#############FGP TV GPU####################")
 # set parameters
 pars = {'algorithm' : FGP_TV, \
         'input' : phantom_noise,\
-        'regularisation_parameter':0.05, \
-        'number_of_iterations' :1500 ,\
-        'tolerance_constant':1e-06,\
+        'regularisation_parameter':0.06, \
+        'number_of_iterations' :930 ,\
+        'tolerance_constant':0.0,\
         'methodTV': 0 ,\
         'nonneg': 0}
 
@@ -194,13 +177,327 @@ Qtools = QualityTools(phantom_tm, fgp_gpu3D)
 RMSE_rof = Qtools.rmse()
 
 # SSIM measure
-Qtools = QualityTools(phantom_tm[128,:,:]*255, fgp_gpu3D[128,:,:]*235)
+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))