From 78a23cff0f17ed7edfcef1bee1048e9df11a8be5 Mon Sep 17 00:00:00 2001 From: jakobsj 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