From 78a23cff0f17ed7edfcef1bee1048e9df11a8be5 Mon Sep 17 00:00:00 2001
From: jakobsj <jakobsj@users.noreply.github.com>
Date: Fri, 1 Jun 2018 14:08:16 +0100
Subject: Astra nexus (#9)
* Add 3D ASTRA operator and nexus demo
* Add astra sophiabeads 3d demo
* 2D astra scaling working simple case, not sophiabeads
* Simplied bp scaling, now seems to work 2D incl sophiabeads
* MC correct astra fp and bp scaling
* Fixed also ASTRA BP scaling in 3D, works sophiabeads 3D
* Minor tidy up
---
Wrappers/Python/wip/demo_astra_mc.py | 4 +-
Wrappers/Python/wip/demo_astra_nexus.py | 171 +++++++++++++++++++++
Wrappers/Python/wip/demo_astra_simple.py | 4 +-
Wrappers/Python/wip/demo_astra_sophiabeads.py | 10 +-
Wrappers/Python/wip/demo_astra_sophiabeads3D.py | 98 ++++++++++++
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, 638 insertions(+), 10 deletions(-)
create mode 100644 Wrappers/Python/wip/demo_astra_nexus.py
create mode 100644 Wrappers/Python/wip/demo_astra_sophiabeads3D.py
create mode 100644 Wrappers/Python/wip/work_out_adjoint.py
create mode 100644 Wrappers/Python/wip/work_out_adjoint3D.py
create mode 100644 Wrappers/Python/wip/work_out_adjoint_sophiabeads.py
(limited to 'Wrappers/Python/wip')
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
--
cgit v1.2.3