summaryrefslogtreecommitdiffstats
path: root/samples/s015_fp_bp.m
blob: 815f4780735294c70ac94fb93736fd46a94e16ab (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
%------------------------------------------------------------------------
% This file is part of the
% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
%
% Copyright: iMinds-Vision Lab, University of Antwerp
% License: Open Source under GPLv3
% Contact: mailto:astra@ua.ac.be
% Website: http://astra.ua.ac.be
%------------------------------------------------------------------------


% 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