-
Notifications
You must be signed in to change notification settings - Fork 2
/
bst_headmodeler.m
2666 lines (2219 loc) · 109 KB
/
bst_headmodeler.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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function [varargout] = bst_headmodeler(varargin);
%BST_HEADMODELER - Solution to the MEG/EEG forward problem
% function [varargout] = bst_headmodeler(varargin);
% Authorized syntax:
% [G, OPTIONS] = bst_headmodeler(StudyFile, OPTIONS);
% [G, OPTIONS] = bst_headmodeler(OPTIONS);
% [OPTIONS] = bst_headmodeler;
%
% --------------------------------- INPUTS -------------------------------------
% INPUTS
%
% [OPTIONS] = bst_headmodeler; Returns the default values for the OPTIONS
% parameter structure
%
% StudyFile: the name of a BrainStorm study file. We suppose the
% corresponding BrainStorm channel file is available in the same folder
% with the conventional file name (e.g., for a study file calles
% meg_brainstormstudy.mat, BST_HEADMODELER expects to find the
% corresponding channel file under the name meg_channel.mat). If the
% channel file were not sitting in the same folder as the study file,
% OPTIONS.ChannelFile enforces computation of the forward model with the
% channel information contained in OPTIONS.ChannelFile.
%
% If no further input arguments are specified, forward modeling is
% completed with the default parameters specified below
%
% OPTIONS a structure where optional parameters may be specified using the
% following fields Note: if no options are specified, BST_HEADMODELER will
% proceed to the computation of the foward problem on a 3D grid of source
% locations that cover the entire head volume See
% OPTIONS.VolumeSourceGridSpacing for default settings
%
% Important notice: there is no need to define all the following fields
% when using the OPTIONS argument. The undefined field(s) will be assigned
% default values.
%
% *
% * Fields Related to forward approach
% *
%
% .Method: is either a character string or a cell array of two strings that
% specifies the kind of approach to be applied to the compuation
% of the foward model In case Method is a cell array, it should
% contain 2 strings, one to specifiy the method to be used for MEG
% and the other for . If only a single string is specified, the
% foward computation will be completed on the set of corresponding
% CHANndx only (i.e. MEG or MEG) Available forward modeling
% methods and corresponding authorized strings for Method:
% - MEG
% 'meg_sphere' (DEFAULT) : Spherical head model designed
% following the Sarvas analytical formulation (i.e.
% considering the true orientation
% of the magnetic field sensors) (see OPTIONS.HeadCenter)
% 'meg_os' : MEG overlapping sphere forward model
% 'meg_bem' : Apply BEM computation (see OPTIONS.BEM for details)
% -
% 'eeg_sphere' : Single-sphere forward modeling (see
% OPTIONS.HeadCenter, OPTIONS.Radii, OPTIONS.Conductivity)
% 'eeg_3sphere' : EEG forward modeling with a set of 3
% concentric spheres (Scalp, Skull, Brain/CSF) (see
% OPTIONS.HeadCenter, OPTIONS.Radii, OPTIONS.Conductivity)
% 'eeg_3sphereBerg' (DEFAULT) : Same as eeg_3sphere with
% correction for possible dipoles outside the sphere
% 'eeg_os' : EEG overlapping sphere head model (see
% OPTIONS.HeadCenter, OPTIONS.Radii, OPTIONS.Conductivity)
% 'eeg_bem' : Apply BEM computation (see OPIONS.BEM for details)
%
% Default is {'meg_sphere','eeg_3sphereBerg'};
%
% .HeadModelName : a character string that specifies the name of the
% headmodel represented by this file, e.g "Spherical",
% "Overlapping Spheres", "Constant Collocation BEM", etc.
% Default is "Default", meaning it will include the the
% name(s) of the method(s) used in the MEG and/or EEG
% forward models
%
% *
% * Fields Related to function's I/O
% *
%
% .HeadModelFile : Specifies the name of the head model file where to store
% the forward model. If set to 'default', the default
% nomenclature for BrainStorm's head model file name is
% used and BST_HEADMODELER creates a file in StudyFile's
% folder.
% Default is empty.
% .ImageGridFile : Specifies the name of the file where to store the full
% cortical gain matrix file If set to 'default', the
% default nomenclature for BrainStorm's head model file
% name is used and BST_HEADMODELER creates a file in
% StudyFile's folder.
% Default is empty.
% .ImageGridBlockSize : Number of sources for which to compute the forward
% model at a time in a block computation routines
% (saves memory space). This option is relevant only
% when some forward modeling on cortical surface is
% requested (i.e. when .Cortex is specified)
% Default is 2000
% .FileNamePrefix : A string that specifies the prefix for all file names
% (Channel, HeadModel, Gain Matrices) when .HeadModelFile
% is set to 'default' and .ChannelFile is empty.
% Default is 'bst_'
% .Verbose : Toggles verbose mode on/off;
% Default is 1 (on)
%
% *
% * Fields Related to Head Geometry *
% *
%
% .Scalp : A structure specifying the Scalp surface envelope to serve
% for parameter adjustment of best-fitting sphere, with
% following fields:
% .FileName : A string specifying the name of the BrainStorm
% tessellation file containing the Scalp
% tessellation (default is 1);
% .iGrid : An integer for the index of the Scalp surface in
% the Faces, Vertices and Comments cell arrays in
% the tessellation file
% Default is empty (Best-fitting sphere is
% computed from the sensor array).
% .HeadCenter: a 3-element vector specifying the coordinates, in the
% sensors coordinate system, of the center of the spheres that
% might be used in the head model.
% Default is estimated from the center of the best-fitting
% sphere to the sensor locations
% .Radii : a 3-element vector containing the radii of the single or 3
% concentric spheres, when needed;
% Order must be the following : [Rcsf, Routerskull, Rscalp];
% Default is estimated from the best-fitting sphere to the
% sensor locations and OPTIONS.Radii is set to: Rscalp [.88
% .93 1]. Rscalp is estimated from the radius of the
% best-fitting sphere;
% .Conductivity : a 3-element vector containing the values for the
% conductivity of the tissues in the following order:
% [Ccsf, Cskull, Cscalp];
% Default is set to [.33 .0042 .33];
% .EEGRef : the NAME (not index of the channel file) of the electrode
% that acts as the reference channel for the EEG. If data is
% referenced to instantaneous average (i.e. so called
% average-reference recording) value is 'AVERAGE REF';
% IMPORTANT NOTICE: When user calls bst_headmodeler with the
% .ChannelLoc option and .ChannelType = 'EEG'and wants the EEG
% reference to be e.g. channel 26, then .EEGRef should be set
% to 'EEG 26'
% Default is 'AVERAGE REF'.
%
% .OS_ComputeParam : if 1, force computation of all sphere parameters when
% choosing a method based on either the MEG or EEG
% overlapping-sphere technique, if 0 and when
% .HeadModelFile is specified, sphere parameters are
% loaded from the pre-computed HeadModel file.
% Default is 1.
%
% .BEM : Structure that specifies the necessary BEM parameters
% .Interpolative : Flag indicating whether exact or
% interpolative approach is used to compute
% the forward solution using BEM.
% if set to 1, exact computation is run on a
% set of points distributed wihtin the inner
% head volume and any subsequent request for
% a forward gain vector (e.g. during a
% volumic source scan using RAP-MUSIC) is
% computed using an interpolation of the
% forward gain vectors of the closest volumic
% grid points. This allows faster computation
% of the BEM solution during source search.
% if set to 0, exact computation is required
% at every source location.
% We recommend to set it to 0 (default) when
% sources have fixed location, e.g.
% constrained on the cortical surface.
% .EnvelopeNames : a cell array of strutures that specifies
% the ORDERED tessellated surfaces to be
% included in the BEM computation.
% .EnvelopeNames{k}.TessFile : A string for
% the name of the tessellation file
% containing the kth surface
% .EnvelopeNames{k}.TessName : A string for
% the name of the surface within the
% tessellation file This string should match
% one of the Comment strings in the
% tessellation file. The chosen surfaces must
% be ordered starting by the the innermost
% surface (e.g. brain or inner skull
% surface) and finishing with the outermost
% layer (e.g. the scalp)
% .Basis : set to either 'constant' or 'linear' (default)
% .Test : set to either 'Galerkin' or 'Collocation' (default)
% .ISA : Isolated-skull approach set to 0 or 1 (default is 1)
% .NVertMax : Maximum number of vertices per envelope,
% therefore leading to decimation of orginal
% surfaces if necessary
% (default is 1000)
% .ForceXferComputation: if set to 1, force recomputation of
% existing transfer matrices in current
% study folder (replace existing
% files);
% Default is 1
%
% *
% * Fields Related to Sensor Information *
% *
%
% .ChannelFile : Specifies the name of the file containing the channel
% information (needs to be a BrainStorm channel file). If
% file does not exists and sensor information is provided in
% .ChannelLoc, a BrainStorm Channl file with name
% .ChannelFile is created
% Default is left blank as this information is extracted
% from the channel file associated to the chosen BrainStorm
% studyfile.
% .Channel : A full BrainStorm channel structure if no channel file is
% specified;
% Default is empty
% .ChannelType : A string specifying the type of channel in ChannelLoc. Can
% be either 'MEG' or 'EEG'. Note that the same channel type
% is assumed for every channel.
% Default is empty
% .ChannelLoc : Specifies the location of the channels where to compute
% the forward model Can be either a 3xNsens (for EEG or
% MEG-magnetometer) or 6xNsens matrix (for the
% MEG-gradiometer case). (for magnetometer or gradiometer
% MEG - channel weights are set to -1 and 1 for each
% magnetometer in the gradiometer respectively) Note that
% in the MEG-gradiometer case, the 3 first (res. last) rows
% of .ChannelLoc stand for each of the magnetometers of the
% gradiometer set. In the case of a mixture of MEG magneto
% and MEG gradio-meters, .ChannelLoc needs to be a 6xNsens
% matrix where the last 3 rows are filled with NaN for
% MEG-magnetometers. If ones wants to handle both EEG and MEG
% sensors, please create a full ChannelFile and use the
% .ChannelFile option.
% Default is empty (Information extracted from the ChannelFile).
% .ChannelOrient : Specifies the orientation of the channels where to
% compute the EEG or MEG forward model Can be either a
% 3xNsens (for EEG and magnetometer MEG) or 6xNsens (for
% gradiometer MEG) matrix or a cell array of such matrices
% (one cell per type of method selected)
% Default is empty (Information extracted from the
% ChannelFile or assume radial orientation when
% .ChannelLoc is filled).
%
% *
% * Fields Related to Source Models *
% *
% .SourceModel : A vector indicating the type of source models to be
% computed; The following code is enfoced:
% -1 : Compute the forward fields of Current Dipole sources
% (available for all forward approaches)
% 1 : 1st-order Current Multipole Sources
% (available for sphere-based MEG approaches only)
% User can set OPTIONS.SourceModel to e.g., [-1 1] to
% compute forward models from both source models.
% Default is -1
%
%
% *
% * Fields Related to Source Localization *
% *
% .Cortex : A structure specifying the Cortex surface
% envelope to serve as an image support with
% following fields.
% .FileName : A string specifying the name of
% the BrainStorm tessellation file
% containing the Cortex
% tessellation;
% .iGrid : An integer for the index of the
% Cortex surface in the Faces,
% Vertices and Comments cell arrays
% in the tessellation file (default
% is 1)
% Default is empty.
% .GridLoc : A 3xNsources matrix that contains the
% locations of the sources at which the forward
% model will be computed
% Default is empty (Information taken from
% OPTIONS.Cortex or OPTIONS.VolumeSourceGrid);
% .GridOrient : a 3xNsources matrix that forces the source
% orientation at every vertex of the .ImageGrid
% cortical surface;
% Defaults is empty; this information being
% extracted from the corresponding tessellated
% surface.
% .ApplyGridOrient : if set to 1, force computation of the forward
% fields by considering the local orientation of
% the cortical surface;
% If set to 0, a set of 3 orthogonal dipoles are
% considered at each vertex location on the
% tessellated surface.
% Default is 1.
% .VolumeSourceGrid : if set to 1, a 3D source grid is designed to
% fit inside the head volume and will serve as a
% source space for scannig techniques such as
% RAP-MUSIC;
% if set to 0, this grid will be computed at the
% first call of e.g. RAP-MUSIC);
% Default is 1
% .VolumeSourceGridSpacing : Spacing in centimeters between two consecutive
% sources in the 3D source grid described above;
% Default is 2 cm.
% .VolumeSourceGridLoc : a 3xN matrix specifying the locations of the
% grid points that will be used to design the
% volumic search grid (see .VolumicSourceGrid)
% Default is empty (locations are estimated
% automatically to cover the estimated inner
% head volume)
% .SourceLoc : a 3xNsources matrix that contains the
% locations of the sources at which the forward
% model will be computed
% Default is empty (Information taken from
% OPTIONS.ImageGrid or
% OPTIONS.VolumeSourceGrid);
% .SourceOrient : a 3xNsources matrix that contains the
% orientations of the sources at which the
% forward model will be computed
% Default is empty.
%
%
% --------------------------------- OUTPUTS ------------------------------------
% OUTPUT
% G if the number of sources (Nsources) if less than
% .ImageGridBlockSize then G is a gain matrix of dimension
% Nsensors x Nsources: Each column of G is the forward field
% created by a dipolar source of unit amplitude. Otherwise, G is
% the name of the binary file containing the gain matrix. This
% file can be read using the READ_GAIN function.
% OPTIONS Returns the OPTIONS structure with updated fields following the
% call to BST_HEADMODELER. Can be useful to obtain a full
% BrainStorm Channel structure when only the .ChannelLoc and
% possibly .ChannelOrient fields were provided.
%<autobegin> ---------------------- 27-Jun-2005 10:43:31 -----------------------
% ------ Automatically Generated Comments Block Using AUTO_COMMENTS_PRE7 -------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\bem_gain.m
% toolbox\bem_xfer.m
% toolbox\berg.m
% toolbox\bst_message_window.m
% toolbox\colnorm.m
% toolbox\get_channel.m
% toolbox\get_user_directory.m
% toolbox\good_channel.m
% toolbox\gridmaker.m
% toolbox\inorcol.m
% toolbox\norlig.m
% toolbox\overlapping_sphere.m
% toolbox\rownorm.m
% toolbox\save_fieldnames.m
% toolbox\source_grids.m
% toolbox\view_surface.m
%
% Subfunctions in this file, in order of occurrence in file:
% BEMGaingridFname = bem_GainGrid(DataType,OPTIONS,BEMChanNdx)
% g = gterm_constant(r,rq)
%
% At Check-in: $Author: Silvin $ $Revision: 68 $ $Date: 12/15/05 4:14a $
%
% This software is part of BrainStorm Toolbox Version 27-June-2005
%
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
% University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
% Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
% CNRS, Hopital de la Salpetriere, Paris, France
%
% See BrainStorm website at http://neuroimage.usc.edu for further information.
%
% Copyright (c) 2005 BrainStorm by the University of Southern California
% This software distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 27-Jun-2005 10:43:31 -----------------------
% /---Script Author--------------------------------------\
% | |
% | *** Sylvain Baillet Ph.D. |
% | Cognitive Neuroscience & Brain Imaging Laboratory |
% | CNRS UPR640 - LENA |
% | Hopital de la Salpetriere, Paris, France |
% | sylvain.baillet@chups.jussieu.fr |
% | |
% \------------------------------------------------------/
%
% Date of creation: March 2002
% ----------------------------- Script History ---------------------------------
%
% SB 07-Aug-2002 Fixed GridLoc saving issues - now GridLoc
% is a cell array of length one which indicates
% the name of the tessellation file that was
% used as an imaging support.
% SB 03-Sep-2002 Fixed minor bugs when function is called from command line
% Now bst_headmodeler can be called even when Brainstorm is not running
% SB 05-Sep-2002 BEM can now be called from command line as well.
% Updated script header
% SB 06-Sep-2002 Updated EEG BEM : added option for arbitrary reference of the
% potentials (takes now OPTIONS.EEGRef into account).
% SB 18-Nov-2002 Updated EEG reference management.
% SB 21-Nov-2002 Minor fixes for command line calls.
% SB 21-Nov-2002bis Fixed bug from the get(*,'VertexNormals') command that
% sometimes returns [0 0 0] at some vertex
% SB 23-Dec-2002 Correct HeadModelFile name is passed to OPTIONS output argument
% JCM 27-May-2004 Cleaning comments
% Found Matlab bug and temp solution.
% Matlab 6.5.0 bug, Technical Solution Number: 1-19NOK
% http://www.mathworks.com/support/solutions/data/1-19NOK.html?solution=1-19NOK
% Must REWIND first before fseek.
% Apparently fixed in 6.5.1
% Example:
% status = fseek(fdest,0,'bof');
% status = fseek(fdest,offset,'bof');
% SB 08-Mar-2005 Fixed bug in BEM computation
% SB 02-Feb-2005 Fixed bug that occured when ChannelLoc field is used to
% passed EEG channel locations (thanks Jeremie Mattout at the FIL)
% ----------------------------- Script History ---------------------------------
% Default options settings--------------------------------------------------------------------------------------------------
DefaultMethod = {'meg_sphere','eeg_3sphereBerg'};
ReducePatchScalpNVerts = 500; % Number of vertices the scalp tessellation will be reduced to in OS computations
BEM_defaults = struct(...
'Basis','linear',...
'Test','galerkin',...
'Interpolative',0,...
'ISA',1,...
'NVertMax',1000,...
'ForceXferComputation', 1, ...
'checksurf',0);
Def_OPTIONS = struct(...
'ApplyGridOrient',1,...
'BEM', BEM_defaults,...
'Channel', [],...
'ChannelFile', '',...
'ChannelLoc', '',...
'ChannelOrient', '',...
'ChannelType', '',...
'Conductivity', [.33 .0042 .33],...
'Cortex',[],...
'EEGRef','',...
'HeadCenter',[],...
'HeadModelFile', '',...
'HeadModelName','Default',...
'ImageGridBlockSize', 2000,...
'ImageGridFile', '',...
'GridOrient',[],...
'Method', {DefaultMethod},...
'OS_ComputeParam', 1,...
'PrefixFileName','bst_',...
'Radii', [],...
'Scalp',[],...
'SourceLoc',[],...
'SourceModel', [-1],...
'SourceOrient',[],...
'StudyFile','',...
'TessellationFile','',...
'VolumeSourceGrid',1,...
'VolumeSourceGridSpacing', 2,...
'VolumeSourceGridLoc', [],...
'Verbose', 1 ...
);
SourceOrderString = {'Current Dipole',...
'Unvalid Order (was Magnetic Dipole Moment)',...
'Current Multipole Expansion',...
'Current Dipole Pairs',...
'Unvalid Order (was Magnetic Dipole Moment PAIRS)',...
'Current Multipole Expansion Pairs'};
if nargin == 0
if nargout > 1
varargout{1} = Def_OPTIONS;
varargout{2} = Def_OPTIONS;
else
varargout{1} = Def_OPTIONS;
end
return
elseif nargin == 1
if ischar(varargin{1}) % User enters only the studyfile name
StudyFile = varargin{1};
OPTIONS = Def_OPTIONS;
elseif isstruct(varargin{1}) % User enters only the OPTION field
OPTIONS = varargin{1};
else
errordlg('Uncorrect input parameter type, please check the help file'), varargout = cell(nargout,1);
return,
end
elseif nargin == 2 % No options were specified, apply default
StudyFile = varargin{1};
OPTIONS = varargin{2};
else
errordlg('Wrong number of arguments when calling head modeler')
varargout = cell(nargout,1);
return
end
% Check field names of passed OPTIONS and fill missing ones with default values
DefFieldNames = fieldnames(Def_OPTIONS);
for k = 1:length(DefFieldNames)
if ~isfield(OPTIONS,DefFieldNames{k}) | strcmp(DefFieldNames{k},'BEM')
if ~isfield(OPTIONS,DefFieldNames{k})
OPTIONS = setfield(OPTIONS,DefFieldNames{k},getfield(Def_OPTIONS,DefFieldNames{k}));
elseif strcmp(DefFieldNames{k},'BEM')
BEM_DefFieldNames = fieldnames(BEM_defaults);
for kk = 1:length(BEM_DefFieldNames)
if ~isfield(OPTIONS.BEM,BEM_DefFieldNames{kk})
OPTIONS.BEM = setfield(OPTIONS.BEM,BEM_DefFieldNames{kk},getfield(BEM_defaults,BEM_DefFieldNames{kk}));
end
end
end
end
end
if isempty(OPTIONS.Conductivity)
OPTIONS.Conductivity = Def_OPTIONS.Conductivity;
end
clear Def_OPTIONS
if isempty(OPTIONS.HeadModelFile) & ~isempty(OPTIONS.ImageGridFile)
% Force creation of a headmodel file
OPTIONS.HeadModelFile = 'default';
end
OPTIONS.HeadModelFileOld = OPTIONS.HeadModelFile;
% What type of forward model (MEG and/or EEG) ?
DataType.MEG = strmatch('meg',OPTIONS.Method); % Rank of the respective forward method in cell array Method (ie could be Method = {'meg_bem','eeg_sphere'} or vice-versa)
DataType.EEG = strmatch('eeg',OPTIONS.Method);
if ~iscell(OPTIONS.Method)
OPTIONS.Method = {OPTIONS.Method};
end
MegMethod = [];
if ~isempty(DataType.MEG)
MegMethod = OPTIONS.Method{DataType.MEG}; % String indicating the forward method selected for MEG (res. EEG)
end
EegMethod = [];
if ~isempty(DataType.EEG)
EegMethod = OPTIONS.Method{DataType.EEG};
end
% Check inputs integrity
%Source models
if ~isempty(find(OPTIONS.SourceModel == 0)) | ~isempty(find(abs(OPTIONS.SourceModel) > 1)) % unValid source models
if ~isempty(DataType.MEG)
errordlg('Valid source model orders for MEG are: -1 (Current Dipole) and 1 (fist-order Current Multipole Expansion)')
end
if ~isempty(DataType.EEG)
errordlg('Valid source model order for EEG is: -1 (Current Dipole) ')
end
varargout = cell(nargout,1);
return
end
% Source locations
if isempty(OPTIONS.SourceLoc) & ~OPTIONS.VolumeSourceGrid & isempty(OPTIONS.Cortex) % No source locations were specified
errordlg('No source locations are specified. Please fill either one of the following fields of OPTIONS: .SourceLoc / .VolumeSourceGrid / .Cortex')
varargout = cell(nargout,1);
return
end
%--------------------------------------------------------------------------------------------------------------------------------------------
%
% HEAD MODELING BEGINS
%
%--------------------------------------------------------------------------------------------------------------------------------------------
if OPTIONS.Verbose,
clockk = fix(clock);
time0 = clock;
bst_message_window({...
' ',...
'__________________________________________',...
sprintf(...
'Head modeling begins (%s - %dH %dmin %ds)',date,clockk(4:6))...
})
end
User = get_user_directory;
if isempty(User) % Function is called from command line: BrainStorm is not running
% Use default folders
User.SUBJECTS = pwd;
User.STUDIES = pwd;
end
try
cd(User.STUDIES)
catch
end
% Get all Subject and Study information------------------------------------------------------------------------------------------------------
if exist('StudyFile','var')
if OPTIONS.Verbose, bst_message_window('Loading Study and Subject Information...'), end
try
load(StudyFile) % Load study information
[StudyPath,tmp,ext] = fileparts(StudyFile); clear tmp
catch
errordlg(['Could not find ',StudyFile,' in current Matlab path'])
varargout = cell(nargout,1);
return
end
if isempty(OPTIONS.StudyFile)
OPTIONS.StudyFile = StudyFile;
end
SubjectFile = BrainStormSubject; clear BrainStormSubject
OPTIONS.rooot = findstr(StudyFile,'brainstormstudy.mat');
if isempty(OPTIONS.rooot)
errordlg('Study file name should be of the form ''*brainstormstudy.mat''')
varargout = cell(nargout,1);
return
end
OPTIONS.rooot = strrep(StudyFile,'brainstormstudy.mat','');
Study = load(StudyFile);
%User = get_user_directory;
try
cd(User.SUBJECTS)
Subject = load(Study.BrainStormSubject);
catch
errordlg(sprintf('Please make sure the subject''s information in %s is available on this plateform',Study.BrainStormSubject)); return
end
if isempty(OPTIONS.TessellationFile) % No Tessellation file was passed to the headmodeler
if isfield(Subject,'Tesselation')
OPTIONS.TessellationFile = Subject.Tesselation; % Take default (OBSOLETE as for BsT MMII because we now consider all tessellation files in Subject's folder)
else
OPTIONS.TessellationFile = '';
%[OPTIONS.TessellationFile,DataPopup, Leader] = find_brainstorm_files('tess',fullfile(Users.SUBJECTS,OPTIONS.subjectpath));
end
end
end
% Get Channel Information ------------------------------------------------------------------------------------------------------------------
if ~isempty(OPTIONS.Channel)
Channel = OPTIONS.Channel;
else
if isempty(OPTIONS.ChannelFile) & isempty(OPTIONS.ChannelLoc) % Load Channel file in current folder
[Channel, ChannelFile] = get_channel(fullfile(User.STUDIES,OPTIONS.StudyFile));
OPTIONS.ChannelFile = fullfile(fileparts(fullfile(User.STUDIES,OPTIONS.StudyFile)),ChannelFile);
if isempty(ChannelFile)
errordlg(sprintf('Channel file %s is missing. Please have it available it the current study folder.',OPTIONS.ChannelFile), 'Missing Channel File')
return
end
elseif isempty(OPTIONS.ChannelLoc) & exist(OPTIONS.ChannelFile,'file') % If no specific channel locations are given and channel file exists, load the proper channel file
OPTIONS.rooot = strrep(lower(OPTIONS.ChannelFile),'channel.mat','');
try
load(OPTIONS.ChannelFile)
catch
cd(User.STUDIES)
load(OPTIONS.ChannelFile)
end
else % Create a dummy Channel structure with Channel Locations (res. Orientations) specified in OPTIONS.ChannelLoc (res .ChannelOrient)
if OPTIONS.Verbose, bst_message_window('Creating the channel structure from information in the OPTIONS fields. . .'), end
% Get Channel Locations
nchan = size(OPTIONS.ChannelLoc,2); % Number of channels
Channel = struct('Loc',[],'Orient',[],'Comment','','Weight',[],'Type','','Name','');
Channel(1:nchan) = deal(Channel);
ChanType = upper(OPTIONS.ChannelType);
if isempty(ChanType)
errordlg('Please specify a channel type (i.e. MEG or EEG) in the ChannelType field of OPTIONS'),
bst_message_window('Please specify a channel type (i.e. MEG or EEG) in the ChannelType field of OPTIONS'),
varargout = cell(nargout,1);
return
end
[Channel(:).Type] = deal(ChanType);
if size(OPTIONS.ChannelLoc,1) == 6 % MEG Gradiometers or mixture gradio/magneto meters
OPTIONS.ChannelLoc = reshape(OPTIONS.ChannelLoc,3,nchan*2);
iGradFlag = 1; % Flag - gradiometer-type sensor set
else
iGradFlag = 0; % Flag - EEG or magnetometer-type sensor set
end
ichan = 1;
for k = 1:nchan
if iGradFlag
Channel(k).Loc = OPTIONS.ChannelLoc(:,ichan:ichan+1);
ichan = ichan+2;
else
if strcmp(ChanType,'MEG')
Channel(k).Loc = OPTIONS.ChannelLoc(:,ichan);
%elseif strcmp(ChanType,'EEG') % Artificially add a dummy column full of zeros to each Channel(k).Loc
%Channel(k).Loc = [OPTIONS.ChannelLoc(:,ichan) [0 0 0]'];
elseif strcmp(ChanType,'EEG') % Artificially add a dummy column full of zeros to each Channel(k).Loc
Channel(k).Loc = [OPTIONS.ChannelLoc(:,ichan)];
end
ichan = ichan+1;
end
Channel(k).Name = sprintf('%s %d',ChanType,k);
Channel(k).Comment = int2str(k);
end
clear ichan k
% Get Channel Orientations
if isempty(OPTIONS.ChannelOrient) & strcmp(ChanType,'MEG') % No channel orientation were specified: use radial sensors
if OPTIONS.Verbose, bst_message_window('Assign radial orientation to all Channels. . .'), end
if isempty(OPTIONS.HeadCenter)
if iGradFlag
Vertices = OPTIONS.ChannelLoc(:,1:2:end)';
else
Vertices = OPTIONS.ChannelLoc';
end
nscalp = size(Vertices,1);
if nscalp > 500 % 500 points is more than enough to compute scalp's best fitting sphere
Vertices = Vertices(unique(round(linspace(1,nscalp,500))),:);
nscalp = size(Vertices,1);
end
nmes = size(Vertices,1);
% Run parameters fit --------------------------------------------------------------------------------------------------------------
mass = mean(Vertices); % center of mass of the scalp vertex locations
R0 = mean(norlig(Vertices - ones(nscalp,1)*mass)); % Average distance between the center of mass and the scalp points
vec0 = [mass,R0];
[minn,brp] = fminsearch('dist_sph',vec0,[],Vertices);
OPTIONS.HeadCenter = minn(1:end-1);
if isempty(OPTIONS.Radii)
OPTIONS.Radii = minn(end);
OPTIONS.Radii = minn(end)*[1/1.14 1/1.08 1];
end
if OPTIONS.Verbose, bst_message_window({...
sprintf('Center of the Sphere : %3.1f %3.1f %3.1f (cm)',100*OPTIONS.HeadCenter),...
sprintf('Radius : %3.1f (cm)',100*OPTIONS.Radii(3)),...
'-> DONE',' '})
end
clear minn brp mass R0 Vertices vec0
end
tmp = [Channel.Loc] - repmat(OPTIONS.HeadCenter',1,nchan*(1+(iGradFlag==1)));
OPTIONS.ChannelOrient = tmp*inorcol(tmp); % Radial orientation for every channel
clear tmp
elseif ~isempty(OPTIONS.ChannelOrient) & strcmp(ChanType,'MEG') & iGradFlag
OPTIONS.ChannelOrient = reshape(OPTIONS.ChannelOrient,3,nchan*2);
end
if strcmp(ChanType,'MEG')
for k=1:nchan
Channel(k).Orient = OPTIONS.ChannelOrient(:,((1+(iGradFlag==1))*k-1):(1+(iGradFlag==1))*k);
end
end
clear k
% Define Weights
if iGradFlag
[Channel(:).Weight] = deal([1 -1]);
else
[Channel(:).Weight] = deal([1]);
end
if strcmp(ChanType,'MEG') % Force no reference channels
Channel(1).irefsens = [];
end
% New Channel stbemructure completed: save as a new channel file
if ~isempty(OPTIONS.ChannelFile)
%OPTIONS.ChannelFile = [OPTIONS.rooot,'channel.mat'];
save(OPTIONS.ChannelFile,'Channel')
if OPTIONS.Verbose, bst_message_window({...
sprintf('Channel file created in : %s',OPTIONS.ChannelFile),...
'-> DONE',' '}), end
end
end
%load(OPTIONS.ChannelFile)
OPTIONS.Channel = Channel;
end
% Find MEG and EEG channel indices
MEGndx = good_channel(Channel,[],'MEG');
EEGndx = good_channel(Channel,[],'EEG');
EEGREFndx = good_channel(Channel,[],'EEG REF');
MEG = ~isempty(MEGndx) & ~isempty(DataType.MEG); % == 1 if MEG is requested and is available
EEG = ~isempty(EEGndx) & ~isempty(DataType.EEG);
if ~EEG & ~MEG
errordlg('Please check that data (MEG or EEG) and channel types are compatible.')
return
end
% Test whether CME models are requested for EEG
if ~isempty(find(OPTIONS.SourceModel == 1)) & EEG
if MEG
if OPTIONS.Verbose
bst_message_window('wrap',...
{'',...
'CME model not available for EEG - Skipping. . .',...
'Keeping it for MEG forward modeling only',...
''...
})
end
else % EEG only is requested
errordlg(...
{'CME model not available for EEG - Skipping. . .',...
'Keeping it for MEG forward modeling only'...
})
return
end
end
% Detect EEG reference - if none is specified in the .Comment or Type fields,
% and if none was passed through OPTIONS.EEGRef
% -> Apply average-referencing of the potentials by default.
if EEG
if isempty(OPTIONS.EEGRef)
% EEG Reference Channel
EEGREFndx = good_channel(Channel,[],'EEG REF');
if isempty(EEGREFndx) % Average EEG reference anyway
[Channel(EEGndx).Comment] = deal('AVERAGE REF');
end
else % EEG Reference is specified
switch(OPTIONS.EEGRef)
case 'AVERAGE REF'
[Channel(:).Comment] = deal('AVERAGE REF');
otherwise
EEGREFndx = strmatch(OPTIONS.EEGRef,char(Channel(:).Name));
if isempty(EEGREFndx)
errordlg(sprintf(...
'No channel named ''%s'' was found amongst available EEG channels. Cannot use it as a reference for EEG.',OPTIONS.EEGRef...
))
return
end
end
end
% if isempty(EEGREFndx) & ~MEG % No reference electrode was defined for the EEG stop if no MEG to compute...
% errordlg('Please make sure you have defined a reference for the EEG')
% return
% elseif isempty(EEGREFndx) & MEG % Skip EEG forward modeling and proceed to MEG
% EEG = 0;
% if OPTIONS.Verbose
% bst_message_window({...
% 'No reference was defined for the EEG',...
% 'Skipping EEG modeling and completing MEG only...',...
% })
% end
% end
end
if OPTIONS.Verbose,
if EEG & MEG
if ~isempty(EEGREFndx)
bst_message_window({'Channel count for current computation:',[int2str(length(MEGndx)), ' MEG Channels'],...
sprintf('%d EEG Channels with Reference: %s',length(EEGndx),Channel(EEGREFndx).Name),' '})
else
bst_message_window({'Channel count for current computation:',[int2str(length(MEGndx)), ' MEG Channels'],...
sprintf('%d EEG Channels with Average Reference',length(EEGndx)),' '})
end
elseif MEG
bst_message_window({'Channel count for current computation:',[int2str(length(MEGndx)), ' MEG Channels'],...
' '})
elseif EEG
if ~isempty(EEGREFndx)
bst_message_window({'Channel count for current computation:',...
sprintf('%d EEG Channels with Reference: %s',length(EEGndx),Channel(EEGREFndx).Name),' '})
else
bst_message_window({'Channel count for current computation:',...
sprintf('%d EEG Channels with Average Reference',length(EEGndx)),' '})
end
end
end
if MEG & isempty(MEGndx) % User wants to compute MEG but no MEG data is available
errordlg('Sorry - No MEG data is available'), return
end
if EEG & isempty(EEGndx) % User wants to compute EEG but no EEG data is available
errordlg('Sorry - No EEG data is available'), return
end
% Computation of parameters of the best-fitting sphere --------------------------------------------------------------------------------------------------------------
if length(findstr('bem',[OPTIONS.Method{:}])) ~= length(OPTIONS.Method)% Only if sphere-based head model is requested in any modality (MEG or EEG)
% Best-fitting sphere parameters --------------------------------------------------------------------------------------------------------------
if isempty(OPTIONS.HeadCenter) | isempty(OPTIONS.Radii)
if OPTIONS.Verbose, bst_message_window('Estimating Center of the Head. . .'), end
if isempty(OPTIONS.Scalp) % Best-fitting sphere is derived from the sensor array
% ans = questdlg('No scalp surface was selected - Do you want to use a spherical approximation of the head derived from the sensor locations instead ?',...
% '','Yes','No','Yes');
if EEG & length(EEGndx) > 9 % If EEG is available but a reasonable number of electrodes, use it to compute the sphere parameters
ndx = [EEGndx]; % Compute Sphere parameters from EEG sensors only
else
ndx = [MEGndx];
end
Vertices = zeros(length(ndx),3);
for k=1:length(ndx)
Vertices(k,:) = Channel(ndx(k)).Loc(:,1)';
end
else % Best-fitting sphere is computed from the set of vertices of the scalp surface enveloppe
try
load(fullfile(User.SUBJECTS,OPTIONS.Scalp.FileName),'Vertices')
catch
load(OPTIONS.Scalp.FileName,'Vertices')
end
if ~isfield(OPTIONS.Scalp,'iGrid') % Apply Default
OPTIONS.Scalp.iGrid = 1;
end
try
Vertices = Vertices{OPTIONS.Scalp.iGrid}';
catch
errordlg(sprintf(...
'Tessellation file %s does not contain %d enveloppes',...
OPTIONS.Scalp.FileName,OPTIONS.Scalp.iGrid))
end
end
nscalp = size(Vertices,1);
nmes = size(Vertices,1);
% Run parameters fit --------------------------------------------------------------------------------------------------------------
mass = mean(Vertices); % center of mass of the scalp vertex locations
R0 = mean(norlig(Vertices - ones(nscalp,1)*mass)); % Average distance between the center of mass and the scalp points
vec0 = [mass,R0];
[SphereParams,brp] = fminsearch('dist_sph',vec0,[],Vertices);
if OPTIONS.Verbose
bst_message_window({...
sprintf('Center of the Sphere : %3.1f %3.1f %3.1f (cm)',100*SphereParams(1:end-1)'),...
sprintf('Radius : %3.1f (cm)',100*SphereParams(end)),...
'-> DONE',' '})
end
end % no head center and sphere radii were specified by user
% Assign default values for sphere parameters
% if none were specified before
if isempty(OPTIONS.Radii)
OPTIONS.Radii = SphereParams(end)*[1/1.14 1/1.08 1];
end
if isempty(OPTIONS.HeadCenter)
OPTIONS.HeadCenter = SphereParams(1:end-1)';
end
if isempty(OPTIONS.Conductivity)