From df5b9b9f93d0dd666d571be6ce2f7afd864fbbf4 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 31 Jul 2018 15:12:25 +0100 Subject: Added working demo of colourbay data load and recon. --- Wrappers/Python/wip/demo_colourbay.py | 137 ++++++++++++++++++++++++++++++++++ 1 file changed, 137 insertions(+) create mode 100644 Wrappers/Python/wip/demo_colourbay.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/demo_colourbay.py b/Wrappers/Python/wip/demo_colourbay.py new file mode 100644 index 0000000..5dbf2e1 --- /dev/null +++ b/Wrappers/Python/wip/demo_colourbay.py @@ -0,0 +1,137 @@ +# This script demonstrates how to load a mat-file with UoM colour-bay data +# into the CIL optimisation framework and run (simple) multichannel +# reconstruction methods. + +# All third-party imports. +import numpy +from scipy.io import loadmat +import matplotlib.pyplot as plt + +# All own imports. +from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.astra.ops import AstraProjectorMC +from ccpi.optimisation.algs import CGLS, FISTA +from ccpi.optimisation.funcs import Norm2sq, Norm1 + +# Load full data and permute to expected ordering. Change path as necessary. +# The loaded X has dims 80x60x80x150, which is pix x angle x pix x channel. +# Permute (numpy.transpose) puts into our default ordering which is +# (channel, angle, vertical, horizontal). + +pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/ColourBay/spectral_data_sets/CarbonPd/' +filename = 'carbonPd_full_sinogram_stripes_removed.mat' + +X = loadmat(pathname + filename) +X = numpy.transpose(X['SS'],(3,1,2,0)) + +# Store geometric variables for reuse +num_channels = X.shape[0] +num_pixels_h = X.shape[3] +num_pixels_v = X.shape[2] +num_angles = X.shape[1] + +# Display a single projection in a single channel +plt.imshow(X[100,5,:,:]) +plt.title('Example of a projection image in one channel' ) +plt.show() + +# Set angles to use +angles = numpy.linspace(-numpy.pi,numpy.pi,num_angles,endpoint=False) + +# Define full 3D acquisition geometry and data container. +# Geometric info is taken from the txt-file in the same dir as the mat-file +ag = AcquisitionGeometry('cone', + '3D', + angles, + pixel_num_h=num_pixels_h, + pixel_size_h=0.25, + pixel_num_v=num_pixels_v, + pixel_size_v=0.25, + dist_source_center=233.0, + dist_center_detector=245.0, + channels=num_channels) +data = AcquisitionData(X, geometry=ag) + +# Reduce to central slice by extracting relevant parameters from data and its +# geometry. Perhaps create function to extract central slice automatically? +data2d = data.subset(vertical=40) +ag2d = AcquisitionGeometry('cone', + '2D', + ag.angles, + pixel_num_h=ag.pixel_num_h, + pixel_size_h=ag.pixel_size_h, + pixel_num_v=1, + pixel_size_v=ag.pixel_size_h, + dist_source_center=ag.dist_source_center, + dist_center_detector=ag.dist_center_detector, + channels=ag.channels) +data2d.geometry = ag2d + +# Set up 2D Image Geometry. +# First need the geometric magnification to scale the voxel size relative +# to the detector pixel size. +mag = (ag.dist_source_center + ag.dist_center_detector)/ag.dist_source_center +ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, + voxel_num_y=ag2d.pixel_num_h, + voxel_size_x=ag2d.pixel_size_h/mag, + voxel_size_y=ag2d.pixel_size_h/mag, + channels=X.shape[0]) + +# Create GPU multichannel projector/backprojector operator with ASTRA. +Aall = AstraProjectorMC(ig2d,ag2d,'gpu') + +# Compute and simple backprojction and display one channel as image. +Xbp = Aall.adjoint(data2d) +plt.imshow(Xbp.subset(channel=100).array) +plt.show() + +# Set initial guess ImageData with zeros for algorithms, and algorithm options. +x_init = ImageData(numpy.zeros((num_channels,num_pixels_v,num_pixels_h)), + geometry=ig2d, + dimension_labels=['channel','horizontal_y','horizontal_x']) +opt_CGLS = {'tol': 1e-4, 'iter': 5} + +# Run CGLS algorithm and display one channel. +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aall, data2d, opt_CGLS) + +plt.imshow(x_CGLS.subset(channel=100).array) +plt.title('CGLS') +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS Criterion vs iterations') +plt.show() + +# Create least squares object instance with projector, test data and a constant +# coefficient of 0.5. Note it is least squares over all channels. +f = Norm2sq(Aall,data2d,c=0.5) + +# Options for FISTA algorithm. +opt = {'tol': 1e-4, 'iter': 100} + +# Run FISTA for least squares without regularization and display one channel +# reconstruction as image. +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) + +plt.imshow(x_fista0.subset(channel=100).array) +plt.title('FISTA LS') +plt.show() + +plt.semilogy(criter0) +plt.title('FISTA LS Criterion vs iterations') +plt.show() + +# Set up 1-norm regularisation (over all channels), solve with FISTA, and +# display one channel of reconstruction. +lam = 0.1 +g0 = Norm1(lam) + +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) + +plt.imshow(x_fista1.subset(channel=100).array) +plt.title('FISTA LS+1') +plt.show() + +plt.semilogy(criter1) +plt.title('FISTA LS+1 Criterion vs iterations') +plt.show() \ No newline at end of file -- cgit v1.2.3 From 2341ba146f16e79b097f55b1efa281bc06eef472 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Tue, 31 Jul 2018 22:37:51 +0100 Subject: Added IMAT white-beam demo loading summed fits files --- Wrappers/Python/wip/demo_imat_whitebeam.py | 128 +++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 Wrappers/Python/wip/demo_imat_whitebeam.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/demo_imat_whitebeam.py b/Wrappers/Python/wip/demo_imat_whitebeam.py new file mode 100644 index 0000000..af3d568 --- /dev/null +++ b/Wrappers/Python/wip/demo_imat_whitebeam.py @@ -0,0 +1,128 @@ +# This script demonstrates how to load IMAT fits data +# into the CIL optimisation framework and run reconstruction methods. +# +# This demo loads the summedImg files which are the non-spectral images +# resulting from summing projections over all spectral channels. + +# needs dxchange: conda install -c conda-forge dxchange +# needs astropy: conda install astropy + + +# All third-party imports. +import numpy +from scipy.io import loadmat +import matplotlib.pyplot as plt +from dxchange.reader import read_fits + +# All own imports. +from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple +from ccpi.optimisation.algs import CGLS, FISTA +from ccpi.optimisation.funcs import Norm2sq, Norm1 + +# Load and display a couple of summed projection as examples +pathname0 = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/' +filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits' + +data0 = read_fits(pathname0 + filename0) + +pathname10 = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/' +filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits' + +data10 = read_fits(pathname10 + filename10) + +# Load a flat field (more are available, should we average over them?) +flat1 = read_fits('/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') + +# Apply flat field and display after flat-field correction and negative log +data0_rel = data0 / flat1 +data10_rel = data10 / flat1 + +plt.imshow(data0_rel) +plt.colorbar() +plt.show() + +plt.imshow(-numpy.log(data0_rel)) +plt.colorbar() +plt.show() + +plt.imshow(data10_rel) +plt.colorbar() +plt.show() + +plt.imshow(-numpy.log(data10_rel)) +plt.colorbar() +plt.show() + +# Set up for loading all summed images at 250 angles. +pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/' +filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits' + +# Dimensions +num_angles = 250 +imsize = 512 + +# Initialise array +data = numpy.zeros((num_angles,imsize,imsize)) + +# Load only 0-249, as 250 is at repetition of zero degrees just like image 0 +for i in range(0,250): + curimfile = (pathname + filename).format(i, i+4675) + data[i,:,:] = read_fits(curimfile) + +# Apply flat field and take negative log +data_rel = -numpy.log(data/flat1) + +# Permute order to get: angles, vertical, horizontal, as default in framework. +data_rel = numpy.transpose(data_rel,(0,2,1)) + +# Set angles to use +angles = numpy.linspace(-numpy.pi,numpy.pi,num_angles,endpoint=False) + +# Create 3D acquisition geometry and acquisition data +ag = AcquisitionGeometry('parallel', + '3D', + angles, + pixel_num_h=imsize, + pixel_num_v=imsize) +b = AcquisitionData(data_rel, geometry=ag) + +# Reduce to single (noncentral) slice by extracting relevant parameters from data and its +# geometry. Perhaps create function to extract central slice automatically? +b2d = b.subset(vertical=128) +ag2d = AcquisitionGeometry('parallel', + '2D', + ag.angles, + pixel_num_h=ag.pixel_num_h) +b2d.geometry = ag2d + +# Create 2D image geometry +ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, + voxel_num_y=ag2d.pixel_num_h) + +# Create GPU projector/backprojector operator with ASTRA. +Aop = AstraProjectorSimple(ig2d,ag2d,'gpu') + +# Demonstrate operator is working by applying simple backprojection. +z = Aop.adjoint(b2d) +plt.imshow(z.array) +plt.title('Simple backprojection') +plt.colorbar() +plt.show() + +# Set initial guess ImageData with zeros for algorithms, and algorithm options. +x_init = ImageData(numpy.zeros((imsize,imsize)), + geometry=ig2d) +opt_CGLS = {'tol': 1e-4, 'iter': 20} + +# Run CGLS algorithm and display reconstruction. +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b2d, opt_CGLS) + +plt.imshow(x_CGLS.array) +plt.title('CGLS') +plt.colorbar() +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS Criterion vs iterations') +plt.show() \ No newline at end of file -- cgit v1.2.3 From ab323e0b38b941d2a7aa5b5e705518b770f1d36b Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 14 Aug 2018 12:45:29 +0100 Subject: corrections to normalization and log with zeroes in flats --- Wrappers/Python/wip/demo_imat_whitebeam.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/demo_imat_whitebeam.py b/Wrappers/Python/wip/demo_imat_whitebeam.py index af3d568..482c1ae 100644 --- a/Wrappers/Python/wip/demo_imat_whitebeam.py +++ b/Wrappers/Python/wip/demo_imat_whitebeam.py @@ -21,22 +21,25 @@ from ccpi.optimisation.algs import CGLS, FISTA from ccpi.optimisation.funcs import Norm2sq, Norm1 # Load and display a couple of summed projection as examples -pathname0 = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/' +pathname0 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle0/' filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits' data0 = read_fits(pathname0 + filename0) -pathname10 = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/' +pathname10 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle10/' filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits' data10 = read_fits(pathname10 + filename10) # Load a flat field (more are available, should we average over them?) -flat1 = read_fits('/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') +flat1 = read_fits('/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') # Apply flat field and display after flat-field correction and negative log -data0_rel = data0 / flat1 -data10_rel = data10 / flat1 +data0_rel = numpy.zeros(numpy.shape(flat1), dtype = float) +nonzero = flat1 > 0 +data0_rel[nonzero] = data0[nonzero] / flat1[nonzero] +data10_rel = numpy.zeros(numpy.shape(flat1), dtype = float) +data10_rel[nonzero] = data10[nonzero] / flat1[nonzero] plt.imshow(data0_rel) plt.colorbar() @@ -55,7 +58,7 @@ plt.colorbar() plt.show() # Set up for loading all summed images at 250 angles. -pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/' +pathname = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle{}/' filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits' # Dimensions @@ -71,7 +74,14 @@ for i in range(0,250): data[i,:,:] = read_fits(curimfile) # Apply flat field and take negative log -data_rel = -numpy.log(data/flat1) +nonzero = flat1 > 0 +for i in range(0,250): + data[i,nonzero] = data[i,nonzero]/flat1[nonzero] + +eqzero = data == 0 +data[eqzero] = 1 + +data_rel = -numpy.log(data) # Permute order to get: angles, vertical, horizontal, as default in framework. data_rel = numpy.transpose(data_rel,(0,2,1)) -- cgit v1.2.3 From c369a2950801fca4db606f67433b7bebc32fbdf1 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 14 Aug 2018 22:26:11 +0100 Subject: IMAT data multichannel script started --- Wrappers/Python/wip/demo_imat_multichan_RGLTK.py | 198 +++++++++++++++++++++++ 1 file changed, 198 insertions(+) create mode 100644 Wrappers/Python/wip/demo_imat_multichan_RGLTK.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py new file mode 100644 index 0000000..dae5f3a --- /dev/null +++ b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py @@ -0,0 +1,198 @@ +# This script demonstrates how to load IMAT fits data +# into the CIL optimisation framework and run reconstruction methods. +# +# This demo loads the summedImg files which are the non-spectral images +# resulting from summing projections over all spectral channels. + +# needs dxchange: conda install -c conda-forge dxchange +# needs astropy: conda install astropy + + +# All third-party imports. +import numpy as np +from scipy.io import loadmat +import matplotlib.pyplot as plt +from dxchange.reader import read_fits + +# All own imports. +from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple +from ccpi.optimisation.algs import CGLS, FISTA +from ccpi.optimisation.funcs import Norm2sq, Norm1 +from ccpi.plugins.regularisers import ROF_TV, FGP_TV, SB_TV + + +pathname0 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle0/' +filenameG = "IMAT00004675_Tomo_test_000_" + +n = 512 +totalAngles = 250 # total number of projection angles +# spectral discretization parameter +num_average = 120 +numChannels = 2970 +totChannels = round(numChannels/num_average) # total no. of averaged channels +Projections_stack = np.zeros((num_average,n,n),dtype='uint16') +ProjAngleChannels = np.zeros((totalAngles,totChannels,n,n),dtype='float32') + +counterT = 0 +for i in range(0,2,1): + for j in range(0,num_average,1): + if counterT < 10: + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathname0,filenameG,'0000',str(counterT)) + if ((counterT >= 10) & (counterT < 100)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathname0,filenameG,'000',str(counterT)) + if ((counterT >= 100) & (counterT < 1000)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathname0,filenameG,'00',str(counterT)) + if ((counterT >= 1000) & (counterT < 10000)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathname0,filenameG,'0',str(counterT)) + Projections_stack[j,:,:] = read_fits(outfile) + counterT = counterT + 1 + AverageProj=np.mean(Projections_stack,axis=0) # averaged projection +ProjAngleChannels[0,i,:,:] = AverageProj + + +filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits' + +data0 = read_fits(pathname0 + filename0) + +pathname10 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle10/' +filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits' + +data10 = read_fits(pathname10 + filename10) + +# Load a flat field (more are available, should we average over them?) +flat1 = read_fits('/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') + +# Apply flat field and display after flat-field correction and negative log +data0_rel = np.zeros(np.shape(flat1), dtype = float) +nonzero = flat1 > 0 +data0_rel[nonzero] = data0[nonzero] / flat1[nonzero] +data10_rel = np.zeros(np.shape(flat1), dtype = float) +data10_rel[nonzero] = data10[nonzero] / flat1[nonzero] + +plt.figure() +plt.imshow(data0_rel) +plt.colorbar() +plt.show() + +plt.figure() +plt.imshow(-np.log(data0_rel)) +plt.colorbar() +plt.show() + +plt.figure() +plt.imshow(data10_rel) +plt.colorbar() +plt.show() + +plt.figure() +plt.imshow(-np.log(data10_rel)) +plt.colorbar() +plt.show() + +# Set up for loading all summed images at 250 angles. +pathname = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle{}/' +filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits' + +# Dimensions +num_angles = 250 +imsize = 512 + +# Initialise array +data = np.zeros((num_angles,imsize,imsize)) + +# Load only 0-249, as 250 is at repetition of zero degrees just like image 0 +for i in range(0,250): + curimfile = (pathname + filename).format(i, i+4675) + data[i,:,:] = read_fits(curimfile) + +# Apply flat field and take negative log +nonzero = flat1 > 0 +for i in range(0,250): + data[i,nonzero] = data[i,nonzero]/flat1[nonzero] + +eqzero = data == 0 +data[eqzero] = 1 + +data_rel = -np.log(data) + +# Permute order to get: angles, vertical, horizontal, as default in framework. +data_rel = np.transpose(data_rel,(0,2,1)) + +# Set angles to use +angles = np.linspace(-np.pi,np.pi,num_angles,endpoint=False) + +# Create 3D acquisition geometry and acquisition data +ag = AcquisitionGeometry('parallel', + '3D', + angles, + pixel_num_h=imsize, + pixel_num_v=imsize) +b = AcquisitionData(data_rel, geometry=ag) + +# Reduce to single (noncentral) slice by extracting relevant parameters from data and its +# geometry. Perhaps create function to extract central slice automatically? +b2d = b.subset(vertical=128) +ag2d = AcquisitionGeometry('parallel', + '2D', + ag.angles, + pixel_num_h=ag.pixel_num_h) +b2d.geometry = ag2d + +# Create 2D image geometry +ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, + voxel_num_y=ag2d.pixel_num_h) + +# Create GPU projector/backprojector operator with ASTRA. +Aop = AstraProjectorSimple(ig2d,ag2d,'gpu') + +# Demonstrate operator is working by applying simple backprojection. +z = Aop.adjoint(b2d) +plt.figure() +plt.imshow(z.array) +plt.title('Simple backprojection') +plt.colorbar() +plt.show() + +# Set initial guess ImageData with zeros for algorithms, and algorithm options. +x_init = ImageData(np.zeros((imsize,imsize)), + geometry=ig2d) +opt_CGLS = {'tol': 1e-4, 'iter': 20} + +# Run CGLS algorithm and display reconstruction. +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b2d, opt_CGLS) +plt.figure() +plt.imshow(x_CGLS.array) +plt.title('CGLS') +plt.colorbar() +plt.show() + +plt.figure() +plt.semilogy(criter_CGLS) +plt.title('CGLS Criterion vs iterations') +plt.show() + + +f = Norm2sq(Aop,b2d,c=0.5) + +opt = {'tol': 1e-4, 'iter': 50} + +lamtv = 1.0 +# Repeat for FGP variant. +g_fgp = FGP_TV(lambdaReg = lamtv, + iterationsTV=50, + tolerance=1e-5, + methodTV=0, + nonnegativity=0, + printing=0, + device='cpu') + +x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp,opt) + +plt.figure() +plt.subplot(121) +plt.imshow(x_fista_fgp.array) +plt.title('FISTA FGP TV') +plt.subplot(122) +plt.semilogy(criter_fgp) +plt.show() -- cgit v1.2.3 From da5901a9fa132091eef500e73eeb9197ff2a3f05 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 22 Aug 2018 21:10:02 +0100 Subject: script to reconstruct multi-channel imat data updated --- Wrappers/Python/wip/demo_imat_multichan_RGLTK.py | 278 ++++++++++------------- 1 file changed, 115 insertions(+), 163 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py index dae5f3a..148aef8 100644 --- a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py +++ b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py @@ -1,8 +1,7 @@ # This script demonstrates how to load IMAT fits data # into the CIL optimisation framework and run reconstruction methods. # -# This demo loads the summedImg files which are the non-spectral images -# resulting from summing projections over all spectral channels. +# Demo to reconstruct energy-discretized channels of IMAT data # needs dxchange: conda install -c conda-forge dxchange # needs astropy: conda install astropy @@ -10,189 +9,142 @@ # All third-party imports. import numpy as np -from scipy.io import loadmat import matplotlib.pyplot as plt from dxchange.reader import read_fits +from astropy.io import fits # All own imports. -from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData +from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageGeometry, ImageData, DataContainer from ccpi.astra.ops import AstraProjectorSimple, AstraProjector3DSimple from ccpi.optimisation.algs import CGLS, FISTA from ccpi.optimisation.funcs import Norm2sq, Norm1 -from ccpi.plugins.regularisers import ROF_TV, FGP_TV, SB_TV - - -pathname0 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle0/' -filenameG = "IMAT00004675_Tomo_test_000_" +from ccpi.plugins.regularisers import FGP_TV +# set main parameters here n = 512 totalAngles = 250 # total number of projection angles # spectral discretization parameter -num_average = 120 -numChannels = 2970 -totChannels = round(numChannels/num_average) # total no. of averaged channels +num_average = 145 # channel discretization frequency (total number of averaged channels) +numChannels = 2970 # 2970 +totChannels = round(numChannels/num_average) # the resulting number of channels Projections_stack = np.zeros((num_average,n,n),dtype='uint16') ProjAngleChannels = np.zeros((totalAngles,totChannels,n,n),dtype='float32') -counterT = 0 -for i in range(0,2,1): - for j in range(0,num_average,1): - if counterT < 10: - outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathname0,filenameG,'0000',str(counterT)) - if ((counterT >= 10) & (counterT < 100)): - outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathname0,filenameG,'000',str(counterT)) - if ((counterT >= 100) & (counterT < 1000)): - outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathname0,filenameG,'00',str(counterT)) - if ((counterT >= 1000) & (counterT < 10000)): - outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathname0,filenameG,'0',str(counterT)) - Projections_stack[j,:,:] = read_fits(outfile) - counterT = counterT + 1 - AverageProj=np.mean(Projections_stack,axis=0) # averaged projection -ProjAngleChannels[0,i,:,:] = AverageProj - - -filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits' - -data0 = read_fits(pathname0 + filename0) - -pathname10 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle10/' -filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits' - -data10 = read_fits(pathname10 + filename10) - -# Load a flat field (more are available, should we average over them?) -flat1 = read_fits('/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') - -# Apply flat field and display after flat-field correction and negative log -data0_rel = np.zeros(np.shape(flat1), dtype = float) +######################################################################### +print ("Loading the data...") +MainPath = '/media/algol/HD-LXU3/DATA_DANIIL/' # path to data +pathname0 = '{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/Sample/') +counterFileName = 4675 +# A main loop over all available angles +for ll in range(0,totalAngles,1): + pathnameData = '{!s}{!s}{!s}/'.format(pathname0,'angle',str(ll)) + filenameCurr = '{!s}{!s}{!s}'.format('IMAT0000',str(counterFileName),'_Tomo_test_000_') + counterT = 0 + # loop over reduced channels (discretized) + for i in range(0,totChannels,1): + sumCount = 0 + # loop over actual channels to obtain averaged one + for j in range(0,num_average,1): + if counterT < 10: + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'0000',str(counterT)) + if ((counterT >= 10) & (counterT < 100)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'000',str(counterT)) + if ((counterT >= 100) & (counterT < 1000)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'00',str(counterT)) + if ((counterT >= 1000) & (counterT < 10000)): + outfile = '{!s}{!s}{!s}{!s}.fits'.format(pathnameData,filenameCurr,'0',str(counterT)) + try: + Projections_stack[j,:,:] = read_fits(outfile) + except: + print("Fits is corrupted, skipping no.", counterT) + sumCount -= 1 + counterT += 1 + sumCount += 1 + AverageProj=np.sum(Projections_stack,axis=0)/sumCount # averaged projection over "num_average" channels + ProjAngleChannels[ll,i,:,:] = AverageProj + print("Angle is processed", ll) + counterFileName += 1 +######################################################################### + +flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')) nonzero = flat1 > 0 -data0_rel[nonzero] = data0[nonzero] / flat1[nonzero] -data10_rel = np.zeros(np.shape(flat1), dtype = float) -data10_rel[nonzero] = data10[nonzero] / flat1[nonzero] - -plt.figure() -plt.imshow(data0_rel) -plt.colorbar() -plt.show() - -plt.figure() -plt.imshow(-np.log(data0_rel)) -plt.colorbar() -plt.show() - -plt.figure() -plt.imshow(data10_rel) -plt.colorbar() -plt.show() - -plt.figure() -plt.imshow(-np.log(data10_rel)) -plt.colorbar() -plt.show() - -# Set up for loading all summed images at 250 angles. -pathname = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle{}/' -filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits' - -# Dimensions -num_angles = 250 -imsize = 512 - -# Initialise array -data = np.zeros((num_angles,imsize,imsize)) - -# Load only 0-249, as 250 is at repetition of zero degrees just like image 0 -for i in range(0,250): - curimfile = (pathname + filename).format(i, i+4675) - data[i,:,:] = read_fits(curimfile) - # Apply flat field and take negative log -nonzero = flat1 > 0 -for i in range(0,250): - data[i,nonzero] = data[i,nonzero]/flat1[nonzero] - -eqzero = data == 0 -data[eqzero] = 1 +for ll in range(0,totalAngles,1): + for i in range(0,totChannels,1): + ProjAngleChannels[ll,i,nonzero] = ProjAngleChannels[ll,i,nonzero]/flat1[nonzero] -data_rel = -np.log(data) - -# Permute order to get: angles, vertical, horizontal, as default in framework. -data_rel = np.transpose(data_rel,(0,2,1)) +eqzero = ProjAngleChannels == 0 +ProjAngleChannels[eqzero] = 1 +ProjAngleChannels_NormLog = -np.log(ProjAngleChannels) # normalised and neg-log data +# extact sinogram over energy channels +selectedVertical_slice = 256 +sino_all_channels = ProjAngleChannels_NormLog[:,:,:,selectedVertical_slice] # Set angles to use -angles = np.linspace(-np.pi,np.pi,num_angles,endpoint=False) +angles = np.linspace(-np.pi,np.pi,totalAngles,endpoint=False) -# Create 3D acquisition geometry and acquisition data +# set the geometry +ig = ImageGeometry(n,n) ag = AcquisitionGeometry('parallel', - '3D', - angles, - pixel_num_h=imsize, - pixel_num_v=imsize) -b = AcquisitionData(data_rel, geometry=ag) - -# Reduce to single (noncentral) slice by extracting relevant parameters from data and its -# geometry. Perhaps create function to extract central slice automatically? -b2d = b.subset(vertical=128) -ag2d = AcquisitionGeometry('parallel', - '2D', - ag.angles, - pixel_num_h=ag.pixel_num_h) -b2d.geometry = ag2d - -# Create 2D image geometry -ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, - voxel_num_y=ag2d.pixel_num_h) - -# Create GPU projector/backprojector operator with ASTRA. -Aop = AstraProjectorSimple(ig2d,ag2d,'gpu') - -# Demonstrate operator is working by applying simple backprojection. -z = Aop.adjoint(b2d) -plt.figure() -plt.imshow(z.array) -plt.title('Simple backprojection') -plt.colorbar() -plt.show() - -# Set initial guess ImageData with zeros for algorithms, and algorithm options. -x_init = ImageData(np.zeros((imsize,imsize)), - geometry=ig2d) -opt_CGLS = {'tol': 1e-4, 'iter': 20} - -# Run CGLS algorithm and display reconstruction. -x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b2d, opt_CGLS) -plt.figure() -plt.imshow(x_CGLS.array) -plt.title('CGLS') -plt.colorbar() -plt.show() - -plt.figure() -plt.semilogy(criter_CGLS) -plt.title('CGLS Criterion vs iterations') -plt.show() - - -f = Norm2sq(Aop,b2d,c=0.5) - -opt = {'tol': 1e-4, 'iter': 50} - -lamtv = 1.0 -# Repeat for FGP variant. -g_fgp = FGP_TV(lambdaReg = lamtv, + '2D', + angles, + n,1) +Aop = AstraProjectorSimple(ig, ag, 'gpu') + + +# loop to reconstruct energy channels +REC_chan = np.zeros((totChannels, n, n), 'float32') +for i in range(0,totChannels,1): + sino_channel = sino_all_channels[:,i,:] # extract a sinogram for i-th channel + + print ("Initial guess") + x_init = ImageData(geometry=ig) + + # Create least squares object instance with projector and data. + print ("Create least squares object instance with projector and data.") + f = Norm2sq(Aop,DataContainer(sino_channel),c=0.5) + + print ("Run FISTA-TV for least squares") + lamtv = 10 + opt = {'tol': 1e-4, 'iter': 200} + g_fgp = FGP_TV(lambdaReg = lamtv, iterationsTV=50, - tolerance=1e-5, + tolerance=1e-6, methodTV=0, nonnegativity=0, printing=0, - device='cpu') - -x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp,opt) - -plt.figure() -plt.subplot(121) -plt.imshow(x_fista_fgp.array) -plt.title('FISTA FGP TV') -plt.subplot(122) -plt.semilogy(criter_fgp) -plt.show() + device='gpu') + + x_fista_fgp, it1, timing1, criter_fgp = FISTA(x_init, f, g_fgp, opt) + REC_chan[i,:,:] = x_fista_fgp.array + """ + plt.figure() + plt.subplot(121) + plt.imshow(x_fista_fgp.array, vmin=0, vmax=0.05) + plt.title('FISTA FGP TV') + plt.subplot(122) + plt.semilogy(criter_fgp) + plt.show() + """ + """ + print ("Run CGLS for least squares") + opt = {'tol': 1e-4, 'iter': 20} + x_init = ImageData(geometry=ig) + x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, DataContainer(sino_channel), opt=opt) + + plt.figure() + plt.imshow(x_CGLS.array,vmin=0, vmax=0.05) + plt.title('CGLS') + plt.show() + """ +# Saving images into fits using astrapy if required +add_val = np.min(REC_chan[:]) +REC_chan += abs(add_val) +REC_chan = REC_chan/np.max(REC_chan[:])*65535 +counter = 0 +filename = 'FISTA_TV_imat_slice' +for i in range(totChannels): + outfile = '{!s}_{!s}_{!s}.fits'.format(filename,str(selectedVertical_slice),str(counter)) + hdu = fits.PrimaryHDU(np.uint16(REC_chan[i,:,:])) + hdu.writeto(outfile, overwrite=True) + counter += 1 \ No newline at end of file -- cgit v1.2.3 From ed0bb129f1b03f1cbfee40a1d3444aefe176c7cf Mon Sep 17 00:00:00 2001 From: algol Date: Thu, 23 Aug 2018 00:21:02 +0100 Subject: some updates to demo --- Wrappers/Python/wip/demo_imat_multichan_RGLTK.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py index 148aef8..0ec116f 100644 --- a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py +++ b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py @@ -32,7 +32,7 @@ ProjAngleChannels = np.zeros((totalAngles,totChannels,n,n),dtype='float32') ######################################################################### print ("Loading the data...") -MainPath = '/media/algol/HD-LXU3/DATA_DANIIL/' # path to data +MainPath = '/media/algol/336F96987817D4B4/DATA/IMAT_DATA/' # path to data pathname0 = '{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/Sample/') counterFileName = 4675 # A main loop over all available angles @@ -105,7 +105,7 @@ for i in range(0,totChannels,1): f = Norm2sq(Aop,DataContainer(sino_channel),c=0.5) print ("Run FISTA-TV for least squares") - lamtv = 10 + lamtv = 5 opt = {'tol': 1e-4, 'iter': 200} g_fgp = FGP_TV(lambdaReg = lamtv, iterationsTV=50, @@ -138,13 +138,14 @@ for i in range(0,totChannels,1): plt.show() """ # Saving images into fits using astrapy if required -add_val = np.min(REC_chan[:]) -REC_chan += abs(add_val) -REC_chan = REC_chan/np.max(REC_chan[:])*65535 counter = 0 filename = 'FISTA_TV_imat_slice' for i in range(totChannels): + im = REC_chan[i,:,:] + add_val = np.min(im[:]) + im += abs(add_val) + im = im/np.max(im[:])*65535 outfile = '{!s}_{!s}_{!s}.fits'.format(filename,str(selectedVertical_slice),str(counter)) - hdu = fits.PrimaryHDU(np.uint16(REC_chan[i,:,:])) + hdu = fits.PrimaryHDU(np.uint16(im)) hdu.writeto(outfile, overwrite=True) counter += 1 \ No newline at end of file -- cgit v1.2.3 From d4b9d039581205a16e313816bb4e6b2934f62015 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Mon, 24 Sep 2018 11:43:47 +0100 Subject: Change local paths --- Wrappers/Python/wip/demo_imat_multichan_RGLTK.py | 7 ++++--- Wrappers/Python/wip/demo_imat_whitebeam.py | 8 ++++---- 2 files changed, 8 insertions(+), 7 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py index 0ec116f..59d634e 100644 --- a/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py +++ b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py @@ -32,8 +32,9 @@ ProjAngleChannels = np.zeros((totalAngles,totChannels,n,n),dtype='float32') ######################################################################### print ("Loading the data...") -MainPath = '/media/algol/336F96987817D4B4/DATA/IMAT_DATA/' # path to data -pathname0 = '{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/Sample/') +#MainPath = '/media/algol/336F96987817D4B4/DATA/IMAT_DATA/' # path to data +MainPath = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/' # path to data +pathname0 = '{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/Sample/') counterFileName = 4675 # A main loop over all available angles for ll in range(0,totalAngles,1): @@ -66,7 +67,7 @@ for ll in range(0,totalAngles,1): counterFileName += 1 ######################################################################### -flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_DATA/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')) +flat1 = read_fits('{!s}{!s}{!s}'.format(MainPath,'PSI_phantom_IMAT/DATA/','OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits')) nonzero = flat1 > 0 # Apply flat field and take negative log for ll in range(0,totalAngles,1): diff --git a/Wrappers/Python/wip/demo_imat_whitebeam.py b/Wrappers/Python/wip/demo_imat_whitebeam.py index 482c1ae..e2ffdb7 100644 --- a/Wrappers/Python/wip/demo_imat_whitebeam.py +++ b/Wrappers/Python/wip/demo_imat_whitebeam.py @@ -21,18 +21,18 @@ from ccpi.optimisation.algs import CGLS, FISTA from ccpi.optimisation.funcs import Norm2sq, Norm1 # Load and display a couple of summed projection as examples -pathname0 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle0/' +pathname0 = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/' filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits' data0 = read_fits(pathname0 + filename0) -pathname10 = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle10/' +pathname10 = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/' filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits' data10 = read_fits(pathname10 + filename10) # Load a flat field (more are available, should we average over them?) -flat1 = read_fits('/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') +flat1 = read_fits('/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') # Apply flat field and display after flat-field correction and negative log data0_rel = numpy.zeros(numpy.shape(flat1), dtype = float) @@ -58,7 +58,7 @@ plt.colorbar() plt.show() # Set up for loading all summed images at 250 angles. -pathname = '/media/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/DATA/Sample/angle{}/' +pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/' filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits' # Dimensions -- cgit v1.2.3 From 94eb6d54fb38f04999b5e8b1d0b2b7b66309b80f Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Thu, 14 Feb 2019 16:16:31 +0000 Subject: Updated paths in IMAT whitebeam script --- Wrappers/Python/wip/demo_imat_whitebeam.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/wip/demo_imat_whitebeam.py b/Wrappers/Python/wip/demo_imat_whitebeam.py index e2ffdb7..e0d213e 100644 --- a/Wrappers/Python/wip/demo_imat_whitebeam.py +++ b/Wrappers/Python/wip/demo_imat_whitebeam.py @@ -21,18 +21,18 @@ from ccpi.optimisation.algs import CGLS, FISTA from ccpi.optimisation.funcs import Norm2sq, Norm1 # Load and display a couple of summed projection as examples -pathname0 = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/' +pathname0 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle0/' filename0 = 'IMAT00004675_Tomo_test_000_SummedImg.fits' data0 = read_fits(pathname0 + filename0) -pathname10 = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/' +pathname10 = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle10/' filename10 = 'IMAT00004685_Tomo_test_000_SummedImg.fits' data10 = read_fits(pathname10 + filename10) # Load a flat field (more are available, should we average over them?) -flat1 = read_fits('/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') +flat1 = read_fits('/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/OpenBeam_aft1/IMAT00004932_Tomo_test_000_SummedImg.fits') # Apply flat field and display after flat-field correction and negative log data0_rel = numpy.zeros(numpy.shape(flat1), dtype = float) @@ -58,7 +58,7 @@ plt.colorbar() plt.show() # Set up for loading all summed images at 250 angles. -pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/' +pathname = '/media/newhd/shared/Data/neutrondata/PSI_phantom_IMAT/DATA/Sample/angle{}/' filename = 'IMAT0000{}_Tomo_test_000_SummedImg.fits' # Dimensions -- cgit v1.2.3