%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 3: Dual-grid MMC example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%% create surface mesh for three concentric spheres with radius = 10mm,
% 23mm and 25mm respectively
[nsph_10,fsph_10]=meshasphere([30 30 30],10,1.0);
[nsph_23,fsph_23]=meshasphere([30 30 30],23,2.0);
[nsph_25,fsph_25]=meshasphere([30 30 30],25,2.0);

%% Define simulation parameters
cfg.nphoton=3e6;            % Define simulated photon count
cfg.seed=1648335518;        % Allocate random number seed (-1 for totally random)
cfg.srcpos=[30.5,30.5,0];   % Default pencil beam source is used here
cfg.srcdir=[0 0 1];         % Source direction
cfg.tstart=0;               % Starting time: 0 
cfg.tend=5e-9;              % Ending time: 5 ns
cfg.tstep=5e-10;            % Time step: 0.5 ns, results in 10 time windows

%== Define optical properties for each tissue type
% simulate the optical properties similar to that 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.02       7.0     0.89   1.37;     %<- tissue label 1: similar to Skull
            0.004      0.009   0.89   1.37;     %<- tissue label 2: similar to CSF
            0.02       9.0     0.89   1.37;     %<- tissue label 3: gray-matter
            0.05       0.0     1.0    1.37];    %<- tissue label 4: non-scattering inclusion

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

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

%== Use the same base simulation parameter for multiple simulations
cfgs=[cfg cfg];

%% generate a coarse CDT volumetric mesh from triangular surfaces of 
% three spheres with an additional bounding box for DMMC

ISO2MESH_SESSION='dmmc1_';

%== Create mesh for a 60x60x60 bounding box
[nbox,ebox]=meshgrid6(0:60:60,0:60:60,0:60:60);
fbox=volface(ebox);

%== Merge three spheres with bounding box
[no,fc]=mergemesh(nsph_10,fsph_10,nsph_23,fsph_23,nsph_25,fsph_25,nbox,fbox);

%== Remove potential isolated nodes (nodes that are not included in any element)
[no,fc]=removeisolatednode(no,fc);

%== Set tetgen option to make the mesh as coarse as possible while achieving good mesh quality
% To learn more about it, please checkout: http://wias-berlin.de/software/tetgen/1.5/doc/manual/manual005.html
ISO2MESH_TETGENOPT='-A -q -Y'

%== Create a relatively coarse volumetric mesh
[cfgs(1).node,cfgs(1).elem]=surf2mesh(no,fc,[0 0 0],[60.1 60.1 60.1],1,100,[1 1 1;30 30 6;30 30 15;30 30 30]);

%== Activate 'DMMC' mode for cfgs(1)
cfgs(1).method='grid';

%% Generate a refined volumetric mesh from triangular surfaces of 
% three spheres with an additional bounding box for MMC
clear ISO2MESH_TETGENOPT
ISO2MESH_SESSION='dmmc2_';

%== Merge surface mesh of three spheres
[no,fc]=mergemesh(nsph_10,fsph_10,nsph_23,fsph_23,nsph_25,fsph_25);
[no,fc]=removeisolatednode(no,fc);

%== Generate volumetric mesh
[cfgs(2).node,cfgs(2).elem,face2]=surf2mesh(no,fc,[0 0 0],[60 60 60],1,30,[1 1 1 3;30 30 6 3;30 30 17 3;30 30 30 3],[],[1 1 1 1 1 1 1 1]);
[cfgs(2).node,cfgs(2).elem]=sortmesh(cfg.srcpos,cfgs(2).node,cfgs(2).elem,1:4);
%--------------------------------------------------------------------------
%%% Hands-on exercise: adjust maximum tetrahedral volume for each tissue   
%%% type and see how this impact mesh densities, simulation speed and 
%%% output accuracy for MMC
%--------------------------------------------------------------------------

%== Use 'havel' ray-tracer for cfgs(2)
cfgs(2).method='havel';

%% Compare DMMC mesh with MMC mesh

figure(1);
subplot(1,2,1);plotmesh(cfgs(1).node,cfgs(1).elem,'y>30');view(3);
subplot(1,2,2);plotmesh(cfgs(2).node,cfgs(2).elem,'y>30');view(3);

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

%% Run mmc simulation
tic;
phimmc=mmclab(cfgs);
toc;

%% Define mcxlab simulation configuration
clear cfg;

cfg.nphoton=1e7;

