From 2d61438ff325b5a3ab26b6389042669a24276914 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Tue, 26 Feb 2019 17:51:13 +0000
Subject: demos_upd
---
.../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 36 ++++++++++++++++++++--
1 file changed, 33 insertions(+), 3 deletions(-)
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 d39d643..aa65cf3 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -23,6 +23,7 @@ import matplotlib.gridspec as gridspec
import numpy as np
import h5py
from ccpi.supp.qualitymetrics import QualityTools
+from scipy.signal import gaussian
# loading the data
h5f = h5py.File('data/TomoSim_data1550671417.h5','r')
@@ -136,12 +137,19 @@ 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)
+#plt.savefig('FBP_phantom.pdf', format='pdf', dpi=1600)
# calculate errors
Qtools = QualityTools(phantom, recFBP)
RMSE_fbp = Qtools.rmse()
print("Root Mean Square Error for FBP is {}".format(RMSE_fbp))
+
+# SSIM measure
+Qtools = QualityTools(phantom[128,:,:]*255, recFBP[128,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_fbp = Qtools.ssim(win2d)
+print("Mean SSIM for FBP is {}".format(ssim_fbp[0]))
#%%
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("Reconstructing with ADMM method using TomoRec software")
@@ -196,6 +204,13 @@ plt.savefig('SBTV_phantom.pdf', format='pdf', dpi=1600)
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))
+
+# SSIM measure
+Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_sbtv[128,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_admm_sbtv = Qtools.ssim(win2d)
+print("Mean SSIM ADMM-SBTV is {}".format(ssim_admm_sbtv[0]))
#%%
print ("Reconstructing with ADMM method using ROFLLT penalty")
RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm,
@@ -236,6 +251,13 @@ plt.show()
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))
+
+# SSIM measure
+Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_rofllt[128,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_admm_rifllt = Qtools.ssim(win2d)
+print("Mean SSIM ADMM-ROFLLT is {}".format(ssim_admm_rifllt[0]))
#%%
print ("Reconstructing with ADMM method using TGV penalty")
RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm,
@@ -244,7 +266,7 @@ RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm,
regularisation = 'TGV', \
regularisation_parameter = optimReg_tgv,\
regularisation_iterations = 600)
-
+#%%
sliceSel = int(0.5*N_size)
max_val = 1
plt.figure(figsize = (20,3))
@@ -269,10 +291,18 @@ plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
plt.title('ADMM-TGV (Y-Z) view', fontsize=19)
plt.colorbar(ax=ax4)
plt.show()
-plt.savefig('TGV_phantom.pdf', format='pdf', dpi=1200)
+#plt.savefig('TGV_phantom.pdf', format='pdf', dpi=1600)
# 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))
+
+# SSIM measure
+#Create a 2d gaussian for the window parameter
+Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_tgv[128,:,:]*235)
+win = np.array([gaussian(11, 1.5)])
+win2d = win * (win.T)
+ssim_admm_tgv = Qtools.ssim(win2d)
+print("Mean SSIM ADMM-TGV is {}".format(ssim_admm_tgv[0]))
#%%
--
cgit v1.2.3