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') 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') 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') 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') 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') 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') 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') 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') 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 From 5ff4c802b06b9e7c1fe557c324bd7041b331f4a8 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 15 Feb 2019 16:09:34 +0000 Subject: add norm to DataContainer (#195) --- Wrappers/Python/ccpi/framework.py | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 9938fb7..318eb92 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -720,6 +720,10 @@ class DataContainer(object): ## reductions def sum(self, out=None, *args, **kwargs): return self.as_array().sum(*args, **kwargs) + def norm(self): + '''return the norm of the DataContainer''' + y = self.as_array() + return numpy.dot(y, y.conjugate()) class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' -- cgit v1.2.3 From 0abb52cd5929d809b3ac1a651d85f465ada80137 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 27 Feb 2019 06:40:14 +0000 Subject: added squared_norm (#204) * added squared_norm closes #203 * added norm and squared_norm closes #203 * Power method uses norm, bugfixes * fix power method --- Wrappers/Python/ccpi/framework.py | 14 ++++++++++---- Wrappers/Python/ccpi/optimisation/algs.py | 5 ++--- Wrappers/Python/ccpi/optimisation/funcs.py | 2 +- Wrappers/Python/ccpi/optimisation/ops.py | 5 +++-- Wrappers/Python/conda-recipe/run_test.py | 28 ++++++++++++++++++++++------ 5 files changed, 38 insertions(+), 16 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 318eb92..e1a2dff 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -3,7 +3,7 @@ # Visual Analytics and Imaging System Group of the Science Technology # Facilities Council, STFC -# Copyright 2018 Edoardo Pasca +# Copyright 2018-2019 Edoardo Pasca # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -26,6 +26,7 @@ import numpy import sys from datetime import timedelta, datetime import warnings +from functools import reduce def find_key(dic, val): """return the key of dictionary dic given the value""" @@ -720,10 +721,15 @@ class DataContainer(object): ## reductions def sum(self, out=None, *args, **kwargs): return self.as_array().sum(*args, **kwargs) - def norm(self): - '''return the norm of the DataContainer''' - y = self.as_array() + def squared_norm(self): + '''return the squared euclidean norm of the DataContainer viewed as a vector''' + shape = self.shape + size = reduce(lambda x,y:x*y, shape, 1) + y = numpy.reshape(self.as_array(), (size, )) return numpy.dot(y, y.conjugate()) + def norm(self): + '''return the euclidean norm of the DataContainer viewed as a vector''' + return numpy.sqrt(self.squared_norm()) class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py index a736277..15638a9 100755 --- a/Wrappers/Python/ccpi/optimisation/algs.py +++ b/Wrappers/Python/ccpi/optimisation/algs.py @@ -158,9 +158,8 @@ def FBPD(x_init, operator=None, constraint=None, data_fidelity=None,\ memopt = opt['memopts'] if 'memopts' in opt.keys() else False # step-sizes - tau = 2 / (data_fidelity.L + 2); - sigma = (1/tau - data_fidelity.L/2) / regulariser.L; - + tau = 2 / (data_fidelity.L + 2) + sigma = (1/tau - data_fidelity.L/2) / regulariser.L inv_sigma = 1/sigma # initialization diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py index c7a6143..9b9fc36 100755 --- a/Wrappers/Python/ccpi/optimisation/funcs.py +++ b/Wrappers/Python/ccpi/optimisation/funcs.py @@ -217,10 +217,10 @@ class ZeroFun(Function): class Norm1(Function): def __init__(self,gamma): + super(Norm1, self).__init__() self.gamma = gamma self.L = 1 self.sign_x = None - super(Norm1, self).__init__() def __call__(self,x,out=None): if out is None: diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py index 450b084..b2f996d 100755 --- a/Wrappers/Python/ccpi/optimisation/ops.py +++ b/Wrappers/Python/ccpi/optimisation/ops.py @@ -209,10 +209,11 @@ def PowerMethodNonsquare(op,numiters , x0=None): # Loop for it in numpy.arange(numiters): x1 = op.adjoint(op.direct(x0)) - x1norm = numpy.sqrt((x1*x1).sum()) + #x1norm = numpy.sqrt((x1*x1).sum()) + x1norm = x1.norm() #print ("x0 **********" ,x0) #print ("x1 **********" ,x1) - s[it] = (x1*x0).sum() / (x0*x0).sum() + s[it] = (x1*x0).sum() / (x0.squared_norm()) x0 = (1.0/x1norm)*x1 return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py index 5bf6538..b52af55 100755 --- a/Wrappers/Python/conda-recipe/run_test.py +++ b/Wrappers/Python/conda-recipe/run_test.py @@ -19,6 +19,7 @@ from ccpi.optimisation.funcs import Norm2 from ccpi.optimisation.ops import LinearOperatorMatrix from ccpi.optimisation.ops import TomoIdentity from ccpi.optimisation.ops import Identity +from ccpi.optimisation.ops import PowerMethodNonsquare import numpy.testing @@ -494,6 +495,9 @@ class TestDataContainer(unittest.TestCase): except AssertionError as err: res = False print(err) + print("expected " , second) + print("actual " , first) + self.assertTrue(res) def test_DataContainerChaining(self): dc = self.create_DataContainer(256,256,256,1) @@ -501,7 +505,13 @@ class TestDataContainer(unittest.TestCase): dc.add(9,out=dc)\ .subtract(1,out=dc) self.assertEqual(1+9-1,dc.as_array().flatten()[0]) - + def test_reduction(self): + print ("test reductions") + dc = self.create_DataContainer(2,2,2,value=1) + sqnorm = dc.squared_norm() + norm = dc.norm() + self.assertEqual(sqnorm, 8.0) + numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7) @@ -643,7 +653,8 @@ class TestAlgorithms(unittest.TestCase): else: self.assertTrue(cvx_not_installable) - def test_FBPD_Norm1_cvx(self): + def skip_test_FBPD_Norm1_cvx(self): + print ("test_FBPD_Norm1_cvx") if not cvx_not_installable: opt = {'memopt': True} # Problem data. @@ -663,12 +674,15 @@ class TestAlgorithms(unittest.TestCase): # Regularization parameter lam = 10 opt = {'memopt': True} + # Initial guess + x_init = DataContainer(np.random.randn(n, 1)) + # Create object instances with the test data A and b. f = Norm2sq(A, b, c=0.5, memopt=True) + f.L = PowerMethodNonsquare(A, 25, x_init)[0] + print ("Lipschitz", f.L) g0 = ZeroFun() - # Initial guess - x_init = DataContainer(np.zeros((n, 1))) # Create 1-norm object instance g1 = Norm1(lam) @@ -710,7 +724,7 @@ class TestAlgorithms(unittest.TestCase): # Set up phantom size NxN by creating ImageGeometry, initialising the # ImageData object with this geometry and empty array and finally put some # data into its array, and display as image. - def test_FISTA_denoise_cvx(self): + def skip_test_FISTA_denoise_cvx(self): if not cvx_not_installable: opt = {'memopt': True} N = 64 @@ -733,6 +747,8 @@ class TestAlgorithms(unittest.TestCase): # Data fidelity term f_denoise = Norm2sq(I, y, c=0.5, memopt=True) + x_init = ImageData(geometry=ig) + f_denoise.L = PowerMethodNonsquare(I, 25, x_init)[0] # 1-norm regulariser lam1_denoise = 1.0 @@ -759,7 +775,7 @@ class TestAlgorithms(unittest.TestCase): # Compare to CVXPY # Construct the problem. - x1_denoise = Variable(N**2, 1) + x1_denoise = Variable(N**2) objective1_denoise = Minimize( 0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1)) prob1_denoise = Problem(objective1_denoise) -- cgit v1.2.3 From 13054c64cff2b506a7affa9d21445f368abbb15b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 28 Feb 2019 18:58:33 +0000 Subject: Proposal of Algorithm class (#179) * initial revision * Removed class members of Algorithm class added update_objective * initial version. Fix inline __idiv__ * First implementation of CompositeOperator/DataContainer * removed __getitem__ added get_item added shape * added CGLS * working unit test, initial tomography test * added reverse multiplication of operator with number * added operators directory * fixed typo * added unittest for CompositeDataContainer * fix TomoIdentity with scalar * check numerical types from numpy * add default stop criterion and run method * add run method * first working implementation of CGLS with CompositeOperator/DataContainer notice problem with _rmul_ and _mul_ methods precedence with numpy. * new Algorithm class and algorithms in separate files Added new Algorithm class and derivatives in different files for GradientDescent, CGLS, FBPD, FISTA * added algorithms and restored CIL_VERSION env variable * removed Algorithms.py * modified run and renamed a few members/methods * uses squared_norm * renamed get_current_objective to get_last_objective update_objective can be issued every N iteration, default 1. fixed run method to run N iterations within the stop criterion. * load class as module files * force py line endings to LF * updates * call super __init__ as first thing * unit tests are now to be found in test directory unit tests are now split in several files in the directory test * install algorithms module * Implementation with Algorithm * skip Reader tests * unittest for linux * commented not needed import Iterable * removed explicit return from __init__ * remove composite operator file --- Wrappers/Python/ccpi/framework.py | 12 +- .../ccpi/optimisation/algorithms/Algorithm.py | 157 ++++ .../Python/ccpi/optimisation/algorithms/CGLS.py | 87 ++ .../Python/ccpi/optimisation/algorithms/FBPD.py | 86 ++ .../Python/ccpi/optimisation/algorithms/FISTA.py | 121 +++ .../optimisation/algorithms/GradientDescent.py | 72 ++ .../ccpi/optimisation/algorithms/__init__.py | 29 + Wrappers/Python/ccpi/optimisation/ops.py | 48 +- Wrappers/Python/conda-recipe/meta.yaml | 9 + Wrappers/Python/conda-recipe/run_test.py | 985 --------------------- Wrappers/Python/setup.py | 8 +- Wrappers/Python/test/TestReader.py | 134 --- Wrappers/Python/test/__init__.py | 0 Wrappers/Python/test/skip_TestReader.py | 134 +++ Wrappers/Python/test/test_DataContainer.py | 499 +++++++++++ Wrappers/Python/test/test_NexusReader.py | 96 ++ Wrappers/Python/test/test_algorithms.py | 123 +++ Wrappers/Python/test/test_run_test.py | 435 +++++++++ 18 files changed, 1903 insertions(+), 1132 deletions(-) create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py create mode 100755 Wrappers/Python/ccpi/optimisation/algorithms/__init__.py delete mode 100755 Wrappers/Python/conda-recipe/run_test.py delete mode 100644 Wrappers/Python/test/TestReader.py create mode 100644 Wrappers/Python/test/__init__.py create mode 100644 Wrappers/Python/test/skip_TestReader.py create mode 100755 Wrappers/Python/test/test_DataContainer.py create mode 100755 Wrappers/Python/test/test_NexusReader.py create mode 100755 Wrappers/Python/test/test_algorithms.py create mode 100755 Wrappers/Python/test/test_run_test.py (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index e1a2dff..0c23628 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -473,7 +473,10 @@ class DataContainer(object): else: raise ValueError('*:Wrong shape: {0} and {1}'.format(self.shape, other.shape)) - elif isinstance(other, (int, float, complex)): + elif isinstance(other, (int, float, complex,\ + numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\ + numpy.float, numpy.float16, numpy.float32, numpy.float64, \ + numpy.complex)): return type(self)(self.as_array() * other, deep_copy=True, dimension_labels=self.dimension_labels, @@ -559,6 +562,8 @@ class DataContainer(object): # __isub__ def __idiv__(self, other): + return self.__itruediv__(other) + def __itruediv__(self, other): if isinstance(other, (int, float)) : numpy.divide(self.array, other, out=self.array) elif issubclass(type(other), DataContainer): @@ -632,6 +637,10 @@ class DataContainer(object): if out is None: if isinstance(x2, (int, float, complex)): out = pwop(self.as_array() , x2 , *args, **kwargs ) + elif isinstance(x2, (numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\ + numpy.float, numpy.float16, numpy.float32, numpy.float64, \ + numpy.complex)): + out = pwop(self.as_array() , x2 , *args, **kwargs ) elif issubclass(type(x2) , DataContainer): out = pwop(self.as_array() , x2.as_array() , *args, **kwargs ) return type(self)(out, @@ -731,6 +740,7 @@ class DataContainer(object): '''return the euclidean norm of the DataContainer viewed as a vector''' return numpy.sqrt(self.squared_norm()) + class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' def __init__(self, diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py new file mode 100755 index 0000000..680b268 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -0,0 +1,157 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import time +from numbers import Integral + +class Algorithm(object): + '''Base class for iterative algorithms + + provides the minimal infrastructure. + Algorithms are iterables so can be easily run in a for loop. They will + stop as soon as the stop cryterion is met. + The user is required to implement the set_up, __init__, update and + and update_objective methods + + A courtesy method run is available to run n iterations. The method accepts + a callback function that receives the current iteration number and the actual objective + value and can be used to trigger print to screens and other user interactions. The run + method will stop when the stopping cryterion is met. + ''' + + def __init__(self): + '''Constructor + + Set the minimal number of parameters: + iteration: current iteration number + max_iteration: maximum number of iterations + memopt: whether to use memory optimisation () + timing: list to hold the times it took to run each iteration + update_objectice_interval: the interval every which we would save the current + objective. 1 means every iteration, 2 every 2 iteration + and so forth. This is by default 1 and should be increased + when evaluating the objective is computationally expensive. + ''' + self.iteration = 0 + self.__max_iteration = 0 + self.__loss = [] + self.memopt = False + self.timing = [] + self.update_objective_interval = 1 + def set_up(self, *args, **kwargs): + '''Set up the algorithm''' + raise NotImplementedError() + def update(self): + '''A single iteration of the algorithm''' + raise NotImplementedError() + + def should_stop(self): + '''default stopping cryterion: number of iterations + + The user can change this in concrete implementatition of iterative algorithms.''' + return self.max_iteration_stop_cryterion() + + def max_iteration_stop_cryterion(self): + '''default stop cryterion for iterative algorithm: max_iteration reached''' + return self.iteration >= self.max_iteration + def __iter__(self): + '''Algorithm is an iterable''' + return self + def next(self): + '''Algorithm is an iterable + + python2 backwards compatibility''' + return self.__next__() + def __next__(self): + '''Algorithm is an iterable + + calling this method triggers update and update_objective + ''' + if self.should_stop(): + raise StopIteration() + else: + time0 = time.time() + self.update() + self.timing.append( time.time() - time0 ) + if self.iteration % self.update_objective_interval == 0: + self.update_objective() + self.iteration += 1 + def get_output(self): + '''Returns the solution found''' + return self.x + def get_last_loss(self): + '''Returns the last stored value of the loss function + + if update_objective_interval is 1 it is the value of the objective at the current + iteration. If update_objective_interval > 1 it is the last stored value. + ''' + return self.__loss[-1] + def get_last_objective(self): + '''alias to get_last_loss''' + return self.get_last_loss() + def update_objective(self): + '''calculates the objective with the current solution''' + raise NotImplementedError() + @property + def loss(self): + '''returns the list of the values of the objective during the iteration + + The length of this list may be shorter than the number of iterations run when + the update_objective_interval > 1 + ''' + return self.__loss + @property + def objective(self): + '''alias of loss''' + return self.loss + @property + def max_iteration(self): + '''gets the maximum number of iterations''' + return self.__max_iteration + @max_iteration.setter + def max_iteration(self, value): + '''sets the maximum number of iterations''' + assert isinstance(value, int) + self.__max_iteration = value + @property + def update_objective_interval(self): + return self.__update_objective_interval + @update_objective_interval.setter + def update_objective_interval(self, value): + if isinstance(value, Integral): + if value >= 1: + self.__update_objective_interval = value + else: + raise ValueError('Update objective interval must be an integer >= 1') + else: + raise ValueError('Update objective interval must be an integer >= 1') + def run(self, iterations, verbose=False, callback=None): + '''run n iterations and update the user with the callback if specified''' + if self.should_stop(): + print ("Stop cryterion has been reached.") + i = 0 + for _ in self: + if verbose: + print ("Iteration {}/{}, objective {}".format(self.iteration, + self.max_iteration, self.get_last_objective()) ) + if callback is not None: + callback(self.iteration, self.get_last_objective()) + i += 1 + if i == iterations: + break + diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py new file mode 100755 index 0000000..7194eb8 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py @@ -0,0 +1,87 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:11:23 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +#from collections.abc import Iterable +class CGLS(Algorithm): + + '''Conjugate Gradient Least Squares algorithm + + Parameters: + x_init: initial guess + operator: operator for forward/backward projections + data: data to operate on + ''' + def __init__(self, **kwargs): + super(CGLS, self).__init__() + self.x = kwargs.get('x_init', None) + self.operator = kwargs.get('operator', None) + self.data = kwargs.get('data', None) + if self.x is not None and self.operator is not None and \ + self.data is not None: + print ("Calling from creator") + self.set_up(x_init =kwargs['x_init'], + operator=kwargs['operator'], + data =kwargs['data']) + + def set_up(self, x_init, operator , data ): + + self.r = data.copy() + self.x = x_init.copy() + + self.operator = operator + self.d = operator.adjoint(self.r) + + + self.normr2 = self.d.squared_norm() + #if isinstance(self.normr2, Iterable): + # self.normr2 = sum(self.normr2) + #self.normr2 = numpy.sqrt(self.normr2) + #print ("set_up" , self.normr2) + + def update(self): + + Ad = self.operator.direct(self.d) + #norm = (Ad*Ad).sum() + #if isinstance(norm, Iterable): + # norm = sum(norm) + norm = Ad.squared_norm() + + alpha = self.normr2/norm + self.x += (self.d * alpha) + self.r -= (Ad * alpha) + s = self.operator.adjoint(self.r) + + normr2_new = s.squared_norm() + #if isinstance(normr2_new, Iterable): + # normr2_new = sum(normr2_new) + #normr2_new = numpy.sqrt(normr2_new) + #print (normr2_new) + + beta = normr2_new/self.normr2 + self.normr2 = normr2_new + self.d = s + beta*self.d + + def update_objective(self): + self.loss.append(self.r.squared_norm()) \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py new file mode 100755 index 0000000..322e9eb --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py @@ -0,0 +1,86 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:09:03 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.funcs import ZeroFun + +class FBPD(Algorithm): + '''FBPD Algorithm + + Parameters: + x_init: initial guess + f: constraint + g: data fidelity + h: regularizer + opt: additional algorithm + ''' + constraint = None + data_fidelity = None + regulariser = None + def __init__(self, **kwargs): + pass + def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\ + regulariser=None, opt=None): + + # default inputs + if constraint is None: + self.constraint = ZeroFun() + else: + self.constraint = constraint + if data_fidelity is None: + data_fidelity = ZeroFun() + else: + self.data_fidelity = data_fidelity + if regulariser is None: + self.regulariser = ZeroFun() + else: + self.regulariser = regulariser + + # algorithmic parameters + + + # step-sizes + self.tau = 2 / (self.data_fidelity.L + 2) + self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L + + self.inv_sigma = 1/self.sigma + + # initialization + self.x = x_init + self.y = operator.direct(self.x) + + + def update(self): + + # primal forward-backward step + x_old = self.x + self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) ) + self.x = self.constraint.prox(self.x, self.tau); + + # dual forward-backward step + self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old); + self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma); + + # time and criterion + self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x)) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py new file mode 100755 index 0000000..bc4489e --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 21 11:07:30 2019 + +@author: ofn77899 +""" + +from ccpi.optimisation.algorithms import Algorithm +from ccpi.optimisation.funcs import ZeroFun +import numpy + +class FISTA(Algorithm): + '''Fast Iterative Shrinkage-Thresholding Algorithm + + Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding + algorithm for linear inverse problems. + SIAM journal on imaging sciences,2(1), pp.183-202. + + Parameters: + x_init: initial guess + f: data fidelity + g: regularizer + h: + opt: additional algorithm + ''' + + def __init__(self, **kwargs): + '''initialisation can be done at creation time if all + proper variables are passed or later with set_up''' + super(FISTA, self).__init__() + self.f = None + self.g = None + self.invL = None + self.t_old = 1 + args = ['x_init', 'f', 'g', 'opt'] + for k,v in kwargs.items(): + if k in args: + args.pop(args.index(k)) + if len(args) == 0: + return self.set_up(kwargs['x_init'], + f=kwargs['f'], + g=kwargs['g'], + opt=kwargs['opt']) + + def set_up(self, x_init, f=None, g=None, opt=None): + + # default inputs + if f is None: + self.f = ZeroFun() + else: + self.f = f + if g is None: + g = ZeroFun() + self.g = g + else: + self.g = g + + # algorithmic parameters + if opt is None: + opt = {'tol': 1e-4, 'memopt':False} + + self.tol = opt['tol'] if 'tol' in opt.keys() else 1e-4 + memopt = opt['memopt'] if 'memopt' in opt.keys() else False + self.memopt = memopt + + # initialization + if memopt: + self.y = x_init.clone() + self.x_old = x_init.clone() + self.x = x_init.clone() + self.u = x_init.clone() + else: + self.x_old = x_init.copy() + self.y = x_init.copy() + + #timing = numpy.zeros(max_iter) + #criter = numpy.zeros(max_iter) + + + self.invL = 1/f.L + + self.t_old = 1 + + def update(self): + # algorithm loop + #for it in range(0, max_iter): + + if self.memopt: + # u = y - invL*f.grad(y) + # store the result in x_old + self.f.gradient(self.y, out=self.u) + self.u.__imul__( -self.invL ) + self.u.__iadd__( self.y ) + # x = g.prox(u,invL) + self.g.proximal(self.u, self.invL, out=self.x) + + self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) + + # y = x + (t_old-1)/t*(x-x_old) + self.x.subtract(self.x_old, out=self.y) + self.y.__imul__ ((self.t_old-1)/self.t) + self.y.__iadd__( self.x ) + + self.x_old.fill(self.x) + self.t_old = self.t + + + else: + u = self.y - self.invL*self.f.grad(self.y) + + self.x = self.g.prox(u,self.invL) + + self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2))) + + self.y = self.x + (self.t_old-1)/self.t*(self.x-self.x_old) + + self.x_old = self.x.copy() + self.t_old = self.t + + def update_objective(self): + self.loss.append( self.f(self.x) + self.g(self.x) ) \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py new file mode 100755 index 0000000..7794b4d --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py @@ -0,0 +1,72 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:05:09 2019 + +@author: ofn77899 +""" +from ccpi.optimisation.algorithms import Algorithm + +class GradientDescent(Algorithm): + '''Implementation of Gradient Descent algorithm + ''' + + def __init__(self, **kwargs): + '''initialisation can be done at creation time if all + proper variables are passed or later with set_up''' + super(GradientDescent, self).__init__() + self.x = None + self.rate = 0 + self.objective_function = None + self.regulariser = None + args = ['x_init', 'objective_function', 'rate'] + for k,v in kwargs.items(): + if k in args: + args.pop(args.index(k)) + if len(args) == 0: + return self.set_up(x_init=kwargs['x_init'], + objective_function=kwargs['objective_function'], + rate=kwargs['rate']) + + def should_stop(self): + '''stopping cryterion, currently only based on number of iterations''' + return self.iteration >= self.max_iteration + + def set_up(self, x_init, objective_function, rate): + '''initialisation of the algorithm''' + self.x = x_init.copy() + if self.memopt: + self.x_update = x_init.copy() + self.objective_function = objective_function + self.rate = rate + self.loss.append(objective_function(x_init)) + self.iteration = 0 + + def update(self): + '''Single iteration''' + if self.memopt: + self.objective_function.gradient(self.x, out=self.x_update) + self.x_update *= -self.rate + self.x += self.x_update + else: + self.x += -self.rate * self.objective_function.grad(self.x) + + def update_objective(self): + self.loss.append(self.objective_function(self.x)) + \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py new file mode 100755 index 0000000..52fe6d7 --- /dev/null +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -0,0 +1,29 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2019 Edoardo Pasca + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Created on Thu Feb 21 11:03:13 2019 + +@author: ofn77899 +""" + +from .Algorithm import Algorithm +from .CGLS import CGLS +from .GradientDescent import GradientDescent +from .FISTA import FISTA +from .FBPD import FBPD \ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py index b2f996d..e9e7f44 100755 --- a/Wrappers/Python/ccpi/optimisation/ops.py +++ b/Wrappers/Python/ccpi/optimisation/ops.py @@ -24,26 +24,49 @@ from ccpi.framework import AcquisitionData from ccpi.framework import ImageData from ccpi.framework import ImageGeometry from ccpi.framework import AcquisitionGeometry - +from numbers import Number # Maybe operators need to know what types they take as inputs/outputs # to not just use generic DataContainer class Operator(object): + '''Operator that maps from a space X -> Y''' + def __init__(self, **kwargs): + self.scalar = 1 + def is_linear(self): + '''Returns if the operator is linear''' + return False def direct(self,x, out=None): - return x - def adjoint(self,x, out=None): - return x + raise NotImplementedError def size(self): # To be defined for specific class raise NotImplementedError - def get_max_sing_val(self): + def norm(self): raise NotImplementedError def allocate_direct(self): + '''Allocates memory on the Y space''' raise NotImplementedError def allocate_adjoint(self): + '''Allocates memory on the X space''' + raise NotImplementedError + def range_dim(self): raise NotImplementedError + def domain_dim(self): + raise NotImplementedError + def __rmul__(self, other): + '''reverse multiplication of Operator with number sets the variable scalar in the Operator''' + assert isinstance(other, Number) + self.scalar = other + return self +class LinearOperator(Operator): + '''Operator that maps from a space X -> Y''' + def is_linear(self): + '''Returns if the operator is linear''' + return True + def adjoint(self,x, out=None): + raise NotImplementedError + class Identity(Operator): def __init__(self): self.s1 = 1.0 @@ -70,21 +93,26 @@ class Identity(Operator): class TomoIdentity(Operator): def __init__(self, geometry, **kwargs): + super(TomoIdentity, self).__init__() self.s1 = 1.0 self.geometry = geometry - super(TomoIdentity, self).__init__() + def direct(self,x,out=None): + if out is None: + if self.scalar != 1: + return x * self.scalar return x.copy() else: + if self.scalar != 1: + out.fill(x * self.scalar) + return out.fill(x) + return def adjoint(self,x, out=None): - if out is None: - return x.copy() - else: - out.fill(x) + return self.direct(x, out) def size(self): return NotImplemented diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index 1b7cae6..8ded429 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -12,6 +12,15 @@ test: requires: - python-wget - cvxpy # [not win] + + source_files: + - ./test # [win] + - ./ccpi/Wrappers/Python/test # [not win] + + commands: + - python -c "import os; print (os.getcwd())" + - python -m unittest discover # [win] + - python -m unittest discover -s ccpi/Wrappers/Python/test # [not win] requirements: build: diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py deleted file mode 100755 index b52af55..0000000 --- a/Wrappers/Python/conda-recipe/run_test.py +++ /dev/null @@ -1,985 +0,0 @@ -import unittest -import numpy -import numpy as np -from ccpi.framework import DataContainer -from ccpi.framework import ImageData -from ccpi.framework import AcquisitionData -from ccpi.framework import ImageGeometry -from ccpi.framework import AcquisitionGeometry -import sys -from timeit import default_timer as timer -from ccpi.optimisation.algs import FISTA -from ccpi.optimisation.algs import FBPD -from ccpi.optimisation.funcs import Norm2sq -from ccpi.optimisation.funcs import ZeroFun -from ccpi.optimisation.funcs import Norm1 -from ccpi.optimisation.funcs import TV2D -from ccpi.optimisation.funcs import Norm2 - -from ccpi.optimisation.ops import LinearOperatorMatrix -from ccpi.optimisation.ops import TomoIdentity -from ccpi.optimisation.ops import Identity -from ccpi.optimisation.ops import PowerMethodNonsquare - - -import numpy.testing -import wget -import os -from ccpi.io.reader import NexusReader - - -try: - from cvxpy import * - cvx_not_installable = False -except ImportError: - cvx_not_installable = True - - -def aid(x): - # This function returns the memory - # block address of an array. - return x.__array_interface__['data'][0] - - -def dt(steps): - return steps[-1] - steps[-2] - - -class TestDataContainer(unittest.TestCase): - def create_simple_ImageData(self): - N = 64 - ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) - Phantom = ImageData(geometry=ig) - - x = Phantom.as_array() - - x[int(round(N/4)):int(round(3*N/4)), - int(round(N/4)):int(round(3*N/4))] = 0.5 - x[int(round(N/8)):int(round(7*N/8)), - int(round(3*N/8)):int(round(5*N/8))] = 1 - - return (ig, Phantom) - - def create_DataContainer(self, X,Y,Z, value=1): - steps = [timer()] - a = value * numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - return ds - - def test_creation_nocopy(self): - shape = (2, 3, 4, 5) - size = shape[0] - for i in range(1, len(shape)): - size = size * shape[i] - #print("a refcount " , sys.getrefcount(a)) - a = numpy.asarray([i for i in range(size)]) - #print("a refcount " , sys.getrefcount(a)) - a = numpy.reshape(a, shape) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z', 'W']) - #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a), 3) - self.assertEqual(ds.dimension_labels, {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'}) - - def testGb_creation_nocopy(self): - X, Y, Z = 512, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - print("test clone") - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a), 3) - ds1 = ds.copy() - self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) - ds1 = ds.clone() - self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) - - def testInlineAlgebra(self): - print("Test Inline Algebra") - X, Y, Z = 1024, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - print(t0) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - #ds.__iadd__( 2 ) - ds += 2 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 3.) - #ds.__isub__( 2 ) - ds -= 2 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - #ds.__imul__( 2 ) - ds *= 2 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 2.) - #ds.__idiv__( 2 ) - ds /= 2 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - - ds1 = ds.copy() - #ds1.__iadd__( 1 ) - ds1 += 1 - #ds.__iadd__( ds1 ) - ds += ds1 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 3.) - #ds.__isub__( ds1 ) - ds -= ds1 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - #ds.__imul__( ds1 ) - ds *= ds1 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 2.) - #ds.__idiv__( ds1 ) - ds /= ds1 - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - - def test_unary_operations(self): - print("Test unary operations") - X, Y, Z = 1024, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = -numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - print(t0) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - - ds.sign(out=ds) - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], -1.) - - ds.abs(out=ds) - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], 1.) - - ds.__imul__(2) - ds.sqrt(out=ds) - steps.append(timer()) - print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0], - numpy.sqrt(2., dtype='float32')) - - def test_binary_operations(self): - self.binary_add() - self.binary_subtract() - self.binary_multiply() - self.binary_divide() - - def binary_add(self): - print("Test binary add") - X, Y, Z = 512, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - ds1 = ds.copy() - - steps.append(timer()) - ds.add(ds1, out=ds) - steps.append(timer()) - t1 = dt(steps) - print("ds.add(ds1, out=ds)", dt(steps)) - steps.append(timer()) - ds2 = ds.add(ds1) - steps.append(timer()) - t2 = dt(steps) - print("ds2 = ds.add(ds1)", dt(steps)) - - self.assertLess(t1, t2) - self.assertEqual(ds.as_array()[0][0][0], 2.) - - ds0 = ds - ds0.add(2, out=ds0) - steps.append(timer()) - print("ds0.add(2,out=ds0)", dt(steps), 3, ds0.as_array()[0][0][0]) - self.assertEqual(4., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) - ds3 = ds0.add(2) - steps.append(timer()) - print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0]) - dt2 = dt(steps) - self.assertLess(dt1, dt2) - - def binary_subtract(self): - print("Test binary subtract") - X, Y, Z = 512, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - ds1 = ds.copy() - - steps.append(timer()) - ds.subtract(ds1, out=ds) - steps.append(timer()) - t1 = dt(steps) - print("ds.subtract(ds1, out=ds)", dt(steps)) - self.assertEqual(0., ds.as_array()[0][0][0]) - - steps.append(timer()) - ds2 = ds.subtract(ds1) - self.assertEqual(-1., ds2.as_array()[0][0][0]) - - steps.append(timer()) - t2 = dt(steps) - print("ds2 = ds.subtract(ds1)", dt(steps)) - - self.assertLess(t1, t2) - - del ds1 - ds0 = ds.copy() - steps.append(timer()) - ds0.subtract(2, out=ds0) - #ds0.__isub__( 2 ) - steps.append(timer()) - print("ds0.subtract(2,out=ds0)", dt( - steps), -2., ds0.as_array()[0][0][0]) - self.assertEqual(-2., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) - ds3 = ds0.subtract(2) - steps.append(timer()) - print("ds3 = ds0.subtract(2)", dt(steps), 0., ds3.as_array()[0][0][0]) - dt2 = dt(steps) - self.assertLess(dt1, dt2) - self.assertEqual(-2., ds0.as_array()[0][0][0]) - self.assertEqual(-4., ds3.as_array()[0][0][0]) - - def binary_multiply(self): - print("Test binary multiply") - X, Y, Z = 1024, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - ds1 = ds.copy() - - steps.append(timer()) - ds.multiply(ds1, out=ds) - steps.append(timer()) - t1 = dt(steps) - print("ds.multiply(ds1, out=ds)", dt(steps)) - steps.append(timer()) - ds2 = ds.multiply(ds1) - steps.append(timer()) - t2 = dt(steps) - print("ds2 = ds.multiply(ds1)", dt(steps)) - - self.assertLess(t1, t2) - - ds0 = ds - ds0.multiply(2, out=ds0) - steps.append(timer()) - print("ds0.multiply(2,out=ds0)", dt( - steps), 2., ds0.as_array()[0][0][0]) - self.assertEqual(2., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) - ds3 = ds0.multiply(2) - steps.append(timer()) - print("ds3 = ds0.multiply(2)", dt(steps), 4., ds3.as_array()[0][0][0]) - dt2 = dt(steps) - self.assertLess(dt1, dt2) - self.assertEqual(4., ds3.as_array()[0][0][0]) - self.assertEqual(2., ds.as_array()[0][0][0]) - - def binary_divide(self): - print("Test binary divide") - X, Y, Z = 1024, 512, 512 - X, Y, Z = 256, 512, 512 - steps = [timer()] - a = numpy.ones((X, Y, Z), dtype='float32') - steps.append(timer()) - t0 = dt(steps) - - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y', 'Z']) - ds1 = ds.copy() - - steps.append(timer()) - ds.divide(ds1, out=ds) - steps.append(timer()) - t1 = dt(steps) - print("ds.divide(ds1, out=ds)", dt(steps)) - steps.append(timer()) - ds2 = ds.divide(ds1) - steps.append(timer()) - t2 = dt(steps) - print("ds2 = ds.divide(ds1)", dt(steps)) - - self.assertLess(t1, t2) - self.assertEqual(ds.as_array()[0][0][0], 1.) - - ds0 = ds - ds0.divide(2, out=ds0) - steps.append(timer()) - print("ds0.divide(2,out=ds0)", dt(steps), 0.5, ds0.as_array()[0][0][0]) - self.assertEqual(0.5, ds0.as_array()[0][0][0]) - - dt1 = dt(steps) - ds3 = ds0.divide(2) - steps.append(timer()) - print("ds3 = ds0.divide(2)", dt(steps), 0.25, ds3.as_array()[0][0][0]) - dt2 = dt(steps) - self.assertLess(dt1, dt2) - self.assertEqual(.25, ds3.as_array()[0][0][0]) - self.assertEqual(.5, ds.as_array()[0][0][0]) - - def test_creation_copy(self): - shape = (2, 3, 4, 5) - size = shape[0] - for i in range(1, len(shape)): - size = size * shape[i] - #print("a refcount " , sys.getrefcount(a)) - a = numpy.asarray([i for i in range(size)]) - #print("a refcount " , sys.getrefcount(a)) - a = numpy.reshape(a, shape) - #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, True, ['X', 'Y', 'Z', 'W']) - #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a), 2) - - def test_subset(self): - shape = (4, 3, 2) - a = [i for i in range(2*3*4)] - a = numpy.asarray(a) - a = numpy.reshape(a, shape) - ds = DataContainer(a, True, ['X', 'Y', 'Z']) - sub = ds.subset(['X']) - res = True - try: - numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([0, 6, 12, 18])) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - sub = ds.subset(['X'], Y=2, Z=0) - res = True - try: - numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([4, 10, 16, 22])) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - sub = ds.subset(['Y']) - try: - numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([0, 2, 4])) - res = True - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - sub = ds.subset(['Z']) - try: - numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([0, 1])) - res = True - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - sub = ds.subset(['Z'], X=1, Y=2) - try: - numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([10, 11])) - res = True - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - print(a) - sub = ds.subset(['X', 'Y'], Z=1) - res = True - try: - numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([[1, 3, 5], - [7, 9, 11], - [13, 15, 17], - [19, 21, 23]])) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def test_ImageData(self): - # create ImageData from geometry - vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2) - vol = ImageData(geometry=vgeometry) - self.assertEqual(vol.shape, (2, 3, 4)) - - vol1 = vol + 1 - self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape)) - - vol1 = vol - 1 - self.assertNumpyArrayEqual(vol1.as_array(), -numpy.ones(vol.shape)) - - vol1 = 2 * (vol + 1) - self.assertNumpyArrayEqual(vol1.as_array(), 2 * numpy.ones(vol.shape)) - - vol1 = (vol + 1) / 2 - self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) / 2) - - vol1 = vol + 1 - self.assertEqual(vol1.sum(), 2*3*4) - vol1 = (vol + 2) ** 2 - self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) * 4) - - self.assertEqual(vol.number_of_dimensions, 3) - - def test_AcquisitionData(self): - sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), - geom_type='parallel', pixel_num_v=3, - pixel_num_h=5, channels=2) - sino = AcquisitionData(geometry=sgeometry) - self.assertEqual(sino.shape, (2, 10, 3, 5)) - - def assertNumpyArrayEqual(self, first, second): - res = True - try: - numpy.testing.assert_array_equal(first, second) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): - res = True - try: - numpy.testing.assert_array_almost_equal(first, second, decimal) - except AssertionError as err: - res = False - print(err) - print("expected " , second) - print("actual " , first) - - self.assertTrue(res) - def test_DataContainerChaining(self): - dc = self.create_DataContainer(256,256,256,1) - - dc.add(9,out=dc)\ - .subtract(1,out=dc) - self.assertEqual(1+9-1,dc.as_array().flatten()[0]) - def test_reduction(self): - print ("test reductions") - dc = self.create_DataContainer(2,2,2,value=1) - sqnorm = dc.squared_norm() - norm = dc.norm() - self.assertEqual(sqnorm, 8.0) - numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7) - - - -class TestAlgorithms(unittest.TestCase): - def assertNumpyArrayEqual(self, first, second): - res = True - try: - numpy.testing.assert_array_equal(first, second) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): - res = True - try: - numpy.testing.assert_array_almost_equal(first, second, decimal) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def test_FISTA_cvx(self): - if not cvx_not_installable: - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0], 1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt': True} - # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) - g0 = ZeroFun() - - # Initial guess - x_init = DataContainer(np.zeros((n, 1))) - - f.grad(x_init) - - # Run FISTA for least squares plus zero function. - x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt) - - # Print solution and final objective/criterion value for comparison - print("FISTA least squares plus zero function solution and objective value:") - print(x_fista0.array) - print(criter0[-1]) - - # Compare to CVXPY - - # Construct the problem. - x0 = Variable(n) - objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0])) - prob0 = Problem(objective0) - - # The optimal objective is returned by prob.solve(). - result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus zero function solution and objective value:") - print(x0.value) - print(objective0.value) - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fista0.array), x0.value, 6) - else: - self.assertTrue(cvx_not_installable) - - def test_FISTA_Norm1_cvx(self): - if not cvx_not_installable: - opt = {'memopt': True} - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0], 1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt': True} - # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) - g0 = ZeroFun() - - # Initial guess - x_init = DataContainer(np.zeros((n, 1))) - - # Create 1-norm object instance - g1 = Norm1(lam) - - g1(x_init) - g1.prox(x_init, 0.02) - - # Combine with least squares and solve using generic FISTA implementation - x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt) - - # Print for comparison - print("FISTA least squares plus 1-norm solution and objective value:") - print(x_fista1.as_array().squeeze()) - print(criter1[-1]) - - # Compare to CVXPY - - # Construct the problem. - x1 = Variable(n) - objective1 = Minimize( - 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) - prob1 = Problem(objective1) - - # The optimal objective is returned by prob.solve(). - result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1.value) - print(objective1.value) - - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fista1.array), x1.value, 6) - else: - self.assertTrue(cvx_not_installable) - - def skip_test_FBPD_Norm1_cvx(self): - print ("test_FBPD_Norm1_cvx") - if not cvx_not_installable: - opt = {'memopt': True} - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0], 1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt': True} - # Initial guess - x_init = DataContainer(np.random.randn(n, 1)) - - # Create object instances with the test data A and b. - f = Norm2sq(A, b, c=0.5, memopt=True) - f.L = PowerMethodNonsquare(A, 25, x_init)[0] - print ("Lipschitz", f.L) - g0 = ZeroFun() - - - # Create 1-norm object instance - g1 = Norm1(lam) - - # Compare to CVXPY - - # Construct the problem. - x1 = Variable(n) - objective1 = Minimize( - 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) - prob1 = Problem(objective1) - - # The optimal objective is returned by prob.solve(). - result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1.value) - print(objective1.value) - - # Now try another algorithm FBPD for same problem: - x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, - Identity(), None, f, g1) - print(x_fbpd1) - print(criterfbpd1[-1]) - - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fbpd1.array), x1.value, 6) - else: - self.assertTrue(cvx_not_installable) - # Plot criterion curve to see both FISTA and FBPD converge to same value. - # Note that FISTA is very efficient for 1-norm minimization so it beats - # FBPD in this test by a lot. But FBPD can handle a larger class of problems - # than FISTA can. - - # Now try 1-norm and TV denoising with FBPD, first 1-norm. - - # Set up phantom size NxN by creating ImageGeometry, initialising the - # ImageData object with this geometry and empty array and finally put some - # data into its array, and display as image. - def skip_test_FISTA_denoise_cvx(self): - if not cvx_not_installable: - opt = {'memopt': True} - N = 64 - ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) - Phantom = ImageData(geometry=ig) - - x = Phantom.as_array() - - x[int(round(N/4)):int(round(3*N/4)), - int(round(N/4)):int(round(3*N/4))] = 0.5 - x[int(round(N/8)):int(round(7*N/8)), - int(round(3*N/8)):int(round(5*N/8))] = 1 - - # Identity operator for denoising - I = TomoIdentity(ig) - - # Data and add noise - y = I.direct(Phantom) - y.array = y.array + 0.1*np.random.randn(N, N) - - # Data fidelity term - f_denoise = Norm2sq(I, y, c=0.5, memopt=True) - x_init = ImageData(geometry=ig) - f_denoise.L = PowerMethodNonsquare(I, 25, x_init)[0] - - # 1-norm regulariser - lam1_denoise = 1.0 - g1_denoise = Norm1(lam1_denoise) - - # Initial guess - x_init_denoise = ImageData(np.zeros((N, N))) - - # Combine with least squares and solve using generic FISTA implementation - x_fista1_denoise, it1_denoise, timing1_denoise, \ - criter1_denoise = \ - FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) - - print(x_fista1_denoise) - print(criter1_denoise[-1]) - - # Now denoise LS + 1-norm with FBPD - x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise,\ - criterfbpd1_denoise = \ - FBPD(x_init_denoise, I, None, f_denoise, g1_denoise) - print(x_fbpd1_denoise) - print(criterfbpd1_denoise[-1]) - - # Compare to CVXPY - - # Construct the problem. - x1_denoise = Variable(N**2) - objective1_denoise = Minimize( - 0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1)) - prob1_denoise = Problem(objective1_denoise) - - # The optimal objective is returned by prob.solve(). - result1_denoise = prob1_denoise.solve( - verbose=False, solver=SCS, eps=1e-12) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1_denoise.value) - print(objective1_denoise.value) - self.assertNumpyArrayAlmostEqual( - x_fista1_denoise.array.flatten(), x1_denoise.value, 5) - - self.assertNumpyArrayAlmostEqual( - x_fbpd1_denoise.array.flatten(), x1_denoise.value, 5) - x1_cvx = x1_denoise.value - x1_cvx.shape = (N, N) - - # Now TV with FBPD - lam_tv = 0.1 - gtv = TV2D(lam_tv) - gtv(gtv.op.direct(x_init_denoise)) - - opt_tv = {'tol': 1e-4, 'iter': 10000} - - x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise,\ - criterfbpdtv_denoise = \ - FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv, opt=opt_tv) - print(x_fbpdtv_denoise) - print(criterfbpdtv_denoise[-1]) - - # Compare to CVXPY - - # Construct the problem. - xtv_denoise = Variable((N, N)) - objectivetv_denoise = Minimize( - 0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise)) - probtv_denoise = Problem(objectivetv_denoise) - - # The optimal objective is returned by prob.solve(). - resulttv_denoise = probtv_denoise.solve( - verbose=False, solver=SCS, eps=1e-12) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(xtv_denoise.value) - print(objectivetv_denoise.value) - - self.assertNumpyArrayAlmostEqual( - x_fbpdtv_denoise.as_array(), xtv_denoise.value, 1) - - else: - self.assertTrue(cvx_not_installable) - - -class TestFunction(unittest.TestCase): - def assertNumpyArrayEqual(self, first, second): - res = True - try: - numpy.testing.assert_array_equal(first, second) - except AssertionError as err: - res = False - print(err) - self.assertTrue(res) - - def create_simple_ImageData(self): - N = 64 - ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) - Phantom = ImageData(geometry=ig) - - x = Phantom.as_array() - - x[int(round(N/4)):int(round(3*N/4)), - int(round(N/4)):int(round(3*N/4))] = 0.5 - x[int(round(N/8)):int(round(7*N/8)), - int(round(3*N/8)):int(round(5*N/8))] = 1 - - return (ig, Phantom) - - def _test_Norm2(self): - print("test Norm2") - opt = {'memopt': True} - ig, Phantom = self.create_simple_ImageData() - x = Phantom.as_array() - print(Phantom) - print(Phantom.as_array()) - - norm2 = Norm2() - v1 = norm2(x) - v2 = norm2(Phantom) - self.assertEqual(v1, v2) - - p1 = norm2.prox(Phantom, 1) - print(p1) - p2 = norm2.prox(x, 1) - self.assertNumpyArrayEqual(p1.as_array(), p2) - - p3 = norm2.proximal(Phantom, 1) - p4 = norm2.proximal(x, 1) - self.assertNumpyArrayEqual(p3.as_array(), p2) - self.assertNumpyArrayEqual(p3.as_array(), p4) - - out = Phantom.copy() - p5 = norm2.proximal(Phantom, 1, out=out) - self.assertEqual(id(p5), id(out)) - self.assertNumpyArrayEqual(p5.as_array(), p3.as_array()) -# ============================================================================= -# def test_upper(self): -# self.assertEqual('foo'.upper(), 'FOO') -# -# def test_isupper(self): -# self.assertTrue('FOO'.isupper()) -# self.assertFalse('Foo'.isupper()) -# -# def test_split(self): -# s = 'hello world' -# self.assertEqual(s.split(), ['hello', 'world']) -# # check that s.split fails when the separator is not a string -# with self.assertRaises(TypeError): -# s.split(2) -# ============================================================================= - -class TestNexusReader(unittest.TestCase): - - def setUp(self): - wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') - self.filename = '24737_fd.nxs' - - def tearDown(self): - os.remove(self.filename) - - - def testGetDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct") - - def testGetProjectionDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct") - - def testLoadProjectionWithoutDimensions(self): - nr = NexusReader(self.filename) - projections = nr.load_projection() - self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") - - def testLoadProjectionWithDimensions(self): - nr = NexusReader(self.filename) - projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) - self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") - - def testLoadProjectionCompareSingle(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) - - def testLoadProjectionCompareMulti(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) - - def testLoadProjectionCompareRandom(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20))) - numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) - - def testLoadProjectionCompareFull(self): - nr = NexusReader(self.filename) - projections_full = nr.load_projection() - projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) - - def testLoadFlatCompareFull(self): - nr = NexusReader(self.filename) - flats_full = nr.load_flat() - flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) - - def testLoadDarkCompareFull(self): - nr = NexusReader(self.filename) - darks_full = nr.load_dark() - darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) - - def testProjectionAngles(self): - nr = NexusReader(self.filename) - angles = nr.get_projection_angles() - self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") - - def test_get_acquisition_data_subset(self): - nr = NexusReader(self.filename) - key = nr.get_image_keys() - sl = nr.get_acquisition_data_subset(0,10) - data = nr.get_acquisition_data().subset(['vertical','horizontal']) - - self.assertTrue(sl.shape , (10,data.shape[1])) - - -if __name__ == '__main__': - unittest.main() - diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index b584344..eaf124b 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -23,12 +23,16 @@ import os import sys -cil_version='0.11.3' +cil_version=os.environ['CIL_VERSION'] +if cil_version == '': + print("Please set the environmental variable CIL_VERSION") + sys.exit(1) setup( name="ccpi-framework", version=cil_version, - packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation'], + packages=['ccpi' , 'ccpi.io', 'ccpi.optimisation', + 'ccpi.optimisation.algorithms'], # Project uses reStructuredText, so ensure that the docutils get # installed or upgraded on the target machine diff --git a/Wrappers/Python/test/TestReader.py b/Wrappers/Python/test/TestReader.py deleted file mode 100644 index 51db052..0000000 --- a/Wrappers/Python/test/TestReader.py +++ /dev/null @@ -1,134 +0,0 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev, Edoardo Pasca and Srikanth Nagella - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -''' -Unit tests for Readers - -@author: Mr. Srikanth Nagella -''' -import unittest - -import numpy.testing -import wget -import os -from ccpi.io.reader import NexusReader -from ccpi.io.reader import XTEKReader -#@unittest.skip -class TestNexusReader(unittest.TestCase): - - def setUp(self): - wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') - self.filename = '24737_fd.nxs' - - def tearDown(self): - os.remove(self.filename) - - - def testGetDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.getSinogramDimensions(), (135, 91, 160), "Sinogram dimensions are not correct") - - def testGetProjectionDimensions(self): - nr = NexusReader(self.filename) - self.assertEqual(nr.getProjectionDimensions(), (91,135,160), "Projection dimensions are not correct") - - def testLoadProjectionWithoutDimensions(self): - nr = NexusReader(self.filename) - projections = nr.loadProjection() - self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") - - def testLoadProjectionWithDimensions(self): - nr = NexusReader(self.filename) - projections = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) - self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") - - def testLoadProjectionCompareSingle(self): - nr = NexusReader(self.filename) - projections_full = nr.loadProjection() - projections_part = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) - - def testLoadProjectionCompareMulti(self): - nr = NexusReader(self.filename) - projections_full = nr.loadProjection() - projections_part = nr.loadProjection((slice(0,3), slice(0,135), slice(0,160))) - numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) - - def testLoadProjectionCompareRandom(self): - nr = NexusReader(self.filename) - projections_full = nr.loadProjection() - projections_part = nr.loadProjection((slice(1,8), slice(5,10), slice(8,20))) - numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) - - def testLoadProjectionCompareFull(self): - nr = NexusReader(self.filename) - projections_full = nr.loadProjection() - projections_part = nr.loadProjection((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) - - def testLoadFlatCompareFull(self): - nr = NexusReader(self.filename) - flats_full = nr.loadFlat() - flats_part = nr.loadFlat((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) - - def testLoadDarkCompareFull(self): - nr = NexusReader(self.filename) - darks_full = nr.loadDark() - darks_part = nr.loadDark((slice(None,None), slice(None,None), slice(None,None))) - numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) - - def testProjectionAngles(self): - nr = NexusReader(self.filename) - angles = nr.getProjectionAngles() - self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") - -class TestXTEKReader(unittest.TestCase): - - def setUp(self): - testpath, filename = os.path.split(os.path.realpath(__file__)) - testpath = os.path.normpath(os.path.join(testpath, '..','..','..')) - self.filename = os.path.join(testpath,'data','SophiaBeads','SophiaBeads_64_averaged.xtekct') - - def tearDown(self): - pass - - def testLoad(self): - xtekReader = XTEKReader(self.filename) - self.assertEqual(xtekReader.geometry.pixel_num_h, 500, "Detector pixel X size is not correct") - self.assertEqual(xtekReader.geometry.pixel_num_v, 500, "Detector pixel Y size is not correct") - self.assertEqual(xtekReader.geometry.dist_source_center, -80.6392412185669, "Distance from source to center is not correct") - self.assertEqual(xtekReader.geometry.dist_center_detector, (1007.006 - 80.6392412185669), "Distance from center to detector is not correct") - - def testReadAngles(self): - xtekReader = XTEKReader(self.filename) - angles = xtekReader.readAngles() - self.assertEqual(angles.shape, (63,), "Angles doesn't match") - self.assertAlmostEqual(angles[46], -085.717, 3, "46th Angles doesn't match") - - def testLoadProjection(self): - xtekReader = XTEKReader(self.filename) - pixels = xtekReader.loadProjection() - self.assertEqual(pixels.shape, (63, 500, 500), "projections shape doesn't match") - - - -if __name__ == "__main__": - #import sys;sys.argv = ['', 'TestNexusReader.testLoad'] - unittest.main() \ No newline at end of file diff --git a/Wrappers/Python/test/__init__.py b/Wrappers/Python/test/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/Wrappers/Python/test/skip_TestReader.py b/Wrappers/Python/test/skip_TestReader.py new file mode 100644 index 0000000..ec7ed58 --- /dev/null +++ b/Wrappers/Python/test/skip_TestReader.py @@ -0,0 +1,134 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev, Edoardo Pasca and Srikanth Nagella + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +''' +Unit tests for Readers + +@author: Mr. Srikanth Nagella +''' +import unittest + +import numpy.testing +import wget +import os +from ccpi.io.reader import NexusReader +from ccpi.io.reader import XTEKReader +#@unittest.skip +class TestNexusReader(unittest.TestCase): + + def setUp(self): + wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + self.filename = '24737_fd.nxs' + + def tearDown(self): + os.remove(self.filename) + + + def testGetDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.getSinogramDimensions(), (135, 91, 160), "Sinogram dimensions are not correct") + + def testGetProjectionDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.getProjectionDimensions(), (91,135,160), "Projection dimensions are not correct") + + def testLoadProjectionWithoutDimensions(self): + nr = NexusReader(self.filename) + projections = nr.loadProjection() + self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionWithDimensions(self): + nr = NexusReader(self.filename) + projections = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) + self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionCompareSingle(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(0,1), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) + + def testLoadProjectionCompareMulti(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(0,3), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) + + def testLoadProjectionCompareRandom(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(1,8), slice(5,10), slice(8,20))) + numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) + + def testLoadProjectionCompareFull(self): + nr = NexusReader(self.filename) + projections_full = nr.loadProjection() + projections_part = nr.loadProjection((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) + + def testLoadFlatCompareFull(self): + nr = NexusReader(self.filename) + flats_full = nr.loadFlat() + flats_part = nr.loadFlat((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) + + def testLoadDarkCompareFull(self): + nr = NexusReader(self.filename) + darks_full = nr.loadDark() + darks_part = nr.loadDark((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) + + def testProjectionAngles(self): + nr = NexusReader(self.filename) + angles = nr.getProjectionAngles() + self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") + +class TestXTEKReader(unittest.TestCase): + + def setUp(self): + testpath, filename = os.path.split(os.path.realpath(__file__)) + testpath = os.path.normpath(os.path.join(testpath, '..','..','..')) + self.filename = os.path.join(testpath,'data','SophiaBeads','SophiaBeads_64_averaged.xtekct') + + def tearDown(self): + pass + + def testLoad(self): + xtekReader = XTEKReader(self.filename) + self.assertEqual(xtekReader.geometry.pixel_num_h, 500, "Detector pixel X size is not correct") + self.assertEqual(xtekReader.geometry.pixel_num_v, 500, "Detector pixel Y size is not correct") + self.assertEqual(xtekReader.geometry.dist_source_center, -80.6392412185669, "Distance from source to center is not correct") + self.assertEqual(xtekReader.geometry.dist_center_detector, (1007.006 - 80.6392412185669), "Distance from center to detector is not correct") + + def testReadAngles(self): + xtekReader = XTEKReader(self.filename) + angles = xtekReader.readAngles() + self.assertEqual(angles.shape, (63,), "Angles doesn't match") + self.assertAlmostEqual(angles[46], -085.717, 3, "46th Angles doesn't match") + + def testLoadProjection(self): + xtekReader = XTEKReader(self.filename) + pixels = xtekReader.loadProjection() + self.assertEqual(pixels.shape, (63, 500, 500), "projections shape doesn't match") + + + +if __name__ == "__main__": + #import sys;sys.argv = ['', 'TestNexusReader.testLoad'] + unittest.main() \ No newline at end of file diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py new file mode 100755 index 0000000..05f3fe8 --- /dev/null +++ b/Wrappers/Python/test/test_DataContainer.py @@ -0,0 +1,499 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:08:21 2019 + +@author: ofn77899 +""" +import sys +import unittest +import numpy +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from timeit import default_timer as timer + + +def dt(steps): + return steps[-1] - steps[-2] +def aid(x): + # This function returns the memory + # block address of an array. + return x.__array_interface__['data'][0] + + + +class TestDataContainer(unittest.TestCase): + def create_simple_ImageData(self): + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + return (ig, Phantom) + + def create_DataContainer(self, X,Y,Z, value=1): + steps = [timer()] + a = value * numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + return ds + + def test_creation_nocopy(self): + shape = (2, 3, 4, 5) + size = shape[0] + for i in range(1, len(shape)): + size = size * shape[i] + #print("a refcount " , sys.getrefcount(a)) + a = numpy.asarray([i for i in range(size)]) + #print("a refcount " , sys.getrefcount(a)) + a = numpy.reshape(a, shape) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z', 'W']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 3) + self.assertEqual(ds.dimension_labels, {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'}) + + def testGb_creation_nocopy(self): + X, Y, Z = 512, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print("test clone") + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 3) + ds1 = ds.copy() + self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) + ds1 = ds.clone() + self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) + + def testInlineAlgebra(self): + print("Test Inline Algebra") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print(t0) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + #ds.__iadd__( 2 ) + ds += 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 3.) + #ds.__isub__( 2 ) + ds -= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + #ds.__imul__( 2 ) + ds *= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 2.) + #ds.__idiv__( 2 ) + ds /= 2 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds1 = ds.copy() + #ds1.__iadd__( 1 ) + ds1 += 1 + #ds.__iadd__( ds1 ) + ds += ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 3.) + #ds.__isub__( ds1 ) + ds -= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + #ds.__imul__( ds1 ) + ds *= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 2.) + #ds.__idiv__( ds1 ) + ds /= ds1 + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + def test_unary_operations(self): + print("Test unary operations") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = -numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + print(t0) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + + ds.sign(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], -1.) + + ds.abs(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds.__imul__(2) + ds.sqrt(out=ds) + steps.append(timer()) + print(dt(steps)) + self.assertEqual(ds.as_array()[0][0][0], + numpy.sqrt(2., dtype='float32')) + + def test_binary_operations(self): + self.binary_add() + self.binary_subtract() + self.binary_multiply() + self.binary_divide() + + def binary_add(self): + print("Test binary add") + X, Y, Z = 512, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.add(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.add(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.add(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.add(ds1)", dt(steps)) + + self.assertLess(t1, t2) + self.assertEqual(ds.as_array()[0][0][0], 2.) + + ds0 = ds + ds0.add(2, out=ds0) + steps.append(timer()) + print("ds0.add(2,out=ds0)", dt(steps), 3, ds0.as_array()[0][0][0]) + self.assertEqual(4., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.add(2) + steps.append(timer()) + print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + + def binary_subtract(self): + print("Test binary subtract") + X, Y, Z = 512, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.subtract(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.subtract(ds1, out=ds)", dt(steps)) + self.assertEqual(0., ds.as_array()[0][0][0]) + + steps.append(timer()) + ds2 = ds.subtract(ds1) + self.assertEqual(-1., ds2.as_array()[0][0][0]) + + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.subtract(ds1)", dt(steps)) + + self.assertLess(t1, t2) + + del ds1 + ds0 = ds.copy() + steps.append(timer()) + ds0.subtract(2, out=ds0) + #ds0.__isub__( 2 ) + steps.append(timer()) + print("ds0.subtract(2,out=ds0)", dt( + steps), -2., ds0.as_array()[0][0][0]) + self.assertEqual(-2., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.subtract(2) + steps.append(timer()) + print("ds3 = ds0.subtract(2)", dt(steps), 0., ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(-2., ds0.as_array()[0][0][0]) + self.assertEqual(-4., ds3.as_array()[0][0][0]) + + def binary_multiply(self): + print("Test binary multiply") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.multiply(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.multiply(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.multiply(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.multiply(ds1)", dt(steps)) + + self.assertLess(t1, t2) + + ds0 = ds + ds0.multiply(2, out=ds0) + steps.append(timer()) + print("ds0.multiply(2,out=ds0)", dt( + steps), 2., ds0.as_array()[0][0][0]) + self.assertEqual(2., ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.multiply(2) + steps.append(timer()) + print("ds3 = ds0.multiply(2)", dt(steps), 4., ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(4., ds3.as_array()[0][0][0]) + self.assertEqual(2., ds.as_array()[0][0][0]) + + def binary_divide(self): + print("Test binary divide") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 + steps = [timer()] + a = numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds1 = ds.copy() + + steps.append(timer()) + ds.divide(ds1, out=ds) + steps.append(timer()) + t1 = dt(steps) + print("ds.divide(ds1, out=ds)", dt(steps)) + steps.append(timer()) + ds2 = ds.divide(ds1) + steps.append(timer()) + t2 = dt(steps) + print("ds2 = ds.divide(ds1)", dt(steps)) + + self.assertLess(t1, t2) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds0 = ds + ds0.divide(2, out=ds0) + steps.append(timer()) + print("ds0.divide(2,out=ds0)", dt(steps), 0.5, ds0.as_array()[0][0][0]) + self.assertEqual(0.5, ds0.as_array()[0][0][0]) + + dt1 = dt(steps) + ds3 = ds0.divide(2) + steps.append(timer()) + print("ds3 = ds0.divide(2)", dt(steps), 0.25, ds3.as_array()[0][0][0]) + dt2 = dt(steps) + self.assertLess(dt1, dt2) + self.assertEqual(.25, ds3.as_array()[0][0][0]) + self.assertEqual(.5, ds.as_array()[0][0][0]) + + def test_creation_copy(self): + shape = (2, 3, 4, 5) + size = shape[0] + for i in range(1, len(shape)): + size = size * shape[i] + #print("a refcount " , sys.getrefcount(a)) + a = numpy.asarray([i for i in range(size)]) + #print("a refcount " , sys.getrefcount(a)) + a = numpy.reshape(a, shape) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, True, ['X', 'Y', 'Z', 'W']) + #print("a refcount " , sys.getrefcount(a)) + self.assertEqual(sys.getrefcount(a), 2) + + def test_subset(self): + shape = (4, 3, 2) + a = [i for i in range(2*3*4)] + a = numpy.asarray(a) + a = numpy.reshape(a, shape) + ds = DataContainer(a, True, ['X', 'Y', 'Z']) + sub = ds.subset(['X']) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([0, 6, 12, 18])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['X'], Y=2, Z=0) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([4, 10, 16, 22])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['Y']) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([0, 2, 4])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + sub = ds.subset(['Z']) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([0, 1])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + sub = ds.subset(['Z'], X=1, Y=2) + try: + numpy.testing.assert_array_equal( + sub.as_array(), numpy.asarray([10, 11])) + res = True + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + print(a) + sub = ds.subset(['X', 'Y'], Z=1) + res = True + try: + numpy.testing.assert_array_equal(sub.as_array(), + numpy.asarray([[1, 3, 5], + [7, 9, 11], + [13, 15, 17], + [19, 21, 23]])) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def test_ImageData(self): + # create ImageData from geometry + vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2) + vol = ImageData(geometry=vgeometry) + self.assertEqual(vol.shape, (2, 3, 4)) + + vol1 = vol + 1 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape)) + + vol1 = vol - 1 + self.assertNumpyArrayEqual(vol1.as_array(), -numpy.ones(vol.shape)) + + vol1 = 2 * (vol + 1) + self.assertNumpyArrayEqual(vol1.as_array(), 2 * numpy.ones(vol.shape)) + + vol1 = (vol + 1) / 2 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) / 2) + + vol1 = vol + 1 + self.assertEqual(vol1.sum(), 2*3*4) + vol1 = (vol + 2) ** 2 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) * 4) + + self.assertEqual(vol.number_of_dimensions, 3) + + def test_AcquisitionData(self): + sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), + geom_type='parallel', pixel_num_v=3, + pixel_num_h=5, channels=2) + sino = AcquisitionData(geometry=sgeometry) + self.assertEqual(sino.shape, (2, 10, 3, 5)) + + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + print("expected " , second) + print("actual " , first) + + self.assertTrue(res) + def test_DataContainerChaining(self): + dc = self.create_DataContainer(256,256,256,1) + + dc.add(9,out=dc)\ + .subtract(1,out=dc) + self.assertEqual(1+9-1,dc.as_array().flatten()[0]) + def test_reduction(self): + print ("test reductions") + dc = self.create_DataContainer(2,2,2,value=1) + sqnorm = dc.squared_norm() + norm = dc.norm() + self.assertEqual(sqnorm, 8.0) + numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7) + + + +if __name__ == '__main__': + unittest.main() + \ No newline at end of file diff --git a/Wrappers/Python/test/test_NexusReader.py b/Wrappers/Python/test/test_NexusReader.py new file mode 100755 index 0000000..55543ba --- /dev/null +++ b/Wrappers/Python/test/test_NexusReader.py @@ -0,0 +1,96 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:05:00 2019 + +@author: ofn77899 +""" + +import unittest +import wget +import os +from ccpi.io.reader import NexusReader +import numpy + + +class TestNexusReader(unittest.TestCase): + + def setUp(self): + wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + self.filename = '24737_fd.nxs' + + def tearDown(self): + os.remove(self.filename) + + + def testGetDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct") + + def testGetProjectionDimensions(self): + nr = NexusReader(self.filename) + self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct") + + def testLoadProjectionWithoutDimensions(self): + nr = NexusReader(self.filename) + projections = nr.load_projection() + self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionWithDimensions(self): + nr = NexusReader(self.filename) + projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) + self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct") + + def testLoadProjectionCompareSingle(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:]) + + def testLoadProjectionCompareMulti(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160))) + numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:]) + + def testLoadProjectionCompareRandom(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20))) + numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20]) + + def testLoadProjectionCompareFull(self): + nr = NexusReader(self.filename) + projections_full = nr.load_projection() + projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:]) + + def testLoadFlatCompareFull(self): + nr = NexusReader(self.filename) + flats_full = nr.load_flat() + flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:]) + + def testLoadDarkCompareFull(self): + nr = NexusReader(self.filename) + darks_full = nr.load_dark() + darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None))) + numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:]) + + def testProjectionAngles(self): + nr = NexusReader(self.filename) + angles = nr.get_projection_angles() + self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct") + + def test_get_acquisition_data_subset(self): + nr = NexusReader(self.filename) + key = nr.get_image_keys() + sl = nr.get_acquisition_data_subset(0,10) + data = nr.get_acquisition_data().subset(['vertical','horizontal']) + + self.assertTrue(sl.shape , (10,data.shape[1])) + + + +if __name__ == '__main__': + unittest.main() + diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py new file mode 100755 index 0000000..b5959b5 --- /dev/null +++ b/Wrappers/Python/test/test_algorithms.py @@ -0,0 +1,123 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 27 15:11:36 2019 + +@author: ofn77899 +""" + +import unittest +import numpy +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from ccpi.optimisation.ops import TomoIdentity +from ccpi.optimisation.funcs import Norm2sq +from ccpi.optimisation.algorithms import GradientDescent +from ccpi.optimisation.algorithms import CGLS +from ccpi.optimisation.algorithms import FISTA +from ccpi.optimisation.algorithms import FBPD + + + + +class TestAlgorithms(unittest.TestCase): + def setUp(self): + #wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs') + #self.filename = '24737_fd.nxs' + # we use TomoIdentity as the operator and solve the simple least squares + # problem for a random-valued ImageData or AcquisitionData b? + # Then we know the minimiser is b itself + + # || I x -b ||^2 + + # create an ImageGeometry + ig = ImageGeometry(12,13,14) + pass + + def tearDown(self): + #os.remove(self.filename) + pass + + def test_GradientDescent(self): + print ("Test GradientDescent") + ig = ImageGeometry(12,13,14) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + norm2sq = Norm2sq(identity, b) + + alg = GradientDescent(x_init=x_init, + objective_function=norm2sq, + rate=0.3) + alg.max_iteration = 20 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + def test_CGLS(self): + print ("Test CGLS") + ig = ImageGeometry(124,153,154) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + alg = CGLS(x_init=x_init, operator=identity, data=b) + alg.max_iteration = 1 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + + def test_FISTA(self): + print ("Test FISTA") + ig = ImageGeometry(127,139,149) + x_init = ImageData(geometry=ig) + b = x_init.copy() + # fill with random numbers + b.fill(numpy.random.random(x_init.shape)) + x_init = ImageData(geometry=ig) + x_init.fill(numpy.random.random(x_init.shape)) + + identity = TomoIdentity(geometry=ig) + + norm2sq = Norm2sq(identity, b) + opt = {'tol': 1e-4, 'memopt':False} + alg = FISTA(x_init=x_init, f=norm2sq, g=None, opt=opt) + alg.max_iteration = 2 + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + alg.run(20, verbose=True) + self.assertNumpyArrayAlmostEqual(alg.x.as_array(), b.as_array()) + + + + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + + + + +if __name__ == '__main__': + unittest.main() + \ No newline at end of file diff --git a/Wrappers/Python/test/test_run_test.py b/Wrappers/Python/test/test_run_test.py new file mode 100755 index 0000000..d0b87f5 --- /dev/null +++ b/Wrappers/Python/test/test_run_test.py @@ -0,0 +1,435 @@ +import unittest +import numpy +import numpy as np +from ccpi.framework import DataContainer +from ccpi.framework import ImageData +from ccpi.framework import AcquisitionData +from ccpi.framework import ImageGeometry +from ccpi.framework import AcquisitionGeometry +from ccpi.optimisation.algs import FISTA +from ccpi.optimisation.algs import FBPD +from ccpi.optimisation.funcs import Norm2sq +from ccpi.optimisation.funcs import ZeroFun +from ccpi.optimisation.funcs import Norm1 +from ccpi.optimisation.funcs import TV2D +from ccpi.optimisation.funcs import Norm2 + +from ccpi.optimisation.ops import LinearOperatorMatrix +from ccpi.optimisation.ops import TomoIdentity +from ccpi.optimisation.ops import Identity +from ccpi.optimisation.ops import PowerMethodNonsquare + + +import numpy.testing + +try: + from cvxpy import * + cvx_not_installable = False +except ImportError: + cvx_not_installable = True + + +def aid(x): + # This function returns the memory + # block address of an array. + return x.__array_interface__['data'][0] + + +def dt(steps): + return steps[-1] - steps[-2] + + + + +class TestAlgorithms(unittest.TestCase): + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): + res = True + try: + numpy.testing.assert_array_almost_equal(first, second, decimal) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def test_FISTA_cvx(self): + if not cvx_not_installable: + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + f.grad(x_init) + + # Run FISTA for least squares plus zero function. + x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt) + + # Print solution and final objective/criterion value for comparison + print("FISTA least squares plus zero function solution and objective value:") + print(x_fista0.array) + print(criter0[-1]) + + # Compare to CVXPY + + # Construct the problem. + x0 = Variable(n) + objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0])) + prob0 = Problem(objective0) + + # The optimal objective is returned by prob.solve(). + result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus zero function solution and objective value:") + print(x0.value) + print(objective0.value) + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fista0.array), x0.value, 6) + else: + self.assertTrue(cvx_not_installable) + + def test_FISTA_Norm1_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + # Create 1-norm object instance + g1 = Norm1(lam) + + g1(x_init) + g1.prox(x_init, 0.02) + + # Combine with least squares and solve using generic FISTA implementation + x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt) + + # Print for comparison + print("FISTA least squares plus 1-norm solution and objective value:") + print(x_fista1.as_array().squeeze()) + print(criter1[-1]) + + # Compare to CVXPY + + # Construct the problem. + x1 = Variable(n) + objective1 = Minimize( + 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) + prob1 = Problem(objective1) + + # The optimal objective is returned by prob.solve(). + result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1.value) + print(objective1.value) + + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fista1.array), x1.value, 6) + else: + self.assertTrue(cvx_not_installable) + + def skip_test_FBPD_Norm1_cvx(self): + print ("test_FBPD_Norm1_cvx") + if not cvx_not_installable: + opt = {'memopt': True} + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Initial guess + x_init = DataContainer(np.random.randn(n, 1)) + + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + f.L = PowerMethodNonsquare(A, 25, x_init)[0] + print ("Lipschitz", f.L) + g0 = ZeroFun() + + + # Create 1-norm object instance + g1 = Norm1(lam) + + # Compare to CVXPY + + # Construct the problem. + x1 = Variable(n) + objective1 = Minimize( + 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) + prob1 = Problem(objective1) + + # The optimal objective is returned by prob.solve(). + result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1.value) + print(objective1.value) + + # Now try another algorithm FBPD for same problem: + x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, + Identity(), None, f, g1) + print(x_fbpd1) + print(criterfbpd1[-1]) + + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fbpd1.array), x1.value, 6) + else: + self.assertTrue(cvx_not_installable) + # Plot criterion curve to see both FISTA and FBPD converge to same value. + # Note that FISTA is very efficient for 1-norm minimization so it beats + # FBPD in this test by a lot. But FBPD can handle a larger class of problems + # than FISTA can. + + # Now try 1-norm and TV denoising with FBPD, first 1-norm. + + # Set up phantom size NxN by creating ImageGeometry, initialising the + # ImageData object with this geometry and empty array and finally put some + # data into its array, and display as image. + def skip_test_FISTA_denoise_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + # Identity operator for denoising + I = TomoIdentity(ig) + + # Data and add noise + y = I.direct(Phantom) + y.array = y.array + 0.1*np.random.randn(N, N) + + # Data fidelity term + f_denoise = Norm2sq(I, y, c=0.5, memopt=True) + x_init = ImageData(geometry=ig) + f_denoise.L = PowerMethodNonsquare(I, 25, x_init)[0] + + # 1-norm regulariser + lam1_denoise = 1.0 + g1_denoise = Norm1(lam1_denoise) + + # Initial guess + x_init_denoise = ImageData(np.zeros((N, N))) + + # Combine with least squares and solve using generic FISTA implementation + x_fista1_denoise, it1_denoise, timing1_denoise, \ + criter1_denoise = \ + FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) + + print(x_fista1_denoise) + print(criter1_denoise[-1]) + + # Now denoise LS + 1-norm with FBPD + x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise,\ + criterfbpd1_denoise = \ + FBPD(x_init_denoise, I, None, f_denoise, g1_denoise) + print(x_fbpd1_denoise) + print(criterfbpd1_denoise[-1]) + + # Compare to CVXPY + + # Construct the problem. + x1_denoise = Variable(N**2) + objective1_denoise = Minimize( + 0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1)) + prob1_denoise = Problem(objective1_denoise) + + # The optimal objective is returned by prob.solve(). + result1_denoise = prob1_denoise.solve( + verbose=False, solver=SCS, eps=1e-12) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1_denoise.value) + print(objective1_denoise.value) + self.assertNumpyArrayAlmostEqual( + x_fista1_denoise.array.flatten(), x1_denoise.value, 5) + + self.assertNumpyArrayAlmostEqual( + x_fbpd1_denoise.array.flatten(), x1_denoise.value, 5) + x1_cvx = x1_denoise.value + x1_cvx.shape = (N, N) + + # Now TV with FBPD + lam_tv = 0.1 + gtv = TV2D(lam_tv) + gtv(gtv.op.direct(x_init_denoise)) + + opt_tv = {'tol': 1e-4, 'iter': 10000} + + x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise,\ + criterfbpdtv_denoise = \ + FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv, opt=opt_tv) + print(x_fbpdtv_denoise) + print(criterfbpdtv_denoise[-1]) + + # Compare to CVXPY + + # Construct the problem. + xtv_denoise = Variable((N, N)) + objectivetv_denoise = Minimize( + 0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise)) + probtv_denoise = Problem(objectivetv_denoise) + + # The optimal objective is returned by prob.solve(). + resulttv_denoise = probtv_denoise.solve( + verbose=False, solver=SCS, eps=1e-12) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(xtv_denoise.value) + print(objectivetv_denoise.value) + + self.assertNumpyArrayAlmostEqual( + x_fbpdtv_denoise.as_array(), xtv_denoise.value, 1) + + else: + self.assertTrue(cvx_not_installable) + + +class TestFunction(unittest.TestCase): + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def create_simple_ImageData(self): + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + return (ig, Phantom) + + def _test_Norm2(self): + print("test Norm2") + opt = {'memopt': True} + ig, Phantom = self.create_simple_ImageData() + x = Phantom.as_array() + print(Phantom) + print(Phantom.as_array()) + + norm2 = Norm2() + v1 = norm2(x) + v2 = norm2(Phantom) + self.assertEqual(v1, v2) + + p1 = norm2.prox(Phantom, 1) + print(p1) + p2 = norm2.prox(x, 1) + self.assertNumpyArrayEqual(p1.as_array(), p2) + + p3 = norm2.proximal(Phantom, 1) + p4 = norm2.proximal(x, 1) + self.assertNumpyArrayEqual(p3.as_array(), p2) + self.assertNumpyArrayEqual(p3.as_array(), p4) + + out = Phantom.copy() + p5 = norm2.proximal(Phantom, 1, out=out) + self.assertEqual(id(p5), id(out)) + self.assertNumpyArrayEqual(p5.as_array(), p3.as_array()) +# ============================================================================= +# def test_upper(self): +# self.assertEqual('foo'.upper(), 'FOO') +# +# def test_isupper(self): +# self.assertTrue('FOO'.isupper()) +# self.assertFalse('Foo'.isupper()) +# +# def test_split(self): +# s = 'hello world' +# self.assertEqual(s.split(), ['hello', 'world']) +# # check that s.split fails when the separator is not a string +# with self.assertRaises(TypeError): +# s.split(2) +# ============================================================================= + + + +if __name__ == '__main__': + unittest.main() + -- cgit v1.2.3