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 ++++++++++++++++++ 16 files changed, 864 insertions(+) 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 (limited to 'samples/matlab') 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); -- cgit v1.2.3