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
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
"""
.. _ex-opm-somatosensory:
Optically pumped magnetometer (OPM) data
========================================
In this dataset, electrical median nerve stimulation was delivered to the
left wrist of the subject. Somatosensory evoked fields were measured using
nine QuSpin SERF OPMs placed over the right-hand side somatomotor area. Here
we demonstrate how to localize these custom OPM data in MNE.
"""
# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
# sphinx_gallery_thumbnail_number = 4
import numpy as np
import mne
data_path = mne.datasets.opm.data_path()
subject = "OPM_sample"
subjects_dir = data_path / "subjects"
raw_fname = data_path / "MEG" / "OPM" / "OPM_SEF_raw.fif"
bem_fname = subjects_dir / subject / "bem" / f"{subject}-5120-5120-5120-bem-sol.fif"
fwd_fname = data_path / "MEG" / "OPM" / "OPM_sample-fwd.fif"
coil_def_fname = data_path / "MEG" / "OPM" / "coil_def.dat"
# %%
# Prepare data for localization
# -----------------------------
# First we filter and epoch the data:
raw = mne.io.read_raw_fif(raw_fname, preload=True)
raw.filter(None, 90, h_trans_bandwidth=10.0)
raw.notch_filter(50.0, notch_widths=1)
# Set epoch rejection threshold a bit larger than for SQUIDs
reject = dict(mag=2e-10)
tmin, tmax = -0.5, 1
# Find median nerve stimulator trigger
event_id = dict(Median=257)
events = mne.find_events(raw, stim_channel="STI101", mask=257, mask_type="and")
picks = mne.pick_types(raw.info, meg=True, eeg=False)
# We use verbose='error' to suppress warning about decimation causing aliasing,
# ideally we would low-pass and then decimate instead
epochs = mne.Epochs(
raw,
events,
event_id,
tmin,
tmax,
verbose="error",
reject=reject,
picks=picks,
proj=False,
decim=10,
preload=True,
)
evoked = epochs.average()
evoked.plot()
cov = mne.compute_covariance(epochs, tmax=0.0)
del epochs, raw
# %%
# Examine our coordinate alignment for source localization and compute a
# forward operator:
#
# .. note:: The Head<->MRI transform is an identity matrix, as the
# co-registration method used equates the two coordinate
# systems. This mis-defines the head coordinate system
# (which should be based on the LPA, Nasion, and RPA)
# but should be fine for these analyses.
bem = mne.read_bem_solution(bem_fname)
trans = mne.transforms.Transform("head", "mri") # identity transformation
# To compute the forward solution, we must
# provide our temporary/custom coil definitions, which can be done as::
#
# with mne.use_coil_def(coil_def_fname):
# fwd = mne.make_forward_solution(
# raw.info, trans, src, bem, eeg=False, mindist=5.0,
# n_jobs=None, verbose=True)
fwd = mne.read_forward_solution(fwd_fname)
# use fixed orientation here just to save memory later
mne.convert_forward_solution(fwd, force_fixed=True, copy=False)
with mne.use_coil_def(coil_def_fname):
fig = mne.viz.plot_alignment(
evoked.info,
trans=trans,
subject=subject,
subjects_dir=subjects_dir,
surfaces=("head", "pial"),
bem=bem,
)
mne.viz.set_3d_view(
figure=fig, azimuth=45, elevation=60, distance=0.4, focalpoint=(0.02, 0, 0.04)
)
# %%
# Perform dipole fitting
# ----------------------
# Fit dipoles on a subset of time points
with mne.use_coil_def(coil_def_fname):
dip_opm, _ = mne.fit_dipole(
evoked.copy().crop(0.040, 0.080), cov, bem, trans, verbose=True
)
idx = np.argmax(dip_opm.gof)
print(
f"Best dipole at t={1000 * dip_opm.times[idx]:0.1f} ms with "
f"{dip_opm.gof[idx]:0.1f}% GOF"
)
# Plot N20m dipole as an example
dip_opm.plot_locations(trans, subject, subjects_dir, mode="orthoview", idx=idx)
# %%
# Perform minimum-norm localization
# ---------------------------------
# Due to the small number of sensors, there will be some leakage of activity
# to areas with low/no sensitivity. Constraining the source space to
# areas we are sensitive to might be a good idea.
inverse_operator = mne.minimum_norm.make_inverse_operator(
evoked.info, fwd, cov, loose=0.0, depth=None
)
del fwd, cov
method = "MNE"
snr = 3.0
lambda2 = 1.0 / snr**2
stc = mne.minimum_norm.apply_inverse(
evoked, inverse_operator, lambda2, method=method, pick_ori=None, verbose=True
)
# Plot source estimate at time of best dipole fit
brain = stc.plot(
hemi="rh",
views="lat",
subjects_dir=subjects_dir,
initial_time=dip_opm.times[idx],
clim=dict(kind="percent", lims=[99, 99.9, 99.99]),
size=(400, 300),
background="w",
)
|
