demo upd

1 files changed, 39 insertions, 105 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 9472a43..8a11f04 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
@@ -32,7 +32,7 @@ darks = h5f['darks'][:]
 angles_rad = h5f['angles_rad'][:]
 h5f.close()
 #%%
-# normalise the data, required format is [detectorsHoriz, Projections, Slices]
+# normalise the data [detectorsVert, Projections, detectorsHoriz]
 data_norm = normaliser(dataRaw, flats, darks, log='log')
 del dataRaw, darks, flats
 
@@ -42,141 +42,75 @@ plt.subplot(131)
 plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max)
 plt.title('2D Projection (analytical)')
 plt.subplot(132)
-plt.imshow(data_norm[:,:,300],vmin=0, vmax=intens_max)
+plt.imshow(data_norm[300,:,:],vmin=0, vmax=intens_max)
 plt.title('Sinogram view')
 plt.subplot(133)
-plt.imshow(data_norm[600,:,:],vmin=0, vmax=intens_max)
+plt.imshow(data_norm[:,:,600],vmin=0, vmax=intens_max)
 plt.title('Tangentogram view')
 plt.show()
 
 
-detectorHoriz = np.size(data_norm,0)
-N_size = 1000
-slice_to_recon = 0 # select which slice to reconstruct
+detectorHoriz = np.size(data_norm,2)
+det_y_crop = [i for i in range(0,detectorHoriz-22)]
+N_size = 950 # reconstruction domain
 #%%
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 from tomorec.methodsDIR import RecToolsDIR
-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
+                    DetectorsDimV = 10,  # 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                   
+                    ObjSize = N_size, # a scalar to define reconstructed object dimensions
                     device='gpu')
 
-FBPrec = RectoolsDIR.FBP(np.transpose(data_norm[det_y_crop,:,slice_to_recon]))
+FBPrec = RectoolsDIR.FBP(data_norm[0:10,:,det_y_crop])
 
 plt.figure()
-plt.imshow(FBPrec[150:550,150:550], vmin=0, vmax=0.005, cmap="gray")
+#plt.imshow(FBPrec[0,150:550,150:550], vmin=0, vmax=0.005, cmap="gray")
+plt.imshow(FBPrec[0,:,:], vmin=0, vmax=0.005, cmap="gray")
 plt.title('FBP reconstruction')
 
 #%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("Reconstructing with ADMM method using TomoRec software")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# initialise TomoRec ITERATIVE reconstruction class ONCE
 from tomorec.methodsIR import RecToolsIR
-# set parameters and initiate a class object
-Rectools = RecToolsIR(DetectorsDimH = np.size(det_y_crop),  # DetectorsDimH # detector dimension (horizontal)
-                    DetectorsDimV = None,  # DetectorsDimV # detector dimension (vertical) for 3D case only
+RectoolsIR = RecToolsIR(DetectorsDimH =  np.size(det_y_crop),  # DetectorsDimH # detector dimension (horizontal)
+                    DetectorsDimV = 5,  # 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='LS',# data fidelity, choose LS, PWLS, GH (wip), Student (wip)
+                    datafidelity='LS',# data fidelity, choose LS, PWLS (wip), 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
+                    OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
                     tolerance = 1e-08, # tolerance to stop outer iterations earlier
                     device='gpu')
-
-lc = Rectools.powermethod() # calculate Lipschitz constant (run once to initilise)
-
-RecFISTA_os_pwls = Rectools.FISTA(np.transpose(data_norm[det_y_crop,:,slice_to_recon]), \
-                             iterationsFISTA = 15, \
-                             lipschitz_const = lc)
-
-fig = plt.figure()
-plt.imshow(RecFISTA_os_pwls[150:550,150:550], vmin=0, vmax=0.003, cmap="gray")
-#plt.imshow(RecFISTA_os_pwls, vmin=0, vmax=0.004, cmap="gray")
-plt.title('FISTA PWLS-OS reconstruction')
-plt.show()
-#fig.savefig('dendr_PWLS.png', dpi=200)
 #%%
-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[det_y_crop,:,slice_to_recon]), \
-                              iterationsFISTA = 15, \
+print ("Reconstructing with ADMM method using ROF-TV penalty")
+
+RecADMM_reg_roftv = RectoolsIR.ADMM(data_norm[0:5,:,det_y_crop],
+                              rho_const = 2000.0, \
+                              iterationsADMM = 3, \
                               regularisation = 'FGP_TV', \
-                              regularisation_parameter = 0.000001,\
-                              regularisation_iterations = 200,\
-                              lipschitz_const = lc)
+                              regularisation_parameter = 0.001,\
+                              regularisation_iterations = 80)
 
