function [sourcemodel, cfg] = ft_prepare_leadfield(cfg, data)

% FT_PREPARE_LEADFIELD computes the forward model for many dipole locations

% on a regular 2D or 3D sourcemodel and stores it for efficient inverse modelling

%

% Use as

% [sourcemodel] = ft_prepare_leadfield(cfg, data)

%

% It is necessary to input the data on which you want to perform the inverse

% computations, since that data generally contain the gradiometer information and

% information about the channels that should be included in the forward model

% computation. The data structure can be either obtained from FT_PREPROCESSING,

% FT_FREQANALYSIS or FT_TIMELOCKANALYSIS. If the data is empty, all channels will be

% included in the forward model.

%

% The configuration should contain

% cfg.channel = Nx1 cell-array with selection of channels (default = 'all'),

% see FT_CHANNELSELECTION for details

%

% The positions of the sources can be specified as a regular 3-D

% sourcemodel that is aligned with the axes of the head coordinate system

% 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) for automatic sourcemodel generation

%

% Alternatively the position of a few sources at locations of interest can

% be 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 sourcemodel point is inside brain (optional)

% cfg.sourcemodel.dim = [Nx Ny Nz] vector with dimensions in case of 3-D sourcemodel (optional)

%

% 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

%

% Optionally, you can modify the leadfields by reducing the rank (i.e. remove the

% weakest orientation), or by normalizing each column.

% 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')

% cfg.normalize = 'yes' or 'no' (default = 'no')

% cfg.normalizeparam = depth normalization parameter (default = 0.5)

% cfg.weight = number or Nx1 vector, weight for each dipole position to compensate

% for the size of the corresponding patch (default = 1)

%

% Depending on the type of headmodel, some additional options may be

% specified.

%

% For OPENMEEG based headmodels:

% cfg.openmeeg.batchsize = scalar (default 1e4), number of dipoles

% for which the leadfield is computed in a

% single call to the low-level code. Trades off

% memory efficiency for speed.

% cfg.openmeeg.dsm = 'no'/'yes', reuse existing DSM if provided

% cfg.openmeeg.keepdsm = 'no'/'yes', option to retain DSM (no by default)

% cfg.openmeeg.nonadaptive = 'no'/'yes'

%

% For SINGLESHELL based headmodels:

% cfg.singleshell.batchsize = scalar or 'all' (default 1), number of dipoles

% for which the leadfield is computed in a

% single call to the low-level code. Trades off

% memory efficiency for speed.

%

% To facilitate data-handling and distributed computing you can use

% cfg.inputfile = ...

% If you specify this option the input data will be read from a *.mat

% file on disk. This mat files should contain only a single variable named 'data',

% corresponding to the input structure.

%

% See also FT_SOURCEANALYSIS, FT_DIPOLEFITTING, FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% Undocumented local options:

% cfg.feedback

% cfg.sel50p = 'no' (default) or 'yes'

% cfg.lbex = 'no' (default) or a number that corresponds with the radius

% cfg.mollify = 'no' (default) or a number that corresponds with the FWHM

% 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 <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 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

% the data can be passed as input arguments or can be read from disk

hasdata = exist('data', 'var');

if ~hasdata

% the data variable will be passed to the prepare_headmodel function below

% where it would be used for channel selection

data = [];

else

% check if the input data is valid for this function

data = ft_checkdata(data);

end

% 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', {'hdmfile', 'headmodel'});

cfg = ft_checkconfig(cfg, 'renamed', {'vol', 'headmodel'});

cfg = ft_checkconfig(cfg, 'renamed', {'grid', 'sourcemodel'});

cfg = ft_checkconfig(cfg, 'renamed', {'om', 'openmeeg'});

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, 'deprecated', {'patchindx', 'patchsize'});

% set the defaults

cfg.lbex = ft_getopt(cfg, 'lbex', 'no');

cfg.sel50p = ft_getopt(cfg, 'sel50p', 'no');

cfg.feedback = ft_getopt(cfg, 'feedback', 'text');

cfg.mollify = ft_getopt(cfg, 'mollify', 'no');

cfg.patchsvd = ft_getopt(cfg, 'patchsvd', 'no');

cfg.checkinside = ft_getopt(cfg, 'checkinside'); % default is set in FT_PREPARE_SOURCEMODEL

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'});

% this code expects the inside to be represented as a logical array

cfg.sourcemodel = ft_checkconfig(cfg.sourcemodel, 'renamed', {'pnt' 'pos'});

cfg = ft_checkconfig(cfg, 'inside2logical', 'yes');

if strcmp(cfg.sel50p, 'yes') && strcmp(cfg.lbex, 'yes')

ft_error('subspace projection with either lbex or sel50p is mutually exclusive');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% collect and preprocess the electrodes/gradiometer and head model

[headmodel, sens, cfg] = prepare_headmodel(cfg, data);

