summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers')
-rw-r--r--Wrappers/Python/wip/demo_imat_multichan_RGLTK.py278
1 files changed, 115 insertions, 163 deletions
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
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-14 22:01:40 +0000