test_3d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test notebook for visualizing in 3D, Neurohackweek 2017\n",
    "\n",
    "- Cloned from the demo introductory_3d.ipynb notebook\n",
    "- Initial functionality based on [IPyvolume](https://ipyvolume.readthedocs.io/en/latest/index.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "load packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import nibabel as nb\n",
    "import ipyvolume\n",
    "import ipyvolume.pylab as p3\n",
    "import numpy as np\n",
    "import os\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Test 3D volume rendering:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ely/Dropbox/SCIENCE/Software/anaconda3/anaconda/lib/python3.6/site-packages/ipyvolume/serialize.py:66: RuntimeWarning: invalid value encountered in true_divide\n",
      "  gradient = gradient / np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e688e3fef21447b5a6a2b76ac13f7917",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# before running, set the global variable NHW to your NHW_visualization git directory\n",
    "# e.g. add 'export NHW=/Users/ely/Desktop/nhw2017/NHW_visualization' to your ~/.bashrc file\n",
    "nhw = os.environ['NHW'] # set the global variable NHW to your NHW_visualization git directory\n",
    "V = nb.load(nhw + '/data/100610/MNINonLinear/T1w_restore_brain.nii.gz').get_data()\n",
    "#V = nb.load('/data/ds000114/derivatives/fmriprep/sub-01/anat/sub-01_t1w_preproc.nii.gz').get_data()\n",
    "ipyvolume.quickvolshow(V, level=[0.25, 0.5, 0.75], opacity=0.03, level_width=0.1, \n",
    "                       width=600, height=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('/tmp/T1.nii.gz', <http.client.HTTPMessage at 0x1274fa710>)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from urllib.request import urlretrieve\n",
    "urlretrieve('https://osf.io/7jbpr/download?version=1', '/tmp/T1.nii.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimensions (176, 256, 256)\n",
      "Pixel dimensions (0.99999803, 1.0, 1.0)\n"
     ]
    }
   ],
   "source": [
    "img = nb.load(nhw + '/data/T1.nii.gz')\n",
    "#img = nb.load('/tmp/T1.nii.gz')\n",
    "#img = nb.load(nhw + '/data/100610/MNINonLinear/T1w_restore_brain.nii.gz')\n",
    "print(\"Dimensions\", img.shape)\n",
    "print(\"Pixel dimensions\", img.header.get_zooms())\n",
    "V = img.get_data()\n",
    "#print(\"Largest dimension\", np.max(img.shape))\n",
    "#print(np.max(img.shape)-img.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ely/Dropbox/SCIENCE/Software/anaconda3/anaconda/lib/python3.6/site-packages/ipyvolume/serialize.py:66: RuntimeWarning: invalid value encountered in true_divide\n",
      "  gradient = gradient / np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "db23180f7d9c46f6adf6f073c30b8f0c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = p3.figure(width=600, height=600)\n",
    "fig.style = {'axes': {'color': 'black',\n",
    "  'label': {'color': 'black'},\n",
    "  'ticklabel': {'color': 'black'},\n",
    "  'visible': False},\n",
    " 'background-color': 'white',\n",
    " 'box': {'visible': False}}\n",
    "fig.camera_fov = 1\n",
    "fig.xlim = (-100, 100)\n",
    "fig.ylim = (-100, 100)\n",
    "fig.zlim = (-100, 100)\n",
    "\n",
    "fig.anglex = 0\n",
    "fig.angley = 0\n",
    "fig.anglez = -np.pi/2\n",
    "vol = p3.volshow(V)\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<bound method HasTraits.trait_names of Figure(anglez=-1.5707963267948966, camera_center=[0.0, 0.0, 0.0], camera_fov=1.0, data_max=1611.9046031311154, height=600, matrix_projection=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], matrix_world=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], style={'axes': {'color': 'black', 'label': {'color': 'black'}, 'ticklabel': {'color': 'black'}, 'visible': False}, 'background-color': 'white', 'box': {'visible': False}}, tf=TransferFunctionWidgetJs3(level3=0.9), volume_data=array([[[ 24.00578153,  25.99806463,  28.99878732, ...,  28.99878732,\n",
       "          24.98962504,   0.        ],\n",
       "        [ 28.99878732,  30.99107042,  33.99179311, ...,  21.98890235,\n",
       "          27.00650422,   0.        ],\n",
       "        [ 33.00794961,  33.00794961,  21.00505884, ...,  30.00722691,\n",
       "          24.98962504,   0.        ],\n",
       "        ..., \n",
       "        [ 16.01205305,  18.00433615,  24.98962504, ...,  36.0086723 ,\n",
       "          21.00505884,   0.        ],\n",
       "        [ 24.00578153,  21.98890235,  16.99589656, ...,  24.00578153,\n",
       "          24.00578153,   0.        ],\n",
       "        [  0.        ,   0.        ,   0.        , ...,   0.        ,\n",
       "           0.        ,   0.        ]],\n",
       "\n",
       "       [[ 18.00433615,  28.99878732,  24.00578153, ...,  30.99107042,\n",
       "          21.98890235,   0.        ],\n",
       "        [ 16.99589656,  30.00722691,  24.00578153, ...,  21.98890235,\n",
       "          25.99806463,   0.        ],\n",
       "        [ 28.99878732,  25.99806463,  15.00361346, ...,  18.00433615,\n",
       "          28.99878732,   0.        ],\n",
       "        ..., \n",
       "        [ 21.00505884,  31.99951001,  22.99734194, ...,  35.0002327 ,\n",
       "          22.99734194,   0.        ],\n",
       "        [ 27.99034773,  19.99661925,  18.00433615, ...,  27.99034773,\n",
       "          16.01205305,   0.        ],\n",
       "        [  0.        ,   0.        ,   0.        , ...,   0.        ,\n",
       "           0.        ,   0.        ]],\n",
       "\n",
       "       [[ 22.99734194,  16.99589656,  36.9925158 , ...,  16.99589656,\n",
       "          24.00578153,   0.        ],\n",
       "        [ 21.98890235,  21.00505884,  39.99323849, ...,  13.99517386,\n",
       "          30.00722691,   0.        ],\n",
       "        [ 30.00722691,  21.98890235,  25.99806463, ...,  21.00505884,\n",
       "          27.99034773,   0.        ],\n",
       "        ..., \n",
       "        [ 19.99661925,  19.99661925,  27.00650422, ...,  27.99034773,\n",
       "          33.99179311,   0.        ],\n",
       "        [ 19.99661925,  30.00722691,  33.00794961, ...,  27.99034773,\n",
       "          24.00578153,   0.        ],\n",
       "        [  0.        ,   0.        ,   0.        , ...,   0.        ,\n",
       "           0.        ,   0.        ]],\n",
       "\n",
       "       ..., \n",
       "       [[ 24.00578153,  22.99734194,  33.00794961, ...,  25.99806463,\n",
       "          18.00433615,   0.        ],\n",
       "        [ 22.99734194,  18.98817965,  16.99589656, ...,  27.00650422,\n",
       "          19.99661925,   0.        ],\n",
       "        [ 18.00433615,  27.99034773,  19.99661925, ...,  24.98962504,\n",
       "          18.00433615,   0.        ],\n",
       "        ..., \n",
       "        [ 15.00361346,  21.98890235,  27.99034773, ...,  25.99806463,\n",
       "          24.98962504,   0.        ],\n",
       "        [ 18.98817965,  21.98890235,  19.99661925, ...,  18.00433615,\n",
       "          27.99034773,   0.        ],\n",
       "        [  0.        ,   0.        ,   0.        , ...,   0.        ,\n",
       "           0.        ,   0.        ]],\n",
       "\n",
       "       [[ 27.00650422,  13.99517386,  19.99661925, ...,  24.98962504,\n",
       "          25.99806463,   0.        ],\n",
       "        [ 21.00505884,  24.00578153,  24.00578153, ...,  18.00433615,\n",
       "          18.98817965,   0.        ],\n",
       "        [ 30.99107042,  30.00722691,  19.99661925, ...,  38.0009554 ,\n",
       "          27.00650422,   0.        ],\n",
       "        ..., \n",
       "        [ 19.99661925,  24.00578153,  21.98890235, ...,  25.99806463,\n",
       "          31.99951001,   0.        ],\n",
       "        [ 24.98962504,  24.98962504,  35.0002327 , ...,  30.99107042,\n",
       "          16.01205305,   0.        ],\n",
       "        [  0.        ,   0.        ,   0.        , ...,   0.        ,\n",
       "           0.        ,   0.        ]],\n",
       "\n",
       "       [[ 39.00939499,  35.0002327 ,  22.99734194, ...,  24.00578153,\n",
       "          22.99734194,   0.        ],\n",
       "        [ 22.99734194,  24.98962504,  12.00289077, ...,  21.98890235,\n",
       "          18.98817965,   0.        ],\n",
       "        [ 33.00794961,  10.01060767,  27.99034773, ...,  19.99661925,\n",
       "          30.99107042,   0.        ],\n",
       "        ..., \n",
       "        [ 19.99661925,  24.98962504,  24.98962504, ...,  16.99589656,\n",
       "          21.00505884,   0.        ],\n",
       "        [ 30.00722691,  30.00722691,  13.99517386, ...,  19.99661925,\n",
       "          16.99589656,   0.        ],\n",
       "        [  0.        ,   0.        ,   0.        , ...,   0.        ,\n",
       "           0.        ,   0.        ]]]), width=600, xlim=[-100.0, 100.0], ylim=[-100.0, 100.0], zlim=[-100.0, 100.0])>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig.trait_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "## Exercise 1:\n",
    "\n",
    "Turn this image into a face\n",
    "\n",
    "- adjust the transfer function using `level`, `opacity`, and `max_opacity`\n",
    "- adjust lighting using `lighting` and `diffuse_coefficient`\n",
    "- turn of controls after you are done adjusting\n",
    "- fix the aspect ratio based on dimensions\n",
    "  The image is always rendered into a cube, so extend the volume with zeros to create a uniform cube\n",
    "- use `downscale` to do faster rendering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ely/Dropbox/SCIENCE/Software/anaconda3/anaconda/lib/python3.6/site-packages/ipyvolume/serialize.py:66: RuntimeWarning: invalid value encountered in true_divide\n",
      "  gradient = gradient / np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bc618d06fda747a8965a21e408a3c35e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = p3.figure(width=600, height=600)\n",
    "fig.style = {'axes': {'color': 'black',\n",
    "  'label': {'color': 'black'},\n",
    "  'ticklabel': {'color': 'black'},\n",
    "  'visible': False},\n",
    " 'background-color': 'white',\n",
    " 'box': {'visible': False}}\n",
    "fig.camera_fov = 1\n",
    "\n",
    "# Set viewing angle\n",
    "fig.anglex = -2.5287287502882\n",
    "fig.angley = 0.432399229463336\n",
    "fig.anglez = -1.7057381623311114\n",
    "\n",
    "# Turn volume into cube\n",
    "Vnew = np.concatenate((np.zeros((256-176, 256, 256)), V), axis=0)\n",
    "#Vnew = np.concatenate((np.zeros((img.shape[0]-176, img.shape[1], img.shape[2])), V), axis=0) #modded to adapt to image size\n",
    "#print(np.max(img.shape)-img.shape[0])\n",
    "#print(np.max(img.shape)-img.shape[1])\n",
    "#print(np.max(img.shape)-img.shape[2])\n",
    "#print(np.zeros((np.max(img.shape)-img.shape[0], np.max(img.shape), np.max(img.shape))))\n",
    "#print(np.zeros((51,1,51)))\n",
    "#Vnew = np.concatenate((np.zeros((np.max(img.shape)-img.shape[0], np.max(img.shape), np.max(img.shape)-img.shape[2])), V), axis=0) #modded to adapt to image size\n",
    "\n",
    "vol = p3.volshow(Vnew, level=[0.21, 0.30], opacity=[0.2, 1.0], max_opacity=1, lighting=True, downscale=2,\n",
    "                 diffuse_coefficient=0.5, controls=False)\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to generate a face is to use the isosurface function. This creates a surface connecting voxels left in the image after thresholding at a given value.\n",
    "\n",
    "Play with the slider to find a value that recovers most of the face."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5414446b8f6d483e92f480b60ecfb680",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = p3.figure(width=600, height=600)\n",
    "fig.style = {'axes': {'color': 'black',\n",
    "  'label': {'color': 'black'},\n",
    "  'ticklabel': {'color': 'black'},\n",
    "  'visible': False},\n",
    " 'background-color': 'white',\n",
    " 'box': {'visible': False}}\n",
    "fig.camera_fov = 1\n",
    "Vnew = np.concatenate((np.zeros((256-176, 256, 256)), V), axis=0)\n",
    "# Set viewing angle\n",
    "fig.anglex = -1.4971217646597959\n",
    "fig.angley = 0.45647865940768106\n",
    "fig.anglez = 3.0730970695793047\n",
    "vol = p3.plot_isosurface(Vnew, 200)\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "## Exercise 2:\n",
    "\n",
    "Can you use this approach to extract the brain from this image? \n",
    "\n",
    "If yes, how? If no, why not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "As can be seen with playing around with the threshold, there is no specific value that separates the skull from the brain in terms of intensities. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IPyvolume can also be used to plot meshes. In fact the isosurface output above is a triangular mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Here is an example from the IPYVolume documentation\n",
    "\n",
    "s = 1/2**0.5\n",
    "\n",
    "# Create 4 vertices for the tetrahedron\n",
    "x = np.array([1.,  -1, 0,  0])\n",
    "y = np.array([0,   0, 1., -1])\n",
    "z = np.array([-s, -s, s,  s])\n",
    "\n",
    "# and 4 surfaces (triangles or faces), where the number refer to the vertex index\n",
    "triangles = [(0, 1, 2), (0, 1, 3), (0, 2, 3), (1,3,2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9ab755da33464892b7da99247c20483a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = p3.figure()\n",
    "# we draw the tetrahedron\n",
    "p3.plot_trisurf(x, y, z, triangles=triangles, color='orange')\n",
    "# and also mark the vertices\n",
    "p3.scatter(x, y, z, marker='sphere', color='blue')\n",
    "p3.xyzlim(-2, 2)\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read the surface geometry of a subject using nibabel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#data_lh = nb.freesurfer.read_geometry(nhw+'/data/100610/T1w/Native/100610.L.pial.native.surf.gii')\n",
    "#data_rh = nb.freesurfer.read_geometry(nhw+'/data/100610/T1w/Native/100610.R.pial.native.surf.gii')\n",
    "data_lh_n = nb.load(nhw+'/data/100610/T1w/Native/100610.L.pial.native.surf.gii')\n",
    "data_lh = nb.load(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.L.midthickness_MSMAll.32k_fs_LR.surf.gii')\n",
    "data_lh_myelin = nb.load(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.L.MyelinMap.32k_fs_LR.func.gii')\n",
    "data_rh = nb.load(nhw+'/data/100610/T1w/Native/100610.R.pial.native.surf.gii')\n",
    "lharray0=data_lh.darrays[0]\n",
    "activation=data_lh_myelin.darrays[0].data\n",
    "#lharray0.data.shape\n",
    "lharray1=data_lh.darrays[1]\n",
    "#data_lh_myelin.darrays[0].shape\n",
    "(data_lh_myelin.darrays[0].data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ely/Dropbox/SCIENCE/Software/anaconda3/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:19: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "93595f341f7c4bc8b933a4a787af0149",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = p3.figure(width=600, height=600)\n",
    "fig.camera_fov = 1\n",
    "fig.style = {'axes': {'color': 'black',\n",
    "  'label': {'color': 'black'},\n",
    "  'ticklabel': {'color': 'black'},\n",
    "  'visible': False},\n",
    " 'background-color': 'white',\n",
    " 'box': {'visible': False}}\n",
    "fig.xlim = (-100, 100)\n",
    "fig.ylim = (-100, 100)\n",
    "fig.zlim = (-100, 100)\n",
    "\n",
    "# plot each hemisphere\n",
    "#x, y, z = data_lh[0][:, 0], data_lh[0][:, 1], data_lh[0][:, 2]\n",
    "#p3.plot_trisurf(x, y, z, triangles=data_lh[1], color='gray')\n",
    "\n",
    "#activation=np.linspace(0,1,x.shape[0]) #np.ones(x.shape[0]) #\n",
    "#activation2 = np.ones(x.shape[0])\n",
    "colors=plt.cm.autumn(np.log(activation))\n",
    "#colors\n",
    "x, y, z = lharray0.data.T\n",
    "surf = p3.plot_trisurf(x, y, z, triangles=lharray1.data, color=colors[:,:3])\n",
    "\n",
    "\n",
    "\n",
    "#x, y, z = data_rh[0][:, 0], data_rh[0][:, 1], data_rh[0][:, 2]\n",
    "#p3.plot_trisurf(x, y, z, triangles=data_rh[1], color='gray')\n",
    "\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9bf73f86823d4a9e97b790e5468630ed",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "color2 = plt.cm.viridis(activation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "surf.color = color2[:, :3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "## Exercise 3: \n",
    "\n",
    "Plot the smoothwm surfaces of subject 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cd24d8fef7e440bfbe05e5aee782ed54",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_lhs = nb.freesurfer.read_geometry('/data/ds000114/derivatives/freesurfer/sub-01/surf/lh.smoothwm')\n",
    "data_rhs = nb.freesurfer.read_geometry('/data/ds000114/derivatives/freesurfer/sub-01/surf/rh.smoothwm')\n",
    "\n",
    "fig = p3.figure(width=600, height=600)\n",
    "fig.camera_fov = 1\n",
    "fig.style = {'axes': {'color': 'black',\n",
    "  'label': {'color': 'black'},\n",
    "  'ticklabel': {'color': 'black'},\n",
    "  'visible': False},\n",
    " 'background-color': 'white',\n",
    " 'box': {'visible': False}}\n",
    "fig.xlim = (-100, 100)\n",
    "fig.ylim = (-100, 100)\n",
    "fig.zlim = (-100, 100)\n",
    "\n",
    "# plot each hemisphere\n",
    "x, y, z = data_lhs[0][:, 0], data_lhs[0][:, 1], data_lhs[0][:, 2]\n",
    "p3.plot_trisurf(x, y, z, triangles=data_lhs[1], color='gray')\n",
    "\n",
    "x, y, z = data_rhs[0][:, 0], data_rhs[0][:, 1], data_rhs[0][:, 2]\n",
    "p3.plot_trisurf(x, y, z, triangles=data_rhs[1], color='gray')\n",
    "\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create solution here\n",
    "\n",
    "# load data\n",
    "\n",
    "# extract and plot each hemisphere"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset we are using also has the surfaces stored in `gifti` format under `derivatives/fmriprep`. In addition, the freesurfer output for these subjects also has the parcellation files. Nibabel can read these files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data_lh = [val.data for val in nb.load('/data/ds000114/derivatives/fmriprep/sub-01/anat/sub-01_t1w_smoothwm.l.surf.gii').darrays]\n",
    "data_rh = [val.data for val in nb.load('/data/ds000114/derivatives/fmriprep/sub-01/anat/sub-01_t1w_smoothwm.r.surf.gii').darrays]\n",
    "\n",
    "# Read the annot files and extract the colors.\n",
    "annot_rh = nb.freesurfer.read_annot('/data/ds000114/derivatives/freesurfer/sub-01/label/rh.aparc.dktatlas.annot')\n",
    "colors_rh = annot_rh[1][annot_rh[0], 0:3]/255.\n",
    "annot_lh = nb.freesurfer.read_annot('/data/ds000114/derivatives/freesurfer/sub-01/label/lh.aparc.dktatlas.annot')\n",
    "colors_lh = annot_lh[1][annot_lh[0], 0:3]/255."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "268981916475422498713befc8649924",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = p3.figure(width=600, height=600)\n",
    "# we draw the tetrahedron\n",
    "x, y, z = data_lh[0][:, 0], data_lh[0][:, 1], data_lh[0][:, 2]\n",
    "p3.plot_trisurf(x, y, z, triangles=data_lh[1], color=colors_lh)\n",
    "x, y, z = data_rh[0][:, 0], data_rh[0][:, 1], data_rh[0][:, 2]\n",
    "p3.plot_trisurf(x, y, z, triangles=data_rh[1], color=colors_rh)\n",
    "fig.camera_fov = 1\n",
    "fig.xlim = (-100, 100)\n",
    "fig.ylim = (-100, 100)\n",
    "fig.zlim = (-100, 100)\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can load the functional time series preprocessed by fmriprep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lh_func = nb.load('/data/ds000114/derivatives/fmriprep/sub-02/ses-test/func/sub-02_ses-test_task-fingerfootlips_bold_space-fsaverage5.l.func.gii')\n",
    "data = np.vstack([val.data for val in lh_func.darrays]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGbRJREFUeJzt3X+w3XWd3/Hnq0HwtwS5WppgE9fsD2SsYtS0tjtWVgi6\nY+gMzITZltRlJlMXt26nO2uoM8tWZQa322WXGWSHlSzBOgbKuiVTYrMZwDqdkR9BFIiIuQKVK6yJ\nG2TdOsLGffeP87l6vJzc+805ufeckOdj5sz5nvf38/2e9/lyuK98f5xzUlVIktTFPxh3A5KkY4eh\nIUnqzNCQJHVmaEiSOjM0JEmdGRqSpM4MDUlSZ4aGJKkzQ0OS1NkJ427gaDv11FNr1apV425Dko4p\n99133/eqamqhcS+40Fi1ahV79uwZdxuSdExJ8n+7jPPwlCSpM0NDktSZoSFJ6szQkCR1ZmhIkjoz\nNCRJnRkakqTODA1JUmeGhiSpswVDI8nWJPuTPDSn/ptJHkmyN8nv99UvSzLd5p3bV1/fatNJtvTV\nVye5O8m+JDclObHVT2qPp9v8VUfjBc9n1ZbbfnKTJD1flz2NG4D1/YUk/xLYALypqt4I/EGrnwFs\nBN7YlvlUkmVJlgHXAOcBZwAXtbEAnwSuqqo1wNPAJa1+CfB0Vb0BuKqNkySN0YKhUVVfAg7OKX8Q\nuLKqnm1j9rf6BmB7VT1bVY8B08Db2226qh6tqueA7cCGJAHeDdzSlt8GnN+3rm1t+hbg7DZekjQm\nw57T+HngX7TDRv87ydtafQXwRN+4mVY7XP3VwPer6tCc+s+sq81/po2XJI3JsN9yewKwHFgHvA24\nOcnrgUF7AsXgcKp5xrPAvJ+RZDOwGeB1r3vdvI1LkoY37J7GDPD56rkH+Hvg1FY/vW/cSuDJeerf\nA05OcsKcOv3LtPmv4vmHyQCoquuqam1VrZ2aWvDr4CVJQxo2NP4HvXMRJPl54ER6AbAD2NiufFoN\nrAHuAe4F1rQrpU6kd7J8R1UVcCdwQVvvJuDWNr2jPabNv6ONlySNyYKHp5J8DngXcGqSGeByYCuw\ntV2G+xywqf1B35vkZuDrwCHg0qr6cVvPh4BdwDJga1XtbU/xEWB7kk8A9wPXt/r1wGeSTNPbw9h4\nFF6vJGkEC4ZGVV10mFn/+jDjrwCuGFDfCewcUH+U3tVVc+s/Ai5cqD9J0tLxE+GSpM4MDUlSZ4aG\nJKkzQ0OS1JmhIUnqzNCQJHVmaEiSOjM0JEmdGRqSpM4MDUlSZ4aGJKkzQ0OS1JmhIUnqzNCQJHVm\naEiSOjM0JEmdGRqSpM4WDI0kW5Psbz/tOnfebyepJKe2x0lydZLpJA8kOatv7KYk+9ptU1/9rUke\nbMtcnSStfkqS3W387iTLj85LliQNq8uexg3A+rnFJKcD7wG+3Vc+D1jTbpuBa9vYU+j9tvg76P20\n6+V9IXBtGzu73OxzbQFur6o1wO3tsSRpjBYMjar6EnBwwKyrgN8Bqq+2Abixeu4CTk5yGnAusLuq\nDlbV08BuYH2b98qq+nJVFXAjcH7fura16W19dUnSmAx1TiPJ+4HvVNXX5sxaATzR93im1earzwyo\nA7y2qp4CaPevmaefzUn2JNlz4MCBIV6RJKmLIw6NJC8FPgr87qDZA2o1RP2IVNV1VbW2qtZOTU0d\n6eKSpI6G2dP4OWA18LUkjwMrga8k+Yf09hRO7xu7EnhygfrKAXWA77bDV7T7/UP0Kkk6io44NKrq\nwap6TVWtqqpV9P7wn1VVfwXsAC5uV1GtA55ph5Z2AeckWd5OgJ8D7GrzfpBkXbtq6mLg1vZUO4DZ\nq6w29dUlSWPS5ZLbzwFfBn4hyUySS+YZvhN4FJgG/hT4DYCqOgh8HLi33T7WagAfBD7dlvkW8IVW\nvxJ4T5J99K7SuvLIXpok6Wg7YaEBVXXRAvNX9U0XcOlhxm0Ftg6o7wHOHFD/a+DshfqTJC0dPxEu\nSerM0JAkdWZoSJI6MzQkSZ0ZGpKkzgwNSVJnhoYkqTNDQ5LUmaEhSerM0JAkdWZoSJI6MzQkSZ0Z\nGpKkzgwNSVJnhoYkqTNDQ5LUWZdf7tuaZH+Sh/pq/yXJN5I8kOQvkpzcN++yJNNJHklybl99fatN\nJ9nSV1+d5O4k+5LclOTEVj+pPZ5u81cdrRctSRpOlz2NG4D1c2q7gTOr6k3AN4HLAJKcAWwE3tiW\n+VSSZUmWAdcA5wFnABe1sQCfBK6qqjXA08Dsz8leAjxdVW8ArmrjJEljtGBoVNWXgINzan9ZVYfa\nw7uAlW16A7C9qp6tqsfo/e7329ttuqoerarngO3AhiQB3g3c0pbfBpzft65tbfoW4Ow2XpI0Jkfj\nnMavA19o0yuAJ/rmzbTa4eqvBr7fF0Cz9Z9ZV5v/TBsvSRqTkUIjyUeBQ8BnZ0sDhtUQ9fnWNaiP\nzUn2JNlz4MCB+ZuWJA1t6NBIsgn4VeDXqmr2j/kMcHrfsJXAk/PUvwecnOSEOfWfWVeb/yrmHCab\nVVXXVdXaqlo7NTU17EuSJC1gqNBIsh74CPD+qvph36wdwMZ25dNqYA1wD3AvsKZdKXUivZPlO1rY\n3Alc0JbfBNzat65NbfoC4I6+cJIkjcEJCw1I8jngXcCpSWaAy+ldLXUSsLudm76rqv5dVe1NcjPw\ndXqHrS6tqh+39XwI2AUsA7ZW1d72FB8Btif5BHA/cH2rXw98Jsk0vT2MjUfh9UqSRrBgaFTVRQPK\n1w+ozY6/ArhiQH0nsHNA/VF6V1fNrf8IuHCh/iRJS8dPhEuSOjM0JEmdGRqSpM4MDUlSZ4aGJKkz\nQ0OS1JmhIUnqzNCQJHVmaEiSOjM0JEmdGRqSpM4MDUlSZ4aGJKkzQ0OS1JmhIUnqzNCQJHW2YGgk\n2Zpkf5KH+mqnJNmdZF+7X97qSXJ1kukkDyQ5q2+ZTW38vvb74rP1tyZ5sC1zddpPAR7uOSRJ49Nl\nT+MGYP2c2hbg9qpaA9zeHgOcR+93wdcAm4FroRcA9H4m9h30fqXv8r4QuLaNnV1u/QLPIUkakwVD\no6q+RO83uvttALa16W3A+X31G6vnLuDkJKcB5wK7q+pgVT0N7AbWt3mvrKovV1UBN85Z16DnkCSN\nybDnNF5bVU8BtPvXtPoK4Im+cTOtNl99ZkB9vueQJI3J0T4RngG1GqJ+ZE+abE6yJ8meAwcOHOni\nkqSOhg2N77ZDS7T7/a0+A5zeN24l8OQC9ZUD6vM9x/NU1XVVtbaq1k5NTQ35kiRJCxk2NHYAs1dA\nbQJu7atf3K6iWgc80w4t7QLOSbK8nQA/B9jV5v0gybp21dTFc9Y16DkkSWNywkIDknwOeBdwapIZ\neldBXQncnOQS4NvAhW34TuC9wDTwQ+ADAFV1MMnHgXvbuI9V1ezJ9Q/Su0LrJcAX2o15nkNadKu2\n3AbA41e+b8ydSJNlwdCoqosOM+vsAWMLuPQw69kKbB1Q3wOcOaD+14OeQ5I0Pn4iXOpo1ZbbfrIH\nIh2vDA1JUmcLHp6SjmfuWUg/yz0N6Qh5mErHM0NDktSZoSFJ6szQkCR15olwqY/nKqT5uachSerM\n0JAkdWZoSJI685yGNKT+8x9+saGOF+5pSJI6MzQkSZ0ZGtJR4FeL6HhhaEiSOjM0JEmdjRQaSf5D\nkr1JHkryuSQvTrI6yd1J9iW5KcmJbexJ7fF0m7+qbz2XtfojSc7tq69vtekkW0bpVZI0uqFDI8kK\n4N8Da6vqTGAZsBH4JHBVVa0BngYuaYtcAjxdVW8ArmrjSHJGW+6NwHrgU0mWJVkGXAOcB5wBXNTG\nSpLGZNTDUycAL0lyAvBS4Cng3cAtbf424Pw2vaE9ps0/O0lafXtVPVtVjwHTwNvbbbqqHq2q54Dt\nbawkaUyGDo2q+g7wB8C36YXFM8B9wPer6lAbNgOsaNMrgCfasofa+Ff31+csc7j68yTZnGRPkj0H\nDhwY9iVJkhYwyuGp5fT+5b8a+EfAy+gdSpqrZhc5zLwjrT+/WHVdVa2tqrVTU1MLtS5JGtIoh6d+\nBXisqg5U1d8Bnwf+GXByO1wFsBJ4sk3PAKcDtPmvAg721+csc7i6JGlMRgmNbwPrkry0nZs4G/g6\ncCdwQRuzCbi1Te9oj2nz76iqavWN7eqq1cAa4B7gXmBNuxrrRHony3eM0K8kaURDf2FhVd2d5Bbg\nK8Ah4H7gOuA2YHuST7Ta9W2R64HPJJmmt4exsa1nb5Kb6QXOIeDSqvoxQJIPAbvoXZm1tar2Dtuv\nJGl0I33LbVVdDlw+p/wovSuf5o79EXDhYdZzBXDFgPpOYOcoPUqSjh6/Gl06ivy6dL3Q+TUikqTO\nDA1JUmeGhiSpM0NDktSZoSFJ6szQkCR1ZmhIkjrzcxo67vnb3lJ37mlIkjozNCRJnRka0iJZteU2\nD33pBcfQkCR1ZmhIkjozNCRJnRkakqTORgqNJCcnuSXJN5I8nOSfJjklye4k+9r98jY2Sa5OMp3k\ngSRn9a1nUxu/L8mmvvpbkzzYlrm6/aysJGlMRt3T+GPgf1XVLwL/BHgY2ALcXlVrgNvbY4Dz6P3+\n9xpgM3AtQJJT6P363zvo/eLf5bNB08Zs7ltu/Yj9SpJGMPQnwpO8Evhl4N8CVNVzwHNJNgDvasO2\nAV8EPgJsAG6sqgLuanspp7Wxu6vqYFvvbmB9ki8Cr6yqL7f6jcD5wBeG7Vnq5+Ww0pEbZU/j9cAB\n4M+S3J/k00leBry2qp4CaPevaeNXAE/0LT/TavPVZwbUJUljMkponACcBVxbVW8B/h8/PRQ1yKDz\nETVE/fkrTjYn2ZNkz4EDB+bvWpI0tFFCYwaYqaq72+Nb6IXId9thJ9r9/r7xp/ctvxJ4coH6ygH1\n56mq66pqbVWtnZqaGuElSZLmM3RoVNVfAU8k+YVWOhv4OrADmL0CahNwa5veAVzcrqJaBzzTDl/t\nAs5JsrydAD8H2NXm/SDJunbV1MV965KOGX6diF5IRv1q9N8EPpvkROBR4AP0gujmJJcA3wYubGN3\nAu8FpoEftrFU1cEkHwfubeM+NntSHPggcAPwEnonwD0JLkljNFJoVNVXgbUDZp09YGwBlx5mPVuB\nrQPqe4AzR+lR6ue/+KXR+IlwSVJnhoYkqTNDQ5LUmaEhSerM0JAkdWZoSJI6MzQkSZ0ZGpKkzgwN\nSVJno36NiHRMmIRPgvf38PiV7xtjJ9Lw3NOQJHVmaEiSOjM0JEmdGRqSpM4MDUlSZ4aGNAb+mp+O\nVSOHRpJlSe5P8j/b49VJ7k6yL8lN7Vf9SHJSezzd5q/qW8dlrf5IknP76utbbTrJllF7lSSN5mjs\naXwYeLjv8SeBq6pqDfA0cEmrXwI8XVVvAK5q40hyBrAReCOwHvhUC6JlwDXAecAZwEVtrCRpTEYK\njSQrgfcBn26PA7wbuKUN2Qac36Y3tMe0+We38RuA7VX1bFU9Ru83xN/ebtNV9WhVPQdsb2MlSWMy\n6ifC/wj4HeAV7fGrge9X1aH2eAZY0aZXAE8AVNWhJM+08SuAu/rW2b/ME3Pq7xixXx1nPG8gHV1D\n72kk+VVgf1Xd118eMLQWmHek9UG9bE6yJ8meAwcOzNO1JGkUoxyeeifw/iSP0zt09G56ex4nJ5nd\ng1kJPNmmZ4DTAdr8VwEH++tzljlc/Xmq6rqqWltVa6empkZ4SZKk+QwdGlV1WVWtrKpV9E5k31FV\nvwbcCVzQhm0Cbm3TO9pj2vw7qqpafWO7umo1sAa4B7gXWNOuxjqxPceOYfuVJpGX3upYsxjfcvsR\nYHuSTwD3A9e3+vXAZ5JM09vD2AhQVXuT3Ax8HTgEXFpVPwZI8iFgF7AM2FpVexehX0lSR0clNKrq\ni8AX2/Sj9K58mjvmR8CFh1n+CuCKAfWdwM6j0aMkaXR+IlyS1JmhIUnqzNCQJHVmaEiSOjM0JEmd\nGRqSpM4MDUlSZ4aGJKkzQ0OS1JmhIUnqzNCQJHW2GF9YKI2V3xorLR5DQ5oA/UH3+JXvG2Mn0vw8\nPCVJ6szQkCR1ZmhIkjozNCRJnQ0dGklOT3JnkoeT7E3y4VY/JcnuJPva/fJWT5Krk0wneSDJWX3r\n2tTG70uyqa/+1iQPtmWuTpJRXqwkaTSj7GkcAv5jVf0SsA64NMkZwBbg9qpaA9zeHgOcB6xpt83A\ntdALGeBy4B30fib28tmgaWM29y23foR+pWPCqi23edmwJtbQoVFVT1XVV9r0D4CHgRXABmBbG7YN\nOL9NbwBurJ67gJOTnAacC+yuqoNV9TSwG1jf5r2yqr5cVQXc2LcuSdIYHJXPaSRZBbwFuBt4bVU9\nBb1gSfKaNmwF8ETfYjOtNl99ZkB90PNvprdHwute97rRXoyOWf7rXFp8I58IT/Jy4M+B36qqv5lv\n6IBaDVF/frHquqpaW1Vrp6amFmpZkjSkkUIjyYvoBcZnq+rzrfzddmiJdr+/1WeA0/sWXwk8uUB9\n5YC6JGlMRrl6KsD1wMNV9Yd9s3YAs1dAbQJu7atf3K6iWgc80w5j7QLOSbK8nQA/B9jV5v0gybr2\nXBf3rUuSNAajnNN4J/BvgAeTfLXV/hNwJXBzkkuAbwMXtnk7gfcC08APgQ8AVNXBJB8H7m3jPlZV\nB9v0B4EbgJcAX2g3SdKYDB0aVfV/GHzeAeDsAeMLuPQw69oKbB1Q3wOcOWyP0rHMLzHUJPIT4ZKk\nzgyNBfhBK0n6KX9PQ8c8Q11aOobGYfiHSJKez8NTkqTODA3pGOC5NU0KD0/pmOQfUGk83NOQJHXm\nnoYmnh9y+6nZbXG8bweNj6GhiTXoEJSHpaTxMjSOMv8leGQMAenYYmh0NMofty5B8kILGw8pLS63\nr8Ylve8RfOFYu3Zt7dmzZ6hlJ+1fvbN/DBbrD8SRBtXc8ZO2vY5XhoaOhiT3VdXahca5pzHBFuOP\n8nzr7Pp8hsVkeaHtpWqyGRrHiaMRFppsHrLSUjA0jjEeIpI0ThMfGknWA38MLAM+XVVXjrmliWBY\naD4estJimejQSLIMuAZ4DzAD3JtkR1V9fbydSccWD13paJno0ADeDkxX1aMASbYDGwBDQ+pgvg9I\nGh4axqSHxgrgib7HM8A7xtSL9ILS9RCn4aJ+kx4aGVB73gdLkmwGNreHf5vkkSGf71Tge0Muu9Ts\ndXHY6xz55FFZjdt18Rytfv9xl0GTHhozwOl9j1cCT84dVFXXAdeN+mRJ9nT5cMsksNfFYa+Lw14X\nz1L3O+lfjX4vsCbJ6iQnAhuBHWPuSZKOWxO9p1FVh5J8CNhF75LbrVW1d8xtSdJxa6JDA6CqdgI7\nl+jpRj7EtYTsdXHY6+Kw18WzpP2+4L6wUJK0eCb9nIYkaYIYGk2S9UkeSTKdZMu4+wFI8niSB5N8\nNcmeVjslye4k+9r98lZPkqtb/w8kOWuRe9uaZH+Sh/pqR9xbkk1t/L4km5aw199L8p22bb+a5L19\n8y5rvT6S5Ny++qK/R5KcnuTOJA8n2Zvkw60+cdt2nl4nbtsmeXGSe5J8rfX6n1t9dZK72za6qV1w\nQ5KT2uPpNn/VQq9hCXq9Icljfdv1za2+tO+Bqjrub/ROsn8LeD1wIvA14IwJ6Otx4NQ5td8HtrTp\nLcAn2/R7gS/Q+2zLOuDuRe7tl4GzgIeG7Q04BXi03S9v08uXqNffA357wNgz2n//k4DV7X2xbKne\nI8BpwFlt+hXAN1tPE7dt5+l14rZt2z4vb9MvAu5u2+tmYGOr/wnwwTb9G8CftOmNwE3zvYYl6vUG\n4IIB45f0PeCeRs9Pvq6kqp4DZr+uZBJtALa16W3A+X31G6vnLuDkJKctVhNV9SXg4Ii9nQvsrqqD\nVfU0sBtYv0S9Hs4GYHtVPVtVjwHT9N4fS/IeqaqnquorbfoHwMP0vhlh4rbtPL0ezti2bds+f9se\nvqjdCng3cEurz92us9v7FuDsJJnnNSxFr4ezpO8BQ6Nn0NeVzPfmXyoF/GWS+9L71DvAa6vqKej9\nTwu8ptUn4TUcaW/j7vlDbXd+6+zhnnl6WvJe2yGRt9D7l+ZEb9s5vcIEbtsky5J8FdhP7w/ot4Dv\nV9WhAc/7k57a/GeAV4+r16qa3a5XtO16VZKT5vY6p6dF6dXQ6On0dSVj8M6qOgs4D7g0yS/PM3ZS\nXwMcvrdx9nwt8HPAm4GngP/a6hPRa5KXA38O/FZV/c18QwfUlrTfAb1O5Latqh9X1ZvpfbPE24Ff\nmud5J6rXJGcClwG/CLyN3iGnj4yjV0Ojp9PXlSy1qnqy3e8H/oLeG/27s4ed2v3+NnwSXsOR9ja2\nnqvqu+1/zL8H/pSfHmIYe69JXkTvj/Bnq+rzrTyR23ZQr5O8bVt/3we+SO/4/8lJZj+v1v+8P+mp\nzX8VvUOc4+p1fTscWFX1LPBnjGm7Gho9E/d1JUleluQVs9PAOcBDra/ZqyA2Abe26R3Axe1KinXA\nM7OHM5bQkfa2CzgnyfJ2COOcVlt0c873/Ct623a2143t6pnVwBrgHpboPdKOm18PPFxVf9g3a+K2\n7eF6ncRtm2Qqyclt+iXAr9A7B3MncEEbNne7zm7vC4A7qnd2+XCvYbF7/UbfPxpC79xL/3ZduvfA\nqGfSXyg3elcgfJPecc6PTkA/r6d3lcbXgL2zPdE7rno7sK/dn1I/veLimtb/g8DaRe7vc/QOPfwd\nvX/RXDJMb8Cv0zuZOA18YAl7/Uzr5YH2P91pfeM/2np9BDhvKd8jwD+ndwjhAeCr7fbeSdy28/Q6\ncdsWeBNwf+vpIeB3+/4/u6dto/8OnNTqL26Pp9v81y/0Gpag1zvadn0I+G/89AqrJX0P+IlwSVJn\nHp6SJHVmaEiSOjM0JEmdGRqSpM4MDUlSZ4aGJKkzQ0OS1JmhIUnq7P8DQ+6y84NRXvEAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcddec65518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(data.flatten(), 128);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10242,) 1542.41 288.133\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/neuro3/lib/python3.6/site-packages/numpy/lib/function_base.py:4011: RuntimeWarning: All-NaN slice encountered\n",
      "  r = func(a, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "data[data<500] = np.nan\n",
    "mean_val = np.nanmean(data)\n",
    "data_change = (data - mean_val)/mean_val\n",
    "mean_epi = np.nanmedian(data, axis=1) - mean_val\n",
    "mean_std = np.nanstd(mean_epi)\n",
    "mean_epi[np.isnan(mean_epi)] = 0.\n",
    "data[np.isnan(data)] = 0.\n",
    "print(mean_epi.shape, mean_val, mean_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFfNJREFUeJzt3X+M5Hd93/Hnq77YLSTgX2vq3B3dc3LQOFFanI1xSoko\n5oftUI62uLKVxifi6pTUpFAahXMtxVEiJEiakKBQowu+YkfExiVQn4JTczEmqFJsvCbGPzDgxTh4\n8eFbamNoaUwd3v1jPhvGe3P7Y2ZnZ+++z4e02u+8v5+Zec93Z+c13x8z31QVkqTu+TuTbkCSNBkG\ngCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUVsm3cByTj/99Jqenp50G5J0TLn7\n7ru/XlVTK43b1AEwPT3N7OzspNuQpGNKkr9azTg3AUlSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLU\nUQaAJHWUASBJHWUASFJHGQDSJje992NM7/3YpNvQccgAkKSOMgAkqaMMAEnqKANAkjrKAJCkjjIA\nJKmjVgyAJPuTHE5y/5L6LyX5QpIHkvxmX/3KJHNt3mv76he02lySvev7MCRJa7WaM4J9APh94PrF\nQpJ/BuwCfryqnk5yRqufDVwC/Cjwg8CfJXlRu9p7gVcD88BdSQ5U1efW64FIktZmxQCoqk8lmV5S\n/kXgnVX1dBtzuNV3ATe2+peTzAHntnlzVfUwQJIb21gDQJImZNh9AC8CXp7kziR/nuQnW30r8Gjf\nuPlWO1pdkjQhw54UfgtwCnAe8JPATUnOAjJgbDE4aGrQDSfZA+wBeOELXzhke5KklQy7BjAPfKR6\nPg18Fzi91bf3jdsGPLZM/QhVta+qZqpqZmpqasj2JEkrGTYA/jvwSoC2k/dE4OvAAeCSJCcl2QHs\nBD4N3AXsTLIjyYn0dhQfGLV5SdLwVtwElOQG4BXA6UnmgauB/cD+dmjod4DdVVXAA0luordz9xng\niqr6m3Y7bwZuBU4A9lfVA2N4PJKkVVrNUUCXHmXWvznK+HcA7xhQvwW4ZU3dSZLGxk8CS1JHGQCS\n1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS\n1FEGgCR1lAEgSR21YgAk2Z/kcDv719J5v5ykkpzeLifJe5LMJbk3yTl9Y3cneaj97F7fhyFJWqvV\nrAF8ALhgaTHJduDVwFf6yhfSOw/wTmAPcE0beyq9U0m+FDgXuDrJKaM0LkkazYoBUFWfAp4YMOvd\nwK8A1VfbBVxfPXcAJyc5E3gtcLCqnqiqJ4GDDAgVSdLGGWofQJLXA1+tqs8umbUVeLTv8nyrHa0+\n6Lb3JJlNMruwsDBMe5KkVVhzACR5DnAV8KuDZg+o1TL1I4tV+6pqpqpmpqam1tqeJGmVhlkD+CFg\nB/DZJI8A24DPJPn79N7Zb+8buw14bJm6JGlC1hwAVXVfVZ1RVdNVNU3vxf2cqvoacAC4rB0NdB7w\nVFUdAm4FXpPklLbz9zWtJkmakNUcBnoD8BfAi5PMJ7l8meG3AA8Dc8AfAP8OoKqeAH4DuKv9/Hqr\nSZImZMtKA6rq0hXmT/dNF3DFUcbtB/avsT9J0pj4SWBJ6igDQJI6ygCQpI4yACSpowwASeooA0CS\nOsoAkKSOMgAkqaMMAEnqKANAkjrKAJCkjjIAJKmjDABJ6igDQJI6ygCQpI5azQlh9ic5nOT+vtpv\nJfl8knuTfDTJyX3zrkwyl+QLSV7bV7+g1eaS7F3/hyJJWovVrAF8ALhgSe0g8GNV9ePAF4ErAZKc\nDVwC/Gi7zn9JckKSE4D3AhcCZwOXtrGSpAlZMQCq6lPAE0tqH6+qZ9rFO+id5B1gF3BjVT1dVV+m\nd2rIc9vPXFU9XFXfAW5sYyVJE7Ie+wB+HvjTNr0VeLRv3nyrHa0uSZqQkQIgyVXAM8AHF0sDhtUy\n9UG3uSfJbJLZhYWFUdqTJC1j6ABIsht4HfCz7WTw0Htnv71v2DbgsWXqR6iqfVU1U1UzU1NTw7Yn\nSVrBUAGQ5ALg7cDrq+rbfbMOAJckOSnJDmAn8GngLmBnkh1JTqS3o/jAaK1LkkaxZaUBSW4AXgGc\nnmQeuJreUT8nAQeTANxRVb9QVQ8kuQn4HL1NQ1dU1d+023kzcCtwArC/qh4Yw+ORJK3SigFQVZcO\nKF+7zPh3AO8YUL8FuGVN3UmSxsZPAktSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJ\nHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdtWIAJNmf5HCS+/tqpyY5mOSh\n9vuUVk+S9ySZS3JvknP6rrO7jX+onU9YkjRBq1kD+ABwwZLaXuC2qtoJ3NYuA1xI7zzAO4E9wDXQ\nCwx6p5J8KXAucPViaEiSJmPFAKiqTwFPLCnvAq5r09cBb+irX189dwAnJzkTeC1wsKqeqKongYMc\nGSqSpA007D6AF1TVIYD2+4xW3wo82jduvtWOVj9Ckj1JZpPMLiwsDNmeJGkl670TOANqtUz9yGLV\nvqqaqaqZqampdW1OkvQ9wwbA423TDu334VafB7b3jdsGPLZMXZI0IcMGwAFg8Uie3cDNffXL2tFA\n5wFPtU1EtwKvSXJK2/n7mlaTJE3IlpUGJLkBeAVwepJ5ekfzvBO4KcnlwFeAi9vwW4CLgDng28Cb\nAKrqiSS/AdzVxv16VS3dsSxJ2kArBkBVXXqUWecPGFvAFUe5nf3A/jV1J0kaGz8JLEkdZQBIUkcZ\nAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZ\nAJLUUSMFQJL/kOSBJPcnuSHJ302yI8mdSR5K8qEkJ7axJ7XLc23+9Ho8AEnScIYOgCRbgX8PzFTV\njwEnAJcA7wLeXVU7gSeBy9tVLgeerKofBt7dxkmSJmTUTUBbgL+XZAvwHOAQ8Ergw23+dcAb2vSu\ndpk2//wkGfH+JUlDGjoAquqrwH+md07gQ8BTwN3AN6rqmTZsHtjaprcCj7brPtPGnzbs/UuSRjPK\nJqBT6L2r3wH8IPBc4MIBQ2vxKsvM67/dPUlmk8wuLCwM254kaQWjbAJ6FfDlqlqoqv8HfAT4J8DJ\nbZMQwDbgsTY9D2wHaPOfDzyx9Earal9VzVTVzNTU1AjtSZKWM0oAfAU4L8lz2rb884HPAbcDb2xj\ndgM3t+kD7TJt/ieq6og1AEnSxhhlH8Cd9Hbmfga4r93WPuDtwNuSzNHbxn9tu8q1wGmt/jZg7wh9\nS5JGtGXlIUdXVVcDVy8pPwycO2DsXwMXj3J/kqT14yeBJamjDABJ6igDQJI6ygCQpI4yACSpowwA\nSeooA0CSOsoAkKSOMgAkqaMMAEnqKANAkjrKAJCkjjIAJKmjDABJ6igDQJI6ygCQpI4aKQCSnJzk\nw0k+n+TBJD+V5NQkB5M81H6f0sYmyXuSzCW5N8k56/MQJEnDGHUN4PeA/1FV/xD4R8CD9E71eFtV\n7QRu43unfrwQ2Nl+9gDXjHjfkqQRDB0ASZ4H/DTtnL9V9Z2q+gawC7iuDbsOeEOb3gVcXz13ACcn\nOXPoziVJIxllDeAsYAH4r0n+Msn7kzwXeEFVHQJov89o47cCj/Zdf77VniXJniSzSWYXFhZGaE+S\ntJxRAmALcA5wTVW9BPg/fG9zzyAZUKsjClX7qmqmqmampqZGaE+StJxRAmAemK+qO9vlD9MLhMcX\nN+2034f7xm/vu/424LER7l+SNIKhA6CqvgY8muTFrXQ+8DngALC71XYDN7fpA8Bl7Wig84CnFjcV\nSZI23pYRr/9LwAeTnAg8DLyJXqjclORy4CvAxW3sLcBFwBzw7TZWkjQhIwVAVd0DzAyYdf6AsQVc\nMcr9SZLWj58ElqSOMgAkqaMMAEnqKANAkjrKAJCkjjIAJKmjDABJ6igDQJI6ygCQpI4yACSpowwA\nSeooA0CSOsoAkKSOMgAkqaMMAEnqqJEDIMkJ7aTwf9Iu70hyZ5KHknyonSyGJCe1y3Nt/vSo9y1J\nGt56rAG8BXiw7/K7gHdX1U7gSeDyVr8ceLKqfhh4dxsnSZqQkQIgyTbgZ4D3t8sBXknvBPEA1wFv\naNO72mXa/PPbeEnSBIy6BvC7wK8A322XTwO+UVXPtMvzwNY2vRV4FKDNf6qNlyRNwNABkOR1wOGq\nuru/PGBorWJe/+3uSTKbZHZhYWHY9iRJKxhlDeBlwOuTPALcSG/Tz+8CJydZPNn8NuCxNj0PbAdo\n858PPLH0RqtqX1XNVNXM1NTUCO1JkpYzdABU1ZVVta2qpoFLgE9U1c8CtwNvbMN2Aze36QPtMm3+\nJ6rqiDUASdLGGMfnAN4OvC3JHL1t/Ne2+rXAaa3+NmDvGO5bkrRKW1YesrKq+iTwyTb9MHDugDF/\nDVy8HvcnSRqdnwSWpI4yACSpowwASeooA0CSOsoAkKSOMgAkqaMMAEnqKANAkjrKAJCkjlqXTwJL\nWn/Tez826RZ0nDMApONUf4A88s6fmWAn2qwMAGkTWeu7fl/kNQoDQDoGDQqKxZpBoNUyAKTjzHJr\nEYaE+hkA0iawmk0/7hTWejMApAmZ5Au6+w4Eo50UfnuS25M8mOSBJG9p9VOTHEzyUPt9SqsnyXuS\nzCW5N8k56/UgJElrN8oHwZ4B/mNV/QhwHnBFkrPpnerxtqraCdzG9079eCGws/3sAa4Z4b4lSSMa\nehNQVR0CDrXpbyV5ENgK7AJe0YZdR+9UkW9v9evbieDvSHJykjPb7Ugao7XuGHZncTesy1dBJJkG\nXgLcCbxg8UW9/T6jDdsKPNp3tflWW3pbe5LMJpldWFhYj/YkDWl678fc+XwcG3kncJLvB/4YeGtV\nfTPJUYcOqNURhap9wD6AmZmZI+ZLGg9f6LtnpABI8n30Xvw/WFUfaeXHFzftJDkTONzq88D2vqtv\nAx4b5f6lzc6jbbSZjXIUUIBrgQer6nf6Zh0Adrfp3cDNffXL2tFA5wFPuf1fkiZnlDWAlwE/B9yX\n5J5W+0/AO4GbklwOfAW4uM27BbgImAO+DbxphPuWjlluatFmMcpRQP+Twdv1Ac4fML6AK4a9P+lY\n4ou8jgV+ElhaR8frC/+gx+U+jWOfZwSTNBQPET32uQYgDWmtH5byxVKbjQEgjajrL+yDDnX18Ndj\ngwEgrUHXX+xX4vI5thgA0ir4wja8pcvONYLNw53AktRRrgFI2lDuH9g8DABpGW760fHMTUCS1FGu\nAUhL+K5fXWEASI0v/BtvuWXu/oHxMwDUeb7wb27uNB6f4zoAPK+pwBeQY5XBPH7HdQBIR+OLizSB\nAEhyAfB7wAnA+6vqneO+T9cENpdhv0RtPcb7wn9sW+5Txf6fr11652nZoDtLTgC+CLya3jmC7wIu\nrarPDRo/MzNTs7OzQ9/fsP/s43gCjWszxGqe9Mvd91q/532lL/5a7jZW8/dY7jaltehyOCS5u6pm\nVhq30WsA5wJzVfUwQJIbgV3AwACYlEFPltW88xj2ftbTam9zPcYNO2/Y+5PWYtBzyaOOnm2j1wDe\nCFxQVf+2Xf454KVV9eZB4ye1BiBJR7OatdRJr31s1jWAQecQflYCJdkD7GkX/3eSL4ypl9OBr4/p\nttebvY6HvY7Hcd1r3jXcmNVcbwVr6fUfrGbQRgfAPLC97/I24LH+AVW1D9g37kaSzK4mITcDex0P\nex0Pex2PcfS60d8FdBewM8mOJCcClwAHNrgHSRIbvAZQVc8keTNwK73DQPdX1QMb2YMkqWfDPwdQ\nVbcAt2z0/Q4w9s1M68hex8Nex8Nex2Pde93Qo4AkSZuH5wOQpI46LgMgycVJHkjy3SQzS+ZdmWQu\nyReSvLavfkGrzSXZ21ffkeTOJA8l+VDbeT2uvj+U5J7280iSe1p9Osn/7Zv3vr7r/ESS+1rf70ky\n6FDbcfT6a0m+2tfTRX3z1rSMN6DX30ry+ST3JvlokpNbfdMt1wG9T2SZLdPP9iS3J3mw/Y+9pdXX\n/HzYoH4faX/He5LMttqpSQ62/+mDSU5p9bS/9Vx7rpyzgX2+uG/Z3ZPkm0neOvblWlXH3Q/wI8CL\ngU8CM331s4HPAicBO4Av0dsZfUKbPgs4sY05u13nJuCSNv0+4Bc36DH8NvCrbXoauP8o4z4N/BS9\nz1j8KXDhBvX3a8AvD6iveRlvQK+vAba06XcB79qsy3VJDxNbZsv0dCZwTpv+AXpf7XL2Wp8PG9jv\nI8DpS2q/Cext03v7ng8Xtb91gPOAOye0jE8AvkbvWP6xLtfjcg2gqh6sqkEfINsF3FhVT1fVl4E5\nel9P8bdfUVFV3wFuBHa1d32vBD7crn8d8IZx99/u918DN6ww7kzgeVX1F9V7Vly/Ef2tYE3LeCMa\nqqqPV9Uz7eId9D5/clSbaLlObJkdTVUdqqrPtOlvAQ8CW5e5ytGeD5O0i97/Mjz7f3oXcH313AGc\n3J4LG+184EtV9VfLjFmX5XpcBsAytgKP9l2eb7Wj1U8DvtH34rFYH7eXA49X1UN9tR1J/jLJnyd5\neattbT0t2qj+Fr25rSrvX1yNZu3LeKP9PL13eYs243JdtFmW2UBJpoGXAHe20lqeDxulgI8nuTu9\nbxkAeEFVHYJeoAFntPqke110Cc9+8ze25XrMBkCSP0ty/4Cf5d4hHe2rKNZaH9oq+76UZz8BDgEv\nrKqXAG8D/ijJ88bR3xp6vQb4IeAft/5+e/FqR+lpkr0ujrkKeAb4YCtNZLmuwWbp4whJvh/4Y+Ct\nVfVN1v582Cgvq6pzgAuBK5L89DJjJ90r6e1jfD3w31pprMv1mD0hTFW9aoirLfdVFIPqX6e3Gril\nrQUc8dUVa7VS30m2AP8S+Im+6zwNPN2m707yJeBF7fH0b84Yub+19NrX8x8Af9IurnUZr4tVLNfd\nwOuA89tmnYkt1zVY8atTJiHJ99F78f9gVX0EoKoe75u/2ufD2FXVY+334SQfpbeZ5PEkZ1bVobaJ\n5/Bm6LW5EPjM4vIc93I9ZtcAhnQAuCTJSUl2ADvp7ewb+BUV7YXiduCN7fq7gZvH3OOrgM9X1d9u\ngkgyld65FEhyVuv74bb6+q0k57X9BpdtQH+LPfVvG/0XwP1tek3LeIN6vQB4O/D6qvp2X33TLdcl\nNt1Xp7TlcS3wYFX9Tl99rc+Hjej1uUl+YHGa3sEA97eedrdh/f/TB4DL0nMe8NTipqIN9Ky1/7Ev\n10ns5R73T1tQ8/Te3T0O3No37yp6e8y/QN+RHfSOAPhim3dVX/2stmDn6K2WnTTm3j8A/MKS2r8C\nHqC31/8zwD/vmzfTnhRfAn6f9uG+DVjGfwjcB9zbnoxnDruMN6DXOXrbS+9pP+/brMt1QO8TWWbL\n9PNP6W1quLdveV40zPNhA3o9q/1tP9v+zle1+mnAbcBD7feprR7gva3X++g7gnCD+n0O8L+A5/fV\nxrpc/SSwJHVU1zYBSZIaA0CSOsoAkKSOMgAkqaMMAEnqKANAkjrKAJCkjjIAJKmj/j+KyeNwDXAV\ncAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcddb7497b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(mean_epi, 128);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def vals2colors(vals):\n",
    "    cmap = np.vstack((plt.cm.viridis(range(256)), plt.cm.cool(range(256))))\n",
    "    indices = 256 + (255 * vals).astype(int)\n",
    "    return cmap[indices, 0:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFhxJREFUeJzt3X+w5XV93/HnKxBJU7UsciHrLmTBWWzRtoveQaaMCQbD\nr3QAW02WmchqaFctpHXSzgjaGR0zTEkaZWSSYNdIgVZBhFC2cS1ZUeN0hlUuBvkpcvkRue6WXcWg\nGdJtwHf/ON+rX3fPvffce849l8v3+Zg5c7/n/f18v9/3/e7d+7rfH+ecVBWSpG76mZVuQJK0cgwB\nSeowQ0CSOswQkKQOMwQkqcMMAUnqMENAkjrMEJCkDlswBJIck+RLSR5K8kCSf9fUj0iyM8kjzdc1\nTT1JrkoyneTeJK9rrWtLM/6RJFuW79uSJA0iC71iOMlaYG1VfT3Jy4C7gfOBdwBPV9UVSS4F1lTV\n+5KcA/w2cA7wBuBjVfWGJEcAU8AkUM16Xl9V359v+0ceeWRt2LBhmO9Rkjrl7rvv/m5VTQwy9tCF\nBlTVHmBPM/3DJA8B64DzgNOaYdcBXwbe19Svr1667EpyeBMkpwE7q+ppgCQ7gbOAG+bb/oYNG5ia\nmhrke5EkAUn+atCxi7omkGQDcBLwVeDoJiBmg+KoZtg64MnWYjNNba66JGmFDBwCSV4K3AK8t6p+\nMN/QPrWap95vW1uTTCWZ2rdv36AtSpIWaaAQSPKz9ALgU1X1p035qeY0z+x1g71NfQY4prX4emD3\nPPWDVNW2qpqsqsmJiYFOa0mSlmCQu4MCfBJ4qKo+2pq1HZi9w2cLcFurfmFzl9ApwDPN6aLbgTOS\nrGnuJDqjqUmSVsiCF4aBU4G3A/cluaepvR+4ArgpyUXAt4G3NfN20LszaBp4FngnQFU9neR3gbua\ncR+evUgsSVoZC94iutImJyfLu4MkaXBJ7q6qyUHG+ophSeowQ0CSOswQkKQOG+TCsKQl2nDp5w6q\nPXHFr61AJ1J/HglIUocZApLUYYaAJHWYISBJHWYISFKHGQKS1GGGgCR1mCEgSR1mCEhShxkCktRh\nhoAkdZghIEkdZghIUocZApLUYYN80Pw1SfYmub9V+0ySe5rHE7OfPZxkQ5K/bc37eGuZ1ye5L8l0\nkquaD7CXJK2gQT5P4FrgD4HrZwtV9Ruz00k+AjzTGv9oVW3qs56rga3ALnofRn8W8PnFtyxJGpUF\njwSq6ivA0/3mNX/N/zpww3zrSLIWeHlV3Vm9T7a/Hjh/8e1KkkZp2GsCbwSeqqpHWrXjkvxlkr9I\n8samtg6YaY2ZaWqSpBU07MdLXsBPHwXsAY6tqu8leT3wP5K8Buh3/r/mWmmSrfROHXHssccO2aIk\naS5LPhJIcijwL4DPzNaqan9Vfa+Zvht4FDiB3l/+61uLrwd2z7XuqtpWVZNVNTkxMbHUFiVJCxjm\ndNCbgW9W1Y9P8ySZSHJIM308sBF4rKr2AD9MckpzHeFC4LYhti1JGoFBbhG9AbgTeHWSmSQXNbM2\nc/AF4V8C7k3yDeBm4N1VNXtR+T3AnwDT9I4QvDNIklbYgtcEquqCOerv6FO7BbhljvFTwGsX2Z8k\naRn5imFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOswQkKQOMwQkqcMMAUnqMENAkjrMEJCkDjME\nJKnDDAFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOmyQzxi+JsneJPe3ah9K8p0k9zSPc1rzLksy\nneThJGe26mc1tekkl47+W5EkLdYgRwLXAmf1qV9ZVZuaxw6AJCfS+wD61zTL/HGSQ5IcAvwRcDZw\nInBBM1aStIIG+aD5ryTZMOD6zgNurKr9wONJpoGTm3nTVfUYQJIbm7EPLrpjSdLIDHNN4JIk9zan\ni9Y0tXXAk60xM01trnpfSbYmmUoytW/fviFalCTNZ6khcDXwKmATsAf4SFNPn7E1T72vqtpWVZNV\nNTkxMbHEFiVJC1nwdFA/VfXU7HSSTwB/1jydAY5pDV0P7G6m56pLklbIko4EkqxtPX0LMHvn0HZg\nc5LDkhwHbAS+BtwFbExyXJKX0Lt4vH3pbUuSRmHBI4EkNwCnAUcmmQE+CJyWZBO9UzpPAO8CqKoH\nktxE74Lvc8DFVfV8s55LgNuBQ4BrquqBkX83kqRFGeTuoAv6lD85z/jLgcv71HcAOxbVnSRpWfmK\nYUnqMENAkjrMEJCkDjMEJKnDDAFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOswQkKQOMwQkqcMM\nAUnqMENAkjrMEJCkDjMEJKnDDAFJ6jBDQJI6bMEQSHJNkr1J7m/V/nOSbya5N8mtSQ5v6huS/G2S\ne5rHx1vLvD7JfUmmk1yVJMvzLUmSBjXIkcC1wFkH1HYCr62qfwJ8C7isNe/RqtrUPN7dql8NbAU2\nNo8D1ylJGrMFQ6CqvgI8fUDtz6vquebpLmD9fOtIshZ4eVXdWVUFXA+cv7SWJUmjMoprAr8FfL71\n/Lgkf5nkL5K8samtA2ZaY2aaWl9JtiaZSjK1b9++EbQoSepnqBBI8gHgOeBTTWkPcGxVnQT8DvDp\nJC8H+p3/r7nWW1XbqmqyqiYnJiaGaVGSNI9Dl7pgki3APwdOb07xUFX7gf3N9N1JHgVOoPeXf/uU\n0Xpg91K3LUkajSUdCSQ5C3gfcG5VPduqTyQ5pJk+nt4F4Meqag/wwySnNHcFXQjcNnT3kqShLHgk\nkOQG4DTgyCQzwAfp3Q10GLCzudNzV3Mn0C8BH07yHPA88O6qmr2o/B56dxr9PXrXENrXESRJK2DB\nEKiqC/qUPznH2FuAW+aYNwW8dlHdSZKWla8YlqQOMwQkqcMMAUnqMENAkjrMEJCkDjMEJKnDDAFJ\n6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOswQkKQOMwQkqcMMAUnqMENAkjrMEJCkDjMEJKnDBgqB\nJNck2Zvk/lbtiCQ7kzzSfF3T1JPkqiTTSe5N8rrWMlua8Y8k2TL6b0eStBiDHglcC5x1QO1S4I6q\n2gjc0TwHOBvY2Dy2AldDLzTofUj9G4CTgQ/OBockaWUMFAJV9RXg6QPK5wHXNdPXAee36tdXzy7g\n8CRrgTOBnVX1dFV9H9jJwcEiSRqjYa4JHF1VewCar0c19XXAk61xM01trvpBkmxNMpVkat++fUO0\nKEmaz3JcGE6fWs1TP7hYta2qJqtqcmJiYqTNSZJ+YpgQeKo5zUPzdW9TnwGOaY1bD+yepy5JWiHD\nhMB2YPYOny3Aba36hc1dQqcAzzSni24HzkiyprkgfEZTkyStkEMHGZTkBuA04MgkM/Tu8rkCuCnJ\nRcC3gbc1w3cA5wDTwLPAOwGq6ukkvwvc1Yz7cFUdeLFZkjRGA4VAVV0wx6zT+4wt4OI51nMNcM3A\n3UmSlpWvGJakDjMEJKnDDAFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOswQkKQOMwQkqcMMAUnq\nMENAkjrMEJCkDjMEJKnDDAFJ6jBDQJI6zBCQpA5bcggkeXWSe1qPHyR5b5IPJflOq35Oa5nLkkwn\neTjJmaP5FiRJSzXQx0v2U1UPA5sAkhwCfAe4ld5nCl9ZVX/QHp/kRGAz8BrglcAXkpxQVc8vtQdJ\n0nBGdTrodODRqvqrecacB9xYVfur6nF6H0R/8oi2L0laglGFwGbghtbzS5Lcm+SaJGua2jrgydaY\nmaYmSVohQ4dAkpcA5wKfbUpXA6+id6poD/CR2aF9Fq851rk1yVSSqX379g3boiRpDqM4Ejgb+HpV\nPQVQVU9V1fNV9SPgE/zklM8McExrufXA7n4rrKptVTVZVZMTExMjaFGS1M8oQuACWqeCkqxtzXsL\ncH8zvR3YnOSwJMcBG4GvjWD7kqQlWvLdQQBJfh74VeBdrfLvJ9lE71TPE7PzquqBJDcBDwLPARd7\nZ5AkrayhQqCqngVecUDt7fOMvxy4fJhtSpJGx1cMS1KHGQKS1GGGgCR1mCEgSR1mCEhShxkCktRh\nhoAkdZghIEkdZghIUocZApLUYYaAJHWYISBJHWYISFKHGQKS1GGGgCR1mCEgSR1mCEhShxkCktRh\nQ4dAkieS3JfkniRTTe2IJDuTPNJ8XdPUk+SqJNNJ7k3yumG3L0laulEdCbypqjZV1WTz/FLgjqra\nCNzRPAc4G9jYPLYCV49o+5KkJViu00HnAdc109cB57fq11fPLuDwJGuXqQdJ0gJGEQIF/HmSu5Ns\nbWpHV9UegObrUU19HfBka9mZpiZJWgGHjmAdp1bV7iRHATuTfHOeselTq4MG9cJkK8Cxxx47ghYl\nSf0MfSRQVbubr3uBW4GTgadmT/M0X/c2w2eAY1qLrwd291nntqqarKrJiYmJYVuUJM1hqBBI8veT\nvGx2GjgDuB/YDmxphm0BbmumtwMXNncJnQI8M3vaSJI0fsOeDjoauDXJ7Lo+XVX/K8ldwE1JLgK+\nDbytGb8DOAeYBp4F3jnk9iVJQxgqBKrqMeCf9ql/Dzi9T72Ai4fZprTabbj0c0te9okrfm2EnUij\nuTAsqTHML/hRbc+g0GL4thGS1GGGgCR1mKeDpFVk3Keb9OLnkYAkdZghIEkdZghIUod5TUBaIs/P\n68XAEJA6wNcTaC6GgPQi4xGKFsNrApLUYYaAJHWYISBJHWYISFKHGQKS1GHeHSQdwNsp1SUeCUhS\nh3kkIHXUgUc8Hu1005JDIMkxwPXALwA/ArZV1ceSfAj418C+Zuj7q2pHs8xlwEXA88C/rarbh+h9\nQR7Wa1R8AZZerIY5EngO+PdV9fUkLwPuTrKzmXdlVf1Be3CSE4HNwGuAVwJfSHJCVT0/RA+SRsQ/\nmrppydcEqmpPVX29mf4h8BCwbp5FzgNurKr9VfU4MA2cvNTtS5KGN5JrAkk2ACcBXwVOBS5JciEw\nRe9o4fv0AmJXa7EZ5g8NaSw81aMuG/ruoCQvBW4B3ltVPwCuBl4FbAL2AB+ZHdpn8ZpjnVuTTCWZ\n2rdvX78hkqQRGOpIIMnP0guAT1XVnwJU1VOt+Z8A/qx5OgMc01p8PbC733qrahuwDWBycrJvUEhL\n4V/90k8b5u6gAJ8EHqqqj7bqa6tqT/P0LcD9zfR24NNJPkrvwvBG4GtL3b6k5efF4he/YY4ETgXe\nDtyX5J6m9n7ggiSb6J3qeQJ4F0BVPZDkJuBBencWXeydQZK0spYcAlX1v+l/nn/HPMtcDly+1G1K\nWnkeHby4+IphSUPz1cerlyGgFy0vAksL8w3kJKnDDAFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowbxHV\ni4K3g76w+IKy1cMQkDQWBsMLk6eDJKnDPBLQquTpH2k0DAFJK8ZTRCvP00GS1GGGgCR1mKeD9ILi\nuX75ttTj5ZGAJHWYRwIdMMzFt0H+KvOvdy2nQX++PGJYmrGHQJKzgI8BhwB/UlVXjLuH5TLoL9tB\nxo1yXYPyl7lWM+80WpqxhkCSQ4A/An4VmAHuSrK9qh4cVw/j/kU36PZG2ddyfo8GhVYTry8sbNxH\nAicD01X1GECSG4HzgLGFwGq2WoJCeqHy1NLBxh0C64AnW89ngDeMuQdJmtdy/5H0QgqZcYdA+tTq\noEHJVmBr8/Rvkjy8xO0dCXx3icuOiz2Ozmro0x5HZzX02bfH/N6yb/cXBx047hCYAY5pPV8P7D5w\nUFVtA7YNu7EkU1U1Oex6lpM9js5q6NMeR2c19Lkaehz36wTuAjYmOS7JS4DNwPYx9yBJaoz1SKCq\nnktyCXA7vVtEr6mqB8bZgyTpJ8b+OoGq2gHsGNPmhj6lNAb2ODqroU97HJ3V0OcLvsdUHXRdVpLU\nEb53kCR12KoPgSRvS/JAkh8lmfMqfJKzkjycZDrJpa36cUm+muSRJJ9pLliPuscjkuxstrEzyZo+\nY96U5J7W4/8mOb+Zd22Sx1vzNq1Ej82451t9bG/Vl30/Dtpnkk1J7mx+Lu5N8hutecu2L+f6GWvN\nP6zZN9PNvtrQmndZU384yZmj6mkJPf5Okgeb/XZHkl9szev7b78CPb4jyb5WL/+qNW9L87PxSJIt\nK9jjla3+vpXkr1vzxrIfB1ZVq/oB/CPg1cCXgck5xhwCPAocD7wE+AZwYjPvJmBzM/1x4D3L0OPv\nA5c205cCv7fA+COAp4Gfb55fC7x1mffjQD0CfzNHfdn346B9AicAG5vpVwJ7gMOXc1/O9zPWGvNv\ngI8305uBzzTTJzbjDwOOa9ZzyAr1+KbWz917Znuc799+BXp8B/CHfZY9Anis+bqmmV6zEj0eMP63\n6d0EM7b9uJjHqj8SqKqHqmqhF5P9+O0qqur/ATcC5yUJ8CvAzc2464Dzl6HN85p1D7qNtwKfr6pn\nl6GXuSy2xx8b436EAfqsqm9V1SPN9G5gLzCxTP3M6vszdsCYdu83A6c3++484Maq2l9VjwPTzfrG\n3mNVfan1c7eL3mt5xmmQ/TiXM4GdVfV0VX0f2Amc9QLo8QLghmXoYyRWfQgMqN/bVawDXgH8dVU9\nd0B91I6uqj0AzdejFhi/mYN/aC5vDtGvTHLYCvb4c0mmkuyaPV3F+PbjYvoEIMnJ9P5ae7RVXo59\nOdfPWN8xzb56ht6+G2TZcfXYdhHw+dbzfv/2ozZoj/+y+Te8OcnsC1BfcPuxOZ12HPDFVnkc+3Fg\nq+LzBJJ8AfiFPrM+UFW3DbKKPrWap75o8/W4yPWsBf4xvddSzLoM+D/0fpltA94HfHiFejy2qnYn\nOR74YpL7gB/0Gbfk285GvC//G7Clqn7UlEeyL/ttrk/twH2w7D+HCxh4O0l+E5gEfrlVPujfvqoe\n7bf8Mvf4P4Ebqmp/knfTO7r6lQGXHYXFbGczcHNVPd+qjWM/DmxVhEBVvXnIVcz1dhXfBQ5Pcmjz\nl1nft7EYtsckTyVZW1V7ml9Me+dZ1a8Dt1bV37XWvaeZ3J/kvwL/YaV6bE6vUFWPJfkycBJwCyPa\nj6PqM8nLgc8B/7GqdrXWPZJ92ccgb4kyO2YmyaHAP6B37Wegt1MZU48keTO9wP3lqto/W5/j337U\nv7wW7LGqvtd6+glg9p14ZoDTDlj2yyPub3Y7g/57bQYubhfGtB8H1pXTQX3frqJ6V2m+RO8cPMAW\nYJAji8Xa3qx7kG0cdP6w+WU3e+79fOD+legxyZrZ0ydJjgROBR4c434ctM+XALcC11fVZw+Yt1z7\ncpC3RGn3/lbgi82+2w5sbu4eOg7YCHxtRH0tqsckJwH/BTi3qva26n3/7Veox7Wtp+cCDzXTtwNn\nNL2uAc7gp4+ox9Zj0+er6V2gvrNVG9d+HNxKX5ke9gG8hV4y7weeAm5v6q8EdrTGnQN8i17ifqBV\nP57ef7hp4LPAYcvQ4yuAO4BHmq9HNPVJep+uNjtuA/Ad4GcOWP6LwH30fmH9d+ClK9Ej8M+aPr7R\nfL1onPtxEX3+JvB3wD2tx6bl3pf9fsbonWo6t5n+uWbfTDf76vjWsh9olnsYOHsZ/78s1OMXmv9H\ns/tt+0L/9ivQ438CHmh6+RLwD1vL/lazf6eBd65Uj83zDwFXHLDc2PbjoA9fMSxJHdaV00GSpD4M\nAUnqMENAkjrMEJCkDjMEJKnDDAFJ6jBDQJI6zBCQpA77/ySjrQsqRB8gAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcdd633f278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "epi_vals = np.minimum(np.maximum(mean_epi/(3 * mean_std), -1), 1)\n",
    "colors = vals2colors(epi_vals)\n",
    "plt.hist(epi_vals, 64);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do not have lh.smoothwm you can retrieve it using the following command\n",
    "\n",
    "```\n",
    "!datalad get /data/ds000114/derivatives/freesurfer/fsaverage5/surf/lh.smoothwm\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "08243ade7eff4bc7be8584dfaa27c104",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_fs = nb.freesurfer.read_geometry('/data/ds000114/derivatives/freesurfer/fsaverage5/surf/lh.smoothwm')\n",
    "\n",
    "fig = p3.figure(width=600, height=600)\n",
    "# we draw the tetrahedron\n",
    "x, y, z = data_fs[0][:, 0], data_fs[0][:, 1], data_fs[0][:, 2]\n",
    "p3.plot_trisurf(x, y, z, triangles=data_fs[1], color=colors)\n",
    "fig.camera_fov = 1\n",
    "fig.anglex = np.pi/2\n",
    "fig.angley = -np.pi/2\n",
    "fig.xlim = (-100, 100)\n",
    "fig.ylim = (-100, 100)\n",
    "fig.zlim = (-100, 100)\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "## Exercise 4:\n",
    "\n",
    "Plot the mean activity of the functional time series (stored as gifti) on the surface of the fsaverage5 mesh. The fmriprep routine generates these surface-based timeseries files.\n",
    "\n",
    "- extract principal components of the data scikit-learn \n",
    "- plot the first two components on the surface of the mesh\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true,
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "foo = PCA(n_components=2).fit_transform(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFhxJREFUeJzt3X+w5XV93/HnKxBJU7UsciHrLmTBWWzRtoveQaaMCQbD\nr3QAW02WmchqaFctpHXSzgjaGR0zTEkaZWSSYNdIgVZBhFC2cS1ZUeN0hlUuBvkpcvkRue6WXcWg\nGdJtwHf/ON+rX3fPvffce849l8v3+Zg5c7/n/f18v9/3/e7d+7rfH+ecVBWSpG76mZVuQJK0cgwB\nSeowQ0CSOswQkKQOMwQkqcMMAUnqMENAkjrMEJCkDlswBJIck+RLSR5K8kCSf9fUj0iyM8kjzdc1\nTT1JrkoyneTeJK9rrWtLM/6RJFuW79uSJA0iC71iOMlaYG1VfT3Jy4C7gfOBdwBPV9UVSS4F1lTV\n+5KcA/w2cA7wBuBjVfWGJEcAU8AkUM16Xl9V359v+0ceeWRt2LBhmO9Rkjrl7rvv/m5VTQwy9tCF\nBlTVHmBPM/3DJA8B64DzgNOaYdcBXwbe19Svr1667EpyeBMkpwE7q+ppgCQ7gbOAG+bb/oYNG5ia\nmhrke5EkAUn+atCxi7omkGQDcBLwVeDoJiBmg+KoZtg64MnWYjNNba66JGmFDBwCSV4K3AK8t6p+\nMN/QPrWap95vW1uTTCWZ2rdv36AtSpIWaaAQSPKz9ALgU1X1p035qeY0z+x1g71NfQY4prX4emD3\nPPWDVNW2qpqsqsmJiYFOa0mSlmCQu4MCfBJ4qKo+2pq1HZi9w2cLcFurfmFzl9ApwDPN6aLbgTOS\nrGnuJDqjqUmSVsiCF4aBU4G3A/cluaepvR+4ArgpyUXAt4G3NfN20LszaBp4FngnQFU9neR3gbua\ncR+evUgsSVoZC94iutImJyfLu4MkaXBJ7q6qyUHG+ophSeowQ0CSOswQkKQOG+TCsKQl2nDp5w6q\nPXHFr61AJ1J/HglIUocZApLUYYaAJHWYISBJHWYISFKHGQKS1GGGgCR1mCEgSR1mCEhShxkCktRh\nhoAkdZghIEkdZghIUocZApLUYYN80Pw1SfYmub9V+0ySe5rHE7OfPZxkQ5K/bc37eGuZ1ye5L8l0\nkquaD7CXJK2gQT5P4FrgD4HrZwtV9Ruz00k+AjzTGv9oVW3qs56rga3ALnofRn8W8PnFtyxJGpUF\njwSq6ivA0/3mNX/N/zpww3zrSLIWeHlV3Vm9T7a/Hjh/8e1KkkZp2GsCbwSeqqpHWrXjkvxlkr9I\n8samtg6YaY2ZaWqSpBU07MdLXsBPHwXsAY6tqu8leT3wP5K8Buh3/r/mWmmSrfROHXHssccO2aIk\naS5LPhJIcijwL4DPzNaqan9Vfa+Zvht4FDiB3l/+61uLrwd2z7XuqtpWVZNVNTkxMbHUFiVJCxjm\ndNCbgW9W1Y9P8ySZSHJIM308sBF4rKr2AD9MckpzHeFC4LYhti1JGoFBbhG9AbgTeHWSmSQXNbM2\nc/AF4V8C7k3yDeBm4N1VNXtR+T3AnwDT9I4QvDNIklbYgtcEquqCOerv6FO7BbhljvFTwGsX2Z8k\naRn5imFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOswQkKQOMwQkqcMMAUnqMENAkjrMEJCkDjME\nJKnDDAFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOmyQzxi+JsneJPe3ah9K8p0k9zSPc1rzLksy\nneThJGe26mc1tekkl47+W5EkLdYgRwLXAmf1qV9ZVZuaxw6AJCfS+wD61zTL/HGSQ5IcAvwRcDZw\nInBBM1aStIIG+aD5ryTZMOD6zgNurKr9wONJpoGTm3nTVfUYQJIbm7EPLrpjSdLIDHNN4JIk9zan\ni9Y0tXXAk60xM01trnpfSbYmmUoytW/fviFalCTNZ6khcDXwKmATsAf4SFNPn7E1T72vqtpWVZNV\nNTkxMbHEFiVJC1nwdFA/VfXU7HSSTwB/1jydAY5pDV0P7G6m56pLklbIko4EkqxtPX0LMHvn0HZg\nc5LDkhwHbAS+BtwFbExyXJKX0Lt4vH3pbUuSRmHBI4EkNwCnAUcmmQE+CJyWZBO9UzpPAO8CqKoH\nktxE74Lvc8DFVfV8s55LgNuBQ4BrquqBkX83kqRFGeTuoAv6lD85z/jLgcv71HcAOxbVnSRpWfmK\nYUnqMENAkjrMEJCkDjMEJKnDDAFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOswQkKQOMwQkqcMM\nAUnqMENAkjrMEJCkDjMEJKnDDAFJ6jBDQJI6bMEQSHJNkr1J7m/V/nOSbya5N8mtSQ5v6huS/G2S\ne5rHx1vLvD7JfUmmk1yVJMvzLUmSBjXIkcC1wFkH1HYCr62qfwJ8C7isNe/RqtrUPN7dql8NbAU2\nNo8D1ylJGrMFQ6CqvgI8fUDtz6vquebpLmD9fOtIshZ4eVXdWVUFXA+cv7SWJUmjMoprAr8FfL71\n/Lgkf5nkL5K8samtA2ZaY2aaWl9JtiaZSjK1b9++EbQoSepnqBBI8gHgOeBTTWkPcGxVnQT8DvDp\nJC8H+p3/r7nWW1XbqmqyqiYnJiaGaVGSNI9Dl7pgki3APwdOb07xUFX7gf3N9N1JHgVOoPeXf/uU\n0Xpg91K3LUkajSUdCSQ5C3gfcG5VPduqTyQ5pJk+nt4F4Meqag/wwySnNHcFXQjcNnT3kqShLHgk\nkOQG4DTgyCQzwAfp3Q10GLCzudNzV3Mn0C8BH07yHPA88O6qmr2o/B56dxr9PXrXENrXESRJK2DB\nEKiqC/qUPznH2FuAW+aYNwW8dlHdSZKWla8YlqQOMwQkqcMMAUnqMENAkjrMEJCkDjMEJKnDDAFJ\n6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOswQkKQOMwQkqcMMAUnqMENAkjrMEJCkDjMEJKnDBgqB\nJNck2Zvk/lbtiCQ7kzzSfF3T1JPkqiTTSe5N8rrWMlua8Y8k2TL6b0eStBiDHglcC5x1QO1S4I6q\n2gjc0TwHOBvY2Dy2AldDLzTofUj9G4CTgQ/OBockaWUMFAJV9RXg6QPK5wHXNdPXAee36tdXzy7g\n8CRrgTOBnVX1dFV9H9jJwcEiSRqjYa4JHF1VewCar0c19XXAk61xM01trvpBkmxNMpVkat++fUO0\nKEmaz3JcGE6fWs1TP7hYta2qJqtqcmJiYqTNSZJ+YpgQeKo5zUPzdW9TnwGOaY1bD+yepy5JWiHD\nhMB2YPYOny3Aba36hc1dQqcAzzSni24HzkiyprkgfEZTkyStkEMHGZTkBuA04MgkM/Tu8rkCuCnJ\nRcC3gbc1w3cA5wDTwLPAOwGq6ukkvwvc1Yz7cFUdeLFZkjRGA4VAVV0wx6zT+4wt4OI51nMNcM3A\n3UmSlpWvGJakDjMEJKnDDAFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowQ0CSOswQkKQOMwQkqcMMAUnq\nMENAkjrMEJCkDjMEJKnDDAFJ6jBDQJI6zBCQpA5bcggkeXWSe1qPHyR5b5IPJflOq35Oa5nLkkwn\neTjJmaP5FiRJSzXQx0v2U1UPA5sAkhwCfAe4ld5nCl9ZVX/QHp/kRGAz8BrglcAXkpxQVc8vtQdJ\n0nBGdTrodODRqvqrecacB9xYVfur6nF6H0R/8oi2L0laglGFwGbghtbzS5Lcm+SaJGua2jrgydaY\nmaYmSVohQ4dAkpcA5wKfbUpXA6+id6poD/CR2aF9Fq851rk1yVSSqX379g3boiRpDqM4Ejgb+HpV\nPQVQVU9V1fNV9SPgE/zklM8McExrufXA7n4rrKptVTVZVZMTExMjaFGS1M8oQuACWqeCkqxtzXsL\ncH8zvR3YnOSwJMcBG4GvjWD7kqQlWvLdQQBJfh74VeBdrfLvJ9lE71TPE7PzquqBJDcBDwLPARd7\nZ5AkrayhQqCqngVecUDt7fOMvxy4fJhtSpJGx1cMS1KHGQKS1GGGgCR1mCEgSR1mCEhShxkCktRh\nhoAkdZghIEkdZghIUocZApLUYYaAJHWYISBJHWYISFKHGQKS1GGGgCR1mCEgSR1mCEhShxkCktRh\nQ4dAkieS3JfkniRTTe2IJDuTPNJ8XdPUk+SqJNNJ7k3yumG3L0laulEdCbypqjZV1WTz/FLgjqra\nCNzRPAc4G9jYPLYCV49o+5KkJViu00HnAdc109cB57fq11fPLuDwJGuXqQdJ0gJGEQIF/HmSu5Ns\nbWpHV9UegObrUU19HfBka9mZpiZJWgGHjmAdp1bV7iRHATuTfHOeselTq4MG9cJkK8Cxxx47ghYl\nSf0MfSRQVbubr3uBW4GTgadmT/M0X/c2w2eAY1qLrwd291nntqqarKrJiYmJYVuUJM1hqBBI8veT\nvGx2GjgDuB/YDmxphm0BbmumtwMXNncJnQI8M3vaSJI0fsOeDjoauDXJ7Lo+XVX/K8ldwE1JLgK+\nDbytGb8DOAeYBp4F3jnk9iVJQxgqBKrqMeCf9ql/Dzi9T72Ai4fZprTabbj0c0te9okrfm2EnUij\nuTAsqTHML/hRbc+g0GL4thGS1GGGgCR1mKeDpFVk3Keb9OLnkYAkdZghIEkdZghIUod5TUBaIs/P\n68XAEJA6wNcTaC6GgPQi4xGKFsNrApLUYYaAJHWYISBJHWYISFKHGQKS1GHeHSQdwNsp1SUeCUhS\nh3kkIHXUgUc8Hu1005JDIMkxwPXALwA/ArZV1ceSfAj418C+Zuj7q2pHs8xlwEXA88C/rarbh+h9\nQR7Wa1R8AZZerIY5EngO+PdV9fUkLwPuTrKzmXdlVf1Be3CSE4HNwGuAVwJfSHJCVT0/RA+SRsQ/\nmrppydcEqmpPVX29mf4h8BCwbp5FzgNurKr9VfU4MA2cvNTtS5KGN5JrAkk2ACcBXwVOBS5JciEw\nRe9o4fv0AmJXa7EZ5g8NaSw81aMuG/ruoCQvBW4B3ltVPwCuBl4FbAL2AB+ZHdpn8ZpjnVuTTCWZ\n2rdvX78hkqQRGOpIIMnP0guAT1XVnwJU1VOt+Z8A/qx5OgMc01p8PbC733qrahuwDWBycrJvUEhL\n4V/90k8b5u6gAJ8EHqqqj7bqa6tqT/P0LcD9zfR24NNJPkrvwvBG4GtL3b6k5efF4he/YY4ETgXe\nDtyX5J6m9n7ggiSb6J3qeQJ4F0BVPZDkJuBBencWXeydQZK0spYcAlX1v+l/nn/HPMtcDly+1G1K\nWnkeHby4+IphSUPz1cerlyGgFy0vAksL8w3kJKnDDAFJ6jBDQJI6zBCQpA4zBCSpwwwBSeowbxHV\ni4K3g76w+IKy1cMQkDQWBsMLk6eDJKnDPBLQquTpH2k0DAFJK8ZTRCvP00GS1GGGgCR1mKeD9ILi\nuX75ttTj5ZGAJHWYRwIdMMzFt0H+KvOvdy2nQX++PGJYmrGHQJKzgI8BhwB/UlVXjLuH5TLoL9tB\nxo1yXYPyl7lWM+80WpqxhkCSQ4A/An4VmAHuSrK9qh4cVw/j/kU36PZG2ddyfo8GhVYTry8sbNxH\nAicD01X1GECSG4HzgLGFwGq2WoJCeqHy1NLBxh0C64AnW89ngDeMuQdJmtdy/5H0QgqZcYdA+tTq\noEHJVmBr8/Rvkjy8xO0dCXx3icuOiz2Ozmro0x5HZzX02bfH/N6yb/cXBx047hCYAY5pPV8P7D5w\nUFVtA7YNu7EkU1U1Oex6lpM9js5q6NMeR2c19Lkaehz36wTuAjYmOS7JS4DNwPYx9yBJaoz1SKCq\nnktyCXA7vVtEr6mqB8bZgyTpJ8b+OoGq2gHsGNPmhj6lNAb2ODqroU97HJ3V0OcLvsdUHXRdVpLU\nEb53kCR12KoPgSRvS/JAkh8lmfMqfJKzkjycZDrJpa36cUm+muSRJJ9pLliPuscjkuxstrEzyZo+\nY96U5J7W4/8mOb+Zd22Sx1vzNq1Ej82451t9bG/Vl30/Dtpnkk1J7mx+Lu5N8hutecu2L+f6GWvN\nP6zZN9PNvtrQmndZU384yZmj6mkJPf5Okgeb/XZHkl9szev7b78CPb4jyb5WL/+qNW9L87PxSJIt\nK9jjla3+vpXkr1vzxrIfB1ZVq/oB/CPg1cCXgck5xhwCPAocD7wE+AZwYjPvJmBzM/1x4D3L0OPv\nA5c205cCv7fA+COAp4Gfb55fC7x1mffjQD0CfzNHfdn346B9AicAG5vpVwJ7gMOXc1/O9zPWGvNv\ngI8305uBzzTTJzbjDwOOa9ZzyAr1+KbWz917Znuc799+BXp8B/CHfZY9Anis+bqmmV6zEj0eMP63\n6d0EM7b9uJjHqj8SqKqHqmqhF5P9+O0qqur/ATcC5yUJ8CvAzc2464Dzl6HN85p1D7qNtwKfr6pn\nl6GXuSy2xx8b436EAfqsqm9V1SPN9G5gLzCxTP3M6vszdsCYdu83A6c3++484Maq2l9VjwPTzfrG\n3mNVfan1c7eL3mt5xmmQ/TiXM4GdVfV0VX0f2Amc9QLo8QLghmXoYyRWfQgMqN/bVawDXgH8dVU9\nd0B91I6uqj0AzdejFhi/mYN/aC5vDtGvTHLYCvb4c0mmkuyaPV3F+PbjYvoEIMnJ9P5ae7RVXo59\nOdfPWN8xzb56ht6+G2TZcfXYdhHw+dbzfv/2ozZoj/+y+Te8OcnsC1BfcPuxOZ12HPDFVnkc+3Fg\nq+LzBJJ8AfiFPrM+UFW3DbKKPrWap75o8/W4yPWsBf4xvddSzLoM+D/0fpltA94HfHiFejy2qnYn\nOR74YpL7gB/0Gbfk285GvC//G7Clqn7UlEeyL/ttrk/twH2w7D+HCxh4O0l+E5gEfrlVPujfvqoe\n7bf8Mvf4P4Ebqmp/knfTO7r6lQGXHYXFbGczcHNVPd+qjWM/DmxVhEBVvXnIVcz1dhXfBQ5Pcmjz\nl1nft7EYtsckTyVZW1V7ml9Me+dZ1a8Dt1bV37XWvaeZ3J/kvwL/YaV6bE6vUFWPJfkycBJwCyPa\nj6PqM8nLgc8B/7GqdrXWPZJ92ccgb4kyO2YmyaHAP6B37Wegt1MZU48keTO9wP3lqto/W5/j337U\nv7wW7LGqvtd6+glg9p14ZoDTDlj2yyPub3Y7g/57bQYubhfGtB8H1pXTQX3frqJ6V2m+RO8cPMAW\nYJAji8Xa3qx7kG0cdP6w+WU3e+79fOD+legxyZrZ0ydJjgROBR4c434ctM+XALcC11fVZw+Yt1z7\ncpC3RGn3/lbgi82+2w5sbu4eOg7YCHxtRH0tqsckJwH/BTi3qva26n3/7Veox7Wtp+cCDzXTtwNn\nNL2uAc7gp4+ox9Zj0+er6V2gvrNVG9d+HNxKX5ke9gG8hV4y7weeAm5v6q8EdrTGnQN8i17ifqBV\nP57ef7hp4LPAYcvQ4yuAO4BHmq9HNPVJep+uNjtuA/Ad4GcOWP6LwH30fmH9d+ClK9Ej8M+aPr7R\nfL1onPtxEX3+JvB3wD2tx6bl3pf9fsbonWo6t5n+uWbfTDf76vjWsh9olnsYOHsZ/78s1OMXmv9H\ns/tt+0L/9ivQ438CHmh6+RLwD1vL/lazf6eBd65Uj83zDwFXHLDc2PbjoA9fMSxJHdaV00GSpD4M\nAUnqMENAkjrMEJCkDjMEJKnDDAFJ6jBDQJI6zBCQpA77/ySjrQsqRB8gAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcdd626db00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "comp_vals = np.minimum(np.maximum(foo[:, 0]/(3 * np.std(foo[:, 0])), -1), 1)\n",
    "colors = vals2colors(comp_vals)\n",
    "plt.hist(epi_vals, 64);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bb64357bb77e4562bf2368b3878e1069",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_fs = nb.freesurfer.read_geometry('/data/ds000114/derivatives/freesurfer/fsaverage5/surf/lh.smoothwm')\n",
    "\n",
    "fig = p3.figure(width=600, height=600)\n",
    "# we draw the tetrahedron\n",
    "x, y, z = data_fs[0][:, 0], data_fs[0][:, 1], data_fs[0][:, 2]\n",
    "p3.plot_trisurf(x, y, z, triangles=data_fs[1], color=colors)\n",
    "fig.camera_fov = 1\n",
    "fig.anglex = np.pi/2\n",
    "fig.angley = -np.pi/2\n",
    "fig.xlim = (-100, 100)\n",
    "fig.ylim = (-100, 100)\n",
    "fig.zlim = (-100, 100)\n",
    "p3.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Look up how to use PCA from scikit learn\n",
    "\n",
    "# Follow the steps from above to plot the two components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}