%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
%
% In this example, we simulate a 4-layer brain model using MCXLAB.
% We will investigate the differences between the solutions with and 
% witout boundary reflections (both external and internal) and show
% you how to display and analyze the resulting data.
%
% This file is part of Monte Carlo eXtreme (MCX) URL:http://mcx.sf.net
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% only clear cfg to avoid accidentally clearing other useful data
clear cfg;

%% preparing the input data
% set seed to make the simulation repeatible
cfg.seed=hex2dec('623F9A9E'); 

cfg.nphoton=5e7;

% define a 4 layer structure
cfg.vol=ones(100,100,50);
cfg.vol(:,:,21:25)=2;
cfg.vol(:,:,26:35)=3;
cfg.vol(:,:,36:end)=4;
cfg.vol=uint8(cfg.vol);

% define the source position
cfg.srcpos=[50,50,0]+1;
cfg.srcdir=[0 0 1];

% use the brain optical properties defined at
% http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/CollinsAtlasMesh
% format: [mua(1/mm) mus(1/mm) g n]

cfg.prop=[0 0 1 1            % medium 0: the environment
   0.019 7.8   0.89 1.37     % medium 1: skin & skull
   0.004 0.009 0.89 1.37     % medium 2: CSF
   0.02  9.0   0.89 1.37     % medium 3: gray matter
   0.08 40.9   0.84 1.37];   % medium 4: white matter

% time-domain simulation parameters
cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=5e-10;

% GPU thread configuration
cfg.autopilot=1;
cfg.gpuid=1;

%% running simulation without boundary reflection
fprintf('running simulation ... this takes about 35 seconds on a GTX 470\n');
cfg.isreflect=0; % disable reflection at exterior boundary
tic;
f1=mcxlabcl(cfg);
toc;

%% running simulation with boundary reflection enabled
fprintf('running simulation ... this takes about 50 seconds on a GTX 470\n');
cfg.isreflect=1; % enable reflection at exterior boundary
cfg.isrefint=1;  % enable reflection at interior boundary too
cfg.issavedet=1; % enable recording partial pathlength of detected photons
cfg.detpos=[31 51 1 2];
tic;
[f2,det2]=mcxlabcl(cfg);
toc;

%% plot the results
figure
subplot(221);
contourf(log10(squeeze(sum(f1.data(:,51,:,:),4))'),1:0.5:8);
hold on
plot([0 100],[21 21],'--',[0 100],[26 26],'--',[0 100],[36 36],'--');
title('Fluence rate with no reflection');
set(gca,'clim',[1 8]);

subplot(222);
contourf(log10(squeeze(sum(f2.data(:,51,:,:),4))'),1:0.5:8);
hold on
plot([0 100],[21 21],'--',[0 100],[26 26],'--',[0 100],[36 36],'--');
title('Fluence rate with reflection at boundaries');
set(gca,'clim',[1 8]);

subplot(223);
hold on;
xi=1e9*((cfg.tstart:cfg.tstep:cfg.tend-cfg.tstep)+0.5*cfg.tstep);
for i=31:49
    semilogy(xi,squeeze(f2.data(51,i,1,:)),'color',[1-(i-25)/25 1-(i-25)/25 1]);
end
xlabel('time (ns)')
ylabel('Fluence rate (1/mm^2/s)')
set(gca,'yscale','log');
title('time point spread functions (TPSF)');

subplot(224);
hold on;
[hs1,c1]=hist(det2.ppath(find(det2.ppath(:,1)),1),200);
[hs2,c2]=hist(det2.ppath(find(det2.ppath(:,2)),2),200);
[hs3,c3]=hist(det2.ppath(find(det2.ppath(:,3)),3),200);
[hs4,c4]=hist(det2.ppath(find(det2.ppath(:,4)),4),200);
bar(c1,hs1,'edgecolor','none','facecolor','r');
bar(c2,hs2,'edgecolor','none','facecolor','b');
bar(c3,hs3,'edgecolor','none','facecolor','g');
bar(c4,hs4,'edgecolor','none','facecolor','k');
legend('tissue 1','tissue 2','tissue 3','tissue 4');
xlabel('partial pathlength (mm)');
title(sprintf('detected %d photons',size(det2.ppath,1)));
