%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Summary:
% This script generates a phantom for heterogeneity analysis, as specified
% by the parameters below.
% See also article DOI: 10.2967/jnumed.116.181859
% !! needs files from phantom_lesions Git !!
%
%   cd('E:\010-work\003-localGit\intra_hetero\genPhantom');
% Matlab version:
% 2015a
% Initially created: 5/15/2017 by Peter Ferjancic
% last updated: 5/15/2017
% Version: 1.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function get_RTset_specific3D
imsize=[32, 32, 16];

bckgrnd = 0;
bckstd  = 20;
smoothR=1;
plot='n';   % plot='y';
% list all the lesion params you want
% y,x,z r,intensity, std of lesion
llist = ...
    [16,16,8, 8, 1000, 50.0;
     16,16,8, 3, 100, 1.0
     16,8,8, 3, 100, 1.0
    ];

% -- no touchy fishy after this line --

%% create canvas (determenistic phantom image)
canvas  = bckgrnd*ones(imsize);
stds    = bckstd *ones(imsize);
contourTarget = zeros(imsize);
contourROI = zeros(imsize);

for i=1:size(llist,1)
    [canvas, stds] = supp_addSphere(canvas, stds, llist(i,1),llist(i,2),llist(i,3),llist(i,4), llist(i,5), llist(i,6));
end

[contourTarget, ~] = supp_addSphere(contourTarget, stds, llist(2,1),llist(2,2),llist(2,3),llist(2,4), 1, 0);
[contourROI, ~] = supp_addSphere(contourROI, stds, llist(3,1),llist(3,2),llist(3,3),llist(3,4), 1, 0);

%% Make it random and smooth it
image=zeros(size(canvas));
for yi=1:size(canvas, 1)
    for xi=1:size(canvas, 2)
        for zi=1:size(canvas, 3)
            image(yi,xi,zi) = canvas(yi,xi,zi) + stds(yi,xi,zi)*randn;
        end
    end
end
imageF=imgaussfilt(image,smoothR);

if plot == 'y'
    %%  plot result
    subplot(2,2,1)
    imagesc(canvas, [0,15])
    title('Signal')
    colorbar

    subplot(2,2,2)
    imagesc(stds, [0,10])
    title('STD')
    colorbar

    subplot(2,2,3)
    imagesc(image)
    title('Added noise')
    colorbar
    colormap('hot')

    subplot(2,2,4)
    imagesc(imageF)
    title('Final phantom')
    colorbar
    colormap('hot')
end

pixelspacing = [0.16,0.16,0.16];
origin = [-1.1,-2.2,-3.3];
encoding = 'ascii';

filename = 'C:\010-work\003_localGit\WiscPlan_v2\PhantomData\test_ct.nrrd';
matrix = imageF;
nrrdWriter(filename, matrix, pixelspacing, origin, encoding);

filename = 'C:\010-work\003_localGit\WiscPlan_v2\PhantomData\test_seg.nrrd';
matrix = contourTarget;
nrrdWriter(filename, matrix, pixelspacing, origin, encoding);

filename = 'C:\010-work\003_localGit\WiscPlan_v2\PhantomData\test_ROI.nrrd';
matrix = contourROI;
nrrdWriter(filename, matrix, pixelspacing, origin, encoding);

end