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": 25,
   "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\n",
    "import sys\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": 26,
   "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",
    "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>Define functions:</b>  \n",
    "(kindly remind us to define oscale and mwmask before running)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def showSurface(surface,overlay=None,frame=0,newfigure=True,colormap='summer',figsize=np.array([600,600]),\n",
    "                figlims=np.array([[-100,100],[-100,100],[-100,100]]),mwall=True):\n",
    "    '''\n",
    "    Displays a surface mesh in gifti or FreeSurfer (FS) surface format with/without an overlay inside  \n",
    "    Jupyter notebook for interactive visualization.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    surface: str, gifti 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",
    "        indice of the frame (timepoint or functional data frame) to show\n",
    "    newfigure: bool\n",
    "        Create a new figure else prints into the last figure (in order to visualize both hemispheres in \n",
    "        one plot), default True\n",
    "    colormap: string\n",
    "        A matplotlib colormap, default summer\n",
    "    figsize: ndarray\n",
    "        Size of the figure to display, default [600,600]\n",
    "    figLims: ndarray\n",
    "        x,y and z limits of the axes, default [[-100,100],[-100,100],[-100,100]])\n",
    "\n",
    "    '''\n",
    "\n",
    "    if isinstance(surface,str):\n",
    "        if not os.path.exists(surface):\n",
    "            error('File does not exist, please provide a valid file path to a gifti or FreeSurfer file.')\n",
    "        filename, file_extension = os.path.splitext(surface)\n",
    "        if file_extension is '.gii':\n",
    "            surface = nb.load(surface)\n",
    "        else:\n",
    "            fsgeometry = nb.freesurfer.read_geometry(surface)\n",
    "            x,y,z = fsgeometry[0].T\n",
    "            vertex_edges=fsgeometry[1]\n",
    "\n",
    "    if isinstance(surface,nb.gifti.gifti.GiftiImage):\n",
    "        try:\n",
    "            vertex_spatial=surface.darrays[0]\n",
    "            vertex_edges=surface.darrays[1]\n",
    "            x, y, z = vertex_spatial.data.T\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid gifti file.')\n",
    "\n",
    "    \n",
    "    if mwall is False:\n",
    "        try:\n",
    "            mkeep,mkill=mwmask(surface,overlay)\n",
    "        except:\n",
    "            raise ValueError('Overlay required for medial wall masking.')\n",
    "    \n",
    "    if not isinstance(frame,int):\n",
    "        ValueError('Please provide a valid integer frame index.')\n",
    "    \n",
    "    if isinstance(overlay,list):\n",
    "        if frame>len(overlay) or frame < 0:\n",
    "            error('Frame index out of bounds, please provide a valid frame index.')\n",
    "        overlay = overlay[frame]\n",
    "    \n",
    "    if isinstance(overlay,str):\n",
    "        if not os.path.exists(overlay):\n",
    "            error('File does not exist, please provide a valid file path to a gifti or FreeSurfer file.')\n",
    "        filename, file_extension = os.path.splitext(overlay)\n",
    "        \n",
    "        if file_extension is '.gii':\n",
    "            overlay = nb.load(overlay)\n",
    "        elif (file_extension in ('.annot','')):\n",
    "            annot = nb.freesurfer.read_annot(overlay)\n",
    "            activation = annot[0]\n",
    "        elif (file_extension in ('.curv','.thickness','.sulc')):\n",
    "            activation = nb.freesurfer.read_morph_data(overlay)\n",
    "        elif (file_extension in ('.dlabel.nii','.dscalar.nii','.dconn.nii','.dtseries.nii')):\n",
    "            left,right,sub=csplit(overlay)\n",
    "            \n",
    "    if isinstance(overlay,nb.gifti.gifti.GiftiImage):\n",
    "        try:\n",
    "            activation=overlay.darrays[0].data\n",
    "        except:\n",
    "            raise ValueError('Please provide a valid gifti file')\n",
    "            \n",
    "    # WARNING: still in development\n",
    "    if isinstance(overlay,nb.cifti2.cifti2.Cifti2Image):\n",
    "        try:\n",
    "            left,right,sub=csplit(overlay)\n",
    "            activation=left.copy() \n",
    "        except:\n",
    "            raise ValueError('Please provide a valid cifti file.')\n",
    "    \n",
    "    \n",
    "    if newfigure:\n",
    "\n",
    "        fig = p3.figure(width=figsize[0], height=figsize[1])\n",
    "        fig.camera_fov = 1\n",
    "        fig.style = {'axes': {'color': 'black',\n",
    "          'label': {'color': 'black'},\n",
    "          'ticklabel': {'color': 'black'},\n",
    "          'visible': False},\n",
    "         'background-color': 'white',\n",
    "         'box': {'visible': False}}\n",
    "        fig.xlim = (figlims[0][0], figlims[0][1])\n",
    "        fig.ylim = (figlims[1][0], figlims[1][1])\n",
    "        fig.zlim = (figlims[2][0], figlims[2][1])\n",
    "\n",
    "        \n",
    "    # plot surface\n",
    "    if overlay is None:\n",
    "        p3.plot_trisurf(x, y, z, triangles=vertex_edges.data)\n",
    "    elif mwall is False:\n",
    "        my_color = plt.cm.get_cmap(colormap)\n",
    "        colors=my_color(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": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mwmask(surf,mask):\n",
    "    '''\n",
    "    mask out medial wall; with some tinkering could be used to mask any vertices of (no) interest\n",
    "    surf and mask should be .gii or cifti structures loaded with nibabel\n",
    "    surf is the target surface\n",
    "    mask is an overlay with no information at medial wall vertices (e.g. myelin map)\n",
    "    '''\n",
    "    keep=(mask.darrays[0].data!=0)\n",
    "    kill=(mask.darrays[0].data==0)\n",
    "    ikeep=[i for i, e in enumerate(keep) if e != 0]\n",
    "    ikill=[i for i, e in enumerate(kill) if e != 0]\n",
    "    killdict={ii:1 for ii in ikill} # fun fact, iterating over a dictionary is ~exponentially faster vs. over a list\n",
    "    mask_kill=np.zeros([surf.darrays[1].data.shape[0]],dtype=bool) # create empty arrays matching surface mesh dimentions\n",
    "    mask_keep=mask_kill.copy()\n",
    "    for ii, row in enumerate(surf.darrays[1].data):\n",
    "        for item in row:\n",
    "            if item in killdict.keys():\n",
    "                mask_kill[ii]=True\n",
    "                continue\n",
    "            else:\n",
    "                mask_keep[ii]=True\n",
    "                continue\n",
    "    return mask_keep,mask_kill"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def oscale(olay):\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",
    "    scaled=np.asarray((olay-min(olay))/(max(olay)-min(olay)))\n",
    "    return scaled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a60c8a052fec4fbba7eb03428864f461",
       "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": "7b94c2ae3e8e4479a1e729be50fb532d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# showSurface example:\n",
    "showSurface(data_rh) # right surface only\n",
    "showSurface(data_rh,data_rh_myelin) # new plot, same surface but with myelin map overlay\n",
    "showSurface(data_lh,data_lh_myelin,newfigure=False,mwall=False) # same plot, added \"very inflated\" left surface, excluded medial wall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>NiWidgets 4D longitudinal or timeseries display:</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "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": 32,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c59ea9e38c3443edb3011504e617bdad",
       "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": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d28780eed5e14b1597dda56169fd3cee",
       "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": "02dc5e2adb0c4fcba9b8c8cc0aa0681d",
       "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>Tangentially related CIFTI tools:</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cread(infile):\n",
    "    # loads a cifti file and stores the data and four metadata objects\n",
    "    # should be called as: cd,ch,cn,cx,cf=cread(input_cifti)\n",
    "    cin=nb.load(infile)\n",
    "    cdata=np.asarray(cin.get_data()).copy()\n",
    "    chead=cin.header.copy()\n",
    "    cnhdr=cin.nifti_header.copy()\n",
    "    cxtra=cin.extra.copy()\n",
    "    cfmap=cin.file_map.copy()\n",
    "    return cdata,chead,cnhdr,cxtra,cfmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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:\n",
    "input_cifti=nhw+'/rfMRI_REST1_LR_Atlas_MSMAll_hp2000_clean.dtseries.nii'\n",
    "cd,ch,cn,cx,cf=cread(input_cifti)\n",
    "print(cd.data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def csave(filename,cdata,chead=None,cnhdr=None,cxtra=None,cfmap=None):\n",
    "    # saves a cifti data structure to file. Output likely unusable without some header info\n",
    "    cout=nb.Cifti2Image(dataobj=cdata,header=chead,nifti_header=cnhdr, extra=cxtra, file_map=cfmap)\n",
    "    cout.to_filename(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "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": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def clthresh(data,thresh,setto=0):\n",
    "    # set all values below threshold to a given value (default 0)\n",
    "    lthr_data=data.copy()\n",
    "    lthr_data[lthr_data<thresh]=setto\n",
    "    return lthr_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cuthresh(data,thresh,setto=0):\n",
    "    # set all values above threshold to a given value (default 0)\n",
    "    uthr_data=data.copy()\n",
    "    uthr_data[uthr_data>thresh]=setto\n",
    "    return uthr_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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 0x1221027f0>"
      ]
     },
     "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": 17,
   "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 0x11b0cd208>"
      ]
     },
     "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,ch,cf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cmask(cdata1,cdata2):\n",
    "    # masks cifti data 1 with non-zero values of cifti data 2\n",
    "    # cdata2 can be a single timepoint or matched to cdata1 dimensions\n",
    "    masked_cdata1=cdata1.copy()\n",
    "    if cdata1.shape != cdata2.shape:\n",
    "        mask=np.tile(cdata2,(cdata1.shape[0],1))\n",
    "    else:\n",
    "        mask=cdata2.copy()\n",
    "    masked_cdata1[mask==0]=0\n",
    "    return masked_cdata1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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 0x11cb5afd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cmask example:\n",
    "plt.hist(cd[np.nonzero(cd)],64,[0,20000]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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 0x11d191828>"
      ]
     },
     "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": 21,
   "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 0x11cb537b8>"
      ]
     },
     "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": 22,
   "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 0x11e443da0>"
      ]
     },
     "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": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cvoxels(infile,xyz=False):\n",
    "    # returns array where each element is an array of voxel coordinates for a particular structure\n",
    "    # setting xyz=True will return MNI coordinate array\n",
    "    cd,ch,cn,cx,cf=cread(infile)\n",
    "    ch1=ch.get_index_map(1)\n",
    "    voxarray=[]\n",
    "    if xyz is True:\n",
    "        from nibabel.affines import apply_affine\n",
    "        warp=ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix\n",
    "    for idx, bm in enumerate(mim.brain_models):\n",
    "        if bm.model_type=='CIFTI_MODEL_TYPE_VOXELS':\n",
    "            if xyz is False:\n",
    "                voxarray.append(bm.voxel_indices_ijk)\n",
    "            else:\n",
    "                voxarray.append(apply_affine(warp,bm.voxel_indices_ijk))\n",
    "        else:\n",
    "            voxarray.append([])\n",
    "    return np.asarray(voxarray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cinfo(infile):\n",
    "    # returns information about input cifti file, particularly header metadata\n",
    "    cd,ch,cn,cx,cf=cread(infile)\n",
    "    ch0=ch.get_index_map(0)\n",
    "    ch1=ch.get_index_map(1)\n",
    "    # useful outputs:\n",
    "    print('cifti data shape = {x}'.format(x=cd.shape))\n",
    "    print('')\n",
    "    if ch.get_index_map(0).indices_map_to_data_type=='CIFTI_INDEX_TYPE_SERIES':\n",
    "        print('cifti data dype = fMRI timeseries')\n",
    "        me=ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.meter_exponent\n",
    "        print('Xdim = {x}x10^{e}m'.format(\n",
    "            x=np.abs(ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix[0,0]),\n",
    "            e=me))\n",
    "        \n",
    "        print('     {x} timepoints, TR = {y} {z}'.format(\n",
    "            x=ch.get_index_map(0).number_of_series_points,\n",
    "            y=ch.get_index_map(0).series_step,\n",
    "            z=ch.get_index_map(0).series_unit))\n",
    "\n",
    "    for idx, bm in enumerate(ch1.brain_models):\n",
    "        print((idx, bm.index_offset, bm.index_count, bm.brain_structure)) # list brain structures, starting indices, and sizes\n",
    "    print('Voxel to MNI transformation matrix:')\n",
    "    print(ch1.volume.transformation_matrix_voxel_indices_ijk_to_xyz.matrix)\n",
    "    \n",
    "    # not returned:\n",
    "#    ch.may_contain_header() # unclear how to use\n",
    "#    bm.model_type # returns 'CIFTI_MODEL_TYPE_SURFACE' or 'CIFTI_MODEL_TYPE_VOXELS'\n",
    "#    bm.surface_number_of_vertices # returns total vertices in surface (includes medial wall)\n",
    "#    print('cifti header number of mapped indices = {x}'.format(x=ch.number_of_mapped_indices))\n",
    "#    print('cifti header mapped indices = {x}'.format(x=ch.mapped_indices))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other stuff still in development:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def csplit(incifti):\n",
    "    # split cifti overlay into separate left surface/right surface/subcortical voxel values\n",
    "    # incifti can be a file location string or loaded cifti object\n",
    "    \n",
    "    if isinstance(incifti,nb.cifti2.cifti2.Cifti2Image):\n",
    "        cd=np.asarray(incifti.get_data()).copy()\n",
    "        ch=incifti.header.copy()\n",
    "    else:\n",
    "        cd,ch,cn,cx,cf=cread(incifti) \n",
    "    \n",
    "    mim = ch.matrix[1]\n",
    "    bm1 = list(mim.brain_models)[0]\n",
    "    bm2 = list(mim.brain_models)[1]\n",
    "    lidx = list(bm1.vertex_indices)\n",
    "    ridx = [bm1.surface_number_of_vertices + val for val in bm2.vertex_indices]\n",
    "    bidx = np.concatenate((lidx, ridx))\n",
    "    # debug:\n",
    "    # print(bidx.shape)\n",
    "    #    for idx, bm in enumerate(mim.brain_models):\n",
    "    #        print((idx, bm.index_offset, bm.brain_structure))\n",
    "    \n",
    "    # split cifti overlay into left surface, right surface, and subcortical volume\n",
    "    clh=cd[:,0:len(lidx)]\n",
    "    crh=cd[:,len(lidx):len(lidx)+len(ridx)+1]\n",
    "    csc=cd[:,len(lidx)+len(ridx)+1:cd.shape[1]]\n",
    "    return clh, crh, csc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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,mwall=False)\n",
    "    \n",
    "#    return 1\n",
    "    "
   ]
  }
 ],
 "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
}