From d79b383e05ca49b833cf21896b2a787f06d7df7f Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 11 Mar 2024 03:30:57 +0000 Subject: [PATCH 01/37] Bump davidslusser/actions_python_bandit from 1.0.0 to 1.0.1 Bumps [davidslusser/actions_python_bandit](https://github.com/davidslusser/actions_python_bandit) from 1.0.0 to 1.0.1. - [Release notes](https://github.com/davidslusser/actions_python_bandit/releases) - [Commits](https://github.com/davidslusser/actions_python_bandit/compare/v1.0.0...v1.0.1) --- updated-dependencies: - dependency-name: davidslusser/actions_python_bandit dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 5abfae43a3b..908555af797 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -30,7 +30,7 @@ jobs: needs: style runs-on: ubuntu-latest steps: - - uses: davidslusser/actions_python_bandit@v1.0.0 + - uses: davidslusser/actions_python_bandit@v1.0.1 with: src: "mne" options: "-c pyproject.toml -ll -r" From 25e39909c57c5bb70b3c7688e1803be3d105b2bd Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 10 Jan 2026 15:00:35 +1000 Subject: [PATCH 02/37] diataxis --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index c76a5c5f569..1e2c69b9ad0 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -6,7 +6,9 @@ ========================================== Here we compute the evoked from raw for the Brainstorm Elekta phantom -tutorial dataset. For comparison, see :footcite:`TadelEtAl2011` and +tutorial dataset. We will compare the actual dipole positions +to the estimated dipole positions. +For comparison, see :footcite:`TadelEtAl2011` and `the original Brainstorm tutorial `__. """ From 96245171fcbc2343c1052b49a6462dbe3c07a0b4 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 24 Jan 2026 15:38:06 +1000 Subject: [PATCH 03/37] first draft using diataxis framework --- .../inverse/80_brainstorm_phantom_elekta.py | 92 ++++++++++++------- 1 file changed, 60 insertions(+), 32 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 1e2c69b9ad0..22f531191d1 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -5,11 +5,12 @@ Brainstorm Elekta phantom dataset tutorial ========================================== -Here we compute the evoked from raw for the Brainstorm Elekta phantom -tutorial dataset. We will compare the actual dipole positions -to the estimated dipole positions. +This tutorial provides a step-by-step guide to +importing and processing Elekta-Neuromag current phantom recordings. +The aim of this tutorial is to show the user how to use phantom recordings to +evaluate source localisation methods by comparing estimated vs real dipole positions. For comparison, see :footcite:`TadelEtAl2011` and -`the original Brainstorm tutorial +`the original Brainstorm tutorial with an explanation of phantom recordings `__. """ # sphinx_gallery_thumbnail_number = 9 @@ -33,41 +34,51 @@ 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`. +# 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, amplitudes can be seen in raw data) +# data are read to construct instances of :class:`mne.io.Raw`. 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: +# Let's first get an idea of the data structure we are working with: +raw.info + +# The data channel array consisted of 204 MEG planor gradiometers, +# 102 axial magnetometers, and 3 stimulus channels. + +# Next, let's look at the events in the phantom data for one stimulus channel: events = find_events(raw, "STI201") raw.plot(events=events) raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels +# we have 32 artificial dipoles stored as events +# the remaining event IDs are not relevant for this tutorial (256, 768 and so on) # %% -# 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: +# Let's explore the phantom data by plotting power spectral density for each sensor. +# Here, we use only the first 30 seconds to save memory raw.compute_psd(tmax=30).plot( average=False, amplitude=False, picks="data", exclude="bads" ) +# We can see that the data has strong line frequency (60 Hz and harmonics) +# and cHPI coil noise (five peaks around 300 Hz). # %% -# Our phantom produces sinusoidal bursts at 20 Hz: +# Next we plot the dipole events raw.plot(events=events) +# We can see that the simulated dipoles produce sinusoidal bursts at 20 Hz +# (can we really see that in the plot?) # %% -# 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. +# Next, we epoch the the data based on the dipoles events (1:32) +# We select 100 ms before and 100ms after the event trigger +# and baseline correct the epochs from -100 ms to - 0.05 before stimulus onset tmin, tmax = -0.1, 0.1 bmax = -0.05 # Avoid capture filter ringing into baseline @@ -75,13 +86,19 @@ epochs = mne.Epochs( raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False ) + +# Here we average the epochs for the first simulated dipole +# and plot the evoked signal epochs["1"].average().plot(time_unit="s") +# we averaged over 640 simulated events for the first dipole +# We can see that the first peak in the data appears close to the trigger onset # %% # .. _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 +# Finally, we source reconstruct the evoked simulated dipole. +# To do this we use a :ref:`sphere head geometry model ` +# and visualise 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 @@ -107,17 +124,22 @@ mri_fiducials=True, subjects_dir=subjects_dir, ) - +# What does this plot tell us? # %% -# 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. +# Let's do some dipole fits. +# To do this we compute the noise covariance for each epoch. +# We plot the whitened data to assess the whitening step. # 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) +# Not sure what we see here TBH + +# Next, we fit the dipoles for the evoked data. +# We choose the timepoint which maximises global field power + data = [] t_peak = 0.036 # true for Elekta phantom for ii in event_id: @@ -129,9 +151,10 @@ dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None) # %% -# 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): +# Let's visualize the explained variance. +# To do this, we need to make sure that the +# data and the residuals are 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) @@ -142,9 +165,10 @@ line.set_color("#98df81") residual.plot(axes=axes) +# again not sure how to interpret those plots # %% -# Now we can compare to the actual locations, taking the difference in mm: - +# Finally, we compare the estimated to the true dipole locations +# To do this, we calculate the difference by ..... actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 100.0 # nAm @@ -170,9 +194,11 @@ ax3.set_xlabel("Dipole index") ax3.set_ylabel("Amplitude error (nAm)") +# We can see that the error magnitude depends on the position of the estimate dipole +# however, the location error is never greater than 5 mm which is good? # %% -# Let's plot the positions and the orientations of the actual and the estimated -# dipoles +# Finally, we can plot the positions and the orientations +# of the actual and the estimated dipoles 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 @@ -190,18 +216,20 @@ subjects_dir=subjects_dir, ) -# Plot the position and the orientation of the actual dipole +# Plot the position and the orientation of the actual 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 ) mne.viz.set_3d_view(figure=fig, azimuth=70, elevation=80, distance=0.5) +# Here the green arrows represent the estimated dipoles, the black arrows +# the true dipole location # %% # References # ---------- From 174d36d014ac787e18aad48e816a3bd50c57ae3c Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 25 Jan 2026 11:19:35 +1000 Subject: [PATCH 04/37] added changelog --- doc/changes/dev/13585.other.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 doc/changes/dev/13585.other.rst diff --git a/doc/changes/dev/13585.other.rst b/doc/changes/dev/13585.other.rst new file mode 100644 index 00000000000..ae5a5a411c8 --- /dev/null +++ b/doc/changes/dev/13585.other.rst @@ -0,0 +1 @@ +Refactoring tutorial on Phantom dipoles using the diataxis framework, by `Carina Forster`_. From aa4ca199312eb76efbe35348640dcd98a5e2e59a Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 25 Jan 2026 11:25:09 +1000 Subject: [PATCH 05/37] changelog update --- doc/changes/dev/13585.other.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/changes/dev/13585.other.rst b/doc/changes/dev/13585.other.rst index ae5a5a411c8..6fd8196d231 100644 --- a/doc/changes/dev/13585.other.rst +++ b/doc/changes/dev/13585.other.rst @@ -1 +1 @@ -Refactoring tutorial on Phantom dipoles using the diataxis framework, by `Carina Forster`_. +Adopting the diataxis framework by improving the phantom dipole tutorial, by `Carina Forster`_. From 37ef2ad06a52d9869645eb9486397fa595a4c198 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 25 Jan 2026 11:35:47 +1000 Subject: [PATCH 06/37] wrong PR --- doc/changes/dev/13584.other.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 doc/changes/dev/13584.other.rst 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`_. From e31c386e52f9748cbc10fa5accefc79d1a1a47e0 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 25 Jan 2026 11:39:22 +1000 Subject: [PATCH 07/37] deleted changelog with wrong PR --- doc/changes/dev/13585.other.rst | 1 - 1 file changed, 1 deletion(-) delete mode 100644 doc/changes/dev/13585.other.rst diff --git a/doc/changes/dev/13585.other.rst b/doc/changes/dev/13585.other.rst deleted file mode 100644 index 6fd8196d231..00000000000 --- a/doc/changes/dev/13585.other.rst +++ /dev/null @@ -1 +0,0 @@ -Adopting the diataxis framework by improving the phantom dipole tutorial, by `Carina Forster`_. From ab6f3c5dd01a9b3b5a0f4f8d0e1457a42df6b2f0 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 14 Feb 2026 15:23:59 +1000 Subject: [PATCH 08/37] added explanations --- .../inverse/80_brainstorm_phantom_elekta.py | 51 +++++++++++-------- 1 file changed, 31 insertions(+), 20 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 22f531191d1..ec965e604ec 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -37,7 +37,7 @@ # 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, amplitudes can be seen in raw data) -# data are read to construct instances of :class:`mne.io.Raw`. +# data are accessed to construct instances of :class:`mne.io.Raw`. data_path = bst_phantom_elekta.data_path(verbose=True) raw_fname = data_path / "kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif" @@ -48,7 +48,7 @@ raw.info -# The data channel array consisted of 204 MEG planor gradiometers, +# The data consists of 204 MEG planor gradiometers, # 102 axial magnetometers, and 3 stimulus channels. # Next, let's look at the events in the phantom data for one stimulus channel: @@ -66,15 +66,14 @@ average=False, amplitude=False, picks="data", exclude="bads" ) -# We can see that the data has strong line frequency (60 Hz and harmonics) -# and cHPI coil noise (five peaks around 300 Hz). +# We can see that the data has strong line frequency (60 Hz, 120 Hz and so on) noise +# and cHPI (continuous head position indicator) coil noise (five peaks around 300 Hz). # %% # Next we plot the dipole events -raw.plot(events=events) +raw.plot(events=events, n_channels=15) # We can see that the simulated dipoles produce sinusoidal bursts at 20 Hz -# (can we really see that in the plot?) # %% # Next, we epoch the the data based on the dipoles events (1:32) # We select 100 ms before and 100ms after the event trigger @@ -92,6 +91,8 @@ epochs["1"].average().plot(time_unit="s") # we averaged over 640 simulated events for the first dipole # We can see that the first peak in the data appears close to the trigger onset +# at around 3ms, with a peak repeating every 3ms. +# Thus the burst repetition rate is 3Hz. # %% # .. _plt_brainstorm_phantom_elekta_eeg_sphere_geometry: @@ -124,21 +125,21 @@ mri_fiducials=True, subjects_dir=subjects_dir, ) -# What does this plot tell us? +# We can see that our head model aligns with the phantom head model. # %% # Let's do some dipole fits. -# To do this we compute the noise covariance for each epoch. -# We plot the whitened data to assess the whitening step. - -# here we can get away with using method='oas' for speed (faster than "shrunk") -# but in general "shrunk" is usually better +# To do this we compute the noise covariance for the window before dipole onset. cov = mne.compute_covariance(epochs, tmax=bmax) +# The covariance captures the sensor noise structure. +# We whiten the data to normalize noise across sensors before fitting dipoles +# tutorial on whitening/covariance estimation for details :ref:`tut-compute-covariance`. mne.viz.plot_evoked_white(epochs["1"].average(), cov) - -# Not sure what we see here TBH +# The plot shows the evoked signal divided by the estimated noise standard deviation. +# After whitening, most baseline activity should fall roughly within ±1 (unit variance). # Next, we fit the dipoles for the evoked data. # We choose the timepoint which maximises global field power +# We have seen in the evoked plot that this is around 3 ms after dipole onset. data = [] t_peak = 0.036 # true for Elekta phantom @@ -147,7 +148,7 @@ 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 +del epochs # save memory dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None) # %% @@ -165,43 +166,53 @@ line.set_color("#98df81") residual.plot(axes=axes) -# again not sure how to interpret those plots +# Here we visualise how well the dipole explains the evoked response (green line). +# The red lines represent the residuals, the leftover noise after dipole fitting. +# A good fit: green lines are strong and residuals are small and roughly flat. # %% # Finally, we compare the estimated to the true dipole locations # To do this, we calculate the difference by ..... +# This is our ground truth. actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 100.0 # nAm fig, (ax1, ax2, ax3) = plt.subplots( nrows=3, ncols=1, figsize=(6, 7), layout="constrained" ) - +# 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)") +# 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 (°)") +# Finally we compare amplitudes by subtracting estimated from true amplitude. amps = actual_amp - 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)") -# We can see that the error magnitude depends on the position of the estimate dipole -# however, the location error is never greater than 5 mm which is good? +# The dipole fits closely match the phantom ground truth. +# We can achieve sub-centimeter accuracy with a mean position error of 2.6 mm. +# This demonstrates that the fitting procedure is accurate. # %% # Finally, we can plot the positions and the orientations # of the actual and the estimated dipoles 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), needed to create Dipole instance dip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof) fig = mne.viz.plot_alignment( From c398ec76fce4f63d2e420cf03ea3c1ef1919f4a0 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 14 Feb 2026 18:12:49 +1000 Subject: [PATCH 09/37] formatting for docs --- .../inverse/80_brainstorm_phantom_elekta.py | 75 ++++++++++--------- 1 file changed, 40 insertions(+), 35 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index ec965e604ec..fbd0afa37c8 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -9,6 +9,7 @@ importing and processing Elekta-Neuromag current phantom recordings. The aim of this tutorial is to show the user how to use phantom recordings to evaluate source localisation methods by comparing estimated vs real dipole positions. + For comparison, see :footcite:`TadelEtAl2011` and `the original Brainstorm tutorial with an explanation of phantom recordings `__. @@ -39,7 +40,6 @@ # Here the medium-amplitude (200 nAm, amplitudes can be seen in raw data) # data are accessed to construct instances of :class:`mne.io.Raw`. 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) @@ -48,12 +48,12 @@ raw.info +# %% # The data consists of 204 MEG planor gradiometers, # 102 axial magnetometers, and 3 stimulus channels. -# Next, let's look at the events in the phantom data for one stimulus channel: +# Next, let's look at the events in the phantom data for one stimulus channel. events = find_events(raw, "STI201") -raw.plot(events=events) raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels # we have 32 artificial dipoles stored as events @@ -66,18 +66,19 @@ average=False, amplitude=False, picks="data", exclude="bads" ) -# We can see that the data has strong line frequency (60 Hz, 120 Hz and so on) noise -# and cHPI (continuous head position indicator) coil noise (five peaks around 300 Hz). +# We can see that the data has strong line frequency (60 Hz, 120 Hz ...) noise +# and cHPI (continuous head position indicator) coil noise (peaks around 300 Hz). # %% # Next we plot the dipole events -raw.plot(events=events, n_channels=15) +raw.plot(events=events, n_channels=10) # We can see that the simulated dipoles produce sinusoidal bursts at 20 Hz # %% -# Next, we epoch the the data based on the dipoles events (1:32) -# We select 100 ms before and 100ms after the event trigger -# and baseline correct the epochs from -100 ms to - 0.05 before stimulus onset +# Next, we epoch the the data based on the dipoles events (1:32). + +# We select 100 ms before and after the event trigger +# and baseline correct the epochs from -100 ms to -0.05 ms before stimulus onset. tmin, tmax = -0.1, 0.1 bmax = -0.05 # Avoid capture filter ringing into baseline @@ -89,8 +90,10 @@ # Here we average the epochs for the first simulated dipole # and plot the evoked signal epochs["1"].average().plot(time_unit="s") -# we averaged over 640 simulated events for the first dipole -# We can see that the first peak in the data appears close to the trigger onset +# %% + +# We averaged over 640 simulated events for the first dipole. +# The first peak in the data appears close to the trigger onset # at around 3ms, with a peak repeating every 3ms. # Thus the burst repetition rate is 3Hz. @@ -125,22 +128,23 @@ mri_fiducials=True, subjects_dir=subjects_dir, ) +# %% # We can see that our head model aligns with the phantom head model. # %% # Let's do some dipole fits. -# To do this we compute the noise covariance for the window before dipole onset. -cov = mne.compute_covariance(epochs, tmax=bmax) + +# First we compute the noise covariance for the baseline window. +# See covariance/whitening tutorial for details :ref:`tut-compute-covariance`. +# %% # The covariance captures the sensor noise structure. -# We whiten the data to normalize noise across sensors before fitting dipoles -# tutorial on whitening/covariance estimation for details :ref:`tut-compute-covariance`. -mne.viz.plot_evoked_white(epochs["1"].average(), cov) +cov = mne.compute_covariance(epochs, tmax=bmax) + # The plot shows the evoked signal divided by the estimated noise standard deviation. -# After whitening, most baseline activity should fall roughly within ±1 (unit variance). +mne.viz.plot_evoked_white(epochs["1"].average(), cov) # Next, we fit the dipoles for the evoked data. # We choose the timepoint which maximises global field power # We have seen in the evoked plot that this is around 3 ms after dipole onset. - data = [] t_peak = 0.036 # true for Elekta phantom for ii in event_id: @@ -150,12 +154,17 @@ evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.0) del epochs # save memory dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None) +# %% +# Whitened global field power (GFP): + +# most baseline activity should fall roughly within ±1 (unit variance). # %% # Let's visualize the explained variance. + # To do this, we need to make sure that the # data and the residuals are on the same scale -# (here the "time points" are the 32 dipole peak values that we fit): +# (here the "time points" are the 32 dipole peak values that we fit). fig, axes = plt.subplots(2, 1) evoked.plot(axes=axes) @@ -165,20 +174,19 @@ for line in ax.lines: line.set_color("#98df81") residual.plot(axes=axes) - +# %% # Here we visualise how well the dipole explains the evoked response (green line). # The red lines represent the residuals, the leftover noise after dipole fitting. # A good fit: green lines are strong and residuals are small and roughly flat. -# %% -# Finally, we compare the estimated to the true dipole locations -# To do this, we calculate the difference by ..... -# This is our ground truth. + +# Finally, we compare the estimated to the true dipole locations. actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 100.0 # nAm fig, (ax1, ax2, ax3) = plt.subplots( nrows=3, ncols=1, figsize=(6, 7), layout="constrained" ) + # 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)) @@ -201,18 +209,16 @@ ax3.bar(event_id, amps) ax3.set_xlabel("Dipole index") ax3.set_ylabel("Amplitude error (nAm)") - -# The dipole fits closely match the phantom ground truth. +# %% +# The dipole fits closely match the true phantom data. # We can achieve sub-centimeter accuracy with a mean position error of 2.6 mm. # This demonstrates that the fitting procedure is accurate. -# %% + # Finally, we can plot the positions and the orientations -# of the actual and the estimated dipoles +# of the estimated and true dipoles. actual_amp = np.ones(len(dip)) # fake amp, needed to create Dipole instance -actual_gof = np.ones( - len(dip) -) # fake goodness-of-fit (GOF), needed to create Dipole instance +actual_gof = np.ones(len(dip)) # fake goodness-of-fit (GOF) dip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof) fig = mne.viz.plot_alignment( @@ -227,7 +233,7 @@ subjects_dir=subjects_dir, ) -# Plot the position and the orientation of the actual dipole in black +# 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 ) @@ -238,9 +244,8 @@ ) mne.viz.set_3d_view(figure=fig, azimuth=70, elevation=80, distance=0.5) - -# Here the green arrows represent the estimated dipoles, the black arrows -# the true dipole location +# %% +# The dipoles overlap and point in the same direction. # %% # References # ---------- From 0be4f483e1f763ccbe21a13224f8ebb81094c746 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 14 Feb 2026 18:45:25 +1000 Subject: [PATCH 10/37] more formatting --- .../inverse/80_brainstorm_phantom_elekta.py | 50 ++++++++----------- 1 file changed, 21 insertions(+), 29 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index fbd0afa37c8..cc772f7bbe4 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -44,7 +44,7 @@ raw = read_raw_fif(raw_fname) # %% -# Let's first get an idea of the data structure we are working with: +# Let's first get an idea of the data structure we are working with. raw.info @@ -53,33 +53,34 @@ # 102 axial magnetometers, and 3 stimulus channels. # Next, let's look at the events in the phantom data for one stimulus channel. +# %% events = find_events(raw, "STI201") raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels - -# we have 32 artificial dipoles stored as events -# the remaining event IDs are not relevant for this tutorial (256, 768 and so on) # %% +# we have 32 artificial dipoles stored as events +# the remaining event IDs are not relevant for this tutorial (256, 768 ... ). + # Let's explore the phantom data by plotting power spectral density for each sensor. -# Here, we use only the first 30 seconds to save memory +# Here, we use only the first 30 seconds to save memory +# %% raw.compute_psd(tmax=30).plot( average=False, amplitude=False, picks="data", exclude="bads" ) - +# %% # We can see that the data has strong line frequency (60 Hz, 120 Hz ...) noise # and cHPI (continuous head position indicator) coil noise (peaks around 300 Hz). -# %% -# Next we plot the dipole events +# Next we plot the dipole events +# %% raw.plot(events=events, n_channels=10) - -# We can see that the simulated dipoles produce sinusoidal bursts at 20 Hz # %% -# Next, we epoch the the data based on the dipoles events (1:32). +# We can see that the simulated dipoles produce sinusoidal bursts at 20 Hz. +# Next, we epoch the the data based on the dipoles events (1:32). +# %% # We select 100 ms before and after the event trigger # and baseline correct the epochs from -100 ms to -0.05 ms before stimulus onset. - tmin, tmax = -0.1, 0.1 bmax = -0.05 # Avoid capture filter ringing into baseline event_id = list(range(1, 33)) @@ -91,15 +92,12 @@ # and plot the evoked signal epochs["1"].average().plot(time_unit="s") # %% - +# .. _plt_brainstorm_phantom_elekta_eeg_sphere_geometry: # We averaged over 640 simulated events for the first dipole. # The first peak in the data appears close to the trigger onset -# at around 3ms, with a peak repeating every 3ms. -# Thus the burst repetition rate is 3Hz. +# at around 3 ms. The burst envelope repeats at approximately 3 Hz. + -# %% -# .. _plt_brainstorm_phantom_elekta_eeg_sphere_geometry: -# # Finally, we source reconstruct the evoked simulated dipole. # To do this we use a :ref:`sphere head geometry model ` # and visualise the coordinate alignment and the sphere location. The phantom @@ -130,9 +128,7 @@ ) # %% # We can see that our head model aligns with the phantom head model. -# %% # Let's do some dipole fits. - # First we compute the noise covariance for the baseline window. # See covariance/whitening tutorial for details :ref:`tut-compute-covariance`. # %% @@ -156,16 +152,12 @@ dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None) # %% # Whitened global field power (GFP): - # most baseline activity should fall roughly within ±1 (unit variance). - -# %% # Let's visualize the explained variance. - # To do this, we need to make sure that the # data and the residuals are 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: @@ -178,8 +170,8 @@ # Here we visualise how well the dipole explains the evoked response (green line). # The red lines represent the residuals, the leftover noise after dipole fitting. # A good fit: green lines are strong and residuals are small and roughly flat. - # Finally, we compare the estimated to the true dipole locations. +# %% actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 100.0 # nAm @@ -210,9 +202,8 @@ ax3.set_xlabel("Dipole index") ax3.set_ylabel("Amplitude error (nAm)") # %% -# The dipole fits closely match the true phantom data. -# We can achieve sub-centimeter accuracy with a mean position error of 2.6 mm. -# This demonstrates that the fitting procedure is accurate. +# The dipole fits closely match the true phantom data, +# achieving sub-centimeter accuracy (mean position error 2.6 mm). # Finally, we can plot the positions and the orientations # of the estimated and true dipoles. @@ -246,6 +237,7 @@ mne.viz.set_3d_view(figure=fig, azimuth=70, elevation=80, distance=0.5) # %% # The dipoles overlap and point in the same direction. + # %% # References # ---------- From c5e8c3f8cb41d3f4f0fd7899ed419e47f39ab5f1 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 15 Feb 2026 14:21:40 +1000 Subject: [PATCH 11/37] final formatting --- .../inverse/80_brainstorm_phantom_elekta.py | 47 ++++++++----------- 1 file changed, 19 insertions(+), 28 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index cc772f7bbe4..3c9dd957561 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -51,34 +51,30 @@ # %% # The data consists of 204 MEG planor gradiometers, # 102 axial magnetometers, and 3 stimulus channels. - +# # Next, let's look at the events in the phantom data for one stimulus channel. -# %% events = find_events(raw, "STI201") raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels # %% -# we have 32 artificial dipoles stored as events -# the remaining event IDs are not relevant for this tutorial (256, 768 ... ). - +# There are 32 artificial dipoles stored as events. +# The remaining event IDs are not relevant for this tutorial (256, 768 ... ). +# # Let's explore the phantom data by plotting power spectral density for each sensor. - -# Here, we use only the first 30 seconds to save memory -# %% +# Here, we use only the first 30 seconds to save memory. raw.compute_psd(tmax=30).plot( average=False, amplitude=False, picks="data", exclude="bads" ) # %% -# We can see that the data has strong line frequency (60 Hz, 120 Hz ...) noise +# We can see that the data has strong line frequency noise (60 Hz, 120 Hz ...) # and cHPI (continuous head position indicator) coil noise (peaks around 300 Hz). - -# Next we plot the dipole events -# %% +# +# Next we plot the dipole events. raw.plot(events=events, n_channels=10) # %% # We can see that the simulated dipoles produce sinusoidal bursts at 20 Hz. - +# # Next, we epoch the the data based on the dipoles events (1:32). -# %% +# # We select 100 ms before and after the event trigger # and baseline correct the epochs from -100 ms to -0.05 ms before stimulus onset. tmin, tmax = -0.1, 0.1 @@ -96,8 +92,7 @@ # We averaged over 640 simulated events for the first dipole. # The first peak in the data appears close to the trigger onset # at around 3 ms. The burst envelope repeats at approximately 3 Hz. - - +# # Finally, we source reconstruct the evoked simulated dipole. # To do this we use a :ref:`sphere head geometry model ` # and visualise the coordinate alignment and the sphere location. The phantom @@ -109,7 +104,6 @@ # dipole locations differ. The phantom_otaniemi scan was aligned to the # phantom's head coordinate frame, so an identity ``trans`` is appropriate # here. - 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) @@ -128,10 +122,11 @@ ) # %% # We can see that our head model aligns with the phantom head model. +# # Let's do some dipole fits. # First we compute the noise covariance for the baseline window. -# See covariance/whitening tutorial for details :ref:`tut-compute-covariance`. -# %% +# Check the covariance/whitening tutorial for details :ref:`tut-compute-covariance`. +# # The covariance captures the sensor noise structure. cov = mne.compute_covariance(epochs, tmax=bmax) @@ -157,7 +152,6 @@ # To do this, we need to make sure that the # data and the residuals are 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: @@ -169,9 +163,10 @@ # %% # Here we visualise how well the dipole explains the evoked response (green line). # The red lines represent the residuals, the leftover noise after dipole fitting. -# A good fit: green lines are strong and residuals are small and roughly flat. +# What is a good fit? The green lines are strong and residuals are small and +# roughly flat. +# # Finally, we compare the estimated to the true dipole locations. -# %% actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 100.0 # nAm @@ -204,10 +199,9 @@ # %% # The dipole fits closely match the true phantom data, # achieving sub-centimeter accuracy (mean position error 2.6 mm). - +# # Finally, we can plot the positions and the orientations # of the estimated and true dipoles. - actual_amp = np.ones(len(dip)) # fake amp, needed to create Dipole instance actual_gof = np.ones(len(dip)) # fake goodness-of-fit (GOF) dip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof) @@ -223,20 +217,17 @@ show_axes=True, subjects_dir=subjects_dir, ) - # 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 in green fig = mne.viz.plot_dipole_locations( dipoles=dip, 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) # %% -# The dipoles overlap and point in the same direction. +# We can see that the dipoles overlap and point in the same direction. # %% # References From 9ccff3f19c89ae2cf550c3fe673821281f102f33 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 15 Feb 2026 14:25:45 +1000 Subject: [PATCH 12/37] add blank lines around directives --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 3c9dd957561..091096a4740 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -232,4 +232,6 @@ # %% # References # ---------- +# # .. footbibliography:: +# From 9fc7c4aba4022c8706f29089e06f8c6afc9963dd Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 15 Feb 2026 14:57:57 +1000 Subject: [PATCH 13/37] add blank line at file end --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 091096a4740..98db3babba2 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -68,10 +68,10 @@ # We can see that the data has strong line frequency noise (60 Hz, 120 Hz ...) # and cHPI (continuous head position indicator) coil noise (peaks around 300 Hz). # -# Next we plot the dipole events. +# Here we plot the dipole events. raw.plot(events=events, n_channels=10) # %% -# We can see that the simulated dipoles produce sinusoidal bursts at 20 Hz. +# The simulated dipoles produce sinusoidal bursts at 20 Hz. # # Next, we epoch the the data based on the dipoles events (1:32). # @@ -148,6 +148,7 @@ # %% # Whitened global field power (GFP): # most baseline activity should fall roughly within ±1 (unit variance). +# # Let's visualize the explained variance. # To do this, we need to make sure that the # data and the residuals are on the same scale @@ -163,6 +164,7 @@ # %% # Here we visualise how well the dipole explains the evoked response (green line). # The red lines represent the residuals, the leftover noise after dipole fitting. +# # What is a good fit? The green lines are strong and residuals are small and # roughly flat. # From e8299e2f1e9b7ca2bb8486a88854e6db2a53a26a Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 16 Feb 2026 19:46:16 +1000 Subject: [PATCH 14/37] fix sphinx error log --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 98db3babba2..5ed0162d5e9 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -237,3 +237,4 @@ # # .. footbibliography:: # +# From 35949f4b7fe43a25a3be75e5f26440dd6e2e71ee Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 21 Feb 2026 15:22:50 +1000 Subject: [PATCH 15/37] fix sphinx formatting error --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 5ed0162d5e9..09686f19493 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -89,6 +89,7 @@ epochs["1"].average().plot(time_unit="s") # %% # .. _plt_brainstorm_phantom_elekta_eeg_sphere_geometry: +# # We averaged over 640 simulated events for the first dipole. # The first peak in the data appears close to the trigger onset # at around 3 ms. The burst envelope repeats at approximately 3 Hz. @@ -132,7 +133,7 @@ # The plot shows the evoked signal divided by the estimated noise standard deviation. mne.viz.plot_evoked_white(epochs["1"].average(), cov) - +# %% # Next, we fit the dipoles for the evoked data. # We choose the timepoint which maximises global field power # We have seen in the evoked plot that this is around 3 ms after dipole onset. @@ -234,7 +235,4 @@ # %% # References # ---------- -# # .. footbibliography:: -# -# From 139d273c449649e204a49684a97d6a34eed7bc16 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 29 Mar 2026 14:51:55 +1000 Subject: [PATCH 16/37] address dan's review --- .../inverse/80_brainstorm_phantom_elekta.py | 173 ++++++++---------- 1 file changed, 72 insertions(+), 101 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 09686f19493..17100e023a7 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -7,11 +7,16 @@ This tutorial provides a step-by-step guide to importing and processing Elekta-Neuromag current phantom recordings. -The aim of this tutorial is to show the user how to use phantom recordings to + +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 learn how to use phantom recordings to evaluate source localisation methods by comparing estimated vs real dipole positions. For comparison, see :footcite:`TadelEtAl2011` and -`the original Brainstorm tutorial with an explanation of phantom recordings +`the original Brainstorm tutorial `__. """ # sphinx_gallery_thumbnail_number = 9 @@ -25,6 +30,7 @@ import matplotlib.pyplot as plt import numpy as np +from scipy.signal import find_peaks import mne from mne import find_events, fit_dipole @@ -37,131 +43,92 @@ # %% # 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, amplitudes can be seen in raw data) -# data are accessed to construct instances of :class:`mne.io.Raw`. +# The dataset has recordings at 3 different current amplitudes (20, 200, and 2000 nAm), +# here we'll load the medium-amplitude one. 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) - -# %% -# Let's first get an idea of the data structure we are working with. - -raw.info - -# %% -# The data consists of 204 MEG planor gradiometers, -# 102 axial magnetometers, and 3 stimulus channels. +raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels # # Next, let's look at the events in the phantom data for one stimulus channel. events = find_events(raw, "STI201") -raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels # %% # There are 32 artificial dipoles stored as events. # The remaining event IDs are not relevant for this tutorial (256, 768 ... ). -# -# Let's explore the phantom data by plotting power spectral density for each sensor. -# Here, we use only the first 30 seconds to save memory. -raw.compute_psd(tmax=30).plot( - average=False, amplitude=False, picks="data", exclude="bads" -) -# %% -# We can see that the data has strong line frequency noise (60 Hz, 120 Hz ...) -# and cHPI (continuous head position indicator) coil noise (peaks around 300 Hz). -# -# Here we plot the dipole events. -raw.plot(events=events, n_channels=10) + # %% -# The simulated dipoles produce sinusoidal bursts at 20 Hz. +# In this data the phantom was set to produce 20 Hz sinusoidal bursts of current. # -# Next, we epoch the the data based on the dipoles events (1:32). -# -# We select 100 ms before and after the event trigger +# Next, we epoch the data based on the dipole events (1:32). # and baseline correct the epochs from -100 ms to -0.05 ms before stimulus onset. -tmin, tmax = -0.1, 0.1 -bmax = -0.05 # Avoid capture filter ringing into baseline +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 + raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=True ) - -# Here we average the epochs for the first simulated dipole -# and plot the evoked signal -epochs["1"].average().plot(time_unit="s") # %% -# .. _plt_brainstorm_phantom_elekta_eeg_sphere_geometry: -# +# Dipole fitting on phantom data +# ------------------------------- +# Next, we use the epoched data to fit dipoles, +# then compare them to the known dipole locations built into the phantom. + +# We need to determine the peak amplitude timepoint we want +# to estimate the dipole at + +# We use the first dipole event to plot the evoked signal +# Skip first and last epoch as they are corrupted +epochs["1"][1:-1].average().plot(time_unit="s") +evoked_tmp = epochs["1"][1:-1].average() # We averaged over 640 simulated events for the first dipole. -# The first peak in the data appears close to the trigger onset -# at around 3 ms. The burst envelope repeats at approximately 3 Hz. -# -# Finally, we source reconstruct the evoked simulated dipole. -# To do this we use a :ref:`sphere head geometry model ` -# and visualise 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. -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, -) -# %% -# We can see that our head model aligns with the phantom head model. -# -# Let's do some dipole fits. -# First we compute the noise covariance for the baseline window. -# Check the covariance/whitening tutorial for details :ref:`tut-compute-covariance`. -# -# The covariance captures the sensor noise structure. -cov = mne.compute_covariance(epochs, tmax=bmax) +# You can see that the burst envelope repeats at approximately 3 Hz. -# The plot shows the evoked signal divided by the estimated noise standard deviation. -mne.viz.plot_evoked_white(epochs["1"].average(), cov) -# %% -# Next, we fit the dipoles for the evoked data. -# We choose the timepoint which maximises global field power -# We have seen in the evoked plot that this is around 3 ms after dipole onset. +# Next we compute Global Field Power (GFP) +gfp = np.std(evoked_tmp.data, axis=0) + +# Determine first peak (after time 0) +times = evoked_tmp.times +post_stim_mask = times > 0 + +# Get index of first peak +peaks, _ = find_peaks(gfp[post_stim_mask]) + +# First peak index in original time vector +first_peak_idx = np.where(post_stim_mask)[0][peaks[0]] +t_peak = times[first_peak_idx] + +print(f"Detected first peak at {t_peak * 1000:.1f} ms") + +# We store the evoked data for each dipole and +# create an evoked array containing all dipoles 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) + +# Next, we need to compute the noise covariance for dipole fitting. +# Check the covariance/whitening tutorial for details :ref:`tut-compute-covariance`. +# +# The covariance captures the sensor noise structure. +cov = mne.compute_covariance(epochs, tmax=bmax) + del epochs # save memory + +# We use a :ref:`sphere head geometry model ` +# which aligns well with our phantom head model. +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) dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None) # %% -# Whitened global field power (GFP): -# most baseline activity should fall roughly within ±1 (unit variance). -# # Let's visualize the explained variance. -# To do this, we need to make sure that the -# data and the residuals are 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) +# The dipole object contains the goodness of fit (GOF) for each dipole. +plt.plot(dip.gof) +plt.xlabel("Dipole") +plt.ylabel("GOF") +plt.show() # %% # Here we visualise how well the dipole explains the evoked response (green line). # The red lines represent the residuals, the leftover noise after dipole fitting. @@ -209,6 +176,9 @@ actual_gof = np.ones(len(dip)) # fake goodness-of-fit (GOF) dip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof) +subject = "phantom_otaniemi" +trans = mne.transforms.Transform("head", "mri", np.eye(4)) + fig = mne.viz.plot_alignment( evoked.info, trans, @@ -230,7 +200,8 @@ ) mne.viz.set_3d_view(figure=fig, azimuth=70, elevation=80, distance=0.5) # %% -# We can see that the dipoles overlap and point in the same direction. +# We can see that the dipoles overlap, have approximately the same magnitude +# and point in the same direction. # %% # References From bd52896c2eb04443148968be4d69a1eb398c32cf Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 29 Mar 2026 15:04:50 +1000 Subject: [PATCH 17/37] fix t_peak --- .../inverse/80_brainstorm_phantom_elekta.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 17100e023a7..5f80337cd9d 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -65,7 +65,7 @@ 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=True + raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False ) # %% # Dipole fitting on phantom data @@ -97,6 +97,7 @@ first_peak_idx = np.where(post_stim_mask)[0][peaks[0]] t_peak = times[first_peak_idx] +t_peak = 0.036 print(f"Detected first peak at {t_peak * 1000:.1f} ms") # We store the evoked data for each dipole and @@ -121,21 +122,18 @@ 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) + +# Let's finally fit some dipoles dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None) # %% # Let's visualize the explained variance. -# The dipole object contains the goodness of fit (GOF) for each dipole. -plt.plot(dip.gof) +# The dipole object stores the goodness of fit (GOF) for each dipole. +plt.plot(dip.gof) # something wrong, GOF of 6? plt.xlabel("Dipole") plt.ylabel("GOF") plt.show() # %% -# Here we visualise how well the dipole explains the evoked response (green line). -# The red lines represent the residuals, the leftover noise after dipole fitting. -# -# What is a good fit? The green lines are strong and residuals are small and -# roughly flat. -# + # Finally, we compare the estimated to the true dipole locations. actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 100.0 # nAm From d39b6e5a644ab2f89a30194a7c1ba38d67950235 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 29 Mar 2026 19:40:26 +1000 Subject: [PATCH 18/37] fixed peak amplitude estimation --- .../inverse/80_brainstorm_phantom_elekta.py | 75 +++++++++---------- 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 5f80337cd9d..aef63890b0e 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -44,22 +44,18 @@ # The data were collected with an Elekta Neuromag VectorView system # at 1000 Hz and low-pass filtered at 330 Hz. # The dataset has recordings at 3 different current amplitudes (20, 200, and 2000 nAm), -# here we'll load the medium-amplitude one. +# here we will load the medium-amplitude current. 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) raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels # -# Next, let's look at the events in the phantom data for one stimulus channel. +# Here we store the events for one stimulus channel. events = find_events(raw, "STI201") -# %% -# There are 32 artificial dipoles stored as events. -# The remaining event IDs are not relevant for this tutorial (256, 768 ... ). +# The first 32 events are the dipole phantoms recorded. +# The remaining event IDs are not relevant for this tutorial. -# %% -# In this data the phantom was set to produce 20 Hz sinusoidal bursts of current. -# -# Next, we epoch the data based on the dipole events (1:32). +# Next, we epoch the data for each dipole event # and baseline correct the epochs from -100 ms to -0.05 ms before stimulus onset. bmax = -0.05 tmin, tmax = -0.1, 0.8 @@ -70,68 +66,71 @@ # %% # Dipole fitting on phantom data # ------------------------------- -# Next, we use the epoched data to fit dipoles, +# We use the epoched data to fit dipoles, # then compare them to the known dipole locations built into the phantom. -# We need to determine the peak amplitude timepoint we want -# to estimate the dipole at - # We use the first dipole event to plot the evoked signal -# Skip first and last epoch as they are corrupted +# and skip first and last epoch as they are corrupted. epochs["1"][1:-1].average().plot(time_unit="s") evoked_tmp = epochs["1"][1:-1].average() -# We averaged over 640 simulated events for the first dipole. +# We averaged over 18 simulated events for the first dipole. +# In this data the phantom was set to produce 20 Hz sinusoidal bursts of current. # You can see that the burst envelope repeats at approximately 3 Hz. -# Next we compute Global Field Power (GFP) +# First, we need to determine the timepoint of the peak amplitude. +# To find this timepoint, we compute Global Field Power (GFP) gfp = np.std(evoked_tmp.data, axis=0) - -# Determine first peak (after time 0) times = evoked_tmp.times -post_stim_mask = times > 0 -# Get index of first peak -peaks, _ = find_peaks(gfp[post_stim_mask]) +# We wamt the peak for the first bursting event +time_mask = (times > 0) & (times <= 0.1) + +# Find peaks in that window +peaks, _ = find_peaks(gfp[time_mask]) -# First peak index in original time vector -first_peak_idx = np.where(post_stim_mask)[0][peaks[0]] -t_peak = times[first_peak_idx] +# Convert to original indices +peak_indices = np.where(time_mask)[0][peaks] -t_peak = 0.036 -print(f"Detected first peak at {t_peak * 1000:.1f} ms") +# Select the strongest peak (max GFP) +strongest_peak_idx = peak_indices[np.argmax(gfp[peak_indices])] +t_peak = times[strongest_peak_idx] -# We store the evoked data for each dipole and -# create an evoked array containing all dipoles +print(f"Strongest peak at {t_peak * 1000:.1f} ms") + +# Here we store the evoked data for each dipole at the peak amplitude. data = [] for ii in event_id: evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak) data.append(evoked.data[:, 0]) +# Finally, we store the evokeds from all dipoles as an evoked, +# treating every dipole as a timepoint. evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.0) -# Next, we need to compute the noise covariance for dipole fitting. -# Check the covariance/whitening tutorial for details :ref:`tut-compute-covariance`. -# -# The covariance captures the sensor noise structure. +# Next, we need to compute the noise covariance to capture the sensor noise structure. +# We use the baseline window to estimate covariance. +# You can explore the covariance tutorial for details :ref:`tut-compute-covariance`. cov = mne.compute_covariance(epochs, tmax=bmax) -del epochs # save memory +del epochs # delete to save memory # We use a :ref:`sphere head geometry model ` -# which aligns well with our phantom head model. +# to fit our phantom head model. 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) -# Let's finally fit some dipoles -dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None) +# Let's finally fit all 32 phantom dipoles. +dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1) # %% # Let's visualize the explained variance. # The dipole object stores the goodness of fit (GOF) for each dipole. plt.plot(dip.gof) # something wrong, GOF of 6? plt.xlabel("Dipole") -plt.ylabel("GOF") +plt.ylabel("GOF %") plt.show() +# We can see that GOF varies between 50 and more than 95 % +# variance explained. # %% # Finally, we compare the estimated to the true dipole locations. @@ -166,7 +165,7 @@ ax3.set_ylabel("Amplitude error (nAm)") # %% # The dipole fits closely match the true phantom data, -# achieving sub-centimeter accuracy (mean position error 2.6 mm). +# achieving sub-centimeter accuracy (mean position error 2.7mm). # # Finally, we can plot the positions and the orientations # of the estimated and true dipoles. From bd754f8515f3a43d227294cfe019efd3ecf57a9f Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 29 Mar 2026 19:48:38 +1000 Subject: [PATCH 19/37] prettified GOF plot --- .../inverse/80_brainstorm_phantom_elekta.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index aef63890b0e..6348bc8c65f 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -125,17 +125,23 @@ # %% # Let's visualize the explained variance. # The dipole object stores the goodness of fit (GOF) for each dipole. -plt.plot(dip.gof) # something wrong, GOF of 6? -plt.xlabel("Dipole") -plt.ylabel("GOF %") +gof = dip.gof # array of GOF values + +# Define colorblind-friendly colors +colors = ["#E69F00" if val < 60 else "#0072B2" for val in gof] +# orange = low, blue = high (from Okabe-Ito palette) + +plt.bar(event_id, gof, color=colors) +plt.xlabel("Phantom dipole estimation") +plt.ylabel("Goodness of fit (%)") plt.show() # We can see that GOF varies between 50 and more than 95 % -# variance explained. +# variance explained for dipoles. # %% # Finally, we compare the estimated to the true dipole locations. actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() -actual_amp = 100.0 # nAm +actual_amp = 200.0 # nAm fig, (ax1, ax2, ax3) = plt.subplots( nrows=3, ncols=1, figsize=(6, 7), layout="constrained" From fa6976a7b83ead3aa00f9f1f9c5a5e77dc298803 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 29 Mar 2026 20:10:53 +1000 Subject: [PATCH 20/37] single evokeds instead of evoked with dipoles as timepoints --- .../inverse/80_brainstorm_phantom_elekta.py | 27 +++++++++++-------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 6348bc8c65f..21fcab312e2 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -98,14 +98,11 @@ print(f"Strongest peak at {t_peak * 1000:.1f} ms") # Here we store the evoked data for each dipole at the peak amplitude. -data = [] +evokeds = [] for ii in event_id: evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak) - data.append(evoked.data[:, 0]) - -# Finally, we store the evokeds from all dipoles as an evoked, -# treating every dipole as a timepoint. -evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.0) + evoked = mne.EvokedArray(np.array(evoked.data), evoked.info, tmin=0.0) + evokeds.append(evoked) # Next, we need to compute the noise covariance to capture the sensor noise structure. # We use the baseline window to estimate covariance. @@ -120,12 +117,16 @@ fetch_phantom("otaniemi", subjects_dir=subjects_dir) sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08) +dip_all, residuals_all = [], [] # Let's finally fit all 32 phantom dipoles. -dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1) +for evoked in evokeds: + dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1) + dip_all.append(dip) + residuals_all.append(residual) # %% # Let's visualize the explained variance. # The dipole object stores the goodness of fit (GOF) for each dipole. -gof = dip.gof # array of GOF values +gof = [dip.gof[0] for dip in dip_all] # array of GOF values # Define colorblind-friendly colors colors = ["#E69F00" if val < 60 else "#0072B2" for val in gof] @@ -142,6 +143,10 @@ # Finally, we compare the estimated to the true dipole locations. actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 200.0 # nAm +# dipole estimations +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" @@ -149,7 +154,7 @@ # 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)) +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") @@ -157,14 +162,14 @@ # 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)))) +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 (°)") # Finally we compare amplitudes by subtracting estimated from true amplitude. -amps = actual_amp - dip.amplitude / 1e-9 +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") From f2a079b6076cb76e0a8447683078bf486184611d Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 29 Mar 2026 20:11:25 +1000 Subject: [PATCH 21/37] add TODO --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 21fcab312e2..f71311d8afc 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -175,6 +175,7 @@ ax3.set_xlabel("Dipole index") ax3.set_ylabel("Amplitude error (nAm)") # %% +# TODO: fix dipole vis # The dipole fits closely match the true phantom data, # achieving sub-centimeter accuracy (mean position error 2.7mm). # From 4e2cdc6ebebb1f9a70ad40fe97717d3ab85b7dd2 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 29 Mar 2026 20:12:02 +1000 Subject: [PATCH 22/37] add formatting TODO --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index f71311d8afc..9c138724ace 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -175,7 +175,7 @@ ax3.set_xlabel("Dipole index") ax3.set_ylabel("Amplitude error (nAm)") # %% -# TODO: fix dipole vis +# TODO: fix dipole vis and format the docs # The dipole fits closely match the true phantom data, # achieving sub-centimeter accuracy (mean position error 2.7mm). # From 5046a47f26557c8662ec59f9584d5a5ed87c8a52 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 6 Apr 2026 13:32:18 +1000 Subject: [PATCH 23/37] removed plotting label --- doc/_includes/forward.rst | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/doc/_includes/forward.rst b/doc/_includes/forward.rst index 40706fb292c..05cd272a6b1 100644 --- a/doc/_includes/forward.rst +++ b/doc/_includes/forward.rst @@ -690,8 +690,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 + :ref:`Brainstorm Elekta phantom dataset tutorial`, and :ref:`tut-source-alignment-without-mri`. When the sphere model is employed, the computation of the EEG solution can be From 93be3a67a995d273467c2702ae395e1b90457710 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 6 Apr 2026 14:33:01 +1000 Subject: [PATCH 24/37] fix 3D plotting (very slow, discuss with Dan) --- .../inverse/80_brainstorm_phantom_elekta.py | 66 +++++++++++-------- 1 file changed, 37 insertions(+), 29 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 9c138724ace..976944620a6 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -22,6 +22,7 @@ # sphinx_gallery_thumbnail_number = 9 # Authors: Eric Larson +# Carina Forster # # License: BSD-3-Clause # Copyright the MNE-Python contributors. @@ -39,8 +40,9 @@ from mne.io import read_raw_fif print(__doc__) - # %% +# Load and prepare the data +# ------------------------------- # The data were collected with an Elekta Neuromag VectorView system # at 1000 Hz and low-pass filtered at 330 Hz. # The dataset has recordings at 3 different current amplitudes (20, 200, and 2000 nAm), @@ -49,13 +51,11 @@ raw_fname = data_path / "kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif" raw = read_raw_fif(raw_fname) raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels -# -# Here we store the events for one stimulus channel. events = find_events(raw, "STI201") # The first 32 events are the dipole phantoms recorded. # The remaining event IDs are not relevant for this tutorial. - -# Next, we epoch the data for each dipole event +# %% +# Here we epoch the data around dipole events # and baseline correct the epochs from -100 ms to -0.05 ms before stimulus onset. bmax = -0.05 tmin, tmax = -0.1, 0.8 @@ -66,59 +66,65 @@ # %% # Dipole fitting on phantom data # ------------------------------- -# We use the epoched data to fit dipoles, -# then compare them to the known dipole locations built into the phantom. +# Next, we fit dipoles based on the epoched phantom recordings and +# compare the estimated to the known dipole locations built into the phantom. # We use the first dipole event to plot the evoked signal -# and skip first and last epoch as they are corrupted. +# and skip the first and the last epoch as they are corrupted. + epochs["1"][1:-1].average().plot(time_unit="s") -evoked_tmp = epochs["1"][1:-1].average() -# We averaged over 18 simulated events for the first dipole. + +# We average over 18 simulated events for the first dipole. # In this data the phantom was set to produce 20 Hz sinusoidal bursts of current. # You can see that the burst envelope repeats at approximately 3 Hz. - +# %% # First, we need to determine the timepoint of the peak amplitude. -# To find this timepoint, we compute Global Field Power (GFP) +# To find this timepoint, we compute Global Field Power (GFP), +# which is the standard deviation of all EEG electrode voltages at a +# specific timepoint. +evoked_tmp = epochs["1"][1:-1].average() gfp = np.std(evoked_tmp.data, axis=0) times = evoked_tmp.times -# We wamt the peak for the first bursting event +# We restrict the time window to find the peak +# for the first bursting event. time_mask = (times > 0) & (times <= 0.1) - -# Find peaks in that window peaks, _ = find_peaks(gfp[time_mask]) - # Convert to original indices peak_indices = np.where(time_mask)[0][peaks] - # Select the strongest peak (max GFP) 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") - +# %% # Here we store the evoked data for each dipole at the peak amplitude. evokeds = [] for ii in event_id: evoked = epochs[str(ii)][1:-1].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 to capture the sensor noise structure. # We use the baseline window to estimate covariance. # You can explore the covariance tutorial for details :ref:`tut-compute-covariance`. + cov = mne.compute_covariance(epochs, tmax=bmax) del epochs # delete to save memory - +# %% # We use a :ref:`sphere head geometry model ` # to fit our phantom head model. 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) +# %% +# We finally fit all 32 phantom dipoles and store them +# as well as the residuals in a list. + dip_all, residuals_all = [], [] -# Let's finally fit all 32 phantom dipoles. + for evoked in evokeds: dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1) dip_all.append(dip) @@ -136,14 +142,16 @@ plt.xlabel("Phantom dipole estimation") plt.ylabel("Goodness of fit (%)") plt.show() +# %% # We can see that GOF varies between 50 and more than 95 % # variance explained for dipoles. # %% - +# Estimated vs true dipole locations +# ------------------------------- # Finally, we compare the estimated to the true dipole locations. actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 200.0 # nAm -# dipole estimations +# 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] @@ -168,14 +176,13 @@ ax2.set_xlabel("Dipole index") ax2.set_ylabel("Angle error (°)") -# Finally we compare amplitudes by subtracting estimated from true amplitude. +# 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)") # %% -# TODO: fix dipole vis and format the docs # The dipole fits closely match the true phantom data, # achieving sub-centimeter accuracy (mean position error 2.7mm). # @@ -203,10 +210,11 @@ 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 in green -fig = mne.viz.plot_dipole_locations( - dipoles=dip, mode="arrow", subject=subject, color=(0.2, 1.0, 0.5), fig=fig -) +for dip in dip_all: + # 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 + ) mne.viz.set_3d_view(figure=fig, azimuth=70, elevation=80, distance=0.5) # %% # We can see that the dipoles overlap, have approximately the same magnitude From 722692e14abc061f9741492dcd0fa8d5e4dd2aa3 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 6 Apr 2026 14:46:39 +1000 Subject: [PATCH 25/37] second dipole object instead of heavy loop --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 976944620a6..99bf8cba6bb 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -145,7 +145,7 @@ # %% # We can see that GOF varies between 50 and more than 95 % # variance explained for dipoles. -# %% + # Estimated vs true dipole locations # ------------------------------- # Finally, we compare the estimated to the true dipole locations. @@ -190,7 +190,9 @@ # of the estimated and true dipoles. actual_amp = np.ones(len(dip)) # fake amp, 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) subject = "phantom_otaniemi" trans = mne.transforms.Transform("head", "mri", np.eye(4)) @@ -210,11 +212,10 @@ fig = mne.viz.plot_dipole_locations( dipoles=dip_true, mode="arrow", subject=subject, color=(0.0, 0.0, 0.0), fig=fig ) -for dip in dip_all: - # 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 - ) +# Plot the position and the orientation of the estimated dipole in green +fig = mne.viz.plot_dipole_locations( + 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) # %% # We can see that the dipoles overlap, have approximately the same magnitude From 01aa3f0eb5c3ae6424da7e58cf672dd174a01bd7 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 6 Apr 2026 15:51:56 +1000 Subject: [PATCH 26/37] final formatting (ask Dan regarding ouputs) --- .../inverse/80_brainstorm_phantom_elekta.py | 99 ++++++++++--------- 1 file changed, 52 insertions(+), 47 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 99bf8cba6bb..6816ae898f0 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -12,8 +12,8 @@ that generates known magnetic signals, allowing validation and benchmarking of MEG system accuracy and analysis methods. -The aim of this tutorial is to learn how to use phantom recordings to -evaluate source localisation methods by comparing estimated vs real dipole positions. +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 @@ -28,7 +28,6 @@ # Copyright the MNE-Python contributors. # %% - import matplotlib.pyplot as plt import numpy as np from scipy.signal import find_peaks @@ -39,60 +38,61 @@ from mne.datasets.brainstorm import bst_phantom_elekta from mne.io import read_raw_fif -print(__doc__) # %% # Load and prepare the data # ------------------------------- + # The data were collected with an Elekta Neuromag VectorView system # at 1000 Hz and low-pass filtered at 330 Hz. -# The dataset has recordings at 3 different current amplitudes (20, 200, and 2000 nAm), -# here we will load the medium-amplitude current. +# +# The dataset contains recordings at three current amplitudes (20, 200, and 2000 nAm). +# Here we load the medium-amplitude condition. 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) -raw.info["bads"] = ["MEG1933", "MEG2421"] # known bad channels + +# Mark known bad channels +raw.info["bads"] = ["MEG1933", "MEG2421"] + events = find_events(raw, "STI201") -# The first 32 events are the dipole phantoms recorded. -# The remaining event IDs are not relevant for this tutorial. + +# The first 32 events correspond to dipole activations. + # %% -# Here we epoch the data around dipole events -# and baseline correct the epochs from -100 ms to -0.05 ms before stimulus onset. +# Epoch the data and plot evokeds +# --------------- +# +# We epoch the data around dipole events and apply baseline correction. + 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 ) -# %% -# Dipole fitting on phantom data -# ------------------------------- -# Next, we fit dipoles based on the epoched phantom recordings and -# compare the estimated to the known dipole locations built into the phantom. - -# We use the first dipole event to plot the evoked signal -# and skip the first and the last epoch as they are corrupted. - +# Plot evoked response for the first clean dipole epochs["1"][1:-1].average().plot(time_unit="s") - -# We average over 18 simulated events for the first dipole. +# %% # In this data the phantom was set to produce 20 Hz sinusoidal bursts of current. # You can see that the burst envelope repeats at approximately 3 Hz. -# %% -# First, we need to determine the timepoint of the peak amplitude. -# To find this timepoint, we compute Global Field Power (GFP), -# which is the standard deviation of all EEG electrode voltages at a -# specific timepoint. +# +# 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. evoked_tmp = epochs["1"][1:-1].average() gfp = np.std(evoked_tmp.data, axis=0) times = evoked_tmp.times -# We restrict the time window to find the peak -# for the first bursting event. +# Restrict to first burst window time_mask = (times > 0) & (times <= 0.1) + peaks, _ = find_peaks(gfp[time_mask]) -# Convert to original indices peak_indices = np.where(time_mask)[0][peaks] -# Select the strongest peak (max GFP) + +# Select the strongest peak strongest_peak_idx = peak_indices[np.argmax(gfp[peak_indices])] t_peak = times[strongest_peak_idx] @@ -107,7 +107,8 @@ # %% # Next, we need to compute the noise covariance to capture the sensor noise structure. # We use the baseline window to estimate covariance. -# You can explore the covariance tutorial for details :ref:`tut-compute-covariance`. +# +# You can explore the covariance tutorial for details: :ref:`tut-compute-covariance`. cov = mne.compute_covariance(epochs, tmax=bmax) @@ -120,9 +121,10 @@ sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08) # %% -# We finally fit all 32 phantom dipoles and store them -# as well as the residuals in a list. +# Fit dipoles +# ----------- +# We fit dipoles for each phantom and store them in a list. dip_all, residuals_all = [], [] for evoked in evokeds: @@ -130,27 +132,27 @@ dip_all.append(dip) residuals_all.append(residual) # %% -# Let's visualize the explained variance. -# The dipole object stores the goodness of fit (GOF) for each dipole. -gof = [dip.gof[0] for dip in dip_all] # array of GOF values +# Evaluate goodness of fit +# ----------------------- -# Define colorblind-friendly colors +# The dipole object stores the goodness of fit (GOF) for each dipole. +gof = [dip.gof[0] for dip in dip_all] colors = ["#E69F00" if val < 60 else "#0072B2" for val in gof] -# orange = low, blue = high (from Okabe-Ito palette) - plt.bar(event_id, gof, color=colors) plt.xlabel("Phantom dipole estimation") plt.ylabel("Goodness of fit (%)") plt.show() + # %% -# We can see that GOF varies between 50 and more than 95 % -# variance explained for dipoles. +# +# We can see that GOF varies between 50 % and up to 95 %. +# +# Compare estimated and true dipoles +# -------------------------------- -# Estimated vs true dipole locations -# ------------------------------- -# Finally, we compare the estimated to the true dipole locations. actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 200.0 # nAm + # estimated dipoles dip_pos = [dip.pos[0] for dip in dip_all] dip_ori = [dip.ori[0] for dip in dip_all] @@ -186,8 +188,9 @@ # The dipole fits closely match the true phantom data, # achieving sub-centimeter accuracy (mean position error 2.7mm). # -# Finally, we can plot the positions and the orientations -# of the estimated and true dipoles. +# Visualise estimated and true dipoles +# ----------------------- + actual_amp = np.ones(len(dip)) # fake amp, needed to create Dipole instance actual_gof = np.ones(len(dip)) # fake goodness-of-fit (GOF) # setup dipole objects for true and estimated dipoles @@ -208,10 +211,12 @@ show_axes=True, subjects_dir=subjects_dir, ) + # 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 in green fig = mne.viz.plot_dipole_locations( dipoles=dip_estimated, mode="arrow", subject=subject, color=(0.2, 1.0, 0.5), fig=fig From e22e2059567369de884f81539489942709a8ad51 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 6 Apr 2026 15:57:06 +1000 Subject: [PATCH 27/37] quick formatting --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 6816ae898f0..5fa186838e7 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -71,6 +71,7 @@ epochs = mne.Epochs( raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False ) + # Plot evoked response for the first clean dipole epochs["1"][1:-1].average().plot(time_unit="s") # %% @@ -106,6 +107,7 @@ evokeds.append(evoked) # %% # Next, we need to compute the noise covariance to capture the sensor noise structure. +# # We use the baseline window to estimate covariance. # # You can explore the covariance tutorial for details: :ref:`tut-compute-covariance`. @@ -144,7 +146,6 @@ plt.show() # %% -# # We can see that GOF varies between 50 % and up to 95 %. # # Compare estimated and true dipoles From b175b501f100e6e49c1d0b91df8ccc5823a4009a Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 6 Apr 2026 16:20:02 +1000 Subject: [PATCH 28/37] removed sphere model reference --- doc/_includes/forward.rst | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/doc/_includes/forward.rst b/doc/_includes/forward.rst index 05cd272a6b1..2ed1ecacaab 100644 --- a/doc/_includes/forward.rst +++ b/doc/_includes/forward.rst @@ -689,8 +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 From 2b1b4a3b9683ba04db53c9d4703697606a2cf9e6 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 6 Apr 2026 16:21:02 +1000 Subject: [PATCH 29/37] sphinx error --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 5fa186838e7..f92f0127f9c 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -40,7 +40,7 @@ # %% # Load and prepare the data -# ------------------------------- +# -------------------------- # The data were collected with an Elekta Neuromag VectorView system # at 1000 Hz and low-pass filtered at 330 Hz. @@ -60,7 +60,7 @@ # %% # Epoch the data and plot evokeds -# --------------- +# -------------------------------- # # We epoch the data around dipole events and apply baseline correction. @@ -79,7 +79,7 @@ # You can see that the burst envelope repeats at approximately 3 Hz. # # 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. @@ -124,7 +124,7 @@ # %% # Fit dipoles -# ----------- +# ------------ # We fit dipoles for each phantom and store them in a list. dip_all, residuals_all = [], [] @@ -149,7 +149,7 @@ # We can see that GOF varies between 50 % and up to 95 %. # # Compare estimated and true dipoles -# -------------------------------- +# ----------------------------------- actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 200.0 # nAm @@ -190,7 +190,7 @@ # achieving sub-centimeter accuracy (mean position error 2.7mm). # # Visualise estimated and true dipoles -# ----------------------- +# ------------------------------------- actual_amp = np.ones(len(dip)) # fake amp, needed to create Dipole instance actual_gof = np.ones(len(dip)) # fake goodness-of-fit (GOF) From e27ac56d086ba3b600443f97e533e02d3401577f Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Mon, 6 Apr 2026 16:57:49 +1000 Subject: [PATCH 30/37] fix doc build sphinx error --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index f92f0127f9c..c0b236abd55 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -40,7 +40,7 @@ # %% # Load and prepare the data -# -------------------------- +# ------------------------- # The data were collected with an Elekta Neuromag VectorView system # at 1000 Hz and low-pass filtered at 330 Hz. @@ -60,7 +60,7 @@ # %% # Epoch the data and plot evokeds -# -------------------------------- +# ------------------------------- # # We epoch the data around dipole events and apply baseline correction. @@ -79,7 +79,7 @@ # You can see that the burst envelope repeats at approximately 3 Hz. # # 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. @@ -124,7 +124,7 @@ # %% # Fit dipoles -# ------------ +# ----------- # We fit dipoles for each phantom and store them in a list. dip_all, residuals_all = [], [] @@ -135,7 +135,7 @@ residuals_all.append(residual) # %% # Evaluate goodness of fit -# ----------------------- +# ------------------------ # The dipole object stores the goodness of fit (GOF) for each dipole. gof = [dip.gof[0] for dip in dip_all] @@ -149,7 +149,7 @@ # We can see that GOF varies between 50 % and up to 95 %. # # Compare estimated and true dipoles -# ----------------------------------- +# ---------------------------------- actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 200.0 # nAm @@ -190,7 +190,7 @@ # achieving sub-centimeter accuracy (mean position error 2.7mm). # # Visualise estimated and true dipoles -# ------------------------------------- +# ------------------------------------ actual_amp = np.ones(len(dip)) # fake amp, needed to create Dipole instance actual_gof = np.ones(len(dip)) # fake goodness-of-fit (GOF) From 111066ef0a9a8aaf74fc3ee45687740fe030c379 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 11 Apr 2026 15:02:29 +1000 Subject: [PATCH 31/37] final formatting --- .../inverse/80_brainstorm_phantom_elekta.py | 35 +++++++------------ 1 file changed, 13 insertions(+), 22 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index c0b236abd55..36880f64233 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -43,9 +43,8 @@ # ------------------------- # The data were collected with an Elekta Neuromag VectorView system -# at 1000 Hz and low-pass filtered at 330 Hz. -# -# The dataset contains recordings at three current amplitudes (20, 200, and 2000 nAm). +# at 1000 Hz, low-pass filtered at 330 Hz and contains recordings +# at three current amplitudes (20, 200, and 2000 nAm). # Here we load the medium-amplitude condition. data_path = bst_phantom_elekta.data_path(verbose=True) raw_fname = data_path / "kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif" @@ -55,14 +54,13 @@ raw.info["bads"] = ["MEG1933", "MEG2421"] events = find_events(raw, "STI201") - # The first 32 events correspond to dipole activations. # %% # Epoch the data and plot evokeds # ------------------------------- # -# We epoch the data around dipole events and apply baseline correction. +# We epoch and baseline correct the data around the dipole events. bmax = -0.05 tmin, tmax = -0.1, 0.8 @@ -72,11 +70,11 @@ raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False ) -# Plot evoked response for the first clean dipole +# Here we plot the evoked response for the first clean dipole epochs["1"][1:-1].average().plot(time_unit="s") # %% # In this data the phantom was set to produce 20 Hz sinusoidal bursts of current. -# You can see that the burst envelope repeats at approximately 3 Hz. +# The burst envelope repeats at approximately 3 Hz. # # Determine peak activation using Global Field Power (GFP) # -------------------------------------------------------- @@ -86,20 +84,16 @@ evoked_tmp = epochs["1"][1:-1].average() gfp = np.std(evoked_tmp.data, axis=0) times = evoked_tmp.times - # Restrict to first burst window time_mask = (times > 0) & (times <= 0.1) - 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") # %% -# Here we store the evoked data for each dipole at the peak amplitude. +# Here we crop the data at the peak amplitude and store the evoked data for each dipole. evokeds = [] for ii in event_id: evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak) @@ -107,13 +101,10 @@ evokeds.append(evoked) # %% # Next, we need to compute the noise covariance to capture the sensor noise structure. -# # We use the baseline window to estimate covariance. -# # You can explore the covariance tutorial for details: :ref:`tut-compute-covariance`. cov = mne.compute_covariance(epochs, tmax=bmax) - del epochs # delete to save memory # %% # We use a :ref:`sphere head geometry model ` @@ -123,16 +114,15 @@ sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08) # %% -# Fit dipoles +# Dipole fitting # ----------- -# We fit dipoles for each phantom and store them in a list. -dip_all, residuals_all = [], [] +# Finally, we fit dipoles for each phantom and store them in a list. +dip_all = [] for evoked in evokeds: dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1) dip_all.append(dip) - residuals_all.append(residual) # %% # Evaluate goodness of fit # ------------------------ @@ -144,13 +134,14 @@ plt.xlabel("Phantom dipole estimation") plt.ylabel("Goodness of fit (%)") plt.show() - +# # %% -# We can see that GOF varies between 50 % and up to 95 %. +# We can see that GOF varies between simulated dipoles from 50 % up to 95 %. # # Compare estimated and true dipoles # ---------------------------------- +# We get the true dipole positions from the phantoms actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 200.0 # nAm @@ -189,7 +180,7 @@ # The dipole fits closely match the true phantom data, # achieving sub-centimeter accuracy (mean position error 2.7mm). # -# Visualise estimated and true dipoles +# Visualise estimated and true dipole fits # ------------------------------------ actual_amp = np.ones(len(dip)) # fake amp, needed to create Dipole instance From b66d7f03c024eb9f44e1ba477e081f5538f22a0f Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 11 Apr 2026 15:24:51 +1000 Subject: [PATCH 32/37] fix docs fail --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 36880f64233..e7c9b10a3ab 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -115,7 +115,7 @@ # %% # Dipole fitting -# ----------- +# -------------- # Finally, we fit dipoles for each phantom and store them in a list. dip_all = [] @@ -181,7 +181,7 @@ # achieving sub-centimeter accuracy (mean position error 2.7mm). # # Visualise estimated and true dipole fits -# ------------------------------------ +# ---------------------------------------- actual_amp = np.ones(len(dip)) # fake amp, needed to create Dipole instance actual_gof = np.ones(len(dip)) # fake goodness-of-fit (GOF) From 7a62a3490564efb1febe17589bed46c8773a1fb5 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sat, 18 Apr 2026 15:40:54 +1000 Subject: [PATCH 33/37] CI: trigger rerun From 1a9931771d41586fb8da1ad8b59ce56807c6f9d4 Mon Sep 17 00:00:00 2001 From: Carina Date: Thu, 23 Apr 2026 10:48:04 +1000 Subject: [PATCH 34/37] Update 80_brainstorm_phantom_elekta.py add review comments from Dan --- .../inverse/80_brainstorm_phantom_elekta.py | 49 ++++++++++--------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index e7c9b10a3ab..47099c767f2 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -41,10 +41,11 @@ # %% # Load and prepare the data # ------------------------- - +# TODO: convert to text from code comment # The data were collected with an Elekta Neuromag VectorView system -# at 1000 Hz, low-pass filtered at 330 Hz and contains recordings +# 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. data_path = bst_phantom_elekta.data_path(verbose=True) raw_fname = data_path / "kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif" @@ -69,19 +70,22 @@ epochs = mne.Epochs( raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False ) - -# Here we plot the evoked response for the first clean dipole -epochs["1"][1:-1].average().plot(time_unit="s") +# We drop the first and last epoch as they contain artefacts +epochs_clean = epochs[1:-1 +# We select the first simulated dipole for visualisation purposes +epochs_firstdip = epochs_clean["1"] +# 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. +epochs_firstdip.average().plot(time_unit="s") + # %% -# In this data the phantom was set to produce 20 Hz sinusoidal bursts of current. -# The burst envelope repeats at approximately 3 Hz. -# # 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. -evoked_tmp = epochs["1"][1:-1].average() +evoked_tmp = epochs_firstdip.average() gfp = np.std(evoked_tmp.data, axis=0) times = evoked_tmp.times # Restrict to first burst window @@ -96,7 +100,7 @@ # Here we crop the data at the peak amplitude and store the evoked data for each dipole. evokeds = [] for ii in event_id: - evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak) + 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) # %% @@ -104,9 +108,10 @@ # We use the baseline window to estimate covariance. # You can explore the covariance tutorial for details: :ref:`tut-compute-covariance`. -cov = mne.compute_covariance(epochs, tmax=bmax) +cov = mne.compute_covariance(epochs_clean, tmax=bmax) del epochs # delete to save memory # %% +# TODO: explain why this head model is used # We use a :ref:`sphere head geometry model ` # to fit our phantom head model. subjects_dir = data_path @@ -127,6 +132,7 @@ # Evaluate goodness of fit # ------------------------ +# TODO: explain drop in GOF at regular intervals # The dipole object stores the goodness of fit (GOF) for each dipole. gof = [dip.gof[0] for dip in dip_all] colors = ["#E69F00" if val < 60 else "#0072B2" for val in gof] @@ -136,14 +142,15 @@ plt.show() # # %% -# We can see that GOF varies between simulated dipoles from 50 % up to 95 %. -# # Compare estimated and true dipoles # ---------------------------------- +# The dipole fits closely match the true phantom data, +# achieving sub-centimeter accuracy (mean position error 2.7mm). + # We get the true dipole positions from the phantoms actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() -actual_amp = 200.0 # nAm +actual_amp = 100.0 # nAm # estimated dipoles dip_pos = [dip.pos[0] for dip in dip_all] @@ -177,12 +184,13 @@ ax3.set_xlabel("Dipole index") ax3.set_ylabel("Amplitude error (nAm)") # %% -# The dipole fits closely match the true phantom data, -# achieving sub-centimeter accuracy (mean position error 2.7mm). -# -# Visualise estimated and true dipole fits -# ---------------------------------------- +# 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 goodness-of-fit (GOF) # setup dipole objects for true and estimated dipoles @@ -214,9 +222,6 @@ 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) -# %% -# We can see that the dipoles overlap, have approximately the same magnitude -# and point in the same direction. # %% # References From a7c42e7896218adee4bafc65b546f776801e22b0 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 26 Apr 2026 14:07:46 +1000 Subject: [PATCH 35/37] addressed Dan's review --- .../inverse/80_brainstorm_phantom_elekta.py | 95 ++++++++++++++----- 1 file changed, 69 insertions(+), 26 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 47099c767f2..be3c1ebd264 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -41,12 +41,12 @@ # %% # Load and prepare the data # ------------------------- -# TODO: convert to text from code comment # 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) @@ -54,50 +54,62 @@ # Mark known bad channels raw.info["bads"] = ["MEG1933", "MEG2421"] -events = find_events(raw, "STI201") # The first 32 events correspond to dipole activations. +events = find_events(raw, "STI201") + # %% # Epoch the data and plot evokeds # ------------------------------- -# # We epoch and baseline correct the data around the dipole events. +# 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 epoch as they contain artefacts -epochs_clean = epochs[1:-1 + +# We drop the first and last event, it can contains dipole-switching artifacts +epochs_clean = epochs[1:-1] + # We select the first simulated dipole for visualisation purposes epochs_firstdip = epochs_clean["1"] + +# %% # 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. + epochs_firstdip.average().plot(time_unit="s") - + # %% # 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. + +# Get the evoked signal of the first dipole evoked_tmp = epochs_firstdip.average() +# Calculate GFP gfp = np.std(evoked_tmp.data, axis=0) -times = evoked_tmp.times + # Restrict to first burst window -time_mask = (times > 0) & (times <= 0.1) +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") # %% -# Here we crop the data at the peak amplitude and store the evoked data for each dipole. +# 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) @@ -111,9 +123,11 @@ cov = mne.compute_covariance(epochs_clean, tmax=bmax) del epochs # delete to save memory # %% -# TODO: explain why this head model is used +# %% # We use a :ref:`sphere head geometry model ` -# to fit our phantom head 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) @@ -121,8 +135,8 @@ # %% # Dipole fitting # -------------- - # Finally, we fit dipoles for each phantom and store them in a list. + dip_all = [] for evoked in evokeds: @@ -131,20 +145,54 @@ # %% # Evaluate goodness of fit # ------------------------ - -# TODO: explain drop in GOF at regular intervals # The dipole object stores the goodness of fit (GOF) for each dipole. +# Some dipoles have lower GOF because... 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() -# + +# %% +# Dipoles with low goodness of fit +# -------------------------------- +# Why do some dipoles have a low (<60) GOF? +# Here we plot the dipole locations of the dipoles with low GOF. +# The 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)) + +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, +) + +# Plot the position and the orientation of the dipoles with low GOF +fig = mne.viz.plot_dipole_locations( + dipoles=low_dips, mode="arrow", subject=subject, color=(1.0, 0.2, 0.2), fig=fig +) # %% # Compare estimated and true dipoles # ---------------------------------- - # The dipole fits closely match the true phantom data, # achieving sub-centimeter accuracy (mean position error 2.7mm). @@ -152,7 +200,7 @@ actual_pos, actual_ori = mne.dipole.get_phantom_dipoles() actual_amp = 100.0 # nAm -# estimated dipoles +# 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] @@ -186,20 +234,15 @@ # %% # 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 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) -subject = "phantom_otaniemi" -trans = mne.transforms.Transform("head", "mri", np.eye(4)) - fig = mne.viz.plot_alignment( evoked.info, trans, From 3b4200fe42d3341f09386677b28faccb8fe03258 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 26 Apr 2026 15:25:31 +1000 Subject: [PATCH 36/37] prettify and spelling --- .../inverse/80_brainstorm_phantom_elekta.py | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index be3c1ebd264..9e58486621c 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -54,7 +54,7 @@ # Mark known bad channels raw.info["bads"] = ["MEG1933", "MEG2421"] -# The first 32 events correspond to dipole activations. +# The first 32 events correspond to dipole activations events = find_events(raw, "STI201") @@ -71,15 +71,16 @@ raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False ) -# We drop the first and last event, it can contains dipole-switching artifacts +# 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"] # %% -# 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. +# 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. epochs_firstdip.average().plot(time_unit="s") @@ -92,6 +93,7 @@ # Get the evoked signal of the first dipole evoked_tmp = epochs_firstdip.average() + # Calculate GFP gfp = np.std(evoked_tmp.data, axis=0) @@ -116,14 +118,12 @@ evoked = mne.EvokedArray(np.array(evoked.data), evoked.info, tmin=0.0) evokeds.append(evoked) # %% -# Next, we need to compute the noise covariance to capture the sensor noise structure. -# We use the baseline window to estimate covariance. -# You can explore the covariance tutorial for details: :ref:`tut-compute-covariance`. +# 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. @@ -146,7 +146,7 @@ # Evaluate goodness of fit # ------------------------ # The dipole object stores the goodness of fit (GOF) for each dipole. -# Some dipoles have lower GOF because... +# 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) @@ -157,9 +157,9 @@ # %% # Dipoles with low goodness of fit # -------------------------------- -# Why do some dipoles have a low (<60) GOF? -# Here we plot the dipole locations of the dipoles with low GOF. -# The dipoles with low GOF are deep in the brain which might explain +# 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 @@ -174,6 +174,7 @@ 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, @@ -186,7 +187,6 @@ subjects_dir=subjects_dir, ) -# Plot the position and the orientation of the dipoles with low GOF fig = mne.viz.plot_dipole_locations( dipoles=low_dips, mode="arrow", subject=subject, color=(1.0, 0.2, 0.2), fig=fig ) @@ -194,7 +194,7 @@ # Compare estimated and true dipoles # ---------------------------------- # The dipole fits closely match the true phantom data, -# achieving sub-centimeter accuracy (mean position error 2.7mm). +# 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() From d18bdeeec12c2c00e51c318ae3b459ecce080166 Mon Sep 17 00:00:00 2001 From: CarinaFo Date: Sun, 26 Apr 2026 16:08:28 +1000 Subject: [PATCH 37/37] last nitpick --- tutorials/inverse/80_brainstorm_phantom_elekta.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tutorials/inverse/80_brainstorm_phantom_elekta.py b/tutorials/inverse/80_brainstorm_phantom_elekta.py index 9e58486621c..129025d10ac 100644 --- a/tutorials/inverse/80_brainstorm_phantom_elekta.py +++ b/tutorials/inverse/80_brainstorm_phantom_elekta.py @@ -158,7 +158,8 @@ # 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 +# 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.