% construct the sourcemodel for which the leadfield will be computed

tmpcfg = keepfields(cfg, {'checkinside', '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);

% find the indices of all sourcemodel points that are inside the brain

insideindx = find(sourcemodel.inside);

% check whether units are equal (NOTE: this was previously not required,

% this check can be removed if the underlying bug is resolved. See

% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387

if ~isfield(headmodel, 'unit') || ~isfield(sourcemodel, 'unit') || ~isfield(sens, 'unit')

ft_warning('cannot determine the units of all geometric objects required for leadfield computation (headmodel, sourcemodel, sensor configuration). THIS CAN LEAD TO WRONG RESULTS! (refer to http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387)');

else

if ~strcmp(headmodel.unit, sourcemodel.unit) || ~strcmp(sourcemodel.unit, sens.unit)

ft_error('geometric objects (headmodel, sourcemodel, sensor configuration) are not expressed in the same units (this used to be allowed, and will be again in the future, but for now there is a bug which prevents a correct leadfield from being computed; see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2387)');

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD

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'));

if ft_headmodeltype(headmodel, 'openmeeg')

ft_hastoolbox('openmeeg', 1); % add to path (if not yet on path)

cfg.openmeeg = ft_getopt(cfg, 'openmeeg', []);

batchsize = ft_getopt(cfg.openmeeg, 'batchsize', 1e4); % number of voxels per DSM batch; set to e.g. 1000 if not much RAM available

keepdsm = ft_getopt(cfg.openmeeg, 'keepdsm', 'no'); % retain DSM

leadfieldopt = ft_setopt(leadfieldopt, 'dsm', ft_getopt(cfg.openmeeg, 'dsm')); % reuse existing DSM if provided

leadfieldopt = ft_setopt(leadfieldopt, 'nonadaptive', ft_getopt(cfg.openmeeg, 'nonadaptive', 'no'));

% repeated system calls to the openmeeg executable makes it rather slow, calling it once is much more efficient

fprintf('calculating leadfield for %d positions at a time, this may take a while...\n', batchsize);

% a dsm in the input cfg currently assumes that the 'content' of the dsm matches

% exactly the positions that were passed in sourcemodel.pos. This is not guaranteed

% of course. If anything, it would make sense to represent the dsm in the input

% sourcemodel, so that the pos, and dsm are bound together. Still no guarantee, but

% better than nothing. This means that the cfg.openmeeg.dsm option should be

% deprecated

ndip = length(insideindx);

numchunks = ceil(ndip/batchsize);

dippos = sourcemodel.pos(insideindx,:);

if(numchunks > 1)

if istrue(keepdsm)

ft_warning('Keeping DSM output not supported when the computation is split into batches')

end

keepdsm = false;

end

% DSM computation is computationally intensive: As it can be reused with same voxel

% sourcemodel (i.e. if voxels are defined in MRI coordinates rather than MEG

% coordinates), optionally save result. Dense voxel grids may require several

% gigabytes of RAM, so optionally split into smaller batches

dsm = ft_getopt(leadfieldopt, 'dsm');

if istrue(keepdsm) && ~isempty(dsm)

% dsm needs to be computed outside ft_compute_leadfield, because it needs to be passed on in the output

dsm = ft_sysmat_openmeeg(dippos, headmodel, sens, ft_getopt(leadfieldopt, 'nonadaptive'));

leadfieldopt = ft_setopt(leadfieldopt, 'dsm', dsm);

end

sourcemodel.leadfield = cell(size(sourcemodel.pos,1),1);

ft_progress('init', cfg.feedback, 'computing leadfield');

for k = 1:numchunks

ft_progress(k/numchunks, 'computing leadfield %d/%d\n', k, numchunks);

diprange = (((k-1)*batchsize + 1):(min(k*batchsize,ndip)));

tmp = ft_compute_leadfield(dippos(diprange,:), sens, headmodel, leadfieldopt{:});

% distribute the columns of the leadfield matrix over the individual dipole positions

% avoid using the options reducerank and backproject, see https://github.com/fieldtrip/fieldtrip/issues/1410#issuecomment-646994620

[m, n] = size(tmp);

sourcemodel.leadfield(insideindx(diprange)) = mat2cell(tmp, m, repmat(n/numel(diprange), 1, numel(diprange)));

end

ft_progress('close');

if istrue(keepdsm)

% retain DSM in cfg if desired -> FIXME this should not be kept in

% the cfg. If anything, it is sourcemodel specific, so it should be

% retained in the output sourcemodel

cfg.openmeeg.dsm = dsm;

end

elseif ft_headmodeltype(headmodel, 'singleshell')

cfg.singleshell = ft_getopt(cfg, 'singleshell', []);

batchsize = ft_getopt(cfg.singleshell, 'batchsize', 1);

if ischar(batchsize) && strcmp(batchsize, 'all')

batchsize = length(insideindx);

end

dippos = sourcemodel.pos(insideindx,:);

ndip = length(insideindx);

numchunks = ceil(ndip/batchsize);

sourcemodel.leadfield = cell(size(sourcemodel.pos,1),1);

ft_progress('init', cfg.feedback, 'computing leadfield');

for k = 1:numchunks

ft_progress(k/numchunks, 'computing leadfield %d/%d\n', k, numchunks);

diprange = (((k-1)*batchsize + 1):(min(k*batchsize,ndip)));

tmp = ft_compute_leadfield(dippos(diprange,:), sens, headmodel, leadfieldopt{:});

% distribute the columns of the leadfield matrix over the individual dipole positions

% avoid using the options reducerank and backproject, see https://github.com/fieldtrip/fieldtrip/issues/1410#issuecomment-646994620

[m, n] = size(tmp);

sourcemodel.leadfield(insideindx(diprange)) = mat2cell(tmp, m, repmat(n/numel(diprange), 1, numel(diprange)));

end

ft_progress('close');

elseif ft_headmodeltype(headmodel, 'duneuro')

% ft_hastoolbox('duneuro', 1); %does not look necessary here? check

% repeated system calls to the duneuro executable makes it rather slow

% calling it once for all dipoles is much more efficient

% find the indices of all grid points that are inside the brain

insideindx = find(sourcemodel.inside);

ft_progress('init', cfg.feedback, 'computing leadfield');

% compute the leadfield on all grid positions inside the brain

lf = ft_compute_leadfield(sourcemodel.pos(insideindx,:), sens, headmodel, 'reducerank', cfg.reducerank, 'normalize', cfg.normalize, 'normalizeparam', cfg.normalizeparam, 'backproject', cfg.backproject);

lf = mat2cell(lf, size(lf,1), repmat(3,1,size(lf,2)/3));

sourcemodel.leadfield(sourcemodel.inside) = lf;

for i=1:length(insideindx)

thisindx = insideindx(i);

if isfield(cfg, 'grid') && isfield(cfg.grid, 'mom')

% multiply with the normalized dipole moment to get the leadfield in the desired orientation

sourcemodel.leadfield{thisindx} = sourcemodel.leadfield{thisindx} * sourcemodel.mom(:,thisindx);

end

end % for all grid locations inside the brain

elseif ft_headmodeltype(headmodel, 'interpolate')

lf = ft_compute_leadfield(sourcemodel.pos(insideindx,:), sens, headmodel, leadfieldopt{:});

sourcemodel.leadfield(insideindx) = mat2cell(lf, 3, 3.*ones(1,numel(insideindx)));

else

ft_progress('init', cfg.feedback, 'computing leadfield');

for i=1:length(insideindx)

% compute the leadfield on all sourcemodel positions inside the brain

ft_progress(i/length(insideindx), 'computing leadfield %d/%d\n', i, length(insideindx));

thisindx = insideindx(i);

sourcemodel.leadfield{thisindx} = ft_compute_leadfield(sourcemodel.pos(thisindx,:), sens, headmodel, leadfieldopt{:});

end % for all sourcemodel locations inside the brain

ft_progress('close');

end

if isfield(cfg, 'sourcemodel') && isfield(cfg.sourcemodel, 'mom')

for i=1:length(insideindx)

% multiply with the normalized dipole moment to get the leadfield in the desired orientation

% FIXME mom and ori seem to be mixed up here, see https://github.com/fieldtrip/fieldtrip/issues/1399

thisindx = insideindx(i);

sourcemodel.leadfield{thisindx} = sourcemodel.leadfield{thisindx} * sourcemodel.mom(:,thisindx);

end

end

% represent the leadfield for positions outside the brain as empty array

sourcemodel.leadfield(~sourcemodel.inside) = {[]};

% add the label of the channels

sourcemodel.label = sens.label;

sourcemodel.leadfielddimord = '{pos}_chan_ori';

% mollify the leadfields

if ~isequal(cfg.mollify, 'no')

sourcemodel = mollify(cfg, sourcemodel);

end

% combine leadfields in patches and do an SVD on them

if ~isequal(cfg.patchsvd, 'no')

sourcemodel = patchsvd(cfg, sourcemodel);

end

% compute the 50 percent channel selection subspace projection

if ~isequal(cfg.sel50p, 'no')

sourcemodel = sel50p(cfg, sourcemodel, sens);

end

% compute the local basis function expansion (LBEX) subspace projection

if ~isequal(cfg.lbex, 'no')

sourcemodel = lbex(cfg, sourcemodel);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% do the general cleanup and bookkeeping at the end of the function

ft_postamble debug

ft_postamble previous data

ft_postamble provenance sourcemodel

ft_postamble history sourcemodel

ft_postamble savevar sourcemodel