%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 2a: Surface-based meshing using SPM data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%% load SPM segmentations for all tissue types (based on a 5 layer model)
data = loadnifti('datasets/1/c2ANTS17-0Years_head.nii.gz');
seg.wm = data.NIFTIData;     % white matter
data = loadnifti('datasets/1/c1ANTS17-0Years_head.nii.gz');
seg.gm = data.NIFTIData;     % gray matter
data = loadnifti('datasets/1/c3ANTS17-0Years_head.nii.gz');
seg.csf = data.NIFTIData;    % CSF
data = loadnifti('datasets/1/c4ANTS17-0Years_head.nii.gz');
seg.skull = data.NIFTIData;  % skull
data = loadnifti('datasets/1/c5ANTS17-0Years_head.nii.gz');
seg.scalp = data.NIFTIData;  % scalp

%== Radius of the Delaunay sphere used in the sampling the surfaces.
% Default values are 1.7, 1.7, 2, 2.5 and 3, respectively (reference values for 1x1x1mm^3)
% Scale proportionally for denser volumes. Lower values correspond to denser, higher
% fidelity surface extraction, but also results in denser meshes.
cfg.samplingwm = 2;
cfg.samplinggm = 2;
cfg.samplingcsf = 2.5;
cfg.samplingskull = 3.5;
cfg.samplingscalp = 3.5;
%--------------------------------------------------------------------------
%%% Hands-on exercise: change Radius of the Delaunay sphere and see how
%%% this impact mesh densities
%--------------------------------------------------------------------------

%== mesh sizing field used at the mesh generation step.
% Higher values can help lowering the final mesh density.
cfg.meshfield = 5;

%== indicates the volumetric maximum size of elements
% Lowering this value helps with obtaining a denser tetrahedral
% mesh. For dense meshes, values close to 3-5 are recommended.
cfg.maxvol = 50;

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

%% Plot volumetric mesh
figure(1);
subplot(1,2,1)
hold on
plotmesh(node,elem(elem(:,5)==5,:),'x>67','FaceColor',[1 1 1],'EdgeAlpha',0.85) %%wm
plotmesh(node,elem(elem(:,5)==4,:),'x>67','FaceColor',[0.35 0.35 0.35],'EdgeAlpha',0.85) %%pial
plotmesh(node,elem(elem(:,5)==3,:),'x>67','FaceColor',[0.2 0.6 1],'EdgeAlpha',0.85) %%csf
plotmesh(node,elem(elem(:,5)==2,:),'x>67','FaceColor',[1 1 0.9],'EdgeAlpha',0.85) %%bone
plotmesh(node,elem(elem(:,5)==1,:),'x>67','FaceColor',[1 0.8 0.7],'EdgeAlpha',0.85) %%scalp
view([270 0])
axis off
subplot(1,2,2)
hold on
plotmesh(node,elem(elem(:,5)==5,:),'y>105','FaceColor',[1 1 1],'EdgeAlpha',0.85) %%wm
plotmesh(node,elem(elem(:,5)==4,:),'y>105','FaceColor',[0.35 0.35 0.35],'EdgeAlpha',0.85) %%pial
plotmesh(node,elem(elem(:,5)==3,:),'y>105','FaceColor',[0.2 0.6 1],'EdgeAlpha',0.85) %%csf
plotmesh(node,elem(elem(:,5)==2,:),'y>105','FaceColor',[1 1 0.9],'EdgeAlpha',0.85) %%bone
plotmesh(node,elem(elem(:,5)==1,:),'y>105','FaceColor',[1 0.8 0.7],'EdgeAlpha',0.85) %%scalp
view([0 0])
axis off

figure(2);
subplot(1,2,1)
plotmesh(node,elem(elem(:,5)==5,:),'FaceColor',[1 1 1],'EdgeAlpha',0.85)
axis off
subplot(1,2,2)
plotmesh(node,elem(elem(:,5)==4,:),'FaceColor',[0.35 0.35 0.35])
axis off