DigitalBreastTomosynthesis/ProjectionNoising.m

%% Author: Aleks Prša
% Date: July, 2022
% alex.prsa@gmail.com
% =========================================================================
%{
% 
%     DESCRIPTION:
% 
%     Program for noising projection images.
%
%     Faculty of Mathematic and Physics
%     University in Ljubljana
%     Ljubljana, Slovenia
%
%}
% =========================================================================
%%                      Choose initial number of photons

close all;clear;clc

It = 3;                  % Anode current - number of photon

%% Load data

path_User = userpath;
path_ProjData = uigetdir(path_User);

addpath(genpath('Functions'));
addpath(genpath('Parameters'));
ParameterSettings_MicroInsertion

% Load projection data  
if(path_ProjData == 0)
    fprintf('Cancelled by user \n');
	return;
else
    userpath(path_ProjData);

    fprintf('Reading DICOMs   ')
    answer = 'Noised';
    [dataProj,infoDicom] = readDicom(path_ProjData,parameter,answer);
    parameter.bitDepth = infoDicom(:,1).BitDepth;
end

%% Projection noising

fprintf('\nProjection noising       ')
for i=1:size(dataProj,3)
    fprintf('\b\b\b\b\b%2d/%2d', i, size(dataProj,3))                           

    attenuation(:,:,i) = (dataProj(:,:,i) - parameter.mu) / (parameter.eo*parameter.Io);
    attenuation(attenuation < 0) = 0;
    attenuation(attenuation > 1) = 1;

    lambda = attenuation(:,:,i) .* (parameter.alpha * It);
    newProjection(:,:,i) = parameter.eo * poissrnd(lambda,parameter.nv,parameter.nu) + parameter.mu;
end

%% Save data

filestring = ['Output',filesep,'Noised',filesep,'Projections'];   
    if(~exist(filestring))
        mkdir(filestring)
    end
    
fprintf('\nWriting DICOM projections   ')
writeDicomProj(newProjection,infoDicom,filestring,parameter);

fprintf('\nSaving projections... \n')
save(['Output',filesep,'Noised',filesep,'NoisedProjections.mat'],'dataProj','-v7.3');
infoDicom = [];

fprintf('\nFINISHED - Data stored in "Noised" folder\n')

beep