%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 2: Define wide-field source and detector in MMCLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear cfg newcfg;
%== 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

%% load and visualize the original mesh 
load Digimouse_Mesh_1L;

%== Extract the outermost surface mesh of the whole mesh domain
external_face=volface(elem(:,1:4));

%== Plot the surface mesh using transparent face color, this allows you to
%see the other tissues inside
figure(1);
skin=plotmesh(node,external_face);
    skin.FaceAlpha=0;
    skin.EdgeAlpha=0.2;
hold on;

%== Extract the tetrahedral element of all tissue types except for skin
tissue_elem=elem(elem(:,5)~=1,:);
%--------------------------------------------------------------------------
%%% Hands-on exercise: adjust the label to visualize your desired tissue, 
%%% e.g. For bone: tissue_elem=elem(elem(:,5)==2,:);   
%--------------------------------------------------------------------------

%== Plot the tissues inside the skin
plotmesh(node,tissue_elem);

set(gca,'XLim',[0 35],'YLim',[0 93],'ZLim',[0 21]);
set(gca,'looseInset',[0 0 0 0]);
axis off;
view([-75.9 4.4]);

%== Create a 3D animation for the mesh 
axis vis3d;
for i=1:360
    view([-75.9+i 4.4]);
    pause(0.025);
end

%% Define mesh domain for MMCLAB configuration
cfg.node = node;
cfg.elem = elem(:,1:4);
cfg.elemprop = elem(:,5);
cfg.face = face(:,1:3);

%== Assume homogeneous tissue for simplicity, all labeled using 1
cfg.elemprop(:) = 1;

%== Define optical properties for each tissue type
%         mua(1/mm) mus(1/mm)  g      n
cfg.prop=[   0          0      1      1;       %<- tissue label 0: background    
           0.0191      6.6    0.9    1.37];    %<- tissue label 1: 

%% Define external wide-field source
cfg.srctype = 'fourier';    % define an SFDI source
cfg.srcpos=[8.5 50 25];     % one corner of the illumination area
kx=3;                       % wave number in the x-dir
ky=0;                       % wave number in the y-dir
xphase=pi/3;                % phase offset in the x-dir, must < 2pi
yphase=0;                   % phase offset in the y-dir, must < 2pi
cfg.srcparam1=[20 0 0 kx+xphase/(2*pi)];    % location offset for the 2nd corner
cfg.srcparam2=[0 20 0 ky+yphase/(2*pi)];    % location offset for the 3rd corner
cfg.srcdir=[0 0 -1];        % source direction
%--------------------------------------------------------------------------
%%% Hands-on exercise: change source position, SFDI source spatial frequence   
%%% or phase shift along either x or y direction and rerun
%--------------------------------------------------------------------------

%% Add external wide-field illumination
newcfg = cfg;

%== Define a structure for the wide-field source
srcdef=struct('srctype',cfg.srctype,'srcpos',cfg.srcpos,'srcdir',cfg.srcdir,...
        'srcparam1',cfg.srcparam1,'srcparam2',cfg.srcparam2);

%== Merge the original mouse mesh with external source domain
[newcfg.node,newcfg.elem] = mmcaddsrc(cfg.node,[cfg.elem cfg.elemprop],srcdef);

figure(2);
subplot(1,2,1);plotmesh(newcfg.node,newcfg.elem);title('After merged with the source domain');
subplot(1,2,2);plotmesh(newcfg.node,newcfg.elem,'y>40');title('After merged with the source domain, y>40');

%% Add external wide-field detection

%== Define a structure for a wide-field planar detector
detdef=struct('srctype','planar','srcpos',[8.5 50 0],'srcdir',[0 0 1],...
                'srcparam1',[20 0 0 0],'srcparam2',[0 20 0 0]);

%== Merge the exteneded mouse mesh with external wide-field detection domain
[newcfg.node,newcfg.elem]=mmcadddet(newcfg.node,newcfg.elem,detdef);

figure(3);
subplot(1,2,1);plotmesh(newcfg.node,newcfg.elem);title('After merged with the detection domain');
subplot(1,2,2);plotmesh(newcfg.node,newcfg.elem,'y>40');title('After merged with the detection domain, y>40');

%== Update the element and tissue label information for merged mesh
newcfg.elemprop=newcfg.elem(:,5);
newcfg.elem=newcfg.elem(:,1:4);

%% time gate information 

newcfg.tstart=0;    % starting time: 0
newcfg.tend=5e-9;   % ending time: 5 ns
newcfg.tstep=1e-10; % time step: 0.1 ns, results in 50 time windows

%% Define other simulation parameters

%== Define simulated photon count
newcfg.nphoton=1e6;         % total number of simulated photons

%== Allocate random number seed (-1 for totally random)
newcfg.seed=12345678;       

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

%== Do reflection at internal & external boundaries (refractive index mismatch)
newcfg.isreflect=1;

%== Save exiting position and direction for detected photons
newcfg.issaveexit=1;

%% Preview the simulation mesh domain
figure(4);mcxpreview(newcfg);title('simulation domain preview');

%% Run mmc simulation

tic;
[mouse,detp,~,~] = mmclab(newcfg);
toc;

%% output data processing

%== extract time-resolved fluence rate
flux = mouse.data;
%--------------------------------------------------------------------------
%%% Hands-on exercise: fluence rate/flux is the default output of mmclab,
%%% other types of output can be generated by changing cfg.outputtype
%%% e.g. cfg.outputtype='energy' for energy deposition
%--------------------------------------------------------------------------

%== Integrate time-resolved fluence rate to continuous-wave fluence
mousecw = sum(flux,2)*newcfg.tstep;

% Plot CW fluence on the original mesh
figure(5);
plotmesh([cfg.node,mousecw(1:length(cfg.node))],cfg.elem,'linestyle','none');colorbar;
title('CW fluence [J/mm^2] on the original mesh')