20170919_Xu_Lab_Presentation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Surface Visualization Demo Notebook, Neurohackweek, 8 Sept 2017\n",
    "\n",
    "- (Mostly) developted by Murat Bilgel, Benjamin Ely, Melanie Ganz, Krisanne Litinas, and Andrea Shafer\n",
    "- Heavily indebted to Satra Ghosh, Chris Holdgraf, Anisha Keshavan, and Tal Yarkoni\n",
    "- Builds on \n",
    "[IPyvolume](https://ipyvolume.readthedocs.io/en/latest/index.html), \n",
    "[NiBabel cifti2 package](http://nipy.org/nibabel/reference/nibabel.cifti2.html), \n",
    "[Michiel Cottaar's cifti tools](https://git.fmrib.ox.ac.uk/ndcn0236/cifti), and [various NHW demo notebooks](https://github.com/neurohackweek)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Load packages etc:</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import nibabel as nb\n",
    "import ipyvolume\n",
    "import ipyvolume.pylab as p3\n",
    "import numpy as np\n",
    "import os, sys, time\n",
    "import niwidgets\n",
    "from niwidgets import NiftiWidget\n",
    "from niwidgets import niwidget_surface\n",
    "import nilearn.plotting as nip\n",
    "from ipywidgets import interact, fixed, IntSlider, FloatSlider\n",
    "from ipywidgets import interact,fixed,IntSlider\n",
    "from nilearn import image\n",
    "from itertools import islice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": [
    "# NB: added 'export NHW=$HOME/Desktop/nhw2017/NHW_visualization' to ~/.bashrc\n",
    "nhw = os.environ['NHW'] \n",
    "# FYI example data primarily from Human Connectome Project subject 100610\n",
    "path_lh = nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.L.very_inflated_MSMAll.32k_fs_LR.surf.gii'\n",
    "data_lh = nb.load(path_lh)\n",
    "path_rh = nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.R.midthickness_MSMAll.32k_fs_LR.surf.gii'\n",
    "data_rh = nb.load(path_rh)\n",
    "path_sc=nhw+'/data/100610/MNINonLinear/ROIs/Atlas_ROIs.2.nii.gz'\n",
    "data_sc = nb.load(path_sc)\n",
    "\n",
    "path_lh_myelin = nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.L.MyelinMap.32k_fs_LR.func.gii'\n",
    "data_lh_myelin = nb.load(path_lh_myelin)\n",
    "path_rh_myelin = nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.R.MyelinMap.32k_fs_LR.func.gii'\n",
    "data_rh_myelin = nb.load(path_rh_myelin)\n",
    "\n",
    "path_lh_corthick = nhw+'/data/100610/MNINonLinear/fsaverage_LR32k/100610.L.corrThickness.32k_fs_LR.shape.gii'\n",
    "data_lh_corthick = nb.load(path_lh_corthick)\n",
    "\n",
    "filename = os.path.join(nhw,'data','subjA','subjA_t0_pet_mni.nii.gz')\n",
    "filename_4D = os.path.join(nhw,'data','subjA','subjA_4D_pet_mni.nii.gz')\n",
    "\n",
    "path_ciftiTS = nhw+'/data/100610/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii'\n",
    "data_ciftiTS = nb.load(path_ciftiTS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Surface display functions:</b>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def showSurface(surface,overlay=None,frame=0,newfigure=True,colormap='summer',figsize=np.array([600,600]),\n",
    "                figlims=np.array([[-100,100],[-100,100],[-100,100]]),showZeroes=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 object\n",
    "        Path to surface file in gifti or FS surface format or an already loaded gifti object of surface\n",
    "    overlay: str, gifti object\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",
    "        index 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",
    "    showZeroes: bool\n",
    "        Display vertices with intensity = 0, default True\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",
    "    \n",
    "    if not isinstance(frame,int):\n",
    "        raise 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",
    "#        elif (file_extension in ('.dlabel.nii','.dscalar.nii','.dconn.nii','.dtseries.nii')):\n",
    "#            left,right,sub=csplit(overlay)\n",
    "\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",
    "    \n",
    "    if showZeroes is False:\n",
    "        try:\n",
    "            mkeep,mkill=zmask(surface,overlay)\n",
    "        except:\n",
    "            raise ValueError('Overlay required for medial wall masking.')\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 showZeroes is False:\n",
    "        my_color = plt.cm.get_cmap(colormap)\n",
    "        colors=my_color(oscale(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(oscale(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": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def zmask(surf,mask):\n",
    "    '''\n",
    "    Masks out vertices with intensity=0 from overlay. Also returns masked-out vertices.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    surf: gifti object\n",
    "        already loaded gifti object of target surface\n",
    "    mask: gifti object\n",
    "        already loaded gifti object of overlay with zeroes at vertices of no interest (e.g. medial wall)\n",
    "    '''\n",
    "    keep=(mask.darrays[0].data!=0) # nonzero values of mask\n",
    "    kill=(mask.darrays[0].data==0) # zero values of mask\n",
    "    ikeep=[i for i, e in enumerate(keep) if e != 0] # indices of nonzero mask values\n",
    "    ikill=[i for i, e in enumerate(kill) if e != 0] # indices of zero mask values\n",
    "    killdict={ii:1 for ii in ikill} # fun fact, iterating over a dictionary is ~exponentially faster vs. over a list\n",
    "    mask_kill=np.zeros([surf.darrays[1].data.shape[0]],dtype=bool) # create empty arrays matching surface mesh dimentions\n",
    "    mask_keep=mask_kill.copy()\n",
    "    for ii, row in enumerate(surf.darrays[1].data):\n",
    "        for item in row:\n",
    "            if item in killdict.keys():\n",
    "                mask_kill[ii]=True\n",
    "                continue\n",
    "            else:\n",
    "                mask_keep[ii]=True\n",
    "                continue\n",
    "    return mask_keep, mask_kill"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mwpad(overlay,template,frame=0):\n",
    "    '''\n",
    "    Adds intensity=0 vertices along medial wall of overlay\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    overlay: gifti object\n",
    "        gifti object overlay to be padded (e.g. hemisphere data array from a cifti file)\n",
    "    template: gifti object\n",
    "        gifti object with zeroes at medial wall vertices, otherwise matched to overlay (e.g. myelin map)\n",
    "    frame: int\n",
    "        index of the timepoint to show, default 0, ignored if overlay only has one timepoint\n",
    "    '''\n",
    "    padolay=nb.gifti.gifti.GiftiImage()\n",
    "    padded=np.zeros(template.darrays[0].data.shape[0])\n",
    "    count=0        \n",
    "    for i, x in enumerate(template.darrays[0].data):\n",
    "        if x!=0:\n",
    "            if overlay.darrays[0].data.ndim > 1:\n",
    "                padded[i]=overlay.darrays[0].data[frame,count]\n",
    "            else:\n",
    "                padded[i]=overlay.darrays[0].data[count]\n",
    "            count=count+1\n",
    "    padarray=nb.gifti.gifti.GiftiDataArray(padded)\n",
    "    padolay.add_gifti_data_array(padarray)\n",
    "    return padolay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def oscale(overlay):\n",
    "    '''\n",
    "    Scales input array values to between zero and one, as needed for iPyVolume colormaps.\n",
    "    For now, olay must be a 1D array (single timepoint).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    overlay: gifti object\n",
    "        gifti object overlay to be scaled (e.g. myelin map)\n",
    "    '''\n",
    "    scaled=np.asarray((overlay-min(overlay))/(max(overlay)-min(overlay)))\n",
    "    return scaled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7b1deb76091148f2968c583db8164996",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# showSurface example 1:\n",
    "# right surface only\n",
    "showSurface(data_rh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "49c5e0d3407d4236b5cd6363d878a799",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# showSurface example 2:\n",
    "# new plot, same surface but with myelin map overlay\n",
    "fig=showSurface(data_rh,data_rh_myelin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# same plot, added \"very inflated\" left surface, excluded medial wall\n",
    "fig=showSurface(data_lh,data_lh_corthick,newfigure=False,showZeroes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>NiWidgets 4D longitudinal or timeseries display:</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 1) get the default plotter in niwidgets and manipulate it.\n",
    "# 2) write function for an spaghetti plot showing the longitudinal trend at selected voxel\n",
    "\n",
    "def plot_timepoint_slices(data, figure, x, y, z, t=0, colormap='viridis', figsize=(15, 5)):\n",
    "    \"\"\"\n",
    "    Plots x,y,z slices at time point t.\n",
    "    \"\"\"\n",
    "    \n",
    "    coords = [x, y, z]\n",
    "    views = ['Sagittal', 'Coronal', 'Axial']\n",
    "    \n",
    "    axes = [figure.add_subplot(231), figure.add_subplot(232), figure.add_subplot(233),\n",
    "            figure.add_subplot(234)]\n",
    "    \n",
    "    \n",
    "    data = data.dataobj.get_unscaled()\n",
    "    print(data.shape)\n",
    "    for subplot in range(3):\n",
    "        slice_obj = 3 * [slice(None)] + [t]\n",
    "        slice_obj[subplot] = coords[subplot]\n",
    "\n",
    "        plt.sca(axes[subplot])\n",
    "        axes[subplot].set_facecolor('black')\n",
    "        axes[subplot].set_title(views[subplot])\n",
    "        axes[subplot].tick_params(\n",
    "            axis='both', which='both', bottom='off', top='off',\n",
    "            labelbottom='off', right='off', left='off', labelleft='off'\n",
    "            )\n",
    "        # fix the axis limits\n",
    "        axis_limits = [limit for i, limit in enumerate(data.shape)\n",
    "                       if i != subplot]\n",
    "        axes[subplot].set_xlim(0, axis_limits[0])\n",
    "        axes[subplot].set_ylim(0, axis_limits[1])\n",
    "\n",
    "        # plot the actual slice\n",
    "        if subplot == 0:\n",
    "            plt.imshow(np.flipud(np.rot90(data[slice_obj], k=1)),\n",
    "                       cmap=colormap)\n",
    "        else:\n",
    "            plt.imshow(np.rot90(data[slice_obj], k=3), cmap=colormap)\n",
    "        # draw guides to show where the other two slices are\n",
    "        guide_positions = [val for i, val in enumerate(coords)\n",
    "                           if i != subplot]\n",
    "        plt.axvline(x=guide_positions[0], color='gray', alpha=0.8)\n",
    "        plt.axhline(y=guide_positions[1], color='gray', alpha=0.8)\n",
    "\n",
    "    #plt.sca(axes[-1])\n",
    "    axes[-1].plot(range(data.shape[3]),data[x,y,z,:],'o-')\n",
    "    axes[-1].plot(t,data[x,y,z,t],'ro')\n",
    "    axes[-1].set_xlabel('Time point')\n",
    "    axes[-1].set_ylabel('Intensity')\n",
    "    \n",
    "    # show the plot\n",
    "    plt.show()\n",
    "    # print the value at that point in case people need to know\n",
    "    print('Value at point {x}, {y}, {z}, {t}: {intensity}'.format(\n",
    "        x=x, y=y, z=z, t=t, intensity=data[x, y, z, t]\n",
    "    ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6ee1e6a4747442eca4ad6e39e5ef734b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# example volumetric 4D data:\n",
    "my_4D_widget = NiftiWidget(filename_4D)\n",
    "my_4D_widget.nifti_plotter(plotting_func=plot_timepoint_slices, \n",
    "                           x=IntSlider(90, min=0, max=181, continuous_update=False),\n",
    "                           y=IntSlider(100, min=0, max=217),\n",
    "                           z=IntSlider(90, min=0, max=181),\n",
    "                           t=IntSlider(min=0, max=3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d327f6bc0c9744108f520eba1888b1c2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d1381aaadb5249619cb3376fda2ffeac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# example surface 4D data:\n",
    "meshfile_lh_inflated = nhw+'/data/006/surf/lh.inflated'\n",
    "overlayfiles_lh = [nhw+'/data/006/surf/lh.thickness',nhw+'/data/006/surf/lh.curv']\n",
    "sw = niwidget_surface.SurfaceWidget(meshfile=meshfile_lh_inflated,\n",
    "                                    overlayfiles = overlayfiles_lh)\n",
    "sw.surface_plotter()\n",
    "sw.fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>CIFTI tools:</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cread(infile):\n",
    "    '''\n",
    "    Loads a cifti file and stores as five objects: data array, cifti header, nifti header, extra metadata, filemap\n",
    "    Example call:\n",
    "        cd,ch,cn,cx,cf=cread(path_to_cifti)\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    infile: str, Cifti2Image\n",
    "        string containing path to cifti file or an already loaded cifti object\n",
    "    '''\n",
    "    if isinstance(infile,nb.cifti2.cifti2.Cifti2Image):\n",
    "        try:\n",
    "            cin=infile\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid cifti file.')\n",
    "    elif isinstance(infile,str):\n",
    "        if not os.path.exists(infile):\n",
    "            raise ValueError('File does not exist, please provide a valid path to a cifti file.')\n",
    "        try:\n",
    "            cin=nb.load(infile) # NB causes apparently benign warning as of NiBabel version 2.2.0dev\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid cifti 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": 21,
   "metadata": {},
   "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": [
      "(1200, 91282)\n"
     ]
    }
   ],
   "source": [
    "# cread example 1:\n",
    "cd,ch,cn,cx,cf=cread(path_ciftiTS)\n",
    "print(cd.data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1200, 91282)\n"
     ]
    }
   ],
   "source": [
    "# cread example 2:\n",
    "cd,ch,cn,cx,cf=cread(data_ciftiTS)\n",
    "print(cd.data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variable                Type                Data/Info\n",
      "-----------------------------------------------------\n",
      "FloatSlider             MetaHasTraits       <class 'ipywidgets.widget<...>idget_float.FloatSlider'>\n",
      "IntSlider               MetaHasTraits       <class 'ipywidgets.widgets.widget_int.IntSlider'>\n",
      "NiftiWidget             type                <class 'niwidgets.niwidget.NiftiWidget'>\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 0x11600b278>\n",
      "cn                      Nifti2Header        <class 'nibabel.nifti2.Ni<...> 0\\nunused_str      : b''\n",
      "cread                   function            <function cread at 0x114051158>\n",
      "cx                      dict                n=0\n",
      "data_ciftiTS            Cifti2Image         <nibabel.cifti2.cifti2.Ci<...>ge object at 0x112cd5898>\n",
      "data_lh                 GiftiImage          <nibabel.gifti.gifti.Gift<...>ge object at 0x112cd56d8>\n",
      "data_lh_corthick        GiftiImage          <nibabel.gifti.gifti.Gift<...>ge object at 0x112cd1f28>\n",
      "data_lh_myelin          GiftiImage          <nibabel.gifti.gifti.Gift<...>ge object at 0x112cd13c8>\n",
      "data_rh                 GiftiImage          <nibabel.gifti.gifti.Gift<...>ge object at 0x112cd5ba8>\n",
      "data_rh_myelin          GiftiImage          <nibabel.gifti.gifti.Gift<...>ge object at 0x112cd1e10>\n",
      "data_sc                 Nifti1Image         <class 'nibabel.nifti1.Ni<...>nmagic           : b'n+1'\n",
      "fig                     int                 1\n",
      "filename                str                 /Users/ely/Desktop/nhw201<...>A/subjA_t0_pet_mni.nii.gz\n",
      "filename_4D             str                 /Users/ely/Desktop/nhw201<...>A/subjA_4D_pet_mni.nii.gz\n",
      "fixed                   MetaHasTraits       <class 'ipywidgets.widgets.interaction.fixed'>\n",
      "image                   module              <module 'nilearn.image' f<...>learn/image/__init__.py'>\n",
      "interact                _InteractFactory    <ipywidgets.widgets.inter<...>ry object at 0x106800828>\n",
      "ipyvolume               module              <module 'ipyvolume' from <...>s/ipyvolume/__init__.py'>\n",
      "islice                  type                <class 'itertools.islice'>\n",
      "meshfile_lh_inflated    str                 /Users/ely/Desktop/nhw201<...>data/006/surf/lh.inflated\n",
      "mwpad                   function            <function mwpad at 0x11364e6a8>\n",
      "my_4D_widget            NiftiWidget         <niwidgets.niwidget.Nifti<...>et object at 0x11683cf28>\n",
      "nb                      module              <module 'nibabel' from '/<...>bel/nibabel/__init__.py'>\n",
      "nhw                     str                 /Users/ely/Desktop/nhw2017/NHW_visualization\n",
      "nip                     module              <module 'nilearn.plotting<...>rn/plotting/__init__.py'>\n",
      "niwidget_surface        module              <module 'niwidgets.niwidg<...>ets/niwidget_surface.py'>\n",
      "niwidgets               module              <module 'niwidgets' from <...>s/niwidgets/__init__.py'>\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",
      "oscale                  function            <function oscale at 0x113679e18>\n",
      "overlayfiles_lh         list                n=2\n",
      "p3                      module              <module 'ipyvolume.pylab'<...>ages/ipyvolume/pylab.py'>\n",
      "path_ciftiTS            str                 /Users/ely/Desktop/nhw201<...>hp2000_clean.dtseries.nii\n",
      "path_lh                 str                 /Users/ely/Desktop/nhw201<...>MSMAll.32k_fs_LR.surf.gii\n",
      "path_lh_corthick        str                 /Users/ely/Desktop/nhw201<...>kness.32k_fs_LR.shape.gii\n",
      "path_lh_myelin          str                 /Users/ely/Desktop/nhw201<...>linMap.32k_fs_LR.func.gii\n",
      "path_rh                 str                 /Users/ely/Desktop/nhw201<...>MSMAll.32k_fs_LR.surf.gii\n",
      "path_rh_myelin          str                 /Users/ely/Desktop/nhw201<...>linMap.32k_fs_LR.func.gii\n",
      "path_sc                 str                 /Users/ely/Desktop/nhw201<...>/ROIs/Atlas_ROIs.2.nii.gz\n",
      "plot_timepoint_slices   function            <function plot_timepoint_slices at 0x112cd0f28>\n",
      "plt                     module              <module 'matplotlib.pyplo<...>es/matplotlib/pyplot.py'>\n",
      "showSurface             function            <function showSurface at 0x113649158>\n",
      "sw                      SurfaceWidget       <niwidgets.niwidget_surfa<...>et object at 0x114052c88>\n",
      "sys                     module              <module 'sys' (built-in)>\n",
      "time                    module              <module 'time' (built-in)>\n",
      "zmask                   function            <function zmask at 0x113649c80>\n"
     ]
    }
   ],
   "source": [
    "%whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def csave(filename,cdata,chead,cnhdr=None,cxtra=None,cfmap=None,returnImg=False):\n",
    "    '''\n",
    "    Saves a cifti image from a collection of cifti substructures to file.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    filename: str\n",
    "        string containing output file path\n",
    "    cdata: numpy ndarray\n",
    "        array containing data to be saved, can be 1D (single timepoint) or 2D (timeseries)\n",
    "    chead: Cifti2Header\n",
    "        cifti header data, generally extracted from a previously read cifti file\n",
    "    cnhdr: Nifti2Header (optional)\n",
    "        nifti header data, also generally extracted from a previously read cifti file\n",
    "    cxtra: dict (optional)\n",
    "        extra cifti metadata, generally empty\n",
    "    cfmap: dict (optional)\n",
    "        file map, generally contains the string 'image'\n",
    "    returnImg: bool\n",
    "        optionally return the compiled cifti object\n",
    "    '''\n",
    "    cout=nb.Cifti2Image(dataobj=cdata,header=chead,nifti_header=cnhdr, extra=cxtra, file_map=cfmap)\n",
    "    cout.to_filename(filename)\n",
    "    if returnImg is True:\n",
    "        return cout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# csave example:\n",
    "output_cifti=nhw+'/save_demo.dtseries.nii'\n",
    "csave(output_cifti,cd,ch,cn,ch,cf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def clthresh(data,thresh,setTo=0):\n",
    "    '''\n",
    "    Set all values below threshold for a cifti file to a given value (default 0).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data: numpy ndarray\n",
    "        array containing data to be thresholded\n",
    "    thresh: float\n",
    "        Threshold value (exclusive)\n",
    "    setTo: float\n",
    "        Value to which voxels/vertices are set if below the threshold, default 0\n",
    "    '''\n",
    "    # \n",
    "    lthr_data=data.copy()\n",
    "    lthr_data[lthr_data<thresh]=setTo\n",
    "    return lthr_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cuthresh(data,thresh,setto=0):\n",
    "    '''\n",
    "    Set all values above threshold for a cifti file to a given value (default 0).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data: numpy ndarray\n",
    "        array containing data to be thresholded\n",
    "    thresh: float\n",
    "        Threshold value (exclusive)\n",
    "    setTo: float\n",
    "        Value to which voxels/vertices are set if above the threshold, default 0\n",
    "    '''\n",
    "    # set all values above threshold to a given value (default 0)\n",
    "    uthr_data=data.copy()\n",
    "    uthr_data[uthr_data>thresh]=setto\n",
    "    return uthr_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAD8CAYAAABQFVIjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFqBJREFUeJzt3X+s39V93/HnqzhQlIYfBs/yDMyO4k4ilZaWKwepadXB\najtpF1jLmKdteBoKmmBTom0qsEijS1IJb1rZ0BYiWhCGpTWMNsLqwjwH2KpJ48d1ICWQMN8EEFgG\nE9tAq22spu/98T0Xff2t7Xvte3y/9/o+H9JX9/M9n3PO5/hzr7+vez6/bqoKSZLm6sfGPQBJ0unB\nQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSepi2bgHMJ8uvPDCWrNmzbiHIUmL\nyu7du39UVStmqrekAmXNmjVMTk6OexiStKgkeXU29TzkJUnqwkCRJHVhoEiSujBQJEldGCiSpC4M\nFElSFwaKJKkLA0WS1IWBIknqYkndKS/1tOaW//znyl65/ZfGMBJpYXCGIknqwkCRJHVhoEiSuphV\noCR5JcnzSZ5LMtnKlifZlWRP+3r+UP1bk0wleSnJxqHyy1o/U0nuTJJWflaSB1v5U0nWDLXZ0rax\nJ8mWofK1re5Ua3vm3HeHJOlknchJ+b9aVT8aen8L8FhV3Z7klvb+5iSXApuBjwN/EfhWkp+sqveB\nu4DPAU8B3wQ2AY8C1wOHqupjSTYDW4G/lWQ5cBswARSwO8mOqjrU6txRVduTfK31cddJ7gfpuI52\nAl7SkeZyyOsqYFtb3gZcPVS+vareq6qXgSlgfZJVwDlV9WRVFXD/SJvpvh4Grmyzl43Arqo62EJk\nF7Cprbui1R3dviRpDGYbKMVgprE7yQ2tbGVV7WvLbwAr2/Jq4LWhtq+3stVtebT8iDZVdRh4B7jg\nOH1dALzd6o72JUkag9ke8vpUVe1N8heAXUm+P7yyqipJ9R/e3LUAvAHgkksuGfNoJOn0NatAqaq9\n7ev+JN8A1gNvJllVVfva4az9rfpe4OKh5he1sr1tebR8uM3rSZYB5wIHWvkvjLT5b23deUmWtVnK\ncF+jY78buBtgYmJiQYaeTh/e7KilbMZDXkk+nOQj08vABuC7wA5g+qqrLcAjbXkHsLldubUWWAc8\n3Q6PvZvk8nYO5LqRNtN9XQM83s6z7AQ2JDm/XUW2AdjZ1j3R6o5uX5I0BrOZoawEvtGu8F0G/E5V\n/ZckzwAPJbkeeBW4FqCqXkjyEPAicBi4qV3hBXAjcB9wNoOrux5t5fcADySZAg4yuEqMqjqY5MvA\nM63el6rqYFu+Gdie5CvAs60PSdKYZPDL/tIwMTFRk5OT4x6GFqG5XDbsIS8tdkl2V9XETPV8OKR0\nFN53Ip04H70iSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhTc2SqfY6E2Sp/rO\neR9QqXFxhiJJ6sIZirSI+YgYLSQGirSIGCBayAwUaZ55jkOnKwNFWgIMMc0HT8pLkrowUCRJXRgo\nkqQuPIciLVGeV1FvBoqWPC/FlfowUKQFyqDTYuM5FElSF85QpAVgocxG5vtBljq9OEORJHVhoEiS\nujBQJEldGCiSpC4MFElSF17lJemYvJteJ8IZiiSpCwNFktTFrAMlyRlJnk3yB+398iS7kuxpX88f\nqntrkqkkLyXZOFR+WZLn27o7k6SVn5XkwVb+VJI1Q222tG3sSbJlqHxtqzvV2p45t10hSZqLE5mh\nfB743tD7W4DHqmod8Fh7T5JLgc3Ax4FNwFeTnNHa3AV8DljXXpta+fXAoar6GHAHsLX1tRy4Dfgk\nsB64bSi4tgJ3tDaHWh+SpDGZVaAkuQj4JeC3h4qvAra15W3A1UPl26vqvap6GZgC1idZBZxTVU9W\nVQH3j7SZ7uth4Mo2e9kI7Kqqg1V1CNgFbGrrrmh1R7cvSRqD2c5Q/i3wa8CfDZWtrKp9bfkNYGVb\nXg28NlTv9Va2ui2Plh/RpqoOA+8AFxynrwuAt1vd0b4kSWMwY6Ak+WVgf1XtPladNuOongPrJckN\nSSaTTL711lvjHo4knbZmM0P5WeCzSV4BtgNXJPmPwJvtMBbt6/5Wfy9w8VD7i1rZ3rY8Wn5EmyTL\ngHOBA8fp6wBwXqs72tcRquruqpqoqokVK1bM4p8rSToZMwZKVd1aVRdV1RoGJ9sfr6q/C+wApq+6\n2gI80pZ3AJvblVtrGZx8f7odHns3yeXtHMh1I22m+7qmbaOAncCGJOe3k/EbgJ1t3ROt7uj2JUlj\nMJc75W8HHkpyPfAqcC1AVb2Q5CHgReAwcFNVvd/a3AjcB5wNPNpeAPcADySZAg4yCC6q6mCSLwPP\ntHpfqqqDbflmYHuSrwDPtj4kSWOSwS/7S8PExERNTk6OexhaYBbKH7daLHz0ytKTZHdVTcxUzzvl\nJUldGCiSpC4MFElSFz6+XtIJ8ZH2OhZnKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdWGgSJK6\nMFAkSV0YKJKkLgwUSVIXBookqQuf5SVpzny+l8AZiiSpEwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJ\nXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sJneWlJOdozpyT14QxFktSFMxRJp8TobNCnD5/+nKFI\nkrqYMVCS/HiSp5N8J8kLSf5lK1+eZFeSPe3r+UNtbk0yleSlJBuHyi9L8nxbd2eStPKzkjzYyp9K\nsmaozZa2jT1JtgyVr211p1rbM/vsEknSyZjNDOU94Iqq+ivAJ4BNSS4HbgEeq6p1wGPtPUkuBTYD\nHwc2AV9Nckbr6y7gc8C69trUyq8HDlXVx4A7gK2tr+XAbcAngfXAbUPBtRW4o7U51PqQJI3JjIFS\nA3/S3n6ovQq4CtjWyrcBV7flq4DtVfVeVb0MTAHrk6wCzqmqJ6uqgPtH2kz39TBwZZu9bAR2VdXB\nqjoE7GIQaAGuaHVHty9JGoNZnUNJckaS54D9DD7gnwJWVtW+VuUNYGVbXg28NtT89Va2ui2Plh/R\npqoOA+8AFxynrwuAt1vd0b5Gx35Dkskkk2+99dZs/rmSpJMwq0Cpqver6hPARQxmGz81sr4YzFoW\nnKq6u6omqmpixYoV4x6OJJ22Tugqr6p6G3iCwbmPN9thLNrX/a3aXuDioWYXtbK9bXm0/Ig2SZYB\n5wIHjtPXAeC8Vne0L0nSGMzmKq8VSc5ry2cDvwh8H9gBTF91tQV4pC3vADa3K7fWMjj5/nQ7PPZu\nksvbOZDrRtpM93UN8Hib9ewENiQ5v52M3wDsbOueaHVHty9JGoPZ3Ni4CtjWrtT6MeChqvqDJP8T\neCjJ9cCrwLUAVfVCkoeAF4HDwE1V9X7r60bgPuBs4NH2ArgHeCDJFHCQwVViVNXBJF8Gnmn1vlRV\nB9vyzcD2JF8Bnm19SJLGJINf9peGiYmJmpycHPcwNEY+y2t8vFN+8Uqyu6omZqrno1ckzYujhbkh\nc3rx0SuSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkL\nn+UlaWx8vtfpxRmKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS\n1IWBIknqwkCRJHVhoEiSuvBpw5IWFJ9AvHg5Q5EkdWGgSJK6MFAkSV0YKJKkLmYMlCQXJ3kiyYtJ\nXkjy+Va+PMmuJHva1/OH2tyaZCrJS0k2DpVfluT5tu7OJGnlZyV5sJU/lWTNUJstbRt7kmwZKl/b\n6k61tmf22SWSpJMxmxnKYeCfVtWlwOXATUkuBW4BHquqdcBj7T1t3Wbg48Am4KtJzmh93QV8DljX\nXpta+fXAoar6GHAHsLX1tRy4DfgksB64bSi4tgJ3tDaHWh+SpDGZMVCqal9Vfbst/zHwPWA1cBWw\nrVXbBlzdlq8CtlfVe1X1MjAFrE+yCjinqp6sqgLuH2kz3dfDwJVt9rIR2FVVB6vqELAL2NTWXdHq\njm5fkjQGJ3QfSjsU9dPAU8DKqtrXVr0BrGzLq4Enh5q93sr+tC2Plk+3eQ2gqg4neQe4YLh8pM0F\nwNtVdfgofY2O+QbgBoBLLrlk1v9WLX5Hu59B0qkz65PySX4C+D3gC1X17vC6NuOozmProqrurqqJ\nqppYsWLFuIcjSaetWQVKkg8xCJOvV9Xvt+I322Es2tf9rXwvcPFQ84ta2d62PFp+RJsky4BzgQPH\n6esAcF6rO9qXJGkMZnOVV4B7gO9V1W8OrdoBTF91tQV4ZKh8c7tyay2Dk+9Pt8Nj7ya5vPV53Uib\n6b6uAR5vs56dwIYk57eT8RuAnW3dE63u6PYlSWMwm3MoPwv8PeD5JM+1sn8O3A48lOR64FXgWoCq\neiHJQ8CLDK4Qu6mq3m/tbgTuA84GHm0vGATWA0mmgIMMrhKjqg4m+TLwTKv3pao62JZvBrYn+Qrw\nbOtDkjQmMwZKVf0PIMdYfeUx2vwG8BtHKZ8Efuoo5f8X+JvH6Ote4N6jlP+QwaXEkqQFwDvlJUld\n+Ph6SQve6CXgPs5+YXKGIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknq\nwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwr/YKGnRGf0LjuBfcVwInKFI\nkrowUCRJXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sJAkSR1YaBIkrrwTnmdNo5297Sk+TPjDCXJ\nvUn2J/nuUNnyJLuS7Glfzx9ad2uSqSQvJdk4VH5ZkufbujuTpJWfleTBVv5UkjVDbba0bexJsmWo\nfG2rO9Xanjn3XSFJmovZHPK6D9g0UnYL8FhVrQMea+9JcimwGfh4a/PVJGe0NncBnwPWtdd0n9cD\nh6rqY8AdwNbW13LgNuCTwHrgtqHg2grc0docan1IksZoxkCpqj8EDo4UXwVsa8vbgKuHyrdX1XtV\n9TIwBaxPsgo4p6qerKoC7h9pM93Xw8CVbfayEdhVVQer6hCwC9jU1l3R6o5uX5I0Jid7Un5lVe1r\ny28AK9vyauC1oXqvt7LVbXm0/Ig2VXUYeAe44Dh9XQC83eqO9iVJGpM5n5SvqkpSPQZzKiS5AbgB\n4JJLLhnzaCSdKj7SfvxOdobyZjuMRfu6v5XvBS4eqndRK9vblkfLj2iTZBlwLnDgOH0dAM5rdUf7\n+nOq6u6qmqiqiRUrVpzgP1OSNFsnGyg7gOmrrrYAjwyVb25Xbq1lcPL96XZ47N0kl7dzINeNtJnu\n6xrg8XaeZSewIcn57WT8BmBnW/dEqzu6fUnSmMx4yCvJ7wK/AFyY5HUGV17dDjyU5HrgVeBagKp6\nIclDwIvAYeCmqnq/dXUjgyvGzgYebS+Ae4AHkkwxOPm/ufV1MMmXgWdavS9V1fTFATcD25N8BXi2\n9SFJGqMMfuFfGiYmJmpycnLcw9Ap4o2NGuU5lD6S7K6qiZnq+egVSVIXBookqQsDRZLUhYEiSerC\nQJEkdWGgSJK68O+haFHyEmFp4TFQJJ22fL7X/PKQlySpCwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJ\nXXjZsBY87zmRFgcDRdKSMvoLivel9GOgaEFxNiItXp5DkSR1YaBIkrrwkJfGykNc0unDQNG8MTyk\n05uBImlJ84nE/RgoOiWcjUhLj4GiOTM8JIGBohkYFpJmy0DRBwwPacDzKifHQDkNGQySxsEbGyVJ\nXThDWeScjUhaKAyURcTwkMbH8yoz85CXJKmLRT1DSbIJ+HfAGcBvV9XtYx7SSXP2IS0+/m2VIy3a\nQElyBvAfgF8EXgeeSbKjql4c78hmZnhIp6elflhs0QYKsB6YqqofAiTZDlwFLPhAkbR0LKWQWcyB\nshp4bej968AnxzSWDzj7kDSTk/2cWOhBtJgDZVaS3ADc0N7+SZKXTrKrC4Ef9RlVV47rxDiuE+O4\nTswpHVe2nnTTuY7rL82m0mIOlL3AxUPvL2plR6iqu4G757qxJJNVNTHXfnpzXCfGcZ0Yx3Vilvq4\nFvNlw88A65KsTXImsBnYMeYxSdKStWhnKFV1OMk/AnYyuGz43qp6YczDkqQla9EGCkBVfRP45jxt\nbs6HzU4Rx3ViHNeJcVwnZkmPK1U1H9uRJJ3mFvM5FEnSAmKgzCDJpiQvJZlKcss8bO/iJE8keTHJ\nC0k+38p/PcneJM+112eG2tzaxvdSko1D5Zcleb6tuzNJ5ji2V1p/zyWZbGXLk+xKsqd9PX8+x5Xk\nLw/tk+eSvJvkC+PYX0nuTbI/yXeHyrrtnyRnJXmwlT+VZM0cxvWvk3w/yR8l+UaS81r5miT/Z2i/\nfe1Ujes4Y+v2veu8zx4cGtMrSZ6bz32WY382jP1n7ANV5esYLwYn+38AfBQ4E/gOcOkp3uYq4Gfa\n8keA/wVcCvw68M+OUv/SNq6zgLVtvGe0dU8DlwMBHgU+PcexvQJcOFL2r4Bb2vItwNb5HtfI9+sN\nBtfMz/v+An4e+Bngu6di/wA3Al9ry5uBB+cwrg3Asra8dWhca4brjfTTdVzHGVu3713PfTay/t8A\n/2I+9xnH/mwY+8/Y9MsZyvF98HiXqvp/wPTjXU6ZqtpXVd9uy38MfI/BUwGO5Spge1W9V1UvA1PA\n+iSrgHOq6ska/HTcD1x9CoZ8FbCtLW8b2sY4xnUl8IOqenWG8Z6ScVXVHwIHj7K9XvtnuK+HgStn\nM4s62riq6r9W1eH29kkG93Ed06kY17HGdhxj3WfTWvtrgd89Xh+9x3Wcz4ax/4xNM1CO72iPdzne\nh3tXbbr508BTregft0MU9w5Na481xtVtebR8Lgr4VpLdGTyBAGBlVe1ry28AK8cwrmmbOfI/+bj3\nF/TdPx+0aWHwDnBBhzH+Awa/pU5b2w7d/PckPze07fkcV6/v3akY288Bb1bVnqGyed1nI58NC+Zn\nzEBZoJL8BPB7wBeq6l3gLgaH3j4B7GMw5Z5vn6qqTwCfBm5K8vPDK9tvO2O5bDCDm1s/C/ynVrQQ\n9tcRxrl/jiXJF4HDwNdb0T7gkvZ9/ifA7yQ5Z56HteC+dyP+Nkf+4jKv++wonw0fGPfPmIFyfLN6\nvEtvST7E4Afm61X1+wBV9WZVvV9Vfwb8FoPDcccb416OPIwx57FX1d72dT/wjTaGN9sUenqKv3++\nx9V8Gvh2Vb3Zxjj2/dX03D8ftEmyDDgXOHCyA0vy94FfBv5O+yCiHR450JZ3Mzju/pPzOa7O37ve\n+2wZ8CvAg0Pjnbd9drTPBhbQz5iBcnzz/niXdrzyHuB7VfWbQ+Wrhqr9DWD66pMdwOZ2dcZaYB3w\ndJsCv5vk8tbndcAjcxjXh5N8ZHqZwUnd77btb2nVtgxtY17GNeSI3xrHvb+G9Nw/w31dAzw+HQQn\nKoM/TvdrwGer6n8Pla/I4G8NkeSjbVw/nK9xte32/N51HRvw14DvV9UHh4zma58d67OBhfQzdiJn\n8JfiC/gMg6spfgB8cR629ykGU9Y/Ap5rr88ADwDPt/IdwKqhNl9s43uJoSuTgAkG/xl/APx72o2s\nJzmujzK4YuQ7wAvT+4LB8dXHgD3At4Dl8zmu1t+HGfwWde5Q2bzvLwaBtg/4UwbHpa/vuX+AH2dw\nSG+KwVU6H53DuKYYHCuf/hmbvrLnV9v39zng28BfP1XjOs7Yun3veu6zVn4f8A9H6s7LPuPYnw1j\n/xmbfnmnvCSpCw95SZK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdfH/AW06crn5\nvnRBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1140a2860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# clthresh/cuthresh example:\n",
    "p1=plt.hist(cd[np.nonzero(cd)], 64,[0,20000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAD8CAYAAABQFVIjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFq9JREFUeJzt3X+s39V93/HnqxAoSgPYxrMsA7OjuJPIpKXlyiA1rTpY\njZN2gXUIedmGp1lBE2xKtE0NLNKcJvkj3rSyoS1UrEQYFmo82girKmOOYes/48d1ICWQMN+UILAM\npraBVptYTd/743su+fqLf1z7Ht/vvb7Ph/TV9/M9n3PO5/hzr+/rns+vm6pCkqTZ+qlxD0CSdHYw\nUCRJXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sJAkSR1YaBIkro4d9wDmEuXXHJJrV69etzDkKQF\nZc+ePX9aVctPVm9RBcrq1auZnJwc9zAkaUFJ8spM6nnIS5LUhYEiSerCQJEkdWGgSJK6MFAkSV0Y\nKJKkLgwUSVIXBookqQsDRZLUxaK6U17q6sF8sOyzNffjkOYJZyiSpC4MFElSFwaKJKmLGQVKkh8n\neT7Jc0kmW9nSJLuS7G3vS4bq35FkKslLSa4bKr+y9TOV5K4kaeXnJ3molT+VZPVQm01tG3uTbBoq\nX9PqTrW2581+d0iSTtepzFD+ZlV9oqom2ufbgd1VtRbY3T6T5ApgI/BxYAPwjSTntDZ3A58D1rbX\nhla+GThcVR8D7gS2tr6WAluAq4B1wJah4NoK3NnaHG59SGfGg/ngS9JRZnPI63pgW1veBtwwVL69\nqt6tqpeBKWBdkpXAhVX1ZFUVcP9Im+m+HgaubbOX64BdVXWoqg4Du4ANbd01re7o9iVJYzDTQCng\nO0n2JLmlla2oqv1t+XVgRVteBbw61Pa1VraqLY+WH9Wmqo4AbwPLTtDXMuCtVne0L0nSGMz0PpRP\nVtW+JH8F2JXkh8Mrq6qSzMsL8FsA3gJw+eWXj3k0knT2mlGgVNW+9n4gybcZnM94I8nKqtrfDmcd\naNX3AZcNNb+0le1ry6Plw21eS3IucBFwsJX/8kib/9HWXZzk3DZLGe5rdOz3APcATExMzMvQ01nE\nmx21iJ30kFeSDyf5yPQysB74PrATmL7qahPwSFveCWxsV26tYXDy/el2eOydJFe3cyA3j7SZ7utG\n4PF2nuUxYH2SJe1k/HrgsbbuiVZ3dPuSpDGYyQxlBfDtdoXvucCDVfXfkjwD7EiyGXgFuAmgql5I\nsgN4ETgC3FZV77W+bgXuAy4AHm0vgHuBB5JMAYcYXCVGVR1K8lXgmVbvK1V1qC1/Edie5GvAs60P\nSdKYZPDL/uIwMTFRk5OT4x6GFqLZXCbsIS8tcEn2DN0yclw+HFI6Fu8zkU6Zj16RJHVhoEiSujBQ\nJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwhsbpTNt9CbJM33nvA+o1Jg4Q5EkdeEMRVrIfESM\n5hEDRVpIDBDNYwaKNNc8x6GzlIEiLQaGmOaAJ+UlSV0YKJKkLgwUSVIXBookqQtPykteiit1YaBI\n85VBpwXGQ16SpC6coUjzwTyZjeQ3jx5HbfFeFc2cMxRJUhcGiiSpCwNFktSFgSJJ6sJAkSR14VVe\nko5r9Kov8MovHZ8zFElSFwaKJKmLGQdKknOSPJvkD9rnpUl2Jdnb3pcM1b0jyVSSl5JcN1R+ZZLn\n27q7kqSVn5/koVb+VJLVQ202tW3sTbJpqHxNqzvV2p43u10hSZqNU5mhfB74wdDn24HdVbUW2N0+\nk+QKYCPwcWAD8I0k57Q2dwOfA9a214ZWvhk4XFUfA+4Etra+lgJbgKuAdcCWoeDaCtzZ2hxufUiS\nxmRGgZLkUuBXgd8ZKr4e2NaWtwE3DJVvr6p3q+plYApYl2QlcGFVPVlVBdw/0ma6r4eBa9vs5Tpg\nV1UdqqrDwC5gQ1t3Tas7un1J0hjMdIby74HfAP5yqGxFVe1vy68DK9ryKuDVoXqvtbJVbXm0/Kg2\nVXUEeBtYdoK+lgFvtbqjfUmSxuCkgZLk14ADVbXneHXajGNeXkuY5JYkk0km33zzzXEPR5LOWjOZ\nofwC8JkkPwa2A9ck+S/AG+0wFu39QKu/D7hsqP2lrWxfWx4tP6pNknOBi4CDJ+jrIHBxqzva11Gq\n6p6qmqiqieXLl8/gnytJOh0nDZSquqOqLq2q1QxOtj9eVf8A2AlMX3W1CXikLe8ENrYrt9YwOPn+\ndDs89k6Sq9s5kJtH2kz3dWPbRgGPAeuTLGkn49cDj7V1T7S6o9uXJI3BbO6U/zqwI8lm4BXgJoCq\neiHJDuBF4AhwW1W919rcCtwHXAA82l4A9wIPJJkCDjEILqrqUJKvAs+0el+pqkNt+YvA9iRfA55t\nfUiSxiSDX/YXh4mJiZqcnBz3MDTfzJM/bjXnPvvB//vHetTKKB+9svgk2VNVEyer553ykqQuDBRJ\nUhcGiiSpCx9fL+mU+Eh7HY8zFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4M\nFElSFwaKJKkLA0WS1IXP8pI0az7fS+AMRZLUiYEiSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIX\nBookqQsDRZLUhYEiSerCQJEkdeGzvLS4PPjBZ05J6sMZiiSpC2coks6I0ScQ+/Ths58zFElSFycN\nlCQ/neTpJN9L8kKS32zlS5PsSrK3vS8ZanNHkqkkLyW5bqj8yiTPt3V3JUkrPz/JQ638qSSrh9ps\natvYm2TTUPmaVneqtT2vzy6RJJ2OmcxQ3gWuqaq/AXwC2JDkauB2YHdVrQV2t88kuQLYCHwc2AB8\nI8k5ra+7gc8Ba9trQyvfDByuqo8BdwJbW19LgS3AVcA6YMtQcG0F7mxtDrc+JEljctJAqYE/bx8/\n1F4FXA9sa+XbgBva8vXA9qp6t6peBqaAdUlWAhdW1ZNVVcD9I22m+3oYuLbNXq4DdlXVoao6DOxi\nEGgBrml1R7cvSRqDGZ1DSXJOkueAAwx+wD8FrKiq/a3K68CKtrwKeHWo+WutbFVbHi0/qk1VHQHe\nBpadoK9lwFut7mhfo2O/Jclkksk333xzJv9cSdJpmFGgVNV7VfUJ4FIGs42/PrK+GMxa5p2quqeq\nJqpqYvny5eMejiSdtU7pKq+qegt4gsG5jzfaYSza+4FWbR9w2VCzS1vZvrY8Wn5UmyTnAhcBB0/Q\n10Hg4lZ3tC9J0hjM5Cqv5UkubssXAL8C/BDYCUxfdbUJeKQt7wQ2tiu31jA4+f50Ozz2TpKr2zmQ\nm0faTPd1I/B4m/U8BqxPsqSdjF8PPNbWPdHqjm5fkjQGM7mxcSWwrV2p9VPAjqr6gyT/C9iRZDPw\nCnATQFW9kGQH8CJwBLitqt5rfd0K3AdcADzaXgD3Ag8kmQIOMbhKjKo6lOSrwDOt3leq6lBb/iKw\nPcnXgGdbH5KkMcngl/3FYWJioiYnJ8c9DI2Tz/L6ic9+8P/+6N3tPXmn/MKVZE9VTZysno9ekTQn\njhVWhszZxUevSJK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdeGNjdIidSbvitfi\n5AxFktSFgSJJ6sJDXpLGxud7nV2coUiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCR\nJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLnzYsaV7xCcQLlzMUSVIXBookqQsDRZLUhYEiSeri\npIGS5LIkTyR5MckLST7fypcm2ZVkb3tfMtTmjiRTSV5Kct1Q+ZVJnm/r7kqSVn5+koda+VNJVg+1\n2dS2sTfJpqHyNa3uVGt7Xp9dIkk6HTOZoRwB/kVVXQFcDdyW5ArgdmB3Va0FdrfPtHUbgY8DG4Bv\nJDmn9XU38DlgbXttaOWbgcNV9THgTmBr62spsAW4ClgHbBkKrq3Ana3N4daHJGlMThooVbW/qr7b\nlv8M+AGwCrge2NaqbQNuaMvXA9ur6t2qehmYAtYlWQlcWFVPVlUB94+0me7rYeDaNnu5DthVVYeq\n6jCwC9jQ1l3T6o5uX5I0Bqd0H0o7FPVzwFPAiqra31a9Dqxoy6uAJ4eavdbK/qItj5ZPt3kVoKqO\nJHkbWDZcPtJmGfBWVR05Rl+jY74FuAXg8ssvn/G/VWeBBz94P4OkM2fGJ+WT/Azwe8AXquqd4XVt\nxjEv7zyqqnuqaqKqJpYvXz7u4UjSWWtGgZLkQwzC5FtV9fut+I12GIv2fqCV7wMuG2p+aSvb15ZH\ny49qk+Rc4CLg4An6Oghc3OqO9iVJGoOZXOUV4F7gB1X1W0OrdgLTV11tAh4ZKt/Yrtxaw+Dk+9Pt\n8Ng7Sa5ufd480ma6rxuBx9us5zFgfZIl7WT8euCxtu6JVnd0+5KkMZjJOZRfAP4h8HyS51rZvwK+\nDuxIshl4BbgJoKpeSLIDeJHBFWK3VdV7rd2twH3ABcCj7QWDwHogyRRwiMFVYlTVoSRfBZ5p9b5S\nVYfa8heB7Um+Bjzb+pAkjUkGv+wvDhMTEzU5OTnuYWiueFL+hLJ33COYOR8OOV5J9lTVxMnqeae8\nJKkLH18vad4bfaS9M5b5yRmKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaK\nJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkL/2KjpAVn9C84gn/FcT5w\nhiJJ6sJAkSR1YaBIkrowUCRJXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sI75XX2ePCDd09Lmjsn\nnaEk+WaSA0m+P1S2NMmuJHvb+5KhdXckmUryUpLrhsqvTPJ8W3dXkrTy85M81MqfSrJ6qM2mto29\nSTYNla9pdada2/NmvyskSbMxk0Ne9wEbRspuB3ZX1Vpgd/tMkiuAjcDHW5tvJDmntbkb+Bywtr2m\n+9wMHK6qjwF3AltbX0uBLcBVwDpgy1BwbQXubG0Otz4kSWN00kCpqj8CDo0UXw9sa8vbgBuGyrdX\n1btV9TIwBaxLshK4sKqerKoC7h9pM93Xw8C1bfZyHbCrqg5V1WFgF7Chrbum1R3dviRpTE73pPyK\nqtrfll8HVrTlVcCrQ/Vea2Wr2vJo+VFtquoI8Daw7AR9LQPeanVH+5IkjcmsT8pXVSWZt8+NTnIL\ncAvA5ZdfPubRSDpTfKT9+J3uDOWNdhiL9n6gle8DLhuqd2kr29eWR8uPapPkXOAi4OAJ+joIXNzq\njvb1AVV1T1VNVNXE8uXLT/GfKUmaqdMNlJ3A9FVXm4BHhso3tiu31jA4+f50Ozz2TpKr2zmQm0fa\nTPd1I/B4O8/yGLA+yZJ2Mn498Fhb90SrO7p9SdKYnPSQV5LfBX4ZuCTJawyuvPo6sCPJZuAV4CaA\nqnohyQ7gReAIcFtVvde6upXBFWMXAI+2F8C9wANJphic/N/Y+jqU5KvAM63eV6pq+uKALwLbk3wN\neLb1IUkaowx+4V8cJiYmanJyctzD0JnijY2nJHvHPYIzz3MofSTZU1UTJ6vno1ckSV0YKJKkLgwU\nSVIXBookqQsDRZLUhYEiSerCv4eihclLhKV5x0CRdNby+V5zy0NekqQuDBRJUhcGiiSpCwNFktSF\ngSJJ6sJAkSR14WXDmv+850RaEAwUSYvK6L0p3pfSj4Gi+cXZiLRgeQ5FktSFgSJJ6sJDXhovD3FJ\nZw0DRXPH8JDOagaKpEXNJxL3Y6DozHA2Ii06Bopmz/CQhIGikzEsJM2QgaKfMDwkwPMqp8tAORsZ\nDJLGwBsbJUldOENZ6JyNSJonDJSFxPCQxsbzKifnIS9JUhcLeoaSZAPwH4BzgN+pqq+PeUinz9mH\ntOD4t1WOtmADJck5wH8CfgV4DXgmyc6qenG8I5sBw0M6Ky32w2ILNlCAdcBUVf0JQJLtwPXA/A8U\nSYvGYgqZhRwoq4BXhz6/Blw1prH8hLMPSSdxrJCZifkeRAs5UGYkyS3ALe3jnyd56TS7ugT40z6j\n6spxnRrHdWoc16k5o+PKl0/7F9bZjuuvzqTSQg6UfcBlQ58vbWVHqap7gHtmu7Ekk1U1Mdt+enNc\np8ZxnRrHdWoW+7gW8mXDzwBrk6xJch6wEdg55jFJ0qK1YGcoVXUkyT8FHmNw2fA3q+qFMQ9Lkhat\nBRsoAFX1h8AfztHmZn3Y7AxxXKfGcZ0ax3VqFvW4UjW/rxqQJC0MC/kciiRpHjFQTiLJhiQvJZlK\ncvscbO+yJE8keTHJC0k+38q/nGRfkufa69NDbe5o43spyXVD5Vcmeb6tuyvJrG6SSfLj1t9zSSZb\n2dIku5Lsbe9L5nJcSf7a0D55Lsk7Sb4wjv2V5JtJDiT5/lBZt/2T5PwkD7Xyp5KsnsW4/m2SHyb5\n4yTfTnJxK1+d5P8O7bffPlPjOsHYun3tOu+zh4bG9OMkz83lPsvxfzaM/XvsfVXl6zgvBif7fwR8\nFDgP+B5wxRne5krg59vyR4D/DVwBfBn4l8eof0Ub1/nAmjbec9q6p4GrgQCPAp+a5dh+DFwyUvZv\ngNvb8u3A1rke18jX63UG18zP+f4Cfgn4eeD7Z2L/ALcCv92WNwIPzWJc64Fz2/LWoXGtHq430k/X\ncZ1gbN2+dj332cj6fwf867ncZxz/Z8PYv8emX85QTuz9x7tU1f8Dph/vcsZU1f6q+m5b/jPgBwye\nCnA81wPbq+rdqnoZmALWJVkJXFhVT9bgu+N+4IYzMOTrgW1tedvQNsYxrmuBH1XVKycZ7xkZV1X9\nEXDoGNvrtX+G+3oYuHYms6hjjauq/ntVHWkfn2RwH9dxnYlxHW9sJzDWfTattb8J+N0T9dF7XCf4\n2TD277FpBsqJHevxLif64d5Vm27+HPBUK/pn7RDFN4emtccb46q2PFo+GwV8J8meDJ5AALCiqva3\n5deBFWMY17SNHP2ffNz7C/run/fbtDB4G1jWYYz/mMFvqdPWtEM3/zPJLw5tey7H1etrdybG9ovA\nG1W1d6hsTvfZyM+GefM9ZqDMU0l+Bvg94AtV9Q5wN4NDb58A9jOYcs+1T1bVJ4BPAbcl+aXhle23\nnbFcNpjBza2fAf5rK5oP++so49w/x5PkS8AR4FutaD9wefs6/3PgwSQXzvGw5t3XbsTf4+hfXOZ0\nnx3jZ8P7xv09ZqCc2Iwe79Jbkg8x+Ib5VlX9PkBVvVFV71XVXwL/mcHhuBONcR9HH8aY9diral97\nPwB8u43hjTaFnp7iH5jrcTWfAr5bVW+0MY59fzU998/7bZKcC1wEHDzdgSX5R8CvAX+//SCiHR45\n2Jb3MDju/rNzOa7OX7ve++xc4NeBh4bGO2f77Fg/G5hH32MGyonN+eNd2vHKe4EfVNVvDZWvHKr2\nd4Dpq092Ahvb1RlrgLXA020K/E6Sq1ufNwOPzGJcH07ykellBid1v9+2v6lV2zS0jTkZ15Cjfmsc\n9/4a0nP/DPd1I/D4dBCcqgz+ON1vAJ+pqv8zVL48g781RJKPtnH9yVyNq22359eu69iAvwX8sKre\nP2Q0V/vseD8bmE/fY6dyBn8xvoBPM7ia4kfAl+Zge59kMGX9Y+C59vo08ADwfCvfCawcavOlNr6X\nGLoyCZhg8J/xR8B/pN3Ieprj+iiDK0a+B7wwvS8YHF/dDewFvgMsnctxtf4+zOC3qIuGyuZ8fzEI\ntP3AXzA4Lr255/4BfprBIb0pBlfpfHQW45picKx8+nts+sqev9u+vs8B3wX+9pka1wnG1u1r13Of\ntfL7gH8yUndO9hnH/9kw9u+x6Zd3ykuSuvCQlySpCwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJXRgo\nkqQuDBRJUhf/H4OJTE49nxD2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1148cce80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tl12000_cd=clthresh(cd,12000)\n",
    "tu12000_cd=cuthresh(cd,12000)\n",
    "p1=plt.hist(tu12000_cd[np.nonzero(tu12000_cd)], 64, [0,20000],color='Orange')\n",
    "p1=plt.hist(tl12000_cd[np.nonzero(tl12000_cd)], 64, [0,20000],color='Green')\n",
    "output_cifti=nhw+'/thresh_demo.dtseries.nii'\n",
    "csave(output_cifti,tl12000_cd,ch,cn,cx,cf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cmask(target,mask):\n",
    "    '''\n",
    "    Masks target cifti file with non-zero values of mask cifti file\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    target: numpy ndarray\n",
    "        numpy array containing cifti data to be masked\n",
    "    mask: numpy ndarray\n",
    "        numpy array containing cifti data to use as mask. Can be 1D or 2D, must match target spatial dimensions.\n",
    "    '''\n",
    "    masked_data=target.copy()\n",
    "    if target.shape != mask.shape:\n",
    "        mask=np.tile(mask,(target.shape[0],1))\n",
    "    masked_data[mask==0]=0 # NB could be sped up by using a data dictionary, not currently high-priority\n",
    "    return masked_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 305,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAD8CAYAAABQFVIjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFqBJREFUeJzt3X+s39V93/HnqzhQlIYfBs/yDMyO4k4ilZaWKwepadXB\najtpF1jLmKdteBoKmmBTom0qsEijS1IJb1rZ0BYiWhCGpTWMNsLqwjwH2KpJ48d1ICWQMN8EEFgG\nE9tAq22spu/98T0Xff2t7Xvte3y/9/o+H9JX9/M9n3PO5/hzr7+vez6/bqoKSZLm6sfGPQBJ0unB\nQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSepi2bgHMJ8uvPDCWrNmzbiHIUmL\nyu7du39UVStmqrekAmXNmjVMTk6OexiStKgkeXU29TzkJUnqwkCRJHVhoEiSujBQJEldGCiSpC4M\nFElSFwaKJKkLA0WS1IWBIknqYkndKS/1tOaW//znyl65/ZfGMBJpYXCGIknqwkCRJHVhoEiSuphV\noCR5JcnzSZ5LMtnKlifZlWRP+3r+UP1bk0wleSnJxqHyy1o/U0nuTJJWflaSB1v5U0nWDLXZ0rax\nJ8mWofK1re5Ua3vm3HeHJOlknchJ+b9aVT8aen8L8FhV3Z7klvb+5iSXApuBjwN/EfhWkp+sqveB\nu4DPAU8B3wQ2AY8C1wOHqupjSTYDW4G/lWQ5cBswARSwO8mOqjrU6txRVduTfK31cddJ7gfpuI52\nAl7SkeZyyOsqYFtb3gZcPVS+vareq6qXgSlgfZJVwDlV9WRVFXD/SJvpvh4Grmyzl43Arqo62EJk\nF7Cprbui1R3dviRpDGYbKMVgprE7yQ2tbGVV7WvLbwAr2/Jq4LWhtq+3stVtebT8iDZVdRh4B7jg\nOH1dALzd6o72JUkag9ke8vpUVe1N8heAXUm+P7yyqipJ9R/e3LUAvAHgkksuGfNoJOn0NatAqaq9\n7ev+JN8A1gNvJllVVfva4az9rfpe4OKh5he1sr1tebR8uM3rSZYB5wIHWvkvjLT5b23deUmWtVnK\ncF+jY78buBtgYmJiQYaeTh/e7KilbMZDXkk+nOQj08vABuC7wA5g+qqrLcAjbXkHsLldubUWWAc8\n3Q6PvZvk8nYO5LqRNtN9XQM83s6z7AQ2JDm/XUW2AdjZ1j3R6o5uX5I0BrOZoawEvtGu8F0G/E5V\n/ZckzwAPJbkeeBW4FqCqXkjyEPAicBi4qV3hBXAjcB9wNoOrux5t5fcADySZAg4yuEqMqjqY5MvA\nM63el6rqYFu+Gdie5CvAs60PSdKYZPDL/tIwMTFRk5OT4x6GFqG5XDbsIS8tdkl2V9XETPV8OKR0\nFN53Ip04H70iSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhTc2SqfY6E2Sp/rO\neR9QqXFxhiJJ6sIZirSI+YgYLSQGirSIGCBayAwUaZ55jkOnKwNFWgIMMc0HT8pLkrowUCRJXRgo\nkqQuPIciLVGeV1FvBoqWPC/FlfowUKQFyqDTYuM5FElSF85QpAVgocxG5vtBljq9OEORJHVhoEiS\nujBQJEldGCiSpC4MFElSF17lJemYvJteJ8IZiiSpCwNFktTFrAMlyRlJnk3yB+398iS7kuxpX88f\nqntrkqkkLyXZOFR+WZLn27o7k6SVn5XkwVb+VJI1Q222tG3sSbJlqHxtqzvV2p45t10hSZqLE5mh\nfB743tD7W4DHqmod8Fh7T5JLgc3Ax4FNwFeTnNHa3AV8DljXXpta+fXAoar6GHAHsLX1tRy4Dfgk\nsB64bSi4tgJ3tDaHWh+SpDGZVaAkuQj4JeC3h4qvAra15W3A1UPl26vqvap6GZgC1idZBZxTVU9W\nVQH3j7SZ7uth4Mo2e9kI7Kqqg1V1CNgFbGrrrmh1R7cvSRqD2c5Q/i3wa8CfDZWtrKp9bfkNYGVb\nXg28NlTv9Va2ui2Plh/RpqoOA+8AFxynrwuAt1vd0b4kSWMwY6Ak+WVgf1XtPladNuOongPrJckN\nSSaTTL711lvjHo4knbZmM0P5WeCzSV4BtgNXJPmPwJvtMBbt6/5Wfy9w8VD7i1rZ3rY8Wn5EmyTL\ngHOBA8fp6wBwXqs72tcRquruqpqoqokVK1bM4p8rSToZMwZKVd1aVRdV1RoGJ9sfr6q/C+wApq+6\n2gI80pZ3AJvblVtrGZx8f7odHns3yeXtHMh1I22m+7qmbaOAncCGJOe3k/EbgJ1t3ROt7uj2JUlj\nMJc75W8HHkpyPfAqcC1AVb2Q5CHgReAwcFNVvd/a3AjcB5wNPNpeAPcADySZAg4yCC6q6mCSLwPP\ntHpfqqqDbflmYHuSrwDPtj4kSWOSwS/7S8PExERNTk6OexhaYBbKH7daLHz0ytKTZHdVTcxUzzvl\nJUldGCiSpC4MFElSFz6+XtIJ8ZH2OhZnKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdWGgSJK6\nMFAkSV0YKJKkLgwUSVIXBookqQuf5SVpzny+l8AZiiSpEwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJ\nXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sJneWlJOdozpyT14QxFktSFMxRJp8TobNCnD5/+nKFI\nkrqYMVCS/HiSp5N8J8kLSf5lK1+eZFeSPe3r+UNtbk0yleSlJBuHyi9L8nxbd2eStPKzkjzYyp9K\nsmaozZa2jT1JtgyVr211p1rbM/vsEknSyZjNDOU94Iqq+ivAJ4BNSS4HbgEeq6p1wGPtPUkuBTYD\nHwc2AV9Nckbr6y7gc8C69trUyq8HDlXVx4A7gK2tr+XAbcAngfXAbUPBtRW4o7U51PqQJI3JjIFS\nA3/S3n6ovQq4CtjWyrcBV7flq4DtVfVeVb0MTAHrk6wCzqmqJ6uqgPtH2kz39TBwZZu9bAR2VdXB\nqjoE7GIQaAGuaHVHty9JGoNZnUNJckaS54D9DD7gnwJWVtW+VuUNYGVbXg28NtT89Va2ui2Plh/R\npqoOA+8AFxynrwuAt1vd0b5Gx35Dkskkk2+99dZs/rmSpJMwq0Cpqver6hPARQxmGz81sr4YzFoW\nnKq6u6omqmpixYoV4x6OJJ22Tugqr6p6G3iCwbmPN9thLNrX/a3aXuDioWYXtbK9bXm0/Ig2SZYB\n5wIHjtPXAeC8Vne0L0nSGMzmKq8VSc5ry2cDvwh8H9gBTF91tQV4pC3vADa3K7fWMjj5/nQ7PPZu\nksvbOZDrRtpM93UN8Hib9ewENiQ5v52M3wDsbOueaHVHty9JGoPZ3Ni4CtjWrtT6MeChqvqDJP8T\neCjJ9cCrwLUAVfVCkoeAF4HDwE1V9X7r60bgPuBs4NH2ArgHeCDJFHCQwVViVNXBJF8Gnmn1vlRV\nB9vyzcD2JF8Bnm19SJLGJINf9peGiYmJmpycHPcwNEY+y2t8vFN+8Uqyu6omZqrno1ckzYujhbkh\nc3rx0SuSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkL\nn+UlaWx8vtfpxRmKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS\n1IWBIknqwkCRJHVhoEiSuvBpw5IWFJ9AvHg5Q5EkdWGgSJK6MFAkSV0YKJKkLmYMlCQXJ3kiyYtJ\nXkjy+Va+PMmuJHva1/OH2tyaZCrJS0k2DpVfluT5tu7OJGnlZyV5sJU/lWTNUJstbRt7kmwZKl/b\n6k61tmf22SWSpJMxmxnKYeCfVtWlwOXATUkuBW4BHquqdcBj7T1t3Wbg48Am4KtJzmh93QV8DljX\nXpta+fXAoar6GHAHsLX1tRy4DfgksB64bSi4tgJ3tDaHWh+SpDGZMVCqal9Vfbst/zHwPWA1cBWw\nrVXbBlzdlq8CtlfVe1X1MjAFrE+yCjinqp6sqgLuH2kz3dfDwJVt9rIR2FVVB6vqELAL2NTWXdHq\njm5fkjQGJ3QfSjsU9dPAU8DKqtrXVr0BrGzLq4Enh5q93sr+tC2Plk+3eQ2gqg4neQe4YLh8pM0F\nwNtVdfgofY2O+QbgBoBLLrlk1v9WLX5Hu59B0qkz65PySX4C+D3gC1X17vC6NuOozmProqrurqqJ\nqppYsWLFuIcjSaetWQVKkg8xCJOvV9Xvt+I322Es2tf9rXwvcPFQ84ta2d62PFp+RJsky4BzgQPH\n6esAcF6rO9qXJGkMZnOVV4B7gO9V1W8OrdoBTF91tQV4ZKh8c7tyay2Dk+9Pt8Nj7ya5vPV53Uib\n6b6uAR5vs56dwIYk57eT8RuAnW3dE63u6PYlSWMwm3MoPwv8PeD5JM+1sn8O3A48lOR64FXgWoCq\neiHJQ8CLDK4Qu6mq3m/tbgTuA84GHm0vGATWA0mmgIMMrhKjqg4m+TLwTKv3pao62JZvBrYn+Qrw\nbOtDkjQmMwZKVf0PIMdYfeUx2vwG8BtHKZ8Efuoo5f8X+JvH6Ote4N6jlP+QwaXEkqQFwDvlJUld\n+Ph6SQve6CXgPs5+YXKGIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknq\nwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwr/YKGnRGf0LjuBfcVwInKFI\nkrowUCRJXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sJAkSR1YaBIkrrwTnmdNo5297Sk+TPjDCXJ\nvUn2J/nuUNnyJLuS7Glfzx9ad2uSqSQvJdk4VH5ZkufbujuTpJWfleTBVv5UkjVDbba0bexJsmWo\nfG2rO9Xanjn3XSFJmovZHPK6D9g0UnYL8FhVrQMea+9JcimwGfh4a/PVJGe0NncBnwPWtdd0n9cD\nh6rqY8AdwNbW13LgNuCTwHrgtqHg2grc0docan1IksZoxkCpqj8EDo4UXwVsa8vbgKuHyrdX1XtV\n9TIwBaxPsgo4p6qerKoC7h9pM93Xw8CVbfayEdhVVQer6hCwC9jU1l3R6o5uX5I0Jid7Un5lVe1r\ny28AK9vyauC1oXqvt7LVbXm0/Ig2VXUYeAe44Dh9XQC83eqO9iVJGpM5n5SvqkpSPQZzKiS5AbgB\n4JJLLhnzaCSdKj7SfvxOdobyZjuMRfu6v5XvBS4eqndRK9vblkfLj2iTZBlwLnDgOH0dAM5rdUf7\n+nOq6u6qmqiqiRUrVpzgP1OSNFsnGyg7gOmrrrYAjwyVb25Xbq1lcPL96XZ47N0kl7dzINeNtJnu\n6xrg8XaeZSewIcn57WT8BmBnW/dEqzu6fUnSmMx4yCvJ7wK/AFyY5HUGV17dDjyU5HrgVeBagKp6\nIclDwIvAYeCmqnq/dXUjgyvGzgYebS+Ae4AHkkwxOPm/ufV1MMmXgWdavS9V1fTFATcD25N8BXi2\n9SFJGqMMfuFfGiYmJmpycnLcw9Ap4o2NGuU5lD6S7K6qiZnq+egVSVIXBookqQsDRZLUhYEiSerC\nQJEkdWGgSJK68O+haFHyEmFp4TFQJJ22fL7X/PKQlySpCwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJ\nXXjZsBY87zmRFgcDRdKSMvoLivel9GOgaEFxNiItXp5DkSR1YaBIkrrwkJfGykNc0unDQNG8MTyk\n05uBImlJ84nE/RgoOiWcjUhLj4GiOTM8JIGBohkYFpJmy0DRBwwPacDzKifHQDkNGQySxsEbGyVJ\nXThDWeScjUhaKAyURcTwkMbH8yoz85CXJKmLRT1DSbIJ+HfAGcBvV9XtYx7SSXP2IS0+/m2VIy3a\nQElyBvAfgF8EXgeeSbKjql4c78hmZnhIp6elflhs0QYKsB6YqqofAiTZDlwFLPhAkbR0LKWQWcyB\nshp4bej968AnxzSWDzj7kDSTk/2cWOhBtJgDZVaS3ADc0N7+SZKXTrKrC4Ef9RlVV47rxDiuE+O4\nTswpHVe2nnTTuY7rL82m0mIOlL3AxUPvL2plR6iqu4G757qxJJNVNTHXfnpzXCfGcZ0Yx3Vilvq4\nFvNlw88A65KsTXImsBnYMeYxSdKStWhnKFV1OMk/AnYyuGz43qp6YczDkqQla9EGCkBVfRP45jxt\nbs6HzU4Rx3ViHNeJcVwnZkmPK1U1H9uRJJ3mFvM5FEnSAmKgzCDJpiQvJZlKcss8bO/iJE8keTHJ\nC0k+38p/PcneJM+112eG2tzaxvdSko1D5Zcleb6tuzNJ5ji2V1p/zyWZbGXLk+xKsqd9PX8+x5Xk\nLw/tk+eSvJvkC+PYX0nuTbI/yXeHyrrtnyRnJXmwlT+VZM0cxvWvk3w/yR8l+UaS81r5miT/Z2i/\nfe1Ujes4Y+v2veu8zx4cGtMrSZ6bz32WY382jP1n7ANV5esYLwYn+38AfBQ4E/gOcOkp3uYq4Gfa\n8keA/wVcCvw68M+OUv/SNq6zgLVtvGe0dU8DlwMBHgU+PcexvQJcOFL2r4Bb2vItwNb5HtfI9+sN\nBtfMz/v+An4e+Bngu6di/wA3Al9ry5uBB+cwrg3Asra8dWhca4brjfTTdVzHGVu3713PfTay/t8A\n/2I+9xnH/mwY+8/Y9MsZyvF98HiXqvp/wPTjXU6ZqtpXVd9uy38MfI/BUwGO5Spge1W9V1UvA1PA\n+iSrgHOq6ska/HTcD1x9CoZ8FbCtLW8b2sY4xnUl8IOqenWG8Z6ScVXVHwIHj7K9XvtnuK+HgStn\nM4s62riq6r9W1eH29kkG93Ed06kY17HGdhxj3WfTWvtrgd89Xh+9x3Wcz4ax/4xNM1CO72iPdzne\nh3tXbbr508BTregft0MU9w5Na481xtVtebR8Lgr4VpLdGTyBAGBlVe1ry28AK8cwrmmbOfI/+bj3\nF/TdPx+0aWHwDnBBhzH+Awa/pU5b2w7d/PckPze07fkcV6/v3akY288Bb1bVnqGyed1nI58NC+Zn\nzEBZoJL8BPB7wBeq6l3gLgaH3j4B7GMw5Z5vn6qqTwCfBm5K8vPDK9tvO2O5bDCDm1s/C/ynVrQQ\n9tcRxrl/jiXJF4HDwNdb0T7gkvZ9/ifA7yQ5Z56HteC+dyP+Nkf+4jKv++wonw0fGPfPmIFyfLN6\nvEtvST7E4Afm61X1+wBV9WZVvV9Vfwb8FoPDcccb416OPIwx57FX1d72dT/wjTaGN9sUenqKv3++\nx9V8Gvh2Vb3Zxjj2/dX03D8ftEmyDDgXOHCyA0vy94FfBv5O+yCiHR450JZ3Mzju/pPzOa7O37ve\n+2wZ8CvAg0Pjnbd9drTPBhbQz5iBcnzz/niXdrzyHuB7VfWbQ+Wrhqr9DWD66pMdwOZ2dcZaYB3w\ndJsCv5vk8tbndcAjcxjXh5N8ZHqZwUnd77btb2nVtgxtY17GNeSI3xrHvb+G9Nw/w31dAzw+HQQn\nKoM/TvdrwGer6n8Pla/I4G8NkeSjbVw/nK9xte32/N51HRvw14DvV9UHh4zma58d67OBhfQzdiJn\n8JfiC/gMg6spfgB8cR629ykGU9Y/Ap5rr88ADwDPt/IdwKqhNl9s43uJoSuTgAkG/xl/APx72o2s\nJzmujzK4YuQ7wAvT+4LB8dXHgD3At4Dl8zmu1t+HGfwWde5Q2bzvLwaBtg/4UwbHpa/vuX+AH2dw\nSG+KwVU6H53DuKYYHCuf/hmbvrLnV9v39zng28BfP1XjOs7Yun3veu6zVn4f8A9H6s7LPuPYnw1j\n/xmbfnmnvCSpCw95SZK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdfH/AW06crn5\nvnRBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12039c1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cmask example:\n",
    "plt.hist(cd[np.nonzero(cd)],64,[0,20000]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAD8CAYAAABQFVIjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFlRJREFUeJzt3X+s3fV93/HnqzhQlAawjWdZBmZH8SaRSUvLlUFqWnWw\n2k7aBbYh5KUbnmYFTbAp0TY1sEijS/JHvGllQ1uoWIkwLNQw2girKvMcw9Z/xo/rQEogYb4pQWAZ\nTG0DrTaxmr73x/lccnxiX1/7fnzPvfbzIR2d7/l8P5/P9+PvPT6v+/n+ODdVhSRJc/VT4x6AJOns\nYKBIkrowUCRJXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sJAkSR1sWTcA5hPl156aa1Zs2bcw5Ck\nRWXv3r1/UlUrTlbvnAqUNWvWMDk5Oe5hSNKikuTV2dTzkJckqQsDRZLUhYEiSerCQJEkdWGgSJK6\nMFAkSV0YKJKkLgwUSVIXBookqYtz6k55qauH8pNln635H4e0QDhDkSR1YaBIkrowUCRJXcwqUJL8\nKMkLSZ5PMtnKliXZnWRfe146VP+OJFNJXk6ycaj8qtbPVJK7k6SVX5Dk4Vb+dJI1Q222tG3sS7Jl\nqHxtqzvV2p4/990hSTpdpzJD+RtV9Ymqmmivbwf2VNU6YE97TZIrgc3Ax4FNwNeTnNfa3AN8DljX\nHpta+VbgSFV9DLgL2Nb6WgbcCVwNrAfuHAqubcBdrc2R1od0ZjyUn3xIOsZcDnldD2xvy9uBG4bK\nd1TVe1X1CjAFrE+yCrioqp6qqgIeGGkz3dejwHVt9rIR2F1Vh6vqCLAb2NTWXdvqjm5fkjQGsw2U\nAr6dZG+SW1rZyqo60JbfAFa25dXAa0NtX29lq9vyaPkxbarqKPAOsHyGvpYDb7e6o31JksZgtveh\nfLKq9if5S8DuJD8YXllVlWRBXoDfAvAWgCuuuGLMo5Gks9esAqWq9rfng0m+xeB8xptJVlXVgXY4\n62Crvh+4fKj5Za1sf1seLR9u83qSJcDFwKFW/ksjbf5HW3dJkiVtljLc1+jY7wXuBZiYmFiQoaez\niDc76hx20kNeST6c5CPTy8AG4HvATmD6qqstwGNteSewuV25tZbByfdn2uGxd5Nc086B3DzSZrqv\nG4En2nmWXcCGJEvbyfgNwK627slWd3T7kqQxmM0MZSXwrXaF7xLgoar6b0meBR5JshV4FbgJoKpe\nTPII8BJwFLitqt5vfd0K3A9cCDzeHgD3AQ8mmQIOM7hKjKo6nOQrwLOt3per6nBb/iKwI8lXgeda\nH5KkMcngl/1zw8TERE1OTo57GFqM5nKZsIe8tMgl2Tt0y8gJ+eWQ0vF4n4l0yvzqFUlSFwaKJKkL\nA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC68sVE600ZvkvTOeZ2lnKFIkrowUCRJXRgokqQu\nPIcizTf/CJfOUs5QJEldOEORFhOvGNMC5gxFktSFMxRpMfN8jBYQZyiSpC6coUj+uV+pC2cokqQu\nnKFIC5UzJy0yzlAkSV04Q5EWAmcjOgs4Q5EkdWGgSJK68JCXdLbxZkeNiTMUSVIXBookqQsDRZLU\nhYEiSepi1oGS5LwkzyX5/fZ6WZLdSfa156VDde9IMpXk5SQbh8qvSvJCW3d3krTyC5I83MqfTrJm\nqM2Wto19SbYMla9tdada2/PntiskSXNxKjOUzwPfH3p9O7CnqtYBe9prklwJbAY+DmwCvp7kvNbm\nHuBzwLr22NTKtwJHqupjwF3AttbXMuBO4GpgPXDnUHBtA+5qbY60PiRJYzKrQElyGfArwG8PFV8P\nbG/L24Ebhsp3VNV7VfUKMAWsT7IKuKiqnqqqAh4YaTPd16PAdW32shHYXVWHq+oIsBvY1NZd2+qO\nbl+SNAaznaH8e+DXgb8YKltZVQfa8hvAyra8GnhtqN7rrWx1Wx4tP6ZNVR0F3gGWz9DXcuDtVne0\nL0nSGJw0UJL8KnCwqvaeqE6bcSzIO6eS3JJkMsnkW2+9Ne7hSNJZazYzlJ8HPpPkR8AO4Nok/wV4\nsx3Goj0fbPX3A5cPtb+sle1vy6Plx7RJsgS4GDg0Q1+HgEta3dG+jlFV91bVRFVNrFixYhb/XEnS\n6ThpoFTVHVV1WVWtYXCy/Ymq+vvATmD6qqstwGNteSewuV25tZbByfdn2uGxd5Nc086B3DzSZrqv\nG9s2CtgFbEiytJ2M3wDsauuebHVHty9JGoO5fJfX14BHkmwFXgVuAqiqF5M8ArwEHAVuq6r3W5tb\ngfuBC4HH2wPgPuDBJFPAYQbBRVUdTvIV4NlW78tVdbgtfxHYkeSrwHOtD0nSmGTwy/65YWJioiYn\nJ8c9DC0058LfIvHLITUHSfZW1cTJ6nmnvCSpCwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJXRgokqQu\nDBRJUhcGiiSpCwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ\n6sJAkSR1sWTcA5A0Dx7KT5Z9tuZ/HDqrOUORJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLLxvW\nueV4l89K6sIZiiSpCwNFktSFgSJJ6uKkgZLkp5M8k+S7SV5M8q9b+bIku5Psa89Lh9rckWQqyctJ\nNg6VX5Xkhbbu7iRp5RckebiVP51kzVCbLW0b+5JsGSpf2+pOtbbn99klkqTTMZsZynvAtVX114FP\nAJuSXAPcDuypqnXAnvaaJFcCm4GPA5uAryc5r/V1D/A5YF17bGrlW4EjVfUx4C5gW+trGXAncDWw\nHrhzKLi2AXe1NkdaH5KkMTlpoNTAn7WXH2qPAq4Htrfy7cANbfl6YEdVvVdVrwBTwPokq4CLquqp\nqirggZE20309ClzXZi8bgd1VdbiqjgC7GQRagGtb3dHtS5LGYFbnUJKcl+R54CCDD/ingZVVdaBV\neQNY2ZZXA68NNX+9la1uy6Plx7SpqqPAO8DyGfpaDrzd6o72NTr2W5JMJpl86623ZvPPlSSdhlkF\nSlW9X1WfAC5jMNv4ayPri8GsZcGpqnuraqKqJlasWDHu4UjSWeuUrvKqqreBJxmc+3izHcaiPR9s\n1fYDlw81u6yV7W/Lo+XHtEmyBLgYODRDX4eAS1rd0b4kSWMwm6u8ViS5pC1fCPwy8ANgJzB91dUW\n4LG2vBPY3K7cWsvg5Psz7fDYu0muaedAbh5pM93XjcATbdazC9iQZGk7Gb8B2NXWPdnqjm5fkjQG\ns/nqlVXA9nal1k8Bj1TV7yf5X8AjSbYCrwI3AVTVi0keAV4CjgK3VdX7ra9bgfuBC4HH2wPgPuDB\nJFPAYQZXiVFVh5N8BXi21ftyVR1uy18EdiT5KvBc60OSNCYZ/LJ/bpiYmKjJyclxD0Pj5Hd5/Zh/\nAlizlGRvVU2crJ53ykuSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQ\nJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkL\nA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknq4qSBkuTyJE8m\neSnJi0k+38qXJdmdZF97XjrU5o4kU0leTrJxqPyqJC+0dXcnSSu/IMnDrfzpJGuG2mxp29iXZMtQ\n+dpWd6q1Pb/PLpEknY7ZzFCOAv+8qq4ErgFuS3IlcDuwp6rWAXvaa9q6zcDHgU3A15Oc1/q6B/gc\nsK49NrXyrcCRqvoYcBewrfW1DLgTuBpYD9w5FFzbgLtamyOtD0nSmJw0UKrqQFV9py3/KfB9YDVw\nPbC9VdsO3NCWrwd2VNV7VfUKMAWsT7IKuKiqnqqqAh4YaTPd16PAdW32shHYXVWHq+oIsBvY1NZd\n2+qObl+SNAZLTqVyOxT1s8DTwMqqOtBWvQGsbMurgaeGmr3eyv68LY+WT7d5DaCqjiZ5B1g+XD7S\nZjnwdlUdPU5fo2O+BbgF4Iorrpj1v1VngYcy7hFI55RZn5RP8jPA7wJfqKp3h9e1GUd1HlsXVXVv\nVU1U1cSKFSvGPRxJOmvNKlCSfIhBmHyzqn6vFb/ZDmPRng+28v3A5UPNL2tl+9vyaPkxbZIsAS4G\nDs3Q1yHgklZ3tC9J0hjM5iqvAPcB36+q3xxatROYvupqC/DYUPnmduXWWgYn359ph8feTXJN6/Pm\nkTbTfd0IPNFmPbuADUmWtpPxG4Bdbd2Tre7o9iVJYzCbcyg/D/wD4IUkz7eyfwl8DXgkyVbgVeAm\ngKp6MckjwEsMrhC7rareb+1uBe4HLgQebw8YBNaDSaaAwwyuEqOqDif5CvBsq/flqjrclr8I7Ejy\nVeC51ockaUwy+GX/3DAxMVGTk5PjHobmiyflZ/bZc+f/vuYmyd6qmjhZPe+UlyR1YaBIkrowUCRJ\nXZzSjY2SziLHO8fkeRXNgTMUSVIXBookqQsDRZLUhYEiSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwU\nSVIXBookqQsDRZLUhYEiSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerC\nQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXS8Y9AKmbhzLuEUjntJPOUJJ8I8nBJN8bKluWZHeSfe15\n6dC6O5JMJXk5ycah8quSvNDW3Z0krfyCJA+38qeTrBlqs6VtY1+SLUPla1vdqdb2/LnvCknSXMzm\nkNf9wKaRstuBPVW1DtjTXpPkSmAz8PHW5utJzmtt7gE+B6xrj+k+twJHqupjwF3AttbXMuBO4Gpg\nPXDnUHBtA+5qbY60PiRJY3TSQKmqPwQOjxRfD2xvy9uBG4bKd1TVe1X1CjAFrE+yCrioqp6qqgIe\nGGkz3dejwHVt9rIR2F1Vh6vqCLAb2NTWXdvqjm5fkjQmp3tSfmVVHWjLbwAr2/Jq4LWheq+3stVt\nebT8mDZVdRR4B1g+Q1/Lgbdb3dG+JEljMuervNqMozqM5YxIckuSySSTb7311riHI0lnrdMNlDfb\nYSza88FWvh+4fKjeZa1sf1seLT+mTZIlwMXAoRn6OgRc0uqO9vUTqureqpqoqokVK1ac4j9TkjRb\npxsoO4Hpq662AI8NlW9uV26tZXDy/Zl2eOzdJNe0cyA3j7SZ7utG4Ik269kFbEiytJ2M3wDsauue\nbHVHty9JGpOT3oeS5HeAXwIuTfI6gyuvvgY8kmQr8CpwE0BVvZjkEeAl4ChwW1W937q6lcEVYxcC\nj7cHwH3Ag0mmGJz839z6OpzkK8Czrd6Xq2r64oAvAjuSfBV4rvUhSRqjDH7hPzdMTEzU5OTkuIeh\nM8UbG+fus+fO54FmL8neqpo4WT2/ekWS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC78eyhanLxE\nWFpwDBRJPzYa1N6XolPgIS9JUhcGiiSpCwNFktSFgSJJ6sJAkSR1YaBIkrrwsmEtfN5zIi0KzlAk\nSV04Q9HC4mxEWrScoUiSujBQJEldeMhL4+UhLumsYaBo/hge0lnNQ16SpC6coejMcDYinXMMFM2d\n4XH2Ot7P1r+RohMwUDQzw0LSLBko+jHDQ9IcGChnI4NB0hh4lZckqQtnKIudsxFJC4SBspgYHloI\nvPJLJ+AhL0lSF4t6hpJkE/AfgPOA366qr415SKfP2YekRW7RBkqS84D/BPwy8DrwbJKdVfXSeEc2\nC4aHzjYeBhOLOFCA9cBUVf0xQJIdwPXAwg8U6VwwGjIGzFlvMQfKauC1odevA1ePaSw/5uxDOj5n\nMWe9xRwos5LkFuCW9vLPkrx8ml1dCvxJn1F15bhOjeM6NWd2XL922r+AnZv76/TNdVx/eTaVFnOg\n7AcuH3p9WSs7RlXdC9w7140lmayqibn205vjOjWO69Q4rlNzro9rMV82/CywLsnaJOcDm4GdYx6T\nJJ2zFu0MpaqOJvknwC4Glw1/o6peHPOwJOmctWgDBaCq/gD4g3na3JwPm50hjuvUOK5T47hOzTk9\nrlR5lYUkae4W8zkUSdICYqCcRJJNSV5OMpXk9nnY3uVJnkzyUpIXk3y+lf9Gkv1Jnm+PTw+1uaON\n7+UkG4fKr0ryQlt3d5I53SST5Eetv+eTTLayZUl2J9nXnpfO57iS/NWhffJ8kneTfGEc+yvJN5Ic\nTPK9obJu+yfJBUkebuVPJ1kzh3H92yQ/SPJHSb6V5JJWvibJ/x3ab791psY1w9i6/ew677OHh8b0\noyTPz+c+y4k/G8b+HvtAVfk4wYPByf4fAh8Fzge+C1x5hre5Cvi5tvwR4H8DVwK/AfyL49S/so3r\nAmBtG+95bd0zwDVAgMeBT81xbD8CLh0p+zfA7W35dmDbfI9r5Of1BoNr5ud9fwG/CPwc8L0zsX+A\nW4HfasubgYfnMK4NwJK2vG1oXGuG643003VcM4yt28+u5z4bWf/vgH81n/uME382jP09Nv1whjKz\nD77epar+HzD99S5nTFUdqKrvtOU/Bb7P4FsBTuR6YEdVvVdVrwBTwPokq4CLquqpGrw7HgBuOAND\nvh7Y3pa3D21jHOO6DvhhVb16kvGekXFV1R8Ch4+zvV77Z7ivR4HrZjOLOt64quq/V9XR9vIpBvdx\nndCZGNeJxjaDse6zaa39TcDvzNRH73HN8Nkw9vfYNANlZsf7epeZPty7atPNnwWebkX/tB2i+MbQ\ntPZEY1zdlkfL56KAbyfZm8E3EACsrKoDbfkNYOUYxjVtM8f+Jx/3/oK+++eDNi0M3gGWdxjjP2Lw\nW+q0te3Qzf9M8gtD257PcfX62Z2Jsf0C8GZV7Rsqm9d9NvLZsGDeYwbKApXkZ4DfBb5QVe8C9zA4\n9PYJ4ACDKfd8+2RVfQL4FHBbkl8cXtl+2xnLZYMZ3Nz6GeC/tqKFsL+OMc79cyJJvgQcBb7Zig4A\nV7Sf8z8DHkpy0TwPa8H97Eb8PY79xWVe99lxPhs+MO73mIEys1l9vUtvST7E4A3zzar6PYCqerOq\n3q+qvwD+M4PDcTONcT/HHsaY89iran97Pgh8q43hzTaFnp7iH5zvcTWfAr5TVW+2MY59fzU9988H\nbZIsAS4GDp3uwJL8Q+BXgV9rH0S0wyOH2vJeBsfd/8p8jqvzz673PlsC/B3g4aHxzts+O95nAwvo\nPWagzGzev96lHa+8D/h+Vf3mUPmqoWp/G5i++mQnsLldnbEWWAc806bA7ya5pvV5M/DYHMb14SQf\nmV5mcFL3e237W1q1LUPbmJdxDTnmt8Zx768hPffPcF83Ak9MB8GpyuCP0/068Jmq+j9D5Ssy+FtD\nJPloG9cfz9e42nZ7/uy6jg34m8APquqDQ0bztc9O9NnAQnqPncoZ/HPxAXyawdUUPwS+NA/b+ySD\nKesfAc+3x6eBB4EXWvlOYNVQmy+18b3M0JVJwASD/4w/BP4j7UbW0xzXRxlcMfJd4MXpfcHg+Ooe\nYB/wbWDZfI6r9fdhBr9FXTxUNu/7i0GgHQD+nMFx6a099w/w0wwO6U0xuErno3MY1xSDY+XT77Hp\nK3v+bvv5Pg98B/hbZ2pcM4yt28+u5z5r5fcD/3ik7rzsM0782TD299j0wzvlJUldeMhLktSFgSJJ\n6sJAkSR1YaBIkrowUCRJXRgokqQuDBRJUhcGiiSpi/8Ppn4Z3eaRiCgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11351c278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "masked3d=cmask(cd,tu12000_cd[0]) # masking 4d timeseries with single 3d timepoint\n",
    "plt.hist(masked3d[np.nonzero(masked3d)],64,[0,20000],color='Orange');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 463,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAD8CAYAAABQFVIjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFi5JREFUeJzt3X+s3fV93/HnqzhQlAawjWdZBmZH8SaRSUvLlUFqWnWw\n2k7XBbYh5GkbnmYFTbAp0TY1sEhzSPJHvGllQ1uoWIkwLK3xaCOsqsxzDFv/GT+uAymBhPmmBIFl\nMLUNtNrEavreH+dz2fEJvr72/fgeX9/nQzo63/P5fj6f78ffc+55nc/3+z3HqSokSZqrnxr3ACRJ\n5wcDRZLUhYEiSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBookqYsl4x7AfLr88strzZo14x6G\nJC0o+/fv/+OqWnGqeosqUNasWcPk5OS4hyFJC0qSV2dTz0NekqQuDBRJUhcGiiSpCwNFktSFgSJJ\n6sJAkSR1YaBIkrowUCRJXRgokqQuFtU35aWecnd+oqy21RhGIp0bnKFIkrowUCRJXRgokqQuZhUo\nSX6c5IUkzyeZbGXLkuxNcqDdLx2qf1eSqSQvJ9k4VH5N62cqyb1J0sovSvJIK386yZqhNlvaNg4k\n2TJUvrbVnWptL5z77pAknanTmaH8tar6VFVNtMd3Avuqah2wrz0mydXAZuCTwCbgG0kuaG3uAz4H\nrGu3Ta18K3Csqj4B3ANsb30tA7YB1wLrgW1DwbUduKe1Odb6kM6K3J2fuEk60VwOed0I7GjLO4Cb\nhsp3VtV7VfUKMAWsT7IKuKSqnqqqAh4aaTPd16PADW32shHYW1VHq+oYsBfY1NZd3+qObl+SNAaz\nDZQCvpNkf5LbWtnKqjrUlt8AVrbl1cBrQ21fb2Wr2/Jo+Qltquo48A6wfIa+lgNvt7qjfUmSxmC2\n30P5dFUdTPIXgL1Jfji8sqoqyTl5AX4LwNsArrrqqjGPRpLOX7MKlKo62O4PJ/k2g/MZbyZZVVWH\n2uGsw636QeDKoeZXtLKDbXm0fLjN60mWAJcCR1r5L420+e9t3WVJlrRZynBfo2O/H7gfYGJi4pwM\nPZ0//LKjFrNTHvJK8tEkH5teBjYA3wd2A9NXXW0BHmvLu4HN7cqttQxOvj/TDo+9m+S6dg7k1pE2\n033dDDzRzrPsATYkWdpOxm8A9rR1T7a6o9uXJI3BbGYoK4Fvtyt8lwC/VVX/NcmzwK4kW4FXgVsA\nqurFJLuAl4DjwB1V9X7r63bgQeBi4PF2A3gAeDjJFHCUwVViVNXRJF8Fnm31vlJVR9vyF4GdSb4G\nPNf6kCSNSQYf9heHiYmJmpycHPcwtADN5TJhD3lpoUuyf+grIyflj0NKH8LvmUinz59ekSR1YaBI\nkrowUCRJXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sIvNkpn2eiXJM/2N+f9gUqNizMUSVIXzlCk\nBcyfiNG5xECRFhADROcyA0WaZ57j0PnKQJEWAUNM88GT8pKkLgwUSVIXBookqQsDRZLUhSflteh5\nKa7Uh4EinaMMOi00HvKSJHXhDEU6Bzgb0fnAGYokqQsDRZLUhYEiSerCQJEkdWGgSJK6MFAkSV0Y\nKJKkLgwUSVIXsw6UJBckeS7J77XHy5LsTXKg3S8dqntXkqkkLyfZOFR+TZIX2rp7k6SVX5TkkVb+\ndJI1Q222tG0cSLJlqHxtqzvV2l44t10hSZqL05mhfB74wdDjO4F9VbUO2Ncek+RqYDPwSWAT8I0k\nF7Q29wGfA9a126ZWvhU4VlWfAO4Btre+lgHbgGuB9cC2oeDaDtzT2hxrfUiSxmRWgZLkCuBvAL85\nVHwjsKMt7wBuGirfWVXvVdUrwBSwPskq4JKqeqqqCnhopM10X48CN7TZy0Zgb1UdrapjwF5gU1t3\nfas7un1J0hjMdoby74BfA/58qGxlVR1qy28AK9vyauC1oXqvt7LVbXm0/IQ2VXUceAdYPkNfy4G3\nW93RviRJY3DKQEnyq8Dhqtp/sjptxlE9B9ZLktuSTCaZfOutt8Y9HEk6b81mhvLzwGeT/BjYCVyf\n5D8Db7bDWLT7w63+QeDKofZXtLKDbXm0/IQ2SZYAlwJHZujrCHBZqzva1wmq6v6qmqiqiRUrVszi\nnytJOhOnDJSququqrqiqNQxOtj9RVX8f2A1MX3W1BXisLe8GNrcrt9YyOPn+TDs89m6S69o5kFtH\n2kz3dXPbRgF7gA1JlraT8RuAPW3dk63u6PYlSWMwl/8P5evAriRbgVeBWwCq6sUku4CXgOPAHVX1\nfmtzO/AgcDHweLsBPAA8nGQKOMoguKiqo0m+Cjzb6n2lqo625S8CO5N8DXiu9SFJGpMMPuwvDhMT\nEzU5OTnuYegcs1j/c6vatnj+9jU3SfZX1cSp6vlNeUlSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiS\nujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaK\nJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1MWScQ9Amk+5O+Me\ngnTecoYiSerCQJEkdWGgSJK6OGWgJPnpJM8k+V6SF5Pc3cqXJdmb5EC7XzrU5q4kU0leTrJxqPya\nJC+0dfcmSSu/KMkjrfzpJGuG2mxp2ziQZMtQ+dpWd6q1vbDPLpEknYnZzFDeA66vqr8KfArYlOQ6\n4E5gX1WtA/a1xyS5GtgMfBLYBHwjyQWtr/uAzwHr2m1TK98KHKuqTwD3ANtbX8uAbcC1wHpg21Bw\nbQfuaW2OtT4kSWNyykCpgT9tDz/SbgXcCOxo5TuAm9ryjcDOqnqvql4BpoD1SVYBl1TVU1VVwEMj\nbab7ehS4oc1eNgJ7q+poVR0D9jIItADXt7qj25ckjcGszqEkuSDJ88BhBm/wTwMrq+pQq/IGsLIt\nrwZeG2r+eitb3ZZHy09oU1XHgXeA5TP0tRx4u9Ud7Wt07LclmUwy+dZbb83mnytJOgOzCpSqer+q\nPgVcwWC28VdG1heDWcs5p6rur6qJqppYsWLFuIcjSeet07rKq6reBp5kcO7jzXYYi3Z/uFU7CFw5\n1OyKVnawLY+Wn9AmyRLgUuDIDH0dAS5rdUf7kiSNwWyu8lqR5LK2fDHwy8APgd3A9FVXW4DH2vJu\nYHO7cmstg5Pvz7TDY+8mua6dA7l1pM10XzcDT7RZzx5gQ5Kl7WT8BmBPW/dkqzu6fUnSGMzmp1dW\nATvalVo/Beyqqt9L8j+BXUm2Aq8CtwBU1YtJdgEvAceBO6rq/dbX7cCDwMXA4+0G8ADwcJIp4CiD\nq8SoqqNJvgo82+p9paqOtuUvAjuTfA14rvUhSRqTDD7sLw4TExM1OTk57mFojPwtr/+vti2ev33N\nTZL9VTVxqnp+U16S1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknq\nwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiS\npC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSF6cMlCRXJnkyyUtJ\nXkzy+Va+LMneJAfa/dKhNnclmUrycpKNQ+XXJHmhrbs3SVr5RUkeaeVPJ1kz1GZL28aBJFuGyte2\nulOt7YV9dokk6UzMZoZyHPjnVXU1cB1wR5KrgTuBfVW1DtjXHtPWbQY+CWwCvpHkgtbXfcDngHXt\ntqmVbwWOVdUngHuA7a2vZcA24FpgPbBtKLi2A/e0NsdaH5KkMTlloFTVoar6blv+E+AHwGrgRmBH\nq7YDuKkt3wjsrKr3quoVYApYn2QVcElVPVVVBTw00ma6r0eBG9rsZSOwt6qOVtUxYC+wqa27vtUd\n3b4kaQyWnE7ldijqZ4GngZVVdaitegNY2ZZXA08NNXu9lf1ZWx4tn27zGkBVHU/yDrB8uHykzXLg\n7ao6/iF9jY75NuA2gKuuumrW/1YtfLk74x6CtKjM+qR8kp8Bfgf4QlW9O7yuzTiq89i6qKr7q2qi\nqiZWrFgx7uFI0nlrVoGS5CMMwuRbVfW7rfjNdhiLdn+4lR8ErhxqfkUrO9iWR8tPaJNkCXApcGSG\nvo4Al7W6o31JksZgNld5BXgA+EFV/frQqt3A9FVXW4DHhso3tyu31jI4+f5MOzz2bpLrWp+3jrSZ\n7utm4Ik269kDbEiytJ2M3wDsaeuebHVHty9JGoPZnEP5eeAfAC8keb6V/Uvg68CuJFuBV4FbAKrq\nxSS7gJcYXCF2R1W939rdDjwIXAw83m4wCKyHk0wBRxlcJUZVHU3yVeDZVu8rVXW0LX8R2Jnka8Bz\nrQ9J0phk8GF/cZiYmKjJyclxD0PzxJPyM6tti+dvX3OTZH9VTZyqnt+UlyR1YaBIkrowUCRJXRgo\nkqQuTuub8pLOHx920YIn6jUXzlAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdWGgSJK6\nMFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBook\nqQsDRZLUhYEiSerCQJEkdWGgSJK6WDLuAUi95O6MewjSonbKGUqSbyY5nOT7Q2XLkuxNcqDdLx1a\nd1eSqSQvJ9k4VH5NkhfaunuTpJVflOSRVv50kjVDbba0bRxIsmWofG2rO9XaXjj3XSFJmovZHPJ6\nENg0UnYnsK+q1gH72mOSXA1sBj7Z2nwjyQWtzX3A54B17Tbd51bgWFV9ArgH2N76WgZsA64F1gPb\nhoJrO3BPa3Os9SFJGqNTBkpV/QFwdKT4RmBHW94B3DRUvrOq3quqV4ApYH2SVcAlVfVUVRXw0Eib\n6b4eBW5os5eNwN6qOlpVx4C9wKa27vpWd3T7kqQxOdOT8iur6lBbfgNY2ZZXA68N1Xu9la1uy6Pl\nJ7SpquPAO8DyGfpaDrzd6o72JUkakzlf5dVmHNVhLGdFktuSTCaZfOutt8Y9HEk6b51poLzZDmPR\n7g+38oPAlUP1rmhlB9vyaPkJbZIsAS4FjszQ1xHgslZ3tK+fUFX3V9VEVU2sWLHiNP+ZkqTZOtNA\n2Q1MX3W1BXhsqHxzu3JrLYOT78+0w2PvJrmunQO5daTNdF83A0+0Wc8eYEOSpe1k/AZgT1v3ZKs7\nun1J0pic8nsoSX4b+CXg8iSvM7jy6uvAriRbgVeBWwCq6sUku4CXgOPAHVX1fuvqdgZXjF0MPN5u\nAA8ADyeZYnDyf3Pr62iSrwLPtnpfqarpiwO+COxM8jXgudaHJGmMMvjAvzhMTEzU5OTkuIehs8Qv\nNs5dbVs87weavST7q2riVPX86RVJUhcGiiSpCwNFktSFgSJJ6sJAkSR1YaBIkrrw/0PRguQlwtK5\nxxmKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC78HorOeX7nRFoYnKFIkrpwhqJzirMR\naeFyhiJJ6sJAkSR14SEvjZWHuKTzh4GieWN4SOc3D3lJkrpwhqKzwtmItPgYKJozw0MSGCg6BcNC\n0mwZKPqA4SFpLgyU85DBIGkcvMpLktSFM5QFztmIpHOFgbKAGB6SzmUe8pIkdbGgZyhJNgH/HrgA\n+M2q+vqYh3TGnH1IWugWbKAkuQD4j8AvA68DzybZXVUvjXdkp2Z4SDofLeRDXuuBqar6o6r6v8BO\n4MYxj0mSFq0FO0MBVgOvDT1+Hbh2TGP5gLMPSYvVQg6UWUlyG3Bbe/inSV4+w64uB/64z6i6clyn\nx3HNIF/+iQ9E58S4PoTjOj1zHddfnE2lhRwoB4Erhx5f0cpOUFX3A/fPdWNJJqtqYq799Oa4To/j\nOj2O6/Qs9nEt5HMozwLrkqxNciGwGdg95jFJ0qK1YGcoVXU8yT8B9jC4bPibVfXimIclSYvWgg0U\ngKr6feD352lzcz5sdpY4rtPjuE6P4zo9i3pcqar52I4k6Ty3kM+hSJLOIQbKKSTZlOTlJFNJ7pyH\n7V2Z5MkkLyV5McnnW/mXkxxM8ny7/cpQm7va+F5OsnGo/JokL7R19yaZ05dkkvy49fd8kslWtizJ\n3iQH2v3S+RxXkr88tE+eT/Juki+MY38l+WaSw0m+P1TWbf8kuSjJI6386SRr5jCuf5Pkh0n+MMm3\nk1zWytck+T9D++03zta4Zhhbt+eu8z57ZGhMP07y/Hzus5z8vWHsr7EPVJW3k9wYnOz/EfBx4ELg\ne8DVZ3mbq4Cfa8sfA/4XcDXwZeBffEj9q9u4LgLWtvFe0NY9A1wHBHgc+Mwcx/Zj4PKRsn8N3NmW\n7wS2z/e4Rp6vNxhcMz/v+wv4ReDngO+fjf0D3A78RlveDDwyh3FtAJa05e1D41ozXG+kn67jmmFs\n3Z67nvtsZP2/Bf7VfO4zTv7eMPbX2PTNGcrM5v3nXarqUFV9ty3/CfADBr8KcDI3Ajur6r2qegWY\nAtYnWQVcUlVP1eDV8RBw01kY8o3Ajra8Y2gb4xjXDcCPqurVU4z3rIyrqv4AOPoh2+u1f4b7ehS4\nYTazqA8bV1X9t6o63h4+xeB7XCd1NsZ1srHNYKz7bFprfwvw2zP10XtcM7w3jP01Ns1AmdmH/bzL\nTG/uXbXp5s8CT7eif9oOUXxzaFp7sjGubsuj5XNRwHeS7M/gFwgAVlbVobb8BrByDOOatpkT/8jH\nvb+g7/75oE0Lg3eA5R3G+I8YfEqdtrYduvkfSX5haNvzOa5ez93ZGNsvAG9W1YGhsnndZyPvDefM\na8xAOUcl+Rngd4AvVNW7wH0MDr19CjjEYMo93z5dVZ8CPgPckeQXh1e2TztjuWwwgy+3fhb4L63o\nXNhfJxjn/jmZJF8CjgPfakWHgKva8/zPgN9Kcsk8D+uce+5G/F1O/OAyr/vsQ94bPjDu15iBMrNZ\n/bxLb0k+wuAF862q+l2Aqnqzqt6vqj8H/hODw3EzjfEgJx7GmPPYq+pguz8MfLuN4c02hZ6e4h+e\n73E1nwG+W1VvtjGOfX81PffPB22SLAEuBY6c6cCS/EPgV4G/196IaIdHjrTl/QyOu/+l+RxX5+eu\n9z5bAvxt4JGh8c7bPvuw9wbOodeYgTKzef95l3a88gHgB1X160Plq4aq/S1g+uqT3cDmdnXGWmAd\n8EybAr+b5LrW563AY3MY10eTfGx6mcFJ3e+37W9p1bYMbWNexjXkhE+N495fQ3run+G+bgaemA6C\n05XBf073a8Bnq+p/D5WvyOD/GiLJx9u4/mi+xtW22/O56zo24K8DP6yqDw4Zzdc+O9l7A+fSa+x0\nzuAvxhvwKwyupvgR8KV52N6nGUxZ/xB4vt1+BXgYeKGV7wZWDbX5UhvfywxdmQRMMPhj/BHwH2hf\nZD3DcX2cwRUj3wNenN4XDI6v7gMOAN8Bls3nuFp/H2XwKerSobJ5318MAu0Q8GcMjktv7bl/gJ9m\ncEhvisFVOh+fw7imGBwrn36NTV/Z83fa8/s88F3gb56tcc0wtm7PXc991sofBP7xSN152Wec/L1h\n7K+x6ZvflJckdeEhL0lSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC7+H2KgFlWB\nwU8hAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x117e117b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "masked4d=cmask(cd,tu12000_cd)    # masking 4d timeseries with 4d timeseries\n",
    "plt.hist(masked4d[np.nonzero(masked4d)],64,[0,20000],color='Green');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 308,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAD8CAYAAABQFVIjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFchJREFUeJzt3X+s3fV93/HnqzhQlAawjWdZhsyOcCfBpNFyZZCaVh1s\ntpN1ga0IudqGp1lBE2xKtE0FFmkOSf6IN61MaAsVGxGGpTUebYRVlXmOYes/48d1ICWQMN+UILAM\ndm0DrTaxmb73x/nc6PgEX1/7fnyP7ft8SEfnez7fz+fzfft7f7zu98c5TlUhSdJc/cy4C5AknR8M\nFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4WjbuA+XT55ZfXqlWrxl2GJJ1T\n9u7d+6dVtexk/RZUoKxatYrJyclxlyFJ55Qkb8ymn6e8JEldGCiSpC4MFElSFwaKJKkLA0WS1IWB\nIknqwkCRJHVhoEiSujBQJEldLKh3yks95b78VFttqTFUIp0dPEKRJHVhoEiSujBQJEldzCpQkvw4\nyctJXkoy2dqWJNmdZF97XjzU/94kU0leS7J+qP26Ns9UkgeSpLVflOTx1v5cklVDYza1bexLsmmo\nfXXrO9XGXjj33SFJOl2ncoTy16vq2qqaaK/vAfZU1RpgT3tNkquBjcA1wAbgG0kuaGMeBD4PrGmP\nDa19M3C0qq4C7ge2trmWAFuA64G1wJah4NoK3N/GHG1zSGdE7stPPSQdby6nvG4GtrXlbcAtQ+3b\nq+qDqnodmALWJlkBXFJVz1ZVAY+OjJme6wngpnb0sh7YXVVHquoosBvY0Nbd2PqObl+SNAazDZQC\nvpNkb5I7WtvyqjrQlt8GlrfllcCbQ2Pfam0r2/Jo+3FjquoY8B6wdIa5lgLvtr6jc0mSxmC270P5\ndFXtT/KXgN1Jfji8sqoqyVl5A34LwDsAPvnJT465Gkk6f80qUKpqf3s+mOTbDK5nvJNkRVUdaKez\nDrbu+4Erh4Zf0dr2t+XR9uExbyVZBFwKHG7tvzoy5r+3dZclWdSOUobnGq39IeAhgImJibMy9HT+\n8M2OWshOesoryceTfGJ6GVgHfB/YCUzfdbUJeLIt7wQ2tju3VjO4+P58Oz32fpIb2jWQ20fGTM91\nK/B0u86yC1iXZHG7GL8O2NXWPdP6jm5fkjQGszlCWQ58u93huwj4nar6r0leAHYk2Qy8AdwGUFWv\nJNkBvAocA+6qqg/bXHcCjwAXA0+1B8DDwGNJpoAjDO4So6qOJPkq8ELr95WqOtKW7wa2J/ka8GKb\nQ5I0Jhn8sb8wTExM1OTk5LjL0DloLrcJe8pL57oke4feMnJCfjik9BF8n4l06vzoFUlSFwaKJKkL\nA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC58Y6N0ho2+SfJMv3PeD6jUuHiEIknqwiMU6Rzm\nR8TobGKgSOcQA0RnMwNFmmde49D5ykCRFgBDTPPBi/KSpC4MFElSFwaKJKkLA0WS1IWBIknqwru8\npLOU7znRucYjFElSFx6hSGcBj0Z0PvAIRZLUhYEiSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIX\nBookqYtZB0qSC5K8mOQP2uslSXYn2deeFw/1vTfJVJLXkqwfar8uyctt3QNJ0tovSvJ4a38uyaqh\nMZvaNvYl2TTUvrr1nWpjL5zbrpAkzcWpHKF8AfjB0Ot7gD1VtQbY016T5GpgI3ANsAH4RpIL2pgH\ngc8Da9pjQ2vfDBytqquA+4Gtba4lwBbgemAtsGUouLYC97cxR9sckqQxmVWgJLkC+FvAfxpqvhnY\n1pa3AbcMtW+vqg+q6nVgClibZAVwSVU9W1UFPDoyZnquJ4Cb2tHLemB3VR2pqqPAbmBDW3dj6zu6\nfUnSGMz2COXfAb8J/MVQ2/KqOtCW3waWt+WVwJtD/d5qbSvb8mj7cWOq6hjwHrB0hrmWAu+2vqNz\nSZLG4KSBkuTXgINVtfdEfdoRR/UsrJckdySZTDJ56NChcZcjSeet2Ryh/BLwuSQ/BrYDNyb5z8A7\n7TQW7flg678fuHJo/BWtbX9bHm0/bkySRcClwOEZ5joMXNb6js51nKp6qKomqmpi2bJls/jnSpJO\nx0kDparuraorqmoVg4vtT1fV3wd2AtN3XW0CnmzLO4GN7c6t1Qwuvj/fTo+9n+SGdg3k9pEx03Pd\n2rZRwC5gXZLF7WL8OmBXW/dM6zu6fUnSGMzl/0P5OrAjyWbgDeA2gKp6JckO4FXgGHBXVX3YxtwJ\nPAJcDDzVHgAPA48lmQKOMAguqupIkq8CL7R+X6mqI235bmB7kq8BL7Y5JEljksEf+wvDxMRETU5O\njrsMnQMWwn94VVsWzs++5ibJ3qqaOFk/3ykvSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBook\nqQsDRZLUhYEiSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdWGg\nSJK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerCQJEkdWGgSJK6MFAkSV0YKJKkLgwUSVIX\nBookqYuTBkqSn03yfJLvJXklyX2tfUmS3Un2tefFQ2PuTTKV5LUk64far0vyclv3QJK09ouSPN7a\nn0uyamjMpraNfUk2DbWvbn2n2tgL++wSSdLpmM0RygfAjVX114BrgQ1JbgDuAfZU1RpgT3tNkquB\njcA1wAbgG0kuaHM9CHweWNMeG1r7ZuBoVV0F3A9sbXMtAbYA1wNrgS1DwbUVuL+NOdrmkCSNyUkD\npQb+vL38WHsUcDOwrbVvA25pyzcD26vqg6p6HZgC1iZZAVxSVc9WVQGPjoyZnusJ4KZ29LIe2F1V\nR6rqKLCbQaAFuLH1Hd2+JGkMZnUNJckFSV4CDjL4Bf8csLyqDrQubwPL2/JK4M2h4W+1tpVtebT9\nuDFVdQx4D1g6w1xLgXdb39G5Rmu/I8lkkslDhw7N5p8rSToNswqUqvqwqq4FrmBwtPFXR9YXg6OW\ns05VPVRVE1U1sWzZsnGXI0nnrVO6y6uq3gWeYXDt4512Gov2fLB12w9cOTTsita2vy2Pth83Jski\n4FLg8AxzHQYua31H55IkjcFs7vJaluSytnwx8DeBHwI7gem7rjYBT7blncDGdufWagYX359vp8fe\nT3JDuwZy+8iY6bluBZ5uRz27gHVJFreL8euAXW3dM63v6PYlSWOw6ORdWAFsa3dq/Qywo6r+IMn/\nBHYk2Qy8AdwGUFWvJNkBvAocA+6qqg/bXHcCjwAXA0+1B8DDwGNJpoAjDO4So6qOJPkq8ELr95Wq\nOtKW7wa2J/ka8GKbQ5I0Jhn8sb8wTExM1OTk5LjL0Dkg92XcJZxxtWXh/OxrbpLsraqJk/XznfKS\npC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWB\nIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEld\nGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujhpoCS5MskzSV5N8kqSL7T2JUl2J9nX\nnhcPjbk3yVSS15KsH2q/LsnLbd0DSdLaL0ryeGt/LsmqoTGb2jb2Jdk01L669Z1qYy/ss0skSadj\nNkcox4B/XlVXAzcAdyW5GrgH2FNVa4A97TVt3UbgGmAD8I0kF7S5HgQ+D6xpjw2tfTNwtKquAu4H\ntra5lgBbgOuBtcCWoeDaCtzfxhxtc0iSxuSkgVJVB6rqu235z4AfACuBm4Ftrds24Ja2fDOwvao+\nqKrXgSlgbZIVwCVV9WxVFfDoyJjpuZ4AbmpHL+uB3VV1pKqOAruBDW3dja3v6PYlSWNwStdQ2qmo\nXwCeA5ZX1YG26m1geVteCbw5NOyt1rayLY+2Hzemqo4B7wFLZ5hrKfBu6zs612jNdySZTDJ56NCh\nU/jXSpJOxawDJcnPAb8HfLGq3h9e1444qnNtXVTVQ1U1UVUTy5YtG3c5knTemlWgJPkYgzD5VlX9\nfmt+p53Goj0fbO37gSuHhl/R2va35dH248YkWQRcChyeYa7DwGWt7+hckqQxmM1dXgEeBn5QVb81\ntGonMH3X1SbgyaH2je3OrdUMLr4/306PvZ/khjbn7SNjpue6FXi6HfXsAtYlWdwuxq8DdrV1z7S+\no9uXJI3BopN34ZeAfwC8nOSl1vYvga8DO5JsBt4AbgOoqleS7ABeZXCH2F1V9WEbdyfwCHAx8FR7\nwCCwHksyBRxhcJcYVXUkyVeBF1q/r1TVkbZ8N7A9ydeAF9sckqQxyeCP/YVhYmKiJicnx12GzgG5\nL+Mu4YyrLQvnZ19zk2RvVU2crN9sjlCk89pCCA9pPvjRK5KkLgwUSVIXBookqQuvoUgL1EddO/JC\nvebCIxRJUhcGiiSpCwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJXRgokqQuDBRJUhcGiiSpCwNFktSF\ngSJJ6sJAkSR1YaBIkrowUCRJXRgokqQuDBRJUhcGiiSpCwNFktSFgSJJ6sJAkSR1YaBIkrowUCRJ\nXRgokqQuDBRJUhcnDZQk30xyMMn3h9qWJNmdZF97Xjy07t4kU0leS7J+qP26JC+3dQ8kSWu/KMnj\nrf25JKuGxmxq29iXZNNQ++rWd6qNvXDuu0KSNBezOUJ5BNgw0nYPsKeq1gB72muSXA1sBK5pY76R\n5II25kHg88Ca9pieczNwtKquAu4Htra5lgBbgOuBtcCWoeDaCtzfxhxtc0iSxuikgVJVfwQcGWm+\nGdjWlrcBtwy1b6+qD6rqdWAKWJtkBXBJVT1bVQU8OjJmeq4ngJva0ct6YHdVHamqo8BuYENbd2Pr\nO7p9SdKYnO41lOVVdaAtvw0sb8srgTeH+r3V2la25dH248ZU1THgPWDpDHMtBd5tfUfnkiSNyZwv\nyrcjjupQyxmR5I4kk0kmDx06NO5yJOm8dbqB8k47jUV7Ptja9wNXDvW7orXtb8uj7ceNSbIIuBQ4\nPMNch4HLWt/RuX5KVT1UVRNVNbFs2bJT/GdKkmbrdANlJzB919Um4Mmh9o3tzq3VDC6+P99Oj72f\n5IZ2DeT2kTHTc90KPN2OenYB65Isbhfj1wG72rpnWt/R7UuSxmTRyTok+V3gV4HLk7zF4M6rrwM7\nkmwG3gBuA6iqV5LsAF4FjgF3VdWHbao7GdwxdjHwVHsAPAw8lmSKwcX/jW2uI0m+CrzQ+n2lqqZv\nDrgb2J7ka8CLbQ5J0hhl8Af/wjAxMVGTk5PjLkNnmdyXcZdw1qgtC+f3gWYvyd6qmjhZP98pL0nq\nwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiS\npC4MFElSFwaKJKkLA0WS1IWBIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKkLA0WS1IWB\nIknqwkCRJHVhoEiSujBQJEldGCiSpC4MFElSFwaKJKmLczpQkmxI8lqSqST3jLseSVrIztlASXIB\n8B+AzwBXA7+R5OrxViVJC9c5GyjAWmCqqv6kqv4vsB24ecw1SdKCdS4HykrgzaHXb7U2SdIYLBp3\nAWdakjuAO9rLP0/y2mlOdTnwp32q6sq6To11zSBfzmjTWVHXR7CuUzPXuv7ybDqdy4GyH7hy6PUV\nre04VfUQ8NBcN5Zksqom5jpPb9Z1aqzr1FjXqVnodZ3Lp7xeANYkWZ3kQmAjsHPMNUnSgnXOHqFU\n1bEk/wTYBVwAfLOqXhlzWZK0YJ2zgQJQVX8I/OE8bW7Op83OEOs6NdZ1aqzr1CzoulJV87EdSdJ5\n7ly+hiJJOosYKCcx3x/vkuTKJM8keTXJK0m+0Nq/nGR/kpfa47NDY+5t9b2WZP1Q+3VJXm7rHkjy\nU/eEnmJtP27zvZRksrUtSbI7yb72vHg+60ryV4b2yUtJ3k/yxXHsryTfTHIwyfeH2rrtnyQXJXm8\ntT+XZNUc6vo3SX6Y5I+TfDvJZa19VZL/M7TffvtM1TVDbd2+dp332eNDNf04yUvzuc9y4t8NY/8e\n+4mq8nGCB4OL/T8CPgVcCHwPuPoMb3MF8Itt+RPA/2Lw0TJfBv7FR/S/utV1EbC61XtBW/c8cAMQ\n4CngM3Os7cfA5SNt/xq4py3fA2yd77pGvl5vM7hnft73F/ArwC8C3z8T+we4E/jttrwReHwOda0D\nFrXlrUN1rRruNzJP17pmqK3b167nPhtZ/2+BfzWf+4wT/24Y+/fY9MMjlJnN+8e7VNWBqvpuW/4z\n4AfM/AkANwPbq+qDqnodmALWJlkBXFJVz9bgu+NR4JYzUPLNwLa2vG1oG+Oo6ybgR1X1xknqPSN1\nVdUfAUc+Ynu99s/wXE8AN83mKOqj6qqq/1ZVx9rLZxm8j+uEzkRdJ6ptBmPdZ9Pa+NuA351pjt51\nzfC7YezfY9MMlJmN9eNd2uHmLwDPtaZ/2k5RfHPosPZENa5sy6Ptc1HAd5LszeATCACWV9WBtvw2\nsHwMdU3byPE/5OPeX9B3//xkTAuD94ClHWr8Rwz+Sp22up26+R9Jfnlo2/NZV6+v3Zmo7ZeBd6pq\n31DbvO6zkd8NZ833mIFylkryc8DvAV+sqveBBxmcersWOMDgkHu+fbqqrmXwCc93JfmV4ZXtr52x\n3DaYwZtbPwf8l9Z0Nuyv44xz/5xIki8Bx4BvtaYDwCfb1/mfAb+T5JJ5Luus+9qN+A2O/8NlXvfZ\nR/xu+Ilxf48ZKDOb1ce79JbkYwy+Yb5VVb8PUFXvVNWHVfUXwH9kcDpuphr3c/xpjDnXXlX72/NB\n4NuthnfaIfT0If7B+a6r+Qzw3ap6p9U49v3V9Nw/PxmTZBFwKXD4dAtL8g+BXwP+XvtFRDs9crgt\n72Vw3v3n57Ouzl+73vtsEfB3gceH6p23ffZRvxs4i77HDJSZzfvHu7TzlQ8DP6iq3xpqXzHU7e8A\n03ef7AQ2trszVgNrgOfbIfD7SW5oc94OPDmHuj6e5BPTywwu6n6/bX9T67ZpaBvzUteQ4/5qHPf+\nGtJz/wzPdSvw9HQQnKokG4DfBD5XVf97qH1ZBv/XEEk+1er6k/mqq22359eua23A3wB+WFU/OWU0\nX/vsRL8bOJu+x07lCv5CfACfZXA3xY+AL83D9j7N4JD1j4GX2uOzwGPAy619J7BiaMyXWn2vMXRn\nEjDB4IfxR8C/p72R9TTr+hSDO0a+B7wyvS8YnF/dA+wDvgMsmc+62nwfZ/BX1KVDbfO+vxgE2gHg\n/zE4L7255/4BfpbBKb0pBnfpfGoOdU0xOFc+/T02fWfPr7ev70vAd4G/fabqmqG2bl+7nvustT8C\n/OORvvOyzzjx74axf49NP3ynvCSpC095SZK6MFAkSV0YKJKkLgwUSVIXBookqQsDRZLUhYEiSerC\nQJEkdfH/Aaj8B7ROi20GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a666278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tl8000_cd=clthresh(cd,8000)\n",
    "range4d=cmask(tl8000_cd,tu12000_cd)    # masking 4d timeseries with 4d timeseries\n",
    "plt.hist(range4d[np.nonzero(range4d)],64,[0,20000],color='Green');\n",
    "output_cifti=nhw+'/range_demo.dtseries.nii'\n",
    "csave(output_cifti,range4d,ch,cn,ch,cf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cinfo(infile):\n",
    "    '''\n",
    "    Returns information about input cifti header\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    infile: str, Cifti2Image\n",
    "        string containing path to cifti file or an already loaded cifti object\n",
    "    verbose: bool\n",
    "        report relatively esoteric cifti header info\n",
    "    '''\n",
    "    if isinstance(infile,nb.cifti2.cifti2.Cifti2Header):\n",
    "        ch=infile\n",
    "    elif isinstance(infile,nb.cifti2.cifti2.Cifti2Image):\n",
    "        cin=infile\n",
    "        cd,ch,cn,cx,cf=cread(infile)\n",
    "    elif isinstance(infile,str):\n",
    "        if not os.path.exists(infile):\n",
    "            raise ValueError('File does not exist, please provide a valid path to a cifti file.')\n",
    "        try:\n",
    "            cd,ch,cn,cx,cf=cread(infile)\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid cifti file.')\n",
    "    ch0=ch.get_index_map(0)\n",
    "    ch1=ch.get_index_map(1)\n",
    "    if ch.get_index_map(0).indices_map_to_data_type=='CIFTI_INDEX_TYPE_SERIES': # support for other cifti types pending\n",
    "        print('cifti data dype = fMRI timeseries')\n",
    "        me=ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.meter_exponent\n",
    "        if me==-3:\n",
    "            units='MM'\n",
    "        else:\n",
    "            units='x10^{e} M'.format(e=me)\n",
    "        print('Xdim = {x} {u} x {xx}'.format(\n",
    "            x=np.abs(ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix[0,0]),\n",
    "            u=units,\n",
    "            xx=ch1.volume.volume_dimensions[0]))\n",
    "        print('Ydim = {y} {u} x {yy}'.format(\n",
    "            y=np.abs(ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix[1,1]),\n",
    "            u=units,\n",
    "            yy=ch1.volume.volume_dimensions[1]))\n",
    "        print('Zdim = {z} {u} x {zz}'.format(\n",
    "            z=np.abs(ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix[2,2]),\n",
    "            u=units,\n",
    "            zz=ch1.volume.volume_dimensions[2]))\n",
    "        print('Tdim = {y} {z} x {t}'.format(\n",
    "            y=ch.get_index_map(0).series_step,\n",
    "            z=ch.get_index_map(0).series_unit,\n",
    "            t=ch.get_index_map(0).number_of_series_points))\n",
    "    disp=[]\n",
    "    voxels=0\n",
    "    print('INDEX      OFFSET     SIZE       STRUCTURE')\n",
    "    for idx, bm in enumerate(ch1.brain_models):\n",
    "        disp.append([idx,bm.index_offset,bm.index_count,bm.brain_structure])\n",
    "        if bm.model_type=='CIFTI_MODEL_TYPE_VOXELS':\n",
    "            voxels=voxels+bm.index_count\n",
    "    col_width = max(len(str(word)) for row in disp for word in row) + 2  # padding\n",
    "    for row in disp:\n",
    "        print(\"{: <10} {: <10} {: <10} {: <20}\".format(*row))\n",
    "    print('total left cortex vertices = {lv} of {lt}'.format(\n",
    "        lv=list(ch1.brain_models)[0].index_count,\n",
    "        lt=list(ch1.brain_models)[0].surface_number_of_vertices))\n",
    "    print('total right cortex vertices = {rv} of {rt}'.format(\n",
    "        rv=list(ch1.brain_models)[1].index_count,\n",
    "        rt=list(ch1.brain_models)[1].surface_number_of_vertices))\n",
    "    print('total subcortical voxels = {v}'.format(v=voxels))\n",
    "    print('total data points = {total}'.format(\n",
    "        total=list(ch1.brain_models)[0].index_count+list(ch1.brain_models)[1].index_count+voxels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": false
   },
   "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 dype = fMRI timeseries\n",
      "Xdim = 2.0 MM x 91\n",
      "Ydim = 2.0 MM x 109\n",
      "Zdim = 2.0 MM x 91\n",
      "Tdim = 0.72 SECOND x 1200\n",
      "INDEX      OFFSET     SIZE       STRUCTURE\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",
      "total left cortex vertices = 29696 of 32492\n",
      "total right cortex vertices = 29716 of 32492\n",
      "total subcortical voxels = 31870\n",
      "total data points = 91282\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "# cinfo examples:\n",
    "print(cinfo(path_ciftiTS))\n",
    "#print(cinfo(data_ciftiTS))\n",
    "#print(cinfo(ch))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cvoxels(header,MNI=True,trim=False):\n",
    "    '''\n",
    "    Returns array where each element is an array of coordinates for a particular cifti brain structure.\n",
    "    Coordinates can be voxel indices (i,j,k) or MNI coordinates (x,y,z; default).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    header: Cifti2Header\n",
    "        cifti header object containing the voxel index mapping and list of brain structures\n",
    "    MNI: bool\n",
    "        True returns MNI coordinates, False returns voxel indices, default=True\n",
    "    trim: bool\n",
    "        removes array elements corresponding to cortical surfaces, default=False\n",
    "    '''\n",
    "    if not isinstance(header,nb.cifti2.cifti2.Cifti2Header):\n",
    "        raise ValueError('Please provide a valid cifti header.')\n",
    "\n",
    "    ch1=header.get_index_map(1)\n",
    "    voxarray=[]\n",
    "    if MNI 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 MNI is False:\n",
    "                voxarray.append(np.asarray(bm.voxel_indices_ijk))\n",
    "            else:\n",
    "                voxarray.append(apply_affine(warp,bm.voxel_indices_ijk))\n",
    "        else:\n",
    "            if trim is False:\n",
    "                voxarray.append([[],[],[]])\n",
    "    return np.asarray(voxarray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 562,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[49 66 28] [ -8.   6. -16.]\n",
      "[50 66 28] [-10.   6. -16.]\n",
      "[48 67 28] [ -6.   8. -16.]\n",
      "[49 67 28] [ -8.   8. -16.]\n",
      "[50 67 28] [-10.   8. -16.]\n",
      "[48 65 29] [ -6.   4. -14.]\n",
      "[49 65 29] [ -8.   4. -14.]\n",
      "[50 65 29] [-10.   4. -14.]\n",
      "[48 66 29] [ -6.   6. -14.]\n",
      "[49 66 29] [ -8.   6. -14.]\n",
      "[50 66 29] [-10.   6. -14.]\n",
      "[48 67 29] [ -6.   8. -14.]\n",
      "[49 67 29] [ -8.   8. -14.]\n",
      "[50 67 29] [-10.   8. -14.]\n",
      "[50 64 30] [-10.   2. -12.]\n",
      "[47 65 30] [ -4.   4. -12.]\n",
      "[48 65 30] [ -6.   4. -12.]\n",
      "[49 65 30] [ -8.   4. -12.]\n",
      "[50 65 30] [-10.   4. -12.]\n",
      "[47 66 30] [ -4.   6. -12.]\n",
      "[48 66 30] [ -6.   6. -12.]\n",
      "[49 66 30] [ -8.   6. -12.]\n",
      "[50 66 30] [-10.   6. -12.]\n",
      "[47 67 30] [ -4.   8. -12.]\n",
      "[48 67 30] [ -6.   8. -12.]\n",
      "[49 67 30] [ -8.   8. -12.]\n",
      "[50 67 30] [-10.   8. -12.]\n",
      "[48 68 30] [ -6.  10. -12.]\n",
      "[49 68 30] [ -8.  10. -12.]\n",
      "[50 68 30] [-10.  10. -12.]\n",
      "[51 68 30] [-12.  10. -12.]\n",
      "[50 69 30] [-10.  12. -12.]\n",
      "[51 69 30] [-12.  12. -12.]\n",
      "[47 65 31] [ -4.   4. -10.]\n",
      "[48 65 31] [ -6.   4. -10.]\n",
      "[49 65 31] [ -8.   4. -10.]\n",
      "[50 65 31] [-10.   4. -10.]\n",
      "[47 66 31] [ -4.   6. -10.]\n",
      "[48 66 31] [ -6.   6. -10.]\n",
      "[49 66 31] [ -8.   6. -10.]\n",
      "[50 66 31] [-10.   6. -10.]\n",
      "[47 67 31] [ -4.   8. -10.]\n",
      "[48 67 31] [ -6.   8. -10.]\n",
      "[49 67 31] [ -8.   8. -10.]\n",
      "[50 67 31] [-10.   8. -10.]\n",
      "[47 68 31] [ -4.  10. -10.]\n",
      "[48 68 31] [ -6.  10. -10.]\n",
      "[49 68 31] [ -8.  10. -10.]\n",
      "[50 68 31] [-10.  10. -10.]\n",
      "[51 68 31] [-12.  10. -10.]\n",
      "[48 69 31] [ -6.  12. -10.]\n",
      "[49 69 31] [ -8.  12. -10.]\n",
      "[50 69 31] [-10.  12. -10.]\n",
      "[51 69 31] [-12.  12. -10.]\n",
      "[49 70 31] [ -8.  14. -10.]\n",
      "[50 70 31] [-10.  14. -10.]\n",
      "[51 70 31] [-12.  14. -10.]\n",
      "[52 70 31] [-14.  14. -10.]\n",
      "[50 71 31] [-10.  16. -10.]\n",
      "[51 71 31] [-12.  16. -10.]\n",
      "[52 71 31] [-14.  16. -10.]\n",
      "[47 65 32] [-4.  4. -8.]\n",
      "[48 65 32] [-6.  4. -8.]\n",
      "[49 65 32] [-8.  4. -8.]\n",
      "[47 66 32] [-4.  6. -8.]\n",
      "[48 66 32] [-6.  6. -8.]\n",
      "[49 66 32] [-8.  6. -8.]\n",
      "[50 66 32] [-10.   6.  -8.]\n",
      "[47 67 32] [-4.  8. -8.]\n",
      "[48 67 32] [-6.  8. -8.]\n",
      "[49 67 32] [-8.  8. -8.]\n",
      "[50 67 32] [-10.   8.  -8.]\n",
      "[47 68 32] [ -4.  10.  -8.]\n",
      "[48 68 32] [ -6.  10.  -8.]\n",
      "[49 68 32] [ -8.  10.  -8.]\n",
      "[50 68 32] [-10.  10.  -8.]\n",
      "[48 69 32] [ -6.  12.  -8.]\n",
      "[49 69 32] [ -8.  12.  -8.]\n",
      "[50 69 32] [-10.  12.  -8.]\n",
      "[51 69 32] [-12.  12.  -8.]\n",
      "[48 70 32] [ -6.  14.  -8.]\n",
      "[49 70 32] [ -8.  14.  -8.]\n",
      "[50 70 32] [-10.  14.  -8.]\n",
      "[51 70 32] [-12.  14.  -8.]\n",
      "[52 70 32] [-14.  14.  -8.]\n",
      "[49 71 32] [ -8.  16.  -8.]\n",
      "[50 71 32] [-10.  16.  -8.]\n",
      "[51 71 32] [-12.  16.  -8.]\n",
      "[52 71 32] [-14.  16.  -8.]\n",
      "[47 65 33] [-4.  4. -6.]\n",
      "[48 65 33] [-6.  4. -6.]\n",
      "[47 66 33] [-4.  6. -6.]\n",
      "[48 66 33] [-6.  6. -6.]\n",
      "[46 67 33] [-2.  8. -6.]\n",
      "[47 67 33] [-4.  8. -6.]\n",
      "[48 67 33] [-6.  8. -6.]\n",
      "[47 68 33] [ -4.  10.  -6.]\n",
      "[48 68 33] [ -6.  10.  -6.]\n",
      "[47 69 33] [ -4.  12.  -6.]\n",
      "[48 69 33] [ -6.  12.  -6.]\n",
      "[48 70 33] [ -6.  14.  -6.]\n",
      "[49 70 33] [ -8.  14.  -6.]\n",
      "[49 71 33] [ -8.  16.  -6.]\n",
      "[50 71 33] [-10.  16.  -6.]\n",
      "[49 72 33] [ -8.  18.  -6.]\n",
      "[50 72 33] [-10.  18.  -6.]\n",
      "[47 65 34] [-4.  4. -4.]\n",
      "[48 65 34] [-6.  4. -4.]\n",
      "[46 66 34] [-2.  6. -4.]\n",
      "[47 66 34] [-4.  6. -4.]\n",
      "[48 66 34] [-6.  6. -4.]\n",
      "[46 67 34] [-2.  8. -4.]\n",
      "[47 67 34] [-4.  8. -4.]\n",
      "[48 67 34] [-6.  8. -4.]\n",
      "[46 68 34] [ -2.  10.  -4.]\n",
      "[47 68 34] [ -4.  10.  -4.]\n",
      "[48 68 34] [ -6.  10.  -4.]\n",
      "[47 69 34] [ -4.  12.  -4.]\n",
      "[48 69 34] [ -6.  12.  -4.]\n",
      "[48 70 34] [ -6.  14.  -4.]\n",
      "[48 71 34] [ -6.  16.  -4.]\n",
      "[49 71 34] [ -8.  16.  -4.]\n",
      "[50 71 34] [-10.  16.  -4.]\n",
      "[49 72 34] [ -8.  18.  -4.]\n",
      "[50 72 34] [-10.  18.  -4.]\n",
      "[47 65 35] [-4.  4. -2.]\n",
      "[46 66 35] [-2.  6. -2.]\n",
      "[47 66 35] [-4.  6. -2.]\n",
      "[46 67 35] [-2.  8. -2.]\n",
      "[47 67 35] [-4.  8. -2.]\n",
      "[47 68 35] [ -4.  10.  -2.]\n",
      "[47 69 35] [ -4.  12.  -2.]\n",
      "[47 70 35] [ -4.  14.  -2.]\n",
      "[48 71 35] [ -6.  16.  -2.]\n",
      "[48 72 35] [ -6.  18.  -2.]\n"
     ]
    }
   ],
   "source": [
    "# cvoxels example:\n",
    "for i in range(0,len(cvoxels(ch,MNI=False,trim=True)[0])):\n",
    "    print(cvoxels(ch,MNI=False,trim=True)[0][i], cvoxels(ch,MNI=True,trim=True)[0][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def csplit(incifti,gii=False):\n",
    "    '''\n",
    "    Split cifti overlay into separate left surface, right surface, and subcortical voxel arrays.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    incifti: str, Cifti2Image\n",
    "        string containing path to cifti file or an already loaded cifti object\n",
    "    gii: bool\n",
    "        returns surfaces as GiftiImage objects, default=False\n",
    "    '''\n",
    "    if isinstance(incifti,nb.cifti2.cifti2.Cifti2Image):\n",
    "        cd=np.asarray(incifti.get_data()).copy()\n",
    "        ch=incifti.header.copy()\n",
    "    elif isinstance(incifti,str):\n",
    "        if not os.path.exists(incifti):\n",
    "            raise ValueError('File does not exist, please provide a valid path to a cifti file.')\n",
    "        try:\n",
    "            cd,ch,cn,cx,cf=cread(incifti)\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid cifti file.')   \n",
    "    mim = ch.matrix[1]\n",
    "    bm1 = list(mim.brain_models)[0]\n",
    "    bm2 = list(mim.brain_models)[1]\n",
    "    lidx = list(bm1.vertex_indices)\n",
    "    ridx = [bm1.surface_number_of_vertices + val for val in bm2.vertex_indices]\n",
    "    bidx = np.concatenate((lidx, ridx))\n",
    "    # split cifti overlay into left surface, right surface, and subcortical volume\n",
    "    clh=cd[:,0:len(lidx)]\n",
    "    crh=cd[:,len(lidx):len(lidx)+len(ridx)+1]\n",
    "    csc=cd[:,len(lidx)+len(ridx)+1:cd.shape[1]]\n",
    "    if gii is True:\n",
    "        glh=nb.gifti.gifti.GiftiImage()\n",
    "        grh=nb.gifti.gifti.GiftiImage()\n",
    "        glh.add_gifti_data_array(nb.gifti.gifti.GiftiDataArray(clh))\n",
    "        grh.add_gifti_data_array(nb.gifti.gifti.GiftiDataArray(crh))\n",
    "        left=glh\n",
    "        right=grh\n",
    "    else:\n",
    "        left=clh\n",
    "        right=crh\n",
    "    sub=csc # option for .nii format subcortex data pending\n",
    "    return left, right, sub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2119488d3bc0416eb2140e6ad82bd90a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# csplit example 1:\n",
    "a,b,c=csplit(data_ciftiTS,gii=True)\n",
    "aa=mwpad(a,data_lh_myelin,frame=0)\n",
    "bb=mwpad(b,data_rh_myelin,frame=0)\n",
    "fig2=showSurface(data_lh,aa,showZeroes=False)\n",
    "fig2=showSurface(data_rh,bb,showZeroes=False,newfigure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2668365002a64af996b9e55094aa9d09",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# csplit example 2:\n",
    "aaa=mwpad(a,data_lh_myelin,frame=0)\n",
    "fig3=showSurface(data_lh,aaa,showZeroes=False)\n",
    "for x in range(1,10):\n",
    "    time.sleep(1)\n",
    "    aaa=mwpad(a,data_lh_myelin,frame=x)\n",
    "    fig3=showSurface(data_lh,aaa,showZeroes=False,newfigure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7a266dbce63b4fc08bbe38f53031349b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# csplit example 2:\n",
    "thrl,thrr,thrs=csplit(nhw+'/thresh_demo.dtseries.nii',gii=True)\n",
    "aaa=mwpad(thrl,data_lh_myelin,frame=0)\n",
    "fig4=showSurface(data_lh,aaa,showZeroes=False)\n",
    "for x in range(1,10):\n",
    "    time.sleep(1)\n",
    "    aaa=mwpad(thrl,data_lh_myelin,frame=x)\n",
    "    fig4=showSurface(data_lh,aaa,showZeroes=False,newfigure=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'range4d' 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-40-9c1ece1cdd42>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0md\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcifti2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcifti2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mCifti2Image\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdataobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mrange4d\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mheader\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mch\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnifti_header\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextra\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfile_map\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mdl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcsplit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mgii\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mddl\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmwpad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdata_lh_myelin\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mfig4\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshowSurface\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata_lh\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maa\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mshowZeroes\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msleep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'range4d' is not defined"
     ]
    }
   ],
   "source": [
    "d=nb.cifti2.cifti2.Cifti2Image(dataobj=range4d, header=ch, nifti_header=cn, extra=cx, file_map=cf)\n",
    "dl,dr,ds=csplit(d,gii=True)\n",
    "ddl=mwpad(dl,data_lh_myelin)\n",
    "fig4=showSurface(data_lh,aa,showZeroes=False)\n",
    "time.sleep(5)\n",
    "fig4=showSurface(data_lh,ddl,showZeroes=False,newfigure=False,colormap='autumn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "03f92a164b814145a02cc5208d3b4ca5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import ciftitools as ct\n",
    "import showSurface_Ely as ss\n",
    "\n",
    "ss.showSurface(data_lh,data_lh_myelin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other stuff still in development:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 314,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cplot(coverlay,lsurf,rsurf,subvol,colormap):\n",
    "    # plots cifti overlay data with specified color on surfaces/subcortex\n",
    "    # lsurf and rsurf should be .surf.gii files matched to the number of vertices in each cifti hemisphere\n",
    "    # subvol should be a HCP-style subcortical volume (e.g. Atlas_ROIs.2.nii.gz)\n",
    "    \n",
    "    # import packages\n",
    "    import nibabel as nb\n",
    "    import matplotlib.pyplot as plt\n",
    "    import ipyvolume\n",
    "    import ipyvolume.pylab as p3\n",
    "    import numpy as np\n",
    "    \n",
    "    # read input files\n",
    "    ls=nb.load(lsurf)\n",
    "    rs=nb.load(rsurf)\n",
    "    sv=nb.load(subvol)\n",
    "    \n",
    "    # parse surface structure\n",
    "    lxyz=ls.darrays[0]\n",
    "    rxyz=rs.darrays[0]\n",
    "    lmesh=ls.darrays[1]\n",
    "    rmesh=rs.darrays[1]\n",
    "    \n",
    "    # plot\n",
    "    showSurface(lsurf,coverlay,showZeroes=False)\n",
    "    \n",
    "#    return 1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 658,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'overlay' 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-658-e59157abfc4d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# WARNING: still in development\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moverlay\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcifti2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcifti2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mCifti2Image\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      3\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m         \u001b[0mleft\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mright\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcsplit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moverlay\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m         \u001b[0mactivation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mleft\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'overlay' is not defined"
     ]
    }
   ],
   "source": [
    "    # WARNING: still in development\n",
    "    if isinstance(overlay,nb.cifti2.cifti2.Cifti2Image):\n",
    "        try:\n",
    "            left,right,sub=csplit(overlay)\n",
    "            activation=left[0]\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid cifti file.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}