diff options
5 files changed, 643 insertions, 0 deletions
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py new file mode 100644 index 0000000..02e7c5c --- /dev/null +++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py @@ -0,0 +1,190 @@ +#!/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, Mark Basham, +Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019 +____________________________________________________________________________ +* Reads real tomographic data (stored at Zenodo) +* 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. TomoPhantom software for data generation + +@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk +GPLv3 license (ASTRA toolbox) +""" +import numpy as np +import matplotlib.pyplot as plt +import scipy.io +from tomorec.supp.suppTools import normaliser + +# load dendritic data +datadict = scipy.io.loadmat('Rawdata_1_frame3D.mat') +# extract data (print(datadict.keys())) +dataRaw = datadict['Rawdata3D'] +angles = datadict['AnglesAr'] +flats = datadict['flats3D'] +darks= datadict['darks3D'] + +dataRaw = np.swapaxes(dataRaw,0,1) +#%% +#flats2 = np.zeros((np.size(flats,0),1, np.size(flats,1)), dtype='float32') +#flats2[:,0,:] = flats[:] +#darks2 = np.zeros((np.size(darks,0),1, np.size(darks,1)), dtype='float32') +#darks2[:,0,:] = darks[:] + +# normalise the data, required format is [detectorsHoriz, Projections, Slices] +data_norm = normaliser(dataRaw, flats, darks, log='log') + +#dataRaw = np.float32(np.divide(dataRaw, np.max(dataRaw).astype(float))) + +intens_max = 70 +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,0) +N_size = 1000 +slice_to_recon = 0 # select which slice to reconstruct +angles_rad = angles*(np.pi/180.0) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +from tomorec.methodsDIR import RecToolsDIR +RectoolsDIR = RecToolsDIR(DetectorsDimH = detectorHoriz, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV = None, # 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(np.transpose(data_norm[:,:,slice_to_recon])) + +plt.figure() +plt.imshow(FBPrec[150:550,150:550], vmin=0, vmax=0.005, cmap="gray") +plt.title('FBP reconstruction') + +from tomorec.methodsIR import RecToolsIR +# set parameters and initiate a class object +Rectools = RecToolsIR(DetectorsDimH = detectorHoriz, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV = None, # 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='PWLS',# data fidelity, choose LS, PWLS, GH (wip), Student (wip) + nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE') + OS_number = 12, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets + tolerance = 1e-08, # tolerance to stop outer iterations earlier + device='gpu') + +lc = Rectools.powermethod(np.transpose(dataRaw[:,:,slice_to_recon])) # calculate Lipschitz constant (run once to initilise) + +RecFISTA_os_pwls = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \ + np.transpose(dataRaw[:,:,slice_to_recon]), \ + iterationsFISTA = 15, \ + lipschitz_const = lc) + +fig = plt.figure() +plt.imshow(RecFISTA_os_pwls[150:550,150:550], vmin=0, vmax=0.003, cmap="gray") +#plt.imshow(RecFISTA_os_pwls, vmin=0, vmax=0.004, cmap="gray") +plt.title('FISTA PWLS-OS reconstruction') +plt.show() +#fig.savefig('dendr_PWLS.png', dpi=200) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("Reconstructing with FISTA PWLS-OS-TV method %%%%%%%%%%%%%%%%") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +# Run FISTA-PWLS-OS reconstrucion algorithm with regularisation +RecFISTA_pwls_os_TV = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \ + np.transpose(dataRaw[:,:,slice_to_recon]), \ + iterationsFISTA = 15, \ + regularisation = 'FGP_TV', \ + regularisation_parameter = 0.000001,\ + regularisation_iterations = 200,\ + lipschitz_const = lc) + +fig = plt.figure() +plt.imshow(RecFISTA_pwls_os_TV[150:550,150:550], vmin=0, vmax=0.003, cmap="gray") +#plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') +plt.title('FISTA PWLS-OS-TV reconstruction') +plt.show() +#fig.savefig('dendr_TV.png', dpi=200) +#%% +""" +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("Reconstructing with FISTA PWLS-OS-Diff4th method %%%%%%%%%%%") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +# Run FISTA-PWLS-OS reconstrucion algorithm with regularisation +RecFISTA_pwls_os_Diff4th = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \ + np.transpose(dataRaw[:,:,slice_to_recon]), \ + iterationsFISTA = 15, \ + regularisation = 'DIFF4th', \ + regularisation_parameter = 0.1,\ + time_marching_parameter = 0.001,\ + edge_param = 0.003,\ + regularisation_iterations = 600,\ + lipschitz_const = lc) + +fig = plt.figure() +plt.imshow(RecFISTA_pwls_os_Diff4th[150:550,150:550], vmin=0, vmax=0.004, cmap="gray") +#plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') +plt.title('FISTA PWLS-OS-Diff4th reconstruction') +plt.show() +#fig.savefig('dendr_Diff4th.png', dpi=200) +""" +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("Reconstructing with FISTA PWLS-OS-ROF_LLT method %%%%%%%%%%%") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +# Run FISTA-PWLS-OS reconstrucion algorithm with regularisation +RecFISTA_pwls_os_rofllt = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \ + np.transpose(dataRaw[:,:,slice_to_recon]), \ + iterationsFISTA = 15, \ + regularisation = 'LLT_ROF', \ + regularisation_parameter = 0.000007,\ + regularisation_parameter2 = 0.0004,\ + regularisation_iterations = 350,\ + lipschitz_const = lc) + +fig = plt.figure() +plt.imshow(RecFISTA_pwls_os_rofllt[150:550,150:550], vmin=0, vmax=0.003, cmap="gray") +#plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') +plt.title('FISTA PWLS-OS-ROF-LLT reconstruction') +plt.show() +#fig.savefig('dendr_ROFLLT.png', dpi=200) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("Reconstructing with FISTA PWLS-OS-TGV method %%%%%%%%%%%%%%%") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +# Run FISTA-PWLS-OS reconstrucion algorithm with regularisation +RecFISTA_pwls_os_tgv = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \ + np.transpose(dataRaw[:,:,slice_to_recon]), \ + iterationsFISTA = 15, \ + regularisation = 'TGV', \ + regularisation_parameter = 0.001,\ + TGV_alpha2 = 0.001,\ + TGV_alpha1 = 2.0,\ + TGV_LipschitzConstant = 24,\ + regularisation_iterations = 300,\ + lipschitz_const = lc) + +fig = plt.figure() +plt.imshow(RecFISTA_pwls_os_tgv[150:550,150:550], vmin=0, vmax=0.003, cmap="gray") +#plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') +plt.title('FISTA PWLS-OS-TGV reconstruction') +plt.show() +#fig.savefig('dendr_TGV.png', dpi=200) diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py new file mode 100644 index 0000000..467e810 --- /dev/null +++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py @@ -0,0 +1,188 @@ +#!/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, Mark Basham, +Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019 +____________________________________________________________________________ +* Reads data which is previosly generated by TomoPhantom software (Zenodo link) +* Optimises for the regularisation parameters which later used in the script: +Demo_SimulData_Recon_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 +3. TomoPhantom software for data generation + +@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk +GPLv3 license (ASTRA toolbox) +""" +import timeit +import os +import matplotlib.pyplot as plt +import numpy as np +import tomophantom +from tomophantom import TomoP3D +from tomophantom.supp.qualitymetrics import QualityTools +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 = 9 # select a model number from the library +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) +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.25*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 = 70 +sliceSel = 150 +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 = 10000 # controls the level of noise +sigma_flats = 0.085 # 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 = 70 +sliceSel = 150 +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() +#%% +# 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 = angles_rad, # array of angles in radians + ObjSize = N_size, # a scalar to define reconstructed object dimensions + device = 'gpu') +#%% +print ("Reconstruction using FBP from TomoRec") +recNumerical= RectoolsDIR.FBP(projData3D_norm) # FBP reconstruction + +sliceSel = int(0.5*N_size) +max_val = 1 +#plt.gray() +plt.figure() +plt.subplot(131) +plt.imshow(recNumerical[sliceSel,:,:],vmin=0, vmax=max_val) +plt.title('3D Reconstruction, axial view') + +plt.subplot(132) +plt.imshow(recNumerical[:,sliceSel,:],vmin=0, vmax=max_val) +plt.title('3D Reconstruction, coronal view') + +plt.subplot(133) +plt.imshow(recNumerical[:,:,sliceSel],vmin=0, vmax=max_val) +plt.title('3D Reconstruction, sagittal view') +plt.show() + +# calculate errors +Qtools = QualityTools(phantom_tm, recNumerical) +RMSE_fbp = Qtools.rmse() +print("Root Mean Square Error for FBP is {}".format(RMSE_fbp)) +#%% +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 = 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') +#%% +# Run ADMM reconstrucion algorithm with 3D regularisation +RecADMM_reg_fgptv = RectoolsIR.ADMM(projData3D_norm, + rho_const = 2000.0, \ + iterationsADMM = 30, \ + regularisation = 'FGP_TV', \ + regularisation_parameter = 0.003,\ + regularisation_iterations = 250) + +sliceSel = int(0.5*N_size) +max_val = 1 +plt.figure() +plt.subplot(131) +plt.imshow(RecADMM_reg_fgptv[sliceSel,:,:],vmin=0, vmax=max_val) +plt.title('3D ADMM-FGP-TV Reconstruction, axial view') + +plt.subplot(132) +plt.imshow(RecADMM_reg_fgptv[:,sliceSel,:],vmin=0, vmax=max_val) +plt.title('3D ADMM-FGP-TV Reconstruction, coronal view') + +plt.subplot(133) +plt.imshow(RecADMM_reg_fgptv[:,:,sliceSel],vmin=0, vmax=max_val) +plt.title('3D ADMM-FGP-TV Reconstruction, sagittal view') +plt.show() + +# calculate errors +Qtools = QualityTools(phantom_tm, RecADMM_reg_fgptv) +RMSE_admm_fgp = Qtools.rmse() +print("Root Mean Square Error for ADMM-FGP-TV is {}".format(RMSE_admm_fgp)) + +#%% diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py new file mode 100644 index 0000000..0a99951 --- /dev/null +++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py @@ -0,0 +1,183 @@ +#!/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, Mark Basham, +Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019 +____________________________________________________________________________ +* Reads data which is previously generated by TomoPhantom software (Zenodo link) +* 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 +3. TomoPhantom software for data generation + +@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk +GPLv3 license (ASTRA toolbox) +""" +import timeit +import os +import matplotlib.pyplot as plt +import numpy as np +from tomophantom.supp.qualitymetrics import QualityTools + +print ("Building 3D phantom using TomoPhantom software") +tic=timeit.default_timer() +model = 9 # select a model number from the library +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) +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.25*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 = 70 +sliceSel = 150 +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 = 10000 # controls the level of noise +sigma_flats = 0.085 # 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 = 70 +sliceSel = 150 +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() +#%% +# 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 = angles_rad, # array of angles in radians + ObjSize = N_size, # a scalar to define reconstructed object dimensions + device = 'gpu') +#%% +print ("Reconstruction using FBP from TomoRec") +recNumerical= RectoolsDIR.FBP(projData3D_norm) # FBP reconstruction + +sliceSel = int(0.5*N_size) +max_val = 1 +#plt.gray() +plt.figure() +plt.subplot(131) +plt.imshow(recNumerical[sliceSel,:,:],vmin=0, vmax=max_val) +plt.title('3D Reconstruction, axial view') + +plt.subplot(132) +plt.imshow(recNumerical[:,sliceSel,:],vmin=0, vmax=max_val) +plt.title('3D Reconstruction, coronal view') + +plt.subplot(133) +plt.imshow(recNumerical[:,:,sliceSel],vmin=0, vmax=max_val) +plt.title('3D Reconstruction, sagittal view') +plt.show() + +# calculate errors +Qtools = QualityTools(phantom_tm, recNumerical) +RMSE_fbp = Qtools.rmse() +print("Root Mean Square Error for FBP is {}".format(RMSE_fbp)) +#%% +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 = 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') +#%% +# Run ADMM reconstrucion algorithm with 3D regularisation +RecADMM_reg_fgptv = RectoolsIR.ADMM(projData3D_norm, + rho_const = 2000.0, \ + iterationsADMM = 30, \ + regularisation = 'FGP_TV', \ + regularisation_parameter = 0.003,\ + regularisation_iterations = 250) + +sliceSel = int(0.5*N_size) +max_val = 1 +plt.figure() +plt.subplot(131) +plt.imshow(RecADMM_reg_fgptv[sliceSel,:,:],vmin=0, vmax=max_val) +plt.title('3D ADMM-FGP-TV Reconstruction, axial view') + +plt.subplot(132) +plt.imshow(RecADMM_reg_fgptv[:,sliceSel,:],vmin=0, vmax=max_val) +plt.title('3D ADMM-FGP-TV Reconstruction, coronal view') + +plt.subplot(133) +plt.imshow(RecADMM_reg_fgptv[:,:,sliceSel],vmin=0, vmax=max_val) +plt.title('3D ADMM-FGP-TV Reconstruction, sagittal view') +plt.show() + +# calculate errors +Qtools = QualityTools(phantom_tm, RecADMM_reg_fgptv) +RMSE_admm_fgp = Qtools.rmse() +print("Root Mean Square Error for ADMM-FGP-TV is {}".format(RMSE_admm_fgp)) + +#%% diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py new file mode 100644 index 0000000..3265d37 --- /dev/null +++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py @@ -0,0 +1,63 @@ +#!/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, Mark Basham, +Martin J. Turner, Philip J. Withers and Alun Ashton; 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 data generation + +@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk +Apache 2.0 license +""" +import numpy as np +import matplotlib.pyplot as plt +import scipy.io +from tomorec.supp.suppTools import normaliser + +# load dendritic data +datadict = scipy.io.loadmat('Rawdata_1_frame3D.mat') +# extract data (print(datadict.keys())) +dataRaw = datadict['Rawdata3D'] +angles = datadict['AnglesAr'] +flats = datadict['flats3D'] +darks= datadict['darks3D'] + +dataRaw = np.swapaxes(dataRaw,0,1) +#%% +#flats2 = np.zeros((np.size(flats,0),1, np.size(flats,1)), dtype='float32') +#flats2[:,0,:] = flats[:] +#darks2 = np.zeros((np.size(darks,0),1, np.size(darks,1)), dtype='float32') +#darks2[:,0,:] = darks[:] + +# normalise the data, required format is [detectorsHoriz, Projections, Slices] +data_norm = normaliser(dataRaw, flats, darks, log='log') + +#dataRaw = np.float32(np.divide(dataRaw, np.max(dataRaw).astype(float))) + +intens_max = 70 +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,0) +N_size = 1000 +slice_to_recon = 0 # select which slice to reconstruct +angles_rad = angles*(np.pi/180.0) +#%% diff --git a/Wrappers/Python/demos/SoftwareX_supp/Readme.md b/Wrappers/Python/demos/SoftwareX_supp/Readme.md new file mode 100644 index 0000000..fa77745 --- /dev/null +++ b/Wrappers/Python/demos/SoftwareX_supp/Readme.md @@ -0,0 +1,19 @@ +# SoftwareX publication [1] supporting files + +## Decription: +The scripts here support publication in SoftwareX journal [1] to ensure +reproducibility of the research. + +## Data: + +## Dependencies: + +## Files description: + + +### References: +[1] "CCPi-Regularisation Toolkit for computed tomographic image reconstruction with proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham, Martin J. Turner, Philip J. Withers and Alun Ashton + + +### 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 |