{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Generate a functional label from source estimates\n\n\nThreshold source estimates and produce a functional label. The label\nis typically the region of interest that contains high values.\nHere we compare the average time course in the anatomical label obtained\nby FreeSurfer segmentation and the average time course from the\nfunctional label. As expected the time course in the functional\nlabel yields higher values.\n\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Author: Luke Bloy <luke.bloy@gmail.com>\n#         Alex Gramfort <alexandre.gramfort@telecom-paristech.fr>\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne.minimum_norm import read_inverse_operator, apply_inverse\nfrom mne.datasets import sample\n\nprint(__doc__)\n\ndata_path = sample.data_path()\nsubjects_dir = data_path + '/subjects'\nfname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'\nfname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif'\nsubjects_dir = data_path + '/subjects'\nsubject = 'sample'\n\nsnr = 3.0\nlambda2 = 1.0 / snr ** 2\nmethod = \"dSPM\"  # use dSPM method (could also be MNE or sLORETA)\n\n# Compute a label/ROI based on the peak power between 80 and 120 ms.\n# The label bankssts-lh is used for the comparison.\naparc_label_name = 'bankssts-lh'\ntmin, tmax = 0.080, 0.120\n\n# Load data\nevoked = mne.read_evokeds(fname_evoked, condition=0, baseline=(None, 0))\ninverse_operator = read_inverse_operator(fname_inv)\nsrc = inverse_operator['src']  # get the source space\n\n# Compute inverse solution\nstc = apply_inverse(evoked, inverse_operator, lambda2, method,\n                    pick_ori='normal')\n\n# Make an STC in the time interval of interest and take the mean\nstc_mean = stc.copy().crop(tmin, tmax).mean()\n\n# use the stc_mean to generate a functional label\n# region growing is halted at 60% of the peak value within the\n# anatomical label / ROI specified by aparc_label_name\nlabel = mne.read_labels_from_annot(subject, parc='aparc',\n                                   subjects_dir=subjects_dir,\n                                   regexp=aparc_label_name)[0]\nstc_mean_label = stc_mean.in_label(label)\ndata = np.abs(stc_mean_label.data)\nstc_mean_label.data[data < 0.6 * np.max(data)] = 0.\n\nfunc_labels, _ = mne.stc_to_label(stc_mean_label, src=src, smooth=True,\n                                  subjects_dir=subjects_dir, connected=True)\n\n# take first as func_labels are ordered based on maximum values in stc\nfunc_label = func_labels[0]\n\n# load the anatomical ROI for comparison\nanat_label = mne.read_labels_from_annot(subject, parc='aparc',\n                                        subjects_dir=subjects_dir,\n                                        regexp=aparc_label_name)[0]\n\n# extract the anatomical time course for each label\nstc_anat_label = stc.in_label(anat_label)\npca_anat = stc.extract_label_time_course(anat_label, src, mode='pca_flip')[0]\n\nstc_func_label = stc.in_label(func_label)\npca_func = stc.extract_label_time_course(func_label, src, mode='pca_flip')[0]\n\n# flip the pca so that the max power between tmin and tmax is positive\npca_anat *= np.sign(pca_anat[np.argmax(np.abs(pca_anat))])\npca_func *= np.sign(pca_func[np.argmax(np.abs(pca_anat))])"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "plot the time courses....\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "plt.figure()\nplt.plot(1e3 * stc_anat_label.times, pca_anat, 'k',\n         label='Anatomical %s' % aparc_label_name)\nplt.plot(1e3 * stc_func_label.times, pca_func, 'b',\n         label='Functional %s' % aparc_label_name)\nplt.legend()\nplt.show()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "plot brain in 3D with PySurfer if available\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "brain = stc_mean.plot(hemi='lh', subjects_dir=subjects_dir)\nbrain.show_view('lateral')\n\n# show both labels\nbrain.add_label(anat_label, borders=True, color='k')\nbrain.add_label(func_label, borders=True, color='b')"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }
  ], 
  "metadata": {
    "kernelspec": {
      "display_name": "Python 2", 
      "name": "python2", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "2.7.13", 
      "pygments_lexer": "ipython2", 
      "codemirror_mode": {
        "version": 2, 
        "name": "ipython"
      }
    }
  }
}