{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Brainstorm Elekta phantom dataset tutorial\n\nHere we compute the evoked from raw for the Brainstorm Elekta phantom\ntutorial dataset. For comparison, see :footcite:`TadelEtAl2011` and\n[the original Brainstorm tutorial](https://neuroimage.usc.edu/brainstorm/Tutorials/PhantomElekta)_.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Eric Larson <larson.eric.d@gmail.com>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport mne\nfrom mne import find_events, fit_dipole\nfrom mne.datasets import fetch_phantom\nfrom mne.datasets.brainstorm import bst_phantom_elekta\nfrom mne.io import read_raw_fif\n\nprint(__doc__)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The data were collected with an Elekta Neuromag VectorView system at 1000 Hz\nand low-pass filtered at 330 Hz. Here the medium-amplitude (200 nAm) data\nare read to construct instances of :class:`mne.io.Raw`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = bst_phantom_elekta.data_path(verbose=True)\n\nraw_fname = data_path / \"kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif\"\nraw = read_raw_fif(raw_fname)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Data channel array consisted of 204 MEG planor gradiometers,\n102 axial magnetometers, and 3 stimulus channels. Let's get the events\nfor the phantom, where each dipole (1-32) gets its own event:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "events = find_events(raw, \"STI201\")\nraw.plot(events=events)\nraw.info[\"bads\"] = [\"MEG1933\", \"MEG2421\"]  # known bad channels"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The data has strong line frequency (60 Hz and harmonics) and cHPI coil\nnoise (five peaks around 300 Hz). Here, we use only the first 30 seconds\nto save memory:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw.compute_psd(tmax=30).plot(\n    average=False, amplitude=False, picks=\"data\", exclude=\"bads\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Our phantom produces sinusoidal bursts at 20 Hz:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw.plot(events=events)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we epoch our data, average it, and look at the first dipole response.\nThe first peak appears around 3 ms. Because we low-passed at 40 Hz,\nwe can also decimate our data to save memory.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tmin, tmax = -0.1, 0.1\nbmax = -0.05  # Avoid capture filter ringing into baseline\nevent_id = list(range(1, 33))\nepochs = mne.Epochs(\n    raw, events, event_id, tmin, tmax, baseline=(None, bmax), preload=False\n)\nepochs[\"1\"].average().plot(time_unit=\"s\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nLet's use a `sphere head geometry model <eeg_sphere_model>`\nand let's see the coordinate alignment and the sphere location. The phantom\nis properly modeled by a single-shell sphere with origin (0., 0., 0.).\n\nEven though this is a VectorView/TRIUX phantom, we can use the Otaniemi\nphantom subject as a surrogate because the \"head\" surface (hemisphere outer\nshell) has the same geometry for both phantoms, even though the internal\ndipole locations differ. The phantom_otaniemi scan was aligned to the\nphantom's head coordinate frame, so an identity ``trans`` is appropriate\nhere.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "subjects_dir = data_path\nfetch_phantom(\"otaniemi\", subjects_dir=subjects_dir)\nsphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08)\nsubject = \"phantom_otaniemi\"\ntrans = mne.transforms.Transform(\"head\", \"mri\", np.eye(4))\nmne.viz.plot_alignment(\n    epochs.info,\n    subject=subject,\n    show_axes=True,\n    bem=sphere,\n    dig=True,\n    surfaces=(\"head-dense\", \"inner_skull\"),\n    trans=trans,\n    mri_fiducials=True,\n    subjects_dir=subjects_dir,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's do some dipole fits. We first compute the noise covariance,\nthen do the fits for each event_id taking the time instant that maximizes\nthe global field power.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# here we can get away with using method='oas' for speed (faster than \"shrunk\")\n# but in general \"shrunk\" is usually better\ncov = mne.compute_covariance(epochs, tmax=bmax)\nmne.viz.plot_evoked_white(epochs[\"1\"].average(), cov)\n\ndata = []\nt_peak = 0.036  # true for Elekta phantom\nfor ii in event_id:\n    # Avoid the first and last trials -- can contain dipole-switching artifacts\n    evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak)\n    data.append(evoked.data[:, 0])\nevoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.0)\ndel epochs\ndip, residual = fit_dipole(evoked, cov, sphere, n_jobs=None)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Do a quick visualization of how much variance we explained, putting the\ndata and residuals on the same scale (here the \"time points\" are the\n32 dipole peak values that we fit):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(2, 1)\nevoked.plot(axes=axes)\nfor ax in axes:\n    for text in list(ax.texts):\n        text.remove()\n    for line in ax.lines:\n        line.set_color(\"#98df81\")\nresidual.plot(axes=axes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can compare to the actual locations, taking the difference in mm:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "actual_pos, actual_ori = mne.dipole.get_phantom_dipoles()\nactual_amp = 100.0  # nAm\n\nfig, (ax1, ax2, ax3) = plt.subplots(\n    nrows=3, ncols=1, figsize=(6, 7), layout=\"constrained\"\n)\n\ndiffs = 1000 * np.sqrt(np.sum((dip.pos - actual_pos) ** 2, axis=-1))\nprint(f\"mean(position error) = {np.mean(diffs):0.1f} mm\")\nax1.bar(event_id, diffs)\nax1.set_xlabel(\"Dipole index\")\nax1.set_ylabel(\"Loc. error (mm)\")\n\nangles = np.rad2deg(np.arccos(np.abs(np.sum(dip.ori * actual_ori, axis=1))))\nprint(f\"mean(angle error) = {np.mean(angles):0.1f}\u00b0\")\nax2.bar(event_id, angles)\nax2.set_xlabel(\"Dipole index\")\nax2.set_ylabel(\"Angle error (\u00b0)\")\n\namps = actual_amp - dip.amplitude / 1e-9\nprint(f\"mean(abs amplitude error) = {np.mean(np.abs(amps)):0.1f} nAm\")\nax3.bar(event_id, amps)\nax3.set_xlabel(\"Dipole index\")\nax3.set_ylabel(\"Amplitude error (nAm)\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's plot the positions and the orientations of the actual and the estimated\ndipoles\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "actual_amp = np.ones(len(dip))  # fake amp, needed to create Dipole instance\nactual_gof = np.ones(len(dip))  # fake GOF, needed to create Dipole instance\ndip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof)\n\nfig = mne.viz.plot_alignment(\n    evoked.info,\n    trans,\n    subject,\n    bem=sphere,\n    surfaces={\"head-dense\": 0.2},\n    coord_frame=\"head\",\n    meg=\"helmet\",\n    show_axes=True,\n    subjects_dir=subjects_dir,\n)\n\n# Plot the position and the orientation of the actual dipole\nfig = mne.viz.plot_dipole_locations(\n    dipoles=dip_true, mode=\"arrow\", subject=subject, color=(0.0, 0.0, 0.0), fig=fig\n)\n\n# Plot the position and the orientation of the estimated dipole\nfig = mne.viz.plot_dipole_locations(\n    dipoles=dip, mode=\"arrow\", subject=subject, color=(0.2, 1.0, 0.5), fig=fig\n)\n\nmne.viz.set_3d_view(figure=fig, azimuth=70, elevation=80, distance=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n.. footbibliography::\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}