diff options
Diffstat (limited to 'Wrappers')
-rwxr-xr-x | Wrappers/Python/ccpi/astra/ops.py | 56 | ||||
-rw-r--r-- | Wrappers/Python/ccpi/astra/processors.py | 168 | ||||
-rwxr-xr-x | Wrappers/Python/ccpi/astra/utils.py | 9 | ||||
-rwxr-xr-x | Wrappers/Python/wip/demo_astra_mc.py | 4 | ||||
-rw-r--r-- | Wrappers/Python/wip/demo_astra_nexus.py | 171 | ||||
-rwxr-xr-x | Wrappers/Python/wip/demo_astra_simple.py | 4 | ||||
-rwxr-xr-x | Wrappers/Python/wip/demo_astra_sophiabeads.py | 10 | ||||
-rw-r--r-- | Wrappers/Python/wip/demo_astra_sophiabeads3D.py | 98 | ||||
-rw-r--r-- | Wrappers/Python/wip/work_out_adjoint.py | 115 | ||||
-rw-r--r-- | Wrappers/Python/wip/work_out_adjoint3D.py | 131 | ||||
-rw-r--r-- | Wrappers/Python/wip/work_out_adjoint_sophiabeads.py | 115 |
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 |