%% Author: Aleks Prša
% Date: July, 2022
% alex.prsa@gmail.com
% =========================================================================
%{
% -------------------------------------------------------------------------
%                           phantom(parameter)
% -------------------------------------------------------------------------
%     DESCRIPTION:
%     This function create a 3D virtual phantom.
%  
%     INPUT:
%     - parameter = Parameter of all geometry
% 
%     OUTPUT:
%     - phantom = Virtual phantom attenuation data.
% 
%     ---------------------------------------------------------------------
%}
% =========================================================================
%% Create virtual phantom
function data3d = phantom(parameter)

data3d = zeros(parameter.nx, parameter.ny, parameter.nz);

xyzr = readmatrix('Mikrokalcinacije.txt');

for i = 1:size(xyzr, 1)
    mu = xyzr(i, 1);                    % Attenuation coefficient of microcalcification (1/mm)
    y0 = xyzr(i, 2);                    % y offset of microcalcification in pixels
    x0 = xyzr(i, 3);                    % x offset of microcalcification in pixels
    z0 = xyzr(i, 4);                    % z offset of microcalcification in slice number of original tomogram
    r = xyzr(i, 5) / 2000;              % Radius of microcalcification (mm)

    z0 = (z0-1) * 10;                   % Transform z offset from slice number to pixel value

    a = r / parameter.dx;               % Radius in x-axis
    b = r / parameter.dy;               % Radius in y-axis
    c = r / parameter.dz;               % Radius in z-axis

    a = a^2;
    b = b^2;
    c = c^2;

    for z = 1:parameter.nz
        for y = 1:parameter.ny
            for x = 1:parameter.nx
                if((x-x0)^2 / a + (y-y0)^2 / b + (z-z0)^2 / c <= 1) 
                    data3d(x,y,z) = data3d(x,y,z) + mu;
                end
            end
        end
    end
end

data3d = data3d - 0.03756;                              % Attenuation correction - microcalcification-H2O
data3d(data3d < 0) = 0;

%data3d(data3d == 0) = 0.0357776;                       % Attenuation coefficient of PMMA (1/mm) surrounding microcalcifications

data3d = flip(data3d, 1);

end