MEG current phantom (CTF)

Authors: Francois Tadel, Elizabeth Bock

This tutorial explains how to import and process CTF 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 Elekta-Neuromag phantom.

License

This tutorial dataset remains a property of its authors: Elizabeth Bock, Francois Tadel and Sylvain Baillet (MEG Lab, McConnell Brain Imaging Center, Montreal Neurological Institute, McGill University, Canada).

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

CTF current phantom

The CTF current dipole phantom is a spherical container filled with a conducting saline solution which contains a current source and sink. This electric dipole is used to simulate brain sources in a conductive medium. The globe has a 130mm inner diameter and small attachment posts to attach the three CTF head localization coils. The position of this dipole can be adjusted within this globe.

The dipole itself is constructed of two gold spheres about 2 mm in diameter, separated by 9.0 mm center to center. The dipole moment can be calculated by the equation m=I.L, where I is the current flow (in Amperes) and L is the length of the dipole (0.009 meters).

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).

Reference

VSM/CTF documentation: PN900-0018, Revision 3.2, 23 November 2006. This document can be found with your full CTF installation at /opt/ctf/docs/Phantom.pdf.

Description of the experiment

Files distributed as part of the CTF phantom tutorial:

Download and installation

Generate anatomy

Access the recordings

Import recordings

Event detection

The sinusoidal signal is generated by the CTF hardware on channel HDAC006. While there are some automatic trigger events generated by the system that can be used for importing, we will have a more precise event average if the events are detected again offline.

Import and average

Noise covariance

Use the empty room recordings.

Source estimation

Compute a forward model and inverse model for a regular grid inside the phantom volume.

Dipole scanning

Scan the for the most significant dipole in the grid of computed dipoles estimate previously.

Dipole fitting with FieldTrip

Perform a non-linear dipole fit with the function ft_dipolefitting from the FieldTrip toolbox.

Results comparison

Condition

Method

Forward model

X

Y

Z

GOF

Nominal location

0

-18

49

mm

200uA

Scanning (5mm)

Single sphere

-1.00

-16.00

44.00

99.80%

Scanning (5mm)

Overlapping spheres

-1.00

-16.00

44.00

99.81%

Scanning (5mm)

OpenMEEG BEM

-1.00

-16.00

44.00

99.76%

Fitting

Single sphere

-1.04

-17.00

43.98

99.95%

Fitting

Overlapping spheres

-1.05

-16.98

44.00

99.95%

CTF software

-1

-17

44

99.95%

MNE software

-0.79

-17.00

43.98

99.9%

20uA

Scanning (5mm)

Single sphere

-1.00

-16.00

44.00

96.94%

Scanning (5mm)

Overlapping spheres

-1.00

-16.00

44.00

96.96%

Scanning (5mm)

OpenMEEG BEM

-1.00

-16.00

44.00

96.90%

Fitting

Single sphere

-1.75

-17.13

44.39

98.25%

Fitting

Overlapping spheres

-1.78

-17.14

44.45

98.25%

CTF software

-1

-17

44

98.38%

MNE software

-1.38

-16.31

44.01

99.1%

The nominal location indicates where the dipole is supposed to be, relative to the center of the sphere. It is measured with rulers with an overall precision of about 5mm. The range of discrepancy we observe between this nominal location and the position of the dipole estimated from the recordings is acceptable.

Advanced

Digitized head points

The head points collected with the Brainstorm digitizer are usually copied to the .ds folders and imported automatically when loading the recordings. We decided not to include them in this example because in the case of this current phantom, there is no ambiguity in the definition of the anatomical fiducials. As this refined registration with the .pos files is not part of the standard CTF workflow, not including it will make it easier to compare the workflow and results with other programs.

For additional testing purposes, the .pos file for the phantom is included in the sample_phantom_ctf.zip package, but you have to add it manually to the recordings. Do not use these points to refine automatically the registration: the fitting algorithm may fail finding the best rotation around the Z axis because the phantom is completely spherical, and the registration is already close to perfection.

Scripting

A script available in the Brainstorm distribution reproduces the different steps of this tutorial:

brainstorm3/toolbox/script/tutorial_phantom_ctf.m

