/
ft_prepare_leadfield.m
372 lines (322 loc) · 17.9 KB
/
ft_prepare_leadfield.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
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 = 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, {'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