%== Define the simulation domain for MCX, which is equivalent to the mesh
% domain
[xi,yi,zi]=ndgrid(0.5:(60-0.5),0.5:(60-0.5),0.5:(60-0.5)); % centroid of each voxel element
dist=(xi-30).^2+(yi-30).^2+(zi-30).^2;                     % distance^2 from each voxel centroid to the center of circle

%== Create MCX simulation domain and update optical property labels to match the mesh domain
cfg.vol=ones(size(xi));                                    
    cfg.vol(dist<625)=2;    % 25^2=625
    cfg.vol(dist<529)=3;    % 23^2=625
    cfg.vol(dist<100)=4;    % 10^2=100 
    cfg.vol=uint8(cfg.vol);

%== Define position, direction of a pencil beam source
cfg.srcpos=[30.5,30.5,0];
cfg.srcdir=[0 0 1];
cfg.issrcfrom0=1;

%== Define optical properties for each tissue type
% simulate the optical properties similar to that 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.02       7.0     0.89   1.37;     %<- tissue label 1: similar to Skull
            0.004      0.009   0.89   1.37;     %<- tissue label 2: similar to CSF
            0.02       9.0     0.89   1.37;     %<- tissue label 3: gray-matter
            0.05       0.0     1.0    1.37];    %<- tissue label 4: non-scattering inclusion

%== Time-domain simulation parameters
cfg.tstart=0;       % starting time: 0
cfg.tend=5e-9;      % ending time: 5 ns
cfg.tstep=5e-10;    % time step: 0.1 ns, results in 10 time windows

%== Define GPU parameters
cfg.gpuid=1;        % use the first GPU
cfg.autopilot=1;    % assign the thread and block automatically

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

%== Grid size in mm
cfg.unitinmm=1; 

%% Running MCX simulations

USE_MCXCL=1;    % activate mcxcl
tic;
phimcx=mcxlab(cfg);
toc;

%% Extract fluence rate results from three simulation results

%== fluence rate output by MCX, integrate to CW fluence
phi_mcx=sum(phimcx.data,4)*cfg.tstep;

%== fluence rate output by MMC, integrate to CW fluence
phi_mmc=sum(phimmc(2).data,2)*cfg.tstep;

%== interpolate MMC node-based CW fluence to the centroid of voxel element 
node2=cfgs(2).node;
phi_mmc_grid=griddata(node2(:,1),node2(:,2),node2(:,3),phi_mmc,xi,yi,zi);

%== fluence rate output by DMMC, integrate to CW fluence
phi_dmmc=sum(phimmc(1).data,4)*cfg.tstep;

%== To ensure the output grid fully covers the mesh domain, voxel dimensions
%are extended along x,y,z directions by 1 grid respectively
phi_dmmc=phi_dmmc(1:end-1,1:end-1,1:end-1); 

%% compare CW fluence using contour plot along plane y=30.5

figure(3);
clines=[-10:0.5:10];
contourf(log10(abs(squeeze(phi_mmc_grid(:,31,:))')),clines,'linewidth',1.5);
hold on
contour(log10(abs(squeeze(phi_dmmc(:,31,:))')),clines,'linestyle','--','color','w','linewidth',1.5);
contour(log10(abs(squeeze(phi_mcx(:,31,:))')),clines,'linestyle','--','color',[0.9100    0.4100    0.1700],'linewidth',1.5);
colorbar('EastOutside');

%== plot the tissue boundaries
% plot a dashed circle with radius of 10
[xcirc,ycirc] = cylinder(10,200);
xcirc=xcirc(1,:)+30.5;
ycirc=ycirc(1,:)+30.5;
plot(xcirc,ycirc,'k--','linewidth',2);

% plot a dashed circle with radius of 23
[xcirc,ycirc] = cylinder(23,200);
xcirc=xcirc(1,:)+30.5;
ycirc=ycirc(1,:)+30.5;
plot(xcirc,ycirc,'k--','linewidth',2);

% plot a dashed circle with radius of 25
[xcirc,ycirc] = cylinder(25,200);
xcirc=xcirc(1,:)+30.5;
ycirc=ycirc(1,:)+30.5;
plot(xcirc,ycirc,'k--','linewidth',2);

axis equal;
colormap;
lg=legend('MMC','DMMC','MCX','Location','northoutside','orientation','horizontal');
set(lg,'Color',[0.5 0.5 0.5]);
set(lg,'box','on');
set(gca,'fontsize',18);