%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  matlab script to produce the validation plots (Fig. 4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (exist('loadjson') ~= 2)
    error('you need to first download and add path to jsonlab (https://github.com/fangq/jsonlab) in order to run this script');
end

c0 = 299792458000;
srcpos = [30 30 1];
detpos = [30 15 10];
twin = [5e-11:1e-10:5e-9];

% ----------------------------------------------------------------
%  Plot time-domain comparison at a given position
% ----------------------------------------------------------------

figure;
hold on;

mcx = loadjson('semi_infinite.jnii', [60 60 60 50]);
mcx = mcx.NIFTIData;
semilogy((1:50) / 10, tddiffusion(0.005, 1, c0, 0, srcpos, detpos, twin), 'r');
semilogy((1:50) / 10, squeeze(mcx(30, 14, 9, :)), 'o');

mcxb = loadjson('semi_infinite_b.jnii', [60 60 60 50]);
mcxb = mcxb.NIFTIData;
semilogy((1:50) / 10, tddiffusion(0.005, 1, c0 / 1.37, 0.493, srcpos, detpos, twin), 'r--');
semilogy((1:50) / 10, squeeze(mcxb(30, 14, 9, :)), '+');

set(gca, 'ylim', [1 10^8.5]);

legend('Diffusion (no reflection)', 'MCX (no reflection)', 'Diffusion (with reflection)', 'MCX (with reflection)');
legend boxoff;
set(gca, 'yscale', 'log');
set(gca, 'fontsize', 18);
xlabel('time (ns)');
ylabel('Fluence in 1/(m^2s)');
box on;
saveas(gcf, 'td_validate.fig');
print -depsc2 td_validate.eps;

% ----------------------------------------------------------------
%  Plot CW comparison along a line
% ----------------------------------------------------------------

figure;
hold on;

cwmcx = sum(mcx, 4);
srcpos = [0 0 0];
detpos = [0.5 + (0:29)', 0.5 * ones(30, 1), 0.5 * ones(30, 1)];
phicw = cwdiffusion(0.005, 1, 0, srcpos, detpos);
phimcx = cwmcx(30, :, 1);

cwmcxb = sum(mcxb, 4);
phicwb = cwdiffusion(0.005, 1, 0.493, srcpos, detpos);
phimcxb = (cwmcxb(30, :, 1));
h = semilogy(1:30, phicw, 'r-', 1:30, phimcx(30:59) * 1e-10, 'o', 1:30, phicwb, 'r--', 1:30, phimcxb(30:59) * 1e-10, 'b+');
box on;
set(gca, 'yscale', 'log');

legend('Diffusion (no reflection)', 'MCX (no reflection)', 'Diffusion (with reflection)', 'MCX(with reflection)');
legend boxoff;
set(gca, 'fontsize', 18);
xlabel('x (mm)');
ylabel('Fluence in 1/(m^2s)');
saveas(gcf, 'td_validate_cw.fig');
print -depsc2 td_validate_cw.eps;

% ----------------------------------------------------------------
%  Plot CW comparison at a cross-section
% ----------------------------------------------------------------

figure;
hold on;

clines = -1.5:-0.5:-4;
srcpos = [29.5 29.5 0.5];
[xi, yi] = meshgrid(0.5:59.5, 0.5:59.5);
detpos = [xi(:) 30 * ones(size(xi(:))) yi(:)];
phicw = reshape(cwdiffusion(0.005, 1.0, 0, srcpos, detpos), size(xi));

[c h2] = contour(xi, yi, log10(max(squeeze(phicw), 1e-8)), clines, 'k-');
[c h1] = contour(log10(max(squeeze(cwmcx(:, 30, :) * 1e-10)', 1e-8)), clines, 'k:');
legend('Diffusion', 'MCX');
legend boxoff;
box on;
set(gca, 'fontsize', 18);
xlabel('x (mm)');
ylabel('z (mm)');
saveas(gcf, 'td_validate_cw_2d.fig');
print -depsc2 td_validate_cw_2d.eps;

% ----------------------------------------------------------------
%  Plot time-domain comparison at a cross-section
% ----------------------------------------------------------------

c0 = 299792458000;
twin = 5e-11:1e-10:5e-9;
mua = 0.005;
musp = 1;
Reff = 0;

srcpos = [29.5 29.5 0];
[xi, yi] = meshgrid(0.5:59.5, 0.5:59.5);

detpos = [xi(:) 29.5 * ones(size(xi(:))) yi(:)];
for i = 1:length(twin)
    phitd(:, :, i) = reshape(tddiffusion(mua, musp, c0, Reff, srcpos, detpos, twin(i)), size(xi));
end

figure;
hold on;

for i = 1:5:21
    halfmax = max(max(phitd(:, :, i))) / 2;
    [c h2] = contour(xi, yi, max(squeeze(phitd(:, :, i)), 1e-12), [halfmax halfmax], 'k-');
    halfmax = max(max(mcx(:, 30, :, i))) / 2;
    [c h2] = contour(max(squeeze(mcx(:, 30, :, i))', 1e-12), [halfmax halfmax], 'k:');
end
legend('Diffusion', 'MCX');
legend boxoff;
box on;
set(gca, 'fontsize', 18);
xlabel('x (mm)');
ylabel('z (mm)');
saveas(gcf, 'td_validate_2d.fig');

print -depsc2 td_validate_2d.eps;
