summaryrefslogtreecommitdiffstats
path: root/Wrappers/Python
diff options
context:
space:
mode:
authorTomas Kulhanek <tmkulhanek@gmail.com>2019-02-28 15:22:10 +0000
committerTomas Kulhanek <tmkulhanek@gmail.com>2019-02-28 15:22:10 +0000
commit68e6f3397e8a450854f39a5d514e1f747b9031a4 (patch)
treeeddf7bc14a998ffabc7e9e01f0cca2ac44b1d88a /Wrappers/Python
parent364a703de9f31b35d4301f3e913f519be9d3a14f (diff)
downloadregularization-68e6f3397e8a450854f39a5d514e1f747b9031a4.tar.gz
regularization-68e6f3397e8a450854f39a5d514e1f747b9031a4.tar.bz2
regularization-68e6f3397e8a450854f39a5d514e1f747b9031a4.tar.xz
regularization-68e6f3397e8a450854f39a5d514e1f747b9031a4.zip
merge
Diffstat (limited to 'Wrappers/Python')
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py231
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py161
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py309
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py117
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Readme.md26
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5bin2408 -> 0 bytes
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5bin2408 -> 0 bytes
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5bin2408 -> 0 bytes
8 files changed, 0 insertions, 844 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
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
--- a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5
+++ /dev/null
Binary files 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
--- a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5
+++ /dev/null
Binary files 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
--- a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5
+++ /dev/null
Binary files differ