diff options
| -rw-r--r-- | Wrappers/Python/wip/demo_imat_multichan_RGLTK.py | 198 |
1 files changed, 198 insertions, 0 deletions
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() |