summaryrefslogtreecommitdiffstats
path: root/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers/Python/wip/demo_astra_sophiabeads3D.py')
-rw-r--r--Wrappers/Python/wip/demo_astra_sophiabeads3D.py103
1 files changed, 0 insertions, 103 deletions
diff --git a/Wrappers/Python/wip/demo_astra_sophiabeads3D.py b/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
deleted file mode 100644
index 624f121..0000000
--- a/Wrappers/Python/wip/demo_astra_sophiabeads3D.py
+++ /dev/null
@@ -1,103 +0,0 @@
-
-# This demo shows how to load a Nikon XTek micro-CT data set and reconstruct
-# a volume of 200 slices using the CGLS method. The SophiaBeads dataset with 64
-# projections is used as test data and can be obtained from here:
-# https://zenodo.org/record/16474
-# The filename with full path to the .xtekct file should be given as string
-# input to NikonDataReader to load in the data.
-
-# Do all imports
-import numpy as np
-import matplotlib.pyplot as plt
-from ccpi.io import NikonDataReader
-from ccpi.framework import ImageGeometry, ImageData
-from ccpi.astra.operators import AstraProjector3DSimple
-from ccpi.optimisation.algorithms import CGLS
-
-# Set up reader object and read in 200 central slices of the data
-datareader= NikonDataReader(xtek_file="REPLACE_THIS_BY_PATH_TO_DATASET/SophiaBeads_64_averaged.xtekct",
- roi=[(901,1101),(0,2000)])
-data = datareader.load_projections()
-
-# Scale and negative-log transform
-data.fill(-np.log(data.as_array()/60000.0))
-
-# Apply centering correction by zero padding, amount found manually
-cor_pad = 30
-data_pad = np.zeros((data.shape[0],data.shape[1],data.shape[2]+cor_pad))
-data_pad[:,:,cor_pad:] = data.array
-data.geometry.pixel_num_h = data.geometry.pixel_num_h + cor_pad
-data.array = data_pad
-
-# 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*np.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)
-
-# Do 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)
-
-# Set tolerance and number of iterations for reconstruction algorithms.
-opt = {'tol': 1e-4, 'iter': 30}
-
-# Run a CGLS reconstruction can be done:
-CGLS_alg = CGLS()
-CGLS_alg.set_up(x_init, Aop, data2)
-CGLS_alg.max_iteration = 2000
-CGLS_alg.run(opt['iter'])
-
-x_CGLS = CGLS_alg.get_output()
-
-# Display ortho slices of reconstruction
-plt.figure()
-plt.imshow(x_CGLS.subset(horizontal_x=1000).as_array())
-plt.show()
-plt.figure()
-plt.imshow(x_CGLS.subset(horizontal_y=1000).as_array())
-plt.show()
-plt.figure()
-plt.imshow(x_CGLS.subset(vertical=100).as_array())
-plt.show()
-
-plt.figure()
-plt.semilogy(CGLS_alg.objective)
-plt.title('CGLS criterion')
-plt.show() \ No newline at end of file