Create BEM headmodel for EEG

This script shows how to make a volume conduction model of the head using the boundary element method (BEM). It starts by reading an anatomical MRI, segments the anatomical MRI and then computes the BEM model using the boundaries betwee the different tissues.

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

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

% read the anatomical MRI: this template MRI is available from 
% http://www.bic.mni.mcgill.ca/brainweb/selection_normal.html
% note that this MRI is aligned with MNI/SPM coordinates, if you have 
% a custom MRI you should ensure that it is aligned with the desired 
% coordinates using the homogeneous transformation matrix
mri = ft_read_mri('t1_icbm_normal_1mm_pn0_rf0.mnc');

% 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 = [];
seg = ft_volumesegment(cfg, mri);

% due to a bug in **[[reference:ft_volumesegment|ft_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

% construct a segmentation of the brain compartment
brain = (seg.gray>0.5 | seg.white>0.5);
s = strel_bol(5);
brain = imclose(brain, s);
% brain = imopen(brain, s);
brain = imdilate(brain, strel_bol(2));
brain = imfill(brain, 'holes');

% construct a segmentation of the skin compartment
skin = (mri.anatomy>300);
% remove the bar at top of head
skin = bwlabeln(skin);
skin(skin==159) = 0;
skin(skin==160) = 0;
skin(skin==161) = 0;
skin = (skin~=0);
% 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;

% construct a segmentation of the skull compartment
s = strel_bol(5);
skin_a  = imerode(skin, s);
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 = '/home/coherence/roboos/src/thom/dipoli'; % where to find the executable
vol = ft_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
% ft_sourceanalysis or ft_dipolefitting as cfg.hdmfile='vol.mat'
save vol.mat vol mri