function [source] = ft_dipolefitting(cfg, data)
% FT_DIPOLEFITTING perform grid search and non-linear fit with one or multiple
% dipoles and try to find the location where the dipole model is best able
% to explain the measured EEG or MEG topography.
%
% This function will initially scan the whole brain with a single dipole on
% a regular coarse grid, and subsequently start at the most optimal location
% with a non-linear search. Alternatively you can specify the initial
% location of the dipole(s) and the non-linear search will start from there.
%
% Use as
% [source] = ft_dipolefitting(cfg, data)
%
% The configuration has the following general fields
% cfg.numdipoles = number, default is 1
% cfg.symmetry = 'x', 'y' or 'z' symmetry for two dipoles, can be empty (default = [])
% cfg.channel = Nx1 cell-array with selection of channels (default = 'all'),
% see FT_CHANNELSELECTION for details
% cfg.gridsearch = 'yes' or 'no', perform global grid search for initial
% guess for the dipole parameters (default = 'yes')
% cfg.nonlinear = 'yes' or 'no', perform nonlinear search for optimal
% dipole parameters (default = 'yes')
%
% If a grid search is performed, a source model needs to be specified. This should either be
% specified as cfg.sourcemodel (see below), or as a set of parameters to define a 3-D regular grid.
% In the latter case, a complete grid is constructed using FT_PREPARE_SOURCEMODEL. The specification
% of a regular 3-D grid, aligned with the axes of the head coordinate system, can be obtained with
% cfg.xgrid = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
% cfg.ygrid = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
% cfg.zgrid = vector (e.g. 0:1:20) or 'auto' (default = 'auto')
% cfg.resolution = number (e.g. 1 cm)
% If the source model destribes a triangulated cortical sheet, it is described as
% cfg.sourcemodel.pos = N*3 matrix with the vertex positions of the cortical sheet
% cfg.sourcemodel.tri = M*3 matrix that describes the triangles connecting the vertices
% Alternatively the position of the dipoles at locations of interest can be
% user-specified, for example obtained from an anatomical or functional MRI
% cfg.sourcemodel.pos = N*3 matrix with position of each source
% cfg.sourcemodel.inside = N*1 vector with boolean value whether grid point is inside brain (optional)
% cfg.sourcemodel.dim = [Nx Ny Nz] vector with dimensions in case of 3-D grid (optional)
%
% If you do not start with a grid search, you have to give a starting location
% for the nonlinear search
% cfg.dip.pos = initial dipole position, matrix of Ndipoles x 3
%
% The conventional approach is to fit dipoles to event-related averages, which
% within FieldTrip can be obtained from the FT_TIMELOCKANALYSIS or from
% the FT_TIMELOCKGRANDAVERAGE function. This has the additional options
% cfg.latency = [begin end] in seconds or 'all' (default = 'all')
% cfg.model = 'moving' or 'regional'
% A moving dipole model has a different position (and orientation) for each
% timepoint, or for each component. A regional dipole model has the same
% position for each timepoint or component, and a different orientation.
%
% You can also fit dipoles to the spatial topographies of an independent
% component analysis, obtained from the FT_COMPONENTANALYSIS function.
% This has the additional options
% cfg.component = array with numbers (can be empty -> all)
%
% You can also fit dipoles to the spatial topographies that are present
% in the data in the frequency domain, which can be obtained using the
% FT_FREQANALYSIS function. This has the additional options
% cfg.frequency = single number (in Hz)
%
% Low level details of the fitting can be specified in the cfg.dipfit structure
% cfg.dipfit.display = level of display, can be 'off', 'iter', 'notify' or 'final' (default = 'iter')
% cfg.dipfit.optimfun = function to use, can be 'fminsearch' or 'fminunc' (default is determined automatic)
% cfg.dipfit.maxiter = maximum number of function evaluations allowed (default depends on the optimfun)
% cfg.dipfit.checkinside = boolean, check that the dipole remains in the source compartment (default = false)
%
% Optionally, you can modify the leadfields by reducing the rank, i.e. remove the weakest orientation
% cfg.reducerank = 'no', or number (default = 3 for EEG, 2 for MEG)
% cfg.backproject = 'yes' or 'no', determines when reducerank is applied whether the
% lower rank leadfield is projected back onto the original linear
% subspace, or not (default = 'yes')
%
% The volume conduction model of the head should be specified as
% cfg.headmodel = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The EEG or MEG sensor positions can be present in the data or can be specified as
% cfg.elec = structure with electrode positions or filename, see FT_READ_SENS
% cfg.grad = structure with gradiometer definition or filename, see FT_READ_SENS
%
% To facilitate data-handling and distributed computing you can use
% cfg.inputfile = ...
% cfg.outputfile = ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_SOURCEANALYSIS, FT_PREPARE_LEADFIELD, FT_PREPARE_HEADMODEL
% TODO change the output format, more suitable would be something like:
% dip.label
% dip.time
% dip.avg (instead of Vdata)
% dip.dip.pos
% dip.dip.mom
% dip.dip.model, or dip.dip.avg
% dip.dimord
% Undocumented local options:
% cfg.dipfit.constr = Source model constraints, depends on cfg.symmetry
% Optionally, you can include a noise covariance structure to sphere the data (is useful when using both
% magnetometers and gradiometers to fit your dipole)
% cfg.dipfit.noisecov = noise covariance matrix, see e.g. FT_TIMELOCK_ANALYSIS
% Copyright (C) 2004-2013, Robert Oostenveld
%
% 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 .
%
% $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 debug
ft_preamble loadvar data
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 data is valid for this function
data = ft_checkdata(data, 'datatype', {'comp', 'timelock', 'freq'}, 'feedback', 'yes');
% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden', {'channels'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed', {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed', {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed', {'optofile', 'opto'});
cfg = ft_checkconfig(cfg, 'renamed', {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'vol', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed', {'grid', 'sourcemodel'});
% get the defaults
cfg.channel = ft_getopt(cfg, 'channel', 'all');
cfg.component = ft_getopt(cfg, 'component', 'all'); % for comp input
cfg.frequency = ft_getopt(cfg, 'frequency'); % for freq input
cfg.latency = ft_getopt(cfg, 'latency', 'all'); % for timelock input
cfg.feedback = ft_getopt(cfg, 'feedback', 'text');
cfg.gridsearch = ft_getopt(cfg, 'gridsearch', 'yes');
cfg.nonlinear = ft_getopt(cfg, 'nonlinear', 'yes');
cfg.symmetry = ft_getopt(cfg, 'symmetry');
cfg.dipfit = ft_getopt(cfg, 'dipfit', []); % the default for this is handled below
cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'}); % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.sourcemodel.unit by the subsequent createsubcfg
% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});
% determine data type
iscomp = ft_datatype(data, 'comp'); % it can also be raw+comp, timelock+comp or freq+comp
isfreq = ft_datatype(data, 'freq'); % it might also be freq+comp, in that case it should be treated as component data
istimelock = ft_datatype(data, 'timelock'); % it might also be timelock+comp, in that case it should be treated as component data
% the default for this depends on the data type
if ~isfield(cfg, 'model')
if iscomp
% each component is fitted independently
cfg.model = 'moving';
elseif isfreq
% fit the data with a dipole at one location
cfg.model = 'regional';
elseif istimelock
% fit the data with a dipole at one location
cfg.model = 'regional';
end
end
if ~isfield(cfg, 'numdipoles')
if isfield(cfg, 'dip')
cfg.numdipoles = size(cfg.dip(1).pos,1);
else
cfg.numdipoles = 1;
end
end
% set up the symmetry constraints
if ~isempty(cfg.symmetry)
if cfg.numdipoles~=2
ft_error('symmetry constraints are only supported for two-dipole models');
elseif strcmp(cfg.symmetry, 'x')
% this structure is passed onto the low-level FT_INVERSE_DIPOLEFIT function
cfg.dipfit.constr.reduce = [1 2 3]; % select the parameters [x1 y1 z1]
cfg.dipfit.constr.expand = [1 2 3 1 2 3]; % repeat them as [x1 y1 z1 x1 y1 z1]
cfg.dipfit.constr.mirror = [1 1 1 -1 1 1]; % multiply each of them with 1 or -1, resulting in [x1 y1 z1 -x1 y1 z1]
elseif strcmp(cfg.symmetry, 'y')
% this structure is passed onto the low-level FT_INVERSE_DIPOLEFIT function
cfg.dipfit.constr.reduce = [1 2 3]; % select the parameters [x1 y1 z1]
cfg.dipfit.constr.expand = [1 2 3 1 2 3]; % repeat them as [x1 y1 z1 x1 y1 z1]
cfg.dipfit.constr.mirror = [1 1 1 1 -1 1]; % multiply each of them with 1 or -1, resulting in [x1 y1 z1 x1 -y1 z1]
elseif strcmp(cfg.symmetry, 'z')
% this structure is passed onto the low-level FT_INVERSE_DIPOLEFIT function
cfg.dipfit.constr.reduce = [1 2 3]; % select the parameters [x1 y1 z1]
cfg.dipfit.constr.expand = [1 2 3 1 2 3]; % repeat them as [x1 y1 z1 x1 y1 z1]
cfg.dipfit.constr.mirror = [1 1 1 1 1 -1]; % multiply each of them with 1 or -1, resulting in [x1 y1 z1 x1 y1 -z1]
else
ft_error('unrecognized symmetry constraint');
end
elseif ~isfield(cfg, 'dipfit') || ~isfield(cfg.dipfit, 'constr')
% no symmetry constraints have been specified
cfg.dipfit.constr = [];
end
if ft_getopt(cfg.dipfit.constr, 'sequential', false) && strcmp(cfg.model, 'moving')
ft_error('the moving dipole model does not combine with the sequential constraint')
% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3119
end
if iscomp
% transform the data into a representation on which the timelocked dipole fit can perform its trick
data = comp2timelock(cfg, data);
% default component selection is all components
if ischar(cfg.component) && strcmp(cfg.component, 'all')
cfg.component = (1:size(data.avg, 2));
end
elseif isfreq
% transform the data into a representation on which the timelocked dipole fit can perform its trick
data = freq2timelock(cfg, data);
elseif istimelock
% no transformation is needed
end
% collect and preprocess the electrodes/gradiometer and head model
% this will also update cfg.channel to match the electrodes/gradiometers
[headmodel, sens, cfg] = prepare_headmodel(cfg, data);
% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD and FT_INVERSE_DIPOLEFIT
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank', ft_getopt(cfg, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject', ft_getopt(cfg, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize', ft_getopt(cfg, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight', ft_getopt(cfg, 'weight'));
% construct the low-level options for the dipole fitting as key-value pairs, these are passed to FT_INVERSE_DIPOLEFIT
dipfitopt = ft_cfg2keyval(cfg.dipfit);
% select the desired channels, ordered according to the sensor structure or configuration
[selcfg, seldata] = match_str(cfg.channel, data.label);
% take the selected channels from the data structure
Vdata = data.avg(seldata, :);
% sphere the date using the noise covariance matrix supplied, if any
% this affects both the gridsearch and the nonlinear optimization
noisecov = ft_getopt(cfg.dipfit, 'noisecov');
if ~isempty(noisecov)
[u, s] = svd(noisecov);
tol = max(size(noisecov)) * eps(norm(s, inf));
s = diag(s);
r1 = sum(s > tol) + 1;
s(1:(r1 - 1)) = 1 ./ sqrt(s(1:(r1 - 1)));
s(r1:end) = 0;
sphere = diag(s) * u';
% apply the sphering to the data
Vdata = sphere * Vdata;
% apply the sphering as a pre-multiplication to the sensor definition
montage = [];
montage.labelold = cfg.channel;
montage.labelnew = cfg.channel;
montage.tra = sphere;
sens = ft_apply_montage(sens, montage, 'balancename', 'sphering');
end
if iscomp
% select the desired component topographies
Vdata = Vdata(:, cfg.component);
elseif isfreq
% the desired frequencies have already been selected
Vdata = Vdata(:, :);
elseif istimelock
% select the desired latencies
if ischar(cfg.latency) && strcmp(cfg.latency, 'all')
cfg.latency = data.time([1 end]);
end
tbeg = nearest(data.time, cfg.latency(1));
tend = nearest(data.time, cfg.latency(end));
cfg.latency = [data.time(tbeg) data.time(tend)];
Vdata = Vdata(:, tbeg:tend);
end
nchans = size(Vdata,1);
ntime = size(Vdata,2);
Vmodel = zeros(nchans, ntime);
ft_info('selected %d channels\n', nchans);
ft_info('selected %d topographies\n', ntime);
if nchans0.001)
ft_warning('the EEG data is not average referenced, correcting this');
end
Vdata = avgref(Vdata);
end
% set to zeros if no initial dipole was specified
if ~isfield(cfg, 'dip')
cfg.dip.pos = zeros(cfg.numdipoles, 3);
cfg.dip.mom = zeros(3*cfg.numdipoles, 1);
end
% set to zeros if no initial dipole position was specified
if ~isfield(cfg.dip, 'pos')
cfg.dip.pos = zeros(cfg.numdipoles, 3);
end
% set to zeros if no initial dipole moment was specified
if ~isfield(cfg.dip, 'mom')
cfg.dip.mom = zeros(3*cfg.numdipoles, 1);
end
% check the specified dipole model
if numel(cfg.dip.pos)~=cfg.numdipoles*3 || numel(cfg.dip.mom)~=cfg.numdipoles*3
ft_error('inconsistent number of dipoles in configuration')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform the dipole scan, this is usefull for generating an initial guess
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(cfg.gridsearch, 'yes')
% test whether we have a valid configuration for dipole scanning
if cfg.numdipoles==1
% this is ok
elseif cfg.numdipoles==2 && ~isempty(cfg.dipfit.constr)
% this is also ok
elseif isfield(cfg.sourcemodel, 'pos') && size(cfg.sourcemodel.pos,2)==cfg.numdipoles*3
% this is also ok
else
ft_error('dipole scanning is only possible for a single dipole or a symmetric dipole pair');
end
if isfield(cfg.sourcemodel, 'leadfield')
ft_notice('using precomputed leadfields for the gridsearch');
sourcemodel = keepfields(cfg.sourcemodel, {'pos', 'tri', 'dim', 'unit', 'coordsys', 'inside', 'leadfield', 'leadfielddimord', 'label'});
% select the channels corresponding to the data and the user configuration
tmpcfg = keepfields(cfg, 'channel');
sourcemodel = ft_selectdata(tmpcfg, sourcemodel);
% sort the channels to be consistent with the data
[dum, chansel] = match_str(data.label, sourcemodel.label);
sourcemodel.label = sourcemodel.label(chansel);
for i=1:numel(sourcemodel.leadfield)
if ~isempty(sourcemodel.leadfield{i})
sourcemodel.leadfield{i} = sourcemodel.leadfield{i}(chansel, :);
end
end
% ensure that the channels are consistent with the data
assert(isequal(sourcemodel.label, cfg.channel), 'cannot match the channels in the sourcemodel to those in the data')
else
ft_notice('computing the leadfields for the gridsearch on the fly');
% construct the dipole positions on which the source reconstruction will be done
tmpcfg = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
tmpcfg.headmodel = headmodel;
if ft_senstype(sens, 'eeg')
tmpcfg.elec = sens;
elseif ft_senstype(sens, 'meg')
tmpcfg.grad = sens;
end
sourcemodel = ft_prepare_sourcemodel(tmpcfg);
end % if precomputed leadfield or not
ngrid = size(sourcemodel.pos,1);
switch cfg.model
case 'regional'
sourcemodel.error = nan(ngrid, 1);
case 'moving'
sourcemodel.error = nan(ngrid, ntime);
otherwise
ft_error('unsupported cfg.model');
end
insideindx = find(sourcemodel.inside);
ft_progress('init', cfg.feedback, 'scanning grid');
for i=1:length(insideindx)
ft_progress(i/length(insideindx), 'scanning grid location %d/%d\n', i, length(insideindx));
thisindx = insideindx(i);
if isfield(sourcemodel, 'leadfield')
% reuse the previously computed leadfield
lf = sourcemodel.leadfield{thisindx};
else
lf = ft_compute_leadfield(sourcemodel.pos(thisindx,:), sens, headmodel, leadfieldopt{:});
end
% the model is V=lf*mom+noise, therefore mom=pinv(lf)*V estimates the
% dipole moment this makes the model potential U=lf*pinv(lf)*V and the
% model error is norm(V-U) = norm(V-lf*pinv(lf)*V) = norm((eye-lf*pinv(lf))*V)
if any(isnan(lf(:)))
% this might happen if one of the dipole locations of the grid is
% outside the brain compartment
lf(:) = 0;
end
switch cfg.model
case 'regional'
% sum the error over all latencies
sourcemodel.error(thisindx,1) = sum(sum(((eye(nchans)-lf*pinv(lf))*Vdata).^2));
case 'moving'
% remember the error for each latency independently
sourcemodel.error(thisindx,:) = sum(((eye(nchans)-lf*pinv(lf))*Vdata).^2);
otherwise
ft_error('unsupported cfg.model');
end % switch model
end % looping over the grid
ft_progress('close');
switch cfg.model
case 'regional'
% find the source position with the minimum error
[err, indx] = min(sourcemodel.error);
dip.pos = sourcemodel.pos(indx,:); % note that for a symmetric dipole pair this results in a vector
dip.pos = reshape(dip.pos,3,cfg.numdipoles)'; % convert to a Nx3 array
dip.mom = zeros(cfg.numdipoles*3,1); % set the dipole moment to zero
if cfg.numdipoles==1
ft_info('found minimum after scanning on grid point [%g %g %g]\n', dip.pos(1), dip.pos(2), dip.pos(3));
elseif cfg.numdipoles==2
ft_info('found minimum after scanning on grid point [%g %g %g; %g %g %g]\n', dip.pos(1,1), dip.pos(1,2), dip.pos(1,3), dip.pos(2,1), dip.pos(2,2), dip.pos(2,3));
end
case 'moving'
for t=1:ntime
% find the source position with the minimum error
[err, indx] = min(sourcemodel.error(:,t));
dip(t).pos = sourcemodel.pos(indx,:); % note that for a symmetric dipole pair this results in a vector
dip(t).pos = reshape(dip(t).pos,3,cfg.numdipoles)'; % convert to a Nx3 array
dip(t).mom = zeros(cfg.numdipoles*3,1); % set the dipole moment to zero
if cfg.numdipoles==1
ft_info('found minimum after scanning for topography %d on grid point [%g %g %g]\n', t, dip(t).pos(1), dip(t).pos(2), dip(t).pos(3));
elseif cfg.numdipoles==2
ft_info('found minimum after scanning for topography %d on grid point [%g %g %g; %g %g %g]\n', t, dip(t).pos(1,1), dip(t).pos(1,2), dip(t).pos(1,3), dip(t).pos(2,1), dip(t).pos(2,2), dip(t).pos(2,3));
end
end
otherwise
ft_error('unsupported cfg.model');
end % switch model
elseif strcmp(cfg.gridsearch, 'no')
% there is no grid needed for dipole scanning
sourcemodel = [];
% use the initial guess supplied in the configuration for the remainder
switch cfg.model
case 'regional'
dip = cfg.dip;
case 'moving'
for t=1:ntime
dip(t) = cfg.dip;
end
otherwise
ft_error('unsupported cfg.model');
end % switch model
end % if gridsearch yes/no
% multiple dipoles can be represented either as a 1x(N*3) vector or as a Nx3 matrix,
% i.e. [x1 y1 z1 x2 y2 z2] or [x1 y1 z1; x2 y2 z2]
switch cfg.model
case 'regional'
dip = fixdipole(dip);
case 'moving'
for t=1:ntime
dip(t) = fixdipole(dip(t));
end
otherwise
ft_error('unsupported cfg.model');
end % switch model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform the non-linear fit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(cfg.nonlinear, 'yes')
switch cfg.model
case 'regional'
% perform the non-linear dipole fit for all latencies together
% catch errors due to non-convergence
try
dip = ft_inverse_dipolefit(dip, sens, headmodel, Vdata, dipfitopt{:}, leadfieldopt{:});
success = 1;
if cfg.numdipoles==1
ft_info('found minimum after non-linear optimization on [%g %g %g]\n', dip.pos(1), dip.pos(2), dip.pos(3));
elseif cfg.numdipoles==2
ft_info('found minimum after non-linear optimization on [%g %g %g; %g %g %g]\n', dip.pos(1,1), dip.pos(1,2), dip.pos(1,3), dip.pos(2,1), dip.pos(2,2), dip.pos(2,3));
end
catch
success = 0;
disp(lasterr);
end
case 'moving'
% perform the non-linear dipole fit for each latency independently
% instead of using dip(t) = ft_inverse_dipolefit(dip(t),...), I am using temporary variables dipin and dipout
% to prevent errors like "Subscripted assignment between dissimilar structures"
dipin = dip;
for t=1:ntime
% catch errors due to non-convergence
try
dipout(t) = ft_inverse_dipolefit(dipin(t), sens, headmodel, Vdata(:,t), dipfitopt{:}, leadfieldopt{:});
success(t) = 1;
if cfg.numdipoles==1
ft_info('found minimum after non-linear optimization for topography %d on [%g %g %g]\n', t, dipout(t).pos(1), dipout(t).pos(2), dipout(t).pos(3));
elseif cfg.numdipoles==2
ft_info('found minimum after non-linear optimization for topography %d on [%g %g %g; %g %g %g]\n', t, dipout(t).pos(1,1), dipout(t).pos(1,2), dipout(t).pos(1,3), dipout(t).pos(2,1), dipout(t).pos(2,2), dipout(t).pos(2,3));
end
catch
% keep the position and moment according to the initial guess
dipout(t).pos = dipin(t).pos;
dipout(t).mom = dipin(t).mom;
success(t) = 0;
disp(lasterr);
end
end
dip = dipout;
clear dipin dipout
otherwise
ft_error('unsupported cfg.model');
end % switch model
end % if nonlinear
if strcmp(cfg.nonlinear, 'no')
% the optimal dipole positions are either obtained from scanning
% or from the initial configured specified by the user
switch cfg.model
case 'regional'
success = 1;
case 'moving'
success = ones(1,ntime);
otherwise
ft_error('unsupported cfg.model');
end % switch model
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the model potential distribution and the residual variance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch cfg.model
case 'regional'
if success
% re-compute the leadfield in order to compute the model potential and dipole moment
lf = ft_compute_leadfield(dip.pos, sens, headmodel, leadfieldopt{:});
if isfield(dip, 'mom') && isfield(dip, 'ampl')
% the orientation and amplitude have already been estimated, this applies to the case of a fixed dipole orientation
dip.pot = (lf * dip.mom) * dip.ampl;
else
% compute all details of the final dipole model using linear estimation
dip.mom = pinv(lf)*Vdata;
dip.pot = lf*dip.mom;
end
dip.rv = rv(Vdata, dip.pot);
Vmodel = dip.pot;
end
case 'moving'
for t=1:ntime
if success(t)
% re-compute the leadfield in order to compute the model potential and dipole moment
lf = ft_compute_leadfield(dip(t).pos, sens, headmodel, leadfieldopt{:});
% compute all details of the final dipole model
dip(t).mom = pinv(lf)*Vdata(:,t);
dip(t).pot = lf*dip(t).mom;
dip(t).rv = rv(Vdata(:,t), dip(t).pot);
Vmodel(:,t) = dip(t).pot;
end
end
otherwise
ft_error('unsupported cfg.model');
end % switch model
switch cfg.model
case 'regional'
if isfreq
% the matrix with the dipole moment is encrypted and cannot be interpreted straight away
% reconstruct the frequency representation of the data at the source level
if isfield(dip, 'mom') && isfield(dip, 'ampl')
% this applies to the case of a fixed dipole orientation
[dip.pow, dip.csd, dip.fourier] = timelock2freq(dip.mom * dip.ampl);
else
[dip.pow, dip.csd, dip.fourier] = timelock2freq(dip.mom);
end
end
case 'moving'
if isfreq
% although this is technically possible so far, it does not make any sense
ft_warning('a moving dipole model in the frequency domain is not supported');
end
otherwise
ft_error('unsupported cfg.model');
end % switch model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% collect the results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
source.label = cfg.channel; % these channels were used in fitting
source.dip = dip;
source.Vdata = Vdata; % FIXME this should be renamed (if possible w.r.t. EEGLAB)
source.Vmodel = Vmodel; % FIXME this should be renamed (if possible w.r.t. EEGLAB)
% the units of the fitted source are the same as the units of the headmodel and the sensor array
for i=1:length(source.dip)
if isfield(headmodel, 'unit')
source.dip(i).unit = headmodel.unit;
elseif isfield(sourcemodel, 'unit')
source.dip(i).unit = sourcemodel.unit;
end
end
% assign a latency, frequeny or component axis to the output
if iscomp
source.component = cfg.component;
% FIXME assign Vdata to an output variable, idem for the model potential
elseif isfreq
source.freq = cfg.frequency;
source.dimord = 'chan_freq';
% FIXME assign Vdata to an output variable, idem for the model potential
elseif istimelock
tbeg = nearest(data.time, cfg.latency(1));
tend = nearest(data.time, cfg.latency(end));
source.time = data.time(tbeg:tend);
source.dimord = 'chan_time';
end
% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous data
ft_postamble provenance source
ft_postamble history source
ft_postamble savevar source