From 5a12eb57a4965dea7241093c1fe7bf50dfac9659 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Tue, 5 Mar 2019 18:38:29 +0000 Subject: software X denoise demo --- demos/SoftwareX_supp/Demo_VolumeDenoise.py | 121 +++++++++++++++++++++++++++++ 1 file changed, 121 insertions(+) create mode 100644 demos/SoftwareX_supp/Demo_VolumeDenoise.py (limited to 'demos/SoftwareX_supp') diff --git a/demos/SoftwareX_supp/Demo_VolumeDenoise.py b/demos/SoftwareX_supp/Demo_VolumeDenoise.py new file mode 100644 index 0000000..2387e94 --- /dev/null +++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py @@ -0,0 +1,121 @@ +#!/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, NDF, Diff4th +#%% +print ("Building 3D phantom using TomoPhantom software") +tic=timeit.default_timer() +model = 9 # select a model number from the library +N_size = 256 # 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 ("#############ROF TV CPU####################") +# set parameters +pars = {'algorithm': ROF_TV, \ + 'input' : phantom_noise,\ + 'regularisation_parameter':0.04,\ + 'number_of_iterations': 100,\ + 'time_marching_parameter': 0.0025 + } + +tic=timeit.default_timer() +rof_cpu3D = ROF_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'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[128,:,:]*255, rof_cpu3D[128,:,:]*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.04,\ + 'number_of_iterations': 600,\ + 'time_marching_parameter': 0.0025 + } + +tic=timeit.default_timer() +rof_gpu3D = ROF_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'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[128,:,:]*255, rof_gpu3D[128,:,:]*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)) + +#%% + -- cgit v1.2.3 From 39baef90c4b209090f006e5308653cb0a3348c4e Mon Sep 17 00:00:00 2001 From: dkazanc Date: Wed, 6 Mar 2019 15:13:58 +0000 Subject: tol work --- demos/SoftwareX_supp/Demo_VolumeDenoise.py | 62 +++++++++++++++++++++++++++++- 1 file changed, 61 insertions(+), 1 deletion(-) (limited to 'demos/SoftwareX_supp') diff --git a/demos/SoftwareX_supp/Demo_VolumeDenoise.py b/demos/SoftwareX_supp/Demo_VolumeDenoise.py index 2387e94..17cdf4d 100644 --- a/demos/SoftwareX_supp/Demo_VolumeDenoise.py +++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py @@ -60,13 +60,14 @@ 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.04,\ - 'number_of_iterations': 100,\ + 'number_of_iterations': 600,\ 'time_marching_parameter': 0.0025 } @@ -116,6 +117,65 @@ 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.04,\ + 'number_of_iterations': 300,\ + 'time_marching_parameter': 0.0025,\ + 'tolerance_constant':1e-05,\ + } + +tic=timeit.default_timer() +fgp_cpu3D = FGP_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'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[128,:,:]*255, fgp_cpu3D[128,:,:]*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.04,\ + 'number_of_iterations': 300,\ + 'time_marching_parameter': 0.0025,\ + 'tolerance_constant':1e-05,\ + } +tic=timeit.default_timer() +fgp_gpu3D = FGP_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'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[128,:,:]*255, fgp_gpu3D[128,:,:]*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)) #%% + + -- cgit v1.2.3 From cfcc4be4413f65a0b9c4ef197687e3a167eff0e8 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 6 Mar 2019 23:34:55 +0000 Subject: cont1 --- demos/SoftwareX_supp/Demo_VolumeDenoise.py | 44 ++++++++++++++++-------------- 1 file changed, 24 insertions(+), 20 deletions(-) (limited to 'demos/SoftwareX_supp') diff --git a/demos/SoftwareX_supp/Demo_VolumeDenoise.py b/demos/SoftwareX_supp/Demo_VolumeDenoise.py index 17cdf4d..6e7ea46 100644 --- a/demos/SoftwareX_supp/Demo_VolumeDenoise.py +++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py @@ -120,19 +120,21 @@ print("ROF-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof, #%% print ("#############FGP TV CPU####################") # set parameters -pars = {'algorithm':FGP_TV, \ +pars = {'algorithm' : FGP_TV, \ 'input' : phantom_noise,\ - 'regularisation_parameter':0.04,\ - 'number_of_iterations': 300,\ - 'time_marching_parameter': 0.0025,\ - 'tolerance_constant':1e-05,\ - } + 'regularisation_parameter':0.05, \ + 'number_of_iterations' :100 ,\ + 'tolerance_constant':1e-04,\ + 'methodTV': 0 ,\ + 'nonneg': 0} tic=timeit.default_timer() -fgp_cpu3D = FGP_TV(pars['input'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['time_marching_parameter'],'cpu') +(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 @@ -149,19 +151,21 @@ print("FGP-TV (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof, #%% print ("#############FGP TV GPU####################") # set parameters -pars = {'algorithm':FGP_TV, \ +pars = {'algorithm' : FGP_TV, \ 'input' : phantom_noise,\ - 'regularisation_parameter':0.04,\ - 'number_of_iterations': 300,\ - 'time_marching_parameter': 0.0025,\ - 'tolerance_constant':1e-05,\ - } + 'regularisation_parameter':0.05, \ + 'number_of_iterations' :80 ,\ + 'tolerance_constant':1e-04,\ + 'methodTV': 0 ,\ + 'nonneg': 0} tic=timeit.default_timer() -fgp_gpu3D = FGP_TV(pars['input'], - pars['regularisation_parameter'], - pars['number_of_iterations'], - pars['time_marching_parameter'],'gpu') +(fgp_gpu3D) = 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 -- cgit v1.2.3 From e2208bfc2ed540065bef2e8e12d914d873464da7 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Sun, 10 Mar 2019 22:23:22 +0000 Subject: all python code updated --- demos/SoftwareX_supp/Demo_VolumeDenoise.py | 55 +++++++++++++++++++++--------- 1 file changed, 38 insertions(+), 17 deletions(-) (limited to 'demos/SoftwareX_supp') diff --git a/demos/SoftwareX_supp/Demo_VolumeDenoise.py b/demos/SoftwareX_supp/Demo_VolumeDenoise.py index 6e7ea46..07e3133 100644 --- a/demos/SoftwareX_supp/Demo_VolumeDenoise.py +++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py @@ -29,7 +29,7 @@ from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, LLT_ROF, NDF, Diff4 #%% print ("Building 3D phantom using TomoPhantom software") tic=timeit.default_timer() -model = 9 # select a model number from the library +model = 16 # select a model number from the library N_size = 256 # Define phantom dimensions using a scalar value (cubic phantom) path = os.path.dirname(tomophantom.__file__) path_library3D = os.path.join(path, "Phantom3DLibrary.dat") @@ -66,16 +66,18 @@ print ("#############ROF TV CPU####################") # set parameters pars = {'algorithm': ROF_TV, \ 'input' : phantom_noise,\ - 'regularisation_parameter':0.04,\ - 'number_of_iterations': 600,\ - 'time_marching_parameter': 0.0025 - } + 'regularisation_parameter':0.02,\ + 'number_of_iterations': 1000,\ + 'time_marching_parameter': 0.001,\ + 'tolerance_constant':0.0} tic=timeit.default_timer() -rof_cpu3D = ROF_TV(pars['input'], +(rof_cpu3D, infcpu) = ROF_TV(pars['input'], pars['regularisation_parameter'], pars['number_of_iterations'], - pars['time_marching_parameter'],'cpu') + pars['time_marching_parameter'], + pars['tolerance_constant'],'cpu') + toc=timeit.default_timer() Run_time_rof = toc - tic @@ -94,28 +96,47 @@ print ("#############ROF TV GPU####################") # set parameters pars = {'algorithm': ROF_TV, \ 'input' : phantom_noise,\ - 'regularisation_parameter':0.04,\ - 'number_of_iterations': 600,\ - 'time_marching_parameter': 0.0025 - } + 'regularisation_parameter':0.06,\ + 'number_of_iterations': 10000,\ + 'time_marching_parameter': 0.00025,\ + 'tolerance_constant':1e-06} tic=timeit.default_timer() -rof_gpu3D = ROF_TV(pars['input'], +(rof_gpu3D, infogpu) = ROF_TV(pars['input'], pars['regularisation_parameter'], pars['number_of_iterations'], - pars['time_marching_parameter'],'gpu') + 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() +sliceNo = 128 # SSIM measure -Qtools = QualityTools(phantom_tm[128,:,:]*255, rof_gpu3D[128,:,:]*235) +Qtools = QualityTools(phantom_tm[sliceNo,:,:]*255, rof_gpu3D[sliceNo,:,:]*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####################") @@ -154,13 +175,13 @@ print ("#############FGP TV GPU####################") pars = {'algorithm' : FGP_TV, \ 'input' : phantom_noise,\ 'regularisation_parameter':0.05, \ - 'number_of_iterations' :80 ,\ - 'tolerance_constant':1e-04,\ + 'number_of_iterations' :1500 ,\ + 'tolerance_constant':1e-06,\ 'methodTV': 0 ,\ 'nonneg': 0} tic=timeit.default_timer() -(fgp_gpu3D) = FGP_TV(pars['input'], +(fgp_gpu3D,infogpu) = FGP_TV(pars['input'], pars['regularisation_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], -- cgit v1.2.3 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/SoftwareX_supp') 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 From 420e71a0dcb42e91e1aa93306c2e2f688b309620 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Tue, 12 Mar 2019 17:29:07 +0000 Subject: cmakelists fixes, matlab wrappers done --- demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 4 ++-- demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 2 +- demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py | 3 +-- 3 files changed, 4 insertions(+), 5 deletions(-) (limited to 'demos/SoftwareX_supp') diff --git a/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py b/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py index 01491d9..ca8f1d2 100644 --- a/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py +++ b/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py @@ -111,7 +111,7 @@ RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # 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") @@ -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 99b9fe8..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 inner iterations earlier + tolerance = 0.0, # tolerance to stop inner (regularisation) iterations earlier device='gpu') #%% print ("Reconstructing with ADMM method using SB-TV penalty") -- cgit v1.2.3 From 1ac06b5ce11b247930489b7aa3afa59215e43c91 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 12 Mar 2019 22:14:27 +0000 Subject: readme updates and demos --- demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) (limited to 'demos/SoftwareX_supp') diff --git a/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py b/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py index ca8f1d2..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,7 +108,7 @@ 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 = 0.0, # tolerance to stop inner (regularisation) iterations earlier @@ -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') -- cgit v1.2.3