/
ft_headmodel_dipoli.m
267 lines (236 loc) · 7.95 KB
/
ft_headmodel_dipoli.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
function headmodel = ft_headmodel_dipoli(mesh, varargin)
% FT_HEADMODEL_DIPOLI creates a volume conduction model of the head
% using the boundary element method (BEM) for EEG. This function takes
% as input the triangulated surfaces that describe the boundaries and
% returns as output a volume conduction model which can be used to
% compute leadfields.
%
% This implements
% Oostendorp TF, van Oosterom A. "Source parameter estimation in
% inhomogeneous volume conductors of arbitrary shape." IEEE Trans
% Biomed Eng. 1989 Mar;36(3):382-91.
%
% The implementation of this function uses an external command-line
% executable with the name "dipoli" which is provided by Thom Oostendorp.
%
% Use as
% headmodel = ft_headmodel_dipoli(mesh, ...)
%
% The mesh is given as a boundary or a struct-array of boundaries (surfaces)
%
% Optional input arguments should be specified in key-value pairs and can
% include
% isolatedsource = string, 'yes' or 'no'
% conductivity = vector, conductivity of each compartment
% tempdir = string, allows you to specify the path for the tempory files (default is automatic)
% tempname = string, allows you to specify the full tempory name including path (default is automatic)
% checkmesh = 'yes' or 'no'
%
% See also FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD
% Copyright (C) 2012, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL
%
% 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$
ft_hastoolbox('dipoli', 1);
% get the optional arguments
isolatedsource = ft_getopt(varargin, 'isolatedsource');
conductivity = ft_getopt(varargin, 'conductivity');
tdir = ft_getopt(varargin, 'tempdir');
tname = ft_getopt(varargin, 'tempname');
checkmesh = ft_getopt(varargin, 'checkmesh', 'yes');
% convert to Boolean value
checkmesh = istrue(checkmesh);
if isfield(mesh, 'bnd')
mesh = mesh.bnd;
end
% replace pnt with pos
mesh = fixpos(mesh);
% start with an empty volume conductor
headmodel = [];
% determine the number of compartments
numboundaries = numel(mesh);
% Dipoli expects surface normals to point inwards;
% this checks and corrects if needed
for i=1:numboundaries
switch surface_orientation(mesh(i).pos, mesh(i).tri)
case 'outward'
ft_warning('flipping mesh %d', i);
mesh(i).tri = fliplr(mesh(i).tri);
case 'inward'
% this is ok
case 'otherwise'
ft_error('incorrect mesh %d', i)
end
end
if isempty(isolatedsource)
if numboundaries>1
% the isolated source compartment is by default the most inner one
isolatedsource = true;
else
isolatedsource = false;
end
else
% convert into a boolean
isolatedsource = istrue(isolatedsource);
end
if isolatedsource
fprintf('using isolated source approach\n');
else
fprintf('not using isolated source approach\n');
end
% determine the desired nesting of the compartments
order = surface_nesting(mesh, 'outsidefirst');
% reorder the boundaries
if numel(mesh)>1
fprintf('reordering the boundaries to: ');
fprintf('%d ', order);
fprintf('\n');
% update the order of the compartments
mesh = mesh(order);
end
if isempty(conductivity)
ft_warning('no conductivity specified, using default values')
if numboundaries == 1
% uniform
conductivity = 1;
elseif numboundaries == 3
% skin/skull/brain
conductivity = [0.3300 0.0042 0.3300];
elseif numboundaries == 4
% skin/skull/csf/brain
conductivity = [0.3300 0.0042 1.7900 0.3300];
else
ft_error('conductivity values must be specified')
end
headmodel.cond = conductivity;
else
if numel(conductivity)~=numboundaries
ft_error('each compartment should have a conductivity value');
end
% reorder the user-specified conductivities
headmodel.cond = conductivity(order);
end
% do some sanity checks on the meshes
if checkmesh
for i=1:numboundaries
ntri = size(mesh(i).tri,1);
npos = size(mesh(i).pos,1);
assert(2*(npos-2)==ntri, 'the number of triangles does not match the number of vertices')
end
% check for each of the vertices that it falls inside the previous surface
for i=2:numboundaries
for j=1:size(mesh(i).pos,1)
inside = surface_inside(mesh(i).pos(j,:), mesh(i-1).pos, mesh(i-1).tri);
assert(inside==true, 'vertex %d of surface %d is outside surface %d', j, i, i-1);
end
end
ft_info('the meshes are closed and properly nested')
end % if checkmesh
headmodel.bnd = mesh;
headmodel.skin_surface = 1;
headmodel.source = numboundaries; % this is now the last one
% just to be sure
clear mesh
if isolatedsource
fprintf('using compartment %d for the isolated source approach\n', headmodel.source);
else
fprintf('not using the isolated source approach\n');
end
% find the location of the dipoli binary
str = which('dipoli.exe');
[p, f, x] = fileparts(str);
dipoli = fullfile(p, f); % without the .exe extension
% the following determines the executable to use, which is in all cases the 32-bit version
switch lower(computer)
case {'maci' 'maci64', 'maca64'}
% apple computer with 32- or 64-bit Intel processor, or Apple ARM processor (M1 or M2)
dipoli = [dipoli '.universal'];
case {'glnx86' 'glnxa64'}
% linux computer
dipoli = [dipoli '.glnx86'];
case {'win32', 'win64'}
% windows computer
dipoli = [dipoli '.exe'];
otherwise
ft_error('there is no dipoli executable for your platform');
end
fprintf('using the executable "%s"\n', dipoli);
% determine the prefix for the temporary files
if isempty(tname)
if isempty(tdir)
if ismac
prefix = tempname('/tmp'); % otherwise the filename with the complete path gets too long
else
prefix = tempname;
end
else
prefix = tempname(tdir);
end
else
prefix = tname;
end
% write the triangulations to file
bndfile = cell(1,numboundaries);
for i=1:numboundaries
bndfile{i} = sprintf('%s_%d.tri', prefix, i);
write_tri(bndfile{i}, headmodel.bnd(i).pos, headmodel.bnd(i).tri);
end
% these will hold the shell script and the inverted system matrix
exefile = [prefix '.sh'];
amafile = [prefix '.ama'];
fid = fopen(exefile, 'w');
fprintf(fid, '#!/bin/sh\n');
fprintf(fid, '\n');
fprintf(fid, '%s -i %s << EOF\n', dipoli, amafile);
for i=1:numboundaries
if isolatedsource && headmodel.source==i
% the isolated potential approach should be applied using this compartment
fprintf(fid, '!%s\n', bndfile{i});
else
fprintf(fid, '%s\n', bndfile{i});
end
fprintf(fid, '%g\n', headmodel.cond(i));
end
fprintf(fid, '\n');
fprintf(fid, '\n');
fprintf(fid, 'EOF\n');
fclose(fid);
% ensure that the temporary shell script can be executed
dos(sprintf('chmod +x %s', exefile));
try
% execute dipoli and read the resulting file
dos(exefile);
ama = loadama(amafile);
headmodel = ama2headmodel(ama);
catch
ft_error('an error ocurred while running the dipoli executable - please look at the screen output');
end
% delete the temporary files
for i=1:numboundaries
delete(bndfile{i})
end
delete(amafile);
delete(exefile);
% maintain the general FieldTrip convention of outward-oriented surfaces
% (also for visualization). The calculation of the head model was done with
% inward-oriented surfaces (sometimes flipped from the original input).
for i=1:numboundaries
headmodel.bnd(i).tri = fliplr(headmodel.bnd(i).tri);
end
% remember that it is a dipoli model
headmodel.type = 'dipoli';