Note
Click here to download the full example code
Here we compute the evoked from raw for the Brainstorm Elekta phantom tutorial dataset. For comparison, see [1] and:
[1] | Tadel F, Baillet S, Mosher JC, Pantazis D, Leahy RM. Brainstorm: A User-Friendly Application for MEG/EEG Analysis. Computational Intelligence and Neuroscience, vol. 2011, Article ID 879716, 13 pages, 2011. doi:10.1155/2011/879716 |
# sphinx_gallery_thumbnail_number = 9
# Authors: Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD (3-clause)
import os.path as op
import numpy as np
import matplotlib.pyplot as plt
import mne
from mne import find_events, fit_dipole
from mne.datasets.brainstorm import bst_phantom_elekta
from mne.io import read_raw_fif
from mayavi import mlab
print(__doc__)
The data were collected with an Elekta Neuromag VectorView system at 1000 Hz
and low-pass filtered at 330 Hz. Here the medium-amplitude (200 nAm) data
are read to construct instances of mne.io.Raw
.
data_path = bst_phantom_elekta.data_path(verbose=True)
raw_fname = op.join(data_path, 'kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif')
raw = read_raw_fif(raw_fname)
Out:
Opening raw data file /home/circleci/mne_data/MNE-brainstorm-data/bst_phantom_elekta/kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif...
Read a total of 13 projection items:
planar-0.0-115.0-PCA-01 (1 x 306) idle
planar-0.0-115.0-PCA-02 (1 x 306) idle
planar-0.0-115.0-PCA-03 (1 x 306) idle
planar-0.0-115.0-PCA-04 (1 x 306) idle
planar-0.0-115.0-PCA-05 (1 x 306) idle
axial-0.0-115.0-PCA-01 (1 x 306) idle
axial-0.0-115.0-PCA-02 (1 x 306) idle
axial-0.0-115.0-PCA-03 (1 x 306) idle
axial-0.0-115.0-PCA-04 (1 x 306) idle
axial-0.0-115.0-PCA-05 (1 x 306) idle
axial-0.0-115.0-PCA-06 (1 x 306) idle
axial-0.0-115.0-PCA-07 (1 x 306) idle
axial-0.0-115.0-PCA-08 (1 x 306) idle
Range : 47000 ... 437999 = 47.000 ... 437.999 secs
Ready.
Current compensation grade : 0
Data channel array consisted of 204 MEG planor gradiometers, 102 axial magnetometers, and 3 stimulus channels. Let’s get the events for the phantom, where each dipole (1-32) gets its own event:
events = find_events(raw, 'STI201')
raw.plot(events=events)
raw.info['bads'] = ['MEG2421']
Out:
645 events found
Event IDs: [ 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 256 768 1792 3840 7936]
The data have strong line frequency (60 Hz and harmonics) and cHPI coil noise (five peaks around 300 Hz). Here we plot only out to 60 seconds to save memory:
raw.plot_psd(tmax=60., average=False)
Out:
Effective window size : 2.048 (s)
Effective window size : 2.048 (s)
Let’s use Maxwell filtering to clean the data a bit. Ideally we would have the fine calibration and cross-talk information for the site of interest, but we don’t, so we just do:
raw.fix_mag_coil_types()
raw = mne.preprocessing.maxwell_filter(raw, origin=(0., 0., 0.))
Out:
101 of 101 T1/T2 magnetometer types replaced with T3.
Maxwell filtering raw data
Loading raw data from disk
Bad MEG channels being reconstructed: ['MEG2421']
Processing 204 gradiometers and 102 magnetometers
Using origin 0.0, 0.0, 0.0 mm in the head frame
Using 87/95 harmonic components for 0.000 (72/80 in, 15/15 out)
Processing 39 data chunks
[done]
We know our phantom produces sinusoidal bursts below 25 Hz, so let’s filter.
raw.filter(None, 40., fir_design='firwin')
raw.plot(events=events)
Out:
Setting up low-pass filter at 40 Hz
h_trans_bandwidth chosen to be 10.0 Hz
Filter length of 331 samples (0.331 sec) selected
Now we epoch our data, average it, and look at the first dipole response. The first peak appears around 3 ms. Because we low-passed at 40 Hz, we can also decimate our data to save memory.
tmin, tmax = -0.1, 0.1
event_id = list(range(1, 33))
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, baseline=(None, -0.01),
decim=3, preload=True)
epochs['1'].average().plot(time_unit='s')
Out:
640 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 640 events and 201 original time points ...
0 bad epochs dropped
Let’s use a sphere head geometry model and let’s see the coordinate alignment and the sphere location. The phantom is properly modeled by a single-shell sphere with origin (0., 0., 0.).
sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=0.08)
mne.viz.plot_alignment(raw.info, subject='sample', show_axes=True,
bem=sphere, dig=True, surfaces='inner_skull')
Out:
Equiv. model fitting -> RV = 0.00347649 %
mu1 = 0.94483 lambda1 = 0.136877
mu2 = 0.667752 lambda2 = 0.683647
mu3 = -0.295407 lambda3 = -0.0101536
Set up EEG sphere model with scalp radius 80.0 mm
Triangle neighbors and vertex normals...
Getting helmet for system 306m
Let’s do some dipole fits. We first compute the noise covariance, then do the fits for each event_id taking the time instant that maximizes the global field power.
# here we can get away with using method='oas' for speed (faster than "shrunk")
# but in general "shrunk" is usually better
cov = mne.compute_covariance(
epochs, tmax=0, method='oas', rank=None)
mne.viz.plot_evoked_white(epochs['1'].average(), cov)
data = []
t_peak = 0.036 # true for Elekta phantom
for ii in event_id:
evoked = epochs[str(ii)].average().crop(t_peak, t_peak)
data.append(evoked.data[:, 0])
evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.)
del epochs, raw
dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1)
Out:
estimated rank (mag + grad): 72
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Reducing data rank from 306 -> 72
Estimating covariance using OAS
Done.
Number of samples used : 21760
[done]
0 projection items activated
estimated rank (mag + grad): 72
0 projection items activated
estimated rank (mag + grad): 72
SSS has been applied to data. Showing mag and grad whitening jointly.
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Created the whitener using a noise covariance matrix with rank 72 (234 small eigenvalues omitted)
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.976295 -0.211976 0.043756 0.29 mm
0.206488 0.972764 0.105326 0.57 mm
-0.064891 -0.093794 0.993475 5.41 mm
0.000000 0.000000 0.000000 1.00
0 bad channels total
Read 306 MEG channels from info
84 coil definitions read
Coordinate transformation: MEG device -> head
0.976295 -0.211976 0.043756 0.29 mm
0.206488 0.972764 0.105326 0.57 mm
-0.064891 -0.093794 0.993475 5.41 mm
0.000000 0.000000 0.000000 1.00
MEG coil definitions created in head coordinates.
Decomposing the sensor noise covariance matrix...
estimated rank (mag + grad): 72
Setting small MEG eigenvalues to zero.
Not doing PCA for MEG.
Created the whitener using a noise covariance matrix with rank 72 (234 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 : 0.0 ms
---- Fitted : 3.0 ms
---- Fitted : 6.0 ms
---- Fitted : 9.0 ms
---- Fitted : 12.0 ms
---- Fitted : 15.0 ms
---- Fitted : 18.0 ms
---- Fitted : 21.0 ms
---- Fitted : 24.0 ms
---- Fitted : 27.0 ms
---- Fitted : 30.0 ms
---- Fitted : 33.0 ms
---- Fitted : 36.0 ms
---- Fitted : 39.0 ms
---- Fitted : 42.0 ms
---- Fitted : 45.0 ms
---- Fitted : 48.0 ms
---- Fitted : 51.0 ms
---- Fitted : 54.0 ms
---- Fitted : 57.0 ms
---- Fitted : 60.0 ms
---- Fitted : 63.0 ms
---- Fitted : 66.0 ms
---- Fitted : 69.0 ms
---- Fitted : 72.0 ms
---- Fitted : 75.0 ms
---- Fitted : 78.0 ms
---- Fitted : 81.0 ms
---- Fitted : 84.0 ms
---- Fitted : 87.0 ms
---- Fitted : 90.0 ms
---- Fitted : 93.0 ms
No projector specified for this dataset. Please consider the method self.add_proj.
32 time points fitted
Do a quick visualization of how much variance we explained, putting the data and residuals on the same scale (here the “time points” are the 32 dipole peak values that we fit):
fig, axes = plt.subplots(2, 1)
evoked.plot(axes=axes)
for ax in axes:
ax.texts = []
for line in ax.lines:
line.set_color('#98df81')
residual.plot(axes=axes)
Now we can compare to the actual locations, taking the difference in mm:
actual_pos, actual_ori = mne.dipole.get_phantom_dipoles()
actual_amp = 100. # nAm
fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1, figsize=(6, 7))
diffs = 1000 * np.sqrt(np.sum((dip.pos - actual_pos) ** 2, axis=-1))
print('mean(position error) = %0.1f mm' % (np.mean(diffs),))
ax1.bar(event_id, diffs)
ax1.set_xlabel('Dipole index')
ax1.set_ylabel('Loc. error (mm)')
angles = np.rad2deg(np.arccos(np.abs(np.sum(dip.ori * actual_ori, axis=1))))
print(u'mean(angle error) = %0.1f°' % (np.mean(angles),))
ax2.bar(event_id, angles)
ax2.set_xlabel('Dipole index')
ax2.set_ylabel(u'Angle error (°)')
amps = actual_amp - dip.amplitude / 1e-9
print('mean(abs amplitude error) = %0.1f nAm' % (np.mean(np.abs(amps)),))
ax3.bar(event_id, amps)
ax3.set_xlabel('Dipole index')
ax3.set_ylabel('Amplitude error (nAm)')
fig.tight_layout()
plt.show()
Out:
mean(position error) = 2.2 mm
mean(angle error) = 1.7°
mean(abs amplitude error) = 13.9 nAm
Let’s plot the positions and the orientations of the actual and the estimated dipoles
def plot_pos_ori(pos, ori, color=(0., 0., 0.), opacity=1.):
x, y, z = pos.T
u, v, w = ori.T
mlab.points3d(x, y, z, scale_factor=0.005, opacity=opacity, color=color)
q = mlab.quiver3d(x, y, z, u, v, w,
scale_factor=0.03, opacity=opacity,
color=color, mode='arrow')
q.glyph.glyph_source.glyph_source.shaft_radius = 0.02
q.glyph.glyph_source.glyph_source.tip_length = 0.1
q.glyph.glyph_source.glyph_source.tip_radius = 0.05
mne.viz.plot_alignment(evoked.info, bem=sphere, surfaces='inner_skull',
coord_frame='head', meg='helmet', show_axes=True)
# Plot the position and the orientation of the actual dipole
plot_pos_ori(actual_pos, actual_ori, color=(0., 0., 0.), opacity=0.5)
# Plot the position and the orientation of the estimated dipole
plot_pos_ori(dip.pos, dip.ori, color=(0.2, 1., 0.5))
mlab.view(70, 80, distance=0.5)
Out:
Triangle neighbors and vertex normals...
Getting helmet for system 306m
Total running time of the script: ( 0 minutes 36.916 seconds)
Estimated memory usage: 932 MB