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
|
"""
.. _tut-phantom-kit:
============================
KIT phantom dataset tutorial
============================
Here we read KIT data obtained from a phantom with 49 dipoles sequentially activated
with 2-cycle 11 Hz sinusoidal bursts to verify source localization accuracy.
"""
# Authors: Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
# %%
import mne_bids
import numpy as np
import mne
data_path = mne.datasets.phantom_kit.data_path()
actual_pos, actual_ori = mne.dipole.get_phantom_dipoles("oyama")
actual_pos, actual_ori = actual_pos[:49], actual_ori[:49] # only 49 of 50 dipoles
bids_path = mne_bids.BIDSPath(
root=data_path,
subject="01",
task="phantom",
run="01",
datatype="meg",
)
# ignore warning about misc units
raw = mne_bids.read_raw_bids(bids_path).load_data()
# Let's apply a little bit of preprocessing (temporal filtering and reference
# regression)
picks_artifact = ["MISC 001", "MISC 002", "MISC 003"]
picks = np.r_[
mne.pick_types(raw.info, meg=True),
mne.pick_channels(raw.info["ch_names"], picks_artifact),
]
raw.filter(None, 40, picks=picks)
mne.preprocessing.regress_artifact(
raw, picks="meg", picks_artifact=picks_artifact, copy=False, proj=False
)
plot_scalings = dict(mag=5e-12) # large-amplitude sinusoids
raw_plot_kwargs = dict(duration=15, n_channels=50, scalings=plot_scalings)
events, event_id = mne.events_from_annotations(raw)
raw.plot(events=events, **raw_plot_kwargs)
n_dip = len(event_id)
assert n_dip == 49 # sanity check
# %%
# We can also look at the power spectral density to see the phantom oscillations at
# 11 Hz plus the expected frequency-domain sinc-like oscillations due to the time-domain
# boxcar windowing of the 11 Hz sinusoid.
spectrum = raw.copy().crop(0, 60).compute_psd(n_fft=10000)
fig = spectrum.plot(amplitude=False)
fig.axes[0].set_xlim(0, 50)
dip_freq = 11.0
fig.axes[0].axvline(dip_freq, color="r", ls="--", lw=2, zorder=4)
# %%
# Now we can figure out our epoching parameters and epoch the data and plot it.
tmin, tmax = -0.08, 0.18
epochs = mne.Epochs(raw, tmin=tmin, tmax=tmax, decim=10, picks="data", preload=True)
del raw
epochs.plot(scalings=plot_scalings)
# %%
# Now we can average the epochs for each dipole, get the activation at the peak time,
# and create an :class:`mne.EvokedArray` from the result.
t_peak = 1.0 / dip_freq / 4.0
data = np.zeros((len(epochs.ch_names), n_dip))
for di in range(n_dip):
data[:, [di]] = epochs[f"dip{di + 1:02d}"].average().crop(t_peak, t_peak).data
evoked = mne.EvokedArray(data, epochs.info, tmin=0, comment="KIT phantom activations")
evoked.plot_joint()
# %%
# Let's fit dipoles at each dipole's peak activation time.
trans = mne.transforms.Transform("head", "mri", np.eye(4))
sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08)
cov = mne.compute_covariance(epochs, tmax=0, method="empirical")
dip, residual = mne.fit_dipole(evoked, cov, sphere, n_jobs=None)
# %%
# Finally let's look at the results.
# sphinx_gallery_thumbnail_number = 5
print(f"Average amplitude: {np.mean(dip.amplitude) * 1e9:0.1f} nAm")
print(f"Average GOF: {np.mean(dip.gof):0.1f}%")
diffs = 1000 * np.sqrt(np.sum((dip.pos - actual_pos) ** 2, axis=-1))
print(f"Average loc error: {np.mean(diffs):0.1f} mm")
angles = np.rad2deg(np.arccos(np.abs(np.sum(dip.ori * actual_ori, axis=1))))
print(f"Average ori error: {np.mean(angles):0.1f}°")
fig = mne.viz.plot_alignment(
evoked.info,
trans,
bem=sphere,
coord_frame="head",
meg="helmet",
show_axes=True,
)
fig = mne.viz.plot_dipole_locations(
dipoles=dip, mode="arrow", color=(0.2, 1.0, 0.5), fig=fig
)
actual_amp = np.ones(len(dip)) # misc amp to create Dipole instance
actual_gof = np.ones(len(dip)) # misc GOF to create Dipole instance
dip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof)
fig = mne.viz.plot_dipole_locations(
dipoles=dip_true, mode="arrow", color=(0.0, 0.0, 0.0), fig=fig
)
mne.viz.set_3d_view(figure=fig, azimuth=90, elevation=90, distance=0.5)
|
