{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# 4D Neuroimaging/BTi phantom dataset tutorial\n\nHere we read 4DBTi epochs data obtained with a spherical phantom\nusing four different dipole locations. For each condition we\ncompute evoked data and compute dipole fits.\n\nData are provided by Jean-Michel Badier from MEG center in Marseille, France.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Alex Gramfort <alexandre.gramfort@inria.fr>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os.path as op\n\nimport numpy as np\n\nimport mne\nfrom mne.datasets import phantom_4dbti"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read data and compute a dipole fit at the peak of the evoked response\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = phantom_4dbti.data_path()\nraw_fname = op.join(data_path, \"{}/e,rfhp1.0Hz\")\n\ndipoles = list()\nsphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.080)\n\nt0 = 0.07  # peak of the response\n\npos = np.empty((4, 3))\nori = np.empty((4, 3))\n\nfor ii in range(4):\n    raw = mne.io.read_raw_bti(\n        raw_fname.format(\n            ii + 1,\n        ),\n        rename_channels=False,\n        preload=True,\n    )\n    raw.info[\"bads\"] = [\"A173\", \"A213\", \"A232\"]\n    events = mne.find_events(raw, \"TRIGGER\", mask=4350, mask_type=\"not_and\")\n    epochs = mne.Epochs(\n        raw, events=events, event_id=8192, tmin=-0.2, tmax=0.4, preload=True\n    )\n    evoked = epochs.average()\n    evoked.plot(time_unit=\"s\")\n    cov = mne.compute_covariance(epochs, tmax=0.0)\n    dip = mne.fit_dipole(evoked.copy().crop(t0, t0), cov, sphere)[0]\n    pos[ii] = dip.pos[0]\n    ori[ii] = dip.ori[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute localisation errors\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "actual_pos = 0.01 * np.array(\n    [[0.16, 1.61, 5.13], [0.17, 1.35, 4.15], [0.16, 1.05, 3.19], [0.13, 0.80, 2.26]]\n)\nactual_pos = np.dot(actual_pos, [[0, 1, 0], [-1, 0, 0], [0, 0, 1]])\n\nerrors = 1e3 * np.linalg.norm(actual_pos - pos, axis=1)\nprint(f\"errors (mm) : {errors}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the dipoles in 3D\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 = mne.Dipole(dip.times, pos, actual_amp, ori, actual_gof)\ndip_true = mne.Dipole(dip.times, actual_pos, actual_amp, ori, actual_gof)\n\nfig = mne.viz.plot_alignment(evoked.info, bem=sphere, surfaces=[])\n\n# Plot the position of the actual dipole\nfig = mne.viz.plot_dipole_locations(\n    dipoles=dip_true, mode=\"sphere\", color=(1.0, 0.0, 0.0), fig=fig\n)\n# Plot the position of the estimated dipole\nfig = mne.viz.plot_dipole_locations(\n    dipoles=dip, mode=\"sphere\", color=(1.0, 1.0, 0.0), fig=fig\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
}