MEG current phantom (Elekta-Neuromag)

Authors: Ken Taylor, Francois Tadel, John Mosher

This tutorial explains how to import and process Elekta-Neuromag current phantom recordings. We decided to release this example for testing and cross-validation purposes. With these datasets, we can evaluate the equivalence of various forward models and dipole fitting methods in the case of simple recordings with single dipoles. The recordings are available in two file formats (native and FIF) to cross-validate the file readers available in Brainstorm and MNE. A similar page exists for the CTF phantom.

Contents

License
The phantom
Description of the experiment
Download and installation
Generate anatomy
Access the recordings
Organize the events
Import the dipole events
Noise covariance
Average trials
Head modeling
Dipole source estimation
Dipole scanning
Results
Scripting

License

This tutorial dataset remains a property of its authors: Ken Taylor, John Mosher (Epilepsy Center, Cleveland Clinic Neurological Institute, Cleveland, OH USA).

If you reference this dataset in your publications, please acknowledge them and cite Brainstorm as indicated on the website. For questions, please contact us through the forum.

The phantom

A current phantom is provided with the Elekta Neuromag for checking the system performance. It contains 32 artificial dipoles and four fixed head-position indicator coils. The phantom is based on the mathematical fact that an equilateral triangular line current produces equivalent magnetic field distribution to that of a tangential current dipole in a spherical conductor, provided that the vertex of the triangle and the origin of the conducting sphere coincide. For a detailed description of how the phantom works, see here.

The phantom dipoles are energized using an internal signal generator which also feeds the HPI coils. An external multiplexer box is used to connect the signal to the individual dipoles. Only one dipole can be activated at a time. The location of the dipole is recorded relative to the center of the sphere (0,0,0)m, where X is positive toward the nasion, Y is positive toward the left ear and Z is positive toward the top of the head (see the CoordinateSystems tutorial for more details).

Use of the phantom is shown below. Note that the uncovered version is the phantom that came with the Neuromag-122, which explicitly shows the wiring. The covered version uses the same principle but somewhat different dipole locations. Further details are available in Section 7.2 of the User's Manual.

References

R.J. Ilmoniemi, M.S. Hämäläinen, and J. Knuutila, The Forward and Inverse Problems in the Spherical Model. In: Biomagnetism: Applications and Theory, eds. H. Weinberg, G. Stroink, T. Katil, Pergamon Press, 1985.
Elekta Neuromag System Hardware User's Manual, Revision G, September 2005.

Description of the experiment

Files distributed as part of the phantom Elekta tutorial:

kojak_all_2000nAm_pp_no_chpi_no_ms_raw.fif
Source = 2000 nAm. This corresponds to an unrealistically strong 1,000 nA-m (2,000 nA-m peak to peak) dipole that gives the highest SNR visualization of the experimental source.
kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif
Source = 200 nAm. This is a weaker dipole, closer to the range of amplitudes we can expect from inter-ictal spikes due to epilepsy, visible in raw data.
kojak_all_20nAm_pp_no_chpi_no_ms_raw.fif
Source = 20 nAm. This represents some of the weakest sources we expect in evoked studies which require averaging to detect and estimate (i.e. generally cannot be seen in single trials).
All recordings were performed by John Mosher in the MEG lab at the Cleveland Clinic with an Elekta phantom.

Download and installation

Requirements: You have followed the introduction tutorials and Brainstorm is installed.
Go to the Download page and download the file: sample_phantom_elekta.zip
Unzip it in a folder that is not in any of the Brainstorm folders (program folder or database folder)
Start Brainstorm (Matlab scripts or stand-alone version)
Select the menu File > Create new protocol. Name it TutorialPhantomElekta and select:
- No, use individual anatomy,
- No, use one channel file per acquisition run (MEG).

Generate anatomy

