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