% this script requires the volplot and triplot functinos, which both can
% be found in the private directory of fieldtrip. Futhermore, it requires
% the image processing toolbox for the segmentations, and the SPM2 toolbox
% for reading the example anatomical data.

% the prepare_bemmodel function is available upon request
% the dipoli stand-alone executable is available upon request

ResultFileName = 'BEMexample01';

% read the individual anatomical MRI
% should be isotropic and have 1*1*1 mm voxels

mri = read_fcdc_mri('iso_example01.hdr')

% check with
sourceplot([],mri);
% Interchange and flip dimensions to match the MNI format
% X left -> right
% Y back -> front
% Z bottom -> top


% construct a segmentation of the brain, i.e. gray+white+csf
% these cannot be directly used for the BEM model, but serve as starting point
cfg = [];
cfg.smooth = 5;
cfg.name = [ResultFileName,'_temp'];
seg = volumesegment(cfg, mri); % answer 'y', 'a', if aligned with MNI coordinates

% due to a bug in volumesegment, the segmentations should be flipped
% in case of an CTF-oriented MRI, but not in case of an SPM-oriented one
% seg.gray  = flipdim(flipdim(flipdim(seg.gray , 3), 2) ,1);
% seg.white = flipdim(flipdim(flipdim(seg.white, 3), 2) ,1);
% seg.csf   = flipdim(flipdim(flipdim(seg.csf  , 3), 2) ,1);

% the construction of the segmentations uses the image processing toolbox
% basically this requires a lot of trial-and-error, with "volplots" in between
% volplot(seg.gray) etc.

% construct a segmentation of the brain compartment
brain = (seg.gray>0.5 | seg.white>0.5  | seg.csf>0.7);
s = strel_bol(5); % 3D sphere with specified radius for 3D image processing
brain = imclose(brain, s);
%brain = imopen(brain, s);
brain = imdilate(brain, strel_bol(2));
brain = imfill(brain, 'holes');

% check for different clusters
brain = bwlabeln(brain,26);
% unique(brain)
% volplot(brain)
brain(brain==2) = 0;
brain = brain~=0;

% construct a segmentation of the skin compartment

% estimate the maximal value for the space around the head
temp = mri.anatomy([1:3,end-2:end],:,:);
border = temp(:);
temp = mri.anatomy(:,[1:3,end-2:end],:);
border = [border;temp(:)];
temp = mri.anatomy(:,:,[1:3,end-2:end]); % check change
border = [border;temp(:)];
 thresEstimate = prctile(border,98.5);
%thresEstimate = max(border)*1.5;

skin = (mri.anatomy>thresEstimate);
[skin,num] = bwlabeln(skin,26); % eliminates "selected" pixel that are in a cluster > 26 connected pixels

%
maxVal=max(skin(:));
tic
for i=1:maxVal
    fprintf('%i/%i\n',i,maxVal)
    len(i)=length(find(skin==i));
end
toc
[val,index]=sort(-len); val=-val;
index=index(1)
skin = (skin==index);

% close the ear holes
% skin(1:30,70:110,1:50) = mri.anatomy(1:30,70:110,1:50)>30;
% skin(150:181,70:110,1:50) = mri.anatomy(150:181,70:110,1:50)>30;

s = strel_bol(5);
skin = imclose(skin, s);
skin(:,:,1) = 1;
skin = imfill(skin, 'holes');
skin(:,:,1) = 0;


% Cut the neck
nCut = 60; % # of voxel to cut from the bottom
skin = skin(:,:,(nCut+1):end);
brain = brain(:,:,(nCut+1):end);
mri.anatomy=mri.anatomy(:,:,(nCut+1):end);
mri.zgrid = mri.zgrid(1:(end-nCut));
mri.transform(3,4) = mri.transform(3,4) + nCut*abs(mri.transform(3,3));
mri.antcomm(3) = mri.antcomm(3) - nCut;
mri.right(3) = mri.right(3) - nCut;
mri.left(3) = mri.left(3) - nCut;
mri.nasion(3) = mri.nasion(3) - nCut;
mri.nCut = nCut;
mri.dim = [size(mri.anatomy,1), size(mri.anatomy,2), size(mri.anatomy,3)];


% construct a segmentation of the skull compartment
s = strel_bol(4); % try 4-6
skin_a  = imerode(skin, s);
s = strel_bol(6); % try 5-7
brain_a = imdilate(brain, s);
skull = (brain_a & skin_a);

% make figures of the segmentations, click around in the figures
volplot(skin  .* mri.anatomy)
volplot(skull .* mri.anatomy)
volplot(brain .* mri.anatomy)
volplot(skin+skull+brain);

% add the BEM segmentation to the anatomical MRI for convenience
% skin = 1, skull = 2, brain = 3
mri.seg = skin+skull+brain;

% construct the triangulated surfaces and compute the BEM model
cfg                = [];
cfg.tissue         = [1 2 3]; % value of each tissue type in the segmentation
cfg.numvertices    = [1000 2000 3000];
cfg.conductivity   = [1 1/80 1];
cfg.isolatedsource = 3;
cfg.method         = 'dipoli';
cfg.dipoli = GetDipoliPath; % e.g.'/projects/EEGBeamformer/eegbem/thom/dipoli'; % where to find the executable
cfg.hdmfile = 'mri/example01/temp_hdm'; % where to write the temporary hdm files
vol = prepare_bemmodel(cfg, mri);

% make figures of the surfaces, note that the brain and skull surface
% should be slightly smoother an additional convolution and threshold of
% their respective segmentations probably would achieve that
figure; triplot(vol.bnd(1).pnt, vol.bnd(1).tri, [], 'faces_skin'); rotate3d
figure; triplot(vol.bnd(2).pnt, vol.bnd(2).tri, [], 'faces_skin'); rotate3d
figure; triplot(vol.bnd(3).pnt, vol.bnd(3).tri, [], 'faces_skin'); rotate3d

% write the BEM model to a matlab file, it can be later specified in
% sourceanalysis or dipolefitting as cfg.hdmfile='vol.mat'
save(sprintf('%s.mat',ResultFileName),'vol','mri')