4D Neuroimaging/BTi phantom dataset tutorial

Here we read 4DBTi epochs data obtained with a spherical phantom using four different dipole locations. For each condition we compute evoked data and compute dipole fits.

Data are provided by Jean-Michel Badier from MEG center in Marseille, France.

# Authors: Alex Gramfort <alexandre.gramfort@inria.fr>
#
# License: BSD (3-clause)

import os.path as op
import numpy as np
from mayavi import mlab
from mne.datasets import phantom_4dbti
import mne

Read data and compute a dipole fit at the peak of the evoked response

data_path = phantom_4dbti.data_path()
raw_fname = op.join(data_path, '%d/e,rfhp1.0Hz')

dipoles = list()
sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=0.080)

t0 = 0.07  # peak of the response

pos = np.empty((4, 3))

for ii in range(4):
    raw = mne.io.read_raw_bti(raw_fname % (ii + 1,),
                              rename_channels=False, preload=True)
    raw.info['bads'] = ['A173', 'A213', 'A232']
    events = mne.find_events(raw, 'TRIGGER', mask=4350, mask_type='not_and')
    epochs = mne.Epochs(raw, events=events, event_id=8192, tmin=-0.2, tmax=0.4,
                        preload=True)
    evoked = epochs.average()
    evoked.plot(time_unit='s')
    cov = mne.compute_covariance(epochs, tmax=0.)
    dip = mne.fit_dipole(evoked.copy().crop(t0, t0), cov, sphere)[0]
    pos[ii] = dip.pos[0]
  • ../../_images/sphx_glr_plot_phantom_4DBTi_001.png
  • ../../_images/sphx_glr_plot_phantom_4DBTi_002.png
  • ../../_images/sphx_glr_plot_phantom_4DBTi_003.png
  • ../../_images/sphx_glr_plot_phantom_4DBTi_004.png

Out:

Equiv. model fitting -> RV = 0.0037285 %
mu1 = 0.943943    lambda1 = 0.139085
mu2 = 0.665517    lambda2 = 0.684836
mu3 = -0.0972643    lambda3 = -0.0135513
Set up EEG sphere model with scalp radius    80.0 mm

Reading 4D PDF file /home/circleci/mne_data/MNE-phantom-4DBTi/1/e,rfhp1.0Hz...
Creating Neuromag info structure ...
... Setting channel info structure.
... putting coil transforms in Neuromag coordinates
... Reading digitization points from /home/circleci/mne_data/MNE-phantom-4DBTi/1/hs_file
... putting digitization points in Neuromag coordinates
... Computing new device to head transform.
Done.
Currently direct inclusion of 4D weight tables is not supported. For critical use cases please take into account the MNE command "mne_create_comp_data" to include weights as printed out by the 4D "print_table" routine.
Current compensation grade : 0
Reading 0 ... 13599  =      0.000 ...    20.052 secs...
Trigger channel has a non-zero initial value of 4350 (consider using initial_event=True to detect this event)
50 events found
Event IDs: [8192]
50 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 50 events and 408 original time points ...
2 bad epochs dropped
Computing data rank from raw with rank=None
    Using tolerance 3.7e-09 (2.2e-16 eps * 245 dim * 6.9e+04  max singular value)
    Estimated rank (mag): 245
    MAG: rank 245 computed from 245 data channels with 0 projectors
Reducing data rank from 245 -> 245
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 6576
[done]
BEM               : <ConductorModel  |  Sphere (3 layers): r0=[0.0, 0.0, 0.0] R=80 mm>
MRI transform     : identity
Sphere model      : origin at (   0.00    0.00    0.00) mm, rad =    0.1 mm
Guess grid        :   20.0 mm
Guess mindist     :    5.0 mm
Guess exclude     :   20.0 mm
Using standard MEG coil definitions.

Coordinate transformation: MRI (surface RAS) -> head
     1.000000  0.000000  0.000000       0.00 mm
     0.000000  1.000000  0.000000       0.00 mm
     0.000000  0.000000  1.000000       0.00 mm
     0.000000  0.000000  0.000000       1.00
Coordinate transformation: MEG device -> head
     0.975564 -0.033891 -0.217085      -1.50 mm
     0.044586  0.998011  0.044560     -34.72 mm
     0.215143 -0.053150  0.975135      -6.31 mm
     0.000000  0.000000  0.000000       1.00
3 bad channels total
Read 245 MEG channels from info
Read  23 MEG compensation channels from info
84 coil definitions read
Coordinate transformation: MEG device -> head
     0.975564 -0.033891 -0.217085      -1.50 mm
     0.044586  0.998011  0.044560     -34.72 mm
     0.215143 -0.053150  0.975135      -6.31 mm
     0.000000  0.000000  0.000000       1.00
