Skip to content
Permalink
release
Switch branches/tags
Go to file
 
 
Cannot retrieve contributors at this time
function [headmodel, cfg] = ft_prepare_headmodel(cfg, data)
% FT_PREPARE_HEADMODEL constructs a volume conduction model from the geometry
% of the head. The volume conduction model specifies how currents that are
% generated by sources in the brain, e.g. dipoles, are propagated through the
% tissue and how these result in externally measureable EEG potentials or MEG
% fields.
%
% FieldTrip implements a variety of forward solutions, partially with internal
% code and some of them using external toolboxes or executables. Each of the
% forward solutions requires a set of configuration options which are listed
% below. This function takes care of all the preparatory steps in the
% construction of the volume conduction model and sets it up so that
% subsequent computations are efficient and fast.
%
% Use as
% headmodel = ft_prepare_headmodel(cfg) or
% headmodel = ft_prepare_headmodel(cfg, mesh) with the output of FT_PREPARE_MESH or FT_READ_HEADSHAPE
% headmodel = ft_prepare_headmodel(cfg, seg) with the output of FT_VOLUMESEGMENT
% headmodel = ft_prepare_headmodel(cfg, elec) with the output of FT_READ_SENS
% headmodel = ft_prepare_headmodel(cfg, sourcemodel) with the output of FT_PREPARE_LEADFIELD
%
% In general the input to this function is a geometrical description of the
% shape of the head and a description of the electrical conductivity. The
% geometrical description can be a set of surface points obtained from
% fT_READ_HEADSHAPE, a surface mesh that was obtained from FT_PREPARE_MESH or
% a segmented anatomical MRI that was obtained from FT_VOLUMESEGMENT.
%
% The cfg argument is a structure that can contain:
% cfg.method = string that specifies the forward solution, see below
% cfg.conductivity = a number or a vector containing the conductivities of the compartments
% cfg.tissue = a string or integer, to be used in combination with a 'seg' for the
% second intput. If 'brain', 'skull', and 'scalp' are fields
% present in 'seg', then cfg.tissue need not be specified, as
% these are defaults, depending on cfg.method. Otherwise,
% cfg.tissue should refer to which field(s) of seg should be used.
%
% For EEG the following methods are available:
% singlesphere analytical single sphere model
% concentricspheres analytical concentric sphere model with up to 4 spheres
% openmeeg boundary element method, based on the OpenMEEG software
% bemcp boundary element method, based on the implementation from Christophe Phillips
% dipoli boundary element method, based on the implementation from Thom Oostendorp
% asa boundary element method, based on the (commercial) ASA software
% simbio finite element method, based on the SimBio software
% duneuro finite element method, based on duneuro software
% fns finite difference method, based on the FNS software
% infinite electric dipole in an infinite homogenous medium
% halfspace infinite homogenous medium on one side, vacuum on the other
% besa finite element leadfield matrix from BESA
% interpolate interpolate the precomputed leadfield
%
% For MEG the following methods are available:
% openmeeg boundary element method, based on the OpenMEEG software
% singlesphere analytical single sphere model
% localspheres local spheres model for MEG, one sphere per channel
% singleshell realisically shaped single shell approximation, based on the implementation from Guido Nolte
% infinite magnetic dipole in an infinite vacuum
%
% Each specific method has its own specific configuration options which are listed below.
%
% BEMCP, DIPOLI, OPENMEEG
% cfg.tissue see above; in combination with 'seg' input
% cfg.isolatedsource (optional)
% cfg.tempdir (optional)
% cfg.tempname (optional)
%
% CONCENTRICSPHERES
% cfg.tissue see above; in combination with 'seg' input
% cfg.order (optional)
% cfg.fitind (optional)
%
% LOCALSPHERES
% cfg.grad
% cfg.tissue see above; in combination with 'seg' input; default options are 'brain' or 'scalp'
% cfg.feedback (optional)
% cfg.radius (optional)
% cfg.maxradius (optional)
% cfg.baseline (optional)
%
% SIMBIO
% cfg.conductivity
%
% DUNEURO
% cfg.conductivity An array with the conductivities must be provided. (see above)
% cfg.grid_filename Alternatively, a filename for the grid and a filename for the conductivities can be passed.
% cfg.tensors_filename "
% cfg.duneuro_settings (optional) Additional settings can be provided for duneuro (see http://www.duneuro.org).
%
% SINGLESHELL
% cfg.tissue see above; in combination with 'seg' input; default options are 'brain' or 'scalp'
% cfg.order (optional)
%
% SINGLESPHERE
% cfg.tissue see above; in combination with 'seg' input; default options are 'brain' or 'scalp'; must be only 1 value
%
% INTERPOLATE
% cfg.outputfile (required) string, filename prefix for the output files
%
% BESA
% cfg.headmodel (required) string, filename of precomputed FEM leadfield
% cfg.elec (required) structure with electrode positions or filename, see FT_READ_SENS
% cfg.outputfile (required) string, filename prefix for the output files
%
% FNS
% cfg.tissue
% cfg.tissueval
% cfg.conductivity
% cfg.elec
% cfg.grad
% cfg.transform
%
% HALFSPACE
% cfg.point
% cfg.submethod (optional)
%
% More details for each of the specific methods can be found in the corresponding
% low-level function which is called FT_HEADMODEL_XXX where XXX is the method
% of choise.
%
% See also FT_PREPARE_MESH, FT_PREPARE_SOURCEMODEL, FT_PREPARE_LEADFIELD,
% FT_HEADMODEL_BEMCP, FT_HEADMODEL_ASA, FT_HEADMODEL_DIPOLI,
% FT_HEADMODEL_SIMBIO, FT_HEADMODEL_FNS, FT_HEADMODEL_HALFSPACE,
% FT_HEADMODEL_INFINITE, FT_HEADMODEL_OPENMEEG, FT_HEADMODEL_SINGLESPHERE,
% FT_HEADMODEL_CONCENTRICSPHERES, FT_HEADMODEL_LOCALSPHERES,
% FT_HEADMODEL_SINGLESHELL, FT_HEADMODEL_INTERPOLATE, FT_HEADMODEL_DUNEURO
% Copyright (C) 2011, Cristiano Micheli
% Copyright (C) 2011-2012, Jan-Mathijs Schoffelen, Robert Oostenveld
% Copyright (C) 2013-2018, Robert Oostenveld, Johanna Zumer
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
% FieldTrip is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% FieldTrip is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$
% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin = nargin;
ft_nargout = nargout;
% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble trackconfig
ft_preamble provenance data
% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
return
end
% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'deprecated', 'geom');
cfg = ft_checkconfig(cfg, 'renamed', {'geom', 'headshape'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'bem_openmeeg', 'openmeeg'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'bem_dipoli', 'dipoli'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'bem_cp', 'bemcp'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'nolte', 'singleshell'});
cfg = ft_checkconfig(cfg, 'renamed', {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'vol', 'headmodel'});
if isfield(cfg, 'headmodel') && ischar(cfg.headmodel) && ~endsWith(cfg.headmodel, '.lft')
cfg.method = 'file'; % FIXME this is not documented, note that it does not apply to BESA headmodels
elseif isfield(cfg, 'headmodel') && isstruct(cfg.headmodel)
cfg.method = 'existing'; % FIXME this is not documented
end
cfg = ft_checkconfig(cfg, 'required', 'method');
% set the general defaults
cfg.headshape = ft_getopt(cfg, 'headshape');
cfg.conductivity = ft_getopt(cfg, 'conductivity');
% volume related options
cfg.tissue = ft_getopt(cfg, 'tissue');
cfg.smooth = ft_getopt(cfg, 'smooth');
cfg.threshold = ft_getopt(cfg, 'threshold');
% other options
cfg.isolatedsource = ft_getopt(cfg, 'isolatedsource'); % used for dipoli and openmeeg
cfg.tempdir = ft_getopt(cfg, 'tempdir'); % used for dipoli
cfg.tempname = ft_getopt(cfg, 'tempname'); % used for dipoli
cfg.point = ft_getopt(cfg, 'point'); % used for halfspace
cfg.submethod = ft_getopt(cfg, 'submethod'); % used for halfspace
cfg.feedback = ft_getopt(cfg, 'feedback');
cfg.radius = ft_getopt(cfg, 'radius');
cfg.maxradius = ft_getopt(cfg, 'maxradius');
cfg.baseline = ft_getopt(cfg, 'baseline');
cfg.singlesphere = ft_getopt(cfg, 'singlesphere');
cfg.grid_filename = ft_getopt(cfg, 'grid_filename'); % used for duneuro
cfg.tensors_filename= ft_getopt(cfg, 'tensors_filename'); % used for duneuro
cfg.duneuro_settings= ft_getopt(cfg, 'duneuro_settings');
cfg.tissueval = ft_getopt(cfg, 'tissueval'); % used for simbio
cfg.transform = ft_getopt(cfg, 'transform');
cfg.siunits = ft_getopt(cfg, 'siunits', 'no'); % yes/no, convert the input and continue with SI units
cfg.unit = ft_getopt(cfg, 'unit');
cfg.smooth = ft_getopt(cfg, 'smooth'); % used for interpolate
cfg.headmodel = ft_getopt(cfg, 'headmodel'); % can contain CTF localspheres model
% the data can be passed as input arguments or can be read from disk
hasdata = exist('data', 'var');
if hasdata
% the data should describe the geometrical mesh
if isfield(data, 'bnd')
data = data.bnd;
end
% check if the input data is valid for this function and ensure that it has the units specified
data = ft_checkdata(data, 'hasunit', 'yes');
% replace pnt by pos
data = fixpos(data);
else
data = [];
end
% convert to SI units
if istrue(cfg.siunits)
cfg.unit = 'm';
end
% convert the geometrical data to the desired units for the source model
if ~isempty(cfg.unit)
if ~isempty(data)
data = ft_convert_units(data, cfg.unit);
end
if isfield(cfg, 'grad') && ~isempty(cfg.grad)
cfg.grad = ft_convert_units(cfg.grad, cfg.unit);
end
if isfield(cfg, 'elec') && ~isempty(cfg.elec)
cfg.elec = ft_convert_units(cfg.elec, cfg.unit);
end
end
% the conductivity in the data overrules cfg.conductivity
if hasdata && isfield(data, 'cond')
cfg.conductivity = data.cond;
end
% boolean variables to manages the different geometrical input data objects
input_mesh = ft_datatype(data, 'mesh');
input_seg = ft_datatype(data, 'segmentation');
input_elec = ft_datatype(data, 'sens');
input_pos = ~input_mesh && isfield(data, 'pos'); % surface points without triangulation
% the construction of the volume conductor model is performed below
switch cfg.method
case 'file'
% read it from file
headmodel = ft_read_headmodel(cfg.headmodel);
cfg = rmfield(cfg, 'method'); % FIXME this is not documented
case 'existing'
% return an existing one
headmodel = cfg.headmodel;
cfg = rmfield(cfg, 'method'); % FIXME this is not documented
case 'interpolate'
% the "data" here represents the output of FT_PREPARE_LEADFIELD, i.e. a regular dipole
% grid with pre-computed leadfields
sens = ft_fetch_sens(cfg, data);
headmodel = ft_headmodel_interpolate(cfg.outputfile, sens, data, 'smooth', cfg.smooth);
case 'besa'
% cfg.headmodel points to the filename of the FEM solution that was computed in BESA
% cfg.elec points to the filename of the corresponding electrode specification
sens = ft_fetch_sens(cfg, data);
headmodel = ft_headmodel_interpolate(cfg.outputfile, sens, cfg.headmodel, 'smooth', cfg.smooth);
case 'asa'
if ~ft_filetype(cfg.headmodel, 'asa_vol')
ft_error('You must supply a valid cfg.headmodel for use with ASA headmodel')
end
headmodel = ft_headmodel_asa(cfg.headmodel);
case {'bemcp' 'dipoli' 'openmeeg'}
% the low-level functions all need a mesh
if isfield(data, 'pos') && isfield(data, 'tri')
if ~isfield(cfg, 'numvertices') || isempty(cfg.numvertices) || isequal(cfg.numvertices, arrayfun(@(x) size(x.pos, 1), data))
% copy the input data
geometry = data;
else
% retriangulate the input to the user-specified number of vertices
tmpcfg.method = 'headshape';
tmpcfg.headshape = data;
tmpcfg.numvertices = cfg.numvertices;
geometry = ft_prepare_mesh(tmpcfg);
end
elseif isfield(data, 'transform') && isfield(data, 'dim')
tmpcfg = keepfields(cfg, {'numvertices', 'tissue', 'spmversion'});
geometry = ft_prepare_mesh(tmpcfg, data);
else
ft_error('Either a segmented MRI or data with closed triangulated mesh is required as data input for the bemcp, dipoli or openmeeg method');
end
if strcmp(cfg.method, 'bemcp')
headmodel = ft_headmodel_bemcp(geometry, 'conductivity', cfg.conductivity);
if any(isnan(headmodel.mat(:)))
% HACK add a little bit of noise, with the NatMEG tutorial data, I discovered that this prevents the warning
% Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
geometry(1).pos = geometry(1).pos + randn(size(geometry(1).pos))*ft_scalingfactor('um', geometry(1).unit);
geometry(2).pos = geometry(2).pos + randn(size(geometry(2).pos))*ft_scalingfactor('um', geometry(2).unit);
geometry(3).pos = geometry(3).pos + randn(size(geometry(3).pos))*ft_scalingfactor('um', geometry(3).unit);
ft_warning('NaN detected, trying once more with slightly different vertex positions');
headmodel = ft_headmodel_bemcp(geometry, 'conductivity', cfg.conductivity);
end
elseif strcmp(cfg.method, 'dipoli')
headmodel = ft_headmodel_dipoli(geometry, 'conductivity', cfg.conductivity, 'isolatedsource', cfg.isolatedsource, 'tempdir', cfg.tempdir, 'tempname', cfg.tempname);
else
headmodel = ft_headmodel_openmeeg(geometry, 'conductivity', cfg.conductivity, 'isolatedsource', cfg.isolatedsource, 'tissue', cfg.tissue);
end
case 'concentricspheres'
cfg.fitind = ft_getopt(cfg, 'fitind');
cfg.order = ft_getopt(cfg, 'order');
% the low-level functions needs surface points, triangles are not needed
if input_mesh || input_pos
geometry = data;
elseif input_seg
tmpcfg = keepfields(cfg, {'numvertices', 'tissue', 'spmversion'});
geometry = ft_prepare_mesh(tmpcfg, data);
elseif input_elec
geometry.pos = data.chanpos;
geometry.unit = data.unit;
elseif ~isempty(cfg.headshape) && ischar(cfg.headshape)
geometry = ft_read_headshape(cfg.headshape);
elseif ~isempty(cfg.headshape)
geometry = fixpos(cfg.headshape);
else
ft_error('You must give a mesh, segmented MRI, sensor data type, or cfg.headshape');
end
headmodel = ft_headmodel_concentricspheres(geometry, 'conductivity', cfg.conductivity, 'fitind', cfg.fitind, 'order', cfg.order);
case 'halfspace'
if input_mesh || input_pos
geometry = data;
else
ft_error('a surface mesh is required as input for the halfspace method');
end
if isempty(cfg.point)
ft_error('cfg.point is required for halfspace method');
end
headmodel = ft_headmodel_halfspace(geometry, cfg.point, 'conductivity', cfg.conductivity, 'sourcemodel', cfg.submethod);
case 'infinite'
% this takes no input arguments
headmodel = ft_headmodel_infinite();
case {'localspheres' 'singlesphere' 'singleshell'}
cfg.grad = ft_getopt(cfg, 'grad'); % used for localspheres
% these three methods all require a single set of surface points
if input_mesh || input_pos
geometry = data;
elseif input_seg
tmpcfg = keepfields(cfg, {'numvertices' 'spmversion'});
if ~isempty(cfg.tissue)
% extract the specified surface
tmpcfg.tissue = cfg.tissue;
geometry = ft_prepare_mesh(tmpcfg, data);
else
% try to extract either the brain or scalp surface
geometry = [];
if isempty(geometry)
try
tmpcfg.tissue = 'brain';
geometry = ft_prepare_mesh(tmpcfg, data);
catch
me = lasterror;
if isequal(me.identifier, 'MATLAB:mex:ErrInvalidMEXFile')
% SPM8 mex file issues are common on macOS, these should not remain invisible
rethrow(me);
end
end
end
if isempty(geometry)
try
tmpcfg.tissue = 'scalp';
geometry = ft_prepare_mesh(tmpcfg, data);
catch
me = lasterror;
if isequal(me.identifier, 'MATLAB:mex:ErrInvalidMEXFile')
% SPM8 mex file issues are common on macOS, these should not remain invisible
rethrow(me);
end
end
end
if isempty(geometry)
ft_error('please specify cfg.tissue and pass an appropriate segmented MRI as input data')
end
end
elseif input_elec
geometry.pos = data.chanpos;
geometry.unit = data.unit;
elseif ~isempty(cfg.headshape)
% get the surface describing the head shape
[geometry.pos, geometry.tri] = headsurface([], [], 'headshape', cfg.headshape);
elseif ~isempty(cfg.headmodel)
% the CTF *.hdm file will be read further down
else
ft_error('this requires a mesh, set of surface points or a segmented mri');
end
switch cfg.method
case 'singlesphere'
if ~isempty(cfg.headmodel)
% read the volume conduction model from a CTF *.hdm file
tmp = ft_read_headmodel(cfg.headmodel);
try
% the single sphere is contained in the "orig" field
headmodel = [];
headmodel.r = tmp.orig.MEG_Sphere.RADIUS;
headmodel.o = [tmp.orig.MEG_Sphere.ORIGIN_X tmp.orig.MEG_Sphere.ORIGIN_Y tmp.orig.MEG_Sphere.ORIGIN_Z];
headmodel.unit = 'cm';
catch
ft_error('the volume conduction model in "%s" is invalid', cfg.headmodel);
end
else
% construct the volume conduction model
headmodel = ft_headmodel_singlesphere(geometry, 'conductivity', cfg.conductivity);
end % headmodel
case 'localspheres'
if ~isempty(cfg.headmodel)
% read the volume conduction model from a CTF *.hdm file
tmp = ft_read_headmodel(cfg.headmodel);
try
headmodel = [];
headmodel.label = tmp.label;
headmodel.r = tmp.r;
headmodel.o = tmp.o;
headmodel.unit = 'cm';
catch
ft_error('the volume conduction model in "%s" is invalid', cfg.headmodel);
end
else
% construct the volume conduction model
cfg.grad = ft_getopt(cfg, 'grad');
if isempty(cfg.grad)
ft_error('for cfg.method = %s, you must also supply cfg.grad', cfg.method);
end
headmodel = ft_headmodel_localspheres(geometry, cfg.grad, 'feedback', cfg.feedback, 'radius', cfg.radius, 'maxradius', cfg.maxradius, 'baseline', cfg.baseline, 'singlesphere', cfg.singlesphere);
end % headmodel
case 'singleshell'
cfg.order = ft_getopt(cfg, 'order');
if ~isfield(geometry, 'tri')
tmpcfg = [];
tmpcfg.headshape = geometry;
geometry = ft_prepare_mesh(tmpcfg);
end
headmodel = ft_headmodel_singleshell(geometry, 'order', cfg.order);
otherwise
ft_error('unsupported method %s', cfg.method);
end % switch method
case {'simbio'}
if input_elec || isfield(data, 'pos') || input_mesh
geometry = data; % more serious checks of validity of the mesh occur inside ft_headmodel_simbio
else
ft_error('you must provide a mesh with tetrahedral or hexahedral elements, where each element has a scalar or tensor conductivity');
end
headmodel = ft_headmodel_simbio(geometry, 'conductivity', cfg.conductivity);
case {'duneuro'}
if input_mesh
geometry = data; % more serious checks of validity of the mesh occur inside ft_headmodel_duneuro
else
error('You must provide a mesh with tetrahedral or hexahedral elements, where each element has a scalar or tensor conductivity');
end
headmodel = ft_headmodel_duneuro(geometry, 'grid_filename', cfg.grid_filename, 'tensors_filename', cfg.tensors_filename,...
'conductivity', cfg.conductivity, 'duneuro_settings', cfg.duneuro_settings);
case {'fns'}
if input_seg
data = ft_datatype_segmentation(data, 'segmentationstyle', 'indexed');
else
ft_error('segmented MRI must be given as data input')
end
sens = ft_fetch_sens(cfg, data);
headmodel = ft_headmodel_fns(data.seg, 'tissue', cfg.tissue, 'tissueval', cfg.tissueval, 'tissuecond', cfg.conductivity, 'sens', sens, 'transform', cfg.transform);
otherwise
ft_error('unsupported method "%s"', cfg.method);
end % switch method
% ensure that the geometrical units are specified
if ~ft_headmodeltype(headmodel, 'infinite')
headmodel = ft_determine_units(headmodel);
end
% do the general cleanup and bookkeeping at the end of the function
ft_postamble trackconfig
ft_postamble provenance
ft_postamble previous data
ft_postamble history headmodel