summaryrefslogtreecommitdiffstats
path: root/Wrappers
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers')
-rwxr-xr-xWrappers/Python/ccpi/astra/ops.py56
-rw-r--r--Wrappers/Python/ccpi/astra/processors.py168
-rwxr-xr-xWrappers/Python/ccpi/astra/utils.py9
-rwxr-xr-xWrappers/Python/wip/demo_astra_mc.py4
-rw-r--r--Wrappers/Python/wip/demo_astra_nexus.py171
-rwxr-xr-xWrappers/Python/wip/demo_astra_simple.py4
-rwxr-xr-xWrappers/Python/wip/demo_astra_sophiabeads.py10
-rw-r--r--Wrappers/Python/wip/demo_astra_sophiabeads3D.py98
-rw-r--r--Wrappers/Python/wip/work_out_adjoint.py115
-rw-r--r--Wrappers/Python/wip/work_out_adjoint3D.py131
-rw-r--r--Wrappers/Python/wip/work_out_adjoint_sophiabeads.py115
11 files changed, 860 insertions, 21 deletions
diff --git a/Wrappers/Python/ccpi/astra/ops.py b/Wrappers/Python/ccpi/astra/ops.py
index cd0ef9e..2e9d6af 100755
--- a/Wrappers/Python/ccpi/astra/ops.py
+++ b/Wrappers/Python/ccpi/astra/ops.py
@@ -17,11 +17,11 @@
from ccpi.optimisation.ops import Operator
import numpy
-import astra
from ccpi.framework import AcquisitionData, ImageData, DataContainer
from ccpi.optimisation.ops import PowerMethodNonsquare
from ccpi.astra.processors import AstraForwardProjector, AstraBackProjector, \
- AstraForwardProjectorMC, AstraBackProjectorMC
+ AstraForwardProjectorMC, AstraBackProjectorMC, AstraForwardProjector3D, \
+ AstraBackProjector3D
class AstraProjectorSimple(Operator):
"""ASTRA projector modified to use DataSet and geometry."""
@@ -74,6 +74,58 @@ class AstraProjectorSimple(Operator):
return DataContainer(numpy.random.randn(inputsize[0],
inputsize[1]))
+class AstraProjector3DSimple(Operator):
+ """ASTRA projector modified to use DataSet and geometry."""
+ def __init__(self, geomv, geomp):
+ super(AstraProjector3DSimple, self).__init__()
+
+ # Store volume and sinogram geometries.
+ self.sinogram_geometry = geomp
+ self.volume_geometry = geomv
+
+ self.fp = AstraForwardProjector3D(volume_geometry=geomv,
+ sinogram_geometry=geomp,
+ output_axes_order=['vertical','angle','horizontal'])
+
+ self.bp = AstraBackProjector3D(volume_geometry=geomv,
+ sinogram_geometry=geomp,
+ output_axes_order=['vertical','horizontal_y','horizontal_x'])
+
+ # Initialise empty for singular value.
+ self.s1 = None
+
+ def direct(self, IM):
+ self.fp.set_input(IM)
+ out = self.fp.get_output()
+ return out
+
+ def adjoint(self, DATA):
+ self.bp.set_input(DATA)
+ out = self.bp.get_output()
+ return out
+
+ #def delete(self):
+ # astra.data2d.delete(self.proj_id)
+
+ def get_max_sing_val(self):
+ self.s1, sall, svec = PowerMethodNonsquare(self,10)
+ return self.s1
+
+ def size(self):
+ # Only implemented for 2D
+ return ( (self.sinogram_geometry.angles.size, \
+ self.sinogram_geometry.pixel_num_h, \
+ self.sinogram_geometry.pixel_num_v,), \
+ (self.volume_geometry.voxel_num_x, \
+ self.volume_geometry.voxel_num_y, \
+ self.volume_geometry.voxel_num_z) )
+
+ def create_image_data(self):
+ inputsize = self.size()[1]
+ return DataContainer(numpy.random.randn(inputsize[2],
+ inputsize[1],
+ inputsize[0]))
+
class AstraProjectorMC(Operator):
"""ASTRA Multichannel projector"""
diff --git a/Wrappers/Python/ccpi/astra/processors.py b/Wrappers/Python/ccpi/astra/processors.py
index db52a6b..855f890 100644
--- a/Wrappers/Python/ccpi/astra/processors.py
+++ b/Wrappers/Python/ccpi/astra/processors.py
@@ -74,8 +74,15 @@ class AstraForwardProjector(DataProcessor):
sinogram_id, DATA.array = astra.create_sino(IM.as_array(),
self.proj_id)
astra.data2d.delete(sinogram_id)
- #return AcquisitionData(array=DATA, geometry=self.sinogram_geometry)
- return DATA
+
+ if self.device == 'cpu':
+ return DATA
+ else:
+ if self.sinogram_geometry.geom_type == 'cone':
+ return DATA
+ else:
+ scaling = 1.0/self.volume_geometry.voxel_size_x
+ return scaling*DATA
class AstraBackProjector(DataProcessor):
'''AstraBackProjector
@@ -145,7 +152,12 @@ class AstraBackProjector(DataProcessor):
rec_id, IM.array = astra.create_backprojection(DATA.as_array(),
self.proj_id)
astra.data2d.delete(rec_id)
- return IM
+
+ if self.device == 'cpu':
+ return IM
+ else:
+ scaling = self.volume_geometry.voxel_size_x**3
+ return scaling*IM
class AstraForwardProjectorMC(AstraForwardProjector):
'''AstraForwardProjector Multi channel
@@ -173,7 +185,15 @@ class AstraForwardProjectorMC(AstraForwardProjector):
sinogram_id, DATA.as_array()[k] = astra.create_sino(IM.as_array()[k],
self.proj_id)
astra.data2d.delete(sinogram_id)
- return DATA
+
+ if self.device == 'cpu':
+ return DATA
+ else:
+ if self.sinogram_geometry.geom_type == 'cone':
+ return DATA
+ else:
+ scaling = (1.0/self.volume_geometry.voxel_size_x)
+ return scaling*DATA
class AstraBackProjectorMC(AstraBackProjector):
'''AstraBackProjector Multi channel
@@ -202,4 +222,142 @@ class AstraBackProjectorMC(AstraBackProjector):
DATA.as_array()[k],
self.proj_id)
astra.data2d.delete(rec_id)
- return IM
+
+ if self.device == 'cpu':
+ return IM
+ else:
+ scaling = self.volume_geometry.voxel_size_x**3
+ return scaling*IM
+
+class AstraForwardProjector3D(DataProcessor):
+ '''AstraForwardProjector3D
+
+ Forward project ImageData to AcquisitionData using ASTRA proj_geom and
+ vol_geom.
+
+ Input: ImageData
+ Parameter: proj_geom, vol_geom
+ Output: AcquisitionData
+ '''
+
+ def __init__(self,
+ volume_geometry=None,
+ sinogram_geometry=None,
+ proj_geom=None,
+ vol_geom=None,
+ output_axes_order=None):
+ kwargs = {
+ 'volume_geometry' : volume_geometry,
+ 'sinogram_geometry' : sinogram_geometry,
+ 'proj_geom' : proj_geom,
+ 'vol_geom' : vol_geom,
+ 'output_axes_order' : output_axes_order
+ }
+
+ #DataProcessor.__init__(self, **kwargs)
+ super(AstraForwardProjector3D, self).__init__(**kwargs)
+
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ # Set up ASTRA Volume and projection geometry, not to be stored in self
+ vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+ self.sinogram_geometry)
+
+ # Also store ASTRA geometries
+ self.vol_geom = vol_geom
+ self.proj_geom = proj_geom
+
+ def check_input(self, dataset):
+ if dataset.number_of_dimensions == 3:
+ return True
+ else:
+ raise ValueError("Expected input dimensions 3, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+ def set_ImageGeometry(self, volume_geometry):
+ self.volume_geometry = volume_geometry
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+ self.sinogram_geometry = sinogram_geometry
+
+ def set_vol_geom(self, vol_geom):
+ self.vol_geom = vol_geom
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+ self.sinogram_geometry = sinogram_geometry
+
+ def process(self):
+ IM = self.get_input()
+ DATA = AcquisitionData(geometry=self.sinogram_geometry,
+ dimension_labels=self.output_axes_order)
+ sinogram_id, DATA.array = astra.create_sino3d_gpu(IM.as_array(),
+ self.proj_geom,
+ self.vol_geom)
+ astra.data3d.delete(sinogram_id)
+ # 3D CUDA FP does not need scaling
+ return DATA
+
+class AstraBackProjector3D(DataProcessor):
+ '''AstraBackProjector3D
+
+ Back project AcquisitionData to ImageData using ASTRA proj_geom, vol_geom.
+
+ Input: AcquisitionData
+ Parameter: proj_geom, vol_geom
+ Output: ImageData
+ '''
+
+ def __init__(self,
+ volume_geometry=None,
+ sinogram_geometry=None,
+ proj_geom=None,
+ vol_geom=None,
+ output_axes_order=None):
+ kwargs = {
+ 'volume_geometry' : volume_geometry,
+ 'sinogram_geometry' : sinogram_geometry,
+ 'proj_geom' : proj_geom,
+ 'vol_geom' : vol_geom,
+ 'output_axes_order' : output_axes_order
+ }
+
+ #DataProcessor.__init__(self, **kwargs)
+ super(AstraBackProjector3D, self).__init__(**kwargs)
+
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ # Set up ASTRA Volume and projection geometry, not to be stored in self
+ vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+ self.sinogram_geometry)
+
+ # Also store ASTRA geometries
+ self.vol_geom = vol_geom
+ self.proj_geom = proj_geom
+
+ def check_input(self, dataset):
+ if dataset.number_of_dimensions == 3:
+ return True
+ else:
+ raise ValueError("Expected input dimensions is 3, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+ def set_ImageGeometry(self, volume_geometry):
+ self.volume_geometry = volume_geometry
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+ self.sinogram_geometry = sinogram_geometry
+
+ def process(self):
+ DATA = self.get_input()
+ IM = ImageData(geometry=self.volume_geometry,
+ dimension_labels=self.output_axes_order)
+ rec_id, IM.array = astra.create_backprojection3d_gpu(DATA.as_array(),
+ self.proj_geom,
+ self.vol_geom)
+ astra.data3d.delete(rec_id)
+
+ # Scaling of 3D ASTRA backprojector, works both parallel and cone.
+ scaling = 1/self.volume_geometry.voxel_size_x**2
+ return scaling*IM
diff --git a/Wrappers/Python/ccpi/astra/utils.py b/Wrappers/Python/ccpi/astra/utils.py
index 9f8fe46..2b19305 100755
--- a/Wrappers/Python/ccpi/astra/utils.py
+++ b/Wrappers/Python/ccpi/astra/utils.py
@@ -1,4 +1,5 @@
import astra
+import numpy as np
def convert_geometry_to_astra(volume_geometry, sinogram_geometry):
# Set up ASTRA Volume and projection geometry, not stored
@@ -20,8 +21,8 @@ def convert_geometry_to_astra(volume_geometry, sinogram_geometry):
sinogram_geometry.pixel_size_h,
sinogram_geometry.pixel_num_h,
sinogram_geometry.angles,
- sinogram_geometry.dist_source_center,
- sinogram_geometry.dist_center_detector)
+ np.abs(sinogram_geometry.dist_source_center),
+ np.abs(sinogram_geometry.dist_center_detector))
else:
NotImplemented
@@ -50,8 +51,8 @@ def convert_geometry_to_astra(volume_geometry, sinogram_geometry):
sinogram_geometry.pixel_num_v,
sinogram_geometry.pixel_num_h,
sinogram_geometry.angles,
- sinogram_geometry.dist_source_center,
- sinogram_geometry.dist_center_detector)
+ np.abs(sinogram_geometry.dist_source_center),
+ np.abs(sinogram_geometry.dist_center_detector))
else:
NotImplemented
diff --git a/Wrappers/Python/wip/demo_astra_mc.py b/Wrappers/Python/wip/demo_astra_mc.py
index 3b78eb3..f09dcb8 100755
--- a/Wrappers/Python/wip/demo_astra_mc.py
+++ b/Wrappers/Python/wip/demo_astra_mc.py
@@ -97,7 +97,7 @@ 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)
+x_init = ImageData(numpy.zeros(x.shape),geometry=ig)
opt = {'tol': 1e-4, 'iter': 200}
# Create least squares object instance with projector, test data and a constant
@@ -120,7 +120,7 @@ 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.
# Again the regulariser is over all channels:
-lam = 0.1
+lam = 10
g0 = Norm1(lam)
# Run FISTA for least squares plus 1-norm function.
diff --git a/Wrappers/Python/wip/demo_astra_nexus.py b/Wrappers/Python/wip/demo_astra_nexus.py
new file mode 100644
index 0000000..1db44c0
--- /dev/null
+++ b/Wrappers/Python/wip/demo_astra_nexus.py
@@ -0,0 +1,171 @@
+
+# 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 and the ASTRA Tomography toolbox.
+#
+# 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.
+
+# 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.astra.ops 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 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()
+
+# 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()
+
+# 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
+
+# Create Acquisition and Image Geometries for setting up projector.
+ag = padded_data2.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 = AstraProjector3DSimple(ig, ag)
+
+# Test backprojection and projection
+z1 = Cop.adjoint(padded_data2)
+z2 = Cop.direct(z1)
+
+plt.imshow(z1.subset(vertical=68).array)
+plt.show()
+
+# Set initial guess
+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")
+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()
+
+# 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)
+
+plt.imshow(x_fista0.subset(horizontal_x=80).array)
+plt.title('FISTA LS+1')
+plt.colorbar()
+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)
+
+plt.imshow(x_fbpd1.subset(horizontal_x=80).array)
+plt.title('FBPD LS+1')
+plt.colorbar()
+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.show()
+
+# Display all reconstructions and decay of objective function
+cols = 4
+rows = 1
+current = 1
+fig = plt.figure()
+
+current = current
+a=fig.add_subplot(rows,cols,current)
+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('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('FBPD LS+1')
+imgplot = plt.imshow(x_fbpd1.subset(horizontal_x=80).as_array())
+
+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())
+
+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')
+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 8529c98..c1dd877 100755
--- a/Wrappers/Python/wip/demo_astra_simple.py
+++ b/Wrappers/Python/wip/demo_astra_simple.py
@@ -139,7 +139,7 @@ plt.show()
# 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': 7000}
+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)
@@ -153,7 +153,7 @@ plt.show()
# The FBPD algorithm can also be used conveniently for TV regularisation:
# Specify TV function object
-lamtv = 1
+lamtv = 10
gtv = TV2D(lamtv)
x_fbpdtv,it_fbpdtv,timing_fbpdtv,criter_fbpdtv=FBPD(x_init,None,f,gtv,opt_FBPD)
diff --git a/Wrappers/Python/wip/demo_astra_sophiabeads.py b/Wrappers/Python/wip/demo_astra_sophiabeads.py
index d8100ea..bcc775e 100755
--- a/Wrappers/Python/wip/demo_astra_sophiabeads.py
+++ b/Wrappers/Python/wip/demo_astra_sophiabeads.py
@@ -26,16 +26,13 @@ cor_pad = 30
sino_pad = np.zeros((sino.shape[0],sino.shape[1]+cor_pad))
sino_pad[:,cor_pad:] = sino
-# Simple beam hardening correction as done in SophiaBeads coda
-sino_pad = sino_pad**2
-
# Extract AcquisitionGeometry for central slice for 2D fanbeam reconstruction
ag2d = AcquisitionGeometry('cone',
'2D',
angles=-np.pi/180*data.geometry.angles,
pixel_num_h=data.geometry.pixel_num_h + cor_pad,
pixel_size_h=data.geometry.pixel_size_h,
- dist_source_center=data.geometry.dist_source_center,
+ dist_source_center=-data.geometry.dist_source_center,
dist_center_detector=data.geometry.dist_center_detector)
# Set up AcquisitionData object for central slice 2D fanbeam
@@ -65,9 +62,10 @@ Aop = AstraProjectorSimple(ig2d, ag2d,"gpu")
x_init = ImageData(np.zeros((N,N)),geometry=ig2d)
# Run 50-iteration CGLS reconstruction
-opt = {'tol': 1e-4, 'iter': 100}
+opt = {'tol': 1e-4, 'iter': 50}
x, it, timing, criter = CGLS(x_init,Aop,data2d,opt=opt)
# Display reconstruction
plt.figure()
-plt.imshow(x.as_array()) \ No newline at end of file
+plt.imshow(x.as_array())
+plt.colorbar() \ No newline at end of file
diff --git a/Wrappers/Python/wip/demo_astra_sophiabeads3D.py b/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
new file mode 100644
index 0000000..edfe1b9
--- /dev/null
+++ b/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
@@ -0,0 +1,98 @@
+
+# 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:
+# 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.
+
+# 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.astra.ops import AstraProjector3DSimple
+from ccpi.optimisation.algs 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
+
+# Scale and negative-log transform
+data.array = -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.geometry.pixel_num_h = data.geometry.pixel_num_h + cor_pad
+data.array = data_pad
+
+# Choose the number of voxels to reconstruct onto as number of detector pixels
+N = data.geometry.pixel_num_h
+
+# Geometric magnification
+mag = (np.abs(data.geometry.dist_center_detector) + \
+ np.abs(data.geometry.dist_source_center)) / \
+ np.abs(data.geometry.dist_source_center)
+
+# Voxel size is detector pixel size divided by mag
+voxel_size_h = data.geometry.pixel_size_h / mag
+
+# Construct the appropriate ImageGeometry
+ig = ImageGeometry(voxel_num_x=N,
+ voxel_num_y=N,
+ voxel_num_z=data.geometry.pixel_num_v,
+ voxel_size_x=voxel_size_h,
+ voxel_size_y=voxel_size_h,
+ voxel_size_z=voxel_size_h)
+
+# 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
+del data
+
+# Extract the Acquisition geometry for setting up projector.
+ag = data2.geometry
+
+# Set up the Projector (AcquisitionModel) using ASTRA on GPU
+Aop = AstraProjector3DSimple(ig, ag)
+
+# So and show simple backprojection
+z = Aop.adjoint(data2)
+plt.figure()
+plt.imshow(z.subset(horizontal_x=1000).as_array())
+plt.show()
+plt.figure()
+plt.imshow(z.subset(horizontal_y=1000).as_array())
+plt.show()
+plt.figure()
+plt.imshow(z.subset(vertical=100).array)
+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)
+
+# Display ortho slices of reconstruction
+plt.figure()
+plt.imshow(x.subset(horizontal_x=1000).as_array())
+plt.show()
+plt.figure()
+plt.imshow(x.subset(horizontal_y=1000).as_array())
+plt.show()
+plt.figure()
+plt.imshow(x.subset(vertical=100).as_array())
+plt.show()
diff --git a/Wrappers/Python/wip/work_out_adjoint.py b/Wrappers/Python/wip/work_out_adjoint.py
new file mode 100644
index 0000000..34d58ff
--- /dev/null
+++ b/Wrappers/Python/wip/work_out_adjoint.py
@@ -0,0 +1,115 @@
+
+# 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.ops 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
new file mode 100644
index 0000000..162e55a
--- /dev/null
+++ b/Wrappers/Python/wip/work_out_adjoint3D.py
@@ -0,0 +1,131 @@
+
+# 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.ops 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
new file mode 100644
index 0000000..2492826
--- /dev/null
+++ b/Wrappers/Python/wip/work_out_adjoint_sophiabeads.py
@@ -0,0 +1,115 @@
+
+# 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.ops 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