1 function tutorial_phantom_ctf(tutorial_dir) 2 % TUTORIAL_PHANTOM_CTF: Script that runs the tests for the CTF current phantom. 3 % 4 % CORRESPONDING ONLINE TUTORIAL: 5 % https://neuroimage.usc.edu/brainstorm/Tutorials/PhantomCtf 6 % 7 % INPUTS: 8 % tutorial_dir: Directory where the sample_phantom.zip file has been unzipped 9 10 % @============================================================================= 11 % This function is part of the Brainstorm software: 12 % https://neuroimage.usc.edu/brainstorm 13 % 14 % Copyright (c) University of Southern California & McGill University 15 % This software is distributed under the terms of the GNU General Public License 16 % as published by the Free Software Foundation. Further details on the GPLv3 17 % license can be found at http://www.gnu.org/copyleft/gpl.html. 18 % 19 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 20 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 21 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 22 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 23 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 24 % 25 % For more information type "brainstorm license" at command prompt. 26 % =============================================================================@ 27 % 28 % Author: Francois Tadel, 2016 29 30 31 % ===== FILES TO IMPORT ===== 32 % You have to specify the folder in which the tutorial dataset is unzipped 33 if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir) 34 error('The first argument must be the full path to the dataset folder.'); 35 end 36 % Build the path of the files to import 37 Run200FileDs = fullfile(tutorial_dir, 'sample_phantom_ctf', 'ds', 'phantom_200uA_20150709_01.ds'); 38 Run20FileDs = fullfile(tutorial_dir, 'sample_phantom_ctf', 'ds', 'phantom_20uA_20150603_03.ds'); 39 NoiseFileDs = fullfile(tutorial_dir, 'sample_phantom_ctf', 'ds', 'emptyroom_20150709_01.ds'); 40 Run200FileFif = fullfile(tutorial_dir, 'sample_phantom_ctf', 'fif', 'phantom_200uA_20150709_01.fif'); 41 Run20FileFif = fullfile(tutorial_dir, 'sample_phantom_ctf', 'fif', 'phantom_20uA_20150603_03.fif'); 42 PosFile = fullfile(tutorial_dir, 'sample_phantom_ctf', 'ds', 'phantom_20160222_01.pos'); 43 % Check if the folder contains the required files 44 if ~file_exist(Run200FileDs) 45 error(['The folder ' tutorial_dir ' does not contain the folder from the file sample_phantom_ctf.zip.']); 46 end 47 48 % ===== CREATE PROTOCOL ===== 49 % The protocol name has to be a valid folder name (no spaces, no weird characters...) 50 ProtocolName = 'TutorialPhantom'; 51 % Start brainstorm without the GUI 52 if ~brainstorm('status') 53 brainstorm nogui 54 end 55 % Delete existing protocol 56 gui_brainstorm('DeleteProtocol', ProtocolName); 57 % Create new protocol 58 gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); 59 % Start a new report 60 bst_report('Start'); 61 62 63 % ===== ANATOMY ===== 64 % Subject name 65 SubjectNameDs = 'PhantomCTF-ds'; 66 SubjectNameFif = 'PhantomCTF-fif'; 67 % Generate the phantom anatomy 68 generate_phantom_ctf(SubjectNameDs); 69 generate_phantom_ctf(SubjectNameFif); 70 71 % ===== LINK CONTINUOUS FILES ===== 72 % Process: Create link to raw files 73 sFilesRun200Ds = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 74 'subjectname', SubjectNameDs, ... 75 'datafile', {Run200FileDs, 'CTF'}, ... 76 'channelalign', 0); 77 sFilesRun20Ds = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 78 'subjectname', SubjectNameDs, ... 79 'datafile', {Run20FileDs, 'CTF'}, ... 80 'channelalign', 0); 81 sFilesNoiseDs = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 82 'subjectname', SubjectNameDs, ... 83 'datafile', {NoiseFileDs, 'CTF'}, ... 84 'channelalign', 0); 85 sFilesRun200Fif = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 86 'subjectname', SubjectNameFif, ... 87 'datafile', {Run200FileFif, 'FIF'}, ... 88 'channelalign', 0); 89 sFilesRun20Fif = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 90 'subjectname', SubjectNameFif, ... 91 'datafile', {Run20FileFif, 'FIF'}, ... 92 'channelalign', 0); 93 sFilesRawDs = [sFilesRun200Ds, sFilesRun20Ds]; 94 sFilesRawFif = [sFilesRun200Fif, sFilesRun20Fif]; 95 sFilesRawAll = [sFilesRawDs, sFilesRawFif]; 96 97 % Process: Convert to continuous (CTF): Continuous 98 bst_process('CallProcess', 'process_ctf_convert', [sFilesRawDs, sFilesNoiseDs], [], ... 99 'rectype', 2); % Continuous 100 101 102 % ===== ADD HEAD POINTS ===== 103 % Process: Add head points 104 bst_process('CallProcess', 'process_headpoints_add', sFilesRawAll, [], ... 105 'channelfile', {PosFile, 'ASCII_NXYZ'}); 106 107 % Process: Snapshot: Sensors/MRI registration 108 bst_process('CallProcess', 'process_snapshot', [sFilesRun200Ds, sFilesRun20Ds], [], ... 109 'target', 1, ... % Sensors/MRI registration 110 'modality', 1, ... % MEG (All) 111 'orient', 1, ... % left 112 'Comment', 'MEG/MRI Registration'); 113 114 115 % ===== DETECT EVENTS ===== 116 % Process: Detect: stim 117 bst_process('CallProcess', 'process_evt_detect_threshold', sFilesRawAll, [], ... 118 'eventname', 'stim', ... 119 'channelname', 'HDAC006', ... 120 'timewindow', [0, 99.99833333], ... 121 'thresholdMAX', 0.5, ... 122 'units', 1, ... % None (10^0) 123 'bandpass', [], ... 124 'isAbsolute', 0, ... 125 'isDCremove', 1); 126 % Process: Convert to simple event 127 bst_process('CallProcess', 'process_evt_simple', sFilesRawAll, [], ... 128 'eventname', 'stim', ... 129 'method', 2); % Keep the middle of the events 130 131 132 % ===== IMPORT EVENTS ===== 133 % Process: Import MEG/EEG: Events 134 sFilesEpochsAll = bst_process('CallProcess', 'process_import_data_event', sFilesRawAll, [], ... 135 'subjectname', [], ... 136 'condition', [], ... 137 'eventname', 'stim', ... 138 'timewindow', [0, 10], ... 139 'epochtime', [-0.07, 0.07], ... 140 'createcond', 0, ... 141 'ignoreshort', 1, ... 142 'usectfcomp', 1, ... 143 'usessp', 1, ... 144 'freq', [], ... 145 'baseline', [-0.07, 0.07]); % Remove DC: Entire file 146 % Process: Average: By folder (subject average) 147 sAvgAll = bst_process('CallProcess', 'process_average', sFilesEpochsAll, [], ... 148 'avgtype', 3, ... % By folder (subject average) 149 'avg_func', 1, ... % Arithmetic average: mean(x) 150 'weighted', 0, ... 151 'keepevents', 0); 152 153 % Process: Snapshot: Recordings time series 154 bst_process('CallProcess', 'process_snapshot', sAvgAll, [], ... 155 'target', 5, ... % Recordings time series 156 'modality', 1); % MEG (All) 157 % Process: Snapshot: Recordings topography (one time) 158 bst_process('CallProcess', 'process_snapshot', sAvgAll, [], ... 159 'target', 6, ... % Recordings topography (one time) 160 'modality', 1, ... % MEG (All) 161 'time', 0); 162 163 164 % ===== NOISE COVARIANCE ===== 165 % Process: Compute noise covariance 166 bst_process('CallProcess', 'process_noisecov', sFilesNoiseDs, [], ... 167 'baseline', [], ... 168 'sensortypes', 'MEG, EEG, SEEG, ECOG', ... 169 'target', 1, ... % Noise covariance 170 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data 171 'identity', 0, ... 172 'copycond', 1, ... 173 'copysubj', 1, ... 174 'replacefile', 1); % Replace 175 % Process: Snapshot: Noise covariance 176 bst_process('CallProcess', 'process_snapshot', sFilesNoiseDs, [], ... 177 'target', 3, ... % Noise covariance 178 'Comment', 'Noise covariance'); 179 180 181 % ===== SOURCE ESTIMATION ===== 182 % Define source grid options 183 VolumeGrid = struct(... 184 'Method', 'isotropic', ... 185 'nLayers', 17, ... 186 'Reduction', 3, ... 187 'nVerticesInit', 4000, ... 188 'Resolution', 0.005); 189 % Process: Compute head model 190 bst_process('CallProcess', 'process_headmodel', sAvgAll, [], ... 191 'Comment', '', ... 192 'sourcespace', 2, ... % MRI volume 193 'volumegrid', VolumeGrid, ... 194 'meg', 2); % Single sphere 195 % Define inverse options 196 InverseOptions = struct(... 197 'Comment', 'Dipoles: Single sphere', ... 198 'InverseMethod', 'gls', ... 199 'InverseMeasure', 'performance', ... 200 'SourceOrient', {{'free'}}, ... 201 'Loose', 0.2, ... 202 'UseDepth', 1, ... 203 'WeightExp', 0.5, ... 204 'WeightLimit', 10, ... 205 'NoiseMethod', 'reg', ... 206 'NoiseReg', 0.1, ... 207 'SnrMethod', 'fixed', ... 208 'SnrRms', 0.001, ... 209 'SnrFixed', 3, ... 210 'ComputeKernel', 1, ... 211 'DataTypes', {{'MEG'}}); 212 % Process: Compute sources [2018] 213 sAvgSrcAll = bst_process('CallProcess', 'process_inverse_2018', sAvgAll, [], ... 214 'output', 2, ... % Kernel only: one per file 215 'inverse', InverseOptions); 216 % Process: Snapshot: Sources (one time) 217 bst_process('CallProcess', 'process_snapshot', sAvgSrcAll, [], ... 218 'target', 8, ... % Sources (one time) 219 'orient', 6, ... % back 220 'time', 0, ... 221 'threshold', 40); 222 223 224 % ===== DIPOLE SCANNING ===== 225 % Process: Dipole scanning 226 sDipScan = bst_process('CallProcess', 'process_dipole_scanning', sAvgSrcAll, [], ... 227 'timewindow', [0, 0], ... 228 'scouts', {}); 229 % Process: Snapshot: Dipoles 230 bst_process('CallProcess', 'process_snapshot', sDipScan, [], ... 231 'target', 13, ... % Dipoles 232 'orient', 3); % top 233 234 235 % ===== FIELDTRIP DIPOLE FITTING ===== 236 % Process: FieldTrip: ft_dipolefitting 237 sDipFit = bst_process('CallProcess', 'process_ft_dipolefitting', sAvgAll, [], ... 238 'filetag', 'Single sphere', ... 239 'timewindow', [0, 0], ... 240 'sensortypes', 'MEG', ... 241 'dipolemodel', 1, ... % Moving dipole 242 'numdipoles', 1, ... 243 'symmetry', 0); 244 % Process: Snapshot: Dipoles 245 bst_process('CallProcess', 'process_snapshot', sDipFit, [], ... 246 'target', 13, ... % Dipoles 247 'orient', 3); % top 248 249 250 % ===== REPORT RESULTS: SINGLE SPHERE ===== 251 strResults = ['CTF current phantom' 10 '==============================================' 10]; 252 % Dipole scanning 253 strResults = [strResults, 'Single sphere / Scanning:' 10]; 254 for i = 1:length(sDipScan) 255 dipMat = load(file_fullpath(sDipScan(i).FileName)); 256 strResults = [strResults, sprintf(' - %s: [%1.2f, %1.2f, %1.2f]mm %1.2f%% gof\n', bst_fileparts(sDipScan(i).FileName), dipMat.Dipole(1).Loc.*1000, dipMat.Dipole(1).Goodness.*100)]; 257 end 258 % Dipole fitting 259 strResults = [strResults, 'Single sphere / Fitting:' 10]; 260 for i = 1:length(sDipFit) 261 dipMat = load(file_fullpath(sDipFit(i).FileName)); 262 strResults = [strResults, sprintf(' - %s: [%1.2f, %1.2f, %1.2f]mm %1.2f%% gof\n', bst_fileparts(sDipFit(i).FileName), dipMat.Dipole(1).Loc.*1000, dipMat.Dipole(1).Goodness.*100)]; 263 end 264 % Report results 265 bst_report('Info', 'process_dipole_scanning', [], strResults); 266 disp([10 '==============================================' 10 strResults]); 267 268 269 % ===== OVERLAPPING SPHERES ===== 270 % Process: Compute head model 271 bst_process('CallProcess', 'process_headmodel', sAvgAll, [], ... 272 'Comment', '', ... 273 'sourcespace', 2, ... % MRI volume 274 'volumegrid', VolumeGrid, ... 275 'meg', 3); % Overlapping spheres 276 % Process: Compute sources [2018] 277 InverseOptions.Comment = 'Dipoles: Overlapping spheres'; 278 sSrcOs = bst_process('CallProcess', 'process_inverse_2018', sAvgAll, [], ... 279 'output', 2, ... % Kernel only: one per file 280 'inverse', InverseOptions); 281 % Process: Snapshot: Sources (one time) 282 bst_process('CallProcess', 'process_snapshot', sSrcOs, [], ... 283 'target', 8, ... % Sources (one time) 284 'orient', 6, ... % back 285 'time', 0, ... 286 'threshold', 40); 287 288 % Process: Dipole scanning 289 sDipScanOs = bst_process('CallProcess', 'process_dipole_scanning', sSrcOs, [], ... 290 'timewindow', [0, 0], ... 291 'scouts', {}); 292 % Process: FieldTrip: ft_dipolefitting 293 sDipFitOs = bst_process('CallProcess', 'process_ft_dipolefitting', sAvgAll, [], ... 294 'filetag', 'Overlapping spheres', ... 295 'timewindow', [0, 0], ... 296 'sensortypes', 'MEG', ... 297 'dipolemodel', 1, ... % Moving dipole 298 'numdipoles', 1, ... 299 'symmetry', 0); 300 301 302 % ===== OVERLAPPING SPHERES: REPORT RESULTS ===== 303 strResults = ''; 304 % Dipole scanning 305 strResults = [strResults, 'Overlapping spheres / Scanning:' 10]; 306 for i = 1:length(sDipScanOs) 307 dipMat = load(file_fullpath(sDipScanOs(i).FileName)); 308 strResults = [strResults, sprintf(' - %s: [%1.2f, %1.2f, %1.2f]mm %1.2f%% gof\n', bst_fileparts(sDipScanOs(i).FileName), dipMat.Dipole(1).Loc.*1000, dipMat.Dipole(1).Goodness.*100)]; 309 end 310 % Dipole fitting 311 strResults = [strResults, 'Overlapping spheres / Fitting:' 10]; 312 for i = 1:length(sDipFitOs) 313 dipMat = load(file_fullpath(sDipFitOs(i).FileName)); 314 strResults = [strResults, sprintf(' - %s: [%1.2f, %1.2f, %1.2f]mm %1.2f%% gof\n', bst_fileparts(sDipFitOs(i).FileName), dipMat.Dipole(1).Loc.*1000, dipMat.Dipole(1).Goodness.*100)]; 315 end 316 % Report results 317 bst_report('Info', 'process_dipole_scanning', [], strResults); 318 disp(strResults); 319 320 321 % ===== OPENMEEG BEM ===== 322 % Process: Generate BEM surfaces (DS) 323 bst_process('CallProcess', 'process_generate_bem', [], [], ... 324 'subjectname', SubjectNameDs, ... 325 'nscalp', 362, ... 326 'nouter', 162, ... 327 'ninner', 162, ... 328 'thickness', 4); 329 % Process: Generate BEM surfaces (FIF) 330 bst_process('CallProcess', 'process_generate_bem', [], [], ... 331 'subjectname', SubjectNameFif, ... 332 'nscalp', 362, ... 333 'nouter', 162, ... 334 'ninner', 162, ... 335 'thickness', 4); 336 337 % Process: Compute head model 338 bst_process('CallProcess', 'process_headmodel', sAvgAll, [], ... 339 'Comment', '', ... 340 'sourcespace', 2, ... % MRI volume 341 'volumegrid', VolumeGrid, ... 342 'meg', 4); % OpenMEEG BEM 343 % Process: Compute sources [2018] 344 InverseOptions.Comment = 'Dipoles: OpenMEEG BEM'; 345 sSrcBem = bst_process('CallProcess', 'process_inverse_2018', sAvgAll, [], ... 346 'output', 2, ... % Kernel only: one per file 347 'inverse', InverseOptions); 348 % Process: Snapshot: Sources (one time) 349 bst_process('CallProcess', 'process_snapshot', sSrcBem, [], ... 350 'target', 8, ... % Sources (one time) 351 'orient', 6, ... % back 352 'time', 0, ... 353 'threshold', 40); 354 % Process: Dipole scanning 355 sDipScanBem = bst_process('CallProcess', 'process_dipole_scanning', sSrcBem, [], ... 356 'timewindow', [0, 0], ... 357 'scouts', {}); 358 359 360 % ===== OPENMEEG BEM: REPORT RESULTS ===== 361 strResults = ''; 362 % Dipole scanning 363 strResults = [strResults, 'OpenMEEG BEM / Scanning:' 10]; 364 for i = 1:length(sDipScanBem) 365 dipMat = load(file_fullpath(sDipScanBem(i).FileName)); 366 strResults = [strResults, sprintf(' - %s: [%1.2f, %1.2f, %1.2f]mm %1.2f%% gof\n', bst_fileparts(sDipScanBem(i).FileName), dipMat.Dipole(1).Loc.*1000, dipMat.Dipole(1).Goodness.*100)]; 367 end 368 % Report results 369 bst_report('Info', 'process_dipole_scanning', [], strResults); 370 disp(strResults); 371 372 373 374 % Save and display report 375 ReportFile = bst_report('Save', []); 376 bst_report('Open', ReportFile); 377 378





Feedback: Comments, bug reports, suggestions, questions
Email address (if you expect an answer):


Tutorials/PhantomCtf (last edited 2021-04-01 09:50:25 by FrancoisTadel)