Newalgo (#22)

* Use new algos, remove FBPD

* Update mc demo to new algs.

* Fix bug overwriting b

* Move work_out_ajoint files

* Update sophiabeads 2d demo to new reader

* Update sophiabeads 3d demo

* Update to generic data path

* Demo read only required part of Nikon data

* Update nexus demo, still missing new reader.

* Nexus demo new loader. Data path still hardcoded.

8 files changed, 325 insertions, 594 deletions

diff --git a/Wrappers/Python/wip/demo_astra_mc.py b/Wrappers/Python/wip/demo_astra_mc.py
index e5999de..fde0120 100755
--- a/Wrappers/Python/wip/demo_astra_mc.py
+++ b/Wrappers/Python/wip/demo_astra_mc.py
@@ -5,9 +5,10 @@
 # regularisation reconstructions.
 
 # Do all imports
-from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry
-from ccpi.optimisation.algs import FISTA
-from ccpi.optimisation.funcs import Norm2sq, Norm1
+from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, \
+                           AcquisitionGeometry
+from ccpi.optimisation.algorithms import FISTA
+from ccpi.optimisation.functions import Norm2Sq, L1Norm
 from ccpi.astra.operators import AstraProjectorMC
 
 import numpy
@@ -97,23 +98,28 @@ plt.show()
 # Using the test data b, different reconstruction methods can now be set up as
 # demonstrated in the rest of this file. In general all methods need an initial 
 # guess and some algorithm options to be set:
-x_init = ImageData(numpy.zeros(x.shape),geometry=ig)
+x_init = ig.allocate(0.0)
 opt = {'tol': 1e-4, 'iter': 200}
 
 # 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(Aop,b,c=0.5)
+f = Norm2Sq(Aop,b,c=0.5)
 
 # Run FISTA for least squares without regularization
-x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
+FISTA_alg = FISTA()
+FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
+FISTA_alg.max_iteration = 2000
+FISTA_alg.run(opt['iter'])
+x_FISTA = FISTA_alg.get_output()
 
-# Display reconstruction and criteration
+# Display reconstruction and criterion
 ff0, axarrf0 = plt.subplots(1,numchannels)
 for k in numpy.arange(3):
-    axarrf0[k].imshow(x_fista0.as_array()[k],vmin=0,vmax=2.5)
+    axarrf0[k].imshow(x_FISTA.as_array()[k],vmin=0,vmax=2.5)
 plt.show()
 
-plt.semilogy(criter0)
+plt.figure()
+plt.semilogy(FISTA_alg.objective)
 plt.title('Criterion vs iterations, least squares')
 plt.show()
 
@@ -121,17 +127,22 @@ plt.show()
 # such as 1-norm regularisation with choice of regularisation parameter lam. 
 # Again the regulariser is over all channels:
 lam = 10
-g0 = Norm1(lam)
+g1 = lam * L1Norm()
 
-# Run FISTA for least squares plus 1-norm function.
-x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
+# Run FISTA for least squares plus 1-norm regularisation.
+FISTA_alg1 = FISTA()
+FISTA_alg1.set_up(x_init=x_init, f=f, g=g1, opt=opt)
+FISTA_alg1.max_iteration = 2000
+FISTA_alg1.run(opt['iter'])
+x_FISTA1 = FISTA_alg1.get_output()
 
-# Display reconstruction and criteration
+# Display reconstruction and criterion
 ff1, axarrf1 = plt.subplots(1,numchannels)
 for k in numpy.arange(3):
-    axarrf1[k].imshow(x_fista1.as_array()[k],vmin=0,vmax=2.5)
+    axarrf1[k].imshow(x_FISTA1.as_array()[k],vmin=0,vmax=2.5)
 plt.show()
 
-plt.semilogy(criter1)
+plt.figure()
+plt.semilogy(FISTA_alg1.objective)
 plt.title('Criterion vs iterations, least squares plus 1-norm regu')
 plt.show()
diff --git a/Wrappers/Python/wip/demo_astra_nexus.py b/Wrappers/Python/wip/demo_astra_nexus.py
index 7892d24..64895dd 100644
--- a/Wrappers/Python/wip/demo_astra_nexus.py
+++ b/Wrappers/Python/wip/demo_astra_nexus.py
@@ -1,6 +1,6 @@
 
 # This script demonstrates how to load a parallel beam data set in Nexus 
-# format, apply dark and flat field correction and reconstruct using the
+# format and reconstruct using the
 # modular optimisation framework and the ASTRA Tomography toolbox.
 # 
 # The data set is available from
@@ -8,67 +8,46 @@
 # and should be downloaded to a local directory to be specified below.
 
 # All own imports
-from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry
-from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-from ccpi.optimisation.funcs import Norm2sq, Norm1
-from ccpi.plugins.ops import CCPiProjectorSimple
-from ccpi.processors import Normalizer, CenterOfRotationFinder 
-from ccpi.plugins.processors import AcquisitionDataPadder
-from ccpi.io.reader import NexusReader
+from ccpi.framework import ImageData, ImageGeometry
+from ccpi.optimisation.algorithms import CGLS, FISTA
+from ccpi.optimisation.functions import Norm2Sq, L1Norm
+from ccpi.processors import CenterOfRotationFinder 
+from ccpi.io import NEXUSDataReader
 from ccpi.astra.operators import AstraProjector3DSimple
 
 # All external imports
 import numpy
 import matplotlib.pyplot as plt
-import os
-
-# Define utility function to average over flat and dark images.
-def avg_img(image):
-    shape = list(numpy.shape(image))
-    l = shape.pop(0)
-    avg = numpy.zeros(shape)
-    for i in range(l):
-        avg += image[i] / l
-    return avg
-    
-# Set up a reader object pointing to the Nexus data set. Revise path as needed.
-reader = NexusReader(os.path.join(".." ,".." ,".." , "..", "CCPi-ReconstructionFramework","data" , "24737_fd.nxs" ))
-
-# Read and print the dimensions of the raw projections
-dims = reader.get_projection_dimensions()
-print (dims)
-
-# Load and average all flat and dark images in preparation for normalising data.
-flat = avg_img(reader.load_flat())
-dark = avg_img(reader.load_dark())
-
-# Set up normaliser object for normalising data by flat and dark images.
-norm = Normalizer(flat_field=flat, dark_field=dark)
-
-# Load the raw projections and pass as input to the normaliser.
-norm.set_input(reader.get_acquisition_data())
 
-# Set up CenterOfRotationFinder object to center data.
-cor = CenterOfRotationFinder()
+## Set up a reader object pointing to the Nexus data set. Revise path as needed.
+# The data is already  corrected for by flat and dark field.
+myreader = NEXUSDataReader(nexus_file="/media/newhd/shared/Data/nexus/24737_fd_normalised.nxs" )
+data= myreader.load_data()
 
+# Negative logarithm transoform.
+data.fill( -numpy.log(data.as_array() ))
+
+# Set up CenterOfRotationFinder object to center data.
 # Set the output of the normaliser as the input and execute to determine center.
-cor.set_input(norm.get_output())
+cor = CenterOfRotationFinder()
+cor.set_input(data)
 center_of_rotation = cor.get_output()
 
-# Set up AcquisitionDataPadder to pad data for centering using the computed 
-# center, set the output of the normaliser as input and execute to produce
-# padded/centered data.
-padder = AcquisitionDataPadder(center_of_rotation=center_of_rotation)
-padder.set_input(norm.get_output())
-padded_data = padder.get_output()
+# From computed center, determine amount of zero-padding to apply, apply
+# and update geometry to wider detector.
+cor_pad = int(2*(center_of_rotation - data.shape[2]/2))
+data_pad = numpy.zeros((data.shape[0],data.shape[1],data.shape[2]+cor_pad))
+data_pad[:,:,:-cor_pad] = data.as_array()
+data.geometry.pixel_num_h = data.geometry.pixel_num_h + cor_pad
+data.array = data_pad
 
 # Permute array and convert angles to radions for ASTRA
-padded_data2 = padded_data.subset(dimensions=['vertical','angle','horizontal'])
-padded_data2.geometry = padded_data.geometry
-padded_data2.geometry.angles = padded_data2.geometry.angles/180*numpy.pi
+padded_data = data.subset(dimensions=['vertical','angle','horizontal'])
+padded_data.geometry = data.geometry
+padded_data.geometry.angles = padded_data.geometry.angles/180*numpy.pi
 
 # Create Acquisition and Image Geometries for setting up projector.
-ag = padded_data2.geometry
+ag = padded_data.geometry
 ig = ImageGeometry(voxel_num_x=ag.pixel_num_h,
                    voxel_num_y=ag.pixel_num_h, 
                    voxel_num_z=ag.pixel_num_v)
@@ -78,94 +57,150 @@ print ("Define projector")
 Cop = AstraProjector3DSimple(ig, ag)
 
 # Test backprojection and projection
-z1 = Cop.adjoint(padded_data2)
+z1 = Cop.adjoint(padded_data)
 z2 = Cop.direct(z1)
 
-plt.imshow(z1.subset(vertical=68).array)
+plt.imshow(z1.subset(horizontal_x=68).array)
 plt.show()
 
-# Set initial guess
+# Set initial guess for reconstruction algorithms.
 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(Cop,padded_data2,c=0.5)
-        
-# Run FISTA reconstruction for least squares without regularization
-print ("Run FISTA for least squares")
+# Set tolerance and number of iterations for reconstruction algorithms.
 opt = {'tol': 1e-4, 'iter': 100}
-x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt)
 
