% This is example code to show how to implement the skull extraction
% library
% Zachary Ip, University of Washington, ipzach@uw.edu, 2019.11.26
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, close all, clc
% Load data
MRI = squeeze(dicomread('T1MPRAGE'));
MRI_info = dicominfo('IM_0004');

% Ensure data is in 3D Matrix
disp(size(MRI));

%% Check which dimensions correspond to coronal, transverse, and saggital
% planes
figure
set(gcf,'Position',[0 0 2000 600])
subplot(1,3,1)
pcolor(flipud(squeeze(MRI(100,:,:)))), shading flat, colormap(bone),axis off
title('1st Dimension')
subplot(1,3,2)
pcolor(flipud(squeeze(MRI(:,200,:)))), shading flat, colormap(bone),axis off
title('2nd Dimension')
subplot(1,3,3)
pcolor(flipud(squeeze(MRI(:,:,40)))), shading flat, colormap(bone),axis off
title('3rd Dimension')

% Dimension 1: Transverse
% Dimension 2: Saggital
% Dimension 3: Coronal

%% Create rudimentary mask
% callibrate mask in your desired dimension, Coronal is reccomended.
[MRI_mask] = CalibrateMask(MRI,3);

 
%% Fill locations for each frame
% Remove 'musculature' region by filling in 'outside' region
% Fill both coronal and transverse sections. Saggital not reccomended.
coronalFill = FillExterior(MRI_mask, 3, false);
transverseFill = FillExterior(MRI_mask, 1, false);

disp('mark')
%% Invert Mask, remove noise
coronalFill = ~coronalFill;
transverseFill = ~transverseFill;
 
coronalClean = RemoveNoise(coronalFill);
transverseClean = RemoveNoise(transverseFill);
 
cleanMask = MergeMasks(coronalClean,transverseClean);

%% View cleaned mask
for i = 1:size(cleanMask,1)
pcolor(flipud(squeeze(cleanMask(i,:,:))))
shading flat, colormap(bone);title(num2str(i))
drawnow
end
%% Scale Image
% Get pixel size info and frame step distance from DICOM header
xsize = MRI_info.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.PixelSpacing(1);
ysize = MRI_info.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.PixelSpacing(2);
zsize = MRI_info.SpacingBetweenSlices;

% Define a scaler to adjust for difference in resolution
scaler = zsize/xsize;

% Interpolate points in the lower resolution dimension to create a uniform
% resolution for printing
[scaledMask, X,Y,Z] = ScaleMask(cleanMask, scaler, 3);
disp('mark')
%% Visualize 3D model
figure;
vis3d = patch(isosurface(scaledMask));
vis3d.FaceColor = 'w';
vis3d.EdgeColor = 'none';
daspect([1 1 1]);
camlight; 
lighting phong
 
%% Locate Craniotomy
% Use anatomical features from MRI to locate Craniotomy origin
% xC will denote x-center-of-craniotomy, and so forth
[xC,yC,zC] = LocateCraniotomy(MRI);

%% Perform virtual craniotomy
% We will make a craniotomy with radius 12.5 mm. xsize will be put in to
% know the resolution of the voxel
[craniotomyMask] = Craniotomy(scaledMask, xC, yC, zC, 12.5, X, Y, Z, xsize);

%% Visualize 3D model
figure;
vis3d = patch(isosurface(craniotomyMask));
vis3d.FaceColor = 'w';
vis3d.EdgeColor = 'none';
daspect([1 1 1]);
camlight; 
lighting phong
%% Crop models
% To reduce file size, printing time, and create a cleaner model, we crop
% our models

scaledMask = scaledMask(1:170,:,:);
craniotomyMask = craniotomyMask(1:170,:,:);

%% Save file as stl
% Make sure you export with a function that encodes resolution.
% This code was taken from https://www.mathworks.com/matlabcentral/fileexchange/68794-make-stl-of-3d-array-optimal-for-3d-printing
make_STL_of_Array('exampleSkull.stl',scaledMask,xsize,xsize,xsize);
make_STL_of_Array('exampleCraniotomySkull.stl',craniotomyMask,xsize,xsize,xsize);
disp('File saved!')