From dc31a371a7f5c66e32226698999f345e350f3049 Mon Sep 17 00:00:00 2001
From: Folkert Bleichrodt <F.Bleichrodt@cwi.nl>
Date: Wed, 8 Apr 2015 16:03:41 +0200
Subject: Added ASTRA-Spot operator opTomo
A Spot operator for the ASTRA projectors. Wraps the forward and
backprojection operation into a Spot operator, which can be used
with matrix-vector syntax in Matlab.
---
matlab/tools/opTomo.m | 280 ++++++++++++++++++++++++++++++++++++
matlab/tools/opTomo_helper_handle.m | 29 ++++
samples/matlab/s017_opTomo.m | 44 ++++++
3 files changed, 353 insertions(+)
create mode 100644 matlab/tools/opTomo.m
create mode 100644 matlab/tools/opTomo_helper_handle.m
create mode 100644 samples/matlab/s017_opTomo.m
diff --git a/matlab/tools/opTomo.m b/matlab/tools/opTomo.m
new file mode 100644
index 0000000..14128d2
--- /dev/null
+++ b/matlab/tools/opTomo.m
@@ -0,0 +1,280 @@
+%OPTOMO Wrapper for ASTRA tomography projector
+%
+% OP = OPTOMO(TYPE, PROJ_GEOM, VOL_GEOM) generates a Spot operator OP for
+% the ASTRA forward and backprojection operations. The string TYPE
+% determines the model used for the projections. Possible choices are:
+% TYPE: * using the CPU
+% 'line' - use a line kernel
+% 'linear' - use a Joseph kernel
+% 'strip' - use the strip kernel
+% * using the GPU
+% 'cuda' - use a Joseph kernel, on the GPU, currently using
+% 'cuda' is the only option in 3D.
+% The PROJ_GEOM and VOL_GEOM structures are projection and volume
+% geometries as used in the ASTRA toolbox.
+%
+% OP = OPTOMO(TYPE, PROJ_GEOM, VOL_GEOM, GPU_INDEX) also specify the
+% index of the GPU that should be used, if multiple GPUs are present in
+% the host system. By default GPU_INDEX is 0.
+%
+% Note: this code depends on the Matlab toolbox
+% "Spot - A Linear-Operator Toolbox" which can be downloaded from
+% http://www.cs.ubc.ca/labs/scl/spot/
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2014-2015, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Author: Folkert Bleichrodt
+% Contact: F.Bleichrodt@cwi.nl
+% Website: http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+% $Id$
+
+classdef opTomo < opSpot
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Properties
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ properties ( Access = private )
+ % multiplication function
+ funHandle
+ % ASTRA identifiers
+ sino_id
+ vol_id
+ fp_alg_id
+ bp_alg_id
+ % ASTRA IDs handle
+ astra_handle
+ % geometries
+ proj_geom;
+ vol_geom;
+ end % properties
+
+ properties ( SetAccess = private, GetAccess = public )
+ proj_size
+ vol_size
+ end % properties
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Methods - public
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ methods
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % opTomo - constructor
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ function op = opTomo(type, proj_geom, vol_geom, gpu_index)
+
+ if nargin < 4 || isempty(gpu_index), gpu_index = 0; end
+
+ proj_size = astra_geom_size(proj_geom);
+ vol_size = astra_geom_size(vol_geom);
+
+ % construct operator
+ op = op@opSpot('opTomo', prod(proj_size), prod(vol_size));
+
+ % determine the dimension
+ is2D = ~isfield(vol_geom, 'GridSliceCount');
+ gpuEnabled = strcmpi(type, 'cuda');
+
+ if is2D
+ % create a projector
+ proj_id = astra_create_projector(type, proj_geom, vol_geom);
+
+ % create a function handle
+ op.funHandle = @opTomo_intrnl2D;
+
+ % Initialize ASTRA data objects.
+ % projection data
+ sino_id = astra_mex_data2d('create', '-sino', proj_geom, 0);
+ % image data
+ vol_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+ % Setup forward and back projection algorithms.
+ if gpuEnabled
+ fp_alg = 'FP_CUDA';
+ bp_alg = 'BP_CUDA';
+ proj_id = [];
+ else
+ fp_alg = 'FP';
+ bp_alg = 'BP';
+ proj_id = astra_create_projector(type, proj_geom, vol_geom);
+ end
+
+ % configuration for ASTRA fp algorithm
+ cfg_fp = astra_struct(fp_alg);
+ cfg_fp.ProjectorId = proj_id;
+ cfg_fp.ProjectionDataId = sino_id;
+ cfg_fp.VolumeDataId = vol_id;
+
+ % configuration for ASTRA bp algorithm
+ cfg_bp = astra_struct(bp_alg);
+ cfg_bp.ProjectionDataId = sino_id;
+ cfg_bp.ProjectorId = proj_id;
+ cfg_bp.ReconstructionDataId = vol_id;
+
+ % set GPU index
+ if gpuEnabled
+ cfg_fp.option.GPUindex = gpu_index;
+ cfg_bp.option.GPUindex = gpu_index;
+ end
+
+ fp_alg_id = astra_mex_algorithm('create', cfg_fp);
+ bp_alg_id = astra_mex_algorithm('create', cfg_bp);
+
+ % Create handle to ASTRA objects, so they will be deleted
+ % if opTomo is deleted.
+ op.astra_handle = opTomo_helper_handle([sino_id, ...
+ vol_id, proj_id, fp_alg_id, bp_alg_id]);
+
+ op.fp_alg_id = fp_alg_id;
+ op.bp_alg_id = bp_alg_id;
+ op.sino_id = sino_id;
+ op.vol_id = vol_id;
+ else
+ % 3D
+ % only gpu/cuda code for 3D
+ if ~gpuEnabled
+ error(['Only type ' 39 'cuda' 39 ' is supported ' ...
+ 'for 3D geometries.'])
+ end
+
+ % create a function handle
+ op.funHandle = @opTomo_intrnl3D;
+ end
+
+
+ % pass object properties
+ op.proj_size = proj_size;
+ op.vol_size = vol_size;
+ op.proj_geom = proj_geom;
+ op.vol_geom = vol_geom;
+ op.cflag = false;
+ op.sweepflag = false;
+
+ end % opTomo - constructor
+
+ end % methods - public
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Methods - protected
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ methods( Access = protected )
+
+ % multiplication
+ function y = multiply(op,x,mode)
+
+ % ASTRA cannot handle sparse vectors
+ if issparse(x)
+ x = full(x);
+ end
+
+ % convert input to single
+ if isa(x, 'single') == false
+ x = single(x);
+ end
+
+ % the multiplication
+ y = op.funHandle(op, x, mode);
+
+ % make sure output is column vector
+ y = y(:);
+
+ end % multiply
+
+ end % methods - protected
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Methods - private
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ methods( Access = private )
+
+ % 2D projection code
+ function y = opTomo_intrnl2D(op,x,mode)
+
+ if mode == 1
+ % X is passed as a vector, reshape it into an image.
+ x = reshape(x, op.vol_size);
+
+ % Matlab data copied to ASTRA data
+ astra_mex_data2d('store', op.vol_id, x);
+
+ % forward projection
+ astra_mex_algorithm('iterate', op.fp_alg_id);
+
+ % retrieve Matlab array
+ y = astra_mex_data2d('get_single', op.sino_id);
+ else
+ % X is passed as a vector, reshape it into a sinogram.
+ x = reshape(x, op.proj_size);
+
+ % Matlab data copied to ASTRA data
+ astra_mex_data2d('store', op.sino_id, x);
+
+ % backprojection
+ astra_mex_algorithm('iterate', op.bp_alg_id);
+
+ % retrieve Matlab array
+ y = astra_mex_data2d('get_single', op.vol_id);
+ end
+ end % opTomo_intrnl2D
+
+
+ % 3D projection code
+ function y = opTomo_intrnl3D(op,x,mode)
+
+ if mode == 1
+ % X is passed as a vector, reshape it into an image
+ x = reshape(x, op.vol_size);
+
+ % initialize output
+ y = zeros(op.proj_size, 'single');
+
+ % link matlab array to ASTRA
+ vol_id = astra_mex_data3d_c('link', '-vol', op.vol_geom, x, 0);
+ sino_id = astra_mex_data3d_c('link', '-sino', op.proj_geom, y, 1);
+
+ % initialize fp algorithm
+ cfg = astra_struct('FP3D_CUDA');
+ cfg.ProjectionDataId = sino_id;
+ cfg.VolumeDataId = vol_id;
+
+ alg_id = astra_mex_algorithm('create', cfg);
+
+ % forward projection
+ astra_mex_algorithm('iterate', alg_id);
+
+ % cleanup
+ astra_mex_data3d('delete', vol_id);
+ astra_mex_data3d('delete', sino_id);
+ else
+ % X is passed as a vector, reshape it into projection data
+ x = reshape(x, op.proj_size);
+
+ % initialize output
+ y = zeros(op.vol_size,'single');
+
+ % link matlab array to ASTRA
+ vol_id = astra_mex_data3d_c('link', '-vol', op.vol_geom, y, 1);
+ sino_id = astra_mex_data3d_c('link', '-sino', op.proj_geom, x, 0);
+
+ % initialize bp algorithm
+ cfg = astra_struct('BP3D_CUDA');
+ cfg.ProjectionDataId = sino_id;
+ cfg.ReconstructionDataId = vol_id;
+
+ alg_id = astra_mex_algorithm('create', cfg);
+
+ % backprojection
+ astra_mex_algorithm('iterate', alg_id);
+
+ % cleanup
+ astra_mex_data3d('delete', vol_id);
+ astra_mex_data3d('delete', sino_id);
+ end
+ end % opTomo_intrnl3D
+
+ end % methods - private
+
+end % classdef
diff --git a/matlab/tools/opTomo_helper_handle.m b/matlab/tools/opTomo_helper_handle.m
new file mode 100644
index 0000000..d9be51f
--- /dev/null
+++ b/matlab/tools/opTomo_helper_handle.m
@@ -0,0 +1,29 @@
+classdef opTomo_helper_handle < handle
+ %ASTRA.OPTOMO_HELPER_HANDLE Handle class around an astra identifier
+ % Automatically deletes the data when object is deleted.
+ % Multiple id's can be passed as an array as input to
+ % the constructor.
+
+ properties
+ id
+ end
+
+ methods
+ function obj = opTomo_helper_handle(id)
+ obj.id = id;
+ end
+ function delete(obj)
+ for i = 1:numel(obj.id)
+ % delete any kind of object
+ astra_mex_data2d('delete', obj.id(i));
+ astra_mex_data3d('delete', obj.id(i));
+ astra_mex_algorithm('delete', obj.id(i));
+ astra_mex_matrix('delete', obj.id(i));
+ astra_mex_projector('delete', obj.id(i));
+ astra_mex_projector3d('delete', obj.id(i))
+ end
+ end
+ end
+
+end
+
diff --git a/samples/matlab/s017_opTomo.m b/samples/matlab/s017_opTomo.m
new file mode 100644
index 0000000..4886cc5
--- /dev/null
+++ b/samples/matlab/s017_opTomo.m
@@ -0,0 +1,44 @@
+% load a phantom image
+im = phantom(256);
+% and flatten it to a vector
+x = im(:);
+
+%% Setting up the geometry
+% projection geometry
+proj_geom = astra_create_proj_geom('parallel', 1, 256, linspace2(0,pi,180));
+% object dimensions
+vol_geom = astra_create_vol_geom(256,256);
+
+%% Generate projection data
+% Create the Spot operator for ASTRA using the GPU.
+W = opTomo('cuda', proj_geom, vol_geom);
+
+p = W*x;
+
+% reshape the vector into a sinogram
+sinogram = reshape(p, W.proj_size);
+imshow(sinogram, []);
+
+
+%% Reconstruction
+% We use a least squares solver lsqr from Matlab to solve the
+% equation W*x = p.
+% Max number of iterations is 100, convergence tolerance of 1e-6.
+y = lsqr(W, p, 1e-6, 100);
+
+% the output is a vector, so we reshape it into an image
+reconstruction = reshape(y, W.vol_size);
+
+subplot(1,3,1);
+imshow(reconstruction, []);
+title('Reconstruction');
+
+subplot(1,3,2);
+imshow(im, []);
+title('Ground truth');
+
+% The transpose of the operator corresponds to the backprojection.
+backProjection = W'*p;
+subplot(1,3,3);
+imshow(reshape(backProjection, W.vol_size), []);
+title('Backprojection');
--
cgit v1.2.3
From e5bbcf730e1f5049b3dfd7f0ee46a7e477042231 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 15 May 2015 13:31:01 +0200
Subject: Add note about Spot to sample script
---
samples/matlab/s017_opTomo.m | 8 ++++++++
1 file changed, 8 insertions(+)
diff --git a/samples/matlab/s017_opTomo.m b/samples/matlab/s017_opTomo.m
index 4886cc5..e23b53d 100644
--- a/samples/matlab/s017_opTomo.m
+++ b/samples/matlab/s017_opTomo.m
@@ -1,3 +1,11 @@
+% This sample illustrates the use of opTomo.
+%
+% opTomo is a wrapper around the FP and BP operations of the ASTRA Toolbox,
+% to allow you to use them as you would a matrix.
+%
+% This class requires the Spot Linear-Operator Toolbox to be installed.
+% You can download this at http://www.cs.ubc.ca/labs/scl/spot/
+
% load a phantom image
im = phantom(256);
% and flatten it to a vector
--
cgit v1.2.3
From a6bdde8d391566972642090cb71cf752da92f69a Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 15 May 2015 13:43:17 +0200
Subject: Add copyright header
---
samples/matlab/s017_opTomo.m | 10 ++++++++++
1 file changed, 10 insertions(+)
diff --git a/samples/matlab/s017_opTomo.m b/samples/matlab/s017_opTomo.m
index e23b53d..891a93d 100644
--- a/samples/matlab/s017_opTomo.m
+++ b/samples/matlab/s017_opTomo.m
@@ -1,3 +1,13 @@
+% -----------------------------------------------------------------------
+% 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 sample illustrates the use of opTomo.
%
% opTomo is a wrapper around the FP and BP operations of the ASTRA Toolbox,
--
cgit v1.2.3
From 732a647d658230b682c6eaf3b61e3ea34af9cdbc Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 15 May 2015 13:56:04 +0200
Subject: Create python sample for OpTomo
---
samples/python/s017_OpTomo.py | 61 +++++++++++++++++++++++++++++++++++++++++++
1 file changed, 61 insertions(+)
create mode 100644 samples/python/s017_OpTomo.py
diff --git a/samples/python/s017_OpTomo.py b/samples/python/s017_OpTomo.py
new file mode 100644
index 0000000..967fa64
--- /dev/null
+++ b/samples/python/s017_OpTomo.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 <http://www.gnu.org/licenses/>.
+#
+#-----------------------------------------------------------------------
+
+import astra
+import numpy as np
+import scipy.sparse.linalg
+
+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('cuda',proj_geom,vol_geom)
+
+# construct the OpTomo object
+W = astra.OpTomo(proj_id)
+
+sinogram = W * P
+sinogram = sinogram.reshape([180, 384])
+
+import pylab
+pylab.gray()
+pylab.figure(1)
+pylab.imshow(P)
+pylab.figure(2)
+pylab.imshow(sinogram)
+
+# Run the lsqr linear solver
+output = scipy.sparse.linalg.lsqr(W, sinogram.flatten(), iter_lim=150)
+rec = output[0].reshape([256, 256])
+
+pylab.figure(3)
+pylab.imshow(rec)
+pylab.show()
+
+# Clean up.
+astra.projector.delete(proj_id)
--
cgit v1.2.3