3 files changed, 236 insertions, 284 deletions
diff --git a/Wrappers/Python/wip/demo_nexus.py b/Wrappers/Python/wip/demo_nexus.py
index 4dcc9f8..03739b1 100755
--- a/Wrappers/Python/wip/demo_nexus.py
+++ b/Wrappers/Python/wip/demo_nexus.py
@@ -1,27 +1,26 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Mar 21 14:26:21 2018
 
-@author: ofn77899
-"""
+# This script demonstrates how to load a parallel beam data set in Nexus 
+# format, apply dark and flat field correction and reconstruct using the
+# modular optimisation framework.
+# 
+# The data set is available from
+# https://github.com/DiamondLightSource/Savu/blob/master/test_data/data/24737_fd.nxs
+# and should be downloaded to a local directory to be specified below.
 
-from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry
+# 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.reconstruction.ccpiops import CCPiProjectorSimple
-from ccpi.reconstruction.parallelbeam import alg as pbalg
-from ccpi.reconstruction.processors import CCPiForwardProjector, CCPiBackwardProjector , \
-Normalizer , CenterOfRotationFinder , AcquisitionDataPadder
-
+from ccpi.plugins.ops import CCPiProjectorSimple
+from ccpi.processors import Normalizer, CenterOfRotationFinder, AcquisitionDataPadder
 from ccpi.io.reader import NexusReader
 
+# All external imports
 import numpy
 import matplotlib.pyplot as plt
-
 import os
-import pickle
-
 
+# Define utility function to average over flat and dark images.
 def avg_img(image):
     shape = list(numpy.shape(image))
     l = shape.pop(0)
@@ -30,116 +29,110 @@ def avg_img(image):
         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" ))
 
-reader = NexusReader(os.path.join(".." ,".." ,".." , "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 the output of the normaliser as the input and execute to determine center.
 cor.set_input(norm.get_output())
 center_of_rotation = cor.get_output()
-voxel_per_pixel = 1
 
+# 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()
 
-pg = padded_data.geometry
-geoms = pbalg.pb_setup_geometry_from_acquisition(padded_data.as_array(),
-                                                pg.angles,
-                                                center_of_rotation,
-                                                voxel_per_pixel )
-vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'],
-                   voxel_num_y=geoms['output_volume_y'], 
-                   voxel_num_z=geoms['output_volume_z'])
-#data = numpy.reshape(reader.getAcquisitionData())
-print ("define projector")
-Cop = CCPiProjectorSimple(vg, pg)
+# Create Acquisition and Image Geometries for setting up projector.
+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)
+
+# Define the projector object
+print ("Define projector")
+Cop = CCPiProjectorSimple(ig, ag)
+
 # 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)
+
+# Set initial guess
 print ("Initial guess")
-# Initial guess
-x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y','vertical'])
+x_init = ImageData(geometry=ig, dimension_labels=['horizontal_x','horizontal_y','vertical'])
         
-#%%
-print ("run FISTA")
-# Run FISTA for least squares without regularization
+# Run FISTA reconstruction for least squares without regularization
+print ("Run FISTA for least squares")
 opt = {'tol': 1e-4, 'iter': 10}
 x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt)
-pickle.dump(x_fista0, open("fista0.pkl", "wb"))
-
 
 plt.imshow(x_fista0.subset(horizontal_x=80).array)
-plt.title('FISTA0')
-#plt.show()
+plt.title('FISTA LS')
+plt.show()
 
-# Now least squares plus 1-norm regularization
+# 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)
 
 # Run FISTA for least squares plus 1-norm function.
 x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0,opt=opt)
-pickle.dump(x_fista1, open("fista1.pkl", "wb"))
 
 plt.imshow(x_fista0.subset(horizontal_x=80).array)
-plt.title('FISTA1')
-#plt.show()
-
-plt.semilogy(criter1)
-#plt.show()
+plt.title('FISTA LS+1')
+plt.show()
 
 # 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)
-pickle.dump(x_fbpd1, open("fbpd1.pkl", "wb"))
 
 plt.imshow(x_fbpd1.subset(horizontal_x=80).array)
-plt.title('FBPD1')
-#plt.show()
-
-plt.semilogy(criter_fbpd1)
-#plt.show()
+plt.title('FBPD LS+1')
+plt.show()
 
-# Run CGLS, which should agree with the FISTA0
+# 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_data, opt=opt)
-pickle.dump(x_CGLS, open("cgls.pkl", "wb"))
 plt.imshow(x_CGLS.subset(horizontal_x=80).array)
 plt.title('CGLS')
-plt.title('CGLS recon, compare FISTA0')
-#plt.show()
-
-plt.semilogy(criter_CGLS)
-plt.title('CGLS criterion')
-#plt.show()
-
+plt.show()
 
+# Display all reconstructions and decay of objective function
 cols = 4
 rows = 1
 current = 1
 fig = plt.figure()
-# projections row
 
 current = current 
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA0')
+a.set_title('FISTA LS')
 imgplot = plt.imshow(x_fista0.subset(horizontal_x=80).as_array())
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA1')
+a.set_title('FISTA LS+1')
 imgplot = plt.imshow(x_fista1.subset(horizontal_x=80).as_array())
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FBPD1')
+a.set_title('FBPD LS+1')
 imgplot = plt.imshow(x_fbpd1.subset(horizontal_x=80).as_array())
 
 current = current + 1
@@ -149,16 +142,12 @@ imgplot = plt.imshow(x_CGLS.subset(horizontal_x=80).as_array())
 
 plt.show()
 
-
-#%%
 fig = plt.figure()
-# projections row
 b=fig.add_subplot(1,1,1)
 b.set_title('criteria')
-imgplot = plt.loglog(criter0 , label='FISTA0')
-imgplot = plt.loglog(criter1 , label='FISTA1')
-imgplot = plt.loglog(criter_fbpd1, label='FBPD1')
+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(criter_fbpdtv, label='FBPD TV')
 b.legend(loc='right')
 plt.show()
 \ No newline at end of file
diff --git a/Wrappers/Python/wip/demo_simple_RGLTK.md b/Wrappers/Python/wip/demo_simple_RGLTK.md
deleted file mode 100644
index 9f0a4c3..0000000
--- a/Wrappers/Python/wip/demo_simple_RGLTK.md
+++ /dev/null
@@ -1,214 +0,0 @@
-
-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.astra.ops import AstraProjectorSimple
-from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-test_case = 1   # 1=parallel2D, 2=cone2D
-
-# Set up phantom
-N = 128
-
-
-vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
-Phantom = ImageData(geometry=vg)
-
-x = Phantom.as_array()
-x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-
-#plt.imshow(x)
-#plt.show()
-
-# Set up measurement geometry
-angles_num = 20; # angles number
-
-if test_case==1:
-    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
-elif test_case==2:
-    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
-else:
-    NotImplemented
-
-det_w = 1.0
-det_num = N
-SourceOrig = 200
-OrigDetec = 0
-
-# Parallelbeam geometry test
-if test_case==1:
-    pg = AcquisitionGeometry('parallel',
-                          '2D',
-                          angles,
-                          det_num,det_w)
-elif test_case==2:
-    pg = AcquisitionGeometry('cone',
-                          '2D',
-                          angles,
-                          det_num,
-                          det_w,
-                          dist_source_center=SourceOrig, 
-                          dist_center_detector=OrigDetec)
-
-# ASTRA operator using volume and sinogram geometries
-Aop = AstraProjectorSimple(vg, pg, 'cpu')
-
-# Unused old astra projector without geometry
-# Aop_old = AstraProjector(det_w, det_num, SourceOrig, 
-#                      OrigDetec, angles, 
-#                      N,'fanbeam','gpu') 
-
-# Try forward and backprojection
-b = Aop.direct(Phantom)
-out2 = Aop.adjoint(b)
-
-#plt.imshow(b.array)
-#plt.show()
-
-#plt.imshow(out2.array)
-#plt.show()
-
-# Create least squares object instance with projector and data.
-f = Norm2sq(Aop,b,c=0.5)
-
-# Initial guess
-x_init = ImageData(np.zeros(x.shape),geometry=vg)
-#%%
-# FISTA with ROF-TV regularisation
-g_rof = _ROF_TV_(lambdaReg = 10.0,iterationsTV=50,tolerance=1e-5,time_marchstep=0.01,device='cpu')
-
-opt = {'tol': 1e-4, 'iter': 100}
-
-x_fista_rof, it1, timing1, criter_rof = FISTA(x_init, f, g_rof,opt)
-
-plt.figure()
-plt.subplot(121)
-plt.imshow(x_fista_rof.array,cmap="BuPu")
-plt.title('FISTA-ROF-TV')
-plt.subplot(122)
-plt.semilogy(criter_rof)
-plt.show()
-#%%
-# FISTA with FGP-TV regularisation
-g_fgp = _FGP_TV_(lambdaReg = 10.0,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,cmap="BuPu")
-plt.title('FISTA-FGP-TV')
-plt.subplot(122)
-plt.semilogy(criter_fgp)
-plt.show()
-#%%
-# 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('FISTA0')
-plt.show()
-#%%
-# Now least squares plus 1-norm regularization
-lam = 0.1
-g0 = Norm1(lam)
-
-# Run FISTA for least squares plus 1-norm function.
-x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0)
-
-plt.imshow(x_fista1.array)
-plt.title('FISTA1')
-plt.show()
-
-plt.semilogy(criter1)
-plt.show()
-#%%
-# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
-opt = {'tol': 1e-4, 'iter': 100}
-x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
-
-plt.imshow(x_fbpd1.array)
-plt.title('FBPD1')
-plt.show()
-
-plt.semilogy(criter_fbpd1)
-plt.show()
-#%%
-opt_FBPD = {'tol': 1e-4, 'iter': 10000}
-# Now FBPD for least squares plus TV
-lamtv = 10.0
-gtv = TV2D(lamtv)
-
-x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt_FBPD)
-
-plt.imshow(x_fbpdtv.array)
-plt.show()
-
-plt.semilogy(criter_fbpdtv)
-plt.show()
-
-
-# Run CGLS, which should agree with the FISTA0
-x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt )
-
-plt.imshow(x_CGLS.array)
-plt.title('CGLS')
-#plt.title('CGLS recon, compare FISTA0')
-plt.show()
-
-plt.semilogy(criter_CGLS)
-plt.title('CGLS criterion')
-plt.show()
-#%%
-
-clims = (0,1)
-cols = 3
-rows = 2
-current = 1
-fig = plt.figure()
-# projections row
-a=fig.add_subplot(rows,cols,current)
-a.set_title('phantom {0}'.format(np.shape(Phantom.as_array())))
-
-imgplot = plt.imshow(Phantom.as_array(),vmin=clims[0],vmax=clims[1])
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA0')
-imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1])
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('FISTA1')
-imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1])
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('FBPD1')
-imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1])
-
-current = current + 1
-a=fig.add_subplot(rows,cols,current)
-a.set_title('CGLS')
-imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1])
-
-#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])
-
-fig = plt.figure()
-# projections row
-b=fig.add_subplot(1,1,1)
-b.set_title('criteria')
-imgplot = plt.loglog(criter0 , label='FISTA0')
-imgplot = plt.loglog(criter1 , label='FISTA1')
-imgplot = plt.loglog(criter_fbpd1, label='FBPD1')
-imgplot = plt.loglog(criter_CGLS, label='CGLS')
-b.legend(loc='right')
-plt.show()
-#%%
diff --git a/Wrappers/Python/wip/demo_simple_RGLTK.py b/Wrappers/Python/wip/demo_simple_RGLTK.py
new file mode 100644
index 0000000..3831603
--- /dev/null
+++ b/Wrappers/Python/wip/demo_simple_RGLTK.py
@@ -0,0 +1,177 @@
+
+# This demo illustrates how the CCPi Regularisation Toolkit can be used 
+# as TV regularisation for use with the FISTA algorithm of the modular 
+# optimisation framework and compares with the FBPD TV implementation.
+
+# All own 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.astra.ops import AstraProjectorSimple
+from ccpi.plugins.regularisers import _ROF_TV_, _FGP_TV_
+
+# All external imports
+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 = 1
+
+# 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 = 128
+ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
+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(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
+
+plt.imshow(x)
+plt.title('Phantom image')
+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 = 20
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
+
+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, 'cpu')
+
+# 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.show()
+
+plt.imshow(z.array)
+plt.title('Backprojected data')
+plt.show()
+
+# Create least squares object instance with projector and data.
+f = Norm2sq(Aop,b,c=0.5)
+
+# Initial guess
+x_init = ImageData(np.zeros(x.shape),geometry=ig)
+
+# Set up FBPD algorithm for TV reconstruction and solve
+opt_FBPD = {'tol': 1e-4, 'iter': 10000}
+
+lamtv = 1.0
+gtv = TV2D(lamtv)
+
+x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,
+                                                         None,
+                                                         f,
+                                                         gtv,
+                                                         opt=opt_FBPD)
+
+plt.figure()
+plt.subplot(121)
+plt.imshow(x_fbpdtv.array)
+plt.title('FBPD TV')
+plt.subplot(122)
+plt.semilogy(criter_fbpdtv)
+plt.show()
+
+# Set up the ROF variant of TV from the CCPi Regularisation Toolkit and run
+# TV-reconstruction using FISTA
+g_rof = _ROF_TV_(lambdaReg = lamtv,
+                 iterationsTV=50,
+                 tolerance=1e-5,
+                 time_marchstep=0.01,
+                 device='cpu')
+
+opt = {'tol': 1e-4, 'iter': 100}
+
+x_fista_rof, it1, timing1, criter_rof = FISTA(x_init, f, g_rof,opt)
+
+plt.figure()
+plt.subplot(121)
+plt.imshow(x_fista_rof.array)
+plt.title('FISTA ROF TV')
+plt.subplot(122)
+plt.semilogy(criter_rof)
+plt.show()
+
+# 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()
+
+# Compare all reconstruction and criteria
+clims = (0,1)
+cols = 3
+rows = 1
+current = 1
+fig = plt.figure()
+
+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])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA ROF TV')
+imgplot = plt.imshow(x_fista_rof.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA FGP TV')
+imgplot = plt.imshow(x_fista_fgp.as_array(),vmin=clims[0],vmax=clims[1])
+
+fig = plt.figure()
+
+b=fig.add_subplot(1,1,1)
+b.set_title('criteria')
+imgplot = plt.loglog(criter_fbpdtv , label='FBPD TV')
+imgplot = plt.loglog(criter_rof , label='ROF TV')
+imgplot = plt.loglog(criter_fgp, label='FGP TV')
+b.legend(loc='right')
+plt.show()