ely_test_3d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ely's 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": 3,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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\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": 4,
   "metadata": {
    "scrolled": true
   },
   "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/ROIs/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": "code",
   "execution_count": 1540,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#important functions:\n",
    "\n",
    "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": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<b>3D Volume Rendering:</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "V = nb.load(nhw + '/data/100610/MNINonLinear/T1w_restore_brain.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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "img = nb.load(nhw + '/data/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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "solution2": "hidden"
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "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",
    "\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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "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": {},
   "source": [
    "\n",
    "\n",
    "<b>Load/Manipulate Ciftis:</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1543,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n"
     ]
    }
   ],
   "source": [
    "# see about importing the cifti more directly\n",
    "\n",
    "cifti_test = nb.load(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.ArealDistortion_MSMAll.32k_fs_LR.dscalar.nii')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1544,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n"
     ]
    }
   ],
   "source": [
    "#print(cifti_test.nifti_header.get_base_affine())\n",
    "#print(cifti_test.nifti_header.get_best_affine())\n",
    "#print(subcort.header.get_base_affine())\n",
    "#cdata=cifti_test.get_data()\n",
    "#cmat=cifti_test.header.matrix.get_index_map(1)\n",
    "#asdf=nb.cifti2.Cifti2Matrix(cifti_test,0)\n",
    "import cifti\n",
    "alt_cifti_test,aa = cifti.read(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.ArealDistortion_MSMAll.32k_fs_LR.dscalar.nii')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1545,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "memmap([[ 0.10375989,  0.02730756,  0.0826102 , ..., -0.68779516,\n",
       "        -0.59200358, -0.56638598]], dtype=float32)"
      ]
     },
     "execution_count": 1545,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alt_cifti_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'fmri_cifti_in' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-2883c0bc3d1c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mthalamus_mask\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkeep_L\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0maffine\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mthal_L\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maffine\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mfMRI_data\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfmri_cifti_in\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0mbm_cortex\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcifti\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mBrainModel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_surface\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m32492\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnvertex\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m32492\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbrain_structure\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'cortex_left'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mbm_thal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcifti\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mBrainModel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_mask\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthalamus_mask\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0many\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maffine\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'thalamus_left'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'fmri_cifti_in' is not defined"
     ]
    }
   ],
   "source": [
    "thal_L=nb.load('/Users/ely/templates/HCP/Subcortical_ROIs/Thalamus_L.nii.gz')\n",
    "keep_L=np.where(thal_L.get_data()!=0)\n",
    "thalamus_mask=np.asarray(keep_L)\n",
    "affine=thal_L.affine\n",
    "fMRI_data=fmri_cifti_in\n",
    "bm_cortex = cifti.BrainModel.from_surface(range(32492), nvertex=32492, brain_structure='cortex_left')\n",
    "bm_thal = cifti.BrainModel.from_mask(thalamus_mask.any(), affine, 'thalamus_left')\n",
    "bm_full = bm_cortex + bm_thal\n",
    "assert bm_full[bm_full.struc == 'CortexLeft'] == bm_cortex\n",
    "assert bm_full[bm_full.struc == 'ThalamusLeft'] == bm_thal\n",
    "parcel = cifti.Parcels([('surface_parcel', bm_cortex[:100]),('volume_parcel', bm_thal)]) #,('combined_parcel', bm_full[0,1])])\n",
    "#series = cifti.Series(start=0, step=100, size=5000)\n",
    "\n",
    "cifti_out=nhw+'/rfMRI1_LR.dtseries.nii'\n",
    "cifti_husk=[series,bm_full]\n",
    "cifti.write(cifti_out, fMRI_data, aa)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n",
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n"
     ]
    }
   ],
   "source": [
    "fmri_cifti='/Users/ely/Desktop/HCP/100610/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii'\n",
    "img = nb.cifti2.Cifti2Image.from_filename(fmri_cifti)\n",
    "img_load = nb.load(fmri_cifti)\n",
    "#%whos\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def cread(file):\n",
    "    # loads a cifti file and stores the data and metadata\n",
    "    cin=nb.load(file)\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": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n"
     ]
    }
   ],
   "source": [
    "cd,ch,cn,cx,cf=cread(fmri_cifti)\n",
    "#%whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1550,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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\n",
    "\n",
    "# clthresh example:\n",
    "# tl12000_cd=clthresh(cd,12000)\n",
    "# plt.hist(cd[np.nonzero(cd)], 32);\n",
    "# plt.hist(tl12000_cd[np.nonzero(tl12000_cd)], 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1551,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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\n",
    "\n",
    "# cuthresh example:\n",
    "# tu12000_cd=cuthresh(cd,12000)\n",
    "# plt.hist(cd[np.nonzero(cd)], 32);\n",
    "# plt.hist(tu12000_cd[np.nonzero(tu12000_cd)], 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "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": 1553,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1200, 91282)\n",
      "(1200, 91282)\n"
     ]
    }
   ],
   "source": [
    "print(cd.shape)\n",
    "print(tl12000_cd.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1554,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22970.512"
      ]
     },
     "execution_count": 1554,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(cd[np.nonzero(cd)])\n",
    "max(cd[np.nonzero(cd)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1557,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fmri_cifti_data=np.asarray(img.get_data()).copy()\n",
    "fmri_cifti_head=img.header.copy()\n",
    "fmri_cifti_nhead=img.nifti_header.copy()\n",
    "fmri_cifti_extra=img.extra.copy()\n",
    "fmri_cifti_fmap=img.file_map.copy()\n",
    "#out_img = nb.cifti2.Cifti2Image(dataobj=fmri_cifti_data, header=fmri_cifti_head,nifti_header=fmri_cifti_nhead, extra=fmri_cifti_extra, file_map=fmri_cifti_fmap)\n",
    "#cifti_out=nhw+'/custom_rfMRI1_LR.dtseries.nii'\n",
    "#img.to_filename(cifti_out)\n",
    "#%whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1558,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#thresh_data=fmri_cifti_data[fmri_cifti_data>10000]\n",
    "thresh_data=fmri_cifti_data.copy()\n",
    "thresh_data[thresh_data<10000]=0\n",
    "#a[a > 10] = 0\n",
    "#thresh_data = [0 if fmri_cifti_data_ > 10000 for fmri_cifti_data_ in fmri_cifti_data]\n",
    "out_img = nb.cifti2.Cifti2Image(dataobj=thresh_data, header=fmri_cifti_head,nifti_header=fmri_cifti_nhead, extra=fmri_cifti_extra, file_map=fmri_cifti_fmap)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.hist(fmri_cifti_data[np.nonzero(fmri_cifti_data)], 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1560,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "stufftosave=[tu12000_cd,ch,cn,cx,cf]\n",
    "out_img = nb.cifti2.Cifti2Image(dataobj=cd,header=ch,nifti_header=cn, extra=cx, file_map=cf)\n",
    "cout=nb.Cifti2Image(dataobj=cd,header=ch,nifti_header=cn, extra=cx, file_map=cf)\n",
    "#nb.save()\n",
    "\n",
    "#%whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1561,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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)\n",
    "#cifti_out=nhw+'/'+thresh_rfMRI1_LR.dtseries.nii'\n",
    "#out_img.to_filename(cifti_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1562,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "csave(nhw+'/test_csave_data_nhdr.dtseries.nii',cd,None,cn)\n",
    "#post_thresh= nb.cifti2.Cifti2Image.from_filename(cifti_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1564,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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",
    "#        return mask.shape\n",
    "    else:\n",
    "        mask=cdata2.copy()\n",
    "#        return mask.shape\n",
    "\n",
    "    masked_cdata1[mask==0]=0\n",
    "    return masked_cdata1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1565,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1200, 91282)"
      ]
     },
     "execution_count": 1565,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.tile(cd[0],(cd.shape[0],1)).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1566,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "masked_cd=cmask(cd,tu12000_cd[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.hist(cd[np.nonzero(cd)], 32);\n",
    "plt.hist(masked_cd[np.nonzero(masked_cd)], 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "csave(nhw+'/test_cmask_cd_singlevol.dtseries.nii',masked_cd,ch,cn,cx,cf)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "post_thresh_data=np.asarray(post_thresh.get_data())\n",
    "\n",
    "plt.hist(post_thresh_data[np.nonzero(post_thresh_data)], 32);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "thresh_data[thresh_data<20000]=0\n",
    "thresh_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#%whos\n",
    "#fMRI_data\n",
    "\n",
    "fmri_cifti_in,aa = cifti.read('/Users/ely/Desktop/HCP/100610/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "aa[0]\n",
    "def write(filename, arr, axes):\n",
    "    img = cifti2.Cifti2Image(np.asarray(arr), axis.to_header(axes))\n",
    "    img.to_filename(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cin=nb.load(fmri_cifti)\n",
    "cdata=cin.get_data()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ktb=ch.get_index_map(1)\n",
    "ktb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ktb.indices_map_to_data_type\n",
    "#nb.cifti2.Cifti2VertexIndices(cdata)\n",
    "fmri_cifti_load=nb.cifti2.load(fmri_cifti)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1569,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<property at 0x117307188>"
      ]
     },
     "execution_count": 1569,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bm=nb.cifti2.cifti2.Cifti2BrainModel(fmri_cifti,ch,cn,cx,cf)\n",
    "#cf[x]\n",
    "nb.fileholders.FileHolder.file_like\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for x in cf:\n",
    "    print (x)\n",
    "    for y in cf[x]:\n",
    "        print (y,':',cf[x][y])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 507,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "b'<BrainModel BrainStructure=\"{}\" IndexCount=\"&lt;nibabel.cifti2.cifti2.Cifti2Header object at 0x127898860&gt;\" IndexOffset=\"[[ 10370.02734375   8975.28808594   8177.61816406 ...,   8267.63867188&#10;    9081.58300781   8428.43261719]&#10; [ 10405.33300781   8885.35449219   8228.50195312 ...,   9814.16601562&#10;   10344.37695312   9132.24707031]&#10; [ 10403.03027344   8885.85253906   8170.07861328 ...,   9360.21777344&#10;    8852.38671875   9805.09863281]&#10; ..., &#10; [ 10651.48046875   8867.24511719   8167.99414062 ...,   8364.47460938&#10;   11047.29199219   9667.34179688]&#10; [ 10547.875        8944.50683594   8189.71582031 ...,  10087.88476562&#10;    9486.22949219   9620.84765625]&#10; [ 10372.42089844   8968.25292969   8143.32861328 ...,   9411.23144531&#10;   10986.9375       8809.45214844]]\" ModelType=\"&lt;class \\'nibabel.nifti2.Nifti2Header\\'&gt; object, endian=\\'&lt;\\'&#10;sizeof_hdr      : 540&#10;magic           : b\\'n+2\\'&#10;eol_check       : [13 10 26 10]&#10;datatype        : float32&#10;bitpix          : 32&#10;dim             : [    2  1200 91282     1     1     1     1     1]&#10;intent_p1       : 0.0&#10;intent_p2       : 0.0&#10;intent_p3       : 0.0&#10;pixdim          : [ 1.  1.  1.  1.  1.  1.  1.  1.]&#10;vox_offset      : 0&#10;scl_slope       : nan&#10;scl_inter       : nan&#10;cal_max         : 0.0&#10;cal_min         : 0.0&#10;slice_duration  : 0.0&#10;toffset         : 0.0&#10;slice_start     : 0&#10;slice_end       : 0&#10;descrip         : b\\'\\'&#10;aux_file        : b\\'\\'&#10;qform_code      : unknown&#10;sform_code      : unknown&#10;quatern_b       : 0.0&#10;quatern_c       : 0.0&#10;quatern_d       : 0.0&#10;qoffset_x       : 0.0&#10;qoffset_y       : 0.0&#10;qoffset_z       : 0.0&#10;srow_x          : [ 0.  0.  0.  0.]&#10;srow_y          : [ 0.  0.  0.  0.]&#10;srow_z          : [ 0.  0.  0.  0.]&#10;slice_code      : unknown&#10;xyzt_units      : 10&#10;intent_code     : dense data series/fiber fans&#10;intent_name     : b\\'ConnDenseSeries\\'&#10;dim_info        : 0&#10;unused_str      : b\\'\\'\" SurfaceNumberOfVertices=\"{\\'image\\': &lt;nibabel.fileholders.FileHolder object at 0x1251d66d8&gt;}\" />'"
      ]
     },
     "execution_count": 507,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bm.to_xml()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1570,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cifti data shape = (1200, 91282)\n",
      "\n",
      "cifti data dype = fMRI timeseries\n",
      "Xdim = 2.0x10^-3m\n",
      "     1200 timepoints, TR = 0.72 SECOND\n",
      "(0, 0, 29696, 'CIFTI_STRUCTURE_CORTEX_LEFT')\n",
      "(1, 29696, 29716, 'CIFTI_STRUCTURE_CORTEX_RIGHT')\n",
      "(2, 59412, 135, 'CIFTI_STRUCTURE_ACCUMBENS_LEFT')\n",
      "(3, 59547, 140, 'CIFTI_STRUCTURE_ACCUMBENS_RIGHT')\n",
      "(4, 59687, 315, 'CIFTI_STRUCTURE_AMYGDALA_LEFT')\n",
      "(5, 60002, 332, 'CIFTI_STRUCTURE_AMYGDALA_RIGHT')\n",
      "(6, 60334, 3472, 'CIFTI_STRUCTURE_BRAIN_STEM')\n",
      "(7, 63806, 728, 'CIFTI_STRUCTURE_CAUDATE_LEFT')\n",
      "(8, 64534, 755, 'CIFTI_STRUCTURE_CAUDATE_RIGHT')\n",
      "(9, 65289, 8709, 'CIFTI_STRUCTURE_CEREBELLUM_LEFT')\n",
      "(10, 73998, 9144, 'CIFTI_STRUCTURE_CEREBELLUM_RIGHT')\n",
      "(11, 83142, 706, 'CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_LEFT')\n",
      "(12, 83848, 712, 'CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_RIGHT')\n",
      "(13, 84560, 764, 'CIFTI_STRUCTURE_HIPPOCAMPUS_LEFT')\n",
      "(14, 85324, 795, 'CIFTI_STRUCTURE_HIPPOCAMPUS_RIGHT')\n",
      "(15, 86119, 297, 'CIFTI_STRUCTURE_PALLIDUM_LEFT')\n",
      "(16, 86416, 260, 'CIFTI_STRUCTURE_PALLIDUM_RIGHT')\n",
      "(17, 86676, 1060, 'CIFTI_STRUCTURE_PUTAMEN_LEFT')\n",
      "(18, 87736, 1010, 'CIFTI_STRUCTURE_PUTAMEN_RIGHT')\n",
      "(19, 88746, 1288, 'CIFTI_STRUCTURE_THALAMUS_LEFT')\n",
      "(20, 90034, 1248, 'CIFTI_STRUCTURE_THALAMUS_RIGHT')\n",
      "Voxel to MNI transformation matrix:\n",
      "[[  -2.    0.    0.   90.]\n",
      " [   0.    2.    0. -126.]\n",
      " [   0.    0.    2.  -72.]\n",
      " [   0.    0.    0.    1.]]\n"
     ]
    }
   ],
   "source": [
    "def cinfo(infile):\n",
    "    # returns information about input cifti file, particularly header metadata\n",
    "    cd,ch,cn,cx,cf=cread(infile)\n",
    "    ch0=ch.get_index_map(0)\n",
    "    ch1=ch.get_index_map(1)\n",
    "    # useful outputs:\n",
    "    print('cifti data shape = {x}'.format(x=cd.shape))\n",
    "    print('')\n",
    "    if ch.get_index_map(0).indices_map_to_data_type=='CIFTI_INDEX_TYPE_SERIES':\n",
    "        print('cifti data dype = fMRI timeseries')\n",
    "        me=ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.meter_exponent\n",
    "        print('Xdim = {x}x10^{e}m'.format(\n",
    "            x=np.abs(ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix[0,0]),\n",
    "            e=me))\n",
    "        \n",
    "        print('     {x} timepoints, TR = {y} {z}'.format(\n",
    "            x=ch.get_index_map(0).number_of_series_points,\n",
    "            y=ch.get_index_map(0).series_step,\n",
    "            z=ch.get_index_map(0).series_unit))\n",
    "\n",
    "    for idx, bm in enumerate(ch1.brain_models):\n",
    "        print((idx, bm.index_offset, bm.index_count, bm.brain_structure)) # list brain structures, starting indices, and sizes\n",
    "    print('Voxel to MNI transformation matrix:')\n",
    "    print(ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix)\n",
    "    \n",
    "    # not returned:\n",
    "#    ch.may_contain_header() # unclear how to use\n",
    "#    bm.model_type # returns 'CIFTI_MODEL_TYPE_SURFACE' or 'CIFTI_MODEL_TYPE_VOXELS'\n",
    "#    bm.surface_number_of_vertices # returns total vertices in surface (includes medial wall)\n",
    "#    print('cifti header number of mapped indices = {x}'.format(x=ch.number_of_mapped_indices))\n",
    "#    print('cifti header mapped indices = {x}'.format(x=ch.mapped_indices))\n",
    "\n",
    "cinfo(fmri_cifti)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1577,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "91"
      ]
     },
     "execution_count": 1577,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mim=ch.get_index_map(1)\n",
    "mim.applies_to_matrix_dimension\n",
    "mim.indices_map_to_data_type\n",
    "mim.number_of_series_points\n",
    "nimm=np.asarray(mim.named_maps)\n",
    "mim.volume.transformation_matrix_voxel_indices_ijk_to_xyz.meter_exponent\n",
    "ch1=ch.get_index_map(1)\n",
    "warp=ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix\n",
    "ch1.volume.volume_dimensions[0]\n",
    "#ch1.to_xml()\n",
    "#for idx, nm in enumerate(mim.named_maps):\n",
    "#    print((idx, nm.map_name, nm.metadata, nm.label_table))\n",
    "#mim.number_of_series_points\n",
    "#mim.series_exponent\n",
    "#mim.series_unit\n",
    "#mim.volume.volume_dimensions[1]\n",
    "#np.zeros((mim.volume.volume_dimensions)).shape\n",
    "#cn.keys()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def cvoxels(infile,xyz=False):\n",
    "    # returns array where each element is an array of voxel coordinates for a particular structure\n",
    "    cd,ch,cn,cx,cf=cread(infile)\n",
    "    ch1=ch.get_index_map(1)\n",
    "    voxarray=[]\n",
    "    if xyz is True:\n",
    "        from nibabel.affines import apply_affine\n",
    "        warp=ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix\n",
    "    for idx, bm in enumerate(ch1.brain_models):\n",
    "        if bm.model_type=='CIFTI_MODEL_TYPE_VOXELS':\n",
    "            if xyz is False:\n",
    "                voxarray.append(bm.voxel_indices_ijk)\n",
    "            else:\n",
    "                voxarray.append(apply_affine(warp,bm.voxel_indices_ijk))\n",
    "        else:\n",
    "            voxarray.append([[],[],[]])\n",
    "    return np.asarray(voxarray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1428,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def cijk2xyz(infile):\n",
    "    # returns array of MNI X/Y/Z values for subcortical structures based on I/J/K voxel indices\n",
    "    from nibabel.affines import apply_affine\n",
    "    cd,ch,cn,cx,cf=cread(infile)\n",
    "    ch1=ch.get_index_map(1)\n",
    "    warp=ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix\n",
    "    ijk=cvoxels(infile)\n",
    "    xyz=[]\n",
    "    for structure in range(0,len(ijk)):\n",
    "        xyz.append(np.asarray(apply_affine(warp,ijk[structure])))\n",
    "    return xyz\n",
    "\n",
    "#xyz=np.zeros((mim.volume.volume_dimensions), dtype = object)\n",
    "#for xx in range (0, subcort_disp.shape[0]):\n",
    "#    for yy in range (0, subcort_disp.shape[1]):\n",
    "#        for zz in range (0, subcort_disp.shape[2]):\n",
    "#            real_coords[xx,yy,zz]=apply_affine(subcort.affine,[xx,yy,zz])\n",
    "#            print(real_coords[xx,yy,zz])\n",
    "#            xyz_sub[i]d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "pixdim[1,2,3] should be non-zero; setting 0 dims to 1\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([list([[], [], []]), list([[], [], []]),\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520d68>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520be0>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520a20>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520860>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c5206a0>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c5204e0>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520320>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520160>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520da0>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520c18>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520a58>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520898>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c5206d8>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520518>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520358>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520198>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c520208>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c8d5550>,\n",
       "       <nibabel.cifti2.cifti2.Cifti2VoxelIndicesIJK object at 0x11c8d5320>], dtype=object)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#from nibabel.affines import apply_affine\n",
    "#warp=mim.volume.transformation_matrix_voxel_indices_ijk_to_xyz.\n",
    "#xyz=np.zeros((mim.volume.volume_dimensions), dtype = object)\n",
    "#real_coords2=apply_affine(warp,varray[0])\n",
    "#vxyz1=cijk2xyz(fmri_cifti)\n",
    "#vijk=cvoxels(fmri_cifti)\n",
    "cvoxels(fmri_cifti,xyz=False)\n",
    "#varray[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(21,)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#a=ch.matrix\n",
    "#b=ch.get_index_map(1)\n",
    "#ch.number_of_mapped_indices\n",
    "#a.number_of_series_points\n",
    "#ab=a.brain_models\n",
    "#enumerate(ab)\n",
    "#whos\n",
    "#cm0=ch.get_index_map(0)\n",
    "#cm1=ch.get_index_map(1)\n",
    "#np.asarray(list(ch.get_index_map(1).brain_models)[3].voxel_indices_ijk)\n",
    "#mim.volume.transformation_matrix_voxel_indices_ijk_to_xyz.meter_exponent\n",
    "vxyz.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1477,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variable           Type                      Data/Info\n",
      "------------------------------------------------------\n",
      "Poly3DCollection   type                      <class 'mpl_toolkits.mplo<...>.art3d.Poly3DCollection'>\n",
      "a                  Cifti2Matrix              <nibabel.cifti2.cifti2.Ci<...>ix object at 0x129235a90>\n",
      "ab                 generator                 <generator object brain_models at 0x129c4e570>\n",
      "apply_affine       function                  <function apply_affine at 0x11730a840>\n",
      "b                  Cifti2MatrixIndicesMap    <nibabel.cifti2.cifti2.Ci<...>ap object at 0x129ed24e0>\n",
      "bb                 int                       91282\n",
      "bidx               ndarray                   59412: 59412 elems, type `int64`, 475296 bytes (464.15625 kb)\n",
      "bm                 Cifti2BrainModel          <nibabel.cifti2.cifti2.Ci<...>el object at 0x12fda45f8>\n",
      "bm1                Cifti2BrainModel          <nibabel.cifti2.cifti2.Ci<...>el object at 0x129ed2438>\n",
      "bm2                Cifti2BrainModel          <nibabel.cifti2.cifti2.Ci<...>el object at 0x129ed22b0>\n",
      "bm_cortex          BrainModel                <cifti.axis.BrainModel object at 0x11c95d668>\n",
      "bm_full            BrainModel                <cifti.axis.BrainModel object at 0x11de0f588>\n",
      "bmx                bytes                     b'<BrainModel IndexOffset<...>t at 0x1251a1390&gt;\" />'\n",
      "c                  int                       0\n",
      "cclh               ndarray                   1200x29696: 35635200 elems, type `float32`, 142540800 bytes (135.9375 Mb)\n",
      "ccrh               ndarray                   1200x29717: 35660400 elems, type `float32`, 142641600 bytes (136.03363037109375 Mb)\n",
      "ccsc               ndarray                   1200x31869: 38242800 elems, type `float32`, 152971200 bytes (145.88470458984375 Mb)\n",
      "cd                 ndarray                   1200x91282: 109538400 elems, type `float32`, 438153600 bytes (417.8558349609375 Mb)\n",
      "cf                 dict                      n=1\n",
      "ch                 Cifti2Header              <nibabel.cifti2.cifti2.Ci<...>er object at 0x1292350b8>\n",
      "cifti              module                    <module 'cifti' from '/Us<...>kages/cifti/__init__.py'>\n",
      "cifti_husk         list                      n=2\n",
      "cifti_out          str                       /Users/ely/Desktop/nhw201<...>00_rfMRI1_LR.dtseries.nii\n",
      "cifti_test         Cifti2Image               <nibabel.cifti2.cifti2.Ci<...>ge object at 0x110f7e8d0>\n",
      "cijk2xyz           function                  <function cijk2xyz at 0x5771338c8>\n",
      "cin                Cifti2Image               <nibabel.cifti2.cifti2.Ci<...>ge object at 0x1251a1390>\n",
      "cinfo              function                  <function cinfo at 0x126039d08>\n",
      "clh                ndarray                   1200x29696: 35635200 elems, type `float32`, 142540800 bytes (135.9375 Mb)\n",
      "clthresh           function                  <function clthresh at 0x127caad08>\n",
      "cm0                Cifti2MatrixIndicesMap    <nibabel.cifti2.cifti2.Ci<...>ap object at 0x129235898>\n",
      "cm1                Cifti2MatrixIndicesMap    <nibabel.cifti2.cifti2.Ci<...>ap object at 0x129ed24e0>\n",
      "cn                 Nifti2Header              <class 'nibabel.nifti2.Ni<...> 0\\nunused_str      : b''\n",
      "color_lh           ndarray                   32492x4: 129968 elems, type `float64`, 1039744 bytes (1015.375 kb)\n",
      "color_mlh          ndarray                   32492x4: 129968 elems, type `float64`, 1039744 bytes (1015.375 kb)\n",
      "color_rh           ndarray                   32492x4: 129968 elems, type `float64`, 1039744 bytes (1015.375 kb)\n",
      "cort_lh            ndarray                   32492: 32492 elems, type `bool`, 32492 bytes\n",
      "cort_rh            ndarray                   32492: 32492 elems, type `bool`, 32492 bytes\n",
      "cort_tri_lh        ndarray                   29696x64980: 1929646080 elems, type `bool`, 1929646080 bytes (1840.25390625 Mb)\n",
      "cort_tri_rh        ndarray                   29696x64980: 1929646080 elems, type `bool`, 1929646080 bytes (1840.25390625 Mb)\n",
      "cort_tri_v2        ndarray                   5833x3: 17499 elems, type `int32`, 69996 bytes\n",
      "cort_tri_v2l       list                      n=29696\n",
      "cort_tri_v3        ndarray                   59628x3: 178884 elems, type `int32`, 715536 bytes (698.765625 kb)\n",
      "count              int                       0\n",
      "cout               Cifti2Image               <nibabel.cifti2.cifti2.Ci<...>ge object at 0x12194a4e0>\n",
      "cout_img           Cifti2Image               <nibabel.cifti2.cifti2.Ci<...>ge object at 0x12780f588>\n",
      "cparse             function                  <function cparse at 0x25c759048>\n",
      "cread              function                  <function cread at 0x12b4f56a8>\n",
      "crh                ndarray                   1200x29717: 35660400 elems, type `float32`, 142641600 bytes (136.03363037109375 Mb)\n",
      "csc                ndarray                   1200x31869: 38242800 elems, type `float32`, 152971200 bytes (145.88470458984375 Mb)\n",
      "cuthresh           function                  <function cuthresh at 0x127caa2f0>\n",
      "cvoxels            function                  <function cvoxels at 0x25986ab70>\n",
      "cx                 dict                      n=0\n",
      "data_ciftiTS       Cifti2Image               <nibabel.cifti2.cifti2.Ci<...>ge object at 0x1212bc208>\n",
      "data_lh            GiftiImage                <nibabel.gifti.gifti.Gift<...>ge object at 0x12780f8d0>\n",
      "data_lh_myelin     GiftiImage                <nibabel.gifti.gifti.Gift<...>ge object at 0x11d15f860>\n",
      "data_rh            GiftiImage                <nibabel.gifti.gifti.Gift<...>ge object at 0x12780f898>\n",
      "data_rh_myelin     GiftiImage                <nibabel.gifti.gifti.Gift<...>ge object at 0x12780ff98>\n",
      "data_sc            Nifti1Image               <class 'nibabel.nifti1.Ni<...>nmagic           : b'n+1'\n",
      "echo               function                  <function echo at 0x136d66378>\n",
      "fig                Figure                    Figure(anglex=0.353035695<...>0], zlim=[-100.0, 100.0])\n",
      "fmri_cifti         str                       /Users/ely/Desktop/HCP/10<...>hp2000_clean.dtseries.nii\n",
      "fmri_cifti_in      memmap                    [[ 10370.02734375   8975.<...>375       8809.45214844]]\n",
      "i                  int                       32491\n",
      "idx                int                       20\n",
      "ii                 int                       64979\n",
      "iii                list                      n=29696\n",
      "img                Cifti2Image               <nibabel.cifti2.cifti2.Ci<...>ge object at 0x1207299e8>\n",
      "indices            list                      n=29696\n",
      "ipyvolume          module                    <module 'ipyvolume' from <...>s/ipyvolume/__init__.py'>\n",
      "item               int32                     21432\n",
      "keep_L             tuple                     n=3\n",
      "killdict           dict                      n=2796\n",
      "ktb                Cifti2MatrixIndicesMap    <nibabel.cifti2.cifti2.Ci<...>ap object at 0x1214f6e48>\n",
      "lidx               list                      n=29696\n",
      "m                  ndarray                   64980x3: 194940 elems, type `bool`, 194940 bytes (190.37109375 kb)\n",
      "md                 dict                      n=29696\n",
      "mesh_lh            GiftiDataArray            <nibabel.gifti.gifti.Gift<...>ay object at 0x12c3f25c0>\n",
      "mesh_rh            GiftiDataArray            <nibabel.gifti.gifti.Gift<...>ay object at 0x135947710>\n",
      "mim                Cifti2MatrixIndicesMap    <nibabel.cifti2.cifti2.Ci<...>ap object at 0x129ed24e0>\n",
      "mkeep              ndarray                   64980: 64980 elems, type `bool`, 64980 bytes\n",
      "mkill              ndarray                   64980: 64980 elems, type `bool`, 64980 bytes\n",
      "mlv                int                       29696\n",
      "mwmask             function                  <function mwmask at 0x12b51b400>\n",
      "nb                 module                    <module 'nibabel' from '/<...>bel/nibabel/__init__.py'>\n",
      "nhw                str                       /Users/ely/Desktop/nhw2017/NHW_visualization\n",
      "nim                generator                 <generator object named_maps at 0x13020aba0>\n",
      "nimm               ndarray                   : 1 elems, type `object`, 8 bytes\n",
      "np                 module                    <module 'numpy' from '/Us<...>kages/numpy/__init__.py'>\n",
      "os                 module                    <module 'os' from '/Users<...>nda/lib/python3.6/os.py'>\n",
      "out_img            Cifti2Image               <nibabel.cifti2.cifti2.Ci<...>ge object at 0x12194ac50>\n",
      "p3                 module                    <module 'ipyvolume.pylab'<...>ages/ipyvolume/pylab.py'>\n",
      "plt                module                    <module 'matplotlib.pyplo<...>es/matplotlib/pyplot.py'>\n",
      "post_thresh        Cifti2Image               <nibabel.cifti2.cifti2.Ci<...>ge object at 0x11ca01898>\n",
      "real_coords        ndarray                   91x109x91: 902629 elems, type `object`, 7221032 bytes (6.886512756347656 Mb)\n",
      "real_coords2       ndarray                   135x3: 405 elems, type `float64`, 3240 bytes\n",
      "ridx               list                      n=29716\n",
      "row                ndarray                   3: 3 elems, type `int32`, 12 bytes\n",
      "series             Series                    <cifti.axis.Series object at 0x11fd172e8>\n",
      "showSurface        function                  <function showSurface at 0x136d6b950>\n",
      "sub_data           ndarray                   91x109x91: 902629 elems, type `float32`, 3610516 bytes (3.443256378173828 Mb)\n",
      "sub_lh             Mesh                      Mesh(color=array([[ 0.124<...>5813141], dtype=float32))\n",
      "sub_rh             ndarray                   32492: 32492 elems, type `bool`, 32492 bytes\n",
      "sub_tri_lh         ndarray                   2796x64980: 181684080 elems, type `bool`, 181684080 bytes (173.26744079589844 Mb)\n",
      "sub_tri_rh         ndarray                   2776x64980: 180384480 elems, type `bool`, 180384480 bytes (172.02804565429688 Mb)\n",
      "subcort            Nifti1Image               <class 'nibabel.nifti1.Ni<...>nmagic           : b'n+1'\n",
      "subcort_disp       ndarray                   91x109x91: 902629 elems, type `bool`, 902629 bytes (881.4736328125 kb)\n",
      "surf_lh            Mesh                      Mesh(color=array([[ 0.124<...>5813141], dtype=float32))\n",
      "surf_mlh           Mesh                      Mesh(color=array([[ 1.   <...>2891235], dtype=float32))\n",
      "surf_rh            Mesh                      Mesh(color=array([[ 0.487<...>741272 ], dtype=float32))\n",
      "sys                module                    <module 'sys' (built-in)>\n",
      "t12000_cd          ndarray                   1200x91282: 109538400 elems, type `float32`, 438153600 bytes (417.8558349609375 Mb)\n",
      "tl12000_cd         ndarray                   1200x91282: 109538400 elems, type `float32`, 438153600 bytes (417.8558349609375 Mb)\n",
      "tu12000_cd         ndarray                   1200x91282: 109538400 elems, type `float32`, 438153600 bytes (417.8558349609375 Mb)\n",
      "uu                 Cifti2MatrixIndicesMap    <nibabel.cifti2.cifti2.Ci<...>ap object at 0x12778b198>\n",
      "v                  GiftiImage                <nibabel.gifti.gifti.Gift<...>ge object at 0x577119b70>\n",
      "v2_ii              list                      n=29696\n",
      "v3_ii              list                      n=2796\n",
      "varray             ndarray                   19: 19 elems, type `object`, 152 bytes\n",
      "vijk               ndarray                   21: 21 elems, type `object`, 168 bytes\n",
      "vv                 list                      n=1\n",
      "vxyz               ndarray                   21: 21 elems, type `object`, 168 bytes\n",
      "vxyz1              list                      n=19\n",
      "warp               ndarray                   4x4: 16 elems, type `float64`, 128 bytes\n",
      "x                  ndarray                   31870: 31870 elems, type `int64`, 254960 bytes (248.984375 kb)\n",
      "xx                 int                       90\n",
      "xyz                ndarray                   91x109x91: 902629 elems, type `object`, 7221032 bytes (6.886512756347656 Mb)\n",
      "xyz_lh             GiftiDataArray            <nibabel.gifti.gifti.Gift<...>ay object at 0x577116630>\n",
      "xyz_rh             GiftiDataArray            <nibabel.gifti.gifti.Gift<...>ay object at 0x57710c198>\n",
      "xyz_sub            list                      n=0\n",
      "y                  ndarray                   31870: 31870 elems, type `int64`, 254960 bytes (248.984375 kb)\n",
      "yy                 int                       108\n",
      "z                  ndarray                   31870: 31870 elems, type `int64`, 254960 bytes (248.984375 kb)\n",
      "zed                ndarray                   32492: 32492 elems, type `float64`, 259936 bytes (253.84375 kb)\n",
      "zz                 int                       90\n"
     ]
    }
   ],
   "source": [
    "%whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1239,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 0, 'CIFTI_STRUCTURE_CORTEX_LEFT')\n",
      "(1, 29696, 'CIFTI_STRUCTURE_CORTEX_RIGHT')\n",
      "(2, 59412, 'CIFTI_STRUCTURE_ACCUMBENS_LEFT')\n",
      "(3, 59547, 'CIFTI_STRUCTURE_ACCUMBENS_RIGHT')\n",
      "(4, 59687, 'CIFTI_STRUCTURE_AMYGDALA_LEFT')\n",
      "(5, 60002, 'CIFTI_STRUCTURE_AMYGDALA_RIGHT')\n",
      "(6, 60334, 'CIFTI_STRUCTURE_BRAIN_STEM')\n",
      "(7, 63806, 'CIFTI_STRUCTURE_CAUDATE_LEFT')\n",
      "(8, 64534, 'CIFTI_STRUCTURE_CAUDATE_RIGHT')\n",
      "(9, 65289, 'CIFTI_STRUCTURE_CEREBELLUM_LEFT')\n",
      "(10, 73998, 'CIFTI_STRUCTURE_CEREBELLUM_RIGHT')\n",
      "(11, 83142, 'CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_LEFT')\n",
      "(12, 83848, 'CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_RIGHT')\n",
      "(13, 84560, 'CIFTI_STRUCTURE_HIPPOCAMPUS_LEFT')\n",
      "(14, 85324, 'CIFTI_STRUCTURE_HIPPOCAMPUS_RIGHT')\n",
      "(15, 86119, 'CIFTI_STRUCTURE_PALLIDUM_LEFT')\n",
      "(16, 86416, 'CIFTI_STRUCTURE_PALLIDUM_RIGHT')\n",
      "(17, 86676, 'CIFTI_STRUCTURE_PUTAMEN_LEFT')\n",
      "(18, 87736, 'CIFTI_STRUCTURE_PUTAMEN_RIGHT')\n",
      "(19, 88746, 'CIFTI_STRUCTURE_THALAMUS_LEFT')\n",
      "(20, 90034, 'CIFTI_STRUCTURE_THALAMUS_RIGHT')\n",
      "(59412,)\n"
     ]
    }
   ],
   "source": [
    "mim = ch.matrix[1]\n",
    "for idx, bm in enumerate(mim.brain_models):\n",
    "    print((idx, bm.index_offset, bm.brain_structure))\n",
    "bm1 = list(mim.brain_models)[0]\n",
    "lidx = list(bm1.vertex_indices)\n",
    "bm2 = list(mim.brain_models)[1]\n",
    "ridx = [bm1.surface_number_of_vertices + val for val in bm2.vertex_indices]\n",
    "bidx = np.concatenate((lidx, ridx))\n",
    "print(bidx.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 694,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "29696"
      ]
     },
     "execution_count": 694,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(lidx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#nb.cifti2.Cifti2BrainModel\n",
    "#%whos \n",
    "data_lh_myelin.darrays[0].data.shape\n",
    "for x in data_lh_myelin.darrays[0].data:\n",
    "    \n",
    "    print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "\n",
    "\n",
    "<b> Surface Rendering: </b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1046,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# import surface mesh and subcortical volume\n",
    "#data_lh = nb.load(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.L.midthickness_MSMAll.32k_fs_LR.surf.gii')\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 = nb.load(nhw+'/data/100610/T1w/Native/100610.L.pial.native.surf.gii')\n",
    "#data_rh = nb.load(nhw+'/data/100610/T1w/Native/100610.R.pial.native.surf.gii')\n",
    "subcort = nb.load('/Users/ely/Desktop/HCP/100610/MNINonLinear/ROIs/Atlas_ROIs.2.nii.gz')\n",
    "\n",
    "xyz_lh=data_lh.darrays[0]\n",
    "xyz_rh=data_rh.darrays[0]\n",
    "mesh_lh=data_lh.darrays[1]\n",
    "mesh_rh=data_rh.darrays[1]\n",
    "cort_lh=(data_lh_myelin.darrays[0].data>0)\n",
    "cort_rh=(data_rh_myelin.darrays[0].data>0)\n",
    "sub_lh=(data_lh_myelin.darrays[0].data==0)\n",
    "sub_rh=(data_rh_myelin.darrays[0].data==0)\n",
    "\n",
    "\n",
    "#cort_olay_lh=data_lh_myelin.darrays[0].data[cort_lh]\n",
    "#cort_olay_rh=data_rh_myelin.darrays[0].data[cort_rh]\n",
    "#sub_olay_lh=data_lh_myelin.darrays[0].data[sub_lh]\n",
    "#sub_olay_rh=data_rh_myelin.darrays[0].data[sub_rh]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1151,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(32492,)"
      ]
     },
     "execution_count": 1151,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_lh_myelin.darrays[0].data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1047,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# import data to overlay\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",
    "olay_lh=data_lh_myelin.darrays[0].data\n",
    "olay_rh=data_rh_myelin.darrays[0].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 991,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# loop through xyz values to create lists of cortical and subcortical triangles\n",
    "# left\n",
    "cort_tri_lh=[]\n",
    "for ii in np.where(cort_lh)[0]:\n",
    "    cort_mask_lha=mesh_lh.data==ii\n",
    "    cort_mask_lh=cort_mask_lha.any(1)\n",
    "    cort_tri_lh.append(cort_mask_lh)\n",
    "    \n",
    "sub_tri_lh=[]\n",
    "for ii in np.where(sub_lh)[0]:\n",
    "    sub_mask_lh=mesh_lh.data==ii\n",
    "    sub_mask_lh=sub_mask_lh.any(1)\n",
    "    sub_tri_lh.append(sub_mask_lh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1048,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cort_mask_kill=np.zeros([mesh_lh.data.shape[0]],dtype=bool)\n",
    "cort_mask_keep=np.zeros([mesh_lh.data.shape[0]],dtype=bool)\n",
    "for ii, row in enumerate(mesh_lh.data):\n",
    "    for item in row:\n",
    "        if item in killdict.keys():\n",
    "            cort_mask_kill[ii]=True\n",
    "            continue\n",
    "        else:\n",
    "            cort_mask_keep[ii]=True\n",
    "            continue\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1013,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(59628, 3)"
      ]
     },
     "execution_count": 1013,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cort_tri_v2=mesh_lh.data[cort_mask_kill]\n",
    "cort_tri_v3=mesh_lh.data[cort_mask_keep]\n",
    "cort_tri_v3.shape\n",
    "#cort_mask_lh.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 965,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "md={}\n",
    "for i in v2_ii:\n",
    "    md[i]=1\n",
    "for ii in killdict:\n",
    "    mm=mesh_lh.data==ii\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1024,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "v2_ii=[i for i, e in enumerate(cort_lh) if e != 0] # indices to keep\n",
    "v3_ii=[i for i, e in enumerate(sub_lh) if e != 0] # indices to keep\n",
    "killdict={ii:1 for ii in v3_ii}\n",
    "\n",
    "#=np.zeros([32492])\n",
    "#cort_mask_v2l=(indices>0)\n",
    "#(maskin>0)\n",
    "#cort_mask_lh\n",
    "#mesh_lh.data.shape\n",
    "#cort_mask_lh.shape\n",
    "#maskin=np.where(mesh_lh.data==v2_ii)\n",
    "#maskin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1061,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 1061,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mkeep,mkill=mwmask(data_lh,data_lh_myelin)\n",
    "np.all(mkill == cort_mask_kill)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1051,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1051-982b1d8eb263>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata_lh_myelin\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mvv\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdarrays\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "v=nb.load(nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.L.MyelinMap.32k_fs_LR.func.gii')\n",
    "vv=data_lh_myelin.\n",
    "vv=v.darrays[1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 861,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   68,    12,     0],\n",
       "       [  180,    12,    68],\n",
       "       [  180,    13,    12],\n",
       "       ..., \n",
       "       [ 8440, 32491, 21432],\n",
       "       [ 8440,  8441, 32491],\n",
       "       [    9,  8440, 21432]], dtype=int32)"
      ]
     },
     "execution_count": 861,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#np.asarray(cort_tri_lh).shape\n",
    "#for ii in np.where(cort_lh)[0]:\n",
    "#    print (ii)\n",
    "m=mesh_lh.data==0\n",
    "mesh_lh.data[:,:3]\n",
    "#for ii in indices:\n",
    "#    cort_mask_v2l=mesh_lh.data==ii\n",
    "#    cort_mask_v2l=cort_mask_lh.any(1)\n",
    "#    cort_tri_lh.append(cort_mask_lh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# right\n",
    "cort_tri_rh=[]\n",
    "for ii in np.where(cort_lh)[0]:\n",
    "    cort_mask_rh=mesh_rh.data==ii\n",
    "    cort_mask_rh=cort_mask_rh.any(1)\n",
    "    cort_tri_rh.append(cort_mask_rh)\n",
    "    \n",
    "sub_tri_rh=[]\n",
    "for ii in np.where(sub_rh)[0]:\n",
    "    sub_mask_rh=mesh_rh.data==ii\n",
    "    sub_mask_rh=sub_mask_rh.any(1)\n",
    "    sub_tri_rh.append(sub_mask_rh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1010,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cort_tri_lh=np.array(cort_tri_lh)\n",
    "sub_tri_lh=np.array(sub_tri_lh)\n",
    "cort_tri_rh=np.array(cort_tri_rh)\n",
    "sub_tri_rh=np.array(sub_tri_rh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1011,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cort_mask_lh=cort_tri_lh.any(0)\n",
    "sub_mask_lh=sub_tri_lh.any(0)\n",
    "cort_mask_rh=cort_tri_rh.any(0)\n",
    "sub_mask_rh=sub_tri_rh.any(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 800,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(64980,)"
      ]
     },
     "execution_count": 800,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cort_mask_lh.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 690,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 690,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "brainmodels=[]\n",
    "for idx, bm in enumerate(mim.brain_models):\n",
    "    np.append(brainmodels,np.asarray([idx, bm.index_offset, bm.brain_structure]))\n",
    "#    np.asarray()\n",
    "#    print((idx, bm.index_offset, bm.brain_structure))\n",
    "brainmodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 721,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1200, 29696)\n",
      "(1200, 29717)\n",
      "(1200, 31869)\n",
      "8412.6\n",
      "6808.48\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 10370.02734375,   8975.28808594,   8177.61816406, ...,\n",
       "          7961.00195312,   8160.44628906,   8412.59570312],\n",
       "       [ 10405.33300781,   8885.35449219,   8228.50195312, ...,\n",
       "          8045.67138672,   8196.06152344,   8447.87597656],\n",
       "       [ 10403.03027344,   8885.85253906,   8170.07861328, ...,\n",
       "          8100.45068359,   8327.37402344,   8508.35644531],\n",
       "       ..., \n",
       "       [ 10651.48046875,   8867.24511719,   8167.99414062, ...,\n",
       "          8208.328125  ,   8336.39746094,   8499.38671875],\n",
       "       [ 10547.875     ,   8944.50683594,   8189.71582031, ...,\n",
       "          8144.38916016,   8298.16308594,   8386.21582031],\n",
       "       [ 10372.42089844,   8968.25292969,   8143.32861328, ...,\n",
       "          8001.57714844,   8115.43945312,   8326.44433594]], dtype=float32)"
      ]
     },
     "execution_count": 721,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cclh=cd[:,0:len(lidx)]\n",
    "ccrh=cd[:,len(lidx):len(lidx)+len(ridx)+1]\n",
    "ccsc=cd[:,len(lidx)+len(ridx)+1:cd.shape[1]]\n",
    "print(cclh.shape)\n",
    "print(ccrh.shape)\n",
    "print(ccsc.shape)\n",
    "print(cclh[0,len(lidx)-1])\n",
    "print(ccrh[0,len(lidx)+1])\n",
    "clh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 559,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1200, 91282)\n",
      "(1200, 29696)\n",
      "(1200, 29717)\n",
      "(1200, 31869)\n",
      "91282\n",
      "29696\n",
      "29717\n",
      "31869\n",
      "91282\n",
      "8412.6\n",
      "10799.1\n"
     ]
    }
   ],
   "source": [
    "clh=cd[:,0:29696]\n",
    "crh=cd[:,29696:59413]\n",
    "csc=cd[:,59413:91282]\n",
    "print(cd.shape)\n",
    "print(clh.shape)\n",
    "print(crh.shape)\n",
    "print(csc.shape)\n",
    "print(len(cd[0]))\n",
    "print(len(clh[0]))\n",
    "print(len(crh[0]))\n",
    "print(len(csc[0]))\n",
    "print(len(clh[0])+len(crh[0])+len(csc[0]))\n",
    "print(clh[0,29695])\n",
    "print(crh[0,0])\n",
    "#cd[0,29692:29698]\n",
    "olay_lh=clh.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 575,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(64980, 3)"
      ]
     },
     "execution_count": 575,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mesh_lh.data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1395,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# compute MNI coordinates from voxel coordinates\n",
    "sub_data=subcort.get_data()\n",
    "subcort_disp=sub_data!=0\n",
    "from nibabel.affines import apply_affine\n",
    "\n",
    "real_coords=np.zeros([subcort_disp.shape[0],subcort_disp.shape[1],subcort_disp.shape[2]], dtype = object)\n",
    "xyz_sub=[]#np.zeros([33259,33259,33259],dtype=object)\n",
    "c=0\n",
    "for xx in range (0, subcort_disp.shape[0]):\n",
    "    for yy in range (0, subcort_disp.shape[1]):\n",
    "        for zz in range (0, subcort_disp.shape[2]):\n",
    "            real_coords[xx,yy,zz]=apply_affine(subcort.affine,[xx,yy,zz])\n",
    "#            print(real_coords[xx,yy,zz])\n",
    "#            xyz_sub[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1398,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[array([  90., -126.,  -72.]), array([  90., -126.,  -70.]),\n",
       "         array([  90., -126.,  -68.]), ..., array([  90., -126.,  104.]),\n",
       "         array([  90., -126.,  106.]), array([  90., -126.,  108.])],\n",
       "        [array([  90., -124.,  -72.]), array([  90., -124.,  -70.]),\n",
       "         array([  90., -124.,  -68.]), ..., array([  90., -124.,  104.]),\n",
       "         array([  90., -124.,  106.]), array([  90., -124.,  108.])],\n",
       "        [array([  90., -122.,  -72.]), array([  90., -122.,  -70.]),\n",
       "         array([  90., -122.,  -68.]), ..., array([  90., -122.,  104.]),\n",
       "         array([  90., -122.,  106.]), array([  90., -122.,  108.])],\n",
       "        ..., \n",
       "        [array([ 90.,  86., -72.]), array([ 90.,  86., -70.]),\n",
       "         array([ 90.,  86., -68.]), ..., array([  90.,   86.,  104.]),\n",
       "         array([  90.,   86.,  106.]), array([  90.,   86.,  108.])],\n",
       "        [array([ 90.,  88., -72.]), array([ 90.,  88., -70.]),\n",
       "         array([ 90.,  88., -68.]), ..., array([  90.,   88.,  104.]),\n",
       "         array([  90.,   88.,  106.]), array([  90.,   88.,  108.])],\n",
       "        [array([ 90.,  90., -72.]), array([ 90.,  90., -70.]),\n",
       "         array([ 90.,  90., -68.]), ..., array([  90.,   90.,  104.]),\n",
       "         array([  90.,   90.,  106.]), array([  90.,   90.,  108.])]],\n",
       "\n",
       "       [[array([  88., -126.,  -72.]), array([  88., -126.,  -70.]),\n",
       "         array([  88., -126.,  -68.]), ..., array([  88., -126.,  104.]),\n",
       "         array([  88., -126.,  106.]), array([  88., -126.,  108.])],\n",
       "        [array([  88., -124.,  -72.]), array([  88., -124.,  -70.]),\n",
       "         array([  88., -124.,  -68.]), ..., array([  88., -124.,  104.]),\n",
       "         array([  88., -124.,  106.]), array([  88., -124.,  108.])],\n",
       "        [array([  88., -122.,  -72.]), array([  88., -122.,  -70.]),\n",
       "         array([  88., -122.,  -68.]), ..., array([  88., -122.,  104.]),\n",
       "         array([  88., -122.,  106.]), array([  88., -122.,  108.])],\n",
       "        ..., \n",
       "        [array([ 88.,  86., -72.]), array([ 88.,  86., -70.]),\n",
       "         array([ 88.,  86., -68.]), ..., array([  88.,   86.,  104.]),\n",
       "         array([  88.,   86.,  106.]), array([  88.,   86.,  108.])],\n",
       "        [array([ 88.,  88., -72.]), array([ 88.,  88., -70.]),\n",
       "         array([ 88.,  88., -68.]), ..., array([  88.,   88.,  104.]),\n",
       "         array([  88.,   88.,  106.]), array([  88.,   88.,  108.])],\n",
       "        [array([ 88.,  90., -72.]), array([ 88.,  90., -70.]),\n",
       "         array([ 88.,  90., -68.]), ..., array([  88.,   90.,  104.]),\n",
       "         array([  88.,   90.,  106.]), array([  88.,   90.,  108.])]],\n",
       "\n",
       "       [[array([  86., -126.,  -72.]), array([  86., -126.,  -70.]),\n",
       "         array([  86., -126.,  -68.]), ..., array([  86., -126.,  104.]),\n",
       "         array([  86., -126.,  106.]), array([  86., -126.,  108.])],\n",
       "        [array([  86., -124.,  -72.]), array([  86., -124.,  -70.]),\n",
       "         array([  86., -124.,  -68.]), ..., array([  86., -124.,  104.]),\n",
       "         array([  86., -124.,  106.]), array([  86., -124.,  108.])],\n",
       "        [array([  86., -122.,  -72.]), array([  86., -122.,  -70.]),\n",
       "         array([  86., -122.,  -68.]), ..., array([  86., -122.,  104.]),\n",
       "         array([  86., -122.,  106.]), array([  86., -122.,  108.])],\n",
       "        ..., \n",
       "        [array([ 86.,  86., -72.]), array([ 86.,  86., -70.]),\n",
       "         array([ 86.,  86., -68.]), ..., array([  86.,   86.,  104.]),\n",
       "         array([  86.,   86.,  106.]), array([  86.,   86.,  108.])],\n",
       "        [array([ 86.,  88., -72.]), array([ 86.,  88., -70.]),\n",
       "         array([ 86.,  88., -68.]), ..., array([  86.,   88.,  104.]),\n",
       "         array([  86.,   88.,  106.]), array([  86.,   88.,  108.])],\n",
       "        [array([ 86.,  90., -72.]), array([ 86.,  90., -70.]),\n",
       "         array([ 86.,  90., -68.]), ..., array([  86.,   90.,  104.]),\n",
       "         array([  86.,   90.,  106.]), array([  86.,   90.,  108.])]],\n",
       "\n",
       "       ..., \n",
       "       [[array([ -86., -126.,  -72.]), array([ -86., -126.,  -70.]),\n",
       "         array([ -86., -126.,  -68.]), ..., array([ -86., -126.,  104.]),\n",
       "         array([ -86., -126.,  106.]), array([ -86., -126.,  108.])],\n",
       "        [array([ -86., -124.,  -72.]), array([ -86., -124.,  -70.]),\n",
       "         array([ -86., -124.,  -68.]), ..., array([ -86., -124.,  104.]),\n",
       "         array([ -86., -124.,  106.]), array([ -86., -124.,  108.])],\n",
       "        [array([ -86., -122.,  -72.]), array([ -86., -122.,  -70.]),\n",
       "         array([ -86., -122.,  -68.]), ..., array([ -86., -122.,  104.]),\n",
       "         array([ -86., -122.,  106.]), array([ -86., -122.,  108.])],\n",
       "        ..., \n",
       "        [array([-86.,  86., -72.]), array([-86.,  86., -70.]),\n",
       "         array([-86.,  86., -68.]), ..., array([ -86.,   86.,  104.]),\n",
       "         array([ -86.,   86.,  106.]), array([ -86.,   86.,  108.])],\n",
       "        [array([-86.,  88., -72.]), array([-86.,  88., -70.]),\n",
       "         array([-86.,  88., -68.]), ..., array([ -86.,   88.,  104.]),\n",
       "         array([ -86.,   88.,  106.]), array([ -86.,   88.,  108.])],\n",
       "        [array([-86.,  90., -72.]), array([-86.,  90., -70.]),\n",
       "         array([-86.,  90., -68.]), ..., array([ -86.,   90.,  104.]),\n",
       "         array([ -86.,   90.,  106.]), array([ -86.,   90.,  108.])]],\n",
       "\n",
       "       [[array([ -88., -126.,  -72.]), array([ -88., -126.,  -70.]),\n",
       "         array([ -88., -126.,  -68.]), ..., array([ -88., -126.,  104.]),\n",
       "         array([ -88., -126.,  106.]), array([ -88., -126.,  108.])],\n",
       "        [array([ -88., -124.,  -72.]), array([ -88., -124.,  -70.]),\n",
       "         array([ -88., -124.,  -68.]), ..., array([ -88., -124.,  104.]),\n",
       "         array([ -88., -124.,  106.]), array([ -88., -124.,  108.])],\n",
       "        [array([ -88., -122.,  -72.]), array([ -88., -122.,  -70.]),\n",
       "         array([ -88., -122.,  -68.]), ..., array([ -88., -122.,  104.]),\n",
       "         array([ -88., -122.,  106.]), array([ -88., -122.,  108.])],\n",
       "        ..., \n",
       "        [array([-88.,  86., -72.]), array([-88.,  86., -70.]),\n",
       "         array([-88.,  86., -68.]), ..., array([ -88.,   86.,  104.]),\n",
       "         array([ -88.,   86.,  106.]), array([ -88.,   86.,  108.])],\n",
       "        [array([-88.,  88., -72.]), array([-88.,  88., -70.]),\n",
       "         array([-88.,  88., -68.]), ..., array([ -88.,   88.,  104.]),\n",
       "         array([ -88.,   88.,  106.]), array([ -88.,   88.,  108.])],\n",
       "        [array([-88.,  90., -72.]), array([-88.,  90., -70.]),\n",
       "         array([-88.,  90., -68.]), ..., array([ -88.,   90.,  104.]),\n",
       "         array([ -88.,   90.,  106.]), array([ -88.,   90.,  108.])]],\n",
       "\n",
       "       [[array([ -90., -126.,  -72.]), array([ -90., -126.,  -70.]),\n",
       "         array([ -90., -126.,  -68.]), ..., array([ -90., -126.,  104.]),\n",
       "         array([ -90., -126.,  106.]), array([ -90., -126.,  108.])],\n",
       "        [array([ -90., -124.,  -72.]), array([ -90., -124.,  -70.]),\n",
       "         array([ -90., -124.,  -68.]), ..., array([ -90., -124.,  104.]),\n",
       "         array([ -90., -124.,  106.]), array([ -90., -124.,  108.])],\n",
       "        [array([ -90., -122.,  -72.]), array([ -90., -122.,  -70.]),\n",
       "         array([ -90., -122.,  -68.]), ..., array([ -90., -122.,  104.]),\n",
       "         array([ -90., -122.,  106.]), array([ -90., -122.,  108.])],\n",
       "        ..., \n",
       "        [array([-90.,  86., -72.]), array([-90.,  86., -70.]),\n",
       "         array([-90.,  86., -68.]), ..., array([ -90.,   86.,  104.]),\n",
       "         array([ -90.,   86.,  106.]), array([ -90.,   86.,  108.])],\n",
       "        [array([-90.,  88., -72.]), array([-90.,  88., -70.]),\n",
       "         array([-90.,  88., -68.]), ..., array([ -90.,   88.,  104.]),\n",
       "         array([ -90.,   88.,  106.]), array([ -90.,   88.,  108.])],\n",
       "        [array([-90.,  90., -72.]), array([-90.,  90., -70.]),\n",
       "         array([-90.,  90., -68.]), ..., array([ -90.,   90.,  104.]),\n",
       "         array([ -90.,   90.,  106.]), array([ -90.,   90.,  108.])]]], dtype=object)"
      ]
     },
     "execution_count": 1398,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "real_coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "subcort_disp=subcort.get_data()!=0\n",
    "np.where(subcort_disp)[2]\n",
    "sub_mask=[]\n",
    "sub_mask0=[]\n",
    "for ii in np.where(subcort_disp)[0]:\n",
    "#    sub_mask0=subcort_disp[ii]\n",
    "#    sub_mask0=sub_mask.any(1)\n",
    "    sub_mask.append(real_coords[ii])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#import IPython\n",
    "np.where(subcort_disp)\n",
    "sub_mask1=np.asarray([np.where(subcort_disp)[0],np.where(subcort_disp)[1],np.where(subcort_disp)[2]])\n",
    "print (sub_mask1.shape)\n",
    "print (xyz_lh.data.shape)\n",
    "np.ones_like(sub_mask1)\n",
    "#%whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 634,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(29696,)"
      ]
     },
     "execution_count": 634,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clh[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print (xyz_lh.data.shape)\n",
    "print (mesh_lh.data.shape)\n",
    "\n",
    "print (xyz_lh.data[0:19])\n",
    "print (mesh_lh.data[0:19])\n",
    "\n",
    "print(subcort_disp[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 648,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(29696,)\n",
      "1.0\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "print(cort_olay_lh.shape)\n",
    "min(olay_lh[0])\n",
    "\n",
    "std_olay_lh=(olay-min(olay))/(max(olay)-min(olay))\n",
    "print(max(std_olay_lh))\n",
    "print(min(std_olay_lh))\n",
    "#mlv=maybe_lverts[0].shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 781,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (1200,29696) (1200,) ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-781-e215821ae1ec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      6\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mscaled\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0mscaled_clh\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0molayscale\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mclh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      9\u001b[0m \u001b[0mscaled_clh\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;31m#scaled_olay=olayscale(clh[0])\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-781-e215821ae1ec>\u001b[0m in \u001b[0;36molayscale\u001b[0;34m(olay)\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0molayscale\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;31m# scales input array values to between zero and one, needed for iPyVolume colormaps\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m     \u001b[0mscaled\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mamin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mamax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mamin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mscaled\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (1200,29696) (1200,) "
     ]
    }
   ],
   "source": [
    "### THIS ONE\n",
    "\n",
    "def olayscale(olay):\n",
    "    # scales input array values to between zero and one, needed for iPyVolume colormaps\n",
    "    scaled=(olay-np.amin(olay,olay.ndim-1))/(np.amax(olay,olay.ndim-1)-np.amin(olay,olay.ndim-1))\n",
    "    return scaled\n",
    "\n",
    "scaled_clh=olayscale(clh)\n",
    "scaled_clh\n",
    "#scaled_olay=olayscale(clh[0])\n",
    "#scaled_olay\n",
    "#print(max(clh[0]))\n",
    "#np.amax(clh,1)\n",
    "#clh.ndim\n",
    "#xx=olay-np.amin(olay,olay.ndim-1)\n",
    "#yy=np.amax(olay,olay.ndim-1)-np.amin(olay,olay.ndim-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 777,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1200, 29696)"
      ]
     },
     "execution_count": 777,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clh.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 733,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-733-4642ad7d08ae>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mscaled_olay\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0molayscale\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0many\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscaled_olay\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0mstd_olay_lh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0mscaled_clh\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0molayscale\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mclh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-728-ae243ad9e789>\u001b[0m in \u001b[0;36molayscale\u001b[0;34m(olay)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0molayscale\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0;31m# scales input array values to between zero and one, needed for iPyVolume colormaps\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m     \u001b[0mscaled\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0molay\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mscaled\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"
     ]
    }
   ],
   "source": [
    "std_olay_lh=np.asarray((olay-min(olay))/(max(olay)-min(olay)))\n",
    "std_olay_lh.shape\n",
    "olay.shape\n",
    "np.corrcoef(std_olay_lh,olay)\n",
    "clh.shape\n",
    "scaled_olay=olayscale(olay)\n",
    "np.any(scaled_olay==std_olay_lh)\n",
    "scaled_clh=olayscale(clh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 782,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "olay=np.zeros([32492])\n",
    "values=clh[0]\n",
    "count=0\n",
    "olay[lidx] = values\n",
    "#for x in 0:len(lidx):\n",
    "#    olay[lidx[x]]=values[x]\n",
    "#    count=count+1\n",
    "    \n",
    "    \n",
    "#for ii in np.where(cort_lh)[0]:\n",
    "#    cort_mask_lh=mesh_lh.data==ii\n",
    "#    cort_mask_lh=cort_mask_lh.any(1)\n",
    "#    cort_tri_lh.append(cort_mask_lh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1065,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "can't assign to keyword (<ipython-input-1065-0ff3e75db027>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-1065-0ff3e75db027>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    [None,a]=mwmask(data_lh,data_lh_myelin)\u001b[0m\n\u001b[0m                                           ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m can't assign to keyword\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1143,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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": 1149,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c4b56f502a384712b3a3fead1dd221f9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 1149,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "showSurface(data_lh,data_lh_myelin,newfigure=True,mwall=True) # comment 0\n",
    "showSurface(data_rh,data_rh_myelin,newfigure=False,mwall=False) # comment 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1063,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d491280e17a54deda097142d696345a4",
       "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",
    "#color_lh=plt.cm.autumn(olay_lh)\n",
    "color_lh=plt.cm.viridis(std_olay_lh)\n",
    "#color_mlh=plt.cm.autumn(np.log(olay_lh+sys.float_info.epsilon))\n",
    "color_rh=plt.cm.viridis(np.log(olay_rh+sys.float_info.epsilon))\n",
    "\n",
    "#colors\n",
    "x, y, z = xyz_lh.data.T\n",
    "#surf_mlh = p3.plot_trisurf(y, z, x, triangles=mesh_lh.data[cort_mask_lh], color=color_mlh[:,:3])\n",
    "\n",
    "surf_lh = p3.plot_trisurf(y,z,x,triangles=mesh_lh.data[mkeep],color=color_lh[:,:3])\n",
    "sub_lh = p3.plot_trisurf(y,z,x, triangles=mesh_lh.data[mkill], color=color_lh[:,:3])\n",
    "x, y, z = xyz_rh.data.T\n",
    "#surf_rh = p3.plot_trisurf(y*1.1, z*1.1, x*1.1, triangles=mesh_rh.data[cort_mask_rh], color=color_rh[:,:3])\n",
    "surf_rh = p3.plot_trisurf(y,z,x,triangles=mesh_rh.data[cort_mask_rh],color=color_rh[:,:3])\n",
    "#sub_rh = p3.plot_trisurf(x, y, z, triangles=mesh_rh.data[sub_mask_rh], color=color_rh[:,:3])\n",
    "\n",
    "x, y, z = sub_mask1\n",
    "#surf_sub = p3.plot_trisurf(x, y, z,triangles=np.ones_like(sub_mask1),color='grey')\n",
    "\n",
    "#x, y, z = xyz_rh.data.T\n",
    "#surf_rh = p3.plot_trisurf(x, y, z, triangles=mesh_rh.data, color=color_rh[:,: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": 652,
   "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": "bbedda88dbac41e997abde8de2d6ac12",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n",
    "from skimage import measure\n",
    "#from skimage.draw import ellipsoid\n",
    "#sub_mask2=np.asarray(sub_mask)\n",
    "#sub_mask2.ndim\n",
    "#verts, faces, normals, values = measure.marching_cubes_lewiner(sub_mask2.T, 0)\n",
    "#%whos\n",
    "#sorted(xyz_lh.data.T[1])\n",
    "#print(max(xyz_lh.data.T[1]))\n",
    "#print(min(xyz_lh.data.T[1]))\n",
    "#xyz_lh.data.T.shape\n",
    "#real_coords.data.shape\n",
    "#p3.plot_wireframe(x,y,z, color='red', wrapx=False, wrapy=False)\n",
    "#x.shape\n",
    "p3.volshow(subcort_disp)\n",
    "#p3.volshow(real_coords)\n",
    "#p3.plot_isosurface(subcort_disp)\n",
    "#p3.plot_isosurface(real_coords)\n",
    "#p3.plot\n",
    "#x, y, z = real_coords.T\n",
    "#p3.plot_trisurf(x, y, z)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "sub_mask2.T.ndim\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "p3.plot_trisurf(x, y, z,triangles=np.ones_like(sub_mask1),color='grey')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 787,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "347b1e9ad6ad44af9d53453526aa8d75",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "echo(fig)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "color2_lh = plt.cm.autumn(np.log(olay_lh+sys.float_info.epsilon))\n",
    "color2_rh = plt.cm.autumn(np.log(olay_rh+sys.float_info.epsilon))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "surf_lh.color = color2_lh[:, :3]\n",
    "surf_rh.color = color2_rh[:, :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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "solution2": "hidden"
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.hist(data.flatten(), 128);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.hist(mean_epi, 128);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "foo = PCA(n_components=2).fit_transform(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "solution2": "hidden"
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false,
    "solution2": "hidden"
   },
   "outputs": [],
   "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()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Look up how to use PCA from scikit learn\n",
    "\n",
    "# Follow the steps from above to plot the two components\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<b> ~~~~~ BONEYARD ~~~~~ </b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from urllib.request import urlretrieve\n",
    "urlretrieve('https://osf.io/7jbpr/download?version=1', '/tmp/T1.nii.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Here is an example of mesh plotting 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": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# code for general 3D rotation\n",
    "from numpy import cross, eye, dot\n",
    "from scipy.linalg import expm3, norm\n",
    "import math\n",
    "def M(axis, theta):\n",
    "    return expm3(cross(eye(3), axis/norm(axis)*theta))\n",
    "\n",
    "v, axis, theta = [3,5,0], [1,0,0], math.pi/2\n",
    "M0 = M(axis, theta)\n",
    "v*M0"
   ]
  }
 ],
 "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
}