apply diataxis to tutorials/inverse/80_brainstorm_phantom_elekta.py

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 <plt_brainstorm_phantom_ctf_eeg_sphere_geometry>`,
-    :ref:`Brainstorm Elekta phantom dataset tutorial
-    <plt_brainstorm_phantom_elekta_eeg_sphere_geometry>`, and
+    dataset tutorial <plt_brainstorm_phantom_ctf_eeg_sphere_geometry>` and
     :ref:`tut-source-alignment-without-mri`.
 
 When the sphere model is employed, the computation of the EEG solution can be

doc/changes/dev/13584.other.rst

@@ -0,0 +1 @@
+Adopting the diataxis framework by improving the phantom dipole tutorial, by `Carina Forster`_.

tutorials/inverse/80_brainstorm_phantom_elekta.py

@@ -5,176 +5,244 @@
 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
 <https://neuroimage.usc.edu/brainstorm/Tutorials/PhantomElekta>`__.
 """
 # sphinx_gallery_thumbnail_number = 9
 
 # Authors: Eric Larson <larson.eric.d@gmail.com>
+#          Carina Forster <carinaforster0611@gmail.com>
 #
 # 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
 from mne.datasets import fetch_phantom
 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 <eeg_sphere_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 <eeg_sphere_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)
 
 # %%