-fig = plt.figure()
-plt.imshow(RecFISTA_pwls_os_TV[150:550,150:550], vmin=0, vmax=0.003, cmap="gray")
-#plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical')
-plt.title('FISTA PWLS-OS-TV reconstruction')
-plt.show()
-#fig.savefig('dendr_TV.png', dpi=200)
-#%%
-"""
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("Reconstructing with FISTA PWLS-OS-Diff4th method %%%%%%%%%%%")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-# Run FISTA-PWLS-OS reconstrucion algorithm with regularisation
-RecFISTA_pwls_os_Diff4th = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \
-                              np.transpose(dataRaw[:,:,slice_to_recon]), \
-                              iterationsFISTA = 15, \
-                              regularisation = 'DIFF4th', \
-                              regularisation_parameter = 0.1,\
-                              time_marching_parameter = 0.001,\
-                              edge_param = 0.003,\
-                              regularisation_iterations = 600,\
-                              lipschitz_const = lc)
 
-fig = plt.figure()
-plt.imshow(RecFISTA_pwls_os_Diff4th[150:550,150:550], vmin=0, vmax=0.004, cmap="gray")
-#plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical')
-plt.title('FISTA PWLS-OS-Diff4th reconstruction')
-plt.show()
-#fig.savefig('dendr_Diff4th.png', dpi=200)
-"""
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("Reconstructing with FISTA PWLS-OS-ROF_LLT method %%%%%%%%%%%")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-# Run FISTA-PWLS-OS reconstrucion algorithm with regularisation
-RecFISTA_pwls_os_rofllt = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \
-                              np.transpose(dataRaw[:,:,slice_to_recon]), \
-                              iterationsFISTA = 15, \
-                              regularisation = 'LLT_ROF', \
-                              regularisation_parameter = 0.000007,\
-                              regularisation_parameter2 = 0.0004,\
-                              regularisation_iterations = 350,\
-                              lipschitz_const = lc)
+sliceSel = 2
+max_val = 0.005
+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')
 
-fig = plt.figure()
-plt.imshow(RecFISTA_pwls_os_rofllt[150:550,150:550], vmin=0, vmax=0.003, cmap="gray")
-#plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical')
-plt.title('FISTA PWLS-OS-ROF-LLT reconstruction')
-plt.show()
-#fig.savefig('dendr_ROFLLT.png', dpi=200)
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("Reconstructing with FISTA PWLS-OS-TGV method %%%%%%%%%%%%%%%")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-# Run FISTA-PWLS-OS reconstrucion algorithm with regularisation
-RecFISTA_pwls_os_tgv = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \
-                              np.transpose(dataRaw[:,:,slice_to_recon]), \
-                              iterationsFISTA = 15, \
-                              regularisation = 'TGV', \
-                              regularisation_parameter = 0.001,\
-                              TGV_alpha2 = 0.001,\
-                              TGV_alpha1 = 2.0,\
-                              TGV_LipschitzConstant = 24,\
-                              regularisation_iterations = 300,\
-                              lipschitz_const = lc)
+plt.subplot(132)
+plt.imshow(RecADMM_reg_roftv[:,sliceSel,:],vmin=0, vmax=max_val)
+plt.title('3D ADMM-ROF-TV Reconstruction, coronal view')
 
-fig = plt.figure()
-plt.imshow(RecFISTA_pwls_os_tgv[150:550,150:550], vmin=0, vmax=0.003, cmap="gray")
-#plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical')
-plt.title('FISTA PWLS-OS-TGV reconstruction')
+plt.subplot(133)
+plt.imshow(RecADMM_reg_roftv[:,:,sliceSel],vmin=0, vmax=max_val)
+plt.title('3D ADMM-ROF-TV Reconstruction, sagittal view')
 plt.show()
-#fig.savefig('dendr_TGV.png', dpi=200)
+#%%
\ No newline at end of file