{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Source localization with equivalent current dipole (ECD) fit\n\nThis shows how to fit a dipole :footcite:`Sarvas1987` using MNE-Python.\n\nFor a comparison of fits between MNE-C and MNE-Python, see\n[this gist](https://gist.github.com/larsoner/ca55f791200fe1dc3dd2)_.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The MNE-Python contributors.\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\nfrom nilearn.datasets import load_mni152_template\nfrom nilearn.plotting import plot_anat\n\nimport mne\nfrom mne.evoked import combine_evoked\nfrom mne.forward import make_forward_dipole\nfrom mne.simulation import simulate_evoked\n\ndata_path = mne.datasets.sample.data_path()\nsubjects_dir = data_path / \"subjects\"\nfname_ave = data_path / \"MEG\" / \"sample\" / \"sample_audvis-ave.fif\"\nfname_cov = data_path / \"MEG\" / \"sample\" / \"sample_audvis-cov.fif\"\nfname_bem = subjects_dir / \"sample\" / \"bem\" / \"sample-5120-bem-sol.fif\"\nfname_trans = data_path / \"MEG\" / \"sample\" / \"sample_audvis_raw-trans.fif\"\nfname_surf_lh = subjects_dir / \"sample\" / \"surf\" / \"lh.white\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's localize the N100m (using MEG only)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "evoked = mne.read_evokeds(fname_ave, condition=\"Right Auditory\", baseline=(None, 0))\nevoked.pick(picks=\"meg\")\nevoked_full = evoked.copy()\nevoked.crop(0.07, 0.08)\n\n# Fit a dipole\ndip = mne.fit_dipole(evoked, fname_cov, fname_bem, fname_trans)[0]\n\n# Plot the result in 3D brain with the MRI image.\ndip.plot_locations(fname_trans, \"sample\", subjects_dir, mode=\"orthoview\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also plot the result using outlines of the head and brain.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "color = [\"k\"] * len(dip)\ncolor[np.argmax(dip.gof)] = \"r\"\ndip.plot_locations(fname_trans, \"sample\", subjects_dir, mode=\"outlines\", color=color)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the result in 3D brain with the MRI image using Nilearn\nIn MRI coordinates and in MNI coordinates (template brain)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "subject = \"sample\"\nmni_pos = dip.to_mni(subject=subject, trans=fname_trans, subjects_dir=subjects_dir)\n\nmri_pos = dip.to_mri(subject=subject, trans=fname_trans, subjects_dir=subjects_dir)\n\n# Find an anatomical label for the best fitted dipole\nbest_dip_idx = dip.gof.argmax()\nlabel = dip.to_volume_labels(\n    fname_trans, subject=subject, subjects_dir=subjects_dir, aseg=\"aparc.a2009s+aseg\"\n)[best_dip_idx]\n\n# Draw dipole position on MRI scan and add anatomical label from parcellation\nt1_fname = subjects_dir / subject / \"mri\" / \"T1.mgz\"\nfig_T1 = plot_anat(t1_fname, cut_coords=mri_pos[0], title=f\"Dipole location: {label}\")\n\ntry:\n    template = load_mni152_template(resolution=1)\nexcept TypeError:  # in nilearn < 0.8.1 this did not exist\n    template = load_mni152_template()\nfig_template = plot_anat(\n    template, cut_coords=mni_pos[0], title=\"Dipole loc. (MNI Space)\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Calculate and visualise magnetic field predicted by dipole with maximum GOF\nand compare to the measured data, highlighting the ipsilateral (right) source\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fwd, stc = make_forward_dipole(dip, fname_bem, evoked.info, fname_trans)\npred_evoked = simulate_evoked(fwd, stc, evoked.info, cov=None, nave=np.inf)\n\n# find time point with highest GOF to plot\nbest_idx = np.argmax(dip.gof)\nbest_time = dip.times[best_idx]\nprint(\n    f\"Highest GOF {dip.gof[best_idx]:0.1f}% at t={best_time * 1000:0.1f} ms with \"\n    f\"confidence volume {dip.conf['vol'][best_idx] * 100**3:0.1f} cm^3\"\n)\n# remember to create a subplot for the colorbar\nfig, axes = plt.subplots(\n    nrows=1,\n    ncols=4,\n    figsize=[10.0, 3.4],\n    gridspec_kw=dict(width_ratios=[1, 1, 1, 0.1], top=0.85),\n    layout=\"constrained\",\n)\nvmin, vmax = -400, 400  # make sure each plot has same colour range\n\n# first plot the topography at the time of the best fitting (single) dipole\nplot_params = dict(times=best_time, ch_type=\"mag\", outlines=\"head\", colorbar=False)\nevoked.plot_topomap(time_format=\"Measured field\", axes=axes[0], **plot_params)\n\n# compare this to the predicted field\npred_evoked.plot_topomap(time_format=\"Predicted field\", axes=axes[1], **plot_params)\n\n# Subtract predicted from measured data (apply equal weights)\ndiff = combine_evoked([evoked, pred_evoked], weights=[1, -1])\nplot_params[\"colorbar\"] = True\ndiff.plot_topomap(time_format=\"Difference\", axes=axes[2:], **plot_params)\nfig.suptitle(\n    f\"Comparison of measured and predicted fields at {best_time * 1000:.0f} ms\",\n    fontsize=16,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Estimate the time course of a single dipole with fixed position and\norientation (the one that maximized GOF) over the entire interval\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dip_fixed = mne.fit_dipole(\n    evoked_full,\n    fname_cov,\n    fname_bem,\n    fname_trans,\n    pos=dip.pos[best_idx],\n    ori=dip.ori[best_idx],\n)[0]\ndip_fixed.plot()"
      ]
    },
    {
      "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
}