From 68e6f3397e8a450854f39a5d514e1f747b9031a4 Mon Sep 17 00:00:00 2001 From: Tomas Kulhanek Date: Thu, 28 Feb 2019 15:22:10 +0000 Subject: merge --- .../demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 231 --------------- .../SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 161 ----------- .../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 309 --------------------- .../demos/SoftwareX_supp/Demo_SimulData_SX.py | 117 -------- Wrappers/Python/demos/SoftwareX_supp/Readme.md | 26 -- .../optim_param/Optim_admm_rofllt.h5 | Bin 2408 -> 0 bytes .../SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 | Bin 2408 -> 0 bytes .../SoftwareX_supp/optim_param/Optim_admm_tgv.h5 | Bin 2408 -> 0 bytes demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 231 +++++++++++++++ .../SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 161 +++++++++++ demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py | 309 +++++++++++++++++++++ demos/SoftwareX_supp/Demo_SimulData_SX.py | 117 ++++++++ demos/SoftwareX_supp/Readme.md | 26 ++ .../optim_param/Optim_admm_rofllt.h5 | Bin 0 -> 2408 bytes .../SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 | Bin 0 -> 2408 bytes demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 | Bin 0 -> 2408 bytes recipe/meta.yaml | 2 + recipe/run_test.py | 2 + test/test_CPU_regularisers.py | 1 + 19 files changed, 849 insertions(+), 844 deletions(-) delete mode 100644 Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py delete mode 100644 Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py delete mode 100644 Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py delete mode 100644 Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py delete mode 100644 Wrappers/Python/demos/SoftwareX_supp/Readme.md delete mode 100644 Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 delete mode 100644 Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 delete mode 100644 Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 create mode 100644 demos/SoftwareX_supp/Demo_RealData_Recon_SX.py create mode 100644 demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py create mode 100644 demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py create mode 100644 demos/SoftwareX_supp/Demo_SimulData_SX.py create mode 100644 demos/SoftwareX_supp/Readme.md create mode 100644 demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 create mode 100644 demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 create mode 100644 demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py deleted file mode 100644 index 01491d9..0000000 --- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py +++ /dev/null @@ -1,231 +0,0 @@ -#!/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 -____________________________________________________________________________ -* Reads real tomographic data (stored at Zenodo) ---- https://doi.org/10.5281/zenodo.2578893 -* Reconstructs using TomoRec software -* Saves reconstructed images -____________________________________________________________________________ ->>>>> Dependencies: <<<<< -1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox -2. TomoRec: conda install -c dkazanc tomorec -or install from https://github.com/dkazanc/TomoRec -3. libtiff if one needs to save tiff images: - install pip install libtiff - -@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk -GPLv3 license (ASTRA toolbox) -""" -import numpy as np -import matplotlib.pyplot as plt -import h5py -from tomorec.supp.suppTools import normaliser -import time - -# load dendritic projection data -h5f = h5py.File('data/DendrData_3D.h5','r') -dataRaw = h5f['dataRaw'][:] -flats = h5f['flats'][:] -darks = h5f['darks'][:] -angles_rad = h5f['angles_rad'][:] -h5f.close() -#%% -# normalise the data [detectorsVert, Projections, detectorsHoriz] -data_norm = normaliser(dataRaw, flats, darks, log='log') -del dataRaw, darks, flats - -intens_max = 2.3 -plt.figure() -plt.subplot(131) -plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max) -plt.title('2D Projection (analytical)') -plt.subplot(132) -plt.imshow(data_norm[300,:,:],vmin=0, vmax=intens_max) -plt.title('Sinogram view') -plt.subplot(133) -plt.imshow(data_norm[:,:,600],vmin=0, vmax=intens_max) -plt.title('Tangentogram view') -plt.show() - -detectorHoriz = np.size(data_norm,2) -det_y_crop = [i for i in range(0,detectorHoriz-22)] -N_size = 950 # reconstruction domain -time_label = int(time.time()) -#%% -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%") -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -from tomorec.methodsDIR import RecToolsDIR - -RectoolsDIR = RecToolsDIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal) - 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 - device='gpu') - -FBPrec = RectoolsDIR.FBP(data_norm[0:100,:,det_y_crop]) - -sliceSel = 50 -max_val = 0.003 -plt.figure() -plt.subplot(131) -plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray") -plt.title('FBP Reconstruction, axial view') - -plt.subplot(132) -plt.imshow(FBPrec[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray") -plt.title('FBP Reconstruction, coronal view') - -plt.subplot(133) -plt.imshow(FBPrec[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray") -plt.title('FBP Reconstruction, sagittal view') -plt.show() - -# saving to tiffs (16bit) -""" -from libtiff import TIFF -FBPrec += np.abs(np.min(FBPrec)) -multiplier = (int)(65535/(np.max(FBPrec))) - -# saving to tiffs (16bit) -for i in range(0,np.size(FBPrec,0)): - tiff = TIFF.open('Dendr_FBP'+'_'+str(i)+'.tiff', mode='w') - tiff.write_image(np.uint16(FBPrec[i,:,:]*multiplier)) - tiff.close() -""" -#%% -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("Reconstructing with ADMM method using TomoRec software") -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -# initialise TomoRec ITERATIVE reconstruction class ONCE -from tomorec.methodsIR import RecToolsIR -RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal) - 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) - 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 - device='gpu') -#%% -print ("Reconstructing with ADMM method using SB-TV penalty") -RecADMM_reg_sbtv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop], - rho_const = 2000.0, \ - iterationsADMM = 15, \ - regularisation = 'SB_TV', \ - regularisation_parameter = 0.00085,\ - regularisation_iterations = 50) - -sliceSel = 50 -max_val = 0.003 -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') - -plt.subplot(132) -plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray") -plt.title('3D ADMM-SB-TV Reconstruction, coronal view') - -plt.subplot(133) -plt.imshow(RecADMM_reg_sbtv[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray") -plt.title('3D ADMM-SB-TV Reconstruction, sagittal view') -plt.show() - - -# saving to tiffs (16bit) -""" -from libtiff import TIFF -multiplier = (int)(65535/(np.max(RecADMM_reg_sbtv))) -for i in range(0,np.size(RecADMM_reg_sbtv,0)): - tiff = TIFF.open('Dendr_ADMM_SBTV'+'_'+str(i)+'.tiff', mode='w') - tiff.write_image(np.uint16(RecADMM_reg_sbtv[i,:,:]*multiplier)) - tiff.close() -""" -# Saving recpnstructed data with a unique time label -np.save('Dendr_ADMM_SBTV'+str(time_label)+'.npy', RecADMM_reg_sbtv) -del RecADMM_reg_sbtv -#%% -print ("Reconstructing with ADMM method using ROF-LLT penalty") -RecADMM_reg_rofllt = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop], - rho_const = 2000.0, \ - iterationsADMM = 15, \ - regularisation = 'LLT_ROF', \ - regularisation_parameter = 0.0009,\ - regularisation_parameter2 = 0.0007,\ - time_marching_parameter = 0.001,\ - regularisation_iterations = 550) - -sliceSel = 50 -max_val = 0.003 -plt.figure() -plt.subplot(131) -plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val) -plt.title('3D ADMM-ROFLLT Reconstruction, axial view') - -plt.subplot(132) -plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val) -plt.title('3D ADMM-ROFLLT Reconstruction, coronal view') - -plt.subplot(133) -plt.imshow(RecADMM_reg_rofllt[:,:,sliceSel],vmin=0, vmax=max_val) -plt.title('3D ADMM-ROFLLT Reconstruction, sagittal view') -plt.show() - -# saving to tiffs (16bit) -""" -from libtiff import TIFF -multiplier = (int)(65535/(np.max(RecADMM_reg_rofllt))) -for i in range(0,np.size(RecADMM_reg_rofllt,0)): - tiff = TIFF.open('Dendr_ADMM_ROFLLT'+'_'+str(i)+'.tiff', mode='w') - tiff.write_image(np.uint16(RecADMM_reg_rofllt[i,:,:]*multiplier)) - tiff.close() -""" - -# Saving recpnstructed data with a unique time label -np.save('Dendr_ADMM_ROFLLT'+str(time_label)+'.npy', RecADMM_reg_rofllt) -del RecADMM_reg_rofllt -#%% -print ("Reconstructing with ADMM method using TGV penalty") -RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop], - rho_const = 2000.0, \ - iterationsADMM = 15, \ - regularisation = 'TGV', \ - regularisation_parameter = 0.01,\ - regularisation_iterations = 500) - -sliceSel = 50 -max_val = 0.003 -plt.figure() -plt.subplot(131) -plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val) -plt.title('3D ADMM-TGV Reconstruction, axial view') - -plt.subplot(132) -plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val) -plt.title('3D ADMM-TGV Reconstruction, coronal view') - -plt.subplot(133) -plt.imshow(RecADMM_reg_tgv[:,:,sliceSel],vmin=0, vmax=max_val) -plt.title('3D ADMM-TGV Reconstruction, sagittal view') -plt.show() - -# saving to tiffs (16bit) -""" -from libtiff import TIFF -multiplier = (int)(65535/(np.max(RecADMM_reg_tgv))) -for i in range(0,np.size(RecADMM_reg_tgv,0)): - tiff = TIFF.open('Dendr_ADMM_TGV'+'_'+str(i)+'.tiff', mode='w') - tiff.write_image(np.uint16(RecADMM_reg_tgv[i,:,:]*multiplier)) - tiff.close() -""" -# 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/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py deleted file mode 100644 index 59ffc0e..0000000 --- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py +++ /dev/null @@ -1,161 +0,0 @@ -#!/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 -____________________________________________________________________________ -* Reads data which is previosly generated by TomoPhantom software (Zenodo link) ---- https://doi.org/10.5281/zenodo.2578893 -* Optimises for the regularisation parameters which later used in the script: -Demo_SimulData_Recon_SX.py -____________________________________________________________________________ ->>>>> Dependencies: <<<<< ->>>>> Dependencies: <<<<< -1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox -2. TomoRec: conda install -c dkazanc tomorec -or install from https://github.com/dkazanc/TomoRec - -@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk -GPLv3 license (ASTRA toolbox) -""" -#import timeit -import matplotlib.pyplot as plt -import numpy as np -import h5py -from ccpi.supp.qualitymetrics import QualityTools - -# loading the data -h5f = h5py.File('data/TomoSim_data1550671417.h5','r') -phantom = h5f['phantom'][:] -projdata_norm = h5f['projdata_norm'][:] -proj_angles = h5f['proj_angles'][:] -h5f.close() - -[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm) -N_size = Vert_det - -sliceSel = 128 -#plt.gray() -plt.figure() -plt.subplot(131) -plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1) -plt.title('3D Phantom, axial view') - -plt.subplot(132) -plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1) -plt.title('3D Phantom, coronal view') - -plt.subplot(133) -plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1) -plt.title('3D Phantom, sagittal view') -plt.show() - -intens_max = 240 -plt.figure() -plt.subplot(131) -plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max) -plt.title('2D Projection (erroneous)') -plt.subplot(132) -plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max) -plt.title('Sinogram view') -plt.subplot(133) -plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max) -plt.title('Tangentogram view') -plt.show() -#%% -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("Reconstructing with ADMM method using TomoRec software") -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -# initialise TomoRec ITERATIVE reconstruction class ONCE -from tomorec.methodsIR import RecToolsIR -RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal) - DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only - AnglesVec = proj_angles, # 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) - 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 - device='gpu') -#%% -param_space = 30 -reg_param_sb_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters -erros_vec_sbtv = np.zeros((param_space)) # a vector of errors - -print ("Reconstructing with ADMM method using SB-TV penalty") -for i in range(0,param_space): - RecADMM_reg_sbtv = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 15, \ - regularisation = 'SB_TV', \ - regularisation_parameter = reg_param_sb_vec[i],\ - regularisation_iterations = 50) - # calculate errors - Qtools = QualityTools(phantom, RecADMM_reg_sbtv) - erros_vec_sbtv[i] = Qtools.rmse() - print("RMSE for regularisation parameter {} for ADMM-SB-TV is {}".format(reg_param_sb_vec[i],erros_vec_sbtv[i])) - -plt.figure() -plt.plot(erros_vec_sbtv) - -# Saving generated data with a unique time label -h5f = h5py.File('Optim_admm_sbtv.h5', 'w') -h5f.create_dataset('reg_param_sb_vec', data=reg_param_sb_vec) -h5f.create_dataset('erros_vec_sbtv', data=erros_vec_sbtv) -h5f.close() -#%% -param_space = 30 -reg_param_rofllt_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters -erros_vec_rofllt = np.zeros((param_space)) # a vector of errors - -print ("Reconstructing with ADMM method using ROF-LLT penalty") -for i in range(0,param_space): - RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 15, \ - regularisation = 'LLT_ROF', \ - regularisation_parameter = reg_param_rofllt_vec[i],\ - regularisation_parameter2 = 0.005,\ - regularisation_iterations = 600) - # calculate errors - Qtools = QualityTools(phantom, RecADMM_reg_rofllt) - erros_vec_rofllt[i] = Qtools.rmse() - print("RMSE for regularisation parameter {} for ADMM-ROF-LLT is {}".format(reg_param_rofllt_vec[i],erros_vec_rofllt[i])) - -plt.figure() -plt.plot(erros_vec_rofllt) - -# Saving generated data with a unique time label -h5f = h5py.File('Optim_admm_rofllt.h5', 'w') -h5f.create_dataset('reg_param_rofllt_vec', data=reg_param_rofllt_vec) -h5f.create_dataset('erros_vec_rofllt', data=erros_vec_rofllt) -h5f.close() -#%% -param_space = 30 -reg_param_tgv_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters -erros_vec_tgv = np.zeros((param_space)) # a vector of errors - -print ("Reconstructing with ADMM method using TGV penalty") -for i in range(0,param_space): - RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 15, \ - regularisation = 'TGV', \ - regularisation_parameter = reg_param_tgv_vec[i],\ - regularisation_iterations = 600) - # calculate errors - Qtools = QualityTools(phantom, RecADMM_reg_tgv) - erros_vec_tgv[i] = Qtools.rmse() - print("RMSE for regularisation parameter {} for ADMM-TGV is {}".format(reg_param_tgv_vec[i],erros_vec_tgv[i])) - -plt.figure() -plt.plot(erros_vec_tgv) - -# Saving generated data with a unique time label -h5f = h5py.File('Optim_admm_tgv.h5', 'w') -h5f.create_dataset('reg_param_tgv_vec', data=reg_param_tgv_vec) -h5f.create_dataset('erros_vec_tgv', data=erros_vec_tgv) -h5f.close() -#%% \ No newline at end of file diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py deleted file mode 100644 index 93b0cef..0000000 --- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py +++ /dev/null @@ -1,309 +0,0 @@ -#!/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 -____________________________________________________________________________ -* Reads data which is previously generated by TomoPhantom software (Zenodo link) ---- https://doi.org/10.5281/zenodo.2578893 -* Reconstruct using optimised regularisation parameters (see Demo_SimulData_ParOptimis_SX.py) -____________________________________________________________________________ ->>>>> Dependencies: <<<<< -1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox -2. TomoRec: conda install -c dkazanc tomorec -or install from https://github.com/dkazanc/TomoRec - -@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk -GPLv3 license (ASTRA toolbox) -""" -#import timeit -import matplotlib.pyplot as plt -import matplotlib.gridspec as gridspec -import numpy as np -import h5py -from ccpi.supp.qualitymetrics import QualityTools -from scipy.signal import gaussian - -# loading the data -h5f = h5py.File('data/TomoSim_data1550671417.h5','r') -phantom = h5f['phantom'][:] -projdata_norm = h5f['projdata_norm'][:] -proj_angles = h5f['proj_angles'][:] -h5f.close() - -[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm) -N_size = Vert_det - -# loading optmisation parameters (the result of running Demo_SimulData_ParOptimis_SX) -h5f = h5py.File('optim_param/Optim_admm_sbtv.h5','r') -reg_param_sb_vec = h5f['reg_param_sb_vec'][:] -erros_vec_sbtv = h5f['erros_vec_sbtv'][:] -h5f.close() - -h5f = h5py.File('optim_param/Optim_admm_rofllt.h5','r') -reg_param_rofllt_vec = h5f['reg_param_rofllt_vec'][:] -erros_vec_rofllt = h5f['erros_vec_rofllt'][:] -h5f.close() - -h5f = h5py.File('optim_param/Optim_admm_tgv.h5','r') -reg_param_tgv_vec = h5f['reg_param_tgv_vec'][:] -erros_vec_tgv = h5f['erros_vec_tgv'][:] -h5f.close() - -index_minSBTV = min(xrange(len(erros_vec_sbtv)), key=erros_vec_sbtv.__getitem__) -index_minROFLLT = min(xrange(len(erros_vec_rofllt)), key=erros_vec_rofllt.__getitem__) -index_minTGV = min(xrange(len(erros_vec_tgv)), key=erros_vec_tgv.__getitem__) -# assign optimal regularisation parameters: -optimReg_sbtv = reg_param_sb_vec[index_minSBTV] -optimReg_rofllt = reg_param_rofllt_vec[index_minROFLLT] -optimReg_tgv = reg_param_tgv_vec[index_minTGV] -#%% -# plot loaded data -sliceSel = 128 -#plt.figure() -fig, (ax1, ax2) = plt.subplots(figsize=(15, 5), ncols=2) -plt.rcParams.update({'xtick.labelsize': 'x-small'}) -plt.rcParams.update({'ytick.labelsize':'x-small'}) -plt.subplot(121) -one = plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1, interpolation='none', cmap="PuOr") -fig.colorbar(one, ax=ax1) -plt.title('3D Phantom, axial (X-Y) view') -plt.subplot(122) -two = plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1,interpolation='none', cmap="PuOr") -fig.colorbar(two, ax=ax2) -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() -#%% -intens_max = 220 -plt.figure() -plt.rcParams.update({'xtick.labelsize': 'x-small'}) -plt.rcParams.update({'ytick.labelsize':'x-small'}) -plt.subplot(131) -plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max, cmap="PuOr") -plt.xlabel('X-detector', fontsize=16) -plt.ylabel('Z-detector', fontsize=16) -plt.title('2D Projection (X-Z) view', fontsize=19) -plt.subplot(132) -plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max, cmap="PuOr") -plt.xlabel('X-detector', fontsize=16) -plt.ylabel('Projection angle', fontsize=16) -plt.title('Sinogram (X-Y) view', fontsize=19) -plt.subplot(133) -plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max, cmap="PuOr") -plt.xlabel('Projection angle', fontsize=16) -plt.ylabel('Z-detector', fontsize=16) -plt.title('Vertical (Y-Z) view', fontsize=19) -plt.show() -#plt.savefig('projdata.pdf', format='pdf', dpi=1200) -#%% -# initialise TomoRec DIRECT reconstruction class ONCE -from tomorec.methodsDIR import RecToolsDIR -RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal) - DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only - AnglesVec = proj_angles, # array of angles in radians - ObjSize = N_size, # a scalar to define reconstructed object dimensions - device = 'gpu') -#%% -print ("Reconstruction using FBP from TomoRec") -recFBP= RectoolsDIR.FBP(projdata_norm) # FBP reconstruction -#%% -x0, y0 = 0, 127 # These are in _pixel_ coordinates!! -x1, y1 = 255, 127 - -sliceSel = int(0.5*N_size) -max_val = 1 -plt.figure(figsize = (20,5)) -gs1 = gridspec.GridSpec(1, 3) -gs1.update(wspace=0.1, hspace=0.05) # set the spacing between axes. -ax1 = plt.subplot(gs1[0]) -plt.imshow(recFBP[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") -ax1.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') -plt.colorbar(ax=ax1) -plt.title('FBP Reconstruction, axial (X-Y) view', fontsize=19) -ax1.set_aspect('equal') -ax3 = plt.subplot(gs1[1]) -plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2) -plt.plot(recFBP[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') -plt.title('Profile', fontsize=19) -ax2 = plt.subplot(gs1[2]) -plt.imshow(recFBP[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") -plt.title('FBP Reconstruction, coronal (Y-Z) view', fontsize=19) -ax2.set_aspect('equal') -plt.show() -#plt.savefig('FBP_phantom.pdf', format='pdf', dpi=1600) - -# calculate errors -Qtools = QualityTools(phantom, recFBP) -RMSE_fbp = Qtools.rmse() -print("Root Mean Square Error for FBP is {}".format(RMSE_fbp)) - -# SSIM measure -Qtools = QualityTools(phantom[128,:,:]*255, recFBP[128,:,:]*235) -win = np.array([gaussian(11, 1.5)]) -win2d = win * (win.T) -ssim_fbp = Qtools.ssim(win2d) -print("Mean SSIM for FBP is {}".format(ssim_fbp[0])) -#%% -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("Reconstructing with ADMM method using TomoRec software") -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -# initialise TomoRec ITERATIVE reconstruction class ONCE -from tomorec.methodsIR import RecToolsIR -RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal) - DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only - AnglesVec = proj_angles, # 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) - 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 - device='gpu') -#%% -print ("Reconstructing with ADMM method using SB-TV penalty") -RecADMM_reg_sbtv = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 25, \ - regularisation = 'SB_TV', \ - regularisation_parameter = optimReg_sbtv,\ - regularisation_iterations = 50) - -sliceSel = int(0.5*N_size) -max_val = 1 -plt.figure(figsize = (20,3)) -gs1 = gridspec.GridSpec(1, 4) -gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. -ax1 = plt.subplot(gs1[0]) -plt.plot(reg_param_sb_vec, erros_vec_sbtv, color='k',linewidth=2) -plt.xlabel('Regularisation parameter', fontsize=16) -plt.ylabel('RMSE value', fontsize=16) -plt.title('Regularisation selection', fontsize=19) -ax2 = plt.subplot(gs1[1]) -plt.imshow(RecADMM_reg_sbtv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") -ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') -plt.title('ADMM-SBTV (X-Y) view', fontsize=19) -#ax2.set_aspect('equal') -ax3 = plt.subplot(gs1[2]) -plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2) -plt.plot(RecADMM_reg_sbtv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') -plt.title('Profile', fontsize=19) -ax4 = plt.subplot(gs1[3]) -plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") -plt.title('ADMM-SBTV (Y-Z) view', fontsize=19) -plt.colorbar(ax=ax4) -plt.show() -plt.savefig('SBTV_phantom.pdf', format='pdf', dpi=1600) - -# calculate errors -Qtools = QualityTools(phantom, RecADMM_reg_sbtv) -RMSE_admm_sbtv = Qtools.rmse() -print("Root Mean Square Error for ADMM-SB-TV is {}".format(RMSE_admm_sbtv)) - -# SSIM measure -Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_sbtv[128,:,:]*235) -win = np.array([gaussian(11, 1.5)]) -win2d = win * (win.T) -ssim_admm_sbtv = Qtools.ssim(win2d) -print("Mean SSIM ADMM-SBTV is {}".format(ssim_admm_sbtv[0])) -#%% -print ("Reconstructing with ADMM method using ROFLLT penalty") -RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 25, \ - regularisation = 'LLT_ROF', \ - regularisation_parameter = optimReg_rofllt,\ - regularisation_parameter2 = 0.0085,\ - regularisation_iterations = 600) - -sliceSel = int(0.5*N_size) -max_val = 1 -plt.figure(figsize = (20,3)) -gs1 = gridspec.GridSpec(1, 4) -gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. -ax1 = plt.subplot(gs1[0]) -plt.plot(reg_param_rofllt_vec, erros_vec_rofllt, color='k',linewidth=2) -plt.xlabel('Regularisation parameter', fontsize=16) -plt.ylabel('RMSE value', fontsize=16) -plt.title('Regularisation selection', fontsize=19) -ax2 = plt.subplot(gs1[1]) -plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") -ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') -plt.title('ADMM-ROFLLT (X-Y) view', fontsize=19) -#ax2.set_aspect('equal') -ax3 = plt.subplot(gs1[2]) -plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2) -plt.plot(RecADMM_reg_rofllt[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') -plt.title('Profile', fontsize=19) -ax4 = plt.subplot(gs1[3]) -plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") -plt.title('ADMM-ROFLLT (Y-Z) view', fontsize=19) -plt.colorbar(ax=ax4) -plt.show() -#plt.savefig('ROFLLT_phantom.pdf', format='pdf', dpi=1600) - -# calculate errors -Qtools = QualityTools(phantom, RecADMM_reg_rofllt) -RMSE_admm_rofllt = Qtools.rmse() -print("Root Mean Square Error for ADMM-ROF-LLT is {}".format(RMSE_admm_rofllt)) - -# SSIM measure -Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_rofllt[128,:,:]*235) -win = np.array([gaussian(11, 1.5)]) -win2d = win * (win.T) -ssim_admm_rifllt = Qtools.ssim(win2d) -print("Mean SSIM ADMM-ROFLLT is {}".format(ssim_admm_rifllt[0])) -#%% -print ("Reconstructing with ADMM method using TGV penalty") -RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 25, \ - regularisation = 'TGV', \ - regularisation_parameter = optimReg_tgv,\ - regularisation_iterations = 600) -#%% -sliceSel = int(0.5*N_size) -max_val = 1 -plt.figure(figsize = (20,3)) -gs1 = gridspec.GridSpec(1, 4) -gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. -ax1 = plt.subplot(gs1[0]) -plt.plot(reg_param_tgv_vec, erros_vec_tgv, color='k',linewidth=2) -plt.xlabel('Regularisation parameter', fontsize=16) -plt.ylabel('RMSE value', fontsize=16) -plt.title('Regularisation selection', fontsize=19) -ax2 = plt.subplot(gs1[1]) -plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") -ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') -plt.title('ADMM-TGV (X-Y) view', fontsize=19) -#ax2.set_aspect('equal') -ax3 = plt.subplot(gs1[2]) -plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2) -plt.plot(RecADMM_reg_tgv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') -plt.title('Profile', fontsize=19) -ax4 = plt.subplot(gs1[3]) -plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") -plt.title('ADMM-TGV (Y-Z) view', fontsize=19) -plt.colorbar(ax=ax4) -plt.show() -#plt.savefig('TGV_phantom.pdf', format='pdf', dpi=1600) - -# calculate errors -Qtools = QualityTools(phantom, RecADMM_reg_tgv) -RMSE_admm_tgv = Qtools.rmse() -print("Root Mean Square Error for ADMM-TGV is {}".format(RMSE_admm_tgv)) - -# SSIM measure -#Create a 2d gaussian for the window parameter -Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_tgv[128,:,:]*235) -win = np.array([gaussian(11, 1.5)]) -win2d = win * (win.T) -ssim_admm_tgv = Qtools.ssim(win2d) -print("Mean SSIM ADMM-TGV is {}".format(ssim_admm_tgv[0])) -#%% \ No newline at end of file diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py deleted file mode 100644 index cdf4325..0000000 --- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py +++ /dev/null @@ -1,117 +0,0 @@ -#!/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 -____________________________________________________________________________ -* Runs TomoPhantom software to simulate tomographic projection data with -some imaging errors and noise -* Saves the data into hdf file to be uploaded in reconstruction scripts -__________________________________________________________________________ - ->>>>> Dependencies: <<<<< -1. TomoPhantom software for phantom and data generation - -@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk -Apache 2.0 license -""" -import timeit -import os -import matplotlib.pyplot as plt -import numpy as np -import tomophantom -from tomophantom import TomoP3D -from tomophantom.supp.flatsgen import flats -from tomophantom.supp.normraw import normaliser_sim - -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) -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)) - -sliceSel = int(0.5*N_size) -#plt.gray() -plt.figure() -plt.subplot(131) -plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1) -plt.title('3D Phantom, axial view') - -plt.subplot(132) -plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1) -plt.title('3D Phantom, coronal view') - -plt.subplot(133) -plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1) -plt.title('3D Phantom, sagittal view') -plt.show() - -# Projection geometry related parameters: -Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal) -Vert_det = N_size # detector row count (vertical) (no reason for it to be > N) -angles_num = int(0.35*np.pi*N_size); # angles number -angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees -angles_rad = angles*(np.pi/180.0) -#%% -print ("Building 3D analytical projection data with TomoPhantom") -projData3D_analyt= TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D) - -intens_max = N_size -sliceSel = int(0.5*N_size) -plt.figure() -plt.subplot(131) -plt.imshow(projData3D_analyt[:,sliceSel,:],vmin=0, vmax=intens_max) -plt.title('2D Projection (analytical)') -plt.subplot(132) -plt.imshow(projData3D_analyt[sliceSel,:,:],vmin=0, vmax=intens_max) -plt.title('Sinogram view') -plt.subplot(133) -plt.imshow(projData3D_analyt[:,:,sliceSel],vmin=0, vmax=intens_max) -plt.title('Tangentogram view') -plt.show() -#%% -print ("Simulate flat fields, add noise and normalise projections...") -flatsnum = 20 # generate 20 flat fields -flatsSIM = flats(Vert_det, Horiz_det, maxheight = 0.1, maxthickness = 3, sigma_noise = 0.2, sigmasmooth = 3, flatsnum=flatsnum) - -plt.figure() -plt.imshow(flatsSIM[0,:,:],vmin=0, vmax=1) -plt.title('A selected simulated flat-field') -#%% -# Apply normalisation of data and add noise -flux_intensity = 60000 # controls the level of noise -sigma_flats = 0.01 # contro the level of noise in flats (higher creates more ring artifacts) -projData3D_norm = normaliser_sim(projData3D_analyt, flatsSIM, sigma_flats, flux_intensity) - -intens_max = N_size -sliceSel = int(0.5*N_size) -plt.figure() -plt.subplot(131) -plt.imshow(projData3D_norm[:,sliceSel,:],vmin=0, vmax=intens_max) -plt.title('2D Projection (erroneous)') -plt.subplot(132) -plt.imshow(projData3D_norm[sliceSel,:,:],vmin=0, vmax=intens_max) -plt.title('Sinogram view') -plt.subplot(133) -plt.imshow(projData3D_norm[:,:,sliceSel],vmin=0, vmax=intens_max) -plt.title('Tangentogram view') -plt.show() -#%% -import h5py -import time -time_label = int(time.time()) -# Saving generated data with a unique time label -h5f = h5py.File('TomoSim_data'+str(time_label)+'.h5', 'w') -h5f.create_dataset('phantom', data=phantom_tm) -h5f.create_dataset('projdata_norm', data=projData3D_norm) -h5f.create_dataset('proj_angles', data=angles_rad) -h5f.close() -#%% \ No newline at end of file diff --git a/Wrappers/Python/demos/SoftwareX_supp/Readme.md b/Wrappers/Python/demos/SoftwareX_supp/Readme.md deleted file mode 100644 index 54e83f1..0000000 --- a/Wrappers/Python/demos/SoftwareX_supp/Readme.md +++ /dev/null @@ -1,26 +0,0 @@ - -# SoftwareX publication [1] supporting files - -## Decription: -The scripts here support publication in SoftwareX journal [1] to ensure reproducibility of the research. The scripts linked with data shared at Zenodo. - -## Data: -Data is shared at Zenodo [here](https://doi.org/10.5281/zenodo.2578893) - -## Dependencies: -1. [ASTRA toolbox](https://github.com/astra-toolbox/astra-toolbox): `conda install -c astra-toolbox astra-toolbox` -2. [TomoRec](https://github.com/dkazanc/TomoRec): `conda install -c dkazanc tomorec` -3. [Tomophantom](https://github.com/dkazanc/TomoPhantom): `conda install tomophantom -c ccpi` - -## Files description: -- `Demo_SimulData_SX.py` - simulates 3D projection data using [Tomophantom](https://github.com/dkazanc/TomoPhantom) software. One can skip this module if the data is taken from [Zenodo](https://doi.org/10.5281/zenodo.2578893) -- `Demo_SimulData_ParOptimis_SX.py` - runs computationally extensive calculations for optimal regularisation parameters, the result are saved into directory `optim_param`. This script can be also skipped. -- `Demo_SimulData_Recon_SX.py` - using established regularisation parameters, one runs iterative reconstruction -- `Demo_RealData_Recon_SX.py` - runs real data reconstructions. Can be quite intense on memory so reduce the size of the reconstructed volume if needed. - -### References: -[1] "CCPi-Regularisation Toolkit for computed tomographic image reconstruction with proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner and Philip J. Withers; SoftwareX, 2019. - -### Acknowledgments: -CCPi-RGL software is a product of the [CCPi](https://www.ccpi.ac.uk/) group, STFC SCD software developers and Diamond Light Source (DLS). Any relevant questions/comments can be e-mailed to Daniil Kazantsev at dkazanc@hotmail.com - diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 deleted file mode 100644 index 63bc4fd..0000000 Binary files a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 and /dev/null differ diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 deleted file mode 100644 index 03c0c14..0000000 Binary files a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 and /dev/null differ diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 deleted file mode 100644 index 056d915..0000000 Binary files a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 and /dev/null differ diff --git a/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py b/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py new file mode 100644 index 0000000..01491d9 --- /dev/null +++ b/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py @@ -0,0 +1,231 @@ +#!/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 +____________________________________________________________________________ +* Reads real tomographic data (stored at Zenodo) +--- https://doi.org/10.5281/zenodo.2578893 +* Reconstructs using TomoRec software +* Saves reconstructed images +____________________________________________________________________________ +>>>>> Dependencies: <<<<< +1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox +2. TomoRec: conda install -c dkazanc tomorec +or install from https://github.com/dkazanc/TomoRec +3. libtiff if one needs to save tiff images: + install pip install libtiff + +@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk +GPLv3 license (ASTRA toolbox) +""" +import numpy as np +import matplotlib.pyplot as plt +import h5py +from tomorec.supp.suppTools import normaliser +import time + +# load dendritic projection data +h5f = h5py.File('data/DendrData_3D.h5','r') +dataRaw = h5f['dataRaw'][:] +flats = h5f['flats'][:] +darks = h5f['darks'][:] +angles_rad = h5f['angles_rad'][:] +h5f.close() +#%% +# normalise the data [detectorsVert, Projections, detectorsHoriz] +data_norm = normaliser(dataRaw, flats, darks, log='log') +del dataRaw, darks, flats + +intens_max = 2.3 +plt.figure() +plt.subplot(131) +plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max) +plt.title('2D Projection (analytical)') +plt.subplot(132) +plt.imshow(data_norm[300,:,:],vmin=0, vmax=intens_max) +plt.title('Sinogram view') +plt.subplot(133) +plt.imshow(data_norm[:,:,600],vmin=0, vmax=intens_max) +plt.title('Tangentogram view') +plt.show() + +detectorHoriz = np.size(data_norm,2) +det_y_crop = [i for i in range(0,detectorHoriz-22)] +N_size = 950 # reconstruction domain +time_label = int(time.time()) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +from tomorec.methodsDIR import RecToolsDIR + +RectoolsDIR = RecToolsDIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal) + 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 + device='gpu') + +FBPrec = RectoolsDIR.FBP(data_norm[0:100,:,det_y_crop]) + +sliceSel = 50 +max_val = 0.003 +plt.figure() +plt.subplot(131) +plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray") +plt.title('FBP Reconstruction, axial view') + +plt.subplot(132) +plt.imshow(FBPrec[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray") +plt.title('FBP Reconstruction, coronal view') + +plt.subplot(133) +plt.imshow(FBPrec[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray") +plt.title('FBP Reconstruction, sagittal view') +plt.show() + +# saving to tiffs (16bit) +""" +from libtiff import TIFF +FBPrec += np.abs(np.min(FBPrec)) +multiplier = (int)(65535/(np.max(FBPrec))) + +# saving to tiffs (16bit) +for i in range(0,np.size(FBPrec,0)): + tiff = TIFF.open('Dendr_FBP'+'_'+str(i)+'.tiff', mode='w') + tiff.write_image(np.uint16(FBPrec[i,:,:]*multiplier)) + tiff.close() +""" +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("Reconstructing with ADMM method using TomoRec software") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +# initialise TomoRec ITERATIVE reconstruction class ONCE +from tomorec.methodsIR import RecToolsIR +RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal) + 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) + 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 + device='gpu') +#%% +print ("Reconstructing with ADMM method using SB-TV penalty") +RecADMM_reg_sbtv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop], + rho_const = 2000.0, \ + iterationsADMM = 15, \ + regularisation = 'SB_TV', \ + regularisation_parameter = 0.00085,\ + regularisation_iterations = 50) + +sliceSel = 50 +max_val = 0.003 +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') + +plt.subplot(132) +plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray") +plt.title('3D ADMM-SB-TV Reconstruction, coronal view') + +plt.subplot(133) +plt.imshow(RecADMM_reg_sbtv[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray") +plt.title('3D ADMM-SB-TV Reconstruction, sagittal view') +plt.show() + + +# saving to tiffs (16bit) +""" +from libtiff import TIFF +multiplier = (int)(65535/(np.max(RecADMM_reg_sbtv))) +for i in range(0,np.size(RecADMM_reg_sbtv,0)): + tiff = TIFF.open('Dendr_ADMM_SBTV'+'_'+str(i)+'.tiff', mode='w') + tiff.write_image(np.uint16(RecADMM_reg_sbtv[i,:,:]*multiplier)) + tiff.close() +""" +# Saving recpnstructed data with a unique time label +np.save('Dendr_ADMM_SBTV'+str(time_label)+'.npy', RecADMM_reg_sbtv) +del RecADMM_reg_sbtv +#%% +print ("Reconstructing with ADMM method using ROF-LLT penalty") +RecADMM_reg_rofllt = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop], + rho_const = 2000.0, \ + iterationsADMM = 15, \ + regularisation = 'LLT_ROF', \ + regularisation_parameter = 0.0009,\ + regularisation_parameter2 = 0.0007,\ + time_marching_parameter = 0.001,\ + regularisation_iterations = 550) + +sliceSel = 50 +max_val = 0.003 +plt.figure() +plt.subplot(131) +plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val) +plt.title('3D ADMM-ROFLLT Reconstruction, axial view') + +plt.subplot(132) +plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val) +plt.title('3D ADMM-ROFLLT Reconstruction, coronal view') + +plt.subplot(133) +plt.imshow(RecADMM_reg_rofllt[:,:,sliceSel],vmin=0, vmax=max_val) +plt.title('3D ADMM-ROFLLT Reconstruction, sagittal view') +plt.show() + +# saving to tiffs (16bit) +""" +from libtiff import TIFF +multiplier = (int)(65535/(np.max(RecADMM_reg_rofllt))) +for i in range(0,np.size(RecADMM_reg_rofllt,0)): + tiff = TIFF.open('Dendr_ADMM_ROFLLT'+'_'+str(i)+'.tiff', mode='w') + tiff.write_image(np.uint16(RecADMM_reg_rofllt[i,:,:]*multiplier)) + tiff.close() +""" + +# Saving recpnstructed data with a unique time label +np.save('Dendr_ADMM_ROFLLT'+str(time_label)+'.npy', RecADMM_reg_rofllt) +del RecADMM_reg_rofllt +#%% +print ("Reconstructing with ADMM method using TGV penalty") +RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop], + rho_const = 2000.0, \ + iterationsADMM = 15, \ + regularisation = 'TGV', \ + regularisation_parameter = 0.01,\ + regularisation_iterations = 500) + +sliceSel = 50 +max_val = 0.003 +plt.figure() +plt.subplot(131) +plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val) +plt.title('3D ADMM-TGV Reconstruction, axial view') + +plt.subplot(132) +plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val) +plt.title('3D ADMM-TGV Reconstruction, coronal view') + +plt.subplot(133) +plt.imshow(RecADMM_reg_tgv[:,:,sliceSel],vmin=0, vmax=max_val) +plt.title('3D ADMM-TGV Reconstruction, sagittal view') +plt.show() + +# saving to tiffs (16bit) +""" +from libtiff import TIFF +multiplier = (int)(65535/(np.max(RecADMM_reg_tgv))) +for i in range(0,np.size(RecADMM_reg_tgv,0)): + tiff = TIFF.open('Dendr_ADMM_TGV'+'_'+str(i)+'.tiff', mode='w') + tiff.write_image(np.uint16(RecADMM_reg_tgv[i,:,:]*multiplier)) + tiff.close() +""" +# 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 new file mode 100644 index 0000000..59ffc0e --- /dev/null +++ b/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py @@ -0,0 +1,161 @@ +#!/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 +____________________________________________________________________________ +* Reads data which is previosly generated by TomoPhantom software (Zenodo link) +--- https://doi.org/10.5281/zenodo.2578893 +* Optimises for the regularisation parameters which later used in the script: +Demo_SimulData_Recon_SX.py +____________________________________________________________________________ +>>>>> Dependencies: <<<<< +>>>>> Dependencies: <<<<< +1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox +2. TomoRec: conda install -c dkazanc tomorec +or install from https://github.com/dkazanc/TomoRec + +@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk +GPLv3 license (ASTRA toolbox) +""" +#import timeit +import matplotlib.pyplot as plt +import numpy as np +import h5py +from ccpi.supp.qualitymetrics import QualityTools + +# loading the data +h5f = h5py.File('data/TomoSim_data1550671417.h5','r') +phantom = h5f['phantom'][:] +projdata_norm = h5f['projdata_norm'][:] +proj_angles = h5f['proj_angles'][:] +h5f.close() + +[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm) +N_size = Vert_det + +sliceSel = 128 +#plt.gray() +plt.figure() +plt.subplot(131) +plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1) +plt.title('3D Phantom, axial view') + +plt.subplot(132) +plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1) +plt.title('3D Phantom, coronal view') + +plt.subplot(133) +plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1) +plt.title('3D Phantom, sagittal view') +plt.show() + +intens_max = 240 +plt.figure() +plt.subplot(131) +plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max) +plt.title('2D Projection (erroneous)') +plt.subplot(132) +plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max) +plt.title('Sinogram view') +plt.subplot(133) +plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max) +plt.title('Tangentogram view') +plt.show() +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("Reconstructing with ADMM method using TomoRec software") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +# initialise TomoRec ITERATIVE reconstruction class ONCE +from tomorec.methodsIR import RecToolsIR +RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only + AnglesVec = proj_angles, # 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) + 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 + device='gpu') +#%% +param_space = 30 +reg_param_sb_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters +erros_vec_sbtv = np.zeros((param_space)) # a vector of errors + +print ("Reconstructing with ADMM method using SB-TV penalty") +for i in range(0,param_space): + RecADMM_reg_sbtv = RectoolsIR.ADMM(projdata_norm, + rho_const = 2000.0, \ + iterationsADMM = 15, \ + regularisation = 'SB_TV', \ + regularisation_parameter = reg_param_sb_vec[i],\ + regularisation_iterations = 50) + # calculate errors + Qtools = QualityTools(phantom, RecADMM_reg_sbtv) + erros_vec_sbtv[i] = Qtools.rmse() + print("RMSE for regularisation parameter {} for ADMM-SB-TV is {}".format(reg_param_sb_vec[i],erros_vec_sbtv[i])) + +plt.figure() +plt.plot(erros_vec_sbtv) + +# Saving generated data with a unique time label +h5f = h5py.File('Optim_admm_sbtv.h5', 'w') +h5f.create_dataset('reg_param_sb_vec', data=reg_param_sb_vec) +h5f.create_dataset('erros_vec_sbtv', data=erros_vec_sbtv) +h5f.close() +#%% +param_space = 30 +reg_param_rofllt_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters +erros_vec_rofllt = np.zeros((param_space)) # a vector of errors + +print ("Reconstructing with ADMM method using ROF-LLT penalty") +for i in range(0,param_space): + RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm, + rho_const = 2000.0, \ + iterationsADMM = 15, \ + regularisation = 'LLT_ROF', \ + regularisation_parameter = reg_param_rofllt_vec[i],\ + regularisation_parameter2 = 0.005,\ + regularisation_iterations = 600) + # calculate errors + Qtools = QualityTools(phantom, RecADMM_reg_rofllt) + erros_vec_rofllt[i] = Qtools.rmse() + print("RMSE for regularisation parameter {} for ADMM-ROF-LLT is {}".format(reg_param_rofllt_vec[i],erros_vec_rofllt[i])) + +plt.figure() +plt.plot(erros_vec_rofllt) + +# Saving generated data with a unique time label +h5f = h5py.File('Optim_admm_rofllt.h5', 'w') +h5f.create_dataset('reg_param_rofllt_vec', data=reg_param_rofllt_vec) +h5f.create_dataset('erros_vec_rofllt', data=erros_vec_rofllt) +h5f.close() +#%% +param_space = 30 +reg_param_tgv_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters +erros_vec_tgv = np.zeros((param_space)) # a vector of errors + +print ("Reconstructing with ADMM method using TGV penalty") +for i in range(0,param_space): + RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm, + rho_const = 2000.0, \ + iterationsADMM = 15, \ + regularisation = 'TGV', \ + regularisation_parameter = reg_param_tgv_vec[i],\ + regularisation_iterations = 600) + # calculate errors + Qtools = QualityTools(phantom, RecADMM_reg_tgv) + erros_vec_tgv[i] = Qtools.rmse() + print("RMSE for regularisation parameter {} for ADMM-TGV is {}".format(reg_param_tgv_vec[i],erros_vec_tgv[i])) + +plt.figure() +plt.plot(erros_vec_tgv) + +# Saving generated data with a unique time label +h5f = h5py.File('Optim_admm_tgv.h5', 'w') +h5f.create_dataset('reg_param_tgv_vec', data=reg_param_tgv_vec) +h5f.create_dataset('erros_vec_tgv', data=erros_vec_tgv) +h5f.close() +#%% \ No newline at end of file diff --git a/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py b/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py new file mode 100644 index 0000000..93b0cef --- /dev/null +++ b/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py @@ -0,0 +1,309 @@ +#!/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 +____________________________________________________________________________ +* Reads data which is previously generated by TomoPhantom software (Zenodo link) +--- https://doi.org/10.5281/zenodo.2578893 +* Reconstruct using optimised regularisation parameters (see Demo_SimulData_ParOptimis_SX.py) +____________________________________________________________________________ +>>>>> Dependencies: <<<<< +1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox +2. TomoRec: conda install -c dkazanc tomorec +or install from https://github.com/dkazanc/TomoRec + +@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk +GPLv3 license (ASTRA toolbox) +""" +#import timeit +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import numpy as np +import h5py +from ccpi.supp.qualitymetrics import QualityTools +from scipy.signal import gaussian + +# loading the data +h5f = h5py.File('data/TomoSim_data1550671417.h5','r') +phantom = h5f['phantom'][:] +projdata_norm = h5f['projdata_norm'][:] +proj_angles = h5f['proj_angles'][:] +h5f.close() + +[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm) +N_size = Vert_det + +# loading optmisation parameters (the result of running Demo_SimulData_ParOptimis_SX) +h5f = h5py.File('optim_param/Optim_admm_sbtv.h5','r') +reg_param_sb_vec = h5f['reg_param_sb_vec'][:] +erros_vec_sbtv = h5f['erros_vec_sbtv'][:] +h5f.close() + +h5f = h5py.File('optim_param/Optim_admm_rofllt.h5','r') +reg_param_rofllt_vec = h5f['reg_param_rofllt_vec'][:] +erros_vec_rofllt = h5f['erros_vec_rofllt'][:] +h5f.close() + +h5f = h5py.File('optim_param/Optim_admm_tgv.h5','r') +reg_param_tgv_vec = h5f['reg_param_tgv_vec'][:] +erros_vec_tgv = h5f['erros_vec_tgv'][:] +h5f.close() + +index_minSBTV = min(xrange(len(erros_vec_sbtv)), key=erros_vec_sbtv.__getitem__) +index_minROFLLT = min(xrange(len(erros_vec_rofllt)), key=erros_vec_rofllt.__getitem__) +index_minTGV = min(xrange(len(erros_vec_tgv)), key=erros_vec_tgv.__getitem__) +# assign optimal regularisation parameters: +optimReg_sbtv = reg_param_sb_vec[index_minSBTV] +optimReg_rofllt = reg_param_rofllt_vec[index_minROFLLT] +optimReg_tgv = reg_param_tgv_vec[index_minTGV] +#%% +# plot loaded data +sliceSel = 128 +#plt.figure() +fig, (ax1, ax2) = plt.subplots(figsize=(15, 5), ncols=2) +plt.rcParams.update({'xtick.labelsize': 'x-small'}) +plt.rcParams.update({'ytick.labelsize':'x-small'}) +plt.subplot(121) +one = plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1, interpolation='none', cmap="PuOr") +fig.colorbar(one, ax=ax1) +plt.title('3D Phantom, axial (X-Y) view') +plt.subplot(122) +two = plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1,interpolation='none', cmap="PuOr") +fig.colorbar(two, ax=ax2) +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() +#%% +intens_max = 220 +plt.figure() +plt.rcParams.update({'xtick.labelsize': 'x-small'}) +plt.rcParams.update({'ytick.labelsize':'x-small'}) +plt.subplot(131) +plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max, cmap="PuOr") +plt.xlabel('X-detector', fontsize=16) +plt.ylabel('Z-detector', fontsize=16) +plt.title('2D Projection (X-Z) view', fontsize=19) +plt.subplot(132) +plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max, cmap="PuOr") +plt.xlabel('X-detector', fontsize=16) +plt.ylabel('Projection angle', fontsize=16) +plt.title('Sinogram (X-Y) view', fontsize=19) +plt.subplot(133) +plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max, cmap="PuOr") +plt.xlabel('Projection angle', fontsize=16) +plt.ylabel('Z-detector', fontsize=16) +plt.title('Vertical (Y-Z) view', fontsize=19) +plt.show() +#plt.savefig('projdata.pdf', format='pdf', dpi=1200) +#%% +# initialise TomoRec DIRECT reconstruction class ONCE +from tomorec.methodsDIR import RecToolsDIR +RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only + AnglesVec = proj_angles, # array of angles in radians + ObjSize = N_size, # a scalar to define reconstructed object dimensions + device = 'gpu') +#%% +print ("Reconstruction using FBP from TomoRec") +recFBP= RectoolsDIR.FBP(projdata_norm) # FBP reconstruction +#%% +x0, y0 = 0, 127 # These are in _pixel_ coordinates!! +x1, y1 = 255, 127 + +sliceSel = int(0.5*N_size) +max_val = 1 +plt.figure(figsize = (20,5)) +gs1 = gridspec.GridSpec(1, 3) +gs1.update(wspace=0.1, hspace=0.05) # set the spacing between axes. +ax1 = plt.subplot(gs1[0]) +plt.imshow(recFBP[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") +ax1.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') +plt.colorbar(ax=ax1) +plt.title('FBP Reconstruction, axial (X-Y) view', fontsize=19) +ax1.set_aspect('equal') +ax3 = plt.subplot(gs1[1]) +plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2) +plt.plot(recFBP[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') +plt.title('Profile', fontsize=19) +ax2 = plt.subplot(gs1[2]) +plt.imshow(recFBP[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") +plt.title('FBP Reconstruction, coronal (Y-Z) view', fontsize=19) +ax2.set_aspect('equal') +plt.show() +#plt.savefig('FBP_phantom.pdf', format='pdf', dpi=1600) + +# calculate errors +Qtools = QualityTools(phantom, recFBP) +RMSE_fbp = Qtools.rmse() +print("Root Mean Square Error for FBP is {}".format(RMSE_fbp)) + +# SSIM measure +Qtools = QualityTools(phantom[128,:,:]*255, recFBP[128,:,:]*235) +win = np.array([gaussian(11, 1.5)]) +win2d = win * (win.T) +ssim_fbp = Qtools.ssim(win2d) +print("Mean SSIM for FBP is {}".format(ssim_fbp[0])) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("Reconstructing with ADMM method using TomoRec software") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +# initialise TomoRec ITERATIVE reconstruction class ONCE +from tomorec.methodsIR import RecToolsIR +RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only + AnglesVec = proj_angles, # 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) + 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 + device='gpu') +#%% +print ("Reconstructing with ADMM method using SB-TV penalty") +RecADMM_reg_sbtv = RectoolsIR.ADMM(projdata_norm, + rho_const = 2000.0, \ + iterationsADMM = 25, \ + regularisation = 'SB_TV', \ + regularisation_parameter = optimReg_sbtv,\ + regularisation_iterations = 50) + +sliceSel = int(0.5*N_size) +max_val = 1 +plt.figure(figsize = (20,3)) +gs1 = gridspec.GridSpec(1, 4) +gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. +ax1 = plt.subplot(gs1[0]) +plt.plot(reg_param_sb_vec, erros_vec_sbtv, color='k',linewidth=2) +plt.xlabel('Regularisation parameter', fontsize=16) +plt.ylabel('RMSE value', fontsize=16) +plt.title('Regularisation selection', fontsize=19) +ax2 = plt.subplot(gs1[1]) +plt.imshow(RecADMM_reg_sbtv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") +ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') +plt.title('ADMM-SBTV (X-Y) view', fontsize=19) +#ax2.set_aspect('equal') +ax3 = plt.subplot(gs1[2]) +plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2) +plt.plot(RecADMM_reg_sbtv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') +plt.title('Profile', fontsize=19) +ax4 = plt.subplot(gs1[3]) +plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") +plt.title('ADMM-SBTV (Y-Z) view', fontsize=19) +plt.colorbar(ax=ax4) +plt.show() +plt.savefig('SBTV_phantom.pdf', format='pdf', dpi=1600) + +# calculate errors +Qtools = QualityTools(phantom, RecADMM_reg_sbtv) +RMSE_admm_sbtv = Qtools.rmse() +print("Root Mean Square Error for ADMM-SB-TV is {}".format(RMSE_admm_sbtv)) + +# SSIM measure +Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_sbtv[128,:,:]*235) +win = np.array([gaussian(11, 1.5)]) +win2d = win * (win.T) +ssim_admm_sbtv = Qtools.ssim(win2d) +print("Mean SSIM ADMM-SBTV is {}".format(ssim_admm_sbtv[0])) +#%% +print ("Reconstructing with ADMM method using ROFLLT penalty") +RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm, + rho_const = 2000.0, \ + iterationsADMM = 25, \ + regularisation = 'LLT_ROF', \ + regularisation_parameter = optimReg_rofllt,\ + regularisation_parameter2 = 0.0085,\ + regularisation_iterations = 600) + +sliceSel = int(0.5*N_size) +max_val = 1 +plt.figure(figsize = (20,3)) +gs1 = gridspec.GridSpec(1, 4) +gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. +ax1 = plt.subplot(gs1[0]) +plt.plot(reg_param_rofllt_vec, erros_vec_rofllt, color='k',linewidth=2) +plt.xlabel('Regularisation parameter', fontsize=16) +plt.ylabel('RMSE value', fontsize=16) +plt.title('Regularisation selection', fontsize=19) +ax2 = plt.subplot(gs1[1]) +plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") +ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') +plt.title('ADMM-ROFLLT (X-Y) view', fontsize=19) +#ax2.set_aspect('equal') +ax3 = plt.subplot(gs1[2]) +plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2) +plt.plot(RecADMM_reg_rofllt[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') +plt.title('Profile', fontsize=19) +ax4 = plt.subplot(gs1[3]) +plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") +plt.title('ADMM-ROFLLT (Y-Z) view', fontsize=19) +plt.colorbar(ax=ax4) +plt.show() +#plt.savefig('ROFLLT_phantom.pdf', format='pdf', dpi=1600) + +# calculate errors +Qtools = QualityTools(phantom, RecADMM_reg_rofllt) +RMSE_admm_rofllt = Qtools.rmse() +print("Root Mean Square Error for ADMM-ROF-LLT is {}".format(RMSE_admm_rofllt)) + +# SSIM measure +Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_rofllt[128,:,:]*235) +win = np.array([gaussian(11, 1.5)]) +win2d = win * (win.T) +ssim_admm_rifllt = Qtools.ssim(win2d) +print("Mean SSIM ADMM-ROFLLT is {}".format(ssim_admm_rifllt[0])) +#%% +print ("Reconstructing with ADMM method using TGV penalty") +RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm, + rho_const = 2000.0, \ + iterationsADMM = 25, \ + regularisation = 'TGV', \ + regularisation_parameter = optimReg_tgv,\ + regularisation_iterations = 600) +#%% +sliceSel = int(0.5*N_size) +max_val = 1 +plt.figure(figsize = (20,3)) +gs1 = gridspec.GridSpec(1, 4) +gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. +ax1 = plt.subplot(gs1[0]) +plt.plot(reg_param_tgv_vec, erros_vec_tgv, color='k',linewidth=2) +plt.xlabel('Regularisation parameter', fontsize=16) +plt.ylabel('RMSE value', fontsize=16) +plt.title('Regularisation selection', fontsize=19) +ax2 = plt.subplot(gs1[1]) +plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") +ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') +plt.title('ADMM-TGV (X-Y) view', fontsize=19) +#ax2.set_aspect('equal') +ax3 = plt.subplot(gs1[2]) +plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2) +plt.plot(RecADMM_reg_tgv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') +plt.title('Profile', fontsize=19) +ax4 = plt.subplot(gs1[3]) +plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") +plt.title('ADMM-TGV (Y-Z) view', fontsize=19) +plt.colorbar(ax=ax4) +plt.show() +#plt.savefig('TGV_phantom.pdf', format='pdf', dpi=1600) + +# calculate errors +Qtools = QualityTools(phantom, RecADMM_reg_tgv) +RMSE_admm_tgv = Qtools.rmse() +print("Root Mean Square Error for ADMM-TGV is {}".format(RMSE_admm_tgv)) + +# SSIM measure +#Create a 2d gaussian for the window parameter +Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_tgv[128,:,:]*235) +win = np.array([gaussian(11, 1.5)]) +win2d = win * (win.T) +ssim_admm_tgv = Qtools.ssim(win2d) +print("Mean SSIM ADMM-TGV is {}".format(ssim_admm_tgv[0])) +#%% \ No newline at end of file diff --git a/demos/SoftwareX_supp/Demo_SimulData_SX.py b/demos/SoftwareX_supp/Demo_SimulData_SX.py new file mode 100644 index 0000000..cdf4325 --- /dev/null +++ b/demos/SoftwareX_supp/Demo_SimulData_SX.py @@ -0,0 +1,117 @@ +#!/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 +____________________________________________________________________________ +* Runs TomoPhantom software to simulate tomographic projection data with +some imaging errors and noise +* Saves the data into hdf file to be uploaded in reconstruction scripts +__________________________________________________________________________ + +>>>>> Dependencies: <<<<< +1. TomoPhantom software for phantom and data generation + +@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk +Apache 2.0 license +""" +import timeit +import os +import matplotlib.pyplot as plt +import numpy as np +import tomophantom +from tomophantom import TomoP3D +from tomophantom.supp.flatsgen import flats +from tomophantom.supp.normraw import normaliser_sim + +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) +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)) + +sliceSel = int(0.5*N_size) +#plt.gray() +plt.figure() +plt.subplot(131) +plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1) +plt.title('3D Phantom, axial view') + +plt.subplot(132) +plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1) +plt.title('3D Phantom, coronal view') + +plt.subplot(133) +plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1) +plt.title('3D Phantom, sagittal view') +plt.show() + +# Projection geometry related parameters: +Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal) +Vert_det = N_size # detector row count (vertical) (no reason for it to be > N) +angles_num = int(0.35*np.pi*N_size); # angles number +angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees +angles_rad = angles*(np.pi/180.0) +#%% +print ("Building 3D analytical projection data with TomoPhantom") +projData3D_analyt= TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D) + +intens_max = N_size +sliceSel = int(0.5*N_size) +plt.figure() +plt.subplot(131) +plt.imshow(projData3D_analyt[:,sliceSel,:],vmin=0, vmax=intens_max) +plt.title('2D Projection (analytical)') +plt.subplot(132) +plt.imshow(projData3D_analyt[sliceSel,:,:],vmin=0, vmax=intens_max) +plt.title('Sinogram view') +plt.subplot(133) +plt.imshow(projData3D_analyt[:,:,sliceSel],vmin=0, vmax=intens_max) +plt.title('Tangentogram view') +plt.show() +#%% +print ("Simulate flat fields, add noise and normalise projections...") +flatsnum = 20 # generate 20 flat fields +flatsSIM = flats(Vert_det, Horiz_det, maxheight = 0.1, maxthickness = 3, sigma_noise = 0.2, sigmasmooth = 3, flatsnum=flatsnum) + +plt.figure() +plt.imshow(flatsSIM[0,:,:],vmin=0, vmax=1) +plt.title('A selected simulated flat-field') +#%% +# Apply normalisation of data and add noise +flux_intensity = 60000 # controls the level of noise +sigma_flats = 0.01 # contro the level of noise in flats (higher creates more ring artifacts) +projData3D_norm = normaliser_sim(projData3D_analyt, flatsSIM, sigma_flats, flux_intensity) + +intens_max = N_size +sliceSel = int(0.5*N_size) +plt.figure() +plt.subplot(131) +plt.imshow(projData3D_norm[:,sliceSel,:],vmin=0, vmax=intens_max) +plt.title('2D Projection (erroneous)') +plt.subplot(132) +plt.imshow(projData3D_norm[sliceSel,:,:],vmin=0, vmax=intens_max) +plt.title('Sinogram view') +plt.subplot(133) +plt.imshow(projData3D_norm[:,:,sliceSel],vmin=0, vmax=intens_max) +plt.title('Tangentogram view') +plt.show() +#%% +import h5py +import time +time_label = int(time.time()) +# Saving generated data with a unique time label +h5f = h5py.File('TomoSim_data'+str(time_label)+'.h5', 'w') +h5f.create_dataset('phantom', data=phantom_tm) +h5f.create_dataset('projdata_norm', data=projData3D_norm) +h5f.create_dataset('proj_angles', data=angles_rad) +h5f.close() +#%% \ No newline at end of file diff --git a/demos/SoftwareX_supp/Readme.md b/demos/SoftwareX_supp/Readme.md new file mode 100644 index 0000000..54e83f1 --- /dev/null +++ b/demos/SoftwareX_supp/Readme.md @@ -0,0 +1,26 @@ + +# SoftwareX publication [1] supporting files + +## Decription: +The scripts here support publication in SoftwareX journal [1] to ensure reproducibility of the research. The scripts linked with data shared at Zenodo. + +## Data: +Data is shared at Zenodo [here](https://doi.org/10.5281/zenodo.2578893) + +## Dependencies: +1. [ASTRA toolbox](https://github.com/astra-toolbox/astra-toolbox): `conda install -c astra-toolbox astra-toolbox` +2. [TomoRec](https://github.com/dkazanc/TomoRec): `conda install -c dkazanc tomorec` +3. [Tomophantom](https://github.com/dkazanc/TomoPhantom): `conda install tomophantom -c ccpi` + +## Files description: +- `Demo_SimulData_SX.py` - simulates 3D projection data using [Tomophantom](https://github.com/dkazanc/TomoPhantom) software. One can skip this module if the data is taken from [Zenodo](https://doi.org/10.5281/zenodo.2578893) +- `Demo_SimulData_ParOptimis_SX.py` - runs computationally extensive calculations for optimal regularisation parameters, the result are saved into directory `optim_param`. This script can be also skipped. +- `Demo_SimulData_Recon_SX.py` - using established regularisation parameters, one runs iterative reconstruction +- `Demo_RealData_Recon_SX.py` - runs real data reconstructions. Can be quite intense on memory so reduce the size of the reconstructed volume if needed. + +### References: +[1] "CCPi-Regularisation Toolkit for computed tomographic image reconstruction with proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner and Philip J. Withers; SoftwareX, 2019. + +### Acknowledgments: +CCPi-RGL software is a product of the [CCPi](https://www.ccpi.ac.uk/) group, STFC SCD software developers and Diamond Light Source (DLS). Any relevant questions/comments can be e-mailed to Daniil Kazantsev at dkazanc@hotmail.com + diff --git a/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 b/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 new file mode 100644 index 0000000..63bc4fd Binary files /dev/null and b/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 differ diff --git a/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 b/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 new file mode 100644 index 0000000..03c0c14 Binary files /dev/null and b/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 differ diff --git a/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 b/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 new file mode 100644 index 0000000..056d915 Binary files /dev/null and b/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 differ diff --git a/recipe/meta.yaml b/recipe/meta.yaml index 527ad32..6f36906 100644 --- a/recipe/meta.yaml +++ b/recipe/meta.yaml @@ -14,6 +14,8 @@ test: requires: - pillow - pillow=4.1.1 # [win] +# command: +# - unittest -d discover .... ../test requirements: build: diff --git a/recipe/run_test.py b/recipe/run_test.py index 21f3216..f551616 100755 --- a/recipe/run_test.py +++ b/recipe/run_test.py @@ -815,5 +815,7 @@ class TestRegularisers(unittest.TestCase): self.assertLess(abs(rms_fgp-rms_fgp_exp) , tolerance) + + if __name__ == '__main__': unittest.main() diff --git a/test/test_CPU_regularisers.py b/test/test_CPU_regularisers.py index 6af4cd4..379b989 100644 --- a/test/test_CPU_regularisers.py +++ b/test/test_CPU_regularisers.py @@ -2,6 +2,7 @@ import unittest import math import os import timeit +import numpy as np from ccpi.filters.regularisers import FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th, ROF_TV from testroutines import * -- cgit v1.2.3