From a344ebca5beddf376229415d63847b3edd545966 Mon Sep 17 00:00:00 2001
From: jakobsj <jakobsj@users.noreply.github.com>
Date: Fri, 14 Jun 2019 12:47:00 +0100
Subject: 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.
---
Wrappers/Python/wip/demo_astra_mc.py | 41 ++--
Wrappers/Python/wip/demo_astra_nexus.py | 233 ++++++++++++---------
Wrappers/Python/wip/demo_astra_simple.py | 132 +++++-------
Wrappers/Python/wip/demo_astra_sophiabeads.py | 93 ++++++--
Wrappers/Python/wip/demo_astra_sophiabeads3D.py | 59 +++---
Wrappers/Python/wip/work_out_adjoint.py | 115 ----------
Wrappers/Python/wip/work_out_adjoint3D.py | 131 ------------
.../Python/wip/work_out_adjoint_sophiabeads.py | 115 ----------
8 files changed, 325 insertions(+), 594 deletions(-)
delete mode 100644 Wrappers/Python/wip/work_out_adjoint.py
delete mode 100644 Wrappers/Python/wip/work_out_adjoint3D.py
delete mode 100644 Wrappers/Python/wip/work_out_adjoint_sophiabeads.py
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
--
cgit v1.2.3