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/matlab/s001_sinogram_par2d.m | 33 +++++++++ samples/matlab/s002_data2d.m | 60 ++++++++++++++++ samples/matlab/s003_gpu_reconstruction.m | 52 ++++++++++++++ samples/matlab/s004_cpu_reconstruction.m | 60 ++++++++++++++++ samples/matlab/s005_3d_geometry.m | 98 +++++++++++++++++++++++++ samples/matlab/s006_3d_data.m | 62 ++++++++++++++++ samples/matlab/s007_3d_reconstruction.m | 53 ++++++++++++++ samples/matlab/s008_gpu_selection.m | 37 ++++++++++ samples/matlab/s009_projection_matrix.m | 45 ++++++++++++ samples/matlab/s010_supersampling.m | 58 +++++++++++++++ samples/matlab/s011_object_info.m | 36 ++++++++++ samples/matlab/s012_masks.m | 60 ++++++++++++++++ samples/matlab/s013_constraints.m | 47 ++++++++++++ samples/matlab/s014_FBP.m | 47 ++++++++++++ samples/matlab/s015_fp_bp.m | 62 ++++++++++++++++ samples/matlab/s016_plots.m | 54 ++++++++++++++ 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 ++++++++++++++++++++++ samples/s001_sinogram_par2d.m | 33 --------- samples/s002_data2d.m | 60 ---------------- samples/s003_gpu_reconstruction.m | 52 -------------- samples/s004_cpu_reconstruction.m | 60 ---------------- samples/s005_3d_geometry.m | 98 ------------------------- samples/s006_3d_data.m | 62 ---------------- samples/s007_3d_reconstruction.m | 53 -------------- samples/s008_gpu_selection.m | 37 ---------- samples/s009_projection_matrix.m | 45 ------------ samples/s010_supersampling.m | 58 --------------- samples/s011_object_info.m | 36 ---------- samples/s012_masks.m | 60 ---------------- samples/s013_constraints.m | 47 ------------ samples/s014_FBP.m | 47 ------------ samples/s015_fp_bp.m | 62 ---------------- samples/s016_plots.m | 54 -------------- 49 files changed, 2106 insertions(+), 864 deletions(-) create mode 100644 samples/matlab/s001_sinogram_par2d.m create mode 100644 samples/matlab/s002_data2d.m create mode 100644 samples/matlab/s003_gpu_reconstruction.m create mode 100644 samples/matlab/s004_cpu_reconstruction.m create mode 100644 samples/matlab/s005_3d_geometry.m create mode 100644 samples/matlab/s006_3d_data.m create mode 100644 samples/matlab/s007_3d_reconstruction.m create mode 100644 samples/matlab/s008_gpu_selection.m create mode 100644 samples/matlab/s009_projection_matrix.m create mode 100644 samples/matlab/s010_supersampling.m create mode 100644 samples/matlab/s011_object_info.m create mode 100644 samples/matlab/s012_masks.m create mode 100644 samples/matlab/s013_constraints.m create mode 100644 samples/matlab/s014_FBP.m create mode 100644 samples/matlab/s015_fp_bp.m create mode 100644 samples/matlab/s016_plots.m 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 delete mode 100644 samples/s001_sinogram_par2d.m delete mode 100644 samples/s002_data2d.m delete mode 100644 samples/s003_gpu_reconstruction.m delete mode 100644 samples/s004_cpu_reconstruction.m delete mode 100644 samples/s005_3d_geometry.m delete mode 100644 samples/s006_3d_data.m delete mode 100644 samples/s007_3d_reconstruction.m delete mode 100644 samples/s008_gpu_selection.m delete mode 100644 samples/s009_projection_matrix.m delete mode 100644 samples/s010_supersampling.m delete mode 100644 samples/s011_object_info.m delete mode 100644 samples/s012_masks.m delete mode 100644 samples/s013_constraints.m delete mode 100644 samples/s014_FBP.m delete mode 100644 samples/s015_fp_bp.m delete mode 100644 samples/s016_plots.m (limited to 'samples') diff --git a/samples/matlab/s001_sinogram_par2d.m b/samples/matlab/s001_sinogram_par2d.m new file mode 100644 index 0000000..4494e7b --- /dev/null +++ b/samples/matlab/s001_sinogram_par2d.m @@ -0,0 +1,33 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +% 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, linspace2(0,pi,180)); + +% Create a 256x256 phantom image using matlab's built-in phantom() function +P = phantom(256); + +% 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. +[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); + +figure(1); imshow(P, []); +figure(2); imshow(sinogram, []); + + +% Free memory +astra_mex_data2d('delete', sinogram_id); diff --git a/samples/matlab/s002_data2d.m b/samples/matlab/s002_data2d.m new file mode 100644 index 0000000..a91071f --- /dev/null +++ b/samples/matlab/s002_data2d.m @@ -0,0 +1,60 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(256, 256); + +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); + + +% Create volumes + +% initialized to zero +v0 = astra_mex_data2d('create', '-vol', vol_geom); + +% initialized to 3.0 +v1 = astra_mex_data2d('create', '-vol', vol_geom, 3.0); + +% initialized to a matrix. A may be a single, double or logical (0/1) array. +A = phantom(256); +v2 = astra_mex_data2d('create', '-vol', vol_geom, A); + + +% Projection data +s0 = astra_mex_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_mex_data2d('store', v0, 0); + +% set to a matrix +astra_mex_data2d('store', v2, A); + + + +% Retrieve data + +R = astra_mex_data2d('get', v2); +imshow(R, []); + + +% Retrieve data as a single array. Since astra internally stores +% data as single precision, this is more efficient: + +R = astra_mex_data2d('get_single', v2); + + +% Free memory +astra_mex_data2d('delete', v0); +astra_mex_data2d('delete', v1); +astra_mex_data2d('delete', v2); +astra_mex_data2d('delete', s0); diff --git a/samples/matlab/s003_gpu_reconstruction.m b/samples/matlab/s003_gpu_reconstruction.m new file mode 100644 index 0000000..efb5c68 --- /dev/null +++ b/samples/matlab/s003_gpu_reconstruction.m @@ -0,0 +1,52 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); + +% As before, create a sinogram from a phantom +P = phantom(256); +[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); +figure(1); imshow(P, []); +figure(2); imshow(sinogram, []); + +astra_mex_data2d('delete', sinogram_id); + +% We now re-create the sinogram data object as we would do when loading +% an external sinogram +sinogram_id = astra_mex_data2d('create', '-sino', proj_geom, sinogram); + +% Create a data object for the reconstruction +rec_id = astra_mex_data2d('create', '-vol', vol_geom); + +% Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra_struct('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_mex_algorithm('create', cfg); + +% Run 150 iterations of the algorithm +astra_mex_algorithm('iterate', alg_id, 150); + +% Get the result +rec = astra_mex_data2d('get', rec_id); +figure(3); imshow(rec, []); + +% Clean up. Note that GPU memory is tied up in the algorithm object, +% and main RAM in the data objects. +astra_mex_algorithm('delete', alg_id); +astra_mex_data2d('delete', rec_id); +astra_mex_data2d('delete', sinogram_id); diff --git a/samples/matlab/s004_cpu_reconstruction.m b/samples/matlab/s004_cpu_reconstruction.m new file mode 100644 index 0000000..f25cd2b --- /dev/null +++ b/samples/matlab/s004_cpu_reconstruction.m @@ -0,0 +1,60 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); + +% 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 +P = phantom(256); +[sinogram_id, sinogram] = astra_create_sino(P, proj_id); +figure(1); imshow(P, []); +figure(2); imshow(sinogram, []); + +astra_mex_data2d('delete', sinogram_id); + +% We now re-create the sinogram data object as we would do when loading +% an external sinogram +sinogram_id = astra_mex_data2d('create', '-sino', proj_geom, sinogram); + +% Create a data object for the reconstruction +rec_id = astra_mex_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_struct('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_mex_algorithm('create', cfg); + +% Run 20 iterations of the algorithm +% This will have a runtime in the order of 10 seconds. +astra_mex_algorithm('iterate', alg_id, 20); + +% Get the result +rec = astra_mex_data2d('get', rec_id); +figure(3); imshow(rec, []); + +% Clean up. +astra_mex_projector('delete', proj_id); +astra_mex_algorithm('delete', alg_id); +astra_mex_data2d('delete', rec_id); +astra_mex_data2d('delete', sinogram_id); diff --git a/samples/matlab/s005_3d_geometry.m b/samples/matlab/s005_3d_geometry.m new file mode 100644 index 0000000..4b7360b --- /dev/null +++ b/samples/matlab/s005_3d_geometry.m @@ -0,0 +1,98 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +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 = linspace2(0, 2*pi, 48); +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 = zeros(numel(angles), 12); +for i = 1:numel(angles) + % ray direction + vectors(i,1) = sin(angles(i)); + vectors(i,2) = -cos(angles(i)); + vectors(i,3) = 0; + + % center of detector + vectors(i,4:6) = 0; + + % vector from detector pixel (0,0) to (0,1) + vectors(i,7) = cos(angles(i)); + vectors(i,8) = sin(angles(i)); + vectors(i,9) = 0; + + % vector from detector pixel (0,0) to (1,0) + vectors(i,10) = 0; + vectors(i,11) = 0; + vectors(i,12) = 1; +end + +% 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 = linspace2(0, 2*pi, 48); +proj_geom = astra_create_proj_geom('cone', 1.0, 1.0, 32, 64, ... + angles, 1000, 0); + +% Free + +vectors = zeros(numel(angles), 12); +for i = 1:numel(angles) + + % source + vectors(i,1) = sin(angles(i)) * 1000; + vectors(i,2) = -cos(angles(i)) * 1000; + vectors(i,3) = 0; + + % center of detector + vectors(i,4:6) = 0; + + % vector from detector pixel (0,0) to (0,1) + vectors(i,7) = cos(angles(i)); + vectors(i,8) = sin(angles(i)); + vectors(i,9) = 0; + + % vector from detector pixel (0,0) to (1,0) + vectors(i,10) = 0; + vectors(i,11) = 0; + vectors(i,12) = 1; +end + +% Parameters: #rows, #columns, vectors +proj_geom = astra_create_proj_geom('cone_vec', 32, 64, vectors); + diff --git a/samples/matlab/s006_3d_data.m b/samples/matlab/s006_3d_data.m new file mode 100644 index 0000000..32d88cc --- /dev/null +++ b/samples/matlab/s006_3d_data.m @@ -0,0 +1,62 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +% 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_mex_data3d('create', '-vol', vol_geom); + +% initialized to 3.0 +v1 = astra_mex_data3d('create', '-vol', vol_geom, 3.0); + +% initialized to a matrix. A may be a single or double array. +% Coordinate order: column, row, slice (x, y, z) +A = zeros(48, 64, 32); +v2 = astra_mex_data3d('create', '-vol', vol_geom, A); + + +% Projection data + +% 2 projection directions, along x and y axis resp. +V = [ 1 0 0 0 0 0 0 1 0 0 0 1 ; ... + 0 1 0 0 0 0 -1 0 0 0 0 1 ]; +% 32 rows (v), 64 columns (u) +proj_geom = astra_create_proj_geom('parallel3d_vec', 32, 64, V); + +s0 = astra_mex_data3d('create', '-proj3d', proj_geom); + +% Initialization to a scalar or zero works exactly as with a volume. + +% Initialized to a matrix: +% Coordinate order: column (u), angle, row (v) +A = zeros(64, 2, 32); +s1 = astra_mex_data3d('create', '-proj3d', proj_geom, A); + + +% Retrieve data: +R = astra_mex_data3d('get', v1); + +% Retrieve data as a single array. Since astra internally stores +% data as single precision, this is more efficient: +R = astra_mex_data3d('get_single', v1); + + + +% Delete all created data objects +astra_mex_data3d('delete', v0); +astra_mex_data3d('delete', v1); +astra_mex_data3d('delete', v2); +astra_mex_data3d('delete', s0); +astra_mex_data3d('delete', s1); diff --git a/samples/matlab/s007_3d_reconstruction.m b/samples/matlab/s007_3d_reconstruction.m new file mode 100644 index 0000000..fc9aca6 --- /dev/null +++ b/samples/matlab/s007_3d_reconstruction.m @@ -0,0 +1,53 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(128, 128, 128); + +angles = linspace2(0, pi, 180); +proj_geom = astra_create_proj_geom('parallel3d', 1.0, 1.0, 128, 192, angles); + +% Create a simple hollow cube phantom +cube = zeros(128,128,128); +cube(17:112,17:112,17:112) = 1; +cube(33:96,33:96,33:96) = 0; + +% Create projection data from this +[proj_id, proj_data] = astra_create_sino3d_cuda(cube, proj_geom, vol_geom); + +% Display a single projection image +figure, imshow(squeeze(proj_data(:,20,:))',[]) + +% Create a data object for the reconstruction +rec_id = astra_mex_data3d('create', '-vol', vol_geom); + +% Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra_struct('SIRT3D_CUDA'); +cfg.ReconstructionDataId = rec_id; +cfg.ProjectionDataId = proj_id; + + +% Create the algorithm object from the configuration structure +alg_id = astra_mex_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_mex_algorithm('iterate', alg_id, 150); + +% Get the result +rec = astra_mex_data3d('get', rec_id); +figure, imshow(squeeze(rec(:,:,65)),[]); + + +% Clean up. Note that GPU memory is tied up in the algorithm object, +% and main RAM in the data objects. +astra_mex_algorithm('delete', alg_id); +astra_mex_data3d('delete', rec_id); +astra_mex_data3d('delete', proj_id); diff --git a/samples/matlab/s008_gpu_selection.m b/samples/matlab/s008_gpu_selection.m new file mode 100644 index 0000000..a9e152d --- /dev/null +++ b/samples/matlab/s008_gpu_selection.m @@ -0,0 +1,37 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); +P = phantom(256); + +% Create a sinogram from a phantom, using GPU #1. (The default is #0) +[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom, 1); + + +% Set up the parameters for a reconstruction algorithm using the GPU +rec_id = astra_mex_data2d('create', '-vol', vol_geom); +cfg = astra_struct('SIRT_CUDA'); +cfg.ReconstructionDataId = rec_id; +cfg.ProjectionDataId = sinogram_id; + +% Use GPU #1 for the reconstruction. (The default is #0.) +cfg.option.GPUindex = 1; + +% Run 150 iterations of the algorithm +alg_id = astra_mex_algorithm('create', cfg); +astra_mex_algorithm('iterate', alg_id, 150); +rec = astra_mex_data2d('get', rec_id); + + +% Clean up. +astra_mex_algorithm('delete', alg_id); +astra_mex_data2d('delete', rec_id); +astra_mex_data2d('delete', sinogram_id); diff --git a/samples/matlab/s009_projection_matrix.m b/samples/matlab/s009_projection_matrix.m new file mode 100644 index 0000000..efda0d2 --- /dev/null +++ b/samples/matlab/s009_projection_matrix.m @@ -0,0 +1,45 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); + +% 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); + +% 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_mex_projector('matrix', proj_id); + +% Get the projection matrix as a Matlab sparse matrix. +W = astra_mex_matrix('get', matrix_id); + + +% Manually use this projection matrix to do a projection: +P = phantom(256)'; +s = W * P(:); +s = reshape(s, [proj_geom.DetectorCount size(proj_geom.ProjectionAngles, 2)])'; +figure(1), imshow(s,[]); + +% Because Matlab's matrices are stored transposed in memory compared to C++, +% reshaping them to a vector doesn't give the right ordering for multiplication +% with W. We have to take the transpose of the input and output to get the same +% results (up to numerical noise) as using the toolbox directly. + +% 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_mex_projector('delete', proj_id); +astra_mex_matrix('delete', matrix_id); diff --git a/samples/matlab/s010_supersampling.m b/samples/matlab/s010_supersampling.m new file mode 100644 index 0000000..80f6f56 --- /dev/null +++ b/samples/matlab/s010_supersampling.m @@ -0,0 +1,58 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 3.0, 128, linspace2(0,pi,180)); +P = phantom(256); + +% Because the astra_create_sino_gpu wrapper does not have support for +% all possible algorithm options, we manually create a sinogram +phantom_id = astra_mex_data2d('create', '-vol', vol_geom, P); +sinogram_id = astra_mex_data2d('create', '-sino', proj_geom); +cfg = astra_struct('FP_CUDA'); +cfg.VolumeDataId = phantom_id; +cfg.ProjectionDataId = sinogram_id; + +% Set up 3 rays per detector element +cfg.option.DetectorSuperSampling = 3; + +alg_id = astra_mex_algorithm('create', cfg); +astra_mex_algorithm('run', alg_id); +astra_mex_algorithm('delete', alg_id); +astra_mex_data2d('delete', phantom_id); + +sinogram3 = astra_mex_data2d('get', sinogram_id); + +figure(1); imshow(P, []); +figure(2); imshow(sinogram3, []); + + +% Create a reconstruction, also using supersampling +rec_id = astra_mex_data2d('create', '-vol', vol_geom); +cfg = astra_struct('SIRT_CUDA'); +cfg.ReconstructionDataId = rec_id; +cfg.ProjectionDataId = sinogram_id; +% Set up 3 rays per detector element +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_mex_algorithm('create', cfg); +astra_mex_algorithm('iterate', alg_id, 150); +astra_mex_algorithm('delete', alg_id); + +rec = astra_mex_data2d('get', rec_id); +figure(3); imshow(rec, []); + diff --git a/samples/matlab/s011_object_info.m b/samples/matlab/s011_object_info.m new file mode 100644 index 0000000..61ecb83 --- /dev/null +++ b/samples/matlab/s011_object_info.m @@ -0,0 +1,36 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +% Create two volume geometries +vol_geom1 = astra_create_vol_geom(256, 256); +vol_geom2 = astra_create_vol_geom(512, 256); + +% Create volumes +v0 = astra_mex_data2d('create', '-vol', vol_geom1); +v1 = astra_mex_data2d('create', '-vol', vol_geom2); +v2 = astra_mex_data2d('create', '-vol', vol_geom2); + +% Show the currently allocated volumes +astra_mex_data2d('info'); + + +astra_mex_data2d('delete', v2); +astra_mex_data2d('info'); + +astra_mex_data2d('clear'); +astra_mex_data2d('info'); + + + +% The same clear and info command also work for other object types: +astra_mex_algorithm('info'); +astra_mex_data3d('info'); +astra_mex_projector('info'); +astra_mex_matrix('info'); diff --git a/samples/matlab/s012_masks.m b/samples/matlab/s012_masks.m new file mode 100644 index 0000000..d3611a6 --- /dev/null +++ b/samples/matlab/s012_masks.m @@ -0,0 +1,60 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + + +% 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 = -127.5:127.5; +[x y] = meshgrid(-127.5:127.5,-127.5:127.5); +mask = (x.^2 + y.^2 < 127.5^2); + +figure(1); imshow(mask, []); + + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,50)); + +% As before, create a sinogram from a phantom +P = phantom(256); +[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); +figure(2); imshow(P, []); +figure(3); imshow(sinogram, []); + +% Create a data object for the reconstruction +rec_id = astra_mex_data2d('create', '-vol', vol_geom); + +% Create a data object for the mask +mask_id = astra_mex_data2d('create', '-vol', vol_geom, mask); + +% Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra_struct('SIRT_CUDA'); +cfg.ReconstructionDataId = rec_id; +cfg.ProjectionDataId = sinogram_id; +cfg.option.ReconstructionMaskId = mask_id; + +% Create the algorithm object from the configuration structure +alg_id = astra_mex_algorithm('create', cfg); + +% Run 150 iterations of the algorithm +astra_mex_algorithm('iterate', alg_id, 150); + +% Get the result +rec = astra_mex_data2d('get', rec_id); +figure(4); imshow(rec, []); + +% Clean up. Note that GPU memory is tied up in the algorithm object, +% and main RAM in the data objects. +astra_mex_algorithm('delete', alg_id); +astra_mex_data2d('delete', mask_id); +astra_mex_data2d('delete', rec_id); +astra_mex_data2d('delete', sinogram_id); diff --git a/samples/matlab/s013_constraints.m b/samples/matlab/s013_constraints.m new file mode 100644 index 0000000..d72195c --- /dev/null +++ b/samples/matlab/s013_constraints.m @@ -0,0 +1,47 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +% 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, linspace2(0,pi,50)); + +% As before, create a sinogram from a phantom +P = phantom(256); +[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); +figure(1); imshow(P, []); +figure(2); imshow(sinogram, []); + +% Create a data object for the reconstruction +rec_id = astra_mex_data2d('create', '-vol', vol_geom); + +% Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra_struct('SIRT_CUDA'); +cfg.ReconstructionDataId = rec_id; +cfg.ProjectionDataId = sinogram_id; +cfg.option.MinConstraint = 0; +cfg.option.MaxConstraint = 1; + +% Create the algorithm object from the configuration structure +alg_id = astra_mex_algorithm('create', cfg); + +% Run 150 iterations of the algorithm +astra_mex_algorithm('iterate', alg_id, 150); + +% Get the result +rec = astra_mex_data2d('get', rec_id); +figure(3); imshow(rec, []); + +% Clean up. Note that GPU memory is tied up in the algorithm object, +% and main RAM in the data objects. +astra_mex_algorithm('delete', alg_id); +astra_mex_data2d('delete', rec_id); +astra_mex_data2d('delete', sinogram_id); diff --git a/samples/matlab/s014_FBP.m b/samples/matlab/s014_FBP.m new file mode 100644 index 0000000..b73149c --- /dev/null +++ b/samples/matlab/s014_FBP.m @@ -0,0 +1,47 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); + +% As before, create a sinogram from a phantom +P = phantom(256); +[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); +figure(1); imshow(P, []); +figure(2); imshow(sinogram, []); + +% Create a data object for the reconstruction +rec_id = astra_mex_data2d('create', '-vol', vol_geom); + +% create configuration +cfg = astra_struct('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_mex_algorithm('create', cfg); +astra_mex_algorithm('run', alg_id); + +% Get the result +rec = astra_mex_data2d('get', rec_id); +figure(3); imshow(rec, []); + +% Clean up. Note that GPU memory is tied up in the algorithm object, +% and main RAM in the data objects. +astra_mex_algorithm('delete', alg_id); +astra_mex_data2d('delete', rec_id); +astra_mex_data2d('delete', sinogram_id); diff --git a/samples/matlab/s015_fp_bp.m b/samples/matlab/s015_fp_bp.m new file mode 100644 index 0000000..8cc417e --- /dev/null +++ b/samples/matlab/s015_fp_bp.m @@ -0,0 +1,62 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + + +% 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. + +function s015_fp_bp + + +% FP/BP wrapper function +function Y = astra_wrap(X,T) + if strcmp(T, 'notransp') + % X is passed as a vector. Reshape it into an image. + [sid, s] = astra_create_sino_cuda(reshape(X,[vol_geom.GridRowCount vol_geom.GridColCount])', proj_geom, vol_geom); + astra_mex_data2d('delete', sid); + % now s is the sinogram. Reshape it back into a vector + Y = reshape(s',[prod(size(s)) 1]); + else + % X is passed as a vector. Reshape it into a sinogram. + v = astra_create_backprojection_cuda(reshape(X, [proj_geom.DetectorCount size(proj_geom.ProjectionAngles,2)])', proj_geom, vol_geom); + % now v is the resulting volume. Reshape it back into a vector + Y = reshape(v',[prod(size(v)) 1]); + end +end + + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); + +% Create a 256x256 phantom image using matlab's built-in phantom() function +P = phantom(256); + +% Create a sinogram using the GPU. +[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); + +% Reshape the sinogram into a vector +b = reshape(sinogram',[prod(size(sinogram)) 1]); + +% Call Matlab's lsqr with ASTRA FP and BP +Y = lsqr(@astra_wrap,b,1e-4,25); + +% Reshape the result into an image +Y = reshape(Y,[vol_geom.GridRowCount vol_geom.GridColCount])'; +imshow(Y,[]); + + +astra_mex_data2d('delete', sinogram_id); + +end diff --git a/samples/matlab/s016_plots.m b/samples/matlab/s016_plots.m new file mode 100644 index 0000000..1455c6d --- /dev/null +++ b/samples/matlab/s016_plots.m @@ -0,0 +1,54 @@ +% ----------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp +% 2014-2015, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@uantwerpen.be +% Website: http://sf.net/projects/astra-toolbox +% ----------------------------------------------------------------------- + +vol_geom = astra_create_vol_geom(256, 256); +proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); + +% As before, create a sinogram from a phantom +P = phantom(256); +[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); +figure(1); imshow(P, []); +figure(2); imshow(sinogram, []); + +% Create a data object for the reconstruction +rec_id = astra_mex_data2d('create', '-vol', vol_geom); + +% Set up the parameters for a reconstruction algorithm using the GPU +cfg = astra_struct('SIRT_CUDA'); +cfg.ReconstructionDataId = rec_id; +cfg.ProjectionDataId = sinogram_id; + +% Create the algorithm object from the configuration structure +alg_id = astra_mex_algorithm('create', cfg); + +% Run 1500 iterations of the algorithm one at a time, keeping track of errors +nIters = 1500; +phantom_error = zeros(1, nIters); +residual_error = zeros(1, nIters); +for i = 1:nIters; + % Run a single iteration + astra_mex_algorithm('iterate', alg_id, 1); + + residual_error(i) = astra_mex_algorithm('get_res_norm', alg_id); + rec = astra_mex_data2d('get', rec_id); + phantom_error(i) = sqrt(sumsqr(rec - P)); +end + +% Get the result +rec = astra_mex_data2d('get', rec_id); +figure(3); imshow(rec, []); + +figure(4); plot(residual_error) +figure(5); plot(phantom_error) + +% Clean up. +astra_mex_algorithm('delete', alg_id); +astra_mex_data2d('delete', rec_id); +astra_mex_data2d('delete', sinogram_id); 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) diff --git a/samples/s001_sinogram_par2d.m b/samples/s001_sinogram_par2d.m deleted file mode 100644 index 4494e7b..0000000 --- a/samples/s001_sinogram_par2d.m +++ /dev/null @@ -1,33 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -% 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, linspace2(0,pi,180)); - -% Create a 256x256 phantom image using matlab's built-in phantom() function -P = phantom(256); - -% 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. -[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); - -figure(1); imshow(P, []); -figure(2); imshow(sinogram, []); - - -% Free memory -astra_mex_data2d('delete', sinogram_id); diff --git a/samples/s002_data2d.m b/samples/s002_data2d.m deleted file mode 100644 index a91071f..0000000 --- a/samples/s002_data2d.m +++ /dev/null @@ -1,60 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(256, 256); - -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); - - -% Create volumes - -% initialized to zero -v0 = astra_mex_data2d('create', '-vol', vol_geom); - -% initialized to 3.0 -v1 = astra_mex_data2d('create', '-vol', vol_geom, 3.0); - -% initialized to a matrix. A may be a single, double or logical (0/1) array. -A = phantom(256); -v2 = astra_mex_data2d('create', '-vol', vol_geom, A); - - -% Projection data -s0 = astra_mex_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_mex_data2d('store', v0, 0); - -% set to a matrix -astra_mex_data2d('store', v2, A); - - - -% Retrieve data - -R = astra_mex_data2d('get', v2); -imshow(R, []); - - -% Retrieve data as a single array. Since astra internally stores -% data as single precision, this is more efficient: - -R = astra_mex_data2d('get_single', v2); - - -% Free memory -astra_mex_data2d('delete', v0); -astra_mex_data2d('delete', v1); -astra_mex_data2d('delete', v2); -astra_mex_data2d('delete', s0); diff --git a/samples/s003_gpu_reconstruction.m b/samples/s003_gpu_reconstruction.m deleted file mode 100644 index efb5c68..0000000 --- a/samples/s003_gpu_reconstruction.m +++ /dev/null @@ -1,52 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); - -% As before, create a sinogram from a phantom -P = phantom(256); -[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); -figure(1); imshow(P, []); -figure(2); imshow(sinogram, []); - -astra_mex_data2d('delete', sinogram_id); - -% We now re-create the sinogram data object as we would do when loading -% an external sinogram -sinogram_id = astra_mex_data2d('create', '-sino', proj_geom, sinogram); - -% Create a data object for the reconstruction -rec_id = astra_mex_data2d('create', '-vol', vol_geom); - -% Set up the parameters for a reconstruction algorithm using the GPU -cfg = astra_struct('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_mex_algorithm('create', cfg); - -% Run 150 iterations of the algorithm -astra_mex_algorithm('iterate', alg_id, 150); - -% Get the result -rec = astra_mex_data2d('get', rec_id); -figure(3); imshow(rec, []); - -% Clean up. Note that GPU memory is tied up in the algorithm object, -% and main RAM in the data objects. -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', rec_id); -astra_mex_data2d('delete', sinogram_id); diff --git a/samples/s004_cpu_reconstruction.m b/samples/s004_cpu_reconstruction.m deleted file mode 100644 index f25cd2b..0000000 --- a/samples/s004_cpu_reconstruction.m +++ /dev/null @@ -1,60 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); - -% 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 -P = phantom(256); -[sinogram_id, sinogram] = astra_create_sino(P, proj_id); -figure(1); imshow(P, []); -figure(2); imshow(sinogram, []); - -astra_mex_data2d('delete', sinogram_id); - -% We now re-create the sinogram data object as we would do when loading -% an external sinogram -sinogram_id = astra_mex_data2d('create', '-sino', proj_geom, sinogram); - -% Create a data object for the reconstruction -rec_id = astra_mex_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_struct('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_mex_algorithm('create', cfg); - -% Run 20 iterations of the algorithm -% This will have a runtime in the order of 10 seconds. -astra_mex_algorithm('iterate', alg_id, 20); - -% Get the result -rec = astra_mex_data2d('get', rec_id); -figure(3); imshow(rec, []); - -% Clean up. -astra_mex_projector('delete', proj_id); -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', rec_id); -astra_mex_data2d('delete', sinogram_id); diff --git a/samples/s005_3d_geometry.m b/samples/s005_3d_geometry.m deleted file mode 100644 index 4b7360b..0000000 --- a/samples/s005_3d_geometry.m +++ /dev/null @@ -1,98 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -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 = linspace2(0, 2*pi, 48); -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 = zeros(numel(angles), 12); -for i = 1:numel(angles) - % ray direction - vectors(i,1) = sin(angles(i)); - vectors(i,2) = -cos(angles(i)); - vectors(i,3) = 0; - - % center of detector - vectors(i,4:6) = 0; - - % vector from detector pixel (0,0) to (0,1) - vectors(i,7) = cos(angles(i)); - vectors(i,8) = sin(angles(i)); - vectors(i,9) = 0; - - % vector from detector pixel (0,0) to (1,0) - vectors(i,10) = 0; - vectors(i,11) = 0; - vectors(i,12) = 1; -end - -% 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 = linspace2(0, 2*pi, 48); -proj_geom = astra_create_proj_geom('cone', 1.0, 1.0, 32, 64, ... - angles, 1000, 0); - -% Free - -vectors = zeros(numel(angles), 12); -for i = 1:numel(angles) - - % source - vectors(i,1) = sin(angles(i)) * 1000; - vectors(i,2) = -cos(angles(i)) * 1000; - vectors(i,3) = 0; - - % center of detector - vectors(i,4:6) = 0; - - % vector from detector pixel (0,0) to (0,1) - vectors(i,7) = cos(angles(i)); - vectors(i,8) = sin(angles(i)); - vectors(i,9) = 0; - - % vector from detector pixel (0,0) to (1,0) - vectors(i,10) = 0; - vectors(i,11) = 0; - vectors(i,12) = 1; -end - -% Parameters: #rows, #columns, vectors -proj_geom = astra_create_proj_geom('cone_vec', 32, 64, vectors); - diff --git a/samples/s006_3d_data.m b/samples/s006_3d_data.m deleted file mode 100644 index 32d88cc..0000000 --- a/samples/s006_3d_data.m +++ /dev/null @@ -1,62 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -% 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_mex_data3d('create', '-vol', vol_geom); - -% initialized to 3.0 -v1 = astra_mex_data3d('create', '-vol', vol_geom, 3.0); - -% initialized to a matrix. A may be a single or double array. -% Coordinate order: column, row, slice (x, y, z) -A = zeros(48, 64, 32); -v2 = astra_mex_data3d('create', '-vol', vol_geom, A); - - -% Projection data - -% 2 projection directions, along x and y axis resp. -V = [ 1 0 0 0 0 0 0 1 0 0 0 1 ; ... - 0 1 0 0 0 0 -1 0 0 0 0 1 ]; -% 32 rows (v), 64 columns (u) -proj_geom = astra_create_proj_geom('parallel3d_vec', 32, 64, V); - -s0 = astra_mex_data3d('create', '-proj3d', proj_geom); - -% Initialization to a scalar or zero works exactly as with a volume. - -% Initialized to a matrix: -% Coordinate order: column (u), angle, row (v) -A = zeros(64, 2, 32); -s1 = astra_mex_data3d('create', '-proj3d', proj_geom, A); - - -% Retrieve data: -R = astra_mex_data3d('get', v1); - -% Retrieve data as a single array. Since astra internally stores -% data as single precision, this is more efficient: -R = astra_mex_data3d('get_single', v1); - - - -% Delete all created data objects -astra_mex_data3d('delete', v0); -astra_mex_data3d('delete', v1); -astra_mex_data3d('delete', v2); -astra_mex_data3d('delete', s0); -astra_mex_data3d('delete', s1); diff --git a/samples/s007_3d_reconstruction.m b/samples/s007_3d_reconstruction.m deleted file mode 100644 index fc9aca6..0000000 --- a/samples/s007_3d_reconstruction.m +++ /dev/null @@ -1,53 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(128, 128, 128); - -angles = linspace2(0, pi, 180); -proj_geom = astra_create_proj_geom('parallel3d', 1.0, 1.0, 128, 192, angles); - -% Create a simple hollow cube phantom -cube = zeros(128,128,128); -cube(17:112,17:112,17:112) = 1; -cube(33:96,33:96,33:96) = 0; - -% Create projection data from this -[proj_id, proj_data] = astra_create_sino3d_cuda(cube, proj_geom, vol_geom); - -% Display a single projection image -figure, imshow(squeeze(proj_data(:,20,:))',[]) - -% Create a data object for the reconstruction -rec_id = astra_mex_data3d('create', '-vol', vol_geom); - -% Set up the parameters for a reconstruction algorithm using the GPU -cfg = astra_struct('SIRT3D_CUDA'); -cfg.ReconstructionDataId = rec_id; -cfg.ProjectionDataId = proj_id; - - -% Create the algorithm object from the configuration structure -alg_id = astra_mex_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_mex_algorithm('iterate', alg_id, 150); - -% Get the result -rec = astra_mex_data3d('get', rec_id); -figure, imshow(squeeze(rec(:,:,65)),[]); - - -% Clean up. Note that GPU memory is tied up in the algorithm object, -% and main RAM in the data objects. -astra_mex_algorithm('delete', alg_id); -astra_mex_data3d('delete', rec_id); -astra_mex_data3d('delete', proj_id); diff --git a/samples/s008_gpu_selection.m b/samples/s008_gpu_selection.m deleted file mode 100644 index a9e152d..0000000 --- a/samples/s008_gpu_selection.m +++ /dev/null @@ -1,37 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); -P = phantom(256); - -% Create a sinogram from a phantom, using GPU #1. (The default is #0) -[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom, 1); - - -% Set up the parameters for a reconstruction algorithm using the GPU -rec_id = astra_mex_data2d('create', '-vol', vol_geom); -cfg = astra_struct('SIRT_CUDA'); -cfg.ReconstructionDataId = rec_id; -cfg.ProjectionDataId = sinogram_id; - -% Use GPU #1 for the reconstruction. (The default is #0.) -cfg.option.GPUindex = 1; - -% Run 150 iterations of the algorithm -alg_id = astra_mex_algorithm('create', cfg); -astra_mex_algorithm('iterate', alg_id, 150); -rec = astra_mex_data2d('get', rec_id); - - -% Clean up. -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', rec_id); -astra_mex_data2d('delete', sinogram_id); diff --git a/samples/s009_projection_matrix.m b/samples/s009_projection_matrix.m deleted file mode 100644 index efda0d2..0000000 --- a/samples/s009_projection_matrix.m +++ /dev/null @@ -1,45 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); - -% 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); - -% 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_mex_projector('matrix', proj_id); - -% Get the projection matrix as a Matlab sparse matrix. -W = astra_mex_matrix('get', matrix_id); - - -% Manually use this projection matrix to do a projection: -P = phantom(256)'; -s = W * P(:); -s = reshape(s, [proj_geom.DetectorCount size(proj_geom.ProjectionAngles, 2)])'; -figure(1), imshow(s,[]); - -% Because Matlab's matrices are stored transposed in memory compared to C++, -% reshaping them to a vector doesn't give the right ordering for multiplication -% with W. We have to take the transpose of the input and output to get the same -% results (up to numerical noise) as using the toolbox directly. - -% 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_mex_projector('delete', proj_id); -astra_mex_matrix('delete', matrix_id); diff --git a/samples/s010_supersampling.m b/samples/s010_supersampling.m deleted file mode 100644 index 80f6f56..0000000 --- a/samples/s010_supersampling.m +++ /dev/null @@ -1,58 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 3.0, 128, linspace2(0,pi,180)); -P = phantom(256); - -% Because the astra_create_sino_gpu wrapper does not have support for -% all possible algorithm options, we manually create a sinogram -phantom_id = astra_mex_data2d('create', '-vol', vol_geom, P); -sinogram_id = astra_mex_data2d('create', '-sino', proj_geom); -cfg = astra_struct('FP_CUDA'); -cfg.VolumeDataId = phantom_id; -cfg.ProjectionDataId = sinogram_id; - -% Set up 3 rays per detector element -cfg.option.DetectorSuperSampling = 3; - -alg_id = astra_mex_algorithm('create', cfg); -astra_mex_algorithm('run', alg_id); -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', phantom_id); - -sinogram3 = astra_mex_data2d('get', sinogram_id); - -figure(1); imshow(P, []); -figure(2); imshow(sinogram3, []); - - -% Create a reconstruction, also using supersampling -rec_id = astra_mex_data2d('create', '-vol', vol_geom); -cfg = astra_struct('SIRT_CUDA'); -cfg.ReconstructionDataId = rec_id; -cfg.ProjectionDataId = sinogram_id; -% Set up 3 rays per detector element -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_mex_algorithm('create', cfg); -astra_mex_algorithm('iterate', alg_id, 150); -astra_mex_algorithm('delete', alg_id); - -rec = astra_mex_data2d('get', rec_id); -figure(3); imshow(rec, []); - diff --git a/samples/s011_object_info.m b/samples/s011_object_info.m deleted file mode 100644 index 61ecb83..0000000 --- a/samples/s011_object_info.m +++ /dev/null @@ -1,36 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -% Create two volume geometries -vol_geom1 = astra_create_vol_geom(256, 256); -vol_geom2 = astra_create_vol_geom(512, 256); - -% Create volumes -v0 = astra_mex_data2d('create', '-vol', vol_geom1); -v1 = astra_mex_data2d('create', '-vol', vol_geom2); -v2 = astra_mex_data2d('create', '-vol', vol_geom2); - -% Show the currently allocated volumes -astra_mex_data2d('info'); - - -astra_mex_data2d('delete', v2); -astra_mex_data2d('info'); - -astra_mex_data2d('clear'); -astra_mex_data2d('info'); - - - -% The same clear and info command also work for other object types: -astra_mex_algorithm('info'); -astra_mex_data3d('info'); -astra_mex_projector('info'); -astra_mex_matrix('info'); diff --git a/samples/s012_masks.m b/samples/s012_masks.m deleted file mode 100644 index d3611a6..0000000 --- a/samples/s012_masks.m +++ /dev/null @@ -1,60 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - - -% 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 = -127.5:127.5; -[x y] = meshgrid(-127.5:127.5,-127.5:127.5); -mask = (x.^2 + y.^2 < 127.5^2); - -figure(1); imshow(mask, []); - - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,50)); - -% As before, create a sinogram from a phantom -P = phantom(256); -[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); -figure(2); imshow(P, []); -figure(3); imshow(sinogram, []); - -% Create a data object for the reconstruction -rec_id = astra_mex_data2d('create', '-vol', vol_geom); - -% Create a data object for the mask -mask_id = astra_mex_data2d('create', '-vol', vol_geom, mask); - -% Set up the parameters for a reconstruction algorithm using the GPU -cfg = astra_struct('SIRT_CUDA'); -cfg.ReconstructionDataId = rec_id; -cfg.ProjectionDataId = sinogram_id; -cfg.option.ReconstructionMaskId = mask_id; - -% Create the algorithm object from the configuration structure -alg_id = astra_mex_algorithm('create', cfg); - -% Run 150 iterations of the algorithm -astra_mex_algorithm('iterate', alg_id, 150); - -% Get the result -rec = astra_mex_data2d('get', rec_id); -figure(4); imshow(rec, []); - -% Clean up. Note that GPU memory is tied up in the algorithm object, -% and main RAM in the data objects. -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', mask_id); -astra_mex_data2d('delete', rec_id); -astra_mex_data2d('delete', sinogram_id); diff --git a/samples/s013_constraints.m b/samples/s013_constraints.m deleted file mode 100644 index d72195c..0000000 --- a/samples/s013_constraints.m +++ /dev/null @@ -1,47 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -% 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, linspace2(0,pi,50)); - -% As before, create a sinogram from a phantom -P = phantom(256); -[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); -figure(1); imshow(P, []); -figure(2); imshow(sinogram, []); - -% Create a data object for the reconstruction -rec_id = astra_mex_data2d('create', '-vol', vol_geom); - -% Set up the parameters for a reconstruction algorithm using the GPU -cfg = astra_struct('SIRT_CUDA'); -cfg.ReconstructionDataId = rec_id; -cfg.ProjectionDataId = sinogram_id; -cfg.option.MinConstraint = 0; -cfg.option.MaxConstraint = 1; - -% Create the algorithm object from the configuration structure -alg_id = astra_mex_algorithm('create', cfg); - -% Run 150 iterations of the algorithm -astra_mex_algorithm('iterate', alg_id, 150); - -% Get the result -rec = astra_mex_data2d('get', rec_id); -figure(3); imshow(rec, []); - -% Clean up. Note that GPU memory is tied up in the algorithm object, -% and main RAM in the data objects. -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', rec_id); -astra_mex_data2d('delete', sinogram_id); diff --git a/samples/s014_FBP.m b/samples/s014_FBP.m deleted file mode 100644 index b73149c..0000000 --- a/samples/s014_FBP.m +++ /dev/null @@ -1,47 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); - -% As before, create a sinogram from a phantom -P = phantom(256); -[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); -figure(1); imshow(P, []); -figure(2); imshow(sinogram, []); - -% Create a data object for the reconstruction -rec_id = astra_mex_data2d('create', '-vol', vol_geom); - -% create configuration -cfg = astra_struct('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_mex_algorithm('create', cfg); -astra_mex_algorithm('run', alg_id); - -% Get the result -rec = astra_mex_data2d('get', rec_id); -figure(3); imshow(rec, []); - -% Clean up. Note that GPU memory is tied up in the algorithm object, -% and main RAM in the data objects. -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', rec_id); -astra_mex_data2d('delete', sinogram_id); diff --git a/samples/s015_fp_bp.m b/samples/s015_fp_bp.m deleted file mode 100644 index 8cc417e..0000000 --- a/samples/s015_fp_bp.m +++ /dev/null @@ -1,62 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - - -% 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. - -function s015_fp_bp - - -% FP/BP wrapper function -function Y = astra_wrap(X,T) - if strcmp(T, 'notransp') - % X is passed as a vector. Reshape it into an image. - [sid, s] = astra_create_sino_cuda(reshape(X,[vol_geom.GridRowCount vol_geom.GridColCount])', proj_geom, vol_geom); - astra_mex_data2d('delete', sid); - % now s is the sinogram. Reshape it back into a vector - Y = reshape(s',[prod(size(s)) 1]); - else - % X is passed as a vector. Reshape it into a sinogram. - v = astra_create_backprojection_cuda(reshape(X, [proj_geom.DetectorCount size(proj_geom.ProjectionAngles,2)])', proj_geom, vol_geom); - % now v is the resulting volume. Reshape it back into a vector - Y = reshape(v',[prod(size(v)) 1]); - end -end - - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); - -% Create a 256x256 phantom image using matlab's built-in phantom() function -P = phantom(256); - -% Create a sinogram using the GPU. -[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); - -% Reshape the sinogram into a vector -b = reshape(sinogram',[prod(size(sinogram)) 1]); - -% Call Matlab's lsqr with ASTRA FP and BP -Y = lsqr(@astra_wrap,b,1e-4,25); - -% Reshape the result into an image -Y = reshape(Y,[vol_geom.GridRowCount vol_geom.GridColCount])'; -imshow(Y,[]); - - -astra_mex_data2d('delete', sinogram_id); - -end diff --git a/samples/s016_plots.m b/samples/s016_plots.m deleted file mode 100644 index 1455c6d..0000000 --- a/samples/s016_plots.m +++ /dev/null @@ -1,54 +0,0 @@ -% ----------------------------------------------------------------------- -% This file is part of the ASTRA Toolbox -% -% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -% 2014-2015, CWI, Amsterdam -% License: Open Source under GPLv3 -% Contact: astra@uantwerpen.be -% Website: http://sf.net/projects/astra-toolbox -% ----------------------------------------------------------------------- - -vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); - -% As before, create a sinogram from a phantom -P = phantom(256); -[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom, vol_geom); -figure(1); imshow(P, []); -figure(2); imshow(sinogram, []); - -% Create a data object for the reconstruction -rec_id = astra_mex_data2d('create', '-vol', vol_geom); - -% Set up the parameters for a reconstruction algorithm using the GPU -cfg = astra_struct('SIRT_CUDA'); -cfg.ReconstructionDataId = rec_id; -cfg.ProjectionDataId = sinogram_id; - -% Create the algorithm object from the configuration structure -alg_id = astra_mex_algorithm('create', cfg); - -% Run 1500 iterations of the algorithm one at a time, keeping track of errors -nIters = 1500; -phantom_error = zeros(1, nIters); -residual_error = zeros(1, nIters); -for i = 1:nIters; - % Run a single iteration - astra_mex_algorithm('iterate', alg_id, 1); - - residual_error(i) = astra_mex_algorithm('get_res_norm', alg_id); - rec = astra_mex_data2d('get', rec_id); - phantom_error(i) = sqrt(sumsqr(rec - P)); -end - -% Get the result -rec = astra_mex_data2d('get', rec_id); -figure(3); imshow(rec, []); - -figure(4); plot(residual_error) -figure(5); plot(phantom_error) - -% Clean up. -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', rec_id); -astra_mex_data2d('delete', sinogram_id); -- cgit v1.2.3