Added working demo of colourbay data load and recon. (#132)

* Added working demo of colourbay data load and recon.

* Added IMAT white-beam demo loading summed fits files

* corrections to normalization and log with zeroes in flats

* IMAT data  multichannel script started

* script to reconstruct multi-channel imat data updated

* some updates to demo

3 files changed, 426 insertions, 0 deletions

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
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..0ec116f
--- /dev/null
+++ b/Wrappers/Python/wip/demo_imat_multichan_RGLTK.py
@@ -0,0 +1,151 @@
+# This script demonstrates how to load IMAT fits data
+# into the CIL optimisation framework and run reconstruction methods.
+#
+# Demo to reconstruct energy-discretized channels of IMAT data
+
+# needs dxchange: conda install -c conda-forge dxchange
+# needs astropy: conda install astropy
+
+
+# All third-party imports.
+import numpy as np
+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, 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 FGP_TV
+
+# set main parameters here
+n = 512 
+totalAngles = 250 # total number of projection angles
+# spectral discretization parameter
+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')
+
+#########################################################################
+print ("Loading the 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 
+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
+# Apply flat field and take negative log
+for ll in range(0,totalAngles,1):
+    for i in range(0,totChannels,1):
+        ProjAngleChannels[ll,i,nonzero] = ProjAngleChannels[ll,i,nonzero]/flat1[nonzero]
+
+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,totalAngles,endpoint=False)
+
+# set the geometry
+ig = ImageGeometry(n,n)
+ag = AcquisitionGeometry('parallel',
+                             '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 = 5
+    opt = {'tol': 1e-4, 'iter': 200}
+    g_fgp = FGP_TV(lambdaReg = lamtv,
+                 iterationsTV=50,
+                 tolerance=1e-6,
+                 methodTV=0,
+                 nonnegativity=0,
+                 printing=0,
+                 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
+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(im))
+    hdu.writeto(outfile, overwrite=True)
+    counter += 1
\ No newline at end of file
diff --git a/Wrappers/Python/wip/demo_imat_whitebeam.py b/Wrappers/Python/wip/demo_imat_whitebeam.py
new file mode 100644
index 0000000..482c1ae
--- /dev/null
+++ b/Wrappers/Python/wip/demo_imat_whitebeam.py
@@ -0,0 +1,138 @@
+# 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/algol/HD-LXU3/DATA_DANIIL/PSI_DATA/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/'
+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 = 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()
+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/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 = 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
+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))
+
+# 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