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 ++++
2 files changed, 309 insertions(+)
create mode 100644 matlab/tools/opTomo.m
create mode 100644 matlab/tools/opTomo_helper_handle.m
(limited to 'matlab/tools')
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
+
--
cgit v1.2.3