demos created

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