%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 4: a replay example using MMCLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
%== add path to the softwares
if ispc
    addpath(genpath('../../Software/win64/MCXStudio'));
elseif isunix
    addpath(genpath('../../Software/linux64/MCXStudio'));
elseif ismac
    addpath(genpath('../../Software/osx64/MCXStudio'));
end

%% Define simulation parameters for the initial MC
clear cfg detp newcfg

cfg.nphoton=1e6;                % Define simulated photon count
cfg.seed=1648335518;            % Allocate random number seed (-1 for totally random)

%== Generate tetrahedral mesh for a homogeneous cubic domain
[cfg.node,cfg.elem]=genT6mesh(0:2:60,0:2:60,0:2:60);

%== Last column of elem is media type
cfg.elemprop=ones(size(cfg.elem,1),1);

%== Pencil beam source position and direction
cfg.srcpos=[30.1,30.2,0];       
cfg.srcdir=[0 0 1];             

%% Define time-gate settings
cfg.tstart=0;       % starting time: 0
cfg.tend=5e-9;      % ending time: 5 ns
cfg.tstep=5e-9;     % time step: 0.1 ns, results a single time window

%== Define optical properties for each tissue type
% simulate the optical properties of skull and gray-matter of a human brain model
% see http://mcx.sourceforge.net/cgi-bin/index.cgi?MMC/Colin27AtlasMesh
%         mua(1/mm) mus(1/mm)  g      n
cfg.prop=[   0          0      1      1;       %<- tissue label 0: background    
           0.005       1.0    0.01    1];      %<- tissue label 1: 

%== "T" to print timing information, "P" to show progress bar
cfg.debuglevel='TP';

%== Disable reflection at internal & external boundaries as there is no refractive index mismatch
cfg.isreflect=0;

%== Define detector position and radius (in length unit [unitinmm], not in mm unit!)
cfg.detpos=[30. 20. 0. 1.
            30. 40. 0. 1.
            20. 30. 0. 1.
            40. 30. 0. 1.];
        
%== Save detected photon exit position and angles
cfg.issaveexit=1;

%== Save detected photon seeds to replay later
cfg.issaveseed=1;

%== preprocessing of the mesh to get the missing fields
newcfg=mmclab(cfg,'prep');

%% Preview the domain
figure(1);mcxpreview(cfg);title('simulation domain preview');

%% Initial MC simulation
[cube, detp, ncfg, seeds]=mmclab(newcfg); 
% detp: detected photon information for detected photon weight and time of flight
% seeds: random number seeds for all detected photons

clear newdetp
%% Set up replay MC for Jmua

%== replay photons detected by det#1, i.e. build Jacobian for
% source and detector #1
newcfg.replaydet=1;     
%--------------------------------------------------------------------------
%%% Hands-on exercise: change the detector ID and rerun to build Jacobian
%%% for other source-detector pairs, such as detector #2,3,4
%--------------------------------------------------------------------------

%== Extract random number seeds for photons detected by detector #1
newcfg.seed=seeds.data(:,detp.data(1,:)==newcfg.replaydet);

%== Extract detected photon output for photon detected by detector #1
newdetp.detid=detp.detid(detp.detid==newcfg.replaydet,:);
    newdetp.nscat=detp.nscat(detp.detid==newcfg.replaydet,:);
    newdetp.ppath=detp.ppath(detp.detid==newcfg.replaydet,:);
    newdetp.p=detp.p(detp.detid==newcfg.replaydet,:);
    newdetp.v=detp.v(detp.detid==newcfg.replaydet,:);
    newdetp.w0=detp.w0(detp.detid==newcfg.replaydet,:);
    newdetp.data=detp.data(:,detp.detid==newcfg.replaydet);
    
%== Calculate the detected photon weight and time of flight using the partial 
% path output and optical properties
newcfg.replayweight=mmcdetweight(newdetp,newcfg.prop);
newcfg.replaytime=mmcdettime(newdetp,newcfg.prop);

%== disable normalization
newcfg.isnormalized=0;

%% Build Jmua (Absorption Jacobian)
newcfg.outputtype='wl';

%== Replay the detected photons to build absorption Jacobian for source and detector #1
[cube2,detp2,~,~]=mmclab(newcfg);

%== the two detected photon arrays should be the same. however, because
% the program uses multi-threading, the orders may be different
if(all(ismember(round(newdetp.data'*1e10)*1e-10,round(detp2.data'*1e10)*1e-10,'rows')))
   disp('replay is successful :-)');
else
   disp('replay failed :-(');
end

%== plot the Jacobian 3D profiles
figure(2);
plotmesh([cfg.node,log10(cube2.data)],cfg.elem(:,1:4),'z=0.2','linestyle','none');
hold on;
plotmesh([cfg.node,log10(cube2.data)],cfg.elem(:,1:4),'x=30','linestyle','none');
title('Absorption Jacobian');
view(3)
set(gca,'xlim',[0 60])
set(gca,'ylim',[0 60])
set(gca,'zlim',[0 60])
shading interp

%% Build Jmus (Scattering Jacobian)
newcfg.outputtype='wp';    

%== Replay the detected photons to build scattering Jacobian for source and detector #1
[cube3,detp3,~,~]=mmclab(newcfg);

%== the two detected photon arrays should be the same. however, because
% the program uses multi-threading, the orders may be different
if(all(ismember(round(newdetp.data'*1e10)*1e-10,round(detp3.data'*1e10)*1e-10,'rows')))
   disp('replay is successful :-)');
else
   disp('replay failed :-(');
end

%== plot the Jacobian 3D profiles
figure(3);
plotmesh([cfg.node,log10(cube3.data)],cfg.elem(:,1:4),'z=0.2','linestyle','none');
hold on;
plotmesh([cfg.node,log10(cube3.data)],cfg.elem(:,1:4),'x=30','linestyle','none');
title('Scattering Jacobian');
view(3)
set(gca,'xlim',[0 60])
set(gca,'ylim',[0 60])
set(gca,'zlim',[0 60])
shading interp