0 compensation data sets in info
MEG coil definitions created in head coordinates.
Decomposing the sensor noise covariance matrix...
Computing data rank from covariance with rank=None
    Using tolerance 3.9e-14 (2.2e-16 eps * 245 dim * 0.72  max singular value)
    Estimated rank (mag): 245
    MAG: rank 245 computed from 245 data channels with 0 projectors
    Setting small MAG eigenvalues to zero (without PCA)
    Created the whitener using a noise covariance matrix with rank 245 (0 small eigenvalues omitted)

---- Computing the forward solution for the guesses...
Making a spherical guess space with radius    72.0 mm...
Filtering (grid =     20 mm)...
Surface CM = (   0.0    0.0    0.0) mm
Surface fits inside a sphere with radius   72.0 mm
Surface extent:
    x =  -72.0 ...   72.0 mm
    y =  -72.0 ...   72.0 mm
    z =  -72.0 ...   72.0 mm
Grid extent:
    x =  -80.0 ...   80.0 mm
    y =  -80.0 ...   80.0 mm
    z =  -80.0 ...   80.0 mm
729 sources before omitting any.
178 sources after omitting infeasible sources.
170 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
Go through all guess source locations...
[done 170 sources]
---- Fitted :    69.3 ms
No projector specified for this dataset. Please consider the method self.add_proj.
1 time points fitted
Reading 4D PDF file /home/circleci/mne_data/MNE-phantom-4DBTi/2/e,rfhp1.0Hz...
Creating Neuromag info structure ...
... Setting channel info structure.
... putting coil transforms in Neuromag coordinates
... Reading digitization points from /home/circleci/mne_data/MNE-phantom-4DBTi/2/hs_file
... putting digitization points in Neuromag coordinates
... Computing new device to head transform.
Done.
Currently direct inclusion of 4D weight tables is not supported. For critical use cases please take into account the MNE command "mne_create_comp_data" to include weights as printed out by the 4D "print_table" routine.
Current compensation grade : 0
Reading 0 ... 13599  =      0.000 ...    20.052 secs...
Trigger channel has a non-zero initial value of 4350 (consider using initial_event=True to detect this event)
50 events found
Event IDs: [8192]
50 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 50 events and 408 original time points ...
2 bad epochs dropped
Computing data rank from raw with rank=None
    Using tolerance 3.5e-09 (2.2e-16 eps * 245 dim * 6.4e+04  max singular value)
    Estimated rank (mag): 245
    MAG: rank 245 computed from 245 data channels with 0 projectors
Reducing data rank from 245 -> 245
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 6576
[done]
BEM               : <ConductorModel  |  Sphere (3 layers): r0=[0.0, 0.0, 0.0] R=80 mm>
MRI transform     : identity
Sphere model      : origin at (   0.00    0.00    0.00) mm, rad =    0.1 mm
Guess grid        :   20.0 mm
Guess mindist     :    5.0 mm
Guess exclude     :   20.0 mm
Using standard MEG coil definitions.

Coordinate transformation: MRI (surface RAS) -> head
     1.000000  0.000000  0.000000       0.00 mm
     0.000000  1.000000  0.000000       0.00 mm
     0.000000  0.000000  1.000000       0.00 mm
     0.000000  0.000000  0.000000       1.00
Coordinate transformation: MEG device -> head
     0.975554 -0.034041 -0.217109      -1.51 mm
     0.044503  0.998063  0.043482     -34.64 mm
     0.215208 -0.052081  0.975178      -6.31 mm
     0.000000  0.000000  0.000000       1.00
3 bad channels total
Read 245 MEG channels from info
Read  23 MEG compensation channels from info
84 coil definitions read
Coordinate transformation: MEG device -> head
     0.975554 -0.034041 -0.217109      -1.51 mm
     0.044503  0.998063  0.043482     -34.64 mm
     0.215208 -0.052081  0.975178      -6.31 mm
     0.000000  0.000000  0.000000       1.00
0 compensation data sets in info
MEG coil definitions created in head coordinates.
Decomposing the sensor noise covariance matrix...
Computing data rank from covariance with rank=None
    Using tolerance 3.4e-14 (2.2e-16 eps * 245 dim * 0.63  max singular value)
    Estimated rank (mag): 245
    MAG: rank 245 computed from 245 data channels with 0 projectors
    Setting small MAG eigenvalues to zero (without PCA)
    Created the whitener using a noise covariance matrix with rank 245 (0 small eigenvalues omitted)

---- Computing the forward solution for the guesses...
Making a spherical guess space with radius    72.0 mm...
Filtering (grid =     20 mm)...
Surface CM = (   0.0    0.0    0.0) mm
Surface fits inside a sphere with radius   72.0 mm
Surface extent:
    x =  -72.0 ...   72.0 mm
    y =  -72.0 ...   72.0 mm
    z =  -72.0 ...   72.0 mm
