From cd2844b6384cb74fae5ee3f7f0fdb91be752653e Mon Sep 17 00:00:00 2001 From: dkazanc Date: Mon, 11 Mar 2019 16:43:05 +0000 Subject: volume denoising complete --- demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py | 2 +- demos/SoftwareX_supp/Demo_VolumeDenoise.py | 363 +++++++++++++++++++++--- 2 files changed, 331 insertions(+), 34 deletions(-) (limited to 'demos') 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)) -- cgit v1.2.3