summaryrefslogtreecommitdiffstats
path: root/Wrappers/Python/demos
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers/Python/demos')
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py2
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py259
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)
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-30 10:52:35 +0000