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.

3 files changed, 353 insertions, 0 deletions

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');