% in this script the spheres in the mathematical phantom are constructed
% they are placed at the location defined for the NEMA quality phantom
function image = constructSpheres(image, imageValue)
    imageSize= size(image);
    radius37 = 18.5;
    center37 = [256+57,256, 90];
    
    radius17 = 8.5;
    center17 = [256-57, 256, 90]

    radius28 = 14;
    center28 = [ 256+28, 256+49 ,90];

    radius22 = 11;
    center22 = [256-28, 256+49 ,90];

    radius13 = 7.5;
    center13 = [256-28, 256-49 ,90];

    radius10 = 5;
    center10 = [256+28, 256-49 ,90];
    % Create a logical image of a circle with specified
    % diameter, center, and image size.
    % First create the image.
    imageSizeX = imageSize(1);
    imageSizeY = imageSize(2);
    imageSizeZ = imageSize(3);
    [columnsInImage rowsInImage depth] = meshgrid(1:imageSizeX, 1:imageSizeY, 1:imageSizeZ);


    circlePixels37 = (rowsInImage - center37(2)).^2 ...
                 + (columnsInImage - center37(1)).^2 ...
                 + (depth - center37(3)).^2 <= radius37.^2;
    circlePixels17 = (rowsInImage - center17(2)).^2 ...
               + (columnsInImage - center17(1)).^2 ...
                 + (depth - center17(3)).^2 <= radius17.^2;
    circlePixels22 = (rowsInImage - center22(2)).^2 ...
                 + (columnsInImage - center22(1)).^2 ...
                 + (depth - center22(3)).^2 <= radius22.^2;

    circlePixels28 = (rowsInImage - center28(2)).^2 ...
                 + (columnsInImage - center28(1)).^2 ...
                 + (depth - center28(3)).^2 <= radius28.^2;        

    circlePixels13 = (rowsInImage - center13(2)).^2 ...
               + (columnsInImage - center13(1)).^2 ...
               + (depth - center13(3)).^2 <= radius13.^2; 

    circlePixels10 = (rowsInImage - center10(2)).^2 ...
               + (columnsInImage - center10(1)).^2 ...
               + (depth - center10(3)).^2 <= radius10.^2; 
    image(circlePixels37 ) = 1;     
    %image(circlePixels37 | circlePixels17 | circlePixels22 | circlePixels28 | circlePixels13 | circlePixels10) = 500;
end