From 96c5d86fe9a9ac9725ae64d1c50ec9af71137608 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Wed, 20 Feb 2019 11:09:21 +0000
Subject: demos created
---
.../demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 190 +++++++++++++++++++++
.../SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 188 ++++++++++++++++++++
.../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 183 ++++++++++++++++++++
.../demos/SoftwareX_supp/Demo_SimulData_SX.py | 63 +++++++
Wrappers/Python/demos/SoftwareX_supp/Readme.md | 19 +++
5 files changed, 643 insertions(+)
create mode 100644 Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
create mode 100644 Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
create mode 100644 Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
create mode 100644 Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
create mode 100644 Wrappers/Python/demos/SoftwareX_supp/Readme.md
(limited to 'Wrappers')
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
new file mode 100644
index 0000000..02e7c5c
--- /dev/null
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
@@ -0,0 +1,190 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This demo scripts support the following publication:
+"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham,
+Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019
+____________________________________________________________________________
+* Reads real tomographic data (stored at Zenodo)
+* Reconstructs using TomoRec software
+* Saves reconstructed images
+____________________________________________________________________________
+>>>>> 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 numpy as np
+import matplotlib.pyplot as plt
+import scipy.io
+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)
+#%%
+#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')
+
+#dataRaw = np.float32(np.divide(dataRaw, np.max(dataRaw).astype(float)))
+
+intens_max = 70
+plt.figure()
+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.title('Sinogram view')
+plt.subplot(133)
+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
+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)
+ 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]))
+
+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)
+ 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)
+ 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)
+
+RecFISTA_os_pwls = Rectools.FISTA(np.transpose(data_norm[:,:,slice_to_recon]), \
+ np.transpose(dataRaw[:,:,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[:,:,slice_to_recon]), \
+ np.transpose(dataRaw[:,:,slice_to_recon]), \
+ iterationsFISTA = 15, \
+ regularisation = 'FGP_TV', \
+ regularisation_parameter = 0.000001,\
+ regularisation_iterations = 200,\
+ lipschitz_const = lc)
+
+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)
+
+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)
+
+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.show()
+#fig.savefig('dendr_TGV.png', dpi=200)
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
new file mode 100644
index 0000000..467e810
--- /dev/null
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This demo scripts support the following publication:
+"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham,
+Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019
+____________________________________________________________________________
+* Reads data which is previosly generated by TomoPhantom software (Zenodo link)
+* Optimises for the regularisation parameters which later used in the script:
+Demo_SimulData_Recon_SX.py
+____________________________________________________________________________
+>>>>> 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 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
+
+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))
+
+sliceSel = int(0.5*N_size)
+#plt.gray()
+plt.figure()
+plt.subplot(131)
+plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1)
+plt.title('3D Phantom, axial view')
+
+plt.subplot(132)
+plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1)
+plt.title('3D Phantom, coronal view')
+
+plt.subplot(133)
+plt.imshow(phantom_tm[:,:,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
+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.title('2D Projection (erroneous)')
+plt.subplot(132)
+plt.imshow(projData3D_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.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 ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# initialise TomoRec ITERATIVE reconstruction class ONCE
+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
+ 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')
+ OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
+ 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)
+
+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()
+
+# 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))
+
+#%%
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
new file mode 100644
index 0000000..0a99951
--- /dev/null
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -0,0 +1,183 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This demo scripts support the following publication:
+"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham,
+Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019
+____________________________________________________________________________
+* Reads data which is previously generated by TomoPhantom software (Zenodo link)
+* Reconstruct using optimised regularisation parameters (see Demo_SimulData_ParOptimis_SX.py)
+____________________________________________________________________________
+>>>>> 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 matplotlib.pyplot as plt
+import numpy as np
+from tomophantom.supp.qualitymetrics import QualityTools
+
+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))
+
+sliceSel = int(0.5*N_size)
+#plt.gray()
+plt.figure()
+plt.subplot(131)
+plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1)
+plt.title('3D Phantom, axial view')
+
+plt.subplot(132)
+plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1)
+plt.title('3D Phantom, coronal view')
+
+plt.subplot(133)
+plt.imshow(phantom_tm[:,:,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
+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.title('2D Projection (erroneous)')
+plt.subplot(132)
+plt.imshow(projData3D_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.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 ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# initialise TomoRec ITERATIVE reconstruction class ONCE
+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
+ 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')
+ OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
+ 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)
+
+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()
+
+# 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))
+
+#%%
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
new file mode 100644
index 0000000..3265d37
--- /dev/null
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This demo scripts support the following publication:
+"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham,
+Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019
+____________________________________________________________________________
+* Runs TomoPhantom software to simulate tomographic projection data with
+some imaging errors and noise
+* Saves the data into hdf file to be uploaded in reconstruction scripts
+____________________________________________________________________________
+
+>>>>> Dependencies: <<<<<
+1. TomoPhantom software for data generation
+
+@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
+Apache 2.0 license
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.io
+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)
+#%%
+#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')
+
+#dataRaw = np.float32(np.divide(dataRaw, np.max(dataRaw).astype(float)))
+
+intens_max = 70
+plt.figure()
+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.title('Sinogram view')
+plt.subplot(133)
+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
+angles_rad = angles*(np.pi/180.0)
+#%%
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Readme.md b/Wrappers/Python/demos/SoftwareX_supp/Readme.md
new file mode 100644
index 0000000..fa77745
--- /dev/null
+++ b/Wrappers/Python/demos/SoftwareX_supp/Readme.md
@@ -0,0 +1,19 @@
+# SoftwareX publication [1] supporting files
+
+## Decription:
+The scripts here support publication in SoftwareX journal [1] to ensure
+reproducibility of the research.
+
+## Data:
+
+## Dependencies:
+
+## Files description:
+
+
+### References:
+[1] "CCPi-Regularisation Toolkit for computed tomographic image reconstruction with proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham, Martin J. Turner, Philip J. Withers and Alun Ashton
+
+
+### Acknowledgments:
+CCPi-RGL software is a product of the [CCPi](https://www.ccpi.ac.uk/) group, STFC SCD software developers and Diamond Light Source (DLS). Any relevant questions/comments can be e-mailed to Daniil Kazantsev at dkazanc@hotmail.com
--
cgit v1.2.3
From 808860a03e254bf52fe5c6b5f3df883ccd62d714 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Wed, 20 Feb 2019 14:31:50 +0000
Subject: progress sim data
---
.../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 75 +++-----------
.../demos/SoftwareX_supp/Demo_SimulData_SX.py | 114 +++++++++++++++------
2 files changed, 100 insertions(+), 89 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 0a99951..7d22bfd 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -13,93 +13,50 @@ ____________________________________________________________________________
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 matplotlib.pyplot as plt
import numpy as np
+import h5py
from tomophantom.supp.qualitymetrics import QualityTools
-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))
+# loading data
+h5f = h5py.File('data/TomoSim_data1550671417.h5','r')
+phantom = h5f['phantom'][:]
+projdata_norm = h5f['projdata_norm'][:]
+proj_angles = h5f['proj_angles'][:]
+h5f.close()
-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
-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
+intens_max = 240
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()
#%%
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
index 3265d37..ce29f0c 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
@@ -12,52 +12,106 @@ some imaging errors and noise
____________________________________________________________________________
>>>>> Dependencies: <<<<<
-1. TomoPhantom software for data generation
+1. TomoPhantom software for phantom and data generation
@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
Apache 2.0 license
"""
-import numpy as np
+import timeit
+import os
import matplotlib.pyplot as plt
-import scipy.io
-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)
-#%%
-#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[:]
+import numpy as np
+import tomophantom
+from tomophantom import TomoP3D
+from tomophantom.supp.flatsgen import flats
+from tomophantom.supp.normraw import normaliser_sim
+
+print ("Building 3D phantom using TomoPhantom software")
+tic=timeit.default_timer()
+model = 16 # select a model number from the library
+N_size = 256 # 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))
+
+sliceSel = int(0.5*N_size)
+#plt.gray()
+plt.figure()
+plt.subplot(131)
+plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1)
+plt.title('3D Phantom, axial view')
+
+plt.subplot(132)
+plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1)
+plt.title('3D Phantom, coronal view')
-# normalise the data, required format is [detectorsHoriz, Projections, Slices]
-data_norm = normaliser(dataRaw, flats, darks, log='log')
+plt.subplot(133)
+plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1)
+plt.title('3D Phantom, sagittal view')
+plt.show()
-#dataRaw = np.float32(np.divide(dataRaw, np.max(dataRaw).astype(float)))
+# 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.35*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
+intens_max = N_size
+sliceSel = int(0.5*N_size)
plt.figure()
plt.subplot(131)
-plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max)
+plt.imshow(projData3D_analyt[:,sliceSel,:],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(projData3D_analyt[sliceSel,:,:],vmin=0, vmax=intens_max)
plt.title('Sinogram view')
plt.subplot(133)
-plt.imshow(data_norm[600,:,:],vmin=0, vmax=intens_max)
+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 = 60000 # controls the level of noise
+sigma_flats = 0.01 # contro the level of noise in flats (higher creates more ring artifacts)
+projData3D_norm = normaliser_sim(projData3D_analyt, flatsSIM, sigma_flats, flux_intensity)
-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)
+intens_max = N_size
+sliceSel = int(0.5*N_size)
+plt.figure()
+plt.subplot(131)
+plt.imshow(projData3D_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.title('Sinogram view')
+plt.subplot(133)
+plt.imshow(projData3D_norm[:,:,sliceSel],vmin=0, vmax=intens_max)
+plt.title('Tangentogram view')
+plt.show()
#%%
+import h5py
+import time
+time_label = int(time.time())
+# Saving generated data with a unique time label
+h5f = h5py.File('TomoSim_data'+str(time_label)+'.h5', 'w')
+h5f.create_dataset('phantom', data=phantom_tm)
+h5f.create_dataset('projdata_norm', data=projData3D_norm)
+h5f.create_dataset('proj_angles', data=angles_rad)
+h5f.close()
+#%%
\ No newline at end of file
--
cgit v1.2.3
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
From 3679e589cde1bdcede504b96d2a0df053289dbe4 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Fri, 22 Feb 2019 10:31:06 +0000
Subject: dem_upd
---
Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py | 6 +++---
1 file changed, 3 insertions(+), 3 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 a9e45b6..5dbd436 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -17,7 +17,7 @@ or install from https://github.com/dkazanc/TomoRec
@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
GPLv3 license (ASTRA toolbox)
"""
-import timeit
+#import timeit
import matplotlib.pyplot as plt
import numpy as np
import h5py
@@ -202,8 +202,8 @@ RecADMM_reg_diff4th = RectoolsIR.ADMM(projdata_norm,
rho_const = 2000.0, \
iterationsADMM = 30, \
regularisation = 'Diff4th', \
- regularisation_parameter = 0.005,\
- regularisation_iterations = 500)
+ regularisation_parameter = 0.0005,\
+ regularisation_iterations = 200)
sliceSel = int(0.5*N_size)
max_val = 1
--
cgit v1.2.3
From 380cb06910578f261ca5a7b3daae6fe42b1f24b0 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Fri, 22 Feb 2019 16:23:30 +0000
Subject: demos_updated2
---
.../demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 48 ++---
.../SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 224 +++++++++------------
2 files changed, 120 insertions(+), 152 deletions(-)
(limited to 'Wrappers')
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
--
cgit v1.2.3
From 3c9af3e44778005948f6727c3e8280fcb6d3a3c9 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Sun, 24 Feb 2019 22:48:55 +0000
Subject: demo upd
---
.../demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 144 ++++++---------------
1 file changed, 39 insertions(+), 105 deletions(-)
(limited to 'Wrappers')
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
--
cgit v1.2.3
From dc5a92771dcf9bfd262ba34b0fc8a1c2df40897d Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Mon, 25 Feb 2019 10:07:50 +0000
Subject: demo upf
---
.../SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 34 ++++------------------
1 file changed, 6 insertions(+), 28 deletions(-)
(limited to 'Wrappers')
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 148fdcc..a79d0a3 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
@@ -79,8 +79,8 @@ RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector d
tolerance = 1e-08, # tolerance to stop outer iterations earlier
device='gpu')
#%%
-param_space = 20
-reg_param_sb_vec = np.linspace(0.001,0.03,param_space,dtype='float32') # a vector of parameters
+param_space = 30
+reg_param_sb_vec = np.linspace(0.03,0.15,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")
@@ -99,8 +99,8 @@ for i in range(0,param_space):
plt.figure()
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
+param_space = 30
+reg_param_rofllt_vec = np.linspace(0.03,0.15,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")
@@ -120,8 +120,8 @@ for i in range(0,param_space):
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
+param_space = 30
+reg_param_tgv_vec = np.linspace(0.03,0.15,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")
@@ -139,26 +139,4 @@ for i in range(0,param_space):
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
--
cgit v1.2.3
From 0587f96cafac66d9ee04005a80b43514c6d2a753 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Mon, 25 Feb 2019 17:13:28 +0000
Subject: some TGV updates and demos
---
.../demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 148 ++++++++++++++++++---
.../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 29 ----
2 files changed, 127 insertions(+), 50 deletions(-)
(limited to 'Wrappers')
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 8a11f04..4cd680e 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
@@ -22,7 +22,8 @@ import numpy as np
import matplotlib.pyplot as plt
import h5py
from tomorec.supp.suppTools import normaliser
-
+import time
+from libtiff import TIFF
# load dendritic projection data
h5f = h5py.File('data/DendrData_3D.h5','r')
@@ -36,7 +37,7 @@ h5f.close()
data_norm = normaliser(dataRaw, flats, darks, log='log')
del dataRaw, darks, flats
-intens_max = 2
+intens_max = 2.3
plt.figure()
plt.subplot(131)
plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max)
@@ -49,29 +50,38 @@ plt.imshow(data_norm[:,:,600],vmin=0, vmax=intens_max)
plt.title('Tangentogram view')
plt.show()
-
detectorHoriz = np.size(data_norm,2)
det_y_crop = [i for i in range(0,detectorHoriz-22)]
N_size = 950 # reconstruction domain
+time_label = int(time.time())
#%%
+"""
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%")
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
from tomorec.methodsDIR import RecToolsDIR
RectoolsDIR = RecToolsDIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal)
- DetectorsDimV = 10, # DetectorsDimV # detector dimension (vertical) for 3D case only
+ DetectorsDimV = 200, # 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(data_norm[0:10,:,det_y_crop])
+FBPrec = RectoolsDIR.FBP(data_norm[20:220,:,det_y_crop])
plt.figure()
-#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')
+FBPrec += np.abs(np.min(FBPrec))
+multiplier = (int)(65535/(np.max(FBPrec)))
+
+# saving to tiffs (16bit)
+for i in range(0,np.size(FBPrec,0)):
+ tiff = TIFF.open('Dendr_FBP'+'_'+str(i)+'.tiff', mode='w')
+ tiff.write_image(np.uint16(FBPrec[i,:,:]*multiplier))
+ tiff.close()
+"""
#%%
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("Reconstructing with ADMM method using TomoRec software")
@@ -79,7 +89,7 @@ print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
# initialise TomoRec ITERATIVE reconstruction class ONCE
from tomorec.methodsIR import RecToolsIR
RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal)
- DetectorsDimV = 5, # DetectorsDimV # detector dimension (vertical) for 3D case only
+ DetectorsDimV = 200, # 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 (wip), GH (wip), Student (wip)
@@ -88,29 +98,125 @@ RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH #
tolerance = 1e-08, # tolerance to stop outer iterations earlier
device='gpu')
#%%
-print ("Reconstructing with ADMM method using ROF-TV penalty")
+print ("Reconstructing with ADMM method using SB-TV penalty")
+RecADMM_reg_sbtv = RectoolsIR.ADMM(data_norm[20:220,:,det_y_crop],
+ rho_const = 2000.0, \
+ iterationsADMM = 15, \
+ regularisation = 'SB_TV', \
+ regularisation_parameter = 0.00085,\
+ regularisation_iterations = 50)
+
+sliceSel = 5
+max_val = 0.003
+plt.figure()
+plt.subplot(131)
+plt.imshow(RecADMM_reg_sbtv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray")
+plt.title('3D ADMM-SB-TV Reconstruction, axial view')
+
+plt.subplot(132)
+plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray")
+plt.title('3D ADMM-SB-TV Reconstruction, coronal view')
+
+plt.subplot(133)
+plt.imshow(RecADMM_reg_sbtv[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray")
+plt.title('3D ADMM-SB-TV Reconstruction, sagittal view')
+plt.show()
+
+multiplier = (int)(65535/(np.max(RecADMM_reg_sbtv)))
+
+# saving to tiffs (16bit)
+for i in range(0,np.size(RecADMM_reg_sbtv,0)):
+ tiff = TIFF.open('Dendr_ADMM_SBTV'+'_'+str(i)+'.tiff', mode='w')
+ tiff.write_image(np.uint16(RecADMM_reg_sbtv[i,:,:]*multiplier))
+ tiff.close()
-RecADMM_reg_roftv = RectoolsIR.ADMM(data_norm[0:5,:,det_y_crop],
+# Saving recpnstructed data with a unique time label
+np.save('Dendr_ADMM_SBTV'+str(time_label)+'.npy', RecADMM_reg_sbtv)
+del RecADMM_reg_sbtv
+#%%
+print ("Reconstructing with ADMM method using ROF-LLT penalty")
+RecADMM_reg_rofllt = RectoolsIR.ADMM(data_norm[20:220,:,det_y_crop],
rho_const = 2000.0, \
- iterationsADMM = 3, \
- regularisation = 'FGP_TV', \
- regularisation_parameter = 0.001,\
- regularisation_iterations = 80)
+ iterationsADMM = 15, \
+ regularisation = 'LLT_ROF', \
+ regularisation_parameter = 0.0009,\
+ regularisation_parameter2 = 0.0007,\
+ time_marching_parameter = 0.001,\
+ regularisation_iterations = 450)
+
+sliceSel = 5
+max_val = 0.003
+plt.figure()
+plt.subplot(131)
+plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val)
+plt.title('3D ADMM-ROFLLT Reconstruction, axial view')
+
+plt.subplot(132)
+plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val)
+plt.title('3D ADMM-ROFLLT Reconstruction, coronal view')
+
+plt.subplot(133)
+plt.imshow(RecADMM_reg_rofllt[:,:,sliceSel],vmin=0, vmax=max_val)
+plt.title('3D ADMM-ROFLLT Reconstruction, sagittal view')
+plt.show()
+
+multiplier = (int)(65535/(np.max(RecADMM_reg_rofllt)))
+
+# saving to tiffs (16bit)
+for i in range(0,np.size(RecADMM_reg_rofllt,0)):
+ tiff = TIFF.open('Dendr_ADMM_ROFLLT'+'_'+str(i)+'.tiff', mode='w')
+ tiff.write_image(np.uint16(RecADMM_reg_rofllt[i,:,:]*multiplier))
+ tiff.close()
+
+
+# Saving recpnstructed data with a unique time label
+np.save('Dendr_ADMM_ROFLLT'+str(time_label)+'.npy', RecADMM_reg_rofllt)
+del RecADMM_reg_rofllt
+#%%
+RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal)
+ 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
+ datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
+ nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
+ OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
+ tolerance = 1e-08, # tolerance to stop outer iterations earlier
+ device='cpu')
+print ("Reconstructing with ADMM method using TGV penalty")
+RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:10,:,det_y_crop],
+ rho_const = 2000.0, \
+ iterationsADMM = 15, \
+ regularisation = 'TGV', \
+ regularisation_parameter = 0.01,\
+ regularisation_iterations = 450)
-sliceSel = 2
-max_val = 0.005
+sliceSel = 7
+max_val = 0.003
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.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_roftv[:,sliceSel,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-ROF-TV Reconstruction, coronal view')
+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_roftv[:,:,sliceSel],vmin=0, vmax=max_val)
-plt.title('3D ADMM-ROF-TV Reconstruction, sagittal view')
+plt.imshow(RecADMM_reg_tgv[:,:,sliceSel],vmin=0, vmax=max_val)
+plt.title('3D ADMM-TGV Reconstruction, sagittal view')
plt.show()
+
+multiplier = (int)(65535/(np.max(RecADMM_reg_tgv)))
+
+# saving to tiffs (16bit)
+for i in range(0,np.size(RecADMM_reg_tgv,0)):
+ tiff = TIFF.open('Dendr_ADMM_TGV'+'_'+str(i)+'.tiff', mode='w')
+ tiff.write_image(np.uint16(RecADMM_reg_tgv[i,:,:]*multiplier))
+ tiff.close()
+
+
+# Saving recpnstructed data with a unique time label
+#np.save('Dendr_ADMM_TGV'+str(time_label)+'.npy', RecADMM_reg_tgv)
+del RecADMM_reg_tgv
#%%
\ No newline at end of file
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 5dbd436..a022ad7 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -197,32 +197,3 @@ 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.0005,\
- regularisation_iterations = 200)
-
-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
From 36cd88670b192e93f611863d58a438b9135b097c Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Mon, 25 Feb 2019 17:24:54 +0000
Subject: optim demo upd
---
.../SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 18 ++++++++++++++++++
1 file changed, 18 insertions(+)
(limited to 'Wrappers')
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 a79d0a3..c4f33ba 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
@@ -98,6 +98,12 @@ for i in range(0,param_space):
plt.figure()
plt.plot(erros_vec_sbtv)
+
+# Saving generated data with a unique time label
+h5f = h5py.File('Optim_admm_sbtv.h5', 'w')
+h5f.create_dataset('reg_param_sb_vec', data=reg_param_sb_vec)
+h5f.create_dataset('erros_vec_sbtv', data=erros_vec_sbtv)
+h5f.close()
#%%
param_space = 30
reg_param_rofllt_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters
@@ -119,6 +125,12 @@ for i in range(0,param_space):
plt.figure()
plt.plot(erros_vec_rofllt)
+
+# Saving generated data with a unique time label
+h5f = h5py.File('Optim_admm_rofllt.h5', 'w')
+h5f.create_dataset('reg_param_rofllt_vec', data=reg_param_rofllt_vec)
+h5f.create_dataset('erros_vec_rofllt', data=erros_vec_rofllt)
+h5f.close()
#%%
param_space = 30
reg_param_tgv_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters
@@ -139,4 +151,10 @@ for i in range(0,param_space):
plt.figure()
plt.plot(erros_vec_tgv)
+
+# Saving generated data with a unique time label
+h5f = h5py.File('Optim_admm_tgv.h5', 'w')
+h5f.create_dataset('reg_param_tgv_vec', data=reg_param_tgv_vec)
+h5f.create_dataset('erros_vec_tgv', data=erros_vec_tgv)
+h5f.close()
#%%
\ No newline at end of file
--
cgit v1.2.3
From 170e4be433e3a2b77e9fcfda47eadf59b4f9f6d4 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Tue, 26 Feb 2019 10:02:25 +0000
Subject: optim parameters added
---
.../SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 | Bin 0 -> 2408 bytes
.../demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 | Bin 0 -> 2408 bytes
.../demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 | Bin 0 -> 2408 bytes
3 files changed, 0 insertions(+), 0 deletions(-)
create mode 100644 Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5
create mode 100644 Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5
create mode 100644 Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5
(limited to 'Wrappers')
diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5
new file mode 100644
index 0000000..63bc4fd
Binary files /dev/null and b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 differ
diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5
new file mode 100644
index 0000000..03c0c14
Binary files /dev/null and b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 differ
diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5
new file mode 100644
index 0000000..056d915
Binary files /dev/null and b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 differ
--
cgit v1.2.3
From e05bcdce76e4e20154228e6ccf56207412340c96 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Tue, 26 Feb 2019 15:16:41 +0000
Subject: optim demo upd, tgv 3d fix
---
.../demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 2 +-
.../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 259 +++++++++++++--------
2 files changed, 168 insertions(+), 93 deletions(-)
(limited to 'Wrappers')
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 4cd680e..dc4eff2 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
@@ -181,7 +181,7 @@ RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH #
nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
tolerance = 1e-08, # tolerance to stop outer iterations earlier
- device='cpu')
+ device='gpu')
print ("Reconstructing with ADMM method using TGV penalty")
RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:10,:,det_y_crop],
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 a022ad7..cf9e797 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -19,6 +19,7 @@ GPLv3 license (ASTRA toolbox)
"""
#import timeit
import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
import numpy as np
import h5py
from ccpi.supp.qualitymetrics import QualityTools
@@ -33,34 +34,73 @@ h5f.close()
[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm)
N_size = Vert_det
-sliceSel = 128
-#plt.gray()
-plt.figure()
-plt.subplot(131)
-plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1)
-plt.title('3D Phantom, axial view')
+# loading optmisation parameters (the result of running Demo_SimulData_ParOptimis_SX)
+h5f = h5py.File('optim_param/Optim_admm_sbtv.h5','r')
+reg_param_sb_vec = h5f['reg_param_sb_vec'][:]
+erros_vec_sbtv = h5f['erros_vec_sbtv'][:]
+h5f.close()
-plt.subplot(132)
-plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1)
-plt.title('3D Phantom, coronal view')
+h5f = h5py.File('optim_param/Optim_admm_rofllt.h5','r')
+reg_param_rofllt_vec = h5f['reg_param_rofllt_vec'][:]
+erros_vec_rofllt = h5f['erros_vec_rofllt'][:]
+h5f.close()
+
+h5f = h5py.File('optim_param/Optim_admm_tgv.h5','r')
+reg_param_tgv_vec = h5f['reg_param_tgv_vec'][:]
+erros_vec_tgv = h5f['erros_vec_tgv'][:]
+h5f.close()
+index_minSBTV = min(xrange(len(erros_vec_sbtv)), key=erros_vec_sbtv.__getitem__)
+index_minROFLLT = min(xrange(len(erros_vec_rofllt)), key=erros_vec_rofllt.__getitem__)
+index_minTGV = min(xrange(len(erros_vec_tgv)), key=erros_vec_tgv.__getitem__)
+# assign optimal regularisation parameters:
+optimReg_sbtv = reg_param_sb_vec[index_minSBTV]
+optimReg_rofllt = reg_param_rofllt_vec[index_minROFLLT]
+optimReg_tgv = reg_param_tgv_vec[index_minTGV]
+#%%
+# plot loaded data
+sliceSel = 128
+#plt.figure()
+fig, (ax1, ax2) = plt.subplots(figsize=(15, 5), ncols=2)
+plt.rcParams.update({'xtick.labelsize': 'x-small'})
+plt.rcParams.update({'ytick.labelsize':'x-small'})
+plt.subplot(121)
+one = plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1, interpolation='none', cmap="PuOr")
+fig.colorbar(one, ax=ax1)
+plt.title('3D Phantom, axial (X-Y) view')
+plt.subplot(122)
+two = plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1,interpolation='none', cmap="PuOr")
+fig.colorbar(two, ax=ax2)
+plt.title('3D Phantom, coronal (Y-Z) view')
+"""
plt.subplot(133)
-plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1)
+plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1, cmap="PuOr")
plt.title('3D Phantom, sagittal view')
-plt.show()
-intens_max = 240
+"""
+plt.show()
+#%%
+intens_max = 220
plt.figure()
+plt.rcParams.update({'xtick.labelsize': 'x-small'})
+plt.rcParams.update({'ytick.labelsize':'x-small'})
plt.subplot(131)
-plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max)
-plt.title('2D Projection (erroneous)')
+plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max, cmap="PuOr")
+plt.xlabel('X-detector', fontsize=16)
+plt.ylabel('Z-detector', fontsize=16)
+plt.title('2D Projection (X-Z) view', fontsize=19)
plt.subplot(132)
-plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max)
-plt.title('Sinogram view')
+plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max, cmap="PuOr")
+plt.xlabel('X-detector', fontsize=16)
+plt.ylabel('Projection angle', fontsize=16)
+plt.title('Sinogram (X-Y) view', fontsize=19)
plt.subplot(133)
-plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max)
-plt.title('Tangentogram view')
+plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max, cmap="PuOr")
+plt.xlabel('Projection angle', fontsize=16)
+plt.ylabel('Z-detector', fontsize=16)
+plt.title('Vertical (Y-Z) view', fontsize=19)
plt.show()
+#plt.savefig('projdata.pdf', format='pdf', dpi=1200)
#%%
# initialise TomoRec DIRECT reconstruction class ONCE
from tomorec.methodsDIR import RecToolsDIR
@@ -72,23 +112,31 @@ RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector
#%%
print ("Reconstruction using FBP from TomoRec")
recFBP= RectoolsDIR.FBP(projdata_norm) # FBP reconstruction
+#%%
+x0, y0 = 0, 127 # These are in _pixel_ coordinates!!
+x1, y1 = 255, 127
sliceSel = int(0.5*N_size)
max_val = 1
-#plt.gray()
-plt.figure()
-plt.subplot(131)
-plt.imshow(recFBP[sliceSel,:,:],vmin=0, vmax=max_val)
-plt.title('3D Reconstruction, axial view')
-
-plt.subplot(132)
-plt.imshow(recFBP[:,sliceSel,:],vmin=0, vmax=max_val)
-plt.title('3D Reconstruction, coronal view')
-
-plt.subplot(133)
-plt.imshow(recFBP[:,:,sliceSel],vmin=0, vmax=max_val)
-plt.title('3D Reconstruction, sagittal view')
+plt.figure(figsize = (20,5))
+gs1 = gridspec.GridSpec(1, 3)
+gs1.update(wspace=0.1, hspace=0.05) # set the spacing between axes.
+ax1 = plt.subplot(gs1[0])
+plt.imshow(recFBP[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
+ax1.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
+plt.colorbar(ax=ax1)
+plt.title('FBP Reconstruction, axial (X-Y) view', fontsize=19)
+ax1.set_aspect('equal')
+ax3 = plt.subplot(gs1[1])
+plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
+plt.plot(recFBP[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
+plt.title('Profile', fontsize=19)
+ax2 = plt.subplot(gs1[2])
+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)
# calculate errors
Qtools = QualityTools(phantom, recFBP)
@@ -110,87 +158,114 @@ RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector d
tolerance = 1e-08, # tolerance to stop outer iterations earlier
device='gpu')
#%%
-print ("Reconstructing with ADMM method using ROF-TV penalty")
-RecADMM_reg_roftv = RectoolsIR.ADMM(projdata_norm,
- rho_const = 2000.0, \
- iterationsADMM = 30, \
- regularisation = 'ROF_TV', \
- regularisation_parameter = 0.003,\
- regularisation_iterations = 500)
+print ("Reconstructing with ADMM method using SB-TV penalty")
+RecADMM_reg_sbtv = RectoolsIR.ADMM(projdata_norm,
+ rho_const = 2000.0, \
+ iterationsADMM = 25, \
+ regularisation = 'SB_TV', \
+ regularisation_parameter = optimReg_sbtv,\
+ regularisation_iterations = 50)
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.figure(figsize = (20,3))
+gs1 = gridspec.GridSpec(1, 4)
+gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes.
+ax1 = plt.subplot(gs1[0])
+plt.plot(reg_param_sb_vec, erros_vec_sbtv, color='k',linewidth=2)
+plt.xlabel('Regularisation parameter', fontsize=16)
+plt.ylabel('RMSE value', fontsize=16)
+plt.title('Regularisation selection', fontsize=19)
+ax2 = plt.subplot(gs1[1])
+plt.imshow(RecADMM_reg_sbtv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
+ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
+plt.title('ADMM-SBTV (X-Y) view', fontsize=19)
+#ax2.set_aspect('equal')
+ax3 = plt.subplot(gs1[2])
+plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
+plt.plot(RecADMM_reg_sbtv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
+plt.title('Profile', fontsize=19)
+ax4 = plt.subplot(gs1[3])
+plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
+plt.title('ADMM-SBTV (Y-Z) view', fontsize=19)
+plt.colorbar(ax=ax4)
plt.show()
+plt.savefig('SBTV_phantom.pdf', format='pdf', dpi=1600)
# 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))
+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))
#%%
-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)
+print ("Reconstructing with ADMM method using ROFLLT penalty")
+RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm,
+ rho_const = 2000.0, \
+ iterationsADMM = 25, \
+ regularisation = 'LLT_ROF', \
+ regularisation_parameter = optimReg_rofllt,\
+ regularisation_parameter2 = 0.0085,\
+ regularisation_iterations = 600)
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.figure(figsize = (20,3))
+gs1 = gridspec.GridSpec(1, 4)
+gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes.
+ax1 = plt.subplot(gs1[0])
+plt.plot(reg_param_rofllt_vec, erros_vec_rofllt, color='k',linewidth=2)
+plt.xlabel('Regularisation parameter', fontsize=16)
+plt.ylabel('RMSE value', fontsize=16)
+plt.title('Regularisation selection', fontsize=19)
+ax2 = plt.subplot(gs1[1])
+plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
+ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
+plt.title('ADMM-ROFLLT (X-Y) view', fontsize=19)
+#ax2.set_aspect('equal')
+ax3 = plt.subplot(gs1[2])
+plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
+plt.plot(RecADMM_reg_rofllt[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
+plt.title('Profile', fontsize=19)
+ax4 = plt.subplot(gs1[3])
+plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
+plt.title('ADMM-ROFLLT (Y-Z) view', fontsize=19)
+plt.colorbar(ax=ax4)
plt.show()
+#plt.savefig('ROFLLT_phantom.pdf', format='pdf', dpi=1600)
# calculate errors
-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))
+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))
#%%
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)
+ rho_const = 2000.0, \
+ iterationsADMM = 25, \
+ regularisation = 'TGV', \
+ regularisation_parameter = optimReg_tgv,\
+ regularisation_iterations = 600)
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.figure(figsize = (20,5))
+gs1 = gridspec.GridSpec(1, 3)
+gs1.update(wspace=0.1, hspace=0.05) # set the spacing between axes.
+ax1 = plt.subplot(gs1[0])
+plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
+ax1.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
+plt.colorbar(ax=ax1)
+plt.title('ADMM-TGV reconstruction, (X-Y) view', fontsize=19)
+ax1.set_aspect('equal')
+ax3 = plt.subplot(gs1[1])
+plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
+plt.plot(RecADMM_reg_tgv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
+plt.title('Profile', fontsize=19)
+ax2 = plt.subplot(gs1[2])
+plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
+plt.title('ADMM-TGV reconstruction, (Y-Z) view', fontsize=19)
+ax2.set_aspect('equal')
plt.show()
+plt.savefig('TGV_phantom.pdf', format='pdf', dpi=1200)
# calculate errors
Qtools = QualityTools(phantom, RecADMM_reg_tgv)
--
cgit v1.2.3
From 57e60c85cb56de6c77a56737af5e6ad57d801420 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Tue, 26 Feb 2019 15:29:08 +0000
Subject: tgv_upd
---
.../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 26 +++++++++++++---------
1 file changed, 15 insertions(+), 11 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 cf9e797..d39d643 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -247,23 +247,27 @@ RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm,
sliceSel = int(0.5*N_size)
max_val = 1
-plt.figure(figsize = (20,5))
-gs1 = gridspec.GridSpec(1, 3)
-gs1.update(wspace=0.1, hspace=0.05) # set the spacing between axes.
+plt.figure(figsize = (20,3))
+gs1 = gridspec.GridSpec(1, 4)
+gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes.
ax1 = plt.subplot(gs1[0])
+plt.plot(reg_param_tgv_vec, erros_vec_tgv, color='k',linewidth=2)
+plt.xlabel('Regularisation parameter', fontsize=16)
+plt.ylabel('RMSE value', fontsize=16)
+plt.title('Regularisation selection', fontsize=19)
+ax2 = plt.subplot(gs1[1])
plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
-ax1.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
-plt.colorbar(ax=ax1)
-plt.title('ADMM-TGV reconstruction, (X-Y) view', fontsize=19)
-ax1.set_aspect('equal')
-ax3 = plt.subplot(gs1[1])
+ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
+plt.title('ADMM-TGV (X-Y) view', fontsize=19)
+#ax2.set_aspect('equal')
+ax3 = plt.subplot(gs1[2])
plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
plt.plot(RecADMM_reg_tgv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
plt.title('Profile', fontsize=19)
-ax2 = plt.subplot(gs1[2])
+ax4 = plt.subplot(gs1[3])
plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
-plt.title('ADMM-TGV reconstruction, (Y-Z) view', fontsize=19)
-ax2.set_aspect('equal')
+plt.title('ADMM-TGV (Y-Z) view', fontsize=19)
+plt.colorbar(ax=ax4)
plt.show()
plt.savefig('TGV_phantom.pdf', format='pdf', dpi=1200)
--
cgit v1.2.3
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(-)
(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 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
From 2930421c6fe94faa81a3871ceebbc2cf154d2b97 Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Wed, 27 Feb 2019 11:39:05 +0000
Subject: demo files updated
---
.../demos/SoftwareX_supp/Demo_RealData_Recon_SX.py | 85 ++++++++++++----------
.../SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py | 5 +-
.../SoftwareX_supp/Demo_SimulData_Recon_SX.py | 7 +-
.../demos/SoftwareX_supp/Demo_SimulData_SX.py | 6 +-
Wrappers/Python/demos/SoftwareX_supp/Readme.md | 17 +++--
5 files changed, 69 insertions(+), 51 deletions(-)
(limited to 'Wrappers')
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 dc4eff2..01491d9 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
@@ -3,10 +3,11 @@
"""
This demo scripts support the following publication:
"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
-proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham,
-Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
+ Philip J. Withers; Software X, 2019
____________________________________________________________________________
* Reads real tomographic data (stored at Zenodo)
+--- https://doi.org/10.5281/zenodo.2578893
* Reconstructs using TomoRec software
* Saves reconstructed images
____________________________________________________________________________
@@ -14,6 +15,8 @@ ____________________________________________________________________________
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. libtiff if one needs to save tiff images:
+ install pip install libtiff
@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
GPLv3 license (ASTRA toolbox)
@@ -23,7 +26,6 @@ import matplotlib.pyplot as plt
import h5py
from tomorec.supp.suppTools import normaliser
import time
-from libtiff import TIFF
# load dendritic projection data
h5f = h5py.File('data/DendrData_3D.h5','r')
@@ -55,24 +57,38 @@ det_y_crop = [i for i in range(0,detectorHoriz-22)]
N_size = 950 # reconstruction domain
time_label = int(time.time())
#%%
-"""
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%")
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
from tomorec.methodsDIR import RecToolsDIR
RectoolsDIR = RecToolsDIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal)
- DetectorsDimV = 200, # DetectorsDimV # detector dimension (vertical) for 3D case only
+ DetectorsDimV = 100, # 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(data_norm[20:220,:,det_y_crop])
+FBPrec = RectoolsDIR.FBP(data_norm[0:100,:,det_y_crop])
-plt.figure()
-plt.imshow(FBPrec[0,:,:], vmin=0, vmax=0.005, cmap="gray")
-plt.title('FBP reconstruction')
+sliceSel = 50
+max_val = 0.003
+plt.figure()
+plt.subplot(131)
+plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray")
+plt.title('FBP Reconstruction, axial view')
+plt.subplot(132)
+plt.imshow(FBPrec[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray")
+plt.title('FBP Reconstruction, coronal view')
+
+plt.subplot(133)
+plt.imshow(FBPrec[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray")
+plt.title('FBP Reconstruction, sagittal view')
+plt.show()
+
+# saving to tiffs (16bit)
+"""
+from libtiff import TIFF
FBPrec += np.abs(np.min(FBPrec))
multiplier = (int)(65535/(np.max(FBPrec)))
@@ -89,7 +105,7 @@ print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
# initialise TomoRec ITERATIVE reconstruction class ONCE
from tomorec.methodsIR import RecToolsIR
RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal)
- DetectorsDimV = 200, # DetectorsDimV # detector dimension (vertical) for 3D case only
+ DetectorsDimV = 100, # 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 (wip), GH (wip), Student (wip)
@@ -99,14 +115,14 @@ RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH #
device='gpu')
#%%
print ("Reconstructing with ADMM method using SB-TV penalty")
-RecADMM_reg_sbtv = RectoolsIR.ADMM(data_norm[20:220,:,det_y_crop],
+RecADMM_reg_sbtv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
rho_const = 2000.0, \
iterationsADMM = 15, \
regularisation = 'SB_TV', \
regularisation_parameter = 0.00085,\
regularisation_iterations = 50)
-sliceSel = 5
+sliceSel = 50
max_val = 0.003
plt.figure()
plt.subplot(131)
@@ -122,29 +138,31 @@ plt.imshow(RecADMM_reg_sbtv[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray")
plt.title('3D ADMM-SB-TV Reconstruction, sagittal view')
plt.show()
-multiplier = (int)(65535/(np.max(RecADMM_reg_sbtv)))
# saving to tiffs (16bit)
+"""
+from libtiff import TIFF
+multiplier = (int)(65535/(np.max(RecADMM_reg_sbtv)))
for i in range(0,np.size(RecADMM_reg_sbtv,0)):
tiff = TIFF.open('Dendr_ADMM_SBTV'+'_'+str(i)+'.tiff', mode='w')
tiff.write_image(np.uint16(RecADMM_reg_sbtv[i,:,:]*multiplier))
tiff.close()
-
+"""
# Saving recpnstructed data with a unique time label
np.save('Dendr_ADMM_SBTV'+str(time_label)+'.npy', RecADMM_reg_sbtv)
del RecADMM_reg_sbtv
#%%
print ("Reconstructing with ADMM method using ROF-LLT penalty")
-RecADMM_reg_rofllt = RectoolsIR.ADMM(data_norm[20:220,:,det_y_crop],
+RecADMM_reg_rofllt = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
rho_const = 2000.0, \
iterationsADMM = 15, \
regularisation = 'LLT_ROF', \
regularisation_parameter = 0.0009,\
regularisation_parameter2 = 0.0007,\
time_marching_parameter = 0.001,\
- regularisation_iterations = 450)
+ regularisation_iterations = 550)
-sliceSel = 5
+sliceSel = 50
max_val = 0.003
plt.figure()
plt.subplot(131)
@@ -160,38 +178,29 @@ plt.imshow(RecADMM_reg_rofllt[:,:,sliceSel],vmin=0, vmax=max_val)
plt.title('3D ADMM-ROFLLT Reconstruction, sagittal view')
plt.show()
-multiplier = (int)(65535/(np.max(RecADMM_reg_rofllt)))
-
# saving to tiffs (16bit)
+"""
+from libtiff import TIFF
+multiplier = (int)(65535/(np.max(RecADMM_reg_rofllt)))
for i in range(0,np.size(RecADMM_reg_rofllt,0)):
tiff = TIFF.open('Dendr_ADMM_ROFLLT'+'_'+str(i)+'.tiff', mode='w')
tiff.write_image(np.uint16(RecADMM_reg_rofllt[i,:,:]*multiplier))
tiff.close()
-
+"""
# Saving recpnstructed data with a unique time label
np.save('Dendr_ADMM_ROFLLT'+str(time_label)+'.npy', RecADMM_reg_rofllt)
del RecADMM_reg_rofllt
#%%
-RectoolsIR = RecToolsIR(DetectorsDimH = np.size(det_y_crop), # DetectorsDimH # detector dimension (horizontal)
- 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
- datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
- nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
- OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
- tolerance = 1e-08, # tolerance to stop outer iterations earlier
- device='gpu')
-
print ("Reconstructing with ADMM method using TGV penalty")
-RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:10,:,det_y_crop],
+RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
rho_const = 2000.0, \
iterationsADMM = 15, \
regularisation = 'TGV', \
regularisation_parameter = 0.01,\
- regularisation_iterations = 450)
+ regularisation_iterations = 500)
-sliceSel = 7
+sliceSel = 50
max_val = 0.003
plt.figure()
plt.subplot(131)
@@ -207,16 +216,16 @@ plt.imshow(RecADMM_reg_tgv[:,:,sliceSel],vmin=0, vmax=max_val)
plt.title('3D ADMM-TGV Reconstruction, sagittal view')
plt.show()
-multiplier = (int)(65535/(np.max(RecADMM_reg_tgv)))
-
# saving to tiffs (16bit)
+"""
+from libtiff import TIFF
+multiplier = (int)(65535/(np.max(RecADMM_reg_tgv)))
for i in range(0,np.size(RecADMM_reg_tgv,0)):
tiff = TIFF.open('Dendr_ADMM_TGV'+'_'+str(i)+'.tiff', mode='w')
tiff.write_image(np.uint16(RecADMM_reg_tgv[i,:,:]*multiplier))
tiff.close()
-
-
+"""
# Saving recpnstructed data with a unique time label
-#np.save('Dendr_ADMM_TGV'+str(time_label)+'.npy', RecADMM_reg_tgv)
+np.save('Dendr_ADMM_TGV'+str(time_label)+'.npy', RecADMM_reg_tgv)
del RecADMM_reg_tgv
#%%
\ No newline at end of file
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 c4f33ba..59ffc0e 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
@@ -3,10 +3,11 @@
"""
This demo scripts support the following publication:
"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
-proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham,
-Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
+ Philip J. Withers; Software X, 2019
____________________________________________________________________________
* Reads data which is previosly generated by TomoPhantom software (Zenodo link)
+--- https://doi.org/10.5281/zenodo.2578893
* Optimises for the regularisation parameters which later used in the script:
Demo_SimulData_Recon_SX.py
____________________________________________________________________________
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 aa65cf3..93b0cef 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
@@ -3,10 +3,11 @@
"""
This demo scripts support the following publication:
"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
-proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham,
-Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
+ Philip J. Withers; Software X, 2019
____________________________________________________________________________
* Reads data which is previously generated by TomoPhantom software (Zenodo link)
+--- https://doi.org/10.5281/zenodo.2578893
* Reconstruct using optimised regularisation parameters (see Demo_SimulData_ParOptimis_SX.py)
____________________________________________________________________________
>>>>> Dependencies: <<<<<
@@ -305,4 +306,4 @@ 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]))
-#%%
+#%%
\ No newline at end of file
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
index ce29f0c..cdf4325 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
+++ b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
@@ -3,13 +3,13 @@
"""
This demo scripts support the following publication:
"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with
-proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham,
-Martin J. Turner, Philip J. Withers and Alun Ashton; Software X, 2019
+proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
+ Philip J. Withers; Software X, 2019
____________________________________________________________________________
* Runs TomoPhantom software to simulate tomographic projection data with
some imaging errors and noise
* Saves the data into hdf file to be uploaded in reconstruction scripts
-____________________________________________________________________________
+__________________________________________________________________________
>>>>> Dependencies: <<<<<
1. TomoPhantom software for phantom and data generation
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Readme.md b/Wrappers/Python/demos/SoftwareX_supp/Readme.md
index fa77745..54e83f1 100644
--- a/Wrappers/Python/demos/SoftwareX_supp/Readme.md
+++ b/Wrappers/Python/demos/SoftwareX_supp/Readme.md
@@ -1,19 +1,26 @@
+
# SoftwareX publication [1] supporting files
## Decription:
-The scripts here support publication in SoftwareX journal [1] to ensure
-reproducibility of the research.
+The scripts here support publication in SoftwareX journal [1] to ensure reproducibility of the research. The scripts linked with data shared at Zenodo.
## Data:
+Data is shared at Zenodo [here](https://doi.org/10.5281/zenodo.2578893)
## Dependencies:
+1. [ASTRA toolbox](https://github.com/astra-toolbox/astra-toolbox): `conda install -c astra-toolbox astra-toolbox`
+2. [TomoRec](https://github.com/dkazanc/TomoRec): `conda install -c dkazanc tomorec`
+3. [Tomophantom](https://github.com/dkazanc/TomoPhantom): `conda install tomophantom -c ccpi`
## Files description:
-
+- `Demo_SimulData_SX.py` - simulates 3D projection data using [Tomophantom](https://github.com/dkazanc/TomoPhantom) software. One can skip this module if the data is taken from [Zenodo](https://doi.org/10.5281/zenodo.2578893)
+- `Demo_SimulData_ParOptimis_SX.py` - runs computationally extensive calculations for optimal regularisation parameters, the result are saved into directory `optim_param`. This script can be also skipped.
+- `Demo_SimulData_Recon_SX.py` - using established regularisation parameters, one runs iterative reconstruction
+- `Demo_RealData_Recon_SX.py` - runs real data reconstructions. Can be quite intense on memory so reduce the size of the reconstructed volume if needed.
### References:
-[1] "CCPi-Regularisation Toolkit for computed tomographic image reconstruction with proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Mark Basham, Martin J. Turner, Philip J. Withers and Alun Ashton
-
+[1] "CCPi-Regularisation Toolkit for computed tomographic image reconstruction with proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner and Philip J. Withers; SoftwareX, 2019.
### Acknowledgments:
CCPi-RGL software is a product of the [CCPi](https://www.ccpi.ac.uk/) group, STFC SCD software developers and Diamond Light Source (DLS). Any relevant questions/comments can be e-mailed to Daniil Kazantsev at dkazanc@hotmail.com
+
--
cgit v1.2.3