In the Matlab command window: type "generate_phantom_elekta".
This creates a new subject Kojak (so named after the 70's TV show) and generates the "anatomy" for this device: one volume and a few surfaces representing the geometry of the phantom.
You can display the MRI and surfaces as presented in the introduction tutorials.

Access the recordings

We can now review one of the raw kojak data sets. These have been generated by sequentially activating each of the 32 phantom dipoles in a single raw file.

Switch to the functional data view by clicking the second button under the protocol name.
Right-click on the subject folder > Review raw file.
- Select the file format: "MEG/EEG: Elekta-Neuromag (*.fif)"
- Select the 200nA file: sample_phantom_elekta/kojak_all_200nAm...fif
Click Event channel, select STI201, and hit OK.
A figure is opened at the end of the import to show the position of the phantom in the MEG:

Organize the events

Double click the link to raw file: it opens the recordings and shows the list of the events.
In the Record tab, rename your groups so that the ordering remains convenient:
- Double-click on event group "1" (or menu Events > Rename group), rename it to 01.
- Repeat this for events #2-#9.
Some unnecessary triggers are generated by the neuromag system, we will delete them:
- Select all the categories with only one event (#256-#7939), press the delete key.
These steps can also be performed from your MATLAB command window:
- Close the file viewer.
- Right-click on the link > File > Export to Matlab > "link"
- Run the lines of code below:
```
for i = 1:32
    STR = str2num(link.F.events(i).label);
    link.F.events(i).label = sprintf('%02.0f',STR);
end
link.F.events(33:end) = [];
```
- Right-click on the link > File > Import from Matlab > "link".
There is also a switching transient in the phantom generator which should be removed:
- Select category 01, click on the first time instant and press the delete key.
Events should now appear as per the image below on the right:
Close the recordings figure to save the modifications.

Import the dipole events

Right click on Link to raw file > Import in database:

Select all the events #01-#32, Epoch time = [-100,300]ms
Uncheck Apply SSP/ICA projectors (we will use noise covariance to stabilize the data)
Check Remove DC offset, time range = [-100,-10]ms
Resample the recordings to 100Hz (change from 1000Hz)
Do not create separate folders
Click Import

Noise covariance

With the epochs imported, we can now calculate the noise covariance from the prestims:

Select all the imported epochs (click on the first group, hold the shift key, click on the last).
Right-click on any epochs group > Noise covariance > Compute from recordings.
Keep the default options, and hit OK.

Average trials

In Process1, select all the 32 trials groups.
Run process Average > Average files > By trial group (folder average).

Head modeling

Right-click on the channel file > Compute head model > MRI volume, Single sphere.
A figure opens showing a sphere fit to the phantom, however we can see that the location is not exactly (0,0,0). We can adjust it manually:
- Click the Edit sphere properties button in the top toolbar.
- Change the center location to [ 0 0 0 ] and hit OK.
Click the [OK] button. A window appears allowing adjustment of the volume source grid.
Select Regular grid (isotropic) with Grid resolution 2.5mm.
This results in an around 56K grid point, which is sufficient.

Dipole source estimation

Right click the Single sphere (volume) and select Compute sources [2016]
Select: Dipole modeling, Unconstrained, MEG GRAD only
Click Show details to show the advanced options, select No covariance regularization
Ignore the other settings, and click OK.

Dipole scanning

In Process1, select the 32 averages and click on [Process sources].
Run process: Sources > Dipole scanning > Time window=[60,60]ms
We would like to group the 32 new dipole files together. Two methods for this:
- Select all the dipoles files from the database explorer (requires that you expand all the average nodes), right-click > Merge dipoles.
- Execute the lines of code below in your Matlab command window:
```
sStudy = bst_get('StudyWithCondition','Kojak/kojak_all_200nAm_pp_no_chpi_no_ms_raw');
DipolesFile = dipoles_merge({sStudy.Dipoles.FileName});
```
As an optional step, we can sign flip the dipoles so that they all point up for convenience:
- Right-click on dipole file Merge: 32 files > File > Export to Matlab > "dip"
- Execute the lines of code below in your Matlab command window:
```
for i = 1:32
    dip.Dipole(i).Amplitude = dip.Dipole(i).Amplitude * ...
    sign(dip.Dipole(i).Amplitude(3));
end
dip.Comment = [dip.Comment ' | flipped'];
```
- Right-click on the folder containing all these files > File > Import from Matlab > "dip".
Merged the flipped dipoles file with the true dipole positions:
- Select "Merge: 32 files | flipped" and "True phantom dipoles" (in folder TrueDipoles)
- Right-click on one of them > Merge dipoles.
- Rename this file into something more meaningful (eg. "comparison").
Right-click on this last file > Display on MRI (3D), or double-click on the file. The true dipole locations are shown in blue (0ms) along with the estimated locations in red (60ms).

Results

 Dipole stats

 Dipole      Loc (mm)      Amp (nA-m)  Gof     Perf       Chi2     RChi2
--------------------------------------------------------------------------
  01 - [ -0.4 -60.4  21.3]    75.3    98.8%    208.6    532 (201)   2.65
  02 - [  2.1 -47.9  23.8]    80.2    97.9%    117.4    303 (201)   1.51
  03 - [ -0.4 -35.4  23.8]    80.1    94.2%     64.5    258 (201)   1.28
  04 - [ -0.4 -25.4  21.3]    75.0    84.3%     36.3    245 (201)   1.22
  05 - [  2.1 -35.4  51.3]    83.8    98.2%    154.4    442 (201)   2.20
  06 - [  2.1 -25.4  43.8]    90.8    96.8%     87.9    258 (201)   1.28
  07 - [  2.1 -12.9  38.8]    94.9    91.9%     53.2    249 (201)   1.24
  08 - [  2.1  -5.4  33.8]    75.7    80.9%     30.6    222 (201)   1.11
  09 - [-60.4  -0.4  21.3]    77.2    98.7%    204.2    560 (201)   2.78
  10 - [-47.9  -0.4  23.8]    82.3    98.1%    120.0    282 (201)   1.40
  11 - [-35.4  -0.4  23.8]    83.2    94.7%     69.0    268 (201)   1.33
  12 - [-25.4  -0.4  21.3]    79.9    88.7%     42.3    227 (201)   1.13
  13 - [-35.4   2.1  51.3]    83.8    98.4%    178.3    513 (201)   2.55
  14 - [-27.9   2.1  46.3]    76.8    96.8%    101.6    336 (201)   1.67
  15 - [-15.4   2.1  38.8]    87.6    93.4%     60.9    261 (201)   1.30
  16 - [ -5.4   2.1  33.8]    73.0    82.9%     35.1    255 (201)   1.27
  17 - [ -0.4  47.1  41.3]    84.9    98.6%    176.4    434 (201)   2.16
  18 - [ -0.4  44.6  31.3]    75.9    98.0%    106.0    229 (201)   1.14
  19 - [ -0.4  39.6  18.8]    79.5    94.4%     68.6    277 (201)   1.38
  20 - [ -0.4  34.6   8.8]    72.1    88.8%     42.3    224 (201)   1.12
  21 - [ -0.4  14.6  61.3]    78.8    98.3%    131.7    305 (201)   1.52
  22 - [  2.1  17.1  48.8]    86.7    96.1%     82.2    272 (201)   1.36
  23 - [  2.1  22.1  36.3]    84.2    92.6%     55.9    251 (201)   1.25
  24 - [  2.1  19.6  23.8]    87.5    83.9%     36.7    259 (201)   1.29
  25 - [ 47.1   2.1  43.8]    80.5    97.1%     94.2    267 (201)   1.33
  26 - [ 42.1  -0.4  31.3]    82.8    93.7%     61.7    255 (201)   1.27
  27 - [ 37.1  -0.4  18.8]    86.2    88.7%     45.6    265 (201)   1.32
  28 - [ 29.6   2.1  13.8]    81.2    87.3%     32.6    155 (201)   0.77
  29 - [ 14.6   2.1  58.8]    91.3    97.6%    129.0    406 (201)   2.02
  30 - [ 17.1   2.1  48.8]    85.6    95.8%     76.7    258 (201)   1.28
  31 - [ 19.6   2.1  38.8]    78.5    91.1%     48.3    227 (201)   1.13
  32 - [ 19.6  -0.4  23.8]    90.2    84.2%     34.2    220 (201)   1.09

 Location errors from true

 Dipole      Loc (mm)            True (mm)         Diff [x y z]     Norm (mm)
------------------------------------------------------------------------------
  01 - [ -0.4 -60.4  21.3] [  0.0 -59.7  22.9] [ -0.4  -0.7  -1.6]   (  1.8)
  02 - [  2.1 -47.9  23.8] [  0.0 -48.6  23.5] [  2.1   0.7   0.3]   (  2.3)
  03 - [ -0.4 -35.4  23.8] [  0.0 -35.8  25.5] [ -0.4   0.4  -1.7]   (  1.8)
  04 - [ -0.4 -25.4  21.3] [  0.0 -24.8  23.1] [ -0.4  -0.6  -1.8]   (  2.0)
  05 - [  2.1 -35.4  51.3] [  0.0 -37.2  52.0] [  2.1   1.8  -0.7]   (  2.9)
  06 - [  2.1 -25.4  43.8] [  0.0 -27.5  46.4] [  2.1   2.1  -2.6]   (  4.0)
  07 - [  2.1 -12.9  38.8] [  0.0 -15.8  41.0] [  2.1   2.9  -2.2]   (  4.2)
  08 - [  2.1  -5.4  33.8] [  0.0  -7.9  33.0] [  2.1   2.5   0.8]   (  3.4)
  09 - [-60.4  -0.4  21.3] [-59.7  -0.0  22.9] [ -0.7  -0.4  -1.6]   (  1.8)
  10 - [-47.9  -0.4  23.8] [-48.6  -0.0  23.5] [  0.7  -0.4   0.3]   (  0.9)
  11 - [-35.4  -0.4  23.8] [-35.8  -0.0  25.5] [  0.4  -0.4  -1.7]   (  1.8)
  12 - [-25.4  -0.4  21.3] [-24.8  -0.0  23.1] [ -0.6  -0.4  -1.8]   (  2.0)
  13 - [-35.4   2.1  51.3] [-37.2  -0.0  52.0] [  1.8   2.1  -0.7]   (  2.9)
  14 - [-27.9   2.1  46.3] [-27.5  -0.0  46.4] [ -0.4   2.1  -0.1]   (  2.1)
  15 - [-15.4   2.1  38.8] [-15.8  -0.0  41.0] [  0.4   2.1  -2.2]   (  3.1)
  16 - [ -5.4   2.1  33.8] [ -7.9  -0.0  33.0] [  2.5   2.1   0.8]   (  3.4)
  17 - [ -0.4  47.1  41.3] [  0.0  46.1  44.4] [ -0.4   1.0  -3.1]   (  3.3)
  18 - [ -0.4  44.6  31.3] [  0.0  41.9  34.0] [ -0.4   2.7  -2.7]   (  3.8)
  19 - [ -0.4  39.6  18.8] [  0.0  38.3  21.6] [ -0.4   1.3  -2.8]   (  3.1)
  20 - [ -0.4  34.6   8.8] [  0.0  31.5  12.7] [ -0.4   3.1  -3.9]   (  5.0)
  21 - [ -0.4  14.6  61.3] [  0.0  13.9  62.4] [ -0.4   0.7  -1.1]   (  1.4)
  22 - [  2.1  17.1  48.8] [  0.0  16.2  51.5] [  2.1   0.9  -2.7]   (  3.6)
  23 - [  2.1  22.1  36.3] [  0.0  20.0  39.1] [  2.1   2.1  -2.8]   (  4.1)
  24 - [  2.1  19.6  23.8] [  0.0  19.3  27.9] [  2.1   0.3  -4.1]   (  4.7)
  25 - [ 47.1   2.1  43.8] [ 46.1  -0.0  44.4] [  1.0   2.1  -0.6]   (  2.4)
  26 - [ 42.1  -0.4  31.3] [ 41.9  -0.0  34.0] [  0.2  -0.4  -2.7]   (  2.8)
  27 - [ 37.1  -0.4  18.8] [ 38.3  -0.0  21.6] [ -1.2  -0.4  -2.8]   (  3.1)
  28 - [ 29.6   2.1  13.8] [ 31.5  -0.0  12.7] [ -1.9   2.1   1.1]   (  3.0)
  29 - [ 14.6   2.1  58.8] [ 13.9  -0.0  62.4] [  0.7   2.1  -3.6]   (  4.2)
  30 - [ 17.1   2.1  48.8] [ 16.2  -0.0  51.5] [  0.9   2.1  -2.7]   (  3.6)
  31 - [ 19.6   2.1  38.8] [ 20.0  -0.0  39.1] [ -0.4   2.1  -0.3]   (  2.1)
  32 - [ 19.6  -0.4  23.8] [ 19.3  -0.0  27.9] [  0.3  -0.4  -4.1]   (  4.2)

 Orientation errors from true

 Dipole  Amp (nA-m)      [X Y Z]              TRUE [X Y Z]     Degrees
------------------------------------------------------------------------
  01 -    ( 75.3)  [-0.01  0.33  0.94] vs [ 0.00  0.36  0.93] (  1.7 )
  02 -    ( 80.2)  [-0.02  0.44  0.90] vs [ 0.00  0.44  0.90] (  1.5 )
  03 -    ( 80.1)  [-0.09  0.56  0.83] vs [ 0.00  0.58  0.81] (  5.4 )
  04 -    ( 75.0)  [-0.09  0.64  0.76] vs [ 0.00  0.68  0.73] (  6.0 )
  05 -    ( 83.8)  [-0.02  0.82  0.57] vs [ 0.00  0.81  0.58] (  1.5 )
  06 -    ( 90.8)  [-0.02  0.86  0.50] vs [ 0.00  0.86  0.51] (  1.1 )
  07 -    ( 94.9)  [-0.01  0.95  0.32] vs [ 0.00  0.93  0.36] (  2.6 )
  08 -    ( 75.7)  [ 0.11  0.98  0.15] vs [ 0.00  0.97  0.23] (  7.8 )
  09 -    ( 77.2)  [ 0.33  0.04  0.94] vs [ 0.36  0.00  0.93] (  2.7 )
  10 -    ( 82.3)  [ 0.44  0.07  0.89] vs [ 0.44  0.00  0.90] (  4.3 )
  11 -    ( 83.2)  [ 0.55  0.10  0.83] vs [ 0.58  0.00  0.81] (  6.1 )
  12 -    ( 79.9)  [ 0.63  0.19  0.75] vs [ 0.68  0.00  0.73] ( 11.2 )
  13 -    ( 83.8)  [ 0.82  0.02  0.57] vs [ 0.81  0.00  0.58] (  1.7 )
  14 -    ( 76.8)  [ 0.86  0.02  0.52] vs [ 0.86  0.00  0.51] (  1.1 )
  15 -    ( 87.6)  [ 0.93  0.10  0.36] vs [ 0.93  0.00  0.36] (  5.8 )
  16 -    ( 73.0)  [ 0.98  0.11  0.15] vs [ 0.97  0.00  0.23] (  8.0 )
  17 -    ( 84.9)  [ 0.00 -0.66  0.75] vs [ 0.00 -0.69  0.72] (  2.7 )
  18 -    ( 75.9)  [-0.05 -0.57  0.82] vs [ 0.00 -0.63  0.78] (  5.1 )
  19 -    ( 79.5)  [-0.14 -0.43  0.89] vs [ 0.00 -0.49  0.87] (  9.2 )
  20 -    ( 72.1)  [-0.15 -0.24  0.96] vs [ 0.00 -0.37  0.93] ( 11.5 )
  21 -    ( 78.8)  [-0.10 -0.97  0.23] vs [ 0.00 -0.98  0.22] (  5.7 )
  22 -    ( 86.7)  [-0.08 -0.94  0.33] vs [ 0.00 -0.95  0.30] (  5.0 )
  23 -    ( 84.2)  [-0.11 -0.85  0.52] vs [ 0.00 -0.89  0.46] (  7.8 )
  24 -    ( 87.5)  [-0.19 -0.75  0.63] vs [ 0.00 -0.82  0.57] ( 12.0 )
  25 -    ( 80.5)  [-0.68  0.03  0.73] vs [-0.69  0.00  0.72] (  1.8 )
  26 -    ( 82.8)  [-0.60 -0.01  0.80] vs [-0.63  0.00  0.78] (  2.6 )
  27 -    ( 86.2)  [-0.45 -0.03  0.89] vs [-0.49  0.00  0.87] (  3.0 )
  28 -    ( 81.2)  [-0.42 -0.05  0.91] vs [-0.37  0.00  0.93] (  4.1 )
  29 -    ( 91.3)  [-0.97 -0.02  0.24] vs [-0.98  0.00  0.22] (  1.8 )
  30 -    ( 85.6)  [-0.94 -0.02  0.33] vs [-0.95  0.00  0.30] (  2.1 )
  31 -    ( 78.5)  [-0.89 -0.00  0.45] vs [-0.89  0.00  0.46] (  0.2 )
  32 -    ( 90.2)  [-0.77 -0.01  0.64] vs [-0.82  0.00  0.57] (  4.9 )

Scripting

The following script from the Brainstorm distribution reproduces the analysis presented in this tutorial page: brainstorm3/toolbox/script/tutorial_phantom_elekta.m

function tutorial_phantom_elekta(tutorial_dir)
% TUTORIAL_PHANTOM_ELEKTA: Script that runs the tests for the Elekta phantom.
%
% CORRESPONDING ONLINE TUTORIAL:
%     https://neuroimage.usc.edu/brainstorm/Tutorials/PhantomElekta
%
% INPUTS: 
%     tutorial_dir: Directory where the sample_phantom.zip file has been unzipped

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
% 
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% 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.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Ken Taylor, John C Mosher, Francois Tadel, 2016

% ===== FILES TO IMPORT =====
% Does the sample folder exist
if ~exist(sprintf('%s',fullfile(tutorial_dir,'sample_phantom_elekta')),'dir'),
    error(sprintf('Missing ''sample_phantom_elekta'' folder in the tutorial directory ''%s''',tutorial_dir));
end
% Build the path of the files to import
FifFile{1} = fullfile(tutorial_dir, 'sample_phantom_elekta', 'kojak_all_2000nAm_pp_no_chpi_no_ms_raw.fif');
FifFile{2} = fullfile(tutorial_dir, 'sample_phantom_elekta', 'kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif');
FifFile{3} = fullfile(tutorial_dir, 'sample_phantom_elekta', 'kojak_all_20nAm_pp_no_chpi_no_ms_raw.fif');
% Select which file is processed with this script: 1=2000nAm, 2=200nAm, 3=20nAm
iFile = 2;
% Check if the folder contains the required files
if ~file_exist(FifFile{iFile})
    error(['The folder ' tutorial_dir ' does not contain the folder from the file sample_phantom_elekta.zip.']);
end


% ===== CREATE PROTOCOL =====
% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = 'TutorialPhantomElekta';
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
% Start a new report
bst_report('Start');


% ===== ANATOMY =====
% Subject name
SubjectName = 'Kojak';
% Generate the phantom anatomy
DipTrueFile = generate_phantom_elekta(SubjectName);

% ===== LINK CONTINUOUS FILES =====
% Process: Create link to raw files
sFilesKojak = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',  SubjectName, ...
    'datafile',     {FifFile{iFile}, 'FIF'}, ...
    'channelalign', 0);

% ===== READ EVENTS =====
% Process: Read from channel
bst_process('CallProcess', 'process_evt_read', sFilesKojak, [], ...
    'stimchan',  'STI201', ...
    'trackmode', 1, ...  % Value: detect the changes of channel value
    'zero',      0);
% Process: Delete spurious other events unrelated to dipoles
bst_process('CallProcess', 'process_evt_delete', sFilesKojak, [], ...
    'eventname', '256, 768, 1792, 3840, 4096, 6144, 7168, 7680, 7936');
% Process: Rename events to have a leading zero, for proper sorting
for i = 1:9
    bst_process('CallProcess', 'process_evt_rename', sFilesKojak, [], ...
        'src',  sprintf('%.0f',i), ...
        'dest', sprintf('%02.0f',i));
end
% Delete the first event of the first category (there is always an artifact)
LinkFile = file_fullpath(sFilesKojak(1).FileName);
LinkMat = load(LinkFile, 'F');
if ~isempty(LinkMat.F.events) && ~isempty(LinkMat.F.events(1).times)
    LinkMat.F.events(1).times(1)   = [];
    LinkMat.F.events(1).epochs(1)  = [];
    if ~isempty(LinkMat.F.events(1).channels)
        LinkMat.F.events(1).channels(1)= [];
    end
    if ~isempty(LinkMat.F.events(1).notes)
        LinkMat.F.events(1).notes(1) = [];
    end
end
bst_save(LinkFile, LinkMat, 'v6', 1);

% ===== IMPORT EVENTS =====
% Process: Import MEG/EEG: Events
%          DC offset correction: [-100ms,-10ms]
%          Resample: 100Hz
sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFilesKojak, [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'eventname',   ['01' sprintf(',%02.0f',2:32)], ...
    'timewindow',  [], ...
    'epochtime',   [-0.1, 0.3], ...
    'createcond',  0, ...
    'ignoreshort', 1, ...
    'usectfcomp',  0, ...
    'usessp',      0, ...
    'freq',        100, ...
    'baseline',    [-0.1, -0.001]);

% Process: Average: By trial group (folder average)
sAvgKojak = bst_process('CallProcess', 'process_average', sFilesEpochs, [], ...
    'avgtype',    5, ...  % By trial group (folder average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

% ===== SOURCE ESTIMATION =====
% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
    'baseline',       [-0.1, -0.01], ...
    'datatimewindow', [0, 0.3], ...
    'sensortypes',    'MEG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       0, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sAvgKojak, [], ...
    'Comment',     '', ...
    'sourcespace', 2, ...  % MRI volume
    'meg',         3, ...  % Overlapping spheres
    'volumegrid',  struct(...
         'Method',        'isotropic', ...
         'nLayers',       17, ...
         'Reduction',     3, ...
         'nVerticesInit', 4000, ...
         'Resolution',    0.0025, ...
         'FileName',      []));
% Process: Compute sources [2018]
sAvgSrcKojak = bst_process('CallProcess', 'process_inverse_2018', sAvgKojak, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment', 'Dipoles: MEG GRAD', ...
         'InverseMethod', 'gls', ...
         'InverseMeasure', 'performance', ...
         'SourceOrient', {{'free'}}, ...
         'Loose', 0.2, ...
         'UseDepth', 0, ...
         'WeightExp', 0.5, ...
         'WeightLimit', 10, ...
         'NoiseMethod', 'none', ...
         'NoiseReg', 0.1, ...
         'SnrMethod', 'rms', ...
         'SnrRms', 0, ...
         'SnrFixed', 3, ...
         'ComputeKernel', 1, ...
         'DataTypes', {{'MEG GRAD'}}));
     
 % ===== DIPOLE SCANNING =====
% Process: Dipole scanning
sDipScan = bst_process('CallProcess', 'process_dipole_scanning', sAvgSrcKojak, [], ...
    'timewindow', [0.06, 0.06], ...
    'scouts',     {});
% Merge all 32 dipoles together
DipMergeFile = dipoles_merge({sDipScan.FileName});

% Flip orientations
dip = load(DipMergeFile);
for i = 1:32
    dip.Dipole(i).Amplitude = dip.Dipole(i).Amplitude * sign(dip.Dipole(i).Amplitude(3));
end
dip.Comment = [dip.Comment ' | flipped'];
DipFlipFile = db_add(sDipScan(1).iStudy, dip);

% Merge with true locations
DipAllFile = dipoles_merge({DipTrueFile, DipFlipFile});
% visualize on the MRI 3D
view_dipoles(DipAllFile, 'Mri3D');


% ===== REPORT =====
% Get the True Dipole Locations
TrueDipoles = load(file_fullpath(DipTrueFile));
true_loc    = [TrueDipoles.Dipole.Loc];
true_orient = [TrueDipoles.Dipole.Amplitude];

% Display stats
fprintf('\n Dipole stats\n\n')
fprintf(' Dipole      Loc (mm)      Amp (nA-m)  Gof     Perf       Chi2     RChi2\n')
fprintf('--------------------------------------------------------------------------\n')
for i = 1:32
    fprintf('  %02.0f - [%5.1f %5.1f %5.1f]   %5.1f   %5.1f%%    %5.1f  %5.0f (%3.0f)   %.2f\n',...
        i, dip.Dipole(i).Loc*1000, norm(dip.Dipole(i).Amplitude)*1e9, ...
        dip.Dipole(i).Goodness*100, dip.Dipole(i).Perform, ...
        dip.Dipole(i).Khi2, dip.Dipole(i).DOF, dip.Dipole(i).Khi2/dip.Dipole(i).DOF),
end

% Compare errors, but dependent on order being correct
fprintf('\n Location errors from true\n\n')
fprintf(' Dipole      Loc (mm)            True (mm)         Diff [x y z]     Norm (mm)\n')
fprintf('------------------------------------------------------------------------------\n')
for i = 1:32
    temp_diff = dip.Dipole(i).Loc-true_loc(:,i);
    fprintf('  %02.0f - [%5.1f %5.1f %5.1f] [%5.1f %5.1f %5.1f] [%5.1f %5.1f %5.1f]   (%5.1f)\n',...
        i, dip.Dipole(i).Loc*1000, true_loc(:,i)*1000, temp_diff*1000, norm(temp_diff)*1000);
end

% Orientation errors
fprintf('\n Orientation errors from true\n\n')
fprintf(' Dipole  Amp (nA-m)      [X Y Z]              TRUE [X Y Z]     Degrees\n')
fprintf('------------------------------------------------------------------------\n')
for i = 1:32
    Amp    = dip.Dipole(i).Amplitude; 
    nAmp   = norm(Amp); 
    Orient = Amp / nAmp;
    fprintf('  %02.0f -    (%5.1f)  [%5.2f %5.2f %5.2f] vs [%5.2f %5.2f %5.2f] (%5.1f )\n',...
        i, nAmp*1e9, Orient, true_orient(:,i), (subspace(Orient,true_orient(:,i)))*180/pi);
end

% Save and display report
ReportFile = bst_report('Save', []);
bst_report('Open', ReportFile);












Feedback: Comments, bug reports, suggestions, questions

Email address (if you expect an answer):