Grid extent:
    x =  -80.0 ...   80.0 mm
    y =  -80.0 ...   80.0 mm
    z =  -80.0 ...   80.0 mm
729 sources before omitting any.
178 sources after omitting infeasible sources.
170 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
Go through all guess source locations...
[done 170 sources]
---- Fitted :    69.3 ms
No projector specified for this dataset. Please consider the method self.add_proj.
1 time points fitted
Reading 4D PDF file /home/circleci/mne_data/MNE-phantom-4DBTi/3/e,rfhp1.0Hz...
Creating Neuromag info structure ...
... Setting channel info structure.
... putting coil transforms in Neuromag coordinates
... Reading digitization points from /home/circleci/mne_data/MNE-phantom-4DBTi/3/hs_file
... putting digitization points in Neuromag coordinates
... Computing new device to head transform.
Done.
Currently direct inclusion of 4D weight tables is not supported. For critical use cases please take into account the MNE command "mne_create_comp_data" to include weights as printed out by the 4D "print_table" routine.
Current compensation grade : 0
Reading 0 ... 13599  =      0.000 ...    20.052 secs...
Trigger channel has a non-zero initial value of 4350 (consider using initial_event=True to detect this event)
50 events found
Event IDs: [8192]
50 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 50 events and 408 original time points ...
2 bad epochs dropped
Computing data rank from raw with rank=None
    Using tolerance 2.6e-09 (2.2e-16 eps * 245 dim * 4.7e+04  max singular value)
    Estimated rank (mag): 245
    MAG: rank 245 computed from 245 data channels with 0 projectors
Reducing data rank from 245 -> 245
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 6576
[done]
BEM               : <ConductorModel  |  Sphere (3 layers): r0=[0.0, 0.0, 0.0] R=80 mm>
MRI transform     : identity
Sphere model      : origin at (   0.00    0.00    0.00) mm, rad =    0.1 mm
Guess grid        :   20.0 mm
Guess mindist     :    5.0 mm
Guess exclude     :   20.0 mm
Using standard MEG coil definitions.

Coordinate transformation: MRI (surface RAS) -> head
     1.000000  0.000000  0.000000       0.00 mm
     0.000000  1.000000  0.000000       0.00 mm
     0.000000  0.000000  1.000000       0.00 mm
     0.000000  0.000000  0.000000       1.00
Coordinate transformation: MEG device -> head
     0.975577 -0.033678 -0.217061      -1.49 mm
     0.044611  0.997960  0.045666     -34.78 mm
     0.215080 -0.054233  0.975089      -6.31 mm
     0.000000  0.000000  0.000000       1.00
3 bad channels total
Read 245 MEG channels from info
Read  23 MEG compensation channels from info
84 coil definitions read
Coordinate transformation: MEG device -> head
     0.975577 -0.033678 -0.217061      -1.49 mm
     0.044611  0.997960  0.045666     -34.78 mm
     0.215080 -0.054233  0.975089      -6.31 mm
     0.000000  0.000000  0.000000       1.00
0 compensation data sets in info
MEG coil definitions created in head coordinates.
Decomposing the sensor noise covariance matrix...
Computing data rank from covariance with rank=None
    Using tolerance 1.9e-14 (2.2e-16 eps * 245 dim * 0.34  max singular value)
    Estimated rank (mag): 245
    MAG: rank 245 computed from 245 data channels with 0 projectors
    Setting small MAG eigenvalues to zero (without PCA)
    Created the whitener using a noise covariance matrix with rank 245 (0 small eigenvalues omitted)

---- Computing the forward solution for the guesses...
Making a spherical guess space with radius    72.0 mm...
Filtering (grid =     20 mm)...
Surface CM = (   0.0    0.0    0.0) mm
Surface fits inside a sphere with radius   72.0 mm
Surface extent:
    x =  -72.0 ...   72.0 mm
    y =  -72.0 ...   72.0 mm
    z =  -72.0 ...   72.0 mm
Grid extent:
    x =  -80.0 ...   80.0 mm
    y =  -80.0 ...   80.0 mm
    z =  -80.0 ...   80.0 mm
729 sources before omitting any.
178 sources after omitting infeasible sources.
170 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
Go through all guess source locations...
[done 170 sources]
---- Fitted :    69.3 ms
No projector specified for this dataset. Please consider the method self.add_proj.
1 time points fitted
Reading 4D PDF file /home/circleci/mne_data/MNE-phantom-4DBTi/4/e,rfhp1.0Hz...
Creating Neuromag info structure ...
... Setting channel info structure.
... putting coil transforms in Neuromag coordinates
... Reading digitization points from /home/circleci/mne_data/MNE-phantom-4DBTi/4/hs_file
... putting digitization points in Neuromag coordinates
... Computing new device to head transform.
Done.
Currently direct inclusion of 4D weight tables is not supported. For critical use cases please take into account the MNE command "mne_create_comp_data" to include weights as printed out by the 4D "print_table" routine.
Current compensation grade : 0
Reading 0 ... 13599  =      0.000 ...    20.052 secs...
Trigger channel has a non-zero initial value of 4350 (consider using initial_event=True to detect this event)
50 events found
Event IDs: [8192]
50 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 50 events and 408 original time points ...
2 bad epochs dropped
Computing data rank from raw with rank=None
    Using tolerance 2.8e-09 (2.2e-16 eps * 245 dim * 5.1e+04  max singular value)
    Estimated rank (mag): 245
    MAG: rank 245 computed from 245 data channels with 0 projectors
