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 between the different tissues.

% this script requires the triplot function, which can
% be found in the private directory of FieldTrip (<directory of Fieldtrip>/private). The ft_prepare_bemmodel function is using functions that are also 
% located in the private directory. Furthermore, the script requires
% the image processing toolbox for the segmentations, and the SPM2 toolbox
% for reading the example anatomical data. 
% To read the mri file correctly, as a windows user, I had to substitute some dll files
% of SPM2 locating in external directory of fieldtrip toolbox by those updated dll files in
% ftp://ftp.fil.ion.ucl.ac.uk/spm/spm2_updates/ 

% 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 = [];
cfg.coordinates='spm';
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. By flipping, the order of the dimensions can be 
% different for different mri data. If the orientation has to be changed with 90 degree, the permute function can
% be used as well.
% 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);

% check whether the "gray", "white and "csf" segmentation worked well
seg.anatomy = mri.anatomy;
cfg = [];
cfg.funparameter = 'gray';
% or
% cfg.funparameter = 'white';
% cfg.funparameter = 'csf';
cfg.interactive = 'yes';
ft_sourceplot(cfg, seg);

% the construction of the segmentations uses the image processing toolbox
% basically this requires a lot of trial-and-error, with ft_sourceplot in between

% construct a segmentation of the brain compartment. It should be whole 3D structure without 
% holes. Hemispheres should be united.
brain = (seg.gray>0.5 | seg.white>0.5);
% Check if this segmentation includes the brain only. And afterwards, fill up
% the holes.
s = strel_bol(5);
brain = imclose(brain, s);
% brain = imopen(brain, s);
brain = imdilate(brain,strel_bol(2));
brain = imfill(brain, 'holes');

% Meanwhile plot the volume. This shows only the brain compartment, and no other parts of the anatomical image.
cfg = [];
cfg.interactive = 'yes';
seg2 = seg;
seg2.anatomy = brain.*seg.anatomy;
% or later:
% seg2.anatomy = skin.*seg.anatomy;
% seg2.anatomy = skull.*seg.anatomy;
% seg2.anatomy = (skin+skull+brain);
figure;
ft_sourceplot(cfg, seg2);

% Sometimes it is difficult to see if the dark parts are holes in the image or not. To see if
% there are remaining holes, you also can use (this will give a black and white picture of the volume):
cfg = [];
cfg.interactive = 'yes';
seg2 = seg;
seg2.anatomy = brain;
% or
% seg2.anatomy = skin;
% seg2.anatomy = skull;
figure;
ft_sourceplot(cfg, seg2);

% You also can check the brain compartment on the top of the entire anatomical data. Use the brain as a mask. The brain compartment will show up on a red scale 
%(you can change the color, see at help ft_sourceplot).
seg2.brain=brain;
seg2.anatomy = mri.anatomy;
% this is needed, for showing the anatomical data outside of mask (=brain)
seg2.gray=mri.anatomy; 
% mask parameter of ft_sourceplot works only if the functional data is also defined because brain will mask the functional data.  Parts of functional data where brain = 0 will be 
% transparent (so the original anatomical data will show up) and the brain compartment (where brain = 1) will show up with a red scaled anatomical structure.

cfg = [];
cfg.interactive = 'yes';
cfg.funparameter = 'gray';
cfg.maskparameter='brain';
figure;
ft_sourceplot(cfg,seg2);


% construct a segmentation of the skin compartment.  
skin = (mri.anatomy>300);
% Use an enough low value that the resulting
% skin volume could consist everything below the head surface, but nothing outside. 
% Check indices of vertices that skin should contain and indices that skin should not with
% ft_sourceplot. (The command line shows the indices when you click on the slices). And get the 
% values from mri.anatomy(i,j,k).


% this mri shows a bar at the top of head. it should be removed. First, locate the bar with
% ft_sourceplot. It is in mri.anatomy(:,:,177:181) where there isn't any part of the head.

skin = bwlabeln(skin); 
unique(skin(:,:,178:181)) %check the non-zero elements
% change each non-zero element to 0
skin(skin==159) = 0;
skin(skin==160) = 0;
skin(skin==161) = 0;
skin(skin==162) = 0;
skin(skin==163) = 0;
skin(skin==164) = 0;
skin = (skin~=0); %skin is logical (type) again

% If it is necessary, close the ear holes. Use ft_sourceplot to locate the earholes.
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;
skin = imclose(skin, s);
skin(:,:,1) = 1; 
skin = imfill(skin, 'holes');
skin(:,:,1) = 0;

% When your mri image is noisy (which is not the case now), it is possible that the skin volume
% will contain some 1's also outside of the actual skin,and they will show up as
% white points in the black outside space.
% This may help on it:
% skin=imopen(skin, strel_bol(5));


% construct a segmentation of the skull compartment
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
cfg = [];
cfg.interactive = 'yes';
seg2 = seg;
seg2.anatomy = brain.*seg.anatomy;
% or
% seg2.anatomy = skin.*seg.anatomy;
% seg2.anatomy = skull.*seg.anatomy;
% seg2.anatomy = (skin+skull+brain);
figure;
ft_sourceplot(cfg, seg2);

% 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. Many of the functions used by ft_prepare bemmodel are located in the private directory. 
%(switch your current directory to the private directory)

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 = true;
cfg.method         = 'dipoli';
%due to a bug you also have to specify both units
cfg.sourceunits = 'mm';
cfg.mriunits = 'mm';
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. (If the 3D picture is
% full of shadows, use Insert... Light from the menu.)
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