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(-) 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