4a-Inverse_source_localization_dipole_fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Source localization with single dipole fit\n",
    "\n",
    "\n",
    "The aim of this lecture is to show you how to do sequential and\n",
    "fixed dipole fitting with MNE-Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import mne\n",
    "\n",
    "mne.set_log_level('WARNING')\n",
    "\n",
    "# Change the following path to where the folder ds000117-practical is on your disk\n",
    "data_path = os.path.expanduser(\"~/work/data/ds000117-practical/\")\n",
    "\n",
    "raw_fname = os.path.join(data_path,\n",
    "    'derivatives/meg_derivatives/sub-01/ses-meg/meg/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg.fif')\n",
    "\n",
    "epochs_fname = raw_fname.replace('_meg.fif', '-epo.fif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjects_dir = os.path.join(data_path, 'derivatives/freesurfer-reconall/')\n",
    "subject = 'sub-01'\n",
    "trans_fname = os.path.join(data_path,\n",
    "    'derivatives/meg_derivatives/sub-01/ses-meg/meg/sub-01-trans.fif')\n",
    "bem_fname = os.path.join(data_path,\n",
    "    'derivatives/meg_derivatives/sub-01/ses-meg/meg/sub-01-bem.fif')\n",
    "fname_surf_lh = os.path.join(subjects_dir, subject, 'surf', 'lh.white')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Check that the head geometry and sensor alignment is correct</li>\n",
    "    </ul>\n",
    "</div>\n",
    "\n",
    "You will need to use the `mne.viz.plot_alignment` function that can take as input a `bem` parameter of type `ConductorModel`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.viz.plot_alignment?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ConductorModel  |  BEM (1 layer)>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bem = mne.bem.read_bem_solution(bem_fname)\n",
    "bem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Could not find the surface for head in the provided BEM model, looking in the subject directory.\n",
      "Using outer_skin.surf for head surface.\n",
      "Getting helmet for system 306m\n"
     ]
    }
   ],
   "source": [
    "# TODO\n",
    "info = mne.io.read_info(epochs_fname)\n",
    "fig = mne.viz.plot_alignment(info, trans_fname, subject=subject, dig=True,\n",
    "                             subjects_dir=subjects_dir, bem=bem, verbose=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's localize the N170m (using MEG only) using dipole fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs = mne.read_epochs(epochs_fname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<EpochsFIF  |   79 events (all good), -0.2 - 2 sec, baseline [-0.2, 0], ~129.1 MB, data loaded,\n",
       " 'face/famous/first': 13\n",
       " 'face/famous/immediate': 3\n",
       " 'face/famous/long': 6\n",
       " 'face/unfamiliar/first': 17\n",
       " 'face/unfamiliar/immediate': 4\n",
       " 'face/unfamiliar/long': 6\n",
       " 'scrambled/first': 15\n",
       " 'scrambled/immediate': 9\n",
       " 'scrambled/long': 6>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epochs.pick_types(meg=True, eeg=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov = mne.compute_covariance(epochs, rank='info')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked_face = epochs['face'].average()\n",
    "evoked_scrambled = epochs['scrambled'].average()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "contrast = mne.combine_evoked([evoked_face, evoked_scrambled], [0.5, -0.5])\n",
    "contrast.crop(None, 0.2)\n",
    "contrast.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Dipole  |  n_times : 6, tmin : 0.173, tmax : 0.190>\n"
     ]
    }
   ],
   "source": [
    "# Fit a dipole using a sequential (time-varying position and orientation) fit\n",
    "contrast_crop = contrast.copy().crop(0.175, 0.190)\n",
    "dip, residual = mne.fit_dipole(contrast_crop, cov, bem_fname,\n",
    "                               trans_fname)\n",
    "print(dip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[33.67613685 36.70061657 35.58746073 34.71198179 38.51052585 41.45988139]\n"
     ]
    }
   ],
   "source": [
    "# Look at our result\n",
    "print(dip.gof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "dip.plot_locations(subject=subject, trans=trans_fname,\n",
    "                   subjects_dir=subjects_dir, mode='orthoview');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Try to improve the dipole fit by using a subselection of channels [0] </li>\n",
    "      <li>What is the improvement in max GOF from using a subselection of channels?</li>     \n",
    "      <li>Should you also try to maximize GOF?</li>\n",
    "    </ul>\n",
    "</div>\n",
    "\n",
    "Tips and tricks:\n",
    "\n",
    "   - [0] Subselect channels (think about ``mne.read_selection(...)``, ``evoked.pick_channels(...)``)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Dipole  |  n_times : 6, tmin : 0.173, tmax : 0.190>\n"
     ]
    }
   ],
   "source": [
    "# TODO\n",
    "selection = mne.read_selection('Left', info=contrast.info)\n",
    "# Fit a dipole using a sequential (time-varying position and orientation) fit\n",
    "dip, residual = \\\n",
    "    mne.fit_dipole(contrast_crop.copy().pick_channels(selection),\n",
    "                   cov, bem_fname, trans_fname)\n",
    "print(dip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[38.15096387 20.45112159 47.88457627 48.48499646 19.37190885 48.05207802]\n"
     ]
    }
   ],
   "source": [
    "print(dip.gof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "dip.plot_locations(subject=subject, trans=trans_fname,\n",
    "                   subjects_dir=subjects_dir, mode='orthoview');"
   ]
  }
 ],
 "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}