%%-----------------------------------------------------------------
%% Comparing widefield MMC and MCX in a homogeneous cubic domain
%  (see MMC paper Fig. 2)
%%-----------------------------------------------------------------
%
% In this example, we compare MMC and MCX using both simple and
% widefield sources.
%
% The cubic domain has a dimension of 60x60x60 mm with optical properties
% mua=0.001, mus=1, n=1.0 and g=0.01. The analytical solution for 
% pencilbeam source can be computed by the cwdiffusion (for CW solution) 
% and tddiffusion (for time-domain solutions) functions in the
% package of Monte Carlo eXtreme (MCX) under the mcx/utils directory.
%
%%-----------------------------------------------------------------

%addpath('/Please/add/path/to/mcx/utils/')
addpath('../../matlab/');
addpath('../../mmclab/');

%%-----------------------------------------------------------------
%% run mmclab for the 60x60x60 homogeneous cubic domain
%%-----------------------------------------------------------------

clear cfg;
clear all;

planarsrc=1; % set to 1 to compare planar src, set to 0 for pencil beam
inclusion=1; % set to 1 to add an inclusion in the domain

cfg.nphoton=3e7;
cfg.seed=1648335518;
if(~inclusion)
    [cfg.node,tmp,cfg.elem]=meshabox([0 0 0],[60 60 60], 2*sqrt(2), 2^3/6);
    cfg.elemprop=ones(size(cfg.elem,1),1);
    cfg.prop=[0 0 1 1;0.005 1.0 0.01 1.0];
else
    [node1,face1]=meshabox([0 0 0],[60 60 60], 2*sqrt(2), 2^3/6);
    [node2,face2]=meshabox([20 20 20],[40 40 40], 2*sqrt(2), 2^3/6);
    [node,face]=mergemesh(node1,face1,node2,face2);
    [cfg.node,cfg.elem]=surf2mesh(node,face,[],[],1,2^3/6,[1 1 1;30 30 30],[],0);
    cfg.elemprop=cfg.elem(:,5);
    cfg.elem=cfg.elem(:,1:4);
    cfg.prop=[0 0 1 1;0.005 1.0 0.01 1.0;0.01 0.2, 0.01 1.37];
end
cfg.srcpos=[30.1,30.2,0];
% cfg.srcdir=[0 0 1 -30]; % to test a diverging beam with a focus of -30
cfg.srcdir=[0 0 1];
cfg.tstart=0;
cfg.tend=5e-9;
cfg.tstep=1e-10;
cfg.debuglevel='TP';
cfg.isreflect=1;
cfg.detpos=[0 30 30 2];

% define wide-field planar source

if(planarsrc)   % when planar src is used, re-tessellate the mesh
    cfg.srctype='planar';
    cfg.srcpos=[10 10 -10];
    cfg.srcparam1=[40 0 0 0];
    cfg.srcparam2=[0 40 0 0];

    % the following block before "end" can be removed as mmclab
    % has this built-in
    srcdef=struct('srctype',cfg.srctype,'srcpos',cfg.srcpos,'srcdir',cfg.srcdir,...
                  'srcparam1',cfg.srcparam1,'srcparam2',cfg.srcparam2);
    if(size(cfg.elem,2)==5)
        cfg.elem=cfg.elem(:,1:4);
    end
    [cfg.node,cfg.elem] = mmcaddsrc(cfg.node,[cfg.elem cfg.elemprop],...
         mmcsrcdomain(srcdef,[min(cfg.node);max(cfg.node)]));
    cfg.elemprop=cfg.elem(:,5);
    cfg.elem=cfg.elem(:,1:4);
end
cfg=mmclab(cfg,'prep');

%%-----------------------------------------------------------------
%% run mmclab
%%-----------------------------------------------------------------

cube=mmclab(cfg);
%cube=mmclab(cfg,'sse'); % this is faster
cube=cube.data;

%%-----------------------------------------------------------------
%% run mcxlab for the same domain
%%-----------------------------------------------------------------

cfgx=cfg;
cfgx=rmfield(cfgx,{'node','elem','elemprop','debuglevel','facenb','evol','e0','isreoriented'});
dim=60;
cfgx.vol=ones(dim,dim,dim);
cfgx.vol=uint8(cfgx.vol);
cfgx.faststep=0;
if(inclusion)
    cfgx.vol(21:40,21:40,21:40)=2;
end
cfgx.nphoton=numel(cfgx.vol)/size(cfg.node,1)*cfg.nphoton; % for mcx, nphoton 
                                                       % needs to be bigger 
                                                       % to match the noise 
                                                       % of mmc
cfgx.autopilot=1;
cfgx.gpuid=1;
cfgx.issrcfrom0=1; % treating the lower-outer corner of the grid as [0 0 0]
                   %to match mmc

cubex=mcxlab(cfgx);
cubex=cubex.data;

%%-----------------------------------------------------------------
%% Make a cross-section
%%-----------------------------------------------------------------

c0=299792458000;

twin=[cfg.tstart+cfg.tstep/2:cfg.tstep:cfg.tend];
gates=length(twin);

cwcb=sum(cube,2);
cwcbx=sum(cubex,4);

[cutpos,cutvalue,facedata]=qmeshcut(cfg.elem(:,1:4),cfg.node,cwcb,'y=30.5');
                                         % slice at y=30.5 to offset half
                                         % grid to match mcx voxels

[xi,yi]=meshgrid(0.5:59.5,0.5:59.5);     % interpolate mmc solution at 
                                         % half-grid because mcx's voxel 
                                         % readings are at the center of 
                                         % the voxels
vi=griddata(cutpos(:,1),cutpos(:,3),cutvalue,xi,yi);
vix=squeeze(cwcbx(:,30,:))';

%%-----------------------------------------------------------------
%% generate a contour plot along y=30.5
%%-----------------------------------------------------------------
figure
hold on
clines = -1.5:-0.25:-8;
if(~planarsrc && ~inclusion)
    srcs=[30.1,30.2,0];
    dets=[xi(:) 30.2*ones(size(xi(:))) yi(:)];
    phicw=reshape(cwdiffusion(0.005, 1.0, 0, srcs, dets),size(xi));
    contour(xi,yi, log10(max(squeeze(phicw),1e-8)), -1.5:-0.5:-4, 'k-' ); 
                 % analytical is semi-infinite, only take the center part
end
contour(xi,yi,log10(abs(vi*cfg.tstep)),clines,'r:')
contour(xi,yi,log10(abs(vix*cfg.tstep)),clines,'b--')

if(inclusion)
    plot([20 20 40 40 20],[20 40 40 20 20],'k--');
end

axis equal
set(gca,'xlim',[1 60]);
set(gca,'fontsize',20);
xlabel('x (mm)');
ylabel('z (mm)');
if(planarsrc || inclusion)
    legend('MMC','MCX');
else
    legend('Diffusion','MMC','MCX');
end
legend boxoff;
box on;
