summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py2
-rw-r--r--demos/SoftwareX_supp/Demo_VolumeDenoise.py363
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))
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 06:55:22 +0000