-plt.imshow(x_fista0.subset(horizontal_x=80).array)
-plt.title('FISTA LS')
-plt.colorbar()
-plt.show()
+# First a CGLS reconstruction can be done:
+CGLS_alg = CGLS()
+CGLS_alg.set_up(x_init, Cop, padded_data)
+CGLS_alg.max_iteration = 2000
+CGLS_alg.run(opt['iter'])
 
-# Set up 1-norm function for FISTA least squares plus 1-norm regularisation
-print ("Run FISTA for least squares plus 1-norm regularisation")
-lam = 0.1
-g0 = Norm1(lam)
+x_CGLS = CGLS_alg.get_output()
 
-# Run FISTA for least squares plus 1-norm function.
-x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0,opt=opt)
+# Fix color and slices for display
+v1 = -0.01
+v2 = 0.13
+hx=80
+hy=80
+v=68
 
-plt.imshow(x_fista0.subset(horizontal_x=80).array)
-plt.title('FISTA LS+1')
-plt.colorbar()
-plt.show()
+# Display ortho slices of reconstruction
+# Display all reconstructions and decay of objective function
+cols = 3
+rows = 1
+fig = plt.figure()
 
-# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
-print ("Run FBPD for least squares plus 1-norm regularisation")
-x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
+current = 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('horizontal_x')
+imgplot = plt.imshow(x_CGLS.subset(horizontal_x=hx).as_array(),vmin=v1,vmax=v2)
 
