summaryrefslogtreecommitdiffstats
path: root/Wrappers/Python/demos
diff options
context:
space:
mode:
authordkazanc <dkazanc@hotmail.com>2019-02-22 16:23:30 +0000
committerdkazanc <dkazanc@hotmail.com>2019-02-22 16:23:30 +0000
commit380cb06910578f261ca5a7b3daae6fe42b1f24b0 (patch)
treeb6093364f42f5b168a1b4f9d49740e4d61f41846 /Wrappers/Python/demos
parent3679e589cde1bdcede504b96d2a0df053289dbe4 (diff)
downloadregularization-380cb06910578f261ca5a7b3daae6fe42b1f24b0.tar.gz
regularization-380cb06910578f261ca5a7b3daae6fe42b1f24b0.tar.bz2
regularization-380cb06910578f261ca5a7b3daae6fe42b1f24b0.tar.xz
regularization-380cb06910578f261ca5a7b3daae6fe42b1f24b0.zip
demos_updated2
Diffstat (limited to 'Wrappers/Python/demos')
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py48
-rw-r--r--Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py224
2 files changed, 120 insertions, 152 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 02e7c5c..9472a43 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
@@ -14,37 +14,29 @@ ____________________________________________________________________________
1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox
2. TomoRec: conda install -c dkazanc tomorec
or install from https://github.com/dkazanc/TomoRec
-3. TomoPhantom software for data generation
@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
GPLv3 license (ASTRA toolbox)
"""
import numpy as np
import matplotlib.pyplot as plt
-import scipy.io
+import h5py
from tomorec.supp.suppTools import normaliser
-# load dendritic data
-datadict = scipy.io.loadmat('Rawdata_1_frame3D.mat')
-# extract data (print(datadict.keys()))
-dataRaw = datadict['Rawdata3D']
-angles = datadict['AnglesAr']
-flats = datadict['flats3D']
-darks= datadict['darks3D']
-dataRaw = np.swapaxes(dataRaw,0,1)
+# load dendritic projection data
+h5f = h5py.File('data/DendrData_3D.h5','r')
+dataRaw = h5f['dataRaw'][:]
+flats = h5f['flats'][:]
+darks = h5f['darks'][:]
+angles_rad = h5f['angles_rad'][:]
+h5f.close()
#%%
-#flats2 = np.zeros((np.size(flats,0),1, np.size(flats,1)), dtype='float32')
-#flats2[:,0,:] = flats[:]
-#darks2 = np.zeros((np.size(darks,0),1, np.size(darks,1)), dtype='float32')
-#darks2[:,0,:] = darks[:]
-
# normalise the data, required format is [detectorsHoriz, Projections, Slices]
data_norm = normaliser(dataRaw, flats, darks, log='log')
+del dataRaw, darks, flats
-#dataRaw = np.float32(np.divide(dataRaw, np.max(dataRaw).astype(float)))
-
-intens_max = 70
+intens_max = 2
plt.figure()
plt.subplot(131)
plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max)
@@ -61,40 +53,41 @@ plt.show()
detectorHoriz = np.size(data_norm,0)
N_size = 1000
slice_to_recon = 0 # select which slice to reconstruct
-angles_rad = angles*(np.pi/180.0)
#%%
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%")
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
from tomorec.methodsDIR import RecToolsDIR
-RectoolsDIR = RecToolsDIR(DetectorsDimH = detectorHoriz, # DetectorsDimH # detector dimension (horizontal)
+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
AnglesVec = angles_rad, # array of angles in radians
ObjSize = N_size, # a scalar to define reconstructed object dimensions
device='gpu')
-FBPrec = RectoolsDIR.FBP(np.transpose(data_norm[:,:,slice_to_recon]))
+FBPrec = RectoolsDIR.FBP(np.transpose(data_norm[det_y_crop,:,slice_to_recon]))
plt.figure()
plt.imshow(FBPrec[150:550,150:550], vmin=0, vmax=0.005, cmap="gray")
plt.title('FBP reconstruction')
+#%%
from tomorec.methodsIR import RecToolsIR
# set parameters and initiate a class object
-Rectools = RecToolsIR(DetectorsDimH = detectorHoriz, # DetectorsDimH # detector dimension (horizontal)
+Rectools = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal)
DetectorsDimV = None, # 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='PWLS',# data fidelity, choose LS, PWLS, GH (wip), Student (wip)
+ datafidelity='LS',# data fidelity, choose LS, PWLS, 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
tolerance = 1e-08, # tolerance to stop outer iterations earlier
device='gpu')
-lc = Rectools.powermethod(np.transpose(dataRaw[:,:,slice_to_recon])) # calculate Lipschitz constant (run once to initilise)
+lc = Rectools.powermethod() # calculate Lipschitz constant (run once to initilise)
-RecFISTA_os_pwls = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \
- np.transpose(dataRaw[:,:,slice_to_recon]), \
+RecFISTA_os_pwls = Rectools.FISTA(np.transpose(data_norm[det_y_crop,:,slice_to_recon]), \
iterationsFISTA = 15, \
lipschitz_const = lc)
@@ -109,8 +102,7 @@ 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[:,:,slice_to_recon]), \
- np.transpose(dataRaw[:,:,slice_to_recon]), \
+RecFISTA_pwls_os_TV = Rectools.FISTA(np.transpose(data_norm[det_y_crop,:,slice_to_recon]), \
iterationsFISTA = 15, \
regularisation = 'FGP_TV', \
regularisation_parameter = 0.000001,\
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
index 467e810..148fdcc 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
@@ -11,136 +11,59 @@ ____________________________________________________________________________
Demo_SimulData_Recon_SX.py
____________________________________________________________________________
>>>>> Dependencies: <<<<<
+>>>>> Dependencies: <<<<<
1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox
2. TomoRec: conda install -c dkazanc tomorec
or install from https://github.com/dkazanc/TomoRec
-3. TomoPhantom software for data generation
@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
GPLv3 license (ASTRA toolbox)
"""
-import timeit
-import os
+#import timeit
import matplotlib.pyplot as plt
import numpy as np
-import tomophantom
-from tomophantom import TomoP3D
-from tomophantom.supp.qualitymetrics import QualityTools
-from tomophantom.supp.flatsgen import flats
-from tomophantom.supp.normraw import normaliser_sim
+import h5py
+from ccpi.supp.qualitymetrics import QualityTools
+
+# 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()
-print ("Building 3D phantom using TomoPhantom software")
-tic=timeit.default_timer()
-model = 9 # select a model number from the library
-N_size = 128 # Define phantom dimensions using a scalar value (cubic phantom)
-path = os.path.dirname(tomophantom.__file__)
-path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
-#This will generate a N_size x N_size x N_size phantom (3D)
-phantom_tm = TomoP3D.Model(model, N_size, path_library3D)
-toc=timeit.default_timer()
-Run_time = toc - tic
-print("Phantom has been built in {} seconds".format(Run_time))
+[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm)
+N_size = Vert_det
-sliceSel = int(0.5*N_size)
+sliceSel = 128
#plt.gray()
plt.figure()
plt.subplot(131)
-plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1)
+plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1)
plt.title('3D Phantom, axial view')
plt.subplot(132)
-plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1)
+plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1)
plt.title('3D Phantom, coronal view')
plt.subplot(133)
-plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1)
+plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1)
plt.title('3D Phantom, sagittal view')
plt.show()
-# Projection geometry related parameters:
-Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal)
-Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
-angles_num = int(0.25*np.pi*N_size); # angles number
-angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees
-angles_rad = angles*(np.pi/180.0)
-#%%
-print ("Building 3D analytical projection data with TomoPhantom")
-projData3D_analyt= TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D)
-
-intens_max = 70
-sliceSel = 150
+intens_max = 240
plt.figure()
plt.subplot(131)
-plt.imshow(projData3D_analyt[:,sliceSel,:],vmin=0, vmax=intens_max)
-plt.title('2D Projection (analytical)')
-plt.subplot(132)
-plt.imshow(projData3D_analyt[sliceSel,:,:],vmin=0, vmax=intens_max)
-plt.title('Sinogram view')
-plt.subplot(133)
-plt.imshow(projData3D_analyt[:,:,sliceSel],vmin=0, vmax=intens_max)
-plt.title('Tangentogram view')
-plt.show()
-#%%
-print ("Simulate flat fields, add noise and normalise projections...")
-flatsnum = 20 # generate 20 flat fields
-flatsSIM = flats(Vert_det, Horiz_det, maxheight = 0.1, maxthickness = 3, sigma_noise = 0.2, sigmasmooth = 3, flatsnum=flatsnum)
-
-plt.figure()
-plt.imshow(flatsSIM[0,:,:],vmin=0, vmax=1)
-plt.title('A selected simulated flat-field')
-#%%
-# Apply normalisation of data and add noise
-flux_intensity = 10000 # controls the level of noise
-sigma_flats = 0.085 # contro the level of noise in flats (higher creates more ring artifacts)
-projData3D_norm = normaliser_sim(projData3D_analyt, flatsSIM, sigma_flats, flux_intensity)
-
-intens_max = 70
-sliceSel = 150
-plt.figure()
-plt.subplot(131)
-plt.imshow(projData3D_norm[:,sliceSel,:],vmin=0, vmax=intens_max)
+plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max)
plt.title('2D Projection (erroneous)')
plt.subplot(132)
-plt.imshow(projData3D_norm[sliceSel,:,:],vmin=0, vmax=intens_max)
+plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max)
plt.title('Sinogram view')
plt.subplot(133)
-plt.imshow(projData3D_norm[:,:,sliceSel],vmin=0, vmax=intens_max)
+plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max)
plt.title('Tangentogram view')
plt.show()
#%%
-# initialise TomoRec DIRECT reconstruction class ONCE
-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
- 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
-
-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.title('3D Reconstruction, axial view')
-
-plt.subplot(132)
-plt.imshow(recNumerical[:,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.title('3D Reconstruction, sagittal view')
-plt.show()
-
-# calculate errors
-Qtools = QualityTools(phantom_tm, recNumerical)
-RMSE_fbp = Qtools.rmse()
-print("Root Mean Square Error for FBP is {}".format(RMSE_fbp))
-#%%
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("Reconstructing with ADMM method using TomoRec software")
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -148,7 +71,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')
@@ -156,33 +79,86 @@ 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,
- rho_const = 2000.0, \
- iterationsADMM = 30, \
- regularisation = 'FGP_TV', \
- regularisation_parameter = 0.003,\
- regularisation_iterations = 250)
+param_space = 20
+reg_param_sb_vec = np.linspace(0.001,0.03,param_space,dtype='float32') # a vector of parameters
+erros_vec_sbtv = np.zeros((param_space)) # a vector of errors
+
+print ("Reconstructing with ADMM method using SB-TV penalty")
+for i in range(0,param_space):
+ RecADMM_reg_sbtv = RectoolsIR.ADMM(projdata_norm,
+ rho_const = 2000.0, \
+ iterationsADMM = 15, \
+ regularisation = 'SB_TV', \
+ regularisation_parameter = reg_param_sb_vec[i],\
+ regularisation_iterations = 50)
+ # calculate errors
+ Qtools = QualityTools(phantom, RecADMM_reg_sbtv)
+ erros_vec_sbtv[i] = Qtools.rmse()
+ print("RMSE for regularisation parameter {} for ADMM-SB-TV is {}".format(reg_param_sb_vec[i],erros_vec_sbtv[i]))
-sliceSel = int(0.5*N_size)
-max_val = 1
plt.figure()
-plt.subplot(131)
-plt.imshow(RecADMM_reg_fgptv[sliceSel,:,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-FGP-TV Reconstruction, axial view')
-
-plt.subplot(132)
-plt.imshow(RecADMM_reg_fgptv[:,sliceSel,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-FGP-TV Reconstruction, coronal view')
-
-plt.subplot(133)
-plt.imshow(RecADMM_reg_fgptv[:,:,sliceSel],vmin=0, vmax=max_val)
-plt.title('3D ADMM-FGP-TV Reconstruction, sagittal view')
-plt.show()
+plt.plot(erros_vec_sbtv)
+#%%
+param_space = 20
+reg_param_rofllt_vec = np.linspace(0.001,0.03,param_space,dtype='float32') # a vector of parameters
+erros_vec_rofllt = np.zeros((param_space)) # a vector of errors
+
+print ("Reconstructing with ADMM method using ROF-LLT penalty")
+for i in range(0,param_space):
+ RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm,
+ rho_const = 2000.0, \
+ iterationsADMM = 15, \
+ regularisation = 'LLT_ROF', \
+ regularisation_parameter = reg_param_rofllt_vec[i],\
+ regularisation_parameter2 = 0.005,\
+ regularisation_iterations = 600)
+ # calculate errors
+ Qtools = QualityTools(phantom, RecADMM_reg_rofllt)
+ erros_vec_rofllt[i] = Qtools.rmse()
+ print("RMSE for regularisation parameter {} for ADMM-ROF-LLT is {}".format(reg_param_rofllt_vec[i],erros_vec_rofllt[i]))
-# calculate errors
-Qtools = QualityTools(phantom_tm, RecADMM_reg_fgptv)
-RMSE_admm_fgp = Qtools.rmse()
-print("Root Mean Square Error for ADMM-FGP-TV is {}".format(RMSE_admm_fgp))
+plt.figure()
+plt.plot(erros_vec_rofllt)
+#%%
+param_space = 20
+reg_param_tgv_vec = np.linspace(0.001,0.03,param_space,dtype='float32') # a vector of parameters
+erros_vec_tgv = np.zeros((param_space)) # a vector of errors
+
+print ("Reconstructing with ADMM method using TGV penalty")
+for i in range(0,param_space):
+ RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm,
+ rho_const = 2000.0, \
+ iterationsADMM = 15, \
+ regularisation = 'TGV', \
+ regularisation_parameter = reg_param_tgv_vec[i],\
+ regularisation_iterations = 600)
+ # calculate errors
+ Qtools = QualityTools(phantom, RecADMM_reg_tgv)
+ erros_vec_tgv[i] = Qtools.rmse()
+ print("RMSE for regularisation parameter {} for ADMM-TGV is {}".format(reg_param_tgv_vec[i],erros_vec_tgv[i]))
+plt.figure()
+plt.plot(erros_vec_tgv)
#%%
+param_space = 1
+reg_param_diff4th = np.linspace(10,100,param_space,dtype='float32') # a vector of parameters
+erros_vec_diff4th = np.zeros((param_space)) # a vector of errors
+
+print ("Reconstructing with ADMM method using Diff4th penalty")
+for i in range(0,param_space):
+ RecADMM_reg_diff4th = RectoolsIR.ADMM(projdata_norm,
+ rho_const = 2000.0, \
+ iterationsADMM = 15, \
+ regularisation = 'Diff4th', \
+ regularisation_parameter = reg_param_diff4th[i],\
+ edge_param = 0.03,
+ time_marching_parameter = 0.001,\
+ regularisation_iterations = 750)
+ # calculate errors
+ Qtools = QualityTools(phantom, RecADMM_reg_diff4th)
+ erros_vec_diff4th[i] = Qtools.rmse()
+ print("RMSE for regularisation parameter {} for ADMM-Diff4th is {}".format(reg_param_diff4th[i],erros_vec_diff4th[i]))
+
+plt.figure()
+plt.plot(erros_vec_diff4th)
+#%% \ No newline at end of file