From 3cae1d138c53a3fd042de3d2c9d9a07cf0650e0f Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Tue, 24 Feb 2015 12:35:45 +0100 Subject: Added Python interface --- samples/python/phantom.mat | Bin 0 -> 5583 bytes samples/python/s001_sinogram_par2d.py | 60 ++++++++++++++++ samples/python/s002_data2d.py | 77 ++++++++++++++++++++ samples/python/s003_gpu_reconstruction.py | 75 ++++++++++++++++++++ samples/python/s004_cpu_reconstruction.py | 81 +++++++++++++++++++++ samples/python/s005_3d_geometry.py | 114 ++++++++++++++++++++++++++++++ samples/python/s006_3d_data.py | 76 ++++++++++++++++++++ samples/python/s007_3d_reconstruction.py | 77 ++++++++++++++++++++ samples/python/s008_gpu_selection.py | 61 ++++++++++++++++ samples/python/s009_projection_matrix.py | 65 +++++++++++++++++ samples/python/s010_supersampling.py | 85 ++++++++++++++++++++++ samples/python/s011_object_info.py | 54 ++++++++++++++ samples/python/s012_masks.py | 92 ++++++++++++++++++++++++ samples/python/s013_constraints.py | 77 ++++++++++++++++++++ samples/python/s014_FBP.py | 76 ++++++++++++++++++++ samples/python/s015_fp_bp.py | 86 ++++++++++++++++++++++ samples/python/s016_plots.py | 86 ++++++++++++++++++++++ 17 files changed, 1242 insertions(+) create mode 100644 samples/python/phantom.mat create mode 100644 samples/python/s001_sinogram_par2d.py create mode 100644 samples/python/s002_data2d.py create mode 100644 samples/python/s003_gpu_reconstruction.py create mode 100644 samples/python/s004_cpu_reconstruction.py create mode 100644 samples/python/s005_3d_geometry.py create mode 100644 samples/python/s006_3d_data.py create mode 100644 samples/python/s007_3d_reconstruction.py create mode 100644 samples/python/s008_gpu_selection.py create mode 100644 samples/python/s009_projection_matrix.py create mode 100644 samples/python/s010_supersampling.py create mode 100644 samples/python/s011_object_info.py create mode 100644 samples/python/s012_masks.py create mode 100644 samples/python/s013_constraints.py create mode 100644 samples/python/s014_FBP.py create mode 100644 samples/python/s015_fp_bp.py create mode 100644 samples/python/s016_plots.py (limited to 'samples/python') diff --git a/samples/python/phantom.mat b/samples/python/phantom.mat new file mode 100644 index 0000000..c465bbe Binary files /dev/null and b/samples/python/phantom.mat differ diff --git a/samples/python/s001_sinogram_par2d.py b/samples/python/s001_sinogram_par2d.py new file mode 100644 index 0000000..009d9b3 --- /dev/null +++ b/samples/python/s001_sinogram_par2d.py @@ -0,0 +1,60 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +# Create a basic 256x256 square volume geometry +vol_geom = astra.create_vol_geom(256, 256) + +# Create a parallel beam geometry with 180 angles between 0 and pi, and +# 384 detector pixels of width 1. +# For more details on available geometries, see the online help of the +# function astra_create_proj_geom . +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) + +# Load a 256x256 phantom image +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] + +# Create a sinogram using the GPU. +# Note that the first time the GPU is accessed, there may be a delay +# of up to 10 seconds for initialization. +proj_id = astra.create_projector('line',proj_geom,vol_geom) +sinogram_id, sinogram = astra.create_sino(P, proj_id,useCUDA=True) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(P) +pylab.figure(2) +pylab.imshow(sinogram) +pylab.show() + + +# Free memory +astra.data2d.delete(sinogram_id) +astra.projector.delete(proj_id) diff --git a/samples/python/s002_data2d.py b/samples/python/s002_data2d.py new file mode 100644 index 0000000..35fb91f --- /dev/null +++ b/samples/python/s002_data2d.py @@ -0,0 +1,77 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(256, 256) + +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) + + +# Create volumes + +# initialized to zero +v0 = astra.data2d.create('-vol', vol_geom) + +# initialized to 3.0 +v1 = astra.data2d.create('-vol', vol_geom, 3.0) + +# initialized to a matrix. A may be a single, double or logical (0/1) array. +import scipy.io +A = scipy.io.loadmat('phantom.mat')['phantom256'] +v2 = astra.data2d.create('-vol', vol_geom, A) + + +# Projection data +s0 = astra.data2d.create('-sino', proj_geom) +# Initialization to a scalar or a matrix also works, exactly as with a volume. + + +# Update data + +# set to zero +astra.data2d.store(v0, 0) + +# set to a matrix +astra.data2d.store(v2, A) + + + +# Retrieve data + +R = astra.data2d.get(v2) +import pylab +pylab.gray() +pylab.imshow(R) +pylab.show() + + +# Free memory +astra.data2d.delete(v0) +astra.data2d.delete(v1) +astra.data2d.delete(v2) +astra.data2d.delete(s0) diff --git a/samples/python/s003_gpu_reconstruction.py b/samples/python/s003_gpu_reconstruction.py new file mode 100644 index 0000000..4f6ec1f --- /dev/null +++ b/samples/python/s003_gpu_reconstruction.py @@ -0,0 +1,75 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) + +# As before, create a sinogram from a phantom +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] +proj_id = astra.create_projector('line',proj_geom,vol_geom) +sinogram_id, sinogram = astra.create_sino(P, proj_id,useCUDA=True) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(P) +pylab.figure(2) +pylab.imshow(sinogram) + +# Create a data object for the reconstruction +rec_id = astra.data2d.create('-vol', vol_geom) + +# Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra.astra_dict('SIRT_CUDA') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = sinogram_id + +# Available algorithms: +# SIRT_CUDA, SART_CUDA, EM_CUDA, FBP_CUDA (see the FBP sample) + + +# Create the algorithm object from the configuration structure +alg_id = astra.algorithm.create(cfg) + +# Run 150 iterations of the algorithm +astra.algorithm.run(alg_id, 150) + +# Get the result +rec = astra.data2d.get(rec_id) +pylab.figure(3) +pylab.imshow(rec) +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) diff --git a/samples/python/s004_cpu_reconstruction.py b/samples/python/s004_cpu_reconstruction.py new file mode 100644 index 0000000..8385cf8 --- /dev/null +++ b/samples/python/s004_cpu_reconstruction.py @@ -0,0 +1,81 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) + +# For CPU-based algorithms, a "projector" object specifies the projection +# model used. In this case, we use the "strip" model. +proj_id = astra.create_projector('strip', proj_geom, vol_geom) + +# Create a sinogram from a phantom +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] +sinogram_id, sinogram = astra.create_sino(P, proj_id) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(P) +pylab.figure(2) +pylab.imshow(sinogram) + +# Create a data object for the reconstruction +rec_id = astra.data2d.create('-vol', vol_geom) + +# Set up the parameters for a reconstruction algorithm using the CPU +# The main difference with the configuration of a GPU algorithm is the +# extra ProjectorId setting. +cfg = astra.astra_dict('SIRT') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = sinogram_id +cfg['ProjectorId'] = proj_id + +# Available algorithms: +# ART, SART, SIRT, CGLS, FBP + + +# Create the algorithm object from the configuration structure +alg_id = astra.algorithm.create(cfg) + +# Run 20 iterations of the algorithm +# This will have a runtime in the order of 10 seconds. +astra.algorithm.run(alg_id, 20) + +# Get the result +rec = astra.data2d.get(rec_id) +pylab.figure(3) +pylab.imshow(rec) +pylab.show() + +# Clean up. +astra.algorithm.delete(alg_id) +astra.data2d.delete(rec_id) +astra.data2d.delete(sinogram_id) +astra.projector.delete(proj_id) diff --git a/samples/python/s005_3d_geometry.py b/samples/python/s005_3d_geometry.py new file mode 100644 index 0000000..f43fc7e --- /dev/null +++ b/samples/python/s005_3d_geometry.py @@ -0,0 +1,114 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +from six.moves import range +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(64, 64, 64) + + +# There are two main 3d projection geometry types: cone beam and parallel beam. +# Each has a regular variant, and a 'vec' variant. +# The 'vec' variants are completely free in the placement of source/detector, +# while the regular variants assume circular trajectories around the z-axis. + + +# ------------- +# Parallel beam +# ------------- + + +# Circular + +# Parameters: width of detector column, height of detector row, #rows, #columns +angles = np.linspace(0, 2*np.pi, 48, False) +proj_geom = astra.create_proj_geom('parallel3d', 1.0, 1.0, 32, 64, angles) + + +# Free + +# We generate the same geometry as the circular one above. +vectors = np.zeros((len(angles), 12)) +for i in range(len(angles)): + # ray direction + vectors[i,0] = np.sin(angles[i]) + vectors[i,1] = -np.cos(angles[i]) + vectors[i,2] = 0 + + # center of detector + vectors[i,3:6] = 0 + + # vector from detector pixel (0,0) to (0,1) + vectors[i,6] = np.cos(angles[i]) + vectors[i,7] = np.sin(angles[i]) + vectors[i,8] = 0; + + # vector from detector pixel (0,0) to (1,0) + vectors[i,9] = 0 + vectors[i,10] = 0 + vectors[i,11] = 1 + +# Parameters: #rows, #columns, vectors +proj_geom = astra.create_proj_geom('parallel3d_vec', 32, 64, vectors) + +# ---------- +# Cone beam +# ---------- + + +# Circular + +# Parameters: width of detector column, height of detector row, #rows, #columns, +# angles, distance source-origin, distance origin-detector +angles = np.linspace(0, 2*np.pi, 48, False) +proj_geom = astra.create_proj_geom('cone', 1.0, 1.0, 32, 64, angles, 1000, 0) + +# Free + +vectors = np.zeros((len(angles), 12)) +for i in range(len(angles)): + # source + vectors[i,0] = np.sin(angles[i]) * 1000 + vectors[i,1] = -np.cos(angles[i]) * 1000 + vectors[i,2] = 0 + + # center of detector + vectors[i,3:6] = 0 + + # vector from detector pixel (0,0) to (0,1) + vectors[i,6] = np.cos(angles[i]) + vectors[i,7] = np.sin(angles[i]) + vectors[i,8] = 0 + + # vector from detector pixel (0,0) to (1,0) + vectors[i,9] = 0 + vectors[i,10] = 0 + vectors[i,11] = 1 + +# Parameters: #rows, #columns, vectors +proj_geom = astra.create_proj_geom('cone_vec', 32, 64, vectors) + diff --git a/samples/python/s006_3d_data.py b/samples/python/s006_3d_data.py new file mode 100644 index 0000000..5178179 --- /dev/null +++ b/samples/python/s006_3d_data.py @@ -0,0 +1,76 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +# Create a 3D volume geometry. +# Parameter order: rows, colums, slices (y, x, z) +vol_geom = astra.create_vol_geom(64, 48, 32) + + +# Create volumes + +# initialized to zero +v0 = astra.data3d.create('-vol', vol_geom) + +# initialized to 3.0 +v1 = astra.data3d.create('-vol', vol_geom, 3.0) + +# initialized to a matrix. A may be a single or double array. +# Coordinate order: slice, row, column (z, y, x) +A = np.zeros((32, 64, 48)) +v2 = astra.data3d.create('-vol', vol_geom, A) + + +# Projection data + +# 2 projection directions, along x and y axis resp. +V = np.array([[ 1,0,0, 0,0,0, 0,1,0, 0,0,1], + [0,1,0, 0,0,0, -1,0,0, 0,0,1]],dtype=np.float) +# 32 rows (v), 64 columns (u) +proj_geom = astra.create_proj_geom('parallel3d_vec', 32, 64, V) + +s0 = astra.data3d.create('-proj3d', proj_geom) + +# Initialization to a scalar or zero works exactly as with a volume. + +# Initialized to a matrix: +# Coordinate order: row (v), angle, column (u) +A = np.zeros((32, 2, 64)) +s1 = astra.data3d.create('-proj3d', proj_geom, A) + + +# Retrieve data: +R = astra.data3d.get(v1) + + +# Delete all created data objects +astra.data3d.delete(v0) +astra.data3d.delete(v1) +astra.data3d.delete(v2) +astra.data3d.delete(s0) +astra.data3d.delete(s1) diff --git a/samples/python/s007_3d_reconstruction.py b/samples/python/s007_3d_reconstruction.py new file mode 100644 index 0000000..40e9556 --- /dev/null +++ b/samples/python/s007_3d_reconstruction.py @@ -0,0 +1,77 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(128, 128, 128) + +angles = np.linspace(0, np.pi, 180,False) +proj_geom = astra.create_proj_geom('parallel3d', 1.0, 1.0, 128, 192, angles) + +# Create a simple hollow cube phantom +cube = np.zeros((128,128,128)) +cube[17:113,17:113,17:113] = 1 +cube[33:97,33:97,33:97] = 0 + +# Create projection data from this +proj_id, proj_data = astra.create_sino3d_gpu(cube, proj_geom, vol_geom) + +# Display a single projection image +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(proj_data[:,20,:]) + +# Create a data object for the reconstruction +rec_id = astra.data3d.create('-vol', vol_geom) + +# Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra.astra_dict('SIRT3D_CUDA') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = proj_id + + +# Create the algorithm object from the configuration structure +alg_id = astra.algorithm.create(cfg) + +# Run 150 iterations of the algorithm +# Note that this requires about 750MB of GPU memory, and has a runtime +# in the order of 10 seconds. +astra.algorithm.run(alg_id, 150) + +# Get the result +rec = astra.data3d.get(rec_id) +pylab.figure(2) +pylab.imshow(rec[:,:,65]) +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.data3d.delete(rec_id) +astra.data3d.delete(proj_id) diff --git a/samples/python/s008_gpu_selection.py b/samples/python/s008_gpu_selection.py new file mode 100644 index 0000000..c42e53b --- /dev/null +++ b/samples/python/s008_gpu_selection.py @@ -0,0 +1,61 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] + +proj_id = astra.create_projector('line',proj_geom,vol_geom) + +# Create a sinogram from a phantom, using GPU #1. (The default is #0) +sinogram_id, sinogram = astra.create_sino(P, proj_id, useCUDA=True, gpuIndex=1) + + +# Set up the parameters for a reconstruction algorithm using the GPU +rec_id = astra.data2d.create('-vol', vol_geom) +cfg = astra.astra_dict('SIRT_CUDA') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = sinogram_id + +# Use GPU #1 for the reconstruction. (The default is #0.) +cfg['option'] = {} +cfg['option']['GPUindex'] = 1 + +# Run 150 iterations of the algorithm +alg_id = astra.algorithm.create(cfg) +astra.algorithm.run(alg_id, 150) +rec = astra.data2d.get(rec_id) + + +# Clean up. +astra.algorithm.delete(alg_id) +astra.data2d.delete(rec_id) +astra.data2d.delete(sinogram_id) +astra.projector.delete(proj_id) diff --git a/samples/python/s009_projection_matrix.py b/samples/python/s009_projection_matrix.py new file mode 100644 index 0000000..c4c4557 --- /dev/null +++ b/samples/python/s009_projection_matrix.py @@ -0,0 +1,65 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) + +# For CPU-based algorithms, a "projector" object specifies the projection +# model used. In this case, we use the "line" model. +proj_id = astra.create_projector('line', proj_geom, vol_geom) + +# Generate the projection matrix for this projection model. +# This creates a matrix W where entry w_{i,j} corresponds to the +# contribution of volume element j to detector element i. +matrix_id = astra.projector.matrix(proj_id) + +# Get the projection matrix as a Scipy sparse matrix. +W = astra.matrix.get(matrix_id) + + +# Manually use this projection matrix to do a projection: +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] +s = W.dot(P.flatten()) +s = np.reshape(s, (len(proj_geom['ProjectionAngles']),proj_geom['DetectorCount'])) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(s) +pylab.show() + +# Each row of the projection matrix corresponds to a detector element. +# Detector t for angle p is for row 1 + t + p*proj_geom.DetectorCount. +# Each column corresponds to a volume pixel. +# Pixel (x,y) corresponds to column 1 + x + y*vol_geom.GridColCount. + + +astra.projector.delete(proj_id) +astra.matrix.delete(matrix_id) diff --git a/samples/python/s010_supersampling.py b/samples/python/s010_supersampling.py new file mode 100644 index 0000000..1a337bc --- /dev/null +++ b/samples/python/s010_supersampling.py @@ -0,0 +1,85 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 3.0, 128, np.linspace(0,np.pi,180,False)) +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] + +# Because the astra.create_sino method does not have support for +# all possible algorithm options, we manually create a sinogram +phantom_id = astra.data2d.create('-vol', vol_geom, P) +sinogram_id = astra.data2d.create('-sino', proj_geom) +cfg = astra.astra_dict('FP_CUDA') +cfg['VolumeDataId'] = phantom_id +cfg['ProjectionDataId'] = sinogram_id + +# Set up 3 rays per detector element +cfg['option'] = {} +cfg['option']['DetectorSuperSampling'] = 3 + +alg_id = astra.algorithm.create(cfg) +astra.algorithm.run(alg_id) +astra.algorithm.delete(alg_id) +astra.data2d.delete(phantom_id) + +sinogram3 = astra.data2d.get(sinogram_id) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(P) +pylab.figure(2) +pylab.imshow(sinogram3) + +# Create a reconstruction, also using supersampling +rec_id = astra.data2d.create('-vol', vol_geom) +cfg = astra.astra_dict('SIRT_CUDA') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = sinogram_id +# Set up 3 rays per detector element +cfg['option'] = {} +cfg['option']['DetectorSuperSampling'] = 3 + +# There is also an option for supersampling during the backprojection step. +# This should be used if your detector pixels are smaller than the voxels. + +# Set up 2 rays per image pixel dimension, for 4 rays total per image pixel. +# cfg['option']['PixelSuperSampling'] = 2 + + +alg_id = astra.algorithm.create(cfg) +astra.algorithm.run(alg_id, 150) +astra.algorithm.delete(alg_id) + +rec = astra.data2d.get(rec_id) +pylab.figure(3) +pylab.imshow(rec) +pylab.show() + diff --git a/samples/python/s011_object_info.py b/samples/python/s011_object_info.py new file mode 100644 index 0000000..02f387a --- /dev/null +++ b/samples/python/s011_object_info.py @@ -0,0 +1,54 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra + +# Create two volume geometries +vol_geom1 = astra.create_vol_geom(256, 256) +vol_geom2 = astra.create_vol_geom(512, 256) + +# Create volumes +v0 = astra.data2d.create('-vol', vol_geom1) +v1 = astra.data2d.create('-vol', vol_geom2) +v2 = astra.data2d.create('-vol', vol_geom2) + +# Show the currently allocated volumes +astra.data2d.info() + + +astra.data2d.delete(v2) +astra.data2d.info() + +astra.data2d.clear() +astra.data2d.info() + + + +# The same clear and info command also work for other object types: +astra.algorithm.info() +astra.data3d.info() +astra.projector.info() +astra.matrix.info() diff --git a/samples/python/s012_masks.py b/samples/python/s012_masks.py new file mode 100644 index 0000000..441d11b --- /dev/null +++ b/samples/python/s012_masks.py @@ -0,0 +1,92 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + + +import astra +import numpy as np + +# In this example we will create a reconstruction in a circular region, +# instead of the usual rectangle. + +# This is done by placing a circular mask on the square reconstruction volume: + +c = np.linspace(-127.5,127.5,256) +x, y = np.meshgrid(c,c) +mask = np.array((x**2 + y**2 < 127.5**2),dtype=np.float) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(mask) + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,50,False)) + +# As before, create a sinogram from a phantom +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] +proj_id = astra.create_projector('line',proj_geom,vol_geom) +sinogram_id, sinogram = astra.create_sino(P, proj_id,useCUDA=True) + +pylab.figure(2) +pylab.imshow(P) +pylab.figure(3) +pylab.imshow(sinogram) + +# Create a data object for the reconstruction +rec_id = astra.data2d.create('-vol', vol_geom) + +# Create a data object for the mask +mask_id = astra.data2d.create('-vol', vol_geom, mask) + +# Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra.astra_dict('SIRT_CUDA') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = sinogram_id +cfg['option'] = {} +cfg['option']['ReconstructionMaskId'] = mask_id + +# Create the algorithm object from the configuration structure +alg_id = astra.algorithm.create(cfg) + +# Run 150 iterations of the algorithm +astra.algorithm.run(alg_id, 150) + +# Get the result +rec = astra.data2d.get(rec_id) + +pylab.figure(4) +pylab.imshow(rec) + +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(mask_id) +astra.data2d.delete(rec_id) +astra.data2d.delete(sinogram_id) +astra.projector.delete(proj_id) diff --git a/samples/python/s013_constraints.py b/samples/python/s013_constraints.py new file mode 100644 index 0000000..009360e --- /dev/null +++ b/samples/python/s013_constraints.py @@ -0,0 +1,77 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +# In this example we will create a reconstruction constrained to +# greyvalues between 0 and 1 + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,50,False)) + +# As before, create a sinogram from a phantom +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] +proj_id = astra.create_projector('line',proj_geom,vol_geom) +sinogram_id, sinogram = astra.create_sino(P, proj_id,useCUDA=True) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(P) +pylab.figure(2) +pylab.imshow(sinogram) + +# Create a data object for the reconstruction +rec_id = astra.data2d.create('-vol', vol_geom) + +# Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra.astra_dict('SIRT_CUDA') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = sinogram_id +cfg['option']={} +cfg['option']['MinConstraint'] = 0 +cfg['option']['MaxConstraint'] = 1 + +# Create the algorithm object from the configuration structure +alg_id = astra.algorithm.create(cfg) + +# Run 150 iterations of the algorithm +astra.algorithm.run(alg_id, 150) + +# Get the result +rec = astra.data2d.get(rec_id) +pylab.figure(3) +pylab.imshow(rec) +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) diff --git a/samples/python/s014_FBP.py b/samples/python/s014_FBP.py new file mode 100644 index 0000000..ef4afc2 --- /dev/null +++ b/samples/python/s014_FBP.py @@ -0,0 +1,76 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +import astra +import numpy as np + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) + +# As before, create a sinogram from a phantom +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] +proj_id = astra.create_projector('line',proj_geom,vol_geom) +sinogram_id, sinogram = astra.create_sino(P, proj_id,useCUDA=True) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(P) +pylab.figure(2) +pylab.imshow(sinogram) + +# 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) + +# Get the result +rec = astra.data2d.get(rec_id) +pylab.figure(3) +pylab.imshow(rec) +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) diff --git a/samples/python/s015_fp_bp.py b/samples/python/s015_fp_bp.py new file mode 100644 index 0000000..10c238d --- /dev/null +++ b/samples/python/s015_fp_bp.py @@ -0,0 +1,86 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + + +# This example demonstrates using the FP and BP primitives with Matlab's lsqr +# solver. Calls to FP (astra_create_sino_cuda) and +# BP (astra_create_backprojection_cuda) are wrapped in a function astra_wrap, +# and a handle to this function is passed to lsqr. + +# Because in this case the inputs/outputs of FP and BP have to be vectors +# instead of images (matrices), the calls require reshaping to and from vectors. + +import astra +import numpy as np + +# FP/BP wrapper class +class astra_wrap(object): + def __init__(self,proj_geom,vol_geom): + self.proj_id = astra.create_projector('line',proj_geom,vol_geom) + self.shape = (proj_geom['DetectorCount']*len(proj_geom['ProjectionAngles']),vol_geom['GridColCount']*vol_geom['GridRowCount']) + self.dtype = np.float + + def matvec(self,v): + sid, s = astra.create_sino(np.reshape(v,(vol_geom['GridRowCount'],vol_geom['GridColCount'])),self.proj_id,useCUDA=True) + astra.data2d.delete(sid) + return s.flatten() + + def rmatvec(self,v): + bid, b = astra.create_backprojection(np.reshape(v,(len(proj_geom['ProjectionAngles']),proj_geom['DetectorCount'],)),self.proj_id,useCUDA=True) + astra.data2d.delete(bid) + return b.flatten() + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) + +# Create a 256x256 phantom image +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] + +# Create a sinogram using the GPU. +proj_id = astra.create_projector('line',proj_geom,vol_geom) +sinogram_id, sinogram = astra.create_sino(P, proj_id,useCUDA=True) + +# Reshape the sinogram into a vector +b = sinogram.flatten() + +# Call lsqr with ASTRA FP and BP +import scipy.sparse.linalg +wrapper = astra_wrap(proj_geom,vol_geom) +result = scipy.sparse.linalg.lsqr(wrapper,b,atol=1e-4,btol=1e-4,iter_lim=25) + +# Reshape the result into an image +Y = np.reshape(result[0],(vol_geom['GridRowCount'], vol_geom['GridColCount'])); + +import pylab +pylab.gray() +pylab.imshow(Y) +pylab.show() + +astra.data2d.delete(sinogram_id) +astra.projector.delete(proj_id) +astra.projector.delete(wrapper.proj_id) + diff --git a/samples/python/s016_plots.py b/samples/python/s016_plots.py new file mode 100644 index 0000000..cd4d98c --- /dev/null +++ b/samples/python/s016_plots.py @@ -0,0 +1,86 @@ +#----------------------------------------------------------------------- +#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 . +# +#----------------------------------------------------------------------- + +from six.moves import range +import astra +import numpy as np + + +vol_geom = astra.create_vol_geom(256, 256) +proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False)) + +# As before, create a sinogram from a phantom +import scipy.io +P = scipy.io.loadmat('phantom.mat')['phantom256'] +proj_id = astra.create_projector('line',proj_geom,vol_geom) +sinogram_id, sinogram = astra.create_sino(P, proj_id,useCUDA=True) + +import pylab +pylab.gray() +pylab.figure(1) +pylab.imshow(P) +pylab.figure(2) +pylab.imshow(sinogram) + +# Create a data object for the reconstruction +rec_id = astra.data2d.create('-vol', vol_geom) + +# Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra.astra_dict('SIRT_CUDA') +cfg['ReconstructionDataId'] = rec_id +cfg['ProjectionDataId'] = sinogram_id + +# Create the algorithm object from the configuration structure +alg_id = astra.algorithm.create(cfg) + +# Run 1500 iterations of the algorithm one at a time, keeping track of errors +nIters = 1500 +phantom_error = np.zeros(nIters) +residual_error = np.zeros(nIters) +for i in range(nIters): + # Run a single iteration + astra.algorithm.run(alg_id, 1) + residual_error[i] = astra.algorithm.get_res_norm(alg_id) + rec = astra.data2d.get(rec_id) + phantom_error[i] = np.sqrt(((rec - P)**2).sum()) + +# Get the result +rec = astra.data2d.get(rec_id) +pylab.figure(3) +pylab.imshow(rec) + +pylab.figure(4) +pylab.plot(residual_error) +pylab.figure(5) +pylab.plot(phantom_error) + +pylab.show() + +# Clean up. +astra.algorithm.delete(alg_id) +astra.data2d.delete(rec_id) +astra.data2d.delete(sinogram_id) +astra.projector.delete(proj_id) -- cgit v1.2.3