-plt.imshow(x_fbpd1.subset(horizontal_x=80).array)
-plt.title('FBPD LS+1')
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('horizontal_y')
+imgplot = plt.imshow(x_CGLS.subset(horizontal_y=hy).as_array(),vmin=v1,vmax=v2)
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('vertical')
+imgplot = plt.imshow(x_CGLS.subset(vertical=v).as_array(),vmin=v1,vmax=v2)
 plt.colorbar()
+
+plt.suptitle('CGLS reconstruction slices')
 plt.show()
 
-# Run CGLS, which should agree with the FISTA least squares
-print ("Run CGLS for least squares")
-x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, padded_data2, opt=opt)
-plt.imshow(x_CGLS.subset(horizontal_x=80).array)
-plt.title('CGLS')
-plt.colorbar()
+plt.figure()
+plt.semilogy(CGLS_alg.objective)
+plt.title('CGLS criterion')
 plt.show()
 
+# Create least squares object instance with projector and data.
+print ("Create least squares object instance with projector and data.")
+f = Norm2Sq(Cop,padded_data,c=0.5)
+
+
+# Run FISTA for least squares without constraints
+FISTA_alg = FISTA()
+FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
+FISTA_alg.max_iteration = 2000
+FISTA_alg.run(opt['iter'])
+x_FISTA = FISTA_alg.get_output()
+
+# Display ortho slices of reconstruction
 # Display all reconstructions and decay of objective function
-cols = 4
-rows = 1
-current = 1
+
 fig = plt.figure()
 
-current = current 
+current = 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA LS')
-imgplot = plt.imshow(x_fista0.subset(horizontal_x=80).as_array())
+a.set_title('horizontal_x')
+imgplot = plt.imshow(x_FISTA.subset(horizontal_x=hx).as_array(),vmin=v1,vmax=v2)
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA LS+1')
-imgplot = plt.imshow(x_fista1.subset(horizontal_x=80).as_array())
+a.set_title('horizontal_y')
+imgplot = plt.imshow(x_FISTA.subset(horizontal_y=hy).as_array(),vmin=v1,vmax=v2)
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FBPD LS+1')
-imgplot = plt.imshow(x_fbpd1.subset(horizontal_x=80).as_array())
+a.set_title('vertical')
+imgplot = plt.imshow(x_FISTA.subset(vertical=v).as_array(),vmin=v1,vmax=v2)
+plt.colorbar()
+
+plt.suptitle('FISTA Least squares reconstruction slices')
+plt.show()
+
+plt.figure()
+plt.semilogy(FISTA_alg.objective)
+plt.title('FISTA Least squares criterion')
+plt.show()
+
+# Create 1-norm function object
+lam = 30.0
+g0 = lam * L1Norm()
+
+# Run FISTA for least squares plus 1-norm function.
+FISTA_alg1 = FISTA()
+FISTA_alg1.set_up(x_init=x_init, f=f, g=g0, opt=opt)
+FISTA_alg1.max_iteration = 2000
+FISTA_alg1.run(opt['iter'])
+x_FISTA1 = FISTA_alg1.get_output()
+
+# Display all reconstructions and decay of objective function
+fig = plt.figure()
+
+current = 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('horizontal_x')
+imgplot = plt.imshow(x_FISTA1.subset(horizontal_x=hx).as_array(),vmin=v1,vmax=v2)
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('CGLS')
-imgplot = plt.imshow(x_CGLS.subset(horizontal_x=80).as_array())
+a.set_title('horizontal_y')
+imgplot = plt.imshow(x_FISTA1.subset(horizontal_y=hy).as_array(),vmin=v1,vmax=v2)
 
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('vertical')
+imgplot = plt.imshow(x_FISTA1.subset(vertical=v).as_array(),vmin=v1,vmax=v2)
+plt.colorbar()
+
+plt.suptitle('FISTA LS+1 reconstruction slices')
+plt.show()
+
+
+plt.figure()
+plt.semilogy(FISTA_alg1.objective)
+plt.title('FISTA LS+1 criterion')
 plt.show()
 
+
 fig = plt.figure()
 b=fig.add_subplot(1,1,1)
 b.set_title('criteria')
-imgplot = plt.loglog(criter0 , label='FISTA LS')
-imgplot = plt.loglog(criter1 , label='FISTA LS+1')
-imgplot = plt.loglog(criter_fbpd1, label='FBPD LS+1')
-imgplot = plt.loglog(criter_CGLS, label='CGLS')
+imgplot = plt.loglog(CGLS_alg.objective , label='CGLS')
+imgplot = plt.loglog(FISTA_alg.objective , label='FISTA LS')
+imgplot = plt.loglog(FISTA_alg1.objective, label='FISTA LS+1')
 b.legend(loc='right')
 plt.show()
diff --git a/Wrappers/Python/wip/demo_astra_simple.py b/Wrappers/Python/wip/demo_astra_simple.py
index 57caaee..925df77 100755
--- a/Wrappers/Python/wip/demo_astra_simple.py
+++ b/Wrappers/Python/wip/demo_astra_simple.py
@@ -2,12 +2,12 @@
 # This demo illustrates how ASTRA 2D projectors can be used with
 # the modular optimisation framework. The demo sets up a 2D test case and 
 # demonstrates reconstruction using CGLS, as well as FISTA for least squares 
-# and 1-norm regularisation and FBPD for 1-norm and TV regularisation.
+# and 1-norm regularisation.
 
 # First make all imports
 from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry
-from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D
+from ccpi.optimisation.algorithms import FISTA, CGLS
+from ccpi.optimisation.functions import Norm2Sq, L1Norm
 from ccpi.astra.operators import AstraProjectorSimple
 
 import numpy as np
@@ -77,22 +77,30 @@ plt.show()
 
 plt.imshow(z.array)
 plt.title('Backprojected data')
+plt.colorbar()
 plt.show()
 
 # Using the test data b, different reconstruction methods can now be set up as
 # demonstrated in the rest of this file. In general all methods need an initial 
 # guess and some algorithm options to be set:
-x_init = ImageData(np.zeros(x.shape),geometry=ig)
-opt = {'tol': 1e-4, 'iter': 1000}
+x_init = ImageData(geometry=ig)
+opt = {'tol': 1e-4, 'iter': 100}
 
 # First a CGLS reconstruction can be done:
-x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt)
+CGLS_alg = CGLS()
+CGLS_alg.set_up(x_init, Aop, b )
+CGLS_alg.max_iteration = 2000
+CGLS_alg.run(opt['iter'])
 
+x_CGLS = CGLS_alg.get_output()
+
+plt.figure()
 plt.imshow(x_CGLS.array)
 plt.title('CGLS')
 plt.show()
 
-plt.semilogy(criter_CGLS)
+plt.figure()
+plt.semilogy(CGLS_alg.objective)
 plt.title('CGLS criterion')
 plt.show()
 
@@ -102,72 +110,56 @@ plt.show()
 
 # Create least squares object instance with projector, test data and a constant 
 # coefficient of 0.5:
-f = Norm2sq(Aop,b,c=0.5)
-
-# Run FISTA for least squares without regularization
-x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None,opt)
-
-plt.imshow(x_fista0.array)
-plt.title('FISTA Least squares')
+f = Norm2Sq(Aop,b,c=0.5)
+
+# Run FISTA for least squares without constraints
+FISTA_alg = FISTA()
+FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
+FISTA_alg.max_iteration = 2000
+FISTA_alg.run(opt['iter'])
+x_FISTA = FISTA_alg.get_output()
+
+plt.figure()
+plt.imshow(x_FISTA.array)
+plt.title('FISTA Least squares reconstruction')
+plt.colorbar()
 plt.show()
 
-plt.semilogy(criter0)
+plt.figure()
+plt.semilogy(FISTA_alg.objective)
 plt.title('FISTA Least squares criterion')
 plt.show()
 
+
 # FISTA can also solve regularised forms by specifying a second function object
 # such as 1-norm regularisation with choice of regularisation parameter lam:
 
 # Create 1-norm function object
-lam = 0.1
-g0 = Norm1(lam)
+lam = 1.0
+g0 = lam * L1Norm()
 
 # Run FISTA for least squares plus 1-norm function.
-x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
-
-plt.imshow(x_fista1.array)
-plt.title('FISTA Least squares plus 1-norm regularisation')
-plt.show()
-
-plt.semilogy(criter1)
-plt.title('FISTA Least squares plus 1-norm regularisation criterion')
-plt.show()
-
-# The least squares plus 1-norm regularisation problem can also be solved by 
-# other algorithms such as the Forward Backward Primal Dual algorithm. This
-# algorithm minimises the sum of three functions and the least squares and 
-# 1-norm functions should be given as the second and third function inputs. 
-# In this test case, this algorithm requires more iterations to converge, so
-# new options are specified.
-opt_FBPD = {'tol': 1e-4, 'iter': 20000}
-x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt_FBPD)
-
-plt.imshow(x_fbpd1.array)
-plt.title('FBPD for least squares plus 1-norm regularisation')
-plt.show()
-
-plt.semilogy(criter_fbpd1)
-plt.title('FBPD for least squares plus 1-norm regularisation criterion')
-plt.show()
-
-# The FBPD algorithm can also be used conveniently for TV regularisation:
-
-# Specify TV function object
-lamtv = 10
-gtv = TV2D(lamtv)
-
-x_fbpdtv,it_fbpdtv,timing_fbpdtv,criter_fbpdtv=FBPD(x_init,None,f,gtv,opt_FBPD)
-
-plt.imshow(x_fbpdtv.array)
+FISTA_alg1 = FISTA()
+FISTA_alg1.set_up(x_init=x_init, f=f, g=g0, opt=opt)
+FISTA_alg1.max_iteration = 2000
+FISTA_alg1.run(opt['iter'])
+x_FISTA1 = FISTA_alg1.get_output()
+
+plt.figure()
+plt.imshow(x_FISTA1.array)
+plt.title('FISTA LS+L1Norm reconstruction')
+plt.colorbar()
 plt.show()
 
-plt.semilogy(criter_fbpdtv)
+plt.figure()
+plt.semilogy(FISTA_alg1.objective)
+plt.title('FISTA LS+L1norm criterion')
 plt.show()
 
 
 # Compare all reconstruction and criteria
 clims = (0,1)
-cols = 3
+cols = 2
 rows = 2
 current = 1
 
@@ -186,34 +178,20 @@ plt.axis('off')
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
 a.set_title('FISTA LS')
-imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1])
+imgplot = plt.imshow(x_FISTA.as_array(),vmin=clims[0],vmax=clims[1])
 plt.axis('off')
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
 a.set_title('FISTA LS+1')
-imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1])
-plt.axis('off')
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('FBPD LS+1')
-imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1])
-plt.axis('off')
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('FBPD TV')
-imgplot = plt.imshow(x_fbpdtv.as_array(),vmin=clims[0],vmax=clims[1])
+imgplot = plt.imshow(x_FISTA1.as_array(),vmin=clims[0],vmax=clims[1])
 plt.axis('off')
 
 fig = plt.figure()
