optim demo upd, tgv 3d fix

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)