diff --git a/doc/_includes/forward.rst b/doc/_includes/forward.rst index 40706fb292c..2ed1ecacaab 100644 --- a/doc/_includes/forward.rst +++ b/doc/_includes/forward.rst @@ -689,9 +689,7 @@ EEG forward solution in the sphere model For examples of using the sphere model when computing the forward model (using :func:`mne.make_forward_solution`), see :ref:`Brainstorm CTF phantom - dataset tutorial `, - :ref:`Brainstorm Elekta phantom dataset tutorial - `, and + dataset tutorial ` and :ref:`tut-source-alignment-without-mri`. When the sphere model is employed, the computation of the EEG solution can be diff --git a/doc/changes/dev/13584.other.rst b/doc/changes/dev/13584.other.rst new file mode 100644 index 00000000000..6fd8196d231 --- /dev/null +++ b/doc/changes/dev/13584.other.rst @@ -0,0 +1 @@ +Adopting the diataxis framework by improving the phantom dipole tutorial, by `Carina Forster`_. diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index c76a5c5f569..129025d10ac 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -5,22 +5,32 @@ Brainstorm Elekta phantom dataset tutorial ========================================== -Here we compute the evoked from raw for the Brainstorm Elekta phantom -tutorial dataset. For comparison, see :footcite:`TadelEtAl2011` and +This tutorial provides a step-by-step guide to +importing and processing Elekta-Neuromag current phantom recordings. + +A phantom recording is a measurement obtained using a device (phantom) +that generates known magnetic signals, +allowing validation and benchmarking of MEG system accuracy and analysis methods. + +The aim of this tutorial is to demonstrate how phantom recordings can be used to +evaluate source localisation methods by comparing estimated and true dipole positions. + +For comparison, see :footcite:`TadelEtAl2011` and `the original Brainstorm tutorial `__. """ # sphinx_gallery_thumbnail_number = 9 # Authors: Eric Larson +# Carina Forster # # License: BSD-3-Clause # Copyright the MNE-Python contributors. # %% - import matplotlib.pyplot as plt import numpy as np +from scipy.signal import find_peaks import mne from mne import find_events, fit_dipole @@ -28,153 +38,211 @@ from mne.datasets.brainstorm import bst_phantom_elekta from mne.io import read_raw_fif -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 :class:`mne.io.Raw`. +# Load and prepare the data +# ------------------------- +# The data were collected with an Elekta Neuromag VectorView system +# at 1000 Hz, low-pass filtered at 330 Hz and contain recordings +# at three current amplitudes (20, 200, and 2000 nAm). +# Here we load the medium-amplitude condition. + +# Load the data data_path = bst_phantom_elekta.data_path(verbose=True) - raw_fname = data_path / "kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif" raw = read_raw_fif(raw_fname) -# %% -# 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: +# Mark known bad channels +raw.info["bads"] = ["MEG1933", "MEG2421"] +# The first 32 events correspond to dipole activations events = find_events(raw, "STI201") -raw.plot(events=events) -raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels + # %% -# The data has strong line frequency (60 Hz and harmonics) and cHPI coil -# noise (five peaks around 300 Hz). Here, we use only the first 30 seconds -# to save memory: +# Epoch the data and plot evokeds +# ------------------------------- +# We epoch and baseline correct the data around the dipole events. -raw.compute_psd(tmax=30).plot( - average=False, amplitude=False, picks="data", exclude="bads" +# Epoch the data +bmax = -0.05 +tmin, tmax = -0.1, 0.8 +event_id = list(range(1, 33)) +epochs = mne.Epochs( + raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False ) +# We drop the first and last event, it can contain dipole-switching artifacts +epochs_clean = epochs[1:-1] + +# We select the first simulated dipole for visualisation purposes +epochs_firstdip = epochs_clean["1"] + # %% -# Our phantom produces sinusoidal bursts at 20 Hz: +# Let's look at the evoked response for the first clean dipole. +# +# We can see that the phantom was set to produce 20 Hz sinusoidal bursts of current +# and the burst envelope repeats at approximately 3 Hz. -raw.plot(events=events) +epochs_firstdip.average().plot(time_unit="s") # %% -# 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. +# Determine peak activation using Global Field Power (GFP) +# -------------------------------------------------------- +# GFP is the standard deviation across sensors at each time +# point, providing a reference-independent measure of signal strength. -tmin, tmax = -0.1, 0.1 -bmax = -0.05 # Avoid capture filter ringing into baseline -event_id = list(range(1, 33)) -epochs = mne.Epochs( - raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False -) -epochs["1"].average().plot(time_unit="s") +# Get the evoked signal of the first dipole +evoked_tmp = epochs_firstdip.average() + +# Calculate GFP +gfp = np.std(evoked_tmp.data, axis=0) +# Restrict to first burst window +times = evoked_tmp.times +time_mask = (times > 0) & (times <= 0.05) + +# Find the peak GFP indices +peaks, _ = find_peaks(gfp[time_mask]) +peak_indices = np.where(time_mask)[0][peaks] + +# Select the strongest peak +strongest_peak_idx = peak_indices[np.argmax(gfp[peak_indices])] +t_peak = times[strongest_peak_idx] +print(f"Strongest peak at {t_peak * 1000:.1f} ms") # %% -# .. _plt_brainstorm_phantom_elekta_eeg_sphere_geometry: -# -# Let's use a :ref:`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.). -# -# Even though this is a VectorView/TRIUX phantom, we can use the Otaniemi -# phantom subject as a surrogate because the "head" surface (hemisphere outer -# shell) has the same geometry for both phantoms, even though the internal -# dipole locations differ. The phantom_otaniemi scan was aligned to the -# phantom's head coordinate frame, so an identity ``trans`` is appropriate -# here. +# Here we select the peak amplitude timepoint and store the evoked data for each dipole. + +evokeds = [] +for ii in event_id: + evoked = epochs_clean[str(ii)].average().crop(t_peak, t_peak) + evoked = mne.EvokedArray(np.array(evoked.data), evoked.info, tmin=0.0) + evokeds.append(evoked) +# %% +# Next, we need to compute the noise covariance in the baseline window +# to capture the sensor noise structure (for details: :ref:`tut-compute-covariance`). + +cov = mne.compute_covariance(epochs_clean, tmax=bmax) +del epochs # delete to save memory +# %% +# We use a :ref:`sphere head geometry model ` +# because the Elekta phantom is designed to approximate a spherical +# conductor with known dipole locations. subjects_dir = data_path fetch_phantom("otaniemi", subjects_dir=subjects_dir) sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08) -subject = "phantom_otaniemi" -trans = mne.transforms.Transform("head", "mri", np.eye(4)) -mne.viz.plot_alignment( - epochs.info, - subject=subject, - show_axes=True, - bem=sphere, - dig=True, - surfaces=("head-dense", "inner_skull"), - trans=trans, - mri_fiducials=True, - subjects_dir=subjects_dir, -) # %% -# 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. +# Dipole fitting +# -------------- +# Finally, we fit dipoles for each phantom and store them in a list. -# 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=bmax) -mne.viz.plot_evoked_white(epochs["1"].average(), cov) +dip_all = [] -data = [] -t_peak = 0.036 # true for Elekta phantom -for ii in event_id: - # Avoid the first and last trials -- can contain dipole-switching artifacts - evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak) - data.append(evoked.data[:, 0]) -evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.0) -del epochs -dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None) +for evoked in evokeds: + dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1) + dip_all.append(dip) +# %% +# Evaluate goodness of fit +# ------------------------ +# The dipole object stores the goodness of fit (GOF) for each dipole. +# Some dipoles have a low GOF (< 60 %). +gof = [dip.gof[0] for dip in dip_all] +colors = ["#E69F00" if val < 60 else "#0072B2" for val in gof] +plt.bar(event_id, gof, color=colors) +plt.xlabel("Phantom dipole estimation") +plt.ylabel("Goodness of fit (%)") +plt.show() # %% -# 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: - for text in list(ax.texts): - text.remove() - for line in ax.lines: - line.set_color("#98df81") -residual.plot(axes=axes) +# Dipoles with low goodness of fit +# -------------------------------- +# Why do some dipoles have a low GOF? +# Here we plot the dipole locations of the dipoles with low GOF. +# +# We can see that dipoles with low GOF are deep in the brain which might explain +# the low GOF. + +# Get indices of low GOF dipoles +low_idx = [i for i, g in enumerate(gof) if g < 60] +low_event_ids = [event_id[i] for i in low_idx] + +print("Low GOF dipoles:", low_event_ids) + +# Let's plot the locations of the dipoles with low GOF. +low_dips = [dip_all[i] for i in low_idx] + +subject = "phantom_otaniemi" +trans = mne.transforms.Transform("head", "mri", np.eye(4)) + +# Plot the position and the orientation of the dipoles with low GOF +fig = mne.viz.plot_alignment( + evoked.info, + trans, + subject, + bem=sphere, + surfaces={"head-dense": 0.2}, + coord_frame="head", + meg="helmet", + show_axes=True, + subjects_dir=subjects_dir, +) +fig = mne.viz.plot_dipole_locations( + dipoles=low_dips, mode="arrow", subject=subject, color=(1.0, 0.2, 0.2), fig=fig +) # %% -# Now we can compare to the actual locations, taking the difference in mm: +# Compare estimated and true dipoles +# ---------------------------------- +# The dipole fits closely match the true phantom data, +# achieving sub-centimeter accuracy (mean position error 2.4mm). +# We get the true dipole positions from the phantoms actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 100.0 # nAm +# Here we store the estimated dipoles +dip_pos = [dip.pos[0] for dip in dip_all] +dip_ori = [dip.ori[0] for dip in dip_all] +dip_amplitude = [dip.amplitude[0] for dip in dip_all] + fig, (ax1, ax2, ax3) = plt.subplots( nrows=3, ncols=1, figsize=(6, 7), layout="constrained" ) -diffs = 1000 * np.sqrt(np.sum((dip.pos - actual_pos) ** 2, axis=-1)) +# Here we calculate the euclidean distance between estimated and true positions. +# We multiply by 1000 to convert from meter to millimeter. +diffs = 1000 * np.sqrt(np.sum((dip_pos - actual_pos) ** 2, axis=-1)) print(f"mean(position error) = {np.mean(diffs):0.1f} mm") 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)))) +# Next we calculate the angle between estimated and true orientation. +# We convert radians to degrees. +angles = np.rad2deg(np.arccos(np.abs(np.sum(dip_ori * actual_ori, axis=1)))) print(f"mean(angle error) = {np.mean(angles):0.1f}°") ax2.bar(event_id, angles) ax2.set_xlabel("Dipole index") ax2.set_ylabel("Angle error (°)") -amps = actual_amp - dip.amplitude / 1e-9 +# Here we compare amplitudes by subtracting estimated from true amplitude. +amps = actual_amp - np.array(dip_amplitude) / 1e-9 print(f"mean(abs amplitude error) = {np.mean(np.abs(amps)):0.1f} nAm") ax3.bar(event_id, amps) ax3.set_xlabel("Dipole index") ax3.set_ylabel("Amplitude error (nAm)") - # %% -# Let's plot the positions and the orientations of the actual and the estimated -# dipoles +# Visualise estimated and true dipole locations +# --------------------------------------------- +# We can see that the dipoles overlap, have approximately the same magnitude +# and point in the same direction. actual_amp = np.ones(len(dip)) # fake amp, needed to create Dipole instance -actual_gof = np.ones(len(dip)) # fake GOF, needed to create Dipole instance +actual_gof = np.ones(len(dip)) # fake goodness-of-fit (GOF) +# setup dipole objects for true and estimated dipoles dip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof) +dip_estimated = mne.Dipole(dip.times, dip_pos, dip_amplitude, dip_ori, actual_gof) fig = mne.viz.plot_alignment( evoked.info, @@ -188,16 +256,15 @@ subjects_dir=subjects_dir, ) -# Plot the position and the orientation of the actual dipole +# Plot the position and the orientation of the true dipole in black fig = mne.viz.plot_dipole_locations( dipoles=dip_true, mode="arrow", subject=subject, color=(0.0, 0.0, 0.0), fig=fig ) -# Plot the position and the orientation of the estimated dipole +# Plot the position and the orientation of the estimated dipole in green fig = mne.viz.plot_dipole_locations( - dipoles=dip, mode="arrow", subject=subject, color=(0.2, 1.0, 0.5), fig=fig + dipoles=dip_estimated, mode="arrow", subject=subject, color=(0.2, 1.0, 0.5), fig=fig ) - mne.viz.set_3d_view(figure=fig, azimuth=70, elevation=80, distance=0.5) # %%