summaryrefslogtreecommitdiffstats
path: root/Wrappers/Python/demos
diff options
context:
space:
mode:
authorDaniil Kazantsev <dkazanc@hotmail.com>2019-02-24 22:48:55 +0000
committerDaniil Kazantsev <dkazanc@hotmail.com>2019-02-24 22:48:55 +0000
commit3c9af3e44778005948f6727c3e8280fcb6d3a3c9 (patch)
treeb916cb3a5039cd735706b127a3933f7ff2a6522f /Wrappers/Python/demos
parent380cb06910578f261ca5a7b3daae6fe42b1f24b0 (diff)
downloadregularization-3c9af3e44778005948f6727c3e8280fcb6d3a3c9.tar.gz
regularization-3c9af3e44778005948f6727c3e8280fcb6d3a3c9.tar.bz2
regularization-3c9af3e44778005948f6727c3e8280fcb6d3a3c9.tar.xz
regularization-3c9af3e44778005948f6727c3e8280fcb6d3a3c9.zip
demo upd
Diffstat (limited to 'Wrappers/Python/demos')
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py144
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