{
 "metadata": {
  "name": "",
  "signature": "sha256:9e42cc11df53c5b9633ffdd7c51838c9133638d7743792a43e8665e57eb1a95f"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compatibility with Python 3\n",
      "from __future__ import print_function  # print('me') instead of print 'me'\n",
      "from __future__ import division  # 1/2 == 0.5, not 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Show figures inside the notebook\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np  # Python array library\n",
      "import matplotlib.pyplot as plt  # Python plotting library"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The major arteries in a T1 MRI image often have high signal (white when displaying in grayscale).\n",
      "\n",
      "Our task is to see if we can pick out the courses of the [vertebral, basilar](http://en.wikipedia.org/wiki/Vertebral_artery#mediaviewer/File:Vertebral_artery_3D_AP.jpg) and [carotid](http://en.wikipedia.org/wiki/Internal_carotid_artery#mediaviewer/File:Magnetic_resonance_angiogram_of_segments_of_the_internal_carotid_artery.jpg) arteries on this image."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This time we are going to load an image using the `nibabel` library.  The image is `ds107_sub001_highres.nii` in the same directory as this notebook."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import nibabel as nib"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Try loading the image using nibabel, to make an image object. Use tab completion on `nib.` or look at the the [day0 introduction notebook](http://nbviewer.ipython.org/github/practical-neuroimaging/pna2015/blob/master/day0/what_is_an_image.ipynb) to work out how to do this."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# img = ?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now get the image array data from the nibabel image object.  Don't forget to use tab completion on the image object if you can remember or don't know the methods of the object."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# data = ?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Try plotting a few slices over the third dimension to see whether you can see the arteries.  For example, to plot the first slice over the third dimension, you might use:\n",
      "\n",
      "`plt.imshow(data[:, :, 0], cmap='gray')`\n",
      "\n",
      "\n",
      "where `data` is your image array data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plotting some slices over the third dimension"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now try looking for a good threshold so that you pick up only the very high signal in the brain.  A good place to start is to use ``plt.hist`` to get an idea of the spread of values within the volume and within the slice."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Here you might try plt.hist or something else to find a threshold"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Try making a binarized image with your threshold and displaying slices with that.  What structures are you picking up?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Maybe display some slices from the data binarized with a threshold"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now try taking a 3D subvolume out of the middle of the image (the approximate middle in all three axes) to pick out a good subvolume of the image that still contains the big arteries."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create a smaller 3D subvolume from the image data that still contains the arteries"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Try binarizing that with some thresholds to see whether you can pick out the arteries without much other stuff.  Hint - you might consider using `np.percentile` or `plt.hist` to find a good threshold."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Try a few plots of binarized slices and other stuff to find a good threshold"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you have a good threshold and a good binarized subset, see if you can see the arterial structure using the fancy function to plot the binarized image with a 3D rendering."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# This function uses matplotlib 3D plotting and sckit-image for rendering\n",
      "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n",
      "from skimage import measure\n",
      "\n",
      "def binarized_surface(binary_array):\n",
      "    \"\"\" Do a 3D plot of the surfaces in a binarized image\n",
      "    \n",
      "    This uses scikit-image and some fancy commands that we don't\n",
      "    need to worry about at the moment, to do the plot.\n",
      "    \"\"\"\n",
      "    verts, faces = measure.marching_cubes(binary_array, 0)\n",
      "    fig = plt.figure(figsize=(10, 12))\n",
      "    ax = fig.add_subplot(111, projection='3d')\n",
      "\n",
      "    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n",
      "    mesh = Poly3DCollection(verts[faces], linewidths=0, alpha=0.5)\n",
      "    ax.add_collection3d(mesh)\n",
      "    ax.set_xlim(0, binary_array.shape[0])\n",
      "    ax.set_ylim(0, binary_array.shape[1])\n",
      "    ax.set_zlim(0, binary_array.shape[2])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For example, let's say you have a binarized subvolume of the original data called ``binarized_subvolume``.  You could do a 3D rendering of this binary image with:\n",
      "\n",
      "`binarized_surface(binarized_subvolume)`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# binarized_surface(binarized_subvolume)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}