{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Head model and forward computation\n\nThe aim of this tutorial is to be a getting started for forward computation.\n\nFor more extensive details and presentation of the general concepts for forward\nmodeling, see `ch_forward`.\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 mne\nfrom mne.datasets import sample\n\ndata_path = sample.data_path()\n\n# the raw file containing the channel location + types\nsample_dir = data_path / \"MEG\" / \"sample\"\nraw_fname = sample_dir / \"sample_audvis_raw.fif\"\n# The paths to Freesurfer reconstructions\nsubjects_dir = data_path / \"subjects\"\nsubject = \"sample\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Computing the forward operator\n\nTo compute a forward operator we need:\n\n   - a ``-trans.fif`` file that contains the coregistration info.\n   - a source space\n   - the :term:`BEM` surfaces\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute and visualize BEM surfaces\n\nThe :term:`BEM` surfaces are the triangulations of the interfaces between\ndifferent tissues needed for forward computation. These surfaces are for\nexample the inner skull surface, the outer skull surface and the outer skin\nsurface, a.k.a. scalp surface.\n\nComputing the BEM surfaces requires FreeSurfer and makes use of\nthe command-line tools `mne watershed_bem` or `mne flash_bem`, or\nthe related functions :func:`mne.bem.make_watershed_bem` or\n:func:`mne.bem.make_flash_bem`.\n\nHere we'll assume it's already computed. It takes a few minutes per subject.\n\nFor EEG we use 3 layers (inner skull, outer skull, and skin) while for\nMEG 1 layer (inner skull) is enough.\n\nLet's look at these surfaces. The function :func:`mne.viz.plot_bem`\nassumes that you have the ``bem`` folder of your subject's FreeSurfer\nreconstruction, containing the necessary surface files. Here we use a smaller\nthan default subset of ``slices`` for speed.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_bem_kwargs = dict(\n    subject=subject,\n    subjects_dir=subjects_dir,\n    brain_surfaces=\"white\",\n    orientation=\"coronal\",\n    slices=[50, 100, 150, 200],\n)\n\nmne.viz.plot_bem(**plot_bem_kwargs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualizing the coregistration\n\nThe coregistration is the operation that allows to position the head and the\nsensors in a common coordinate system. In the MNE software the transformation\nto align the head and the sensors in stored in a so-called **trans file**.\nIt is a FIF file that ends with ``-trans.fif``. It can be obtained with\n:func:`mne.gui.coregistration` (or its convenient command line\nequivalent `mne coreg`), or mrilab if you're using a Neuromag\nsystem.\n\nHere we assume the coregistration is done, so we just visually check the\nalignment with the following code. See `creating-trans` for instructions\non creating the ``-trans.fif`` file interactively.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# The transformation file obtained by coregistration\ntrans = sample_dir / \"sample_audvis_raw-trans.fif\"\n\ninfo = mne.io.read_info(raw_fname)\n# Here we look at the dense head, which isn't used for BEM computations but\n# is useful for coregistration.\nmne.viz.plot_alignment(\n    info,\n    trans,\n    subject=subject,\n    dig=True,\n    meg=[\"helmet\", \"sensors\"],\n    subjects_dir=subjects_dir,\n    surfaces=\"head-dense\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n## Compute Source Space\n\nThe source space defines the position and orientation of the candidate source\nlocations. There are two types of source spaces:\n\n- **surface-based** source space when the candidates are confined to a\n  surface.\n\n- **volumetric or discrete** source space when the candidates are discrete,\n  arbitrarily located source points bounded by the surface.\n\n**Surface-based** source space is computed using\n:func:`mne.setup_source_space`, while **volumetric** source space is computed\nusing :func:`mne.setup_volume_source_space`.\n\nWe will now compute a surface-based source space with an ``'oct4'``\nresolution. See `setting_up_source_space` for details on source space\ndefinition and spacing parameter.\n\n<div class=\"alert alert-danger\"><h4>Warning</h4><p>``'oct4'`` is used here just for speed, for real analyses the recommended\n    spacing is ``'oct6'``.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "src = mne.setup_source_space(\n    subject, spacing=\"oct4\", add_dist=\"patch\", subjects_dir=subjects_dir\n)\nprint(src)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The surface based source space ``src`` contains two parts, one for the left\nhemisphere (258 locations) and one for the right hemisphere (258\nlocations). Sources can be visualized on top of the BEM surfaces in purple.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mne.viz.plot_bem(src=src, **plot_bem_kwargs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To compute a volume based source space defined with a grid of candidate\ndipoles inside a sphere of radius 90mm centered at (0.0, 0.0, 40.0) mm\nyou can use the following code.\nObviously here, the sphere is not perfect. It is not restricted to the\nbrain and it can miss some parts of the cortex.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sphere = (0.0, 0.0, 0.04, 0.09)\nvol_src = mne.setup_volume_source_space(\n    subject,\n    subjects_dir=subjects_dir,\n    sphere=sphere,\n    sphere_units=\"m\",\n    add_interpolator=False,\n)  # just for speed!\nprint(vol_src)\n\nmne.viz.plot_bem(src=vol_src, **plot_bem_kwargs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To compute a volume based source space defined with a grid of candidate\ndipoles inside the brain (requires the :term:`BEM` surfaces) you can use the\nfollowing.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "surface = subjects_dir / subject / \"bem\" / \"inner_skull.surf\"\nvol_src = mne.setup_volume_source_space(\n    subject, subjects_dir=subjects_dir, surface=surface, add_interpolator=False\n)  # Just for speed!\nprint(vol_src)\n\nmne.viz.plot_bem(src=vol_src, **plot_bem_kwargs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>Some sources may appear to be outside the BEM inner skull contour.\n          This is because the ``slices`` are decimated for plotting here.\n          Each slice in the figure actually represents several MRI slices,\n          but only the MRI voxels and BEM boundaries for a single (midpoint\n          of the given slice range) slice are shown, whereas the source space\n          points plotted on that midpoint slice consist of all points\n          for which that slice (out of all slices shown) was the closest.</p></div>\n\nNow let's see how to view all sources in 3D.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = mne.viz.plot_alignment(\n    subject=subject,\n    subjects_dir=subjects_dir,\n    surfaces=\"white\",\n    coord_frame=\"mri\",\n    src=src,\n)\nmne.viz.set_3d_view(\n    fig,\n    azimuth=173.78,\n    elevation=101.75,\n    distance=0.30,\n    focalpoint=(-0.03, -0.01, 0.03),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n## Compute forward solution\n\nWe can now compute the forward solution.\nTo reduce computation we'll just compute a single layer BEM (just inner\nskull) that can then be used for MEG (not EEG).\nWe specify if we want a one-layer or a three-layer BEM using the\n``conductivity`` parameter.\nThe BEM solution requires a BEM model which describes the geometry\nof the head the conductivities of the different tissues.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "conductivity = (0.3,)  # for single layer\n# conductivity = (0.3, 0.006, 0.3)  # for three layers\nmodel = mne.make_bem_model(\n    subject=\"sample\", ico=4, conductivity=conductivity, subjects_dir=subjects_dir\n)\nbem = mne.make_bem_solution(model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that the :term:`BEM` does not involve any use of the trans file. The BEM\nonly depends on the head geometry and conductivities.\nIt is therefore independent from the MEG data and the head position.\n\nLet's now compute the forward operator, commonly referred to as the\ngain or leadfield matrix.\nSee :func:`mne.make_forward_solution` for details on the meaning of each\nparameter.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fwd = mne.make_forward_solution(\n    raw_fname,\n    trans=trans,\n    src=src,\n    bem=bem,\n    meg=True,\n    eeg=False,\n    mindist=5.0,\n    n_jobs=None,\n    verbose=True,\n)\nprint(fwd)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-danger\"><h4>Warning</h4><p>Forward computation can remove vertices that are too close to (or outside)\n   the inner skull surface. For example, here we have gone from 516 to 474\n   vertices in use. For many functions, such as\n   :func:`mne.compute_source_morph`, it is important to pass ``fwd['src']``\n   or ``inv['src']`` so that this removal is adequately accounted for.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"Before: {src}\")\nprint(f'After:  {fwd[\"src\"]}')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can explore the content of ``fwd`` to access the numpy array that contains\nthe gain matrix.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "leadfield = fwd[\"sol\"][\"data\"]\nprint(f\"Leadfield size : {leadfield.shape[0]} sensors x {leadfield.shape[1]} dipoles\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To extract the numpy array containing the forward operator corresponding to\nthe source space ``fwd['src']`` with cortical orientation constraint\nwe can use the following:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fwd_fixed = mne.convert_forward_solution(\n    fwd, surf_ori=True, force_fixed=True, use_cps=True\n)\nleadfield = fwd_fixed[\"sol\"][\"data\"]\nprint(f\"Leadfield size : {leadfield.shape[0]} sensors x {leadfield.shape[1]} dipoles\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is equivalent to the following code that explicitly applies the\nforward operator to a source estimate composed of the identity operator\n(which we omit here because it uses a lot of memory)::\n\n    >>> import numpy as np\n    >>> n_dipoles = leadfield.shape[1]\n    >>> vertices = [src_hemi['vertno'] for src_hemi in fwd_fixed['src']]\n    >>> stc = mne.SourceEstimate(1e-9 * np.eye(n_dipoles), vertices)\n    >>> leadfield = mne.apply_forward(fwd_fixed, stc, info).data / 1e-9\n\nTo save to disk a forward solution you can use\n:func:`mne.write_forward_solution` and to read it back from disk\n:func:`mne.read_forward_solution`. Don't forget that FIF files containing\nforward solution should end with :file:`-fwd.fif`.\n\nTo get a fixed-orientation forward solution, use\n:func:`mne.convert_forward_solution` to convert the free-orientation\nsolution to (surface-oriented) fixed orientation.\n\n## Exercise\n\nBy looking at `ex-sensitivity-maps`\nplot the sensitivity maps for EEG and compare it with the MEG, can you\njustify the claims that:\n\n  - MEG is not sensitive to radial sources\n  - EEG is more sensitive to deep sources\n\nHow will the MEG sensitivity maps and histograms change if you use a free\ninstead if a fixed/surface oriented orientation?\n\nTry this changing the mode parameter in :func:`mne.sensitivity_map`\naccordingly. Why don't we see any dipoles on the gyri?\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
}