From 15040fcdb083db496dd42dbee216e983726620bf Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Wed, 20 Feb 2019 17:22:33 +0000
Subject: rec_demo
---
.../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 114 ++++++++++++++++++---
1 file changed, 101 insertions(+), 13 deletions(-)
(limited to 'Wrappers')
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 7d22bfd..a9e45b6 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -21,15 +21,17 @@ import timeit
import matplotlib.pyplot as plt
import numpy as np
import h5py
-from tomophantom.supp.qualitymetrics import QualityTools
+from ccpi.supp.qualitymetrics import QualityTools
-# loading data
+# loading the data
h5f = h5py.File('data/TomoSim_data1550671417.h5','r')
phantom = h5f['phantom'][:]
projdata_norm = h5f['projdata_norm'][:]
proj_angles = h5f['proj_angles'][:]
h5f.close()
+[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm)
+N_size = Vert_det
sliceSel = 128
#plt.gray()
@@ -64,32 +66,32 @@ plt.show()
from tomorec.methodsDIR import RecToolsDIR
RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal)
DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only
- AnglesVec = angles_rad, # array of angles in radians
+ AnglesVec = proj_angles, # array of angles in radians
ObjSize = N_size, # a scalar to define reconstructed object dimensions
device = 'gpu')
#%%
print ("Reconstruction using FBP from TomoRec")
-recNumerical= RectoolsDIR.FBP(projData3D_norm) # FBP reconstruction
+recFBP= RectoolsDIR.FBP(projdata_norm) # FBP reconstruction
sliceSel = int(0.5*N_size)
max_val = 1
#plt.gray()
plt.figure()
plt.subplot(131)
-plt.imshow(recNumerical[sliceSel,:,:],vmin=0, vmax=max_val)
+plt.imshow(recFBP[sliceSel,:,:],vmin=0, vmax=max_val)
plt.title('3D Reconstruction, axial view')
plt.subplot(132)
-plt.imshow(recNumerical[:,sliceSel,:],vmin=0, vmax=max_val)
+plt.imshow(recFBP[:,sliceSel,:],vmin=0, vmax=max_val)
plt.title('3D Reconstruction, coronal view')
plt.subplot(133)
-plt.imshow(recNumerical[:,:,sliceSel],vmin=0, vmax=max_val)
+plt.imshow(recFBP[:,:,sliceSel],vmin=0, vmax=max_val)
plt.title('3D Reconstruction, sagittal view')
plt.show()
# calculate errors
-Qtools = QualityTools(phantom_tm, recNumerical)
+Qtools = QualityTools(phantom, recFBP)
RMSE_fbp = Qtools.rmse()
print("Root Mean Square Error for FBP is {}".format(RMSE_fbp))
#%%
@@ -100,7 +102,7 @@ print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
from tomorec.methodsIR import RecToolsIR
RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal)
DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only
- AnglesVec = angles_rad, # array of angles in radians
+ AnglesVec = proj_angles, # array of angles in radians
ObjSize = N_size, # a scalar to define reconstructed object dimensions
datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
@@ -108,12 +110,41 @@ RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector d
tolerance = 1e-08, # tolerance to stop outer iterations earlier
device='gpu')
#%%
-# Run ADMM reconstrucion algorithm with 3D regularisation
-RecADMM_reg_fgptv = RectoolsIR.ADMM(projData3D_norm,
+print ("Reconstructing with ADMM method using ROF-TV penalty")
+RecADMM_reg_roftv = RectoolsIR.ADMM(projdata_norm,
rho_const = 2000.0, \
iterationsADMM = 30, \
- regularisation = 'FGP_TV', \
+ regularisation = 'ROF_TV', \
regularisation_parameter = 0.003,\
+ regularisation_iterations = 500)
+
+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.show()
+
+# 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))
+#%%
+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)
sliceSel = int(0.5*N_size)
@@ -133,8 +164,65 @@ plt.title('3D ADMM-FGP-TV Reconstruction, sagittal view')
plt.show()
# calculate errors
-Qtools = QualityTools(phantom_tm, RecADMM_reg_fgptv)
+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))
+#%%
+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)
+
+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.show()
+
+# calculate errors
+Qtools = QualityTools(phantom, RecADMM_reg_tgv)
+RMSE_admm_tgv = Qtools.rmse()
+print("Root Mean Square Error for ADMM-TGV is {}".format(RMSE_admm_tgv))
+#%%
+print ("Reconstructing with ADMM method using Diff4th penalty")
+RecADMM_reg_diff4th = RectoolsIR.ADMM(projdata_norm,
+ rho_const = 2000.0, \
+ iterationsADMM = 30, \
+ regularisation = 'Diff4th', \
+ regularisation_parameter = 0.005,\
+ regularisation_iterations = 500)
+sliceSel = int(0.5*N_size)
+max_val = 1
+plt.figure()
+plt.subplot(131)
+plt.imshow(RecADMM_reg_diff4th[sliceSel,:,:],vmin=0, vmax=max_val)
+plt.title('3D ADMM-Diff4th Reconstruction, axial view')
+
+plt.subplot(132)
+plt.imshow(RecADMM_reg_diff4th[:,sliceSel,:],vmin=0, vmax=max_val)
+plt.title('3D ADMM-Diff4th Reconstruction, coronal view')
+
+plt.subplot(133)
+plt.imshow(RecADMM_reg_diff4th[:,:,sliceSel],vmin=0, vmax=max_val)
+plt.title('3D ADMM-Diff4th Reconstruction, sagittal view')
+plt.show()
+
+# calculate errors
+Qtools = QualityTools(phantom, RecADMM_reg_diff4th)
+RMSE_admm_diff4th = Qtools.rmse()
+print("Root Mean Square Error for ADMM-Diff4th is {}".format(RMSE_admm_diff4th))
#%%
--
cgit v1.2.3