In [ ]:
%matplotlib inline
In [ ]:
# coding: utf-8
import datetime
import time
import sys
import astra
import numpy as np
import pylab as plt
In [ ]:
def create_test_cube(size):
'''
Create a simple hollow cube phantom
:param size: Cube size length (int >0)
:return: numpy mmaped array cube
'''
# cube = np.zeros((size,size,size), dtype='float32')
cube = np.memmap('/diskmnt/fast/makov/cube.tmp',
dtype='float32', mode='w+', shape=(size, size, size))
x0 = int(128.*size/1024)
x1 = int(895.*size/1024)
y0 = int(256.*size/1024)
y1 = int(767.*size/1024)
cube[x0:x1, x0:x1, x0:x1] = 1
cube[y0:y1, y0:y1, y0:y1, ] = 0
return cube
def build_geometry(angles, rec_size):
"""
Build astra circular geometry for square detector. See example #5 from python samples
:params angles_count:List of angles in radians
:param rec_size: Size of reconstruction volume
"""
# All distances in [pixels] (MMC1 dataset configuration)
pixel_size = 2.82473e-3
os_distance = 56.135 / pixel_size # object-sample distance
ds_distance = 225.082 / pixel_size # detector-sample distance
detector_size = rec_size
# Circular
# Parameters: width of detector column, height of detector row, #rows, #columns,
# angles, distance source-origin, distance origin-detector
proj_geom = astra.create_proj_geom('cone', ds_distance / os_distance, ds_distance / os_distance,
detector_size, detector_size, angles,
os_distance, (ds_distance - os_distance))
return proj_geom
def create_sino3d_gpu_mm(data, proj_geom, vol_geom, returnData=True, gpuIndex=None, sinogram_volume=None):
"""Create a forward projection of an image (3D).
:param data: Image data or ID.
:type data: :class:`numpy.ndarray` or :class:`int`
:param proj_geom: Projection geometry.
:type proj_geom: :class:`dict`
:param vol_geom: Volume geometry.
:type vol_geom: :class:`dict`
:param returnData: If False, only return the ID of the forward projection.
:type returnData: :class:`bool`
:param gpuIndex: Optional GPU index.
:type gpuIndex: :class:`int`
:returns: :class:`int` or (:class:`int`, :class:`numpy.ndarray`) -- If ``returnData=False``, returns the ID of the forward projection. Otherwise, returns a tuple containing the ID of the forward projection and the forward projection itself, in that order.
"""
if isinstance(data, np.ndarray):
volume_id = astra.data3d.link('-vol', vol_geom, data)
else:
volume_id = data
if isinstance(sinogram_volume, np.ndarray):
sino_id = astra.data3d.link('-sino', proj_geom, sinogram_volume)
else:
sino_id = astra.data3d.create('-sino', proj_geom, 0)
algString = 'FP3D_CUDA'
cfg = astra.astra_dict(algString)
if gpuIndex is not None:
cfg['option'] = {'GPUindex': gpuIndex}
cfg['ProjectionDataId'] = sino_id
cfg['VolumeDataId'] = volume_id
alg_id = astra.algorithm.create(cfg)
astra.algorithm.run(alg_id)
astra.algorithm.delete(alg_id)
if isinstance(data, np.ndarray):
astra.data3d.delete(volume_id)
if returnData:
return sino_id, astra.data3d.get(sino_id)
else:
return sino_id
def test_reconstruction_cone_fdk(size, angles_count,
algorithm ='FDK_CUDA', iterations=100, gpus_list=[0, 1]):
"""
Build test object, sinogram and reconstrucion
See http://www.astra-toolbox.com/docs/index.html for other odcs
:param size: Size of reconstruction volume
:params angles_count: Number of angles in range 0-PI
:params gpus_list: list of GPU to use for reconstruction, usually [0,] or [1,] or [0,1]
"""
# Set up multi-GPU usage.
# This only works for 3D GPU forward projection and back projection.
astra.astra.set_gpu_index(gpus_list)
angles = np.linspace(0, np.pi, angles_count, False)
# angles = angles-(angles.max()+angles.min())/2
# print('Creating phantom')
t = time.time()
# Create a simple hollow cube phantom
cube = create_test_cube(size)
# Create projection geometry
proj_geom = build_geometry(angles, cube.shape[0])
# Create volume desscription
vol_geom = astra.create_vol_geom(cube.shape)
# Upload data to ASTRA
cube_id = astra.data3d.link('-vol', vol_geom, cube)
# print(time.time() - t)
# print('Creating sinogram')
# t = time.time()
# Create mmapfile for sinogram
sinogram_volume = np.memmap('/diskmnt/fast/makov/sinogram.tmp',
dtype='float32', mode='w+',
shape=(cube.shape[0], len(angles), cube.shape[1])) # TODO: check for non cube
# Create sinogram from cube and NOT send data to scipt
proj_id = create_sino3d_gpu_mm(cube_id, proj_geom, vol_geom,
returnData=False, sinogram_volume=sinogram_volume)
astra.data3d.delete(cube_id)
del cube
# print(time.time() - t)
# print('Creating volume')
# t = time.time()
# To get data by id
# projection_data = astra.data3d.get(proj_id)
rec_volume = np.memmap('/diskmnt/fast/makov/rec.tmp',
dtype='float32', mode='w+',
shape=(vol_geom['GridRowCount'], # TODO: check for non cube
vol_geom['GridColCount'],
vol_geom['GridSliceCount'])
)
rec_id = astra.data3d.link('-vol', vol_geom, rec_volume)
# Set up the parameters for a reconstruction algorithm using the GPU
# Complete list of suporteed algoritms can be found
# http://www.astra-toolbox.com/docs/algs/index.html
# print(time.time() - t)
# print('Start reconstruction')
cfg = astra.astra_dict(algorithm)
# cfg = astra.astra_dict('CGLS3D_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = proj_id
cfg['option'] = {'ShortScan': True}
# Create the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
t = time.time()
# For FDK number of iterationa my be 1
astra.algorithm.run(alg_id, iterations)
# print('Stop reconstruction')
rec_time = time.time() - t
# rec = astra.data3d.get(rec_id)
# plt.figure()
# plt.imshow(rec[:,:,int(rec.shape[-1]/2)], cmap=plt.cm.jet)
# plt.show()
# Clean up. Note that GPU memory is tied up in the algorithm object,
# and main RAM in the data objects.
# astra.data3d.info()
astra.algorithm.delete(alg_id)
astra.data3d.delete(proj_id)
astra.data3d.delete(rec_id)
del rec_volume
del sinogram_volume
return rec_time
In [ ]:
#warm up
rec_time = test_reconstruction_cone_fdk(
algorithm='FDK_CUDA',
size=300, angles_count=300, gpus_list=[0, 1])
for method in ['FDK_CUDA', 'CGLS3D_CUDA']:
rec_time = test_reconstruction_cone_fdk(
algorithm=method,
size=300, angles_count=300, gpus_list=[0, 1])
print(method, rec_time)
# load and show reconstruction cut
rec_volume = np.memmap('/diskmnt/fast/makov/rec.tmp',
dtype='float32', mode='r',
shape=(300, 300, 300)
)
plt.figure()
plt.imshow(rec_volume[:, :, int(rec_volume.shape[-1] / 2)],
cmap=plt.cm.gray_r)
plt.colorbar()
plt.savefig('{}_{}_{}.png'.format(method, astra.__version__, datetime.datetime.now()))
plt.show()
del rec_volume
with open('astra_bech.log', 'a') as f:
f.write('{}\t{}\t{}\n'.format(datetime.datetime.now(), method, rec_time))
In [ ]:
!less astra_bech.log
In [ ]:
phantom = create_test_cube(size=200)[128]
In [ ]:
#-----------------------------------------------------------------------
#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam
#
#Author: Daniel M. Pelt
#Contact: D.M.Pelt@cwi.nl
#Website: http://dmpelt.github.io/pyastratoolbox/
#
#
#This file is part of the Python interface to the
#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
#
#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#
#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#
#-----------------------------------------------------------------------
import astra
import numpy as np
pixel_size = 11.000435e-3
os_distance = 2*54. / pixel_size
ds_distance = 3*225.315 / pixel_size
vol_geom = astra.create_vol_geom(200, 200)
angles = np.linspace(0,np.pi,90,False)
# proj_geom = astra.create_proj_geom('parallel', 1.0, 384, angles)
proj_geom = astra.create_proj_geom('fanflat', ds_distance / os_distance, 384, angles,
os_distance, (ds_distance - os_distance))
# As before, create a sinogram from a phantom
import scipy.io
P = phantom
proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
sinogram_id, sinogram = astra.create_sino(P, proj_id)
import pylab
pylab.gray()
pylab.figure(1)
pylab.imshow(P)
plt.colorbar()
pylab.figure(2)
pylab.imshow(sinogram)
plt.colorbar()
# Create a data object for the reconstruction
rec_id = astra.data2d.create('-vol', vol_geom)
# create configuration
cfg = astra.astra_dict('FBP_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['FilterType'] = 'Ram-Lak'
# possible values for FilterType:
# none, ram-lak, shepp-logan, cosine, hamming, hann, tukey, lanczos,
# triangular, gaussian, barlett-hann, blackman, nuttall, blackman-harris,
# blackman-nuttall, flat-top, kaiser, parzen
# Create and run the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
astra.algorithm.run(alg_id,10)
# Get the result
rec = astra.data2d.get(rec_id)
k = sinogram.shape[0]*(os_distance/ds_distance)**2/(np.pi/2)
rec*=k
pylab.figure(3)
pylab.imshow(rec)
plt.colorbar()
pylab.show()
# Clean up. Note that GPU memory is tied up in the algorithm object,
# and main RAM in the data objects.
astra.algorithm.delete(alg_id)
astra.data2d.delete(rec_id)
astra.data2d.delete(sinogram_id)
astra.projector.delete(proj_id)
print(phantom.mean(), rec.mean(), 1./(phantom.mean()/ rec.mean()))
In [ ]: