diff options
| -rw-r--r-- | Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 48 | ||||
| -rw-r--r-- | Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 224 |
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 |