%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear cfg;

%% preparing the input data
% set seed to make the simulation repeatible
rand_seed = randi([1 2^31-1], 1, 1);
cfg.seed=rand_seed; 

cfg.nphoton=5e7;

% load brain image
load('brain.mat')
brain = uint8(brain);
cfg.vol = brain;

% define the source position
cfg.srcpos=[85 87 204];
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]
% [background; gray matter; white matter; CSF; bone; skull]
cfg.prop = [0 0 1 1;
    0.02 9.0 0.89 1.37;
    0.08 40.9 0.84 1.37;
    0.004 0.009 0.89 1.37;
    0.019 7.8   0.89 1.37;
    0.019 7.8   0.89 1.37];

% 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 with boundary reflection enabled
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]=mcxlab(cfg);
toc;

%% Plot
n = 1;
caxis = [0 8];
brain_slice = squeeze(brain(85,:,:));
brain_slice = rot90(brain_slice);
data = squeeze(f2.data(85,:,:,n));
data = rot90(data);
figure,imagesc(squeeze(log10(single(brain_slice)+data)),caxis),colormap jet
axis off
axis equal
title(['time=' num2str(0.5*n) ' ns'])

n = 2;
data = squeeze(f2.data(85,:,:,n));
data = rot90(data);
figure,imagesc(squeeze(log10(single(brain_slice)+data)),caxis),colormap jet
axis off
axis equal
title(['time=' num2str(0.5*n) ' ns'])

n = 3;
data = squeeze(f2.data(85,:,:,n));
data = rot90(data);
figure,imagesc(squeeze(log10(single(brain_slice)+data)),caxis),colormap jet
axis off
axis equal
title(['time=' num2str(0.5*n) ' ns'])