From 96c5d86fe9a9ac9725ae64d1c50ec9af71137608 Mon Sep 17 00:00:00 2001 From: dkazanc 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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 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/Python') 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