-b=fig.add_subplot(1,1,1)
-b.set_title('criteria')
-imgplot = plt.loglog(criter_CGLS, label='CGLS')
-imgplot = plt.loglog(criter0 , label='FISTA LS')
-imgplot = plt.loglog(criter1 , label='FISTA LS+1')
-imgplot = plt.loglog(criter_fbpd1, label='FBPD LS+1')
-imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV')
-b.legend(loc='lower left')
+a=fig.add_subplot(1,1,1)
+a.set_title('criteria')
+imgplot = plt.loglog(CGLS_alg.objective, label='CGLS')
+imgplot = plt.loglog(FISTA_alg.objective , label='FISTA LS')
+imgplot = plt.loglog(FISTA_alg1.objective , label='FISTA LS+1')
+a.legend(loc='lower left')
 plt.show()
diff --git a/Wrappers/Python/wip/demo_astra_sophiabeads.py b/Wrappers/Python/wip/demo_astra_sophiabeads.py
index b72a5f9..dce3d91 100755
--- a/Wrappers/Python/wip/demo_astra_sophiabeads.py
+++ b/Wrappers/Python/wip/demo_astra_sophiabeads.py
@@ -1,25 +1,27 @@
 
 # This demo shows how to load a Nikon XTek micro-CT data set and reconstruct
-# the central slice using the CGLS method. The SophiaBeads dataset with 256 
-# projections is used as test data and can be obtained from here:
+# the central slice using the CGLS and FISTA methods. The SophiaBeads dataset 
+# with 64 projections is used as test data and can be obtained from here:
 # https://zenodo.org/record/16474
 # The filename with full path to the .xtekct file should be given as string 
-# input to XTEKReader to  load in the data.
+# input to NikonDataReader to  load in the data.
 
 # Do all imports
-from ccpi.io.reader import XTEKReader
 import numpy as np
 import matplotlib.pyplot as plt
+from ccpi.io import NikonDataReader
 from ccpi.framework import ImageGeometry, AcquisitionGeometry, AcquisitionData, ImageData
 from ccpi.astra.operators import AstraProjectorSimple
-from ccpi.optimisation.algs import CGLS
+from ccpi.optimisation.algorithms import FISTA, CGLS
+from ccpi.optimisation.functions import Norm2Sq, L1Norm
 
-# Set up reader object and read the data
-datareader = XTEKReader("REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_256_averaged.xtekct")
-data = datareader.get_acquisition_data()
+# Set up reader object and read in central slice the data
+datareader = NikonDataReader(xtek_file="REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_64_averaged.xtekct",
+                             roi=[(1000,1001),(0,2000)])
+data = datareader.load_projections()
 
 # Extract central slice, scale and negative-log transform
-sino = -np.log(data.as_array()[:,:,1000]/60000.0)
+sino = -np.log(data.as_array()[:,0,:]/60000.0)
 
 # Apply centering correction by zero padding, amount found manually
 cor_pad = 30
