%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 2b: 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
tissues = {'FAT','FAT2','GM','MARROW','SKIN-MUSCLES','SKULL','VESSELS',...
    'WM','CSF','MUSCLES','DURA'};
%1. Fat 2. Fat2 3. GM 4. Marrow 5. Skin-muscles 6. Skull 7. Vessels
%8. WM 9. CSF 10. Muscles 11. Dura
for i = 1:length(tissues)
    data(i) = loadnifti(sprintf('datasets/2b/%s.nii.gz',tissues{i}));
    data(i).NIFTIData = mat2gray(data(i).NIFTIData);
end

%% Simplified to a 5-layer model
seg.wm = data(8).NIFTIData; %WM
seg.gm = data(3).NIFTIData; %GM
seg.csf = data(9).NIFTIData  + data(7).NIFTIData + data(11).NIFTIData; % CSF + Vessels + Dura
seg.skull = data(6).NIFTIData + data(4).NIFTIData; % Skull + Marrow
seg.scalp = data(5).NIFTIData + data(10).NIFTIData + data(1).NIFTIData + data(2).NIFTIData; %Skin-muscles + muscles + fat/fat2

%== 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 = 3.4;
cfg.samplinggm = 3.4;
cfg.samplingcsf = 4;
cfg.samplingskull = 5;
cfg.samplingscalp = 6;

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

%== 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 = 100;

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

%% Plot volumetric mesh
subplot(1,2,1)
hold on
plotmesh(node,elem(elem(:,5)==5,:),'x>137','FaceColor',[1 1 1],'EdgeAlpha',0.85) %%wm
plotmesh(node,elem(elem(:,5)==4,:),'x>137','FaceColor',[0.35 0.35 0.35],'EdgeAlpha',0.85) %%pial
plotmesh(node,elem(elem(:,5)==3,:),'x>137','FaceColor',[0.2 0.6 1],'EdgeAlpha',0.85) %%csf
plotmesh(node,elem(elem(:,5)==2,:),'x>137','FaceColor',[1 1 0.9],'EdgeAlpha',0.85) %%bone
plotmesh(node,elem(elem(:,5)==1,:),'x>137','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>185','FaceColor',[1 1 1],'EdgeAlpha',0.85) %%wm
plotmesh(node,elem(elem(:,5)==4,:),'y>185','FaceColor',[0.35 0.35 0.35],'EdgeAlpha',0.85) %%pial
plotmesh(node,elem(elem(:,5)==3,:),'y>185','FaceColor',[0.2 0.6 1],'EdgeAlpha',0.85) %%csf
plotmesh(node,elem(elem(:,5)==2,:),'y>185','FaceColor',[1 1 0.9],'EdgeAlpha',0.85) %%bone
plotmesh(node,elem(elem(:,5)==1,:),'y>185','FaceColor',[1 0.8 0.7],'EdgeAlpha',0.85) %%scalp
view([0 0])
axis off

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