%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 3: MMC simulation of Brain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%% Reading in data from c1-c5 SPM segmentations using a 4D array
for i = 1:5
    A = loadnifti(sprintf('datasets/3/c%iANTS19-5Years_head.nii.gz',i));
    dim = size(A.NIFTIData);
    data(:,:,:,i) = A.NIFTIData;
end

%% load SPM segmentations for all tissue types (based on a 5 layer model)
seg.wm = data(:,:,:,2); 
seg.gm = data(:,:,:,1);
seg.csf = data(:,:,:,3); 
seg.skull = data(:,:,:,4); 
seg.scalp = data(:,:,:,5); 

%% Run brain2mesh to generate volumetric mesh for the SPM segmentation
[node,elem,face] = brain2mesh(seg);

%% prepare cfg for MMCLAB simulation
clear cfg

%== Define simulated photon count
cfg.nphoton=3e6; %takes ~120s on Intel i7-8700K(12 threads)

%== Define the mesh domain for MMC simulation
cfg.node=node;
    cfg.elem=elem;
    
%== medium labels:0-ambient air,1-air cavities,2-scalp,3-skull,4-csf,5-gray matter,6-white matter 
cfg.elem(:,5)=cfg.elem(:,5)+1;

%== tetrahedral element label which represents different tissue types
cfg.elemprop=cfg.elem(:,5);

%== 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          0      1      1;       %<- tissue label 1: air cavity
            0.019      7.8    0.89   1.37;     %<- tissue label 2: scalp
            0.019      7.8    0.89   1.37;     %<- tissue label 3: skull
            0.004      0.009  0.89   1.37;     %<- tissue label 4: CSF
            0.02       9.0    0.89   1.37;     %<- tissue label 5: gray matter
            0.08      40.9    0.84   1.37];    %<- tissue label 6: white matter

%== Define light source: pencil beam
cfg.srctype='pencil';

%== Pencil beam source placed at EEG 10-5 landmark:"C4h", calculated using Mesh2EEG
cfg.srcpos=[133.5370,90.1988,200.0700];

%== Source direction along node norm, calculated using nodesurfnorm (https://github.com/fangq/metch/blob/master/nodesurfnorm.m)
cfg.srcdir=[-0.5086,-0.1822,-0.8415]; 
    cfg.srcpos=cfg.srcpos+0.001*cfg.srcdir; % ensure source is inside the mesh

%% Define time-gate settings
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 other simulation parameters

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

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

%== dual-grid mode
cfg.method='grid'; 

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

%% Run mmc simulation
[flux]=mmclab(cfg);

%% post-simulation data processing and visualization
% convert time-resolved fluence to CW fluence
CWfluence=sum(flux.data*cfg.tstep,4); 

% coronal plane selected for fluence plot
y_plane=90.5;   

% fix spatial offset between the coordinate origins of DMMC AABB and the mesh
DMMC_AABB_origin=min(cfg.node); %spatial offset of the DMMC AABB
DMMC_voxel_origin=DMMC_AABB_origin+0.5; %centroid of the 1st voxel for DMMC AABB
DMMC_AABB_dimension=size(CWfluence);
CWfluence_DMMC=squeeze(CWfluence(:,ceil(y_plane-DMMC_AABB_origin(1,2)),:)); %along y-axis, DMMC AABB happens to have an offset very close to 3(2.9955)
[xx,zz]=meshgrid(DMMC_voxel_origin(1,1):(DMMC_voxel_origin(1,1)+DMMC_AABB_dimension(1,1)-1),DMMC_voxel_origin(1,3):(DMMC_voxel_origin(1,3)+DMMC_AABB_dimension(1,3)-1));

figure;
clines=-20:0.5:0;
contourf(xx,zz,log10(abs(CWfluence_DMMC')),clines,'linestyle','--','color','w','linewidth',1.5,'DisplayName','DMMC');
hold on;axis equal;
colorbar('EastOutside');

%== plot tissue boundaries, auxiliary functions are part of iso2mesh(http://iso2mesh.sf.net)
plane=[0 y_plane 0; 0 y_plane 1; 1 y_plane 0];
label=unique(cfg.elemprop);
face=[];
for i=1:2:size(label,1)
    newface=volface(cfg.elem(cfg.elemprop==label(i),1:4));
    face=[face;newface];
end
[cutpos,cutvalue,cutedges]=qmeshcut(face(:,1:3),node,node(:,1),plane);
[cutpos,cutedges]=removedupnodes(cutpos,cutedges);
cutloop=extractloops(cutedges);
[nanidx]=find(isnan(cutloop));

for i=1:size(nanidx,2)
    if(i==9) 
        continue;
    end
    if(i==1)
        plot(cutpos(cutloop(1:(nanidx(i)-1)),1),cutpos(cutloop(1:(nanidx(i)-1)),3),'linestyle','--','color','k','LineWidth',1.25,'HandleVisibility','off');
    elseif(i==2)
        plot(cutpos(cutloop((nanidx(i-1)+1):(nanidx(i)-1)),1),cutpos(cutloop((nanidx(i-1)+1):(nanidx(i)-1)),3),'linestyle','--','color','k','LineWidth',1.25,'HandleVisibility','off');
    else
        plot(cutpos(cutloop((nanidx(i-1)+1):(nanidx(i)-1)),1),cutpos(cutloop((nanidx(i-1)+1):(nanidx(i)-1)),3),'linestyle','--','color','k','LineWidth',1.25,'HandleVisibility','off');
    end
end

%== plot source, legend, etc.
plot(cfg.srcpos(1,1),cfg.srcpos(1,3),'o','MarkerEdgeColor','r','MarkerFaceColor','r','MarkerSize',10,'DisplayName','source');
lg=legend('Location','northeast');
set(lg,'color','[0.85 0.85 0.85]');
set(lg,'box','on');
set(gca,'ylim', [160 225]);ylabel('z [mm]');
set(gca,'xlim', [10 165]);xlabel('x [mm]');
set(gca,'clim',[-12 0]);
set(gca,'fontsize',18);
set(gca, 'FontName', 'Times New Roman');