software X denoise demo

1 files changed, 121 insertions, 0 deletions

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))
+
+#%%
+