summaryrefslogtreecommitdiffstats
path: root/matlab/tools/astra_add_noise_to_sino.m
blob: 957b5d01698d0a41ee2dbbab774ecada922d3682 (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
function sinogram_out = astra_add_noise_to_sino(sinogram_in,I0)

%--------------------------------------------------------------------------
% sinogram_out = astra_add_noise_to_sino(sinogram_in,I0)
%
% Add poisson noise to a sinogram.
%
% sinogram_in: input sinogram, can be either MATLAB-data or an
% astra-identifier.  In the latter case, this operation is inplace and the
% result will also be stored in this data object.
% I0: background intensity, used to set noise level, lower equals more
% noise
% sinogram_out: output sinogram in MATLAB-data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the ASTRA Toolbox
% 
% Copyright: 2010-2018, imec Vision Lab, University of Antwerp
%            2014-2018, CWI, Amsterdam
% License: Open Source under GPLv3
% Contact: astra@astra-toolbox.com
% Website: http://www.astra-toolbox.com/
%--------------------------------------------------------------------------

if numel(sinogram_in) == 1
	sinogramRaw = astra_mex_data2d('get', sinogram_in);
else
	sinogramRaw = sinogram_in;
end

% scale to [0,1]
max_sinogramRaw = max(sinogramRaw(:));
sinogramRawScaled = sinogramRaw ./ max_sinogramRaw;
% to detector count
sinogramCT = I0 * exp(-sinogramRawScaled);
% add poison noise
sinogramCT_A = sinogramCT * 1e-12;
sinogramCT_B = double(imnoise(sinogramCT_A, 'poisson'));
sinogramCT_C = sinogramCT_B * 1e12;
% to density
sinogramCT_D = sinogramCT_C / I0;
sinogram_out = -max_sinogramRaw * log(sinogramCT_D);

if numel(sinogram_in) == 1
	astra_mex_data2d('store', sinogram_in, sinogram_out);
end