%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo 1: SPM to Mesh using Cgalmesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%% Extracting the tissue probability maps for each tissue
filename = 'ANTS17-0Years_head.nii.gz';

%== Reading in dataand storing into 4D array
for i = 1:6
    A= loadnifti(sprintf('datasets/1/c%i%s',i,filename));
    dim = size(A.NIFTIData);
    seg(:,:,:,i) = A.NIFTIData;
end

%== Find indices with max probability
[~,I] = max(seg,[],4);
I(:,:,1:70,:) =0;

%== Removing background and non-brain related areas
I(I==6)=0;

%== max radius of the Delaunay sphere
opt.radbound = 2;

%== maximum deviation from the specified isosurfaces
opt.distbound = 2;

%% Run automated image-to-mesh conversion using cgalmesh method with v2m
[node,elem,face] = v2m(uint8(I),[],opt,50,'cgalmesh');
%--------------------------------------------------------------------------
%%% Hands-on exercise: change maxvol and rerun, compare the mesh densities
%--------------------------------------------------------------------------

%% mesh visualization
figure(1);
subplot(2,2,1);plotmesh(node(:,1:3),elem);

%== Plot specific color and transparency and label
subplot(2,2,2);plotmesh(node(:,1:3),elem(elem(:,5)==2,:),'FaceColor',[0.2 0.6 1],'EdgeAlpha',0.85,'FaceAlpha',0.3)

%== Limit plotting to a certain region
subplot(2,2,3);plotmesh(node(:,1:3),elem,'x>78 | y<110 | z>100')
 
%== Combine all parameters:
subplot(2,2,4);plotmesh(node(:,1:3),elem(elem(:,5)==2,:),'x>78| y<110','FaceColor',[0.35 0.35 0.35],'EdgeAlpha',0.85)

%== To plot specific element labels, e.g. labels 1 and 2:
figure(2);plotmesh(node(:,1:3),elem(elem(:,5)==1 | elem(:,5)==2,:))
