Maxwell filter data with movement compensation

Demonstrate movement compensation on simulated data. The simulated data contains bilateral activation of auditory cortices, repeated over 14 different head rotations (head center held fixed). See the following for details:

mne-tools/mne-misc-data

# Authors: Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD-3-Clause

import mne
from mne.preprocessing import maxwell_filter

print(__doc__)

data_path = mne.datasets.misc.data_path(verbose=True) / "movement"

head_pos = mne.chpi.read_head_pos(data_path / "simulated_quats.pos")
raw = mne.io.read_raw_fif(data_path / "simulated_movement_raw.fif")
raw_stat = mne.io.read_raw_fif(data_path / "simulated_stationary_raw.fif")

Opening raw data file /home/circleci/mne_data/MNE-misc-data/movement/simulated_movement_raw.fif...
    Range : 25800 ... 34208 =     42.956 ...    56.955 secs
Ready.
Opening raw data file /home/circleci/mne_data/MNE-misc-data/movement/simulated_stationary_raw.fif...
    Range : 25800 ... 34208 =     42.956 ...    56.955 secs
Ready.

Visualize the “subject” head movements. By providing the measurement information, the distance to the nearest sensor in each direction (e.g., left/right for the X direction, forward/backward for Y) can be shown in blue, and the destination (if given) shown in red.

mne.viz.plot_head_positions(
    head_pos, mode="traces", destination=raw.info["dev_head_t"], info=raw.info
)

This can also be visualized using a quiver.

mne.viz.plot_head_positions(
    head_pos, mode="field", destination=raw.info["dev_head_t"], info=raw.info
)

Getting helmet for system 306m

Process our simulated raw data (taking into account head movements).

# extract our resulting events
events = mne.find_events(raw, stim_channel="STI 014")
events[:, 2] = 1
raw.plot(events=events)

topo_kwargs = dict(times=[0, 0.1, 0.2], ch_type="mag", vlim=(-500, 500))

14 events found
Event IDs: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14]

First, take the average of stationary data (bilateral auditory patterns).

evoked_stat = mne.Epochs(raw_stat, events, 1, -0.2, 0.8).average()
fig = evoked_stat.plot_topomap(**topo_kwargs)
fig.suptitle("Stationary")

Not setting metadata
14 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] s
Applying baseline correction (mode: mean)
0 projection items activated
NOTE: pick_channels() is a legacy function. New code should use inst.pick(...).

Second, take a naive average, which averages across epochs that have been simulated to have different head positions and orientations, thereby spatially smearing the activity.

epochs = mne.Epochs(raw, events, 1, -0.2, 0.8)
evoked = epochs.average()
fig = evoked.plot_topomap(**topo_kwargs)
fig.suptitle("Moving: naive average")

Not setting metadata
14 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] s
Applying baseline correction (mode: mean)
0 projection items activated
NOTE: pick_channels() is a legacy function. New code should use inst.pick(...).

Third, use raw movement compensation (restores pattern).

raw_sss = maxwell_filter(raw, head_pos=head_pos)
evoked_raw_mc = mne.Epochs(raw_sss, events, 1, -0.2, 0.8).average()
fig = evoked_raw_mc.plot_topomap(**topo_kwargs)
fig.suptitle("Moving: movement compensated (raw)")

Maxwell filtering raw data
    No bad MEG channels
    Processing 203 gradiometers and 102 magnetometers
    Automatic origin fit: head of radius 91.0 mm
    Using origin -4.1, 16.0, 51.7 mm in the head frame
        Using 90/95 harmonic components for    0.000  (75/80 in, 15/15 out)
    Appending head position result channels and loading raw data from disk
    Processing 1 data chunk
        Using 87/95 harmonic components for    0.000  (72/80 in, 15/15 out)
        Using 88/95 harmonic components for    1.001  (73/80 in, 15/15 out)
        Using 90/95 harmonic components for    2.000  (75/80 in, 15/15 out)
        Using 88/95 harmonic components for    3.000  (73/80 in, 15/15 out)
        Using 88/95 harmonic components for    3.999  (73/80 in, 15/15 out)
        Using 88/95 harmonic components for    5.000  (73/80 in, 15/15 out)
        Using 89/95 harmonic components for    6.001  (74/80 in, 15/15 out)
        Using 93/95 harmonic components for    6.999  (78/80 in, 15/15 out)
        Using 88/95 harmonic components for    8.000  (73/80 in, 15/15 out)
        Using 91/95 harmonic components for    9.001  (76/80 in, 15/15 out)
        Using 93/95 harmonic components for   10.000  (78/80 in, 15/15 out)
        Using 93/95 harmonic components for   11.000  (78/80 in, 15/15 out)
        Using 89/95 harmonic components for   11.999  (74/80 in, 15/15 out)
        Using 88/95 harmonic components for   13.000  (73/80 in, 15/15 out)
        Used  14 head positions for    0.000 -   13.999 s (#1/1)
[done]
Not setting metadata
14 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] s
Applying baseline correction (mode: mean)
0 projection items activated
NOTE: pick_channels() is a legacy function. New code should use inst.pick(...).

Fourth, use evoked movement compensation. For these data, which contain very large rotations, it does not as cleanly restore the pattern.

evoked_evo_mc = mne.epochs.average_movements(epochs, head_pos=head_pos)
fig = evoked_evo_mc.plot_topomap(**topo_kwargs)
fig.suptitle("Moving: movement compensated (evoked)")

    Automatic origin fit: head of radius 91.0 mm
Aligning and averaging up to 14 epochs
    No bad MEG channels
    Processing 203 gradiometers and 102 magnetometers
    Processing epoch 1 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 2 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 3 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 4 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 5 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 6 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 7 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 8 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 9 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 10 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 11 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 12 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 13 (device location: 0.0, 0.0, 40.0 mm)
    Processing epoch 14 (device location: 0.0, 0.0, 40.0 mm)
NOTE: pick_channels() is a legacy function. New code should use inst.pick(...).
Created Evoked dataset from 14 epochs

Estimated memory usage: 10 MB