diff options
| -rw-r--r-- | Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 2 | ||||
| -rw-r--r-- | Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py | 259 |
2 files changed, 168 insertions, 93 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 4cd680e..dc4eff2 100644 --- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py +++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py @@ -181,7 +181,7 @@ RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # 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='cpu') + device='gpu') print ("Reconstructing with ADMM method using TGV penalty") RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:10,:,det_y_crop], diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py index a022ad7..cf9e797 100644 --- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py +++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py @@ -19,6 +19,7 @@ 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 @@ -33,34 +34,73 @@ 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') +# 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() -plt.subplot(132) -plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1) -plt.title('3D Phantom, coronal view') +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) +plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1, cmap="PuOr") plt.title('3D Phantom, sagittal view') -plt.show() -intens_max = 240 +""" +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) -plt.title('2D Projection (erroneous)') +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) -plt.title('Sinogram view') +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) -plt.title('Tangentogram view') +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 @@ -72,23 +112,31 @@ RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector #%% 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.gray() -plt.figure() -plt.subplot(131) -plt.imshow(recFBP[sliceSel,:,:],vmin=0, vmax=max_val) -plt.title('3D Reconstruction, axial view') - -plt.subplot(132) -plt.imshow(recFBP[:,sliceSel,:],vmin=0, vmax=max_val) -plt.title('3D Reconstruction, coronal view') - -plt.subplot(133) -plt.imshow(recFBP[:,:,sliceSel],vmin=0, vmax=max_val) -plt.title('3D Reconstruction, sagittal view') +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=1200) # calculate errors Qtools = QualityTools(phantom, recFBP) @@ -110,87 +158,114 @@ RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector d tolerance = 1e-08, # tolerance to stop outer iterations earlier device='gpu') #%% -print ("Reconstructing with ADMM method using ROF-TV penalty") -RecADMM_reg_roftv = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 30, \ - regularisation = 'ROF_TV', \ - regularisation_parameter = 0.003,\ - regularisation_iterations = 500) +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() -plt.subplot(131) -plt.imshow(RecADMM_reg_roftv[sliceSel,:,:],vmin=0, vmax=max_val) -plt.title('3D ADMM-ROF-TV Reconstruction, axial view') - -plt.subplot(132) -plt.imshow(RecADMM_reg_roftv[:,sliceSel,:],vmin=0, vmax=max_val) -plt.title('3D ADMM-ROF-TV Reconstruction, coronal view') - -plt.subplot(133) -plt.imshow(RecADMM_reg_roftv[:,:,sliceSel],vmin=0, vmax=max_val) -plt.title('3D ADMM-ROF-TV Reconstruction, sagittal view') +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_roftv) -RMSE_admm_rof = Qtools.rmse() -print("Root Mean Square Error for ADMM-ROF-TV is {}".format(RMSE_admm_rof)) +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)) #%% -print ("Reconstructing with ADMM method using FGP-TV penalty") -RecADMM_reg_fgptv = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 30, \ - regularisation = 'FGP_TV', \ - regularisation_parameter = 0.005,\ - regularisation_iterations = 250) +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() -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.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_fgptv) -RMSE_admm_fgp = Qtools.rmse() -print("Root Mean Square Error for ADMM-FGP-TV is {}".format(RMSE_admm_fgp)) +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)) #%% print ("Reconstructing with ADMM method using TGV penalty") RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm, - rho_const = 2000.0, \ - iterationsADMM = 30, \ - regularisation = 'TGV', \ - regularisation_parameter = 0.005,\ - regularisation_iterations = 350) + 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() -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.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(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr") +ax1.plot([x0, x1], [y0, y1], 'ko-', linestyle='--') +plt.colorbar(ax=ax1) +plt.title('ADMM-TGV reconstruction, (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(RecADMM_reg_tgv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g') +plt.title('Profile', fontsize=19) +ax2 = plt.subplot(gs1[2]) +plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr") +plt.title('ADMM-TGV reconstruction, (Y-Z) view', fontsize=19) +ax2.set_aspect('equal') plt.show() +plt.savefig('TGV_phantom.pdf', format='pdf', dpi=1200) # calculate errors Qtools = QualityTools(phantom, RecADMM_reg_tgv) |