Reducing data rank from 245 -> 245
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 6576
[done]
BEM               : <ConductorModel  |  Sphere (3 layers): r0=[0.0, 0.0, 0.0] R=80 mm>
MRI transform     : identity
Sphere model      : origin at (   0.00    0.00    0.00) mm, rad =    0.1 mm
Guess grid        :   20.0 mm
Guess mindist     :    5.0 mm
Guess exclude     :   20.0 mm
Using standard MEG coil definitions.

Coordinate transformation: MRI (surface RAS) -> head
     1.000000  0.000000  0.000000       0.00 mm
     0.000000  1.000000  0.000000       0.00 mm
     0.000000  0.000000  1.000000       0.00 mm
     0.000000  0.000000  0.000000       1.00
Coordinate transformation: MEG device -> head
     0.975557 -0.033946 -0.217110      -1.50 mm
     0.044391  0.998071  0.043409     -34.60 mm
     0.215218 -0.051986  0.975181      -6.32 mm
     0.000000  0.000000  0.000000       1.00
3 bad channels total
Read 245 MEG channels from info
Read  23 MEG compensation channels from info
84 coil definitions read
Coordinate transformation: MEG device -> head
     0.975557 -0.033946 -0.217110      -1.50 mm
     0.044391  0.998071  0.043409     -34.60 mm
     0.215218 -0.051986  0.975181      -6.32 mm
     0.000000  0.000000  0.000000       1.00
0 compensation data sets in info
MEG coil definitions created in head coordinates.
Decomposing the sensor noise covariance matrix...
Computing data rank from covariance with rank=None
    Using tolerance 2.1e-14 (2.2e-16 eps * 245 dim * 0.39  max singular value)
    Estimated rank (mag): 245
    MAG: rank 245 computed from 245 data channels with 0 projectors
    Setting small MAG eigenvalues to zero (without PCA)
    Created the whitener using a noise covariance matrix with rank 245 (0 small eigenvalues omitted)

---- Computing the forward solution for the guesses...
Making a spherical guess space with radius    72.0 mm...
Filtering (grid =     20 mm)...
Surface CM = (   0.0    0.0    0.0) mm
Surface fits inside a sphere with radius   72.0 mm
Surface extent:
    x =  -72.0 ...   72.0 mm
    y =  -72.0 ...   72.0 mm
    z =  -72.0 ...   72.0 mm
Grid extent:
    x =  -80.0 ...   80.0 mm
    y =  -80.0 ...   80.0 mm
    z =  -80.0 ...   80.0 mm
729 sources before omitting any.
178 sources after omitting infeasible sources.
170 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
Go through all guess source locations...
[done 170 sources]
---- Fitted :    69.3 ms
No projector specified for this dataset. Please consider the method self.add_proj.
1 time points fitted

Compute localisation errors

actual_pos = 0.01 * np.array([[0.16, 1.61, 5.13],
                              [0.17, 1.35, 4.15],
                              [0.16, 1.05, 3.19],
                              [0.13, 0.80, 2.26]])
actual_pos = np.dot(actual_pos, [[0, 1, 0], [-1, 0, 0], [0, 0, 1]])

errors = 1e3 * np.linalg.norm(actual_pos - pos, axis=1)
print("errors (mm) : %s" % errors)

Out:

errors (mm) : [1.374363 1.377723 1.257478 1.28907 ]

Plot the dipoles in 3D

def plot_pos(pos, color=(0., 0., 0.)):
    mlab.points3d(pos[:, 0], pos[:, 1], pos[:, 2], scale_factor=0.005,
                  color=color)


mne.viz.plot_alignment(evoked.info, bem=sphere, surfaces=[])
# Plot the position of the actual dipole
plot_pos(actual_pos, color=(1., 0., 0.))
# Plot the position of the estimated dipole
plot_pos(pos, color=(1., 1., 0.))
../../_images/sphx_glr_plot_phantom_4DBTi_005.png

Out:

Getting helmet for system Magnes_3600wh

Total running time of the script: ( 0 minutes 9.315 seconds)

Estimated memory usage: 8 MB

Gallery generated by Sphinx-Gallery