@@ -59,13 +61,74 @@ ig2d = ImageGeometry(voxel_num_x=N,
 Aop = AstraProjectorSimple(ig2d, ag2d,"gpu")
 
 # Set initial guess for CGLS reconstruction
-x_init = ImageData(np.zeros((N,N)),geometry=ig2d)
+x_init = ImageData(geometry=ig2d)
 
-# Run 50-iteration CGLS reconstruction
+# Set tolerance and number of iterations for reconstruction algorithms.
 opt = {'tol': 1e-4, 'iter': 50}
-x, it, timing, criter = CGLS(x_init,Aop,data2d,opt=opt)
 
-# Display reconstruction
+# First a CGLS reconstruction can be done:
+CGLS_alg = CGLS()
+CGLS_alg.set_up(x_init, Aop, data2d)
+CGLS_alg.max_iteration = 2000
+CGLS_alg.run(opt['iter'])
+
+x_CGLS = CGLS_alg.get_output()
+
+plt.figure()
+plt.imshow(x_CGLS.array)
+plt.title('CGLS')
+plt.colorbar()
+plt.show()
+
+plt.figure()
+plt.semilogy(CGLS_alg.objective)
+plt.title('CGLS criterion')
+plt.show()
+
+# CGLS solves the simple least-squares problem. The same problem can be solved 
+# by FISTA by setting up explicitly a least squares function object and using 
+# no regularisation:
+
+# Create least squares object instance with projector, test data and a constant 
+# coefficient of 0.5:
+f = Norm2Sq(Aop,data2d)
+
+# Run FISTA for least squares without constraints
+FISTA_alg = FISTA()
+FISTA_alg.set_up(x_init=x_init, f=f, opt=opt)
+FISTA_alg.max_iteration = 2000
+FISTA_alg.run(opt['iter'])
+x_FISTA = FISTA_alg.get_output()
+
+plt.figure()
+plt.imshow(x_FISTA.array)
+plt.title('FISTA Least squares reconstruction')
+plt.colorbar()
+plt.show()
+
+plt.figure()
+plt.semilogy(FISTA_alg.objective)
+plt.title('FISTA Least squares criterion')
+plt.show()
+
+# Create 1-norm function object
+lam = 1.0
+g0 = lam * L1Norm()
+
+# Run FISTA for least squares plus 1-norm function.
+FISTA_alg1 = FISTA()
+FISTA_alg1.set_up(x_init=x_init, f=f, g=g0, opt=opt)
+FISTA_alg1.max_iteration = 2000
+FISTA_alg1.run(opt['iter'])
+x_FISTA1 = FISTA_alg1.get_output()
+
+plt.figure()
+plt.imshow(x_FISTA1.array)
+plt.title('FISTA LS+L1Norm reconstruction')
+plt.colorbar()
+plt.show()
+
 plt.figure()
-plt.imshow(x.as_array())
-plt.colorbar()
\ No newline at end of file
+plt.semilogy(FISTA_alg1.objective)
+plt.title('FISTA LS+L1norm criterion')
+plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/demo_astra_sophiabeads3D.py b/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
index 8c43657..624f121 100644
--- a/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
+++ b/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
@@ -1,38 +1,31 @@
 
 # This demo shows how to load a Nikon XTek micro-CT data set and reconstruct
-# the central slice using the CGLS method. The SophiaBeads dataset with 256 
+# a volume of 200 slices using the CGLS method. The SophiaBeads dataset with 64 
 # projections is used as test data and can be obtained from here:
 # https://zenodo.org/record/16474
 # The filename with full path to the .xtekct file should be given as string 
-# input to XTEKReader to  load in the data.
+# input to NikonDataReader to  load in the data.
 
 # Do all imports
-from ccpi.io.reader import XTEKReader
 import numpy as np
 import matplotlib.pyplot as plt
-from ccpi.framework import ImageGeometry, AcquisitionData, ImageData
+from ccpi.io import NikonDataReader
+from ccpi.framework import ImageGeometry, ImageData
 from ccpi.astra.operators import AstraProjector3DSimple
-from ccpi.optimisation.algs import CGLS
+from ccpi.optimisation.algorithms import CGLS
 
-import numpy
-
-# Set up reader object and read the data
-datareader = XTEKReader("REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_256_averaged.xtekct")
-data = datareader.get_acquisition_data()
-
-# Crop data and fix dimension labels
-data = AcquisitionData(data.array[:,:,901:1101],
-                            geometry=data.geometry,
-                            dimension_labels=['angle','horizontal','vertical'])
-data.geometry.pixel_num_v = 200
+# Set up reader object and read in 200 central slices of the data
+datareader= NikonDataReader(xtek_file="REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_64_averaged.xtekct",
+                            roi=[(901,1101),(0,2000)])
+data = datareader.load_projections()
 
 # Scale and negative-log transform
-data.array = -np.log(data.as_array()/60000.0)
+data.fill(-np.log(data.as_array()/60000.0))
 
 # Apply centering correction by zero padding, amount found manually
 cor_pad = 30
-data_pad = np.zeros((data.shape[0],data.shape[1]+cor_pad,data.shape[2]))
-data_pad[:,cor_pad:,:] = data.array
+data_pad = np.zeros((data.shape[0],data.shape[1],data.shape[2]+cor_pad))
+data_pad[:,:,cor_pad:] = data.array
 data.geometry.pixel_num_h = data.geometry.pixel_num_h + cor_pad
 data.array = data_pad
 
@@ -58,7 +51,7 @@ ig = ImageGeometry(voxel_num_x=N,
 # Permute array and convert angles to radions for ASTRA; delete old data.
 data2 = data.subset(dimensions=['vertical','angle','horizontal'])
 data2.geometry = data.geometry
-data2.geometry.angles = -data2.geometry.angles/180*numpy.pi
+data2.geometry.angles = -data2.geometry.angles/180*np.pi
 del data
 
 # Extract the Acquisition geometry for setting up projector.
@@ -67,7 +60,7 @@ ag = data2.geometry
 # Set up the Projector (AcquisitionModel) using ASTRA on GPU
 Aop = AstraProjector3DSimple(ig, ag)
 
-# So and show simple backprojection
+# Do and show simple backprojection
 z = Aop.adjoint(data2)
 plt.figure()
 plt.imshow(z.subset(horizontal_x=1000).as_array())
@@ -82,17 +75,29 @@ plt.show()
 # Set initial guess for CGLS reconstruction
 x_init = ImageData(geometry=ig)
 
-# Run 50-iteration CGLS reconstruction
-opt = {'tol': 1e-4, 'iter': 20}
-x, it, timing, criter = CGLS(x_init,Aop,data2,opt=opt)
+# Set tolerance and number of iterations for reconstruction algorithms.
+opt = {'tol': 1e-4, 'iter': 30}
+
+# Run a CGLS reconstruction can be done:
+CGLS_alg = CGLS()
+CGLS_alg.set_up(x_init, Aop, data2)
+CGLS_alg.max_iteration = 2000
+CGLS_alg.run(opt['iter'])
+
+x_CGLS = CGLS_alg.get_output()
 
 # Display ortho slices of reconstruction
 plt.figure()
-plt.imshow(x.subset(horizontal_x=1000).as_array())
+plt.imshow(x_CGLS.subset(horizontal_x=1000).as_array())
 plt.show()
 plt.figure()
-plt.imshow(x.subset(horizontal_y=1000).as_array())
+plt.imshow(x_CGLS.subset(horizontal_y=1000).as_array())
 plt.show()
 plt.figure()
-plt.imshow(x.subset(vertical=100).as_array())
+plt.imshow(x_CGLS.subset(vertical=100).as_array())
 plt.show()
+
+plt.figure()
+plt.semilogy(CGLS_alg.objective)
+plt.title('CGLS criterion')
+plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/work_out_adjoint.py b/Wrappers/Python/wip/work_out_adjoint.py
deleted file mode 100644
index 276fb76..0000000
--- a/Wrappers/Python/wip/work_out_adjoint.py
+++ /dev/null
@@ -1,115 +0,0 @@
-
-# This demo illustrates how ASTRA 2D projectors can be used with
-# the modular optimisation framework. The demo sets up a 2D test case and 
-# demonstrates reconstruction using CGLS, as well as FISTA for least squares 
-# and 1-norm regularisation and FBPD for 1-norm and TV regularisation.
-
-# First make all imports
-from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry, AcquisitionData
-from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D
-from ccpi.astra.operators import AstraProjectorSimple
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
-test_case = 2
-
-# 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.
-N = 300
-x1 = -1
-x2 = 1
-dx = (x2-x1)/N
-ig = ImageGeometry(voxel_num_x=N,
-                   voxel_num_y=N,
-                   voxel_size_x=dx,
-                   voxel_size_y=dx)
-Phantom = ImageData(geometry=ig)
-
-x = Phantom.as_array()
-x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-x[:,round(3*N/8):round(5*N/8)] = 1
-
-plt.imshow(x)
-plt.title('Phantom image')
-plt.colorbar()
-plt.show()
-
-# Set up AcquisitionGeometry object to hold the parameters of the measurement
-# setup geometry: # Number of angles, the actual angles from 0 to 
-# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector 
-# pixel relative to an object pixel, the number of detector pixels, and the 
-# source-origin and origin-detector distance (here the origin-detector distance 
-# set to 0 to simulate a "virtual detector" with same detector pixel size as
-# object pixel size).
-angles_num = 40
-det_num = 400
-
-SourceOrig = 20
-OrigDetec = 100
-
-geo_mag = (SourceOrig+OrigDetec)/SourceOrig
-
-det_w = geo_mag*dx*1
-
-if test_case==1:
-    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
-    ag = AcquisitionGeometry('parallel',
-                             '2D',
-                             angles,
-                             det_num,det_w)
-elif test_case==2:
-    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
-    ag = AcquisitionGeometry('cone',
-                             '2D',
-                             angles,
-                             det_num,
-                             det_w,
-                             dist_source_center=SourceOrig, 
-                             dist_center_detector=OrigDetec)
-else:
-    NotImplemented
-
-# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
-# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
-Aop = AstraProjectorSimple(ig, ag, 'gpu')
-
-# Forward and backprojection are available as methods direct and adjoint. Here 
-# generate test data b and do simple backprojection to obtain z.
-b = Aop.direct(Phantom)
-z = Aop.adjoint(b)
-
-plt.imshow(b.array)
-plt.title('Simulated data')
-plt.colorbar()
-plt.show()
-
-plt.imshow(z.array)
-plt.title('Backprojected data')
-plt.colorbar()
-plt.show()
-
-
-One = ImageData(geometry=ig)
-xOne = One.as_array()
-xOne[:,:] = 1.0
-
-OneD = AcquisitionData(geometry=ag)
-y1 = OneD.as_array()
-y1[:,:] = 1.0
-
-s1 = (OneD*(Aop.direct(One))).sum()
-s2 = (One*(Aop.adjoint(OneD))).sum()
-print(s1)
-print(s2)
-print(s2/s1)
-
-print((b*b).sum())
-print((z*Phantom).sum())
-print((z*Phantom).sum() / (b*b).sum())
-
-#print(N/det_num)
-#print(0.5*det_w/dx)
\ No newline at end of file
diff --git a/Wrappers/Python/wip/work_out_adjoint3D.py b/Wrappers/Python/wip/work_out_adjoint3D.py
deleted file mode 100644
index 0e4f847..0000000
--- a/Wrappers/Python/wip/work_out_adjoint3D.py
+++ /dev/null
@@ -1,131 +0,0 @@
-
-# This demo illustrates how ASTRA 2D projectors can be used with
-# the modular optimisation framework. The demo sets up a 2D test case and 
-# demonstrates reconstruction using CGLS, as well as FISTA for least squares 
-# and 1-norm regularisation and FBPD for 1-norm and TV regularisation.
-
-# First make all imports
-from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry, AcquisitionData
-from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D
-from ccpi.astra.operators import AstraProjector3DSimple
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
-test_case = 2
-
-# 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.
-N = 100
-x1 = -1
-x2 = 1
-dx = (x2-x1)/N
-ig = ImageGeometry(voxel_num_x=N,
-                   voxel_num_y=N,
-                   voxel_num_z=N,
-                   voxel_size_x=dx,
-                   voxel_size_y=dx,
-                   voxel_size_z=dx)
-Phantom = ImageData(geometry=ig)
-
-x = Phantom.as_array()
-x[:,round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-x[:,round(3*N/8):round(5*N/8),:] = 1
-
-plt.imshow(x[90,:,:])
-plt.title('Phantom image')
-plt.colorbar()
-plt.show()
-
-# Set up AcquisitionGeometry object to hold the parameters of the measurement
-# setup geometry: # Number of angles, the actual angles from 0 to 
-# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector 
-# pixel relative to an object pixel, the number of detector pixels, and the 
-# source-origin and origin-detector distance (here the origin-detector distance 
-# set to 0 to simulate a "virtual detector" with same detector pixel size as
-# object pixel size).
-angles_num = 200
-det_num = 100
-
-det_w = dx
-
-if test_case==1:
-    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
-    ag = AcquisitionGeometry('parallel',
-                             '3D',
-                             angles,
-                             det_num,
-                             det_w,
-                             det_num,
-                             det_w)
-elif test_case==2:
-    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
-    
-    SourceOrig = 20
-    OrigDetec = 100
-
-    geo_mag = (SourceOrig+OrigDetec)/SourceOrig
-
-    det_w = geo_mag*dx
-    
-    ag = AcquisitionGeometry('cone',
-                             '3D',
-                             angles,
-                             det_num,
-                             det_w,
-                             det_num,
-                             det_w,
-                             dist_source_center=SourceOrig, 
-                             dist_center_detector=OrigDetec)
-else:
-    NotImplemented
-
-# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
-# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
-Aop = AstraProjector3DSimple(ig, ag)
-
-# Forward and backprojection are available as methods direct and adjoint. Here 
-# generate test data b and do simple backprojection to obtain z.
-b = Aop.direct(Phantom)
-z = Aop.adjoint(b)
-
-print((b*b).sum())
-print((z*Phantom).sum())
-print((z*Phantom).sum() / (b*b).sum())
-
-plt.imshow(b.array[:,0,:])
-plt.title('Simulated data')
-plt.colorbar()
-plt.show()
-plt.imshow(b.array[:,1,:])
-plt.title('Simulated data')
-plt.colorbar()
-plt.show()
-
-plt.imshow(z.array[50,:,:])
-plt.title('Backprojected data')
-plt.colorbar()
-plt.show()
-
-
-One = ImageData(geometry=ig)
-xOne = One.as_array()
-xOne[:,:,:] = 1.0
-
-OneD = b.clone()
-y1 = OneD.as_array()
-y1[:,:,:] = 1.0
-
-s1 = (OneD*(Aop.direct(One))).sum()
-s2 = (One*(Aop.adjoint(OneD))).sum()
-print(s1)
-print(s2)
-print(s2/s1)
-
-
-
-#print(N/det_num)
-#print(0.5*det_w/dx)
\ No newline at end of file
diff --git a/Wrappers/Python/wip/work_out_adjoint_sophiabeads.py b/Wrappers/Python/wip/work_out_adjoint_sophiabeads.py
deleted file mode 100644
index 7f35a16..0000000
--- a/Wrappers/Python/wip/work_out_adjoint_sophiabeads.py
+++ /dev/null
@@ -1,115 +0,0 @@
-
-# This demo illustrates how ASTRA 2D projectors can be used with
-# the modular optimisation framework. The demo sets up a 2D test case and 
-# demonstrates reconstruction using CGLS, as well as FISTA for least squares 
-# and 1-norm regularisation and FBPD for 1-norm and TV regularisation.
-
-# First make all imports
-from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry, AcquisitionData
-from ccpi.optimisation.algs import FISTA, FBPD, CGLS
-from ccpi.optimisation.funcs import Norm2sq, Norm1, TV2D
-from ccpi.astra.operators import AstraProjectorSimple
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
-test_case = 2
-
-# 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.
-N = 2000
-x1 = -16.015690364093633
-x2 =  16.015690364093633
-dx = (x2-x1)/N
-ig = ImageGeometry(voxel_num_x=N,
-                   voxel_num_y=N,
-                   voxel_size_x=dx,
-                   voxel_size_y=dx)
-Phantom = ImageData(geometry=ig)
-
-x = Phantom.as_array()
-x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-x[:,round(3*N/8):round(5*N/8)] = 1
-
-plt.imshow(x)
-plt.title('Phantom image')
-plt.colorbar()
-plt.show()
-
-# Set up AcquisitionGeometry object to hold the parameters of the measurement
-# setup geometry: # Number of angles, the actual angles from 0 to 
-# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector 
-# pixel relative to an object pixel, the number of detector pixels, and the 
-# source-origin and origin-detector distance (here the origin-detector distance 
-# set to 0 to simulate a "virtual detector" with same detector pixel size as
-# object pixel size).
-angles_num = 4
-det_num = 2500
-
-SourceOrig = 80.6392412185669
-OrigDetec = 926.3637587814331
-
-geo_mag = (SourceOrig+OrigDetec)/SourceOrig
-
-det_w = geo_mag*dx*1
-
-if test_case==1:
-    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
-    ag = AcquisitionGeometry('parallel',
-                             '2D',
-                             angles,
-                             det_num,det_w)
-elif test_case==2:
-    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
-    ag = AcquisitionGeometry('cone',
-                             '2D',
-                             angles,
-                             det_num,
-                             det_w,
-                             dist_source_center=SourceOrig, 
-                             dist_center_detector=OrigDetec)
-else:
-    NotImplemented
-
-# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
-# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
-Aop = AstraProjectorSimple(ig, ag, 'gpu')
-
-# Forward and backprojection are available as methods direct and adjoint. Here 
-# generate test data b and do simple backprojection to obtain z.
-b = Aop.direct(Phantom)
-z = Aop.adjoint(b)
-
-plt.imshow(b.array)
-plt.title('Simulated data')
-plt.colorbar()
-plt.show()
-
-plt.imshow(z.array)
-plt.title('Backprojected data')
-plt.colorbar()
-plt.show()
-
-
-One = ImageData(geometry=ig)
-xOne = One.as_array()
-xOne[:,:] = 1.0
-
-OneD = AcquisitionData(geometry=ag)
-y1 = OneD.as_array()
-y1[:,:] = 1.0
-
-s1 = (OneD*(Aop.direct(One))).sum()
-s2 = (One*(Aop.adjoint(OneD))).sum()
-print(s1)
-print(s2)
-print(s2/s1)
-
-print((b*b).sum())
-print((z*Phantom).sum())
-print((z*Phantom).sum() / (b*b).sum())
-
-#print(N/det_num)
-#print(0.5*det_w/dx)
\ No newline at end of file