demos_updated2

2 files changed, 120 insertions, 152 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
index 02e7c5c..9472a43 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
@@ -14,37 +14,29 @@ ____________________________________________________________________________
 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
+import h5py
 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)
+# 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()
 #%%
-#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')
+del dataRaw, darks, flats
 
-#dataRaw = np.float32(np.divide(dataRaw, np.max(dataRaw).astype(float)))
-
-intens_max = 70
+intens_max = 2
 plt.figure() 
 plt.subplot(131)
 plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max)
@@ -61,40 +53,41 @@ 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)
+det_y_crop = [i for i in range(0,detectorHoriz-22)]
+
+RectoolsDIR = RecToolsDIR(DetectorsDimH = np.size(det_y_crop),  # 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]))
+FBPrec = RectoolsDIR.FBP(np.transpose(data_norm[det_y_crop,:,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)
+Rectools = RecToolsIR(DetectorsDimH = np.size(det_y_crop),  # 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)
+                    datafidelity='LS',# 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)
+lc = Rectools.powermethod() # 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]), \
+RecFISTA_os_pwls = Rectools.FISTA(np.transpose(data_norm[det_y_crop,:,slice_to_recon]), \
                              iterationsFISTA = 15, \
                              lipschitz_const = lc)
 
@@ -109,8 +102,7 @@ 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]), \
+RecFISTA_pwls_os_TV = Rectools.FISTA(np.transpose(data_norm[det_y_crop,:,slice_to_recon]), \
                               iterationsFISTA = 15, \
                               regularisation = 'FGP_TV', \
                               regularisation_parameter = 0.000001,\
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
index 467e810..148fdcc 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
@@ -11,136 +11,59 @@ ____________________________________________________________________________
 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
-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 timeit
 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
+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()
 
-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))
+[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm)
+N_size = Vert_det
 
-sliceSel = int(0.5*N_size)
+sliceSel = 128
 #plt.gray()
 plt.figure() 
 plt.subplot(131)
-plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1)
+plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1)
 plt.title('3D Phantom, axial view')
 
 plt.subplot(132)
-plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1)
+plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1)
 plt.title('3D Phantom, coronal view')
 
 plt.subplot(133)
-plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1)
+plt.imshow(phantom[:,:,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
+intens_max = 240
 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.imshow(projdata_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.imshow(projdata_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.imshow(projdata_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 ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -148,7 +71,7 @@ print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 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
+                    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')
@@ -156,33 +79,86 @@ RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector d
                     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)
+param_space = 20
+reg_param_sb_vec = np.linspace(0.001,0.03,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]))
 
-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()
+plt.plot(erros_vec_sbtv)
+#%%
+param_space = 20
+reg_param_rofllt_vec = np.linspace(0.001,0.03,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]))
 
-# 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))
+plt.figure() 
+plt.plot(erros_vec_rofllt)
+#%%
+param_space = 20
+reg_param_tgv_vec = np.linspace(0.001,0.03,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)
 #%%
+param_space = 1
+reg_param_diff4th = np.linspace(10,100,param_space,dtype='float32') # a vector of parameters
+erros_vec_diff4th = np.zeros((param_space)) # a vector of errors
+
+print ("Reconstructing with ADMM method using Diff4th penalty")
+for i in range(0,param_space):
+    RecADMM_reg_diff4th = RectoolsIR.ADMM(projdata_norm,
+                                  rho_const = 2000.0, \
+                                  iterationsADMM = 15, \
+                                  regularisation = 'Diff4th', \
+                                  regularisation_parameter = reg_param_diff4th[i],\
+                                  edge_param = 0.03,
+                                  time_marching_parameter = 0.001,\
+                                  regularisation_iterations = 750)
+    # calculate errors 
+    Qtools = QualityTools(phantom, RecADMM_reg_diff4th)
+    erros_vec_diff4th[i] = Qtools.rmse()
+    print("RMSE for regularisation parameter {} for ADMM-Diff4th is {}".format(reg_param_diff4th[i],erros_vec_diff4th[i]))
+
+plt.figure() 
+plt.plot(erros_vec_diff4th)
+#%%
\ No newline at end of file