ely_good.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Putatively Useful Code from Ely's Test Notebook, Neurohackweek 2017\n",
    "\n",
    "- Cloned from the demo introductory_3d.ipynb notebook\n",
    "- Builds on \n",
    "[IPyvolume](https://ipyvolume.readthedocs.io/en/latest/index.html), \n",
    "[NiBabel cifti2 package](http://nipy.org/nibabel/reference/nibabel.cifti2.html), \n",
    "[Michiel Cottaar's cifti tools](https://git.fmrib.ox.ac.uk/ndcn0236/cifti), and [various NHW demo notebooks](https://github.com/neurohackweek)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load packages etc\n",
    "%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\n",
    "\n",
    "# be sure to set the global variable NHW to your NHW_visualization git directory\n",
    "# e.g. add 'export NHW=$HOME/Desktop/nhw2017/NHW_visualization' to your ~/.bashrc file\n",
    "nhw = os.environ['NHW'] # set the global variable NHW to your NHW_visualization git directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n"
     ]
    }
   ],
   "source": [
    "# specify some example data\n",
    "data_lh = nb.load(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.L.very_inflated_MSMAll.32k_fs_LR.surf.gii')\n",
    "data_rh = nb.load(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.R.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_myelin = nb.load(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.R.MyelinMap.32k_fs_LR.func.gii')\n",
    "data_sc = nb.load(nhw+'/data/100610/MNINonLinear/Atlas_ROIs.2.nii.gz')\n",
    "data_ciftiTS = nb.load(nhw+'/data/100610/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<b>CIFTI Tools:</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cread(infile):\n",
    "    # loads a cifti file and stores the data and four metadata objects\n",
    "    # should be called as: cd,ch,cn,cx,cf=cread(input_cifti)\n",
    "    cin=nb.load(infile)\n",
    "    cdata=np.asarray(cin.get_data()).copy()\n",
    "    chead=cin.header.copy()\n",
    "    cnhdr=cin.nifti_header.copy()\n",
    "    cxtra=cin.extra.copy()\n",
    "    cfmap=cin.file_map.copy()\n",
    "    return cdata,chead,cnhdr,cxtra,cfmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n"
     ]
    }
   ],
   "source": [
    "# cread example:\n",
    "input_cifti=nhw+'/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii'\n",
    "cd,ch,cn,cx,cf=cread(input_cifti)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def clthresh(data,thresh,setto=0):\n",
    "    # set all values below threshold to a given value (default 0)\n",
    "    lthr_data=data.copy()\n",
    "    lthr_data[lthr_data<thresh]=setto\n",
    "    return lthr_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cuthresh(data,thresh,setto=0):\n",
    "    # set all values above threshold to a given value (default 0)\n",
    "    uthr_data=data.copy()\n",
    "    uthr_data[uthr_data>thresh]=setto\n",
    "    return uthr_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEDCAYAAADOc0QpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAETdJREFUeJzt3X+MZWV9x/H3p+v6I4UW7U50s7AuNqStGit0shA0htjQ\n8iuhTWiytlFrTTZYNJrYP6gm/uhf2KQ24iqbrRKlseDP6kaXGqwkYlLAYbOsLFRdlYYlVEYMrBup\ndem3f9wDXocZ7pmZO3Nnnnm/kps595xnzvmek5vPPvOc555NVSFJasuvTboASdL4Ge6S1CDDXZIa\nZLhLUoMMd0lqkOEuSQ2aaLgnuT7Jw0nu6dH2H5Mc6l7fTfLoatQoSetRJjnPPclrgBPADVX18kX8\n3tuAs6vqr1asOElaxybac6+qbwA/GV6X5LeT/FuSu5LcluR35/nV1wE3rkqRkrQOPWvSBcxjH3Bl\nVX0vybnAR4HXPrkxyYuBM4GvT6g+SVrz1lS4JzkFOB/4bJInVz9nTrNdwOeq6onVrE2S1pM1Fe4M\nhokerapXPkObXcBVq1SPJK1La2oqZFUdB36Y5M8AMvD7T27vxt+fD/zHhEqUpHVh0lMhb2QQ1L+T\n5FiSNwN/Abw5yd3AEeDyoV/ZBdxUPspSkp7RRKdCSpJWxsiee5LnJrkzyd1JjiR5/zxtLkjy2NCX\njN6zMuVKkvroc0P158Brq+pEks3AN5PcXFW3z2l3W1Vd1vfAW7ZsqR07diyiVEnSXXfd9eOqmhrV\nbmS4d+PbJ7q3m7vXssdyduzYwczMzHJ3I0kbSpL/6tOu1w3VJJuSHAIeBm6pqjvmaXZ+ksNJbk7y\nsgX2szvJTJKZ2dnZPoeWJC1Br3Cvqie6ueenAzuTzH0OzEFge1W9Avgw8MUF9rOvqqaranpqauRf\nFZKkJVrUVMiqehS4FbhozvrjVXWiWz4AbE6yZWxVSpIWpc9smakkp3XLzwMuBP5zTpsXpXteQJKd\n3X4fGX+5kqQ++syW2Qp8MskmBqH9mar6cpIrAapqL3AF8JYkJ4HHgV1+0UiSJqfPbJnDwNnzrN87\ntLwH2DPe0iRJS7Wmni0jSRoPw12SGmS4S1KD1trz3KWn7Lj6KxM79v3XXDqxY0vjYM9dkhpkuEtS\ngwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXI\ncJekBhnuktQgw12SGmS4S1KDRoZ7kucmuTPJ3UmOJHn/PG2S5NokR5McTnLOypQrSeqjz/+h+nPg\ntVV1Islm4JtJbq6q24faXAyc1b3OBa7rfkqSJmBkz70GTnRvN3evmtPscuCGru3twGlJto63VElS\nX73G3JNsSnIIeBi4parumNNkG/DA0Ptj3bq5+9mdZCbJzOzs7FJrliSN0Cvcq+qJqnolcDqwM8nL\nl3KwqtpXVdNVNT01NbWUXUiSeljUbJmqehS4FbhozqYHgTOG3p/erZMkTcDIG6pJpoBfVNWjSZ4H\nXAh8YE6z/cBbk9zE4EbqY1X10NirlVbJjqu/MpHj3n/NpRM5rtrTZ7bMVuCTSTYx6Ol/pqq+nORK\ngKraCxwALgGOAj8D3rRC9UqSehgZ7lV1GDh7nvV7h5YLuGq8pUmSlspvqEpSgwx3SWqQ4S5JDTLc\nJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3q8+AwbXCTekKipKWz5y5JDTLcJalB\nhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQSPDPckZSW5Ncm+SI0nePk+bC5I8\nluRQ93rPypQrSeqjz7NlTgLvrKqDSU4F7kpyS1XdO6fdbVV12fhLlCQt1siee1U9VFUHu+WfAvcB\n21a6MEnS0i1qzD3JDuBs4I55Np+f5HCSm5O8bIHf351kJsnM7OzsoouVJPXTO9yTnAJ8HnhHVR2f\ns/kgsL2qXgF8GPjifPuoqn1VNV1V01NTU0utWZI0Qq9wT7KZQbB/qqq+MHd7VR2vqhPd8gFgc5It\nY61UktRbn9kyAT4O3FdVH1ygzYu6diTZ2e33kXEWKknqr89smVcBrwe+neRQt+5dwHaAqtoLXAG8\nJclJ4HFgV1XVCtQrSephZLhX1TeBjGizB9gzrqIkScvjN1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpk\nuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7\nJDXIcJekBhnuktQgw12SGmS4S1KDRoZ7kjOS3Jrk3iRHkrx9njZJcm2So0kOJzlnZcqVJPXxrB5t\nTgLvrKqDSU4F7kpyS1XdO9TmYuCs7nUucF33U5I0ASN77lX1UFUd7JZ/CtwHbJvT7HLghhq4HTgt\nydaxVytJ6mVRY+5JdgBnA3fM2bQNeGDo/TGe/g8ASXYnmUkyMzs7u7hKJUm99Q73JKcAnwfeUVXH\nl3KwqtpXVdNVNT01NbWUXUiSeugV7kk2Mwj2T1XVF+Zp8iBwxtD707t1kqQJ6DNbJsDHgfuq6oML\nNNsPvKGbNXMe8FhVPTTGOiVJi9BntsyrgNcD305yqFv3LmA7QFXtBQ4AlwBHgZ8Bbxp/qZKkvkaG\ne1V9E8iINgVcNa6iJEnL4zdUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWp\nQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho0\nMtyTXJ/k4ST3LLD9giSPJTnUvd4z/jIlSYvxrB5tPgHsAW54hja3VdVlY6lIkrRsI3vuVfUN4Cer\nUIskaUzGNeZ+fpLDSW5O8rIx7VOStER9hmVGOQhsr6oTSS4BvgicNV/DJLuB3QDbt28fw6ElSfNZ\nds+9qo5X1Ylu+QCwOcmWBdruq6rpqpqemppa7qElSQtYdrgneVGSdMs7u30+stz9SpKWbuSwTJIb\ngQuALUmOAe8FNgNU1V7gCuAtSU4CjwO7qqpWrGJJ0kgjw72qXjdi+x4GUyUlSWuE31CVpAYZ7pLU\nIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGjeOpkJLGZMfVX5nIce+/5tKJHFcrx567\nJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho0MtyT\nXJ/k4ST3LLA9Sa5NcjTJ4STnjL9MSdJi9Om5fwK46Bm2Xwyc1b12A9ctvyxJ0nKMDPeq+gbwk2do\ncjlwQw3cDpyWZOu4CpQkLd44xty3AQ8MvT/WrXuaJLuTzCSZmZ2dHcOhJUnzWdUbqlW1r6qmq2p6\nampqNQ8tSRvKOML9QeCMofend+skSRMyjnDfD7yhmzVzHvBYVT00hv1KkpZo5P+hmuRG4AJgS5Jj\nwHuBzQBVtRc4AFwCHAV+BrxppYqVJPUzMtyr6nUjthdw1dgqkiQtm99QlaQGGe6S1CDDXZIaZLhL\nUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1\nyHCXpAaN/J+YpJV06u9d/bR1P73vmglUIrXFnrskNcieu1bNfL30vu3szUuLY89dkhpkz13rwmr0\n5v2LQS3pFe5JLgI+BGwCPlZV18zZfgHwJeCH3aovVNXfjbFOrTN9h2BWw1qqRVotI8M9ySbgI8CF\nwDHgW0n2V9W9c5reVlWXrUCN0rzszUsL69Nz3wkcraofACS5CbgcmBvu0sTZS5cG+txQ3QY8MPT+\nWLdurvOTHE5yc5KXzbejJLuTzCSZmZ2dXUK5kqQ+xnVD9SCwvapOJLkE+CJw1txGVbUP2AcwPT1d\nYzq2tKruf+6f/8r7Hf/zLxOqRFpYn3B/EDhj6P3p3bqnVNXxoeUDST6aZEtV/Xg8ZUpr19ywBwNf\nk9dnWOZbwFlJzkzybGAXsH+4QZIXJUm3vLPb7yPjLlaS1M/InntVnUzyVuCrDKZCXl9VR5Jc2W3f\nC1wBvCXJSeBxYFdVOeyiDcuhG01arzH3qjoAHJizbu/Q8h5gz3hLkyQtlY8fkKQG+fgBLZtzy0fz\npqtWmz13SWqQPXdpQuzNayXZc5ekBtlzXyd2XP2VSZegVeAUSo2LPXdJapDhLkkNclhGWsO86aql\nsucuSQ2y5y6tM950VR+GuxbFb6NK64PDMpLUIHvu0jrnTVfNx567JDXInvsi+U1RSeuB4S41yBk1\nMtwlTfQv0vuvuXRix26Z4a4FOe1RWr8Md2kDcEbNxmO4C7CXLrXGcN+ADHKBN11b1yvck1wEfAjY\nBHysqq6Zsz3d9kuAnwF/WVUHx1yrRjC0tRwO3bRlZLgn2QR8BLgQOAZ8K8n+qrp3qNnFwFnd61zg\nuu7nilhLc80NVLXM3v361afnvhM4WlU/AEhyE3A5MBzulwM3VFUBtyc5LcnWqnpo7BUvk2EsLd18\nvfs+/Edh9fUJ923AA0Pvj/H0Xvl8bbYBvxLuSXYDu7u3J5J8Z1HV/qotwI+X8fut8DoMrNp1yGoc\nZOnW6OfhsgW35AMrdtA1ei2W7cV9Gq3qDdWq2gfsG8e+ksxU1fQ49rWeeR0GvA4DXodf2ujXos+D\nwx4Ezhh6f3q3brFtJEmrpE+4fws4K8mZSZ4N7AL2z2mzH3hDBs4DHluL4+2StFGMHJapqpNJ3gp8\nlcFUyOur6kiSK7vte4EDDKZBHmUwFfJNK1fyU8YyvNMAr8OA12HA6/BLG/paZDDBRZLUEv+zDklq\nkOEuSQ1al+Ge5KIk30lyNElz30pKcn+Sbyc5lGSmW/eCJLck+V738/lD7f+2uxbfSfLHQ+v/oNvP\n0STXdo+JWNOSXJ/k4ST3DK0b27kneU6ST3fr70iyYzXPr68FrsP7kjzYfS4OJblkaFtz1yHJGUlu\nTXJvkiNJ3t6t33CfhyWpqnX1YnBT9/vAS4BnA3cDL510XWM+x/uBLXPW/T1wdbd8NfCBbvml3TV4\nDnBmd202ddvuBM5j8L2bm4GLJ31uPc79NcA5wD0rce7AXwN7u+VdwKcnfc6LuA7vA/5mnrZNXgdg\nK3BOt3wq8N3uXDfc52Epr/XYc3/qcQhV9b/Ak49DaN3lwCe75U8CfzK0/qaq+nlV/ZDBjKWdSbYC\nv1FVt9fgk3vD0O+sWVX1DeAnc1aP89yH9/U54A/X4l80C1yHhTR5HarqoeoeQFhVPwXuY/DN9w33\neViK9RjuCz3qoCUFfC3JXd0jGwBeWL/87sB/Ay/slhe6Htu65bnr16NxnvtTv1NVJ4HHgN9ambJX\nxNuSHO6GbZ4cjmj+OnTDJWcDd+DnoZf1GO4bwaur6pUMnrZ5VZLXDG/seh8bcg7rRj53Bk9bfQnw\nSgbPbfqHyZazOpKcAnweeEdVHR/etsE/D89oPYZ78486qKoHu58PA//KYCjqR92fl3Q/H+6aL3Q9\nHuyW565fj8Z57k/9TpJnAb8JPLJilY9RVf2oqp6oqv8D/onB5wIavg5JNjMI9k9V1Re61X4eeliP\n4d7ncQjrVpJfT3Lqk8vAHwH3MDjHN3bN3gh8qVveD+zq7vqfyeCZ+nd2f7YeT3JeN4b4hqHfWW/G\nee7D+7oC+HrX+1vzngy0zp8y+FxAo9ehq/njwH1V9cGhTX4e+pj0Hd2lvBg86uC7DO6Gv3vS9Yz5\n3F7C4I7/3cCRJ8+PwTjgvwPfA74GvGDod97dXYvvMDQjBphmEADfB/bQfSN5Lb+AGxkMOfyCwdjo\nm8d57sBzgc8yuNl2J/CSSZ/zIq7DPwPfBg4zCKWtLV8H4NUMhlwOA4e61yUb8fOwlJePH5CkBq3H\nYRlJ0giGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQ/wPADRk7Qk0WuQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x14b873b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# clthresh/cuthresh example:\n",
    "tl12000_cd=clthresh(cd,12000)\n",
    "tu12000_cd=cuthresh(cd,12000)\n",
    "plt.hist(cd[np.nonzero(cd)], 32);\n",
    "plt.hist(tl12000_cd[np.nonzero(tl12000_cd)], 32);\n",
    "plt.hist(tu12000_cd[np.nonzero(tu12000_cd)], 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def csave(filename,cdata,chead=None,cnhdr=None,cxtra=None,cfmap=None):\n",
    "    # saves a cifti data structure to file. Output likely unusable without some header info\n",
    "    cout=nb.Cifti2Image(dataobj=cdata,header=chead,nifti_header=cnhdr, extra=cxtra, file_map=cfmap)\n",
    "    cout.to_filename(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# csave example:\n",
    "output_cifti=nhw+'/save_demo.dtseries.nii'\n",
    "csave(output_cifti,cd,ch,cn,ch,cf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cmask(cdata1,cdata2):\n",
    "    # masks cifti data 1 with non-zero values of cifti data 2\n",
    "    # cdata2 can be a single timepoint or matched to cdata1 dimensions\n",
    "    masked_cdata1=cdata1.copy()\n",
    "    if cdata1.shape != cdata2.shape:\n",
    "        mask=np.tile(cdata2,(cdata1.shape[0],1))\n",
    "    else:\n",
    "        mask=cdata2.copy()\n",
    "    masked_cdata1[mask==0]=0\n",
    "    return masked_cdata1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEDCAYAAADOc0QpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEVFJREFUeJzt3WusZWV9x/Hvr4M4acUrU0MHcIYGFVK8niJpjJca6wBN\npiY2GTDeipnQgrEvmjCNqZL4ZqyxrUZwMqUEaaoYldZpGaW1tsXEogwNt4GAI0xlRivj/RZrR/59\nsdfA5sycOeucs87Z5zzn+0l2Zq+1nr33f61sfjznWWs9O1WFJKktvzTpAiRJwzPcJalBhrskNchw\nl6QGGe6S1CDDXZIaNNFwT3JtkkeS3NOj7V8muaN7PJDk+0tRoyStRJnkde5JXgH8GLi+qn5jDq97\nB/DiqvqDRStOklawifbcq+oW4Lvj65L8epLPJbk9yReTPP8YL70I+PiSFClJK9AJky7gGHYCl1bV\nV5O8DLga+O0jG5M8B9gIfGFC9UnSsreswj3JU4DfAj6Z5MjqJ09rtgX4VFX9Yilrk6SVZFmFO6Nh\nou9X1YuO02YLcNkS1SNJK9KyuhSyqn4IPJTk9wEy8sIj27vx92cA/zmhEiVpRZj0pZAfZxTUz0ty\nIMklwBuBS5LcCewFNo+9ZAtwQzmVpSQd10QvhZQkLY5lNSwjSRrGxE6onnzyybVhw4ZJfbwkrUi3\n3377t6tq3WztJhbuGzZsYM+ePZP6eElakZL8d592DstIUoMMd0lq0KzhPtvMjUnemOSuJHcn+dL4\ndemSpMno03O/Dth0nO0PAa+sqnOA9zKaG0aSNEGznlCtqluSbDjO9i+NLd4KnLrwsiRJCzH0mPsl\nwGdn2phka5I9SfYcOnRo4I+WJB0xWLgneTWjcL9ipjZVtbOqpqpqat26WS/TlCTN0yDXuSd5AXAN\ncH5VfWeI95Qkzd+Ce+5JTgduBN5UVQ8svCRJ0kLN2nPvZm58FXBykgPAe4AnAVTVDuDdwLOAq7sf\n2DhcVVOLVbDasmHbTbO22b/9wiWoRGpLn6tlLppl+9uBtw9WkSRpwbxDVZIaZLhLUoMMd0lqkOEu\nSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGD/FiHNF2fqXwlLR57\n7pLUIMNdkhpkuEtSgxxz17LXd/zen+OTHme4a1XxfxRaLRyWkaQGGe6S1CDDXZIaZLhLUoM8oSod\nQ58Tr5501XJmz12SGmS4S1KDZg33JNcmeSTJPTNsT5IPJdmX5K4kLxm+TEnSXPTpuV8HbDrO9vOB\nM7vHVuAjCy9LkrQQs55Qrapbkmw4TpPNwPVVVcCtSZ6e5JSq+uZANUq9OM2w9LghxtzXAw+PLR/o\n1h0lydYke5LsOXTo0AAfLUk6liU9oVpVO6tqqqqm1q1bt5QfLUmryhDhfhA4bWz51G6dJGlChgj3\nXcCbu6tmzgN+4Hi7JE3WrCdUk3wceBVwcpIDwHuAJwFU1Q5gN3ABsA/4KfC2xSpWktRPn6tlLppl\newGXDVaRJGnBnFtGmid/+EPLmdMPSFKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpk\nuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoP8JSYtuf1rL55x24affWwJ\nK5HaZc9dkhpkuEtSgxyW0bLikI00DHvuktQgw12SGmS4S1KDDHdJalCvcE+yKcn9SfYl2XaM7U9L\n8o9J7kyyN8nbhi9VktTXrOGeZA1wFXA+cDZwUZKzpzW7DLi3ql4IvAr4QJITB65VktRTn0shzwX2\nVdWDAEluADYD9461KeCkJAGeAnwXODxwrdKKtGHbTbO22b/9wiWoRKtJn2GZ9cDDY8sHunXjPgyc\nBXwDuBt4Z1U9Ov2NkmxNsifJnkOHDs2zZEnSbIY6ofo64A7g14AXAR9O8tTpjapqZ1VNVdXUunXr\nBvpoSdJ0fcL9IHDa2PKp3bpxbwNurJF9wEPA84cpUZI0V33G3G8DzkyykVGobwGm3yP+deA1wBeT\nPBt4HvDgkIVqZTneNAKSFt+s4V5Vh5NcDtwMrAGuraq9SS7ttu8A3gtcl+RuIMAVVfXtRaxbq5Dz\nzkj99Zo4rKp2A7unrdsx9vwbwO8MW5okab68Q1WSGmS4S1KDnM9dzXOsXquRPXdJapDhLkkNMtwl\nqUGGuyQ1yBOqasJ874j1ZKtaZc9dkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDD\nXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIKX+lZWDDtpt6tdu//cJFrkStsOcuSQ0y3CWpQQ7L\naKLO2Xj6UevufujrE6hEakuvnnuSTUnuT7IvybYZ2rwqyR1J9ib5j2HLlCTNxaw99yRrgKuA1wIH\ngNuS7Kqqe8faPB24GthUVV9P8quLVbAkaXZ9eu7nAvuq6sGq+jlwA7B5WpuLgRur6usAVfXIsGVK\nkuaiz5j7euDhseUDwMumtXku8KQk/w6cBHywqq6f/kZJtgJbAU4//eixVrXtWOPrfds5Di/NzVBX\ny5wAvBS4EHgd8GdJnju9UVXtrKqpqppat27dQB8tSZquT8/9IHDa2PKp3bpxB4DvVNVPgJ8kuQV4\nIfDAIFVqWelzw83+tUtQiKQZ9em53wacmWRjkhOBLcCuaW0+A7w8yQlJfpnRsM19w5YqSepr1p57\nVR1OcjlwM7AGuLaq9ia5tNu+o6ruS/I54C7gUeCaqrpnMQuXJM2s101MVbUb2D1t3Y5py+8H3j9c\naZKk+XL6AUlqkNMPaEXw8khpbgx3rUp9rrk/iW386L7tS1CNNDzDXYui7w1LkhaH4a7m+T8arUaG\nu1aspRiHP+msoydBdahGK4FXy0hSgwx3SWqQ4S5JDXLMXU3x5Kk0Ys9dkhpkuEtSgxyW0bztX3vx\njNvOweERaZLsuUtSgwx3SWqQ4S5JDTLcpTk63rkGabkw3CWpQYa7JDXISyGlFWTDtptmbbN/+4VL\nUImWO3vuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUG9wj3JpiT3J9mX5OhfDH683W8mOZzkDcOV\nKEmaq1nDPcka4CrgfOBs4KIkZ8/Q7n3APw9dpCRpbvr03M8F9lXVg1X1c+AGYPMx2r0D+DTwyID1\nSZLmoU+4rwceHls+0K17TJL1wOuBjwxXmiRpvoY6ofpXwBVV9ejxGiXZmmRPkj2HDh0a6KMlSdP1\nmVvmIHDa2PKp3bpxU8ANSQBOBi5Icriq/mG8UVXtBHYCTE1N1XyL1vJyzkZ/Uk9abvqE+23AmUk2\nMgr1LcATJrSuqo1Hnie5Dvin6cEuSVo6s4Z7VR1OcjlwM7AGuLaq9ia5tNu+Y5FrlCTNUa8pf6tq\nN7B72rpjhnpVvXXhZUmSFsI7VCWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwl\nqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa\ndMKkC5A0rA3bburVbv/2Cxe5Ek2SPXdJapDhLkkNMtwlqUGGuyQ1yHCXpAb1ulomySbgg8Aa4Jqq\n2j5t+xuBK4AAPwL+sKruHLhWLQPnfPQcTjqre87pky1G0oxmDfcka4CrgNcCB4DbkuyqqnvHmj0E\nvLKqvpfkfGAn8LLFKFhL7MqnPXF5o4EurQR9hmXOBfZV1YNV9XPgBmDzeIOq+lJVfa9bvBU4ddgy\nJUlz0Sfc1wMPjy0f6NbN5BLgs8fakGRrkj1J9hw6dKh/lZKkORn0hGqSVzMK9yuOtb2qdlbVVFVN\nrVu3bsiPliSN6XNC9SBw2tjyqd26J0jyAuAa4Pyq+s4w5WkpHeu29f1rJ1CIpAXr03O/DTgzycYk\nJwJbgF3jDZKcDtwIvKmqHhi+TEnSXMzac6+qw0kuB25mdCnktVW1N8ml3fYdwLuBZwFXJwE4XFVT\ni1e2JOl4el3nXlW7gd3T1u0Ye/524O3DliZJmi/vUJWkBhnuktQgw12SGmS4S1KDDHdJapDhLkkN\nMtwlqUG9rnOX9ET7114847YNP/vYElYiHZs9d0lqkD13zeiks7b5a0sNO9ZEcdPt337hElSixWC4\nrxJ9/kOW1A7DXTOOH9trl1Yux9wlqUGGuyQ1yHCXpAY55r7C9T1RerzrsiW1x3AXAOds9OSp1BLD\nfZny0kVJC2G4N8ShF0lHGO4TMOleuUMwUvsM9wENFdpD9sANci1E7xP2TlOw7BjuDTHIJR1huE/I\nQnvnBrmk4zHce5rpz9PxkN6/duGfY2ivfM71ruXAcF8EBrSkSVv14d5rTuu1F8/YKzfIJS1HvcI9\nySbgg8Aa4Jqq2j5te7rtFwA/Bd5aVf81cK2Lqs8YuEEuHZs//LH8zBruSdYAVwGvBQ4AtyXZVVX3\njjU7Hzize7wM+Ej37/Jy5dOOWrV/7Si0nbtcS8HxeC2VPj33c4F9VfUgQJIbgM3AeLhvBq6vqgJu\nTfL0JKdU1TcHr3g2YwF+VE/bnreWsdaDf8ib9/wrYHZ9wn098PDY8gGO7pUfq8164AnhnmQrsLVb\n/HGS++dU7ROdDHz7+E3uWcDbrxg9jsOqsGTHIUvxIUf53b4NV8X3Ie/r1azVY/GcPo2W9IRqVe0E\ndg7xXkn2VNXUEO+1knkcRjwOIx6Hx632Y9HnxzoOAqeNLZ/arZtrG0nSEukT7rcBZybZmOREYAuw\na1qbXcCbM3Ie8IOJjLdLkoAewzJVdTjJ5cDNjC6FvLaq9ia5tNu+A9jN6DLIfYwuhXzb4pX8mEGG\ndxrgcRjxOIx4HB63qo9FRhe4SJJa4g9kS1KDDHdJatCKDPckm5Lcn2Rfkm2TrmdoSfYnuTvJHUn2\ndOuemeRfkny1+/cZY+3/tDsW9yd53dj6l3bvsy/Jh7ppIpa1JNcmeSTJPWPrBtv3JE9O8olu/ZeT\nbFjK/etrhuNwZZKD3ffijiQXjG1r7jgkOS3JvyW5N8neJO/s1q+678O8VNWKejA6qfs14AzgROBO\n4OxJ1zXwPu4HTp627s+Bbd3zbcD7uudnd8fgycDG7tis6bZ9BTiP0X03nwXOn/S+9dj3VwAvAe5Z\njH0H/gjY0T3fAnxi0vs8h+NwJfAnx2jb5HEATgFe0j0/CXig29dV932Yz2Ml9twfmw6hqn4OHJkO\noXWbgY92zz8K/N7Y+huq6n+r6iFGVyydm+QU4KlVdWuNvrnXj71m2aqqW4DvTls95L6Pv9engNcs\nx79oZjgOM2nyOFTVN6ubgLCqfgTcx+jO91X3fZiPlRjuM0110JICPp/k9m7KBoBn1+P3DvwP8Ozu\n+UzHY333fPr6lWjIfX/sNVV1GPgB8KzFKXtRvCPJXd2wzZHhiOaPQzdc8mLgy/h96GUlhvtq8PKq\nehGj2TYvS/KK8Y1d72NVXsO6mved0WyrZwAvYjRv0wcmW87SSPIU4NPAH1fVD8e3rfLvw3GtxHBv\nfqqDqjrY/fsI8PeMhqK+1f15SffvI13zmY7Hwe759PUr0ZD7/thrkpwAPA34zqJVPqCq+lZV/aKq\nHgX+mtH3Aho+DkmexCjY/66qbuxW+33oYSWGe5/pEFasJL+S5KQjz4HfYTS95S7gLV2ztwCf6Z7v\nArZ0Z/03MppT/yvdn60/THJeN4b45rHXrDRD7vv4e70B+ELX+1v2jgRa5/U8Pu1pk8ehq/lvgPuq\n6i/GNvl96GPSZ3Tn82A01cEDjM6Gv2vS9Qy8b2cwOuN/J7D3yP4xGgf8V+CrwOeBZ4695l3dsbif\nsStigClGAfA14MN0dyQv5wfwcUZDDv/HaGz0kiH3HVgLfJLRybavAGdMep/ncBz+FrgbuItRKJ3S\n8nEAXs5oyOUu4I7uccFq/D7M5+H0A5LUoJU4LCNJmoXhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNd\nkhr0/9KbvFXCDIsEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f8e85c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cmask example:\n",
    "masked3d=cmask(cd,tu12000_cd[0]) # masking 4d timeseries with single 3d timepoint\n",
    "masked4d=cmask(cd,tu12000_cd)    # masking 4d timeseries with 4d timeseries\n",
    "plt.hist(cd[np.nonzero(cd)], 32);\n",
    "plt.hist(masked3d[np.nonzero(masked3d)], 32);\n",
    "plt.hist(masked4d[np.nonzero(masked4d)], 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def oscale(olay):\n",
    "    # scales input array values to between zero and one, needed for iPyVolume colormaps\n",
    "    # for now, olay must be a 1D array (single timepoint)\n",
    "    scaled=np.asarray((olay-min(olay))/(max(olay)-min(olay)))\n",
    "    return scaled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mwmask(surf,mask):\n",
    "    # mask out medial wall; with some tinkering can be used to mask any vertices of (no) interest\n",
    "    # surf and mask should be .gii structures loaded with nibabel\n",
    "    # surf is the target surface\n",
    "    # mask is an overlay with no information at medial wall vertices (e.g. myelin map)\n",
    "    keep=(mask.darrays[0].data!=0)\n",
    "    kill=(mask.darrays[0].data==0)\n",
    "    ikeep=[i for i, e in enumerate(keep) if e != 0]\n",
    "    ikill=[i for i, e in enumerate(kill) if e != 0]\n",
    "    killdict={ii:1 for ii in ikill} # fun fact, iterating over a dictionary is ~exponentially faster vs. over a list\n",
    "    mask_kill=np.zeros([surf.darrays[1].data.shape[0]],dtype=bool) # create empty arrays matching surface mesh dimentions\n",
    "    mask_keep=mask_kill.copy()\n",
    "    for ii, row in enumerate(surf.darrays[1].data):\n",
    "        for item in row:\n",
    "            if item in killdict.keys():\n",
    "                mask_kill[ii]=True\n",
    "                continue\n",
    "            else:\n",
    "                mask_keep[ii]=True\n",
    "                continue\n",
    "    return mask_keep,mask_kill"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def showSurface(surface,overlay=None,frame=0,newfigure=True,colormap='summer',figsize=np.array([600,600]),\n",
    "                figlims=np.array([[-100,100],[-100,100],[-100,100]]),mwall=True):\n",
    "    '''\n",
    "    Displays a surface mesh in gifti or FreeSurfer (FS) surface format with/without an overlay inside  \n",
    "    Jupyter notebook for interactive visualization.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    surface: str, gifti opject\n",
    "        Path to surface file in gifti or FS surface format or an already loaded gifti object of surface\n",
    "    overlay: str, gifti opject\n",
    "        Path to overlay file in gifti or FS annot or anaotimcal (.curv,.sulc,.thickness) format or an already loaded \n",
    "        gifti object of overlay, default None\n",
    "    frame: int\n",
    "        indice of the frame (timepoint or functional data frame) to show\n",
    "    newfigure: bool\n",
    "        Create a new figure else prints into the last figure (in order to visualize both hemispheres in \n",
    "        one plot), default True\n",
    "    colormap: string\n",
    "        A matplotlib colormap, default summer\n",
    "    figsize: ndarray\n",
    "        Size of the figure to display, default [600,600]\n",
    "    figLims: ndarray\n",
    "        x,y and z limits of the axes, default [[-100,100],[-100,100],[-100,100]])\n",
    "\n",
    "    '''\n",
    "\n",
    "    if isinstance(surface,str):\n",
    "        if not os.path.exists(surface):\n",
    "            error('File does not exist, please provide a valid file path to a gifti or FreeSurfer file.')\n",
    "        filename, file_extension = os.path.splitext(surface)\n",
    "        if file_extension is '.gii':\n",
    "            surface = nb.load(surface)\n",
    "        else:\n",
    "            fsgeometry = nb.freesurfer.read_geometry(surface)\n",
    "            x,y,z = fsgeometry[0].T\n",
    "            vertex_edges=fsgeometry[1]\n",
    "\n",
    "    if isinstance(surface,nb.gifti.gifti.GiftiImage):\n",
    "        try:\n",
    "            vertex_spatial=surface.darrays[0]\n",
    "            vertex_edges=surface.darrays[1]\n",
    "            x, y, z = vertex_spatial.data.T\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid gifti file.')\n",
    "    \n",
    "    if mwall is False:\n",
    "        try:\n",
    "            mkeep,mkill=mwmask(surface,overlay)\n",
    "        except:\n",
    "            raise ValueError('Overlay required for medial wall masking.')\n",
    "    \n",
    "    if not isinstance(frame,int):\n",
    "        ValueError('Please provide a valid integer frame index.')\n",
    "    \n",
    "    if isinstance(overlay,list):\n",
    "        if frame>len(overlay) or frame < 0:\n",
    "            error('Frame index out of bounds, please provide a valid frame index.')\n",
    "        overlay = overlay[frame]\n",
    "    \n",
    "    if isinstance(overlay,str):\n",
    "        if not os.path.exists(overlay):\n",
    "            error('File does not exist, please provide a valid file path to a gifti or FreeSurfer file.')\n",
    "        filename, file_extension = os.path.splitext(overlay)\n",
    "        \n",
    "        if file_extension is '.gii':\n",
    "            overlay = nb.load(overlay)\n",
    "        elif (file_extension in ('.annot','')):\n",
    "            annot = nb.freesurfer.read_annot(overlay)\n",
    "            activation = annot[0]\n",
    "        elif (file_extension in ('.curv','.thickness','.sulc')):\n",
    "            activation = nb.freesurfer.read_morph_data(overlay)\n",
    "            \n",
    "    if isinstance(overlay,nb.gifti.gifti.GiftiImage):\n",
    "        try:\n",
    "            activation=overlay.darrays[0].data\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid gifti file')\n",
    "    \n",
    "    \n",
    "    if newfigure:\n",
    "\n",
    "        fig = p3.figure(width=figsize[0], height=figsize[1])\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 = (figlims[0][0], figlims[0][1])\n",
    "        fig.ylim = (figlims[1][0], figlims[1][1])\n",
    "        fig.zlim = (figlims[2][0], figlims[2][1])\n",
    "\n",
    "        \n",
    "    # plot surface\n",
    "    if overlay is None:\n",
    "        p3.plot_trisurf(x, y, z, triangles=vertex_edges.data)\n",
    "    elif mwall is False:\n",
    "        my_color = plt.cm.get_cmap(colormap)\n",
    "        colors=my_color((activation-min(activation))/(max(activation)-min(activation)))\n",
    "        p3.plot_trisurf(x, y, z, triangles=vertex_edges.data[mkeep], color=colors[:,:3])\n",
    "    else:\n",
    "        my_color = plt.cm.get_cmap(colormap)\n",
    "        colors=my_color((activation-min(activation))/(max(activation)-min(activation)))\n",
    "        p3.plot_trisurf(x, y, z, triangles=vertex_edges.data, color=colors[:,:3])\n",
    "\n",
    "    if newfigure:\n",
    "        p3.show()\n",
    "        \n",
    "    return 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4c2c58f77f97459290900061ef5db62b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4eb54b8cb1b7461dbe42016c57d045a7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# showSurface example:\n",
    "showSurface(data_rh) # right surface only\n",
    "showSurface(data_rh,data_rh_myelin) # new plot, same surface but with myelin map overlay\n",
    "showSurface(data_lh,data_lh_myelin,newfigure=False,mwall=False) # same plot, added \"very inflated\" left surface, excluded medial wall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cplot(coverlay,lsurf,rsurf,subvol,colormap):\n",
    "    # plots cifti overlay data with specified color on surfaces/subcortex\n",
    "    # lsurf and rsurf should be .surf.gii files matched to the number of vertices in each cifti hemisphere\n",
    "    # subvol should be a HCP-style subcortical volume (e.g. Atlas_ROIs.2.nii.gz)\n",
    "    \n",
    "    # import packages\n",
    "    import nibabel as nb\n",
    "    import matplotlib.pyplot as plt\n",
    "    import ipyvolume\n",
    "    import ipyvolume.pylab as p3\n",
    "    import numpy as np\n",
    "    \n",
    "    # read input files\n",
    "    cd,ch,cn,cx,cf=cread(coverlay) \n",
    "    ls=nb.load(lsurf)\n",
    "    rs=nb.load(rsurf)\n",
    "    sv=nb.load(subvol)\n",
    "    \n",
    "    # parse surface structure\n",
    "    lxyz=ls.darrays[0]\n",
    "    rxyz=rs.darrays[0]\n",
    "    lmesh=ls.darrays[1]\n",
    "    rmesh=rs.darrays[1]\n",
    "    \n",
    "    # parse cifti header\n",
    "    mim = ch.matrix[1]\n",
    "    bm1 = list(mim.brain_models)[0]\n",
    "    bm2 = list(mim.brain_models)[1]\n",
    "    lidx = list(bm1.vertex_indices)\n",
    "    ridx = [bm1.surface_number_of_vertices + val for val in bm2.vertex_indices]\n",
    "    bidx = np.concatenate((lidx, ridx))\n",
    "    # debug:\n",
    "    # print(bidx.shape)\n",
    "    #    for idx, bm in enumerate(mim.brain_models):\n",
    "    #        print((idx, bm.index_offset, bm.brain_structure))\n",
    "    \n",
    "    # split cifti overlay into left surface, right surface, and subcortical volume\n",
    "    clh=cd[:,0:len(lidx)]\n",
    "    crh=cd[:,len(lidx):len(lidx)+len(ridx)+1]\n",
    "    csc=cd[:,len(lidx)+len(ridx)+1:cd.shape[1]]\n",
    "    \n",
    "    # scale overlay values to between 0 and 1\n",
    "    slh=oscale(clh)\n",
    "    srh=oscale(crh)\n",
    "    scs=oscale(scs)\n",
    "    \n",
    "    # mask out medial wall vertices\n",
    "    lolay=np.zeros([lxyz.data.shape[0]])\n",
    "    rolay=np.zeros([rxyz.data.shape[0]])\n",
    "    lolay[lidx]=slh.copy()\n",
    "    rolay[ridx]=srh.copy()\n",
    "    \n",
    "    return lolay, rolay\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cplot(data_ciftiTS,data_lh,data_rh,data_sc,None)"
   ]
  }
 ],
 "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
}