{
 "metadata": {
  "name": "",
  "signature": "sha256:88b23cc92653863fbc531463bccd04c70f867edb44aa6f230bc19bcea11c6577"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Working with four dimensional images, with some revision"
     ]
    },
    {
     "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  # the Python array package\n",
      "import matplotlib.pyplot as plt  # the Python plotting package"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Reading text files"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We've been reading values from text files in the exercises.\n",
      "\n",
      "Here is some revision on how to do that, going from the crude to the elegant way."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I first write a little text file out to disk:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "numbers = [1.2, 2.3, 3.4, 4.5]\n",
      "fobj = open('some_numbers.txt', 'wt')\n",
      "for number in numbers:\n",
      "    fobj.write(str(number) + '\\n')  # String version of number, plus end-of-line character\n",
      "fobj.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now I read it back again.  First, I will just read the strings:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fobj = open('some_numbers.txt', 'rt')\n",
      "lines = fobj.readlines()\n",
      "len(lines)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lines[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next I will read the file, but converting each number to a float:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fobj = open('some_numbers.txt', 'rt')\n",
      "numbers_again = []\n",
      "for line in fobj.readlines():\n",
      "    numbers_again.append(float(line))\n",
      "numbers_again"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In fact we can do this even more simply by using `np.loadtxt`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.loadtxt('some_numbers.txt')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Loading data with nibabel"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This isn't really revision, because we haven't covered this properly yet."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Nibabel allows us to `load` images.  First we need to `import` the nibabel library:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import nibabel as nib"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First I `load` the image, go give me an \"image\" object:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img = nib.load('ds114_sub009_t2r1.nii')\n",
      "type(img)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can explore the image object with ``img.`` followed by TAB in the notebook."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Because images can have large arrays, nibabel does not load the image array when you `load` the image, in order to save time and memory.  The best way to get the image array data is with the `get_data()` method."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data = img.get_data()\n",
      "type(data)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `memmap` is a special type of array that saves memory, but otherwise behaves the same as any normal array.\n",
      "\n",
      "I recommend you disable the use of the 'memmap' special arrays by using the `mmap` *keyword argument* when you load the image."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img = nib.load('ds114_sub009_t2r1.nii', mmap=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this case you will get a normal array from `get_data`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data = img.get_data()\n",
      "type(data)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Four dimensional arrays - space + time"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The image we have just loaded is a four-dimensional image, with a four-dimensional array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first three axes represent three dimensional space.  The last axis represents time.  Here the last (time) axis has length 173.  This means that, for each of these 173 elements, there is one whole three dimensional image.  We often call the three-dimensional images *volumes*.  So we could say that this 4D image contains 173 volumes."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We have previously been taking slices over the third axis of a three-dimensional image.  We can now take slices over the fourth axis of this four-dimensional image:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "first_vol = data[:, :, :, 0]  # A slice over the final (time) axis"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This slice selects the first three-dimensional volume (3D image) from the 4D array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "first_vol.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can use the ellipsis `...` when slicing an array.  The ellipsis is a short cut for \"everything on the previous axes\".  For example, these two have exactly the same meaning:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "first_vol = data[:, :, :, 0]\n",
      "first_vol_again = data[..., 0]  # Using the ellipsis"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`first_vol` is a 3D image just like the 3D images you have already seen:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.imshow(first_vol[:, :, 14], cmap='gray')  # A slice over the third dimension"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Numpy operations work on the whole array by default"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Numpy operatons like `min`, and `max` and `std` operate on the whole numpy array by default, ignoring any array shape. For example, here is the maximum value for the whole 4D array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.max(data)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is exactly the same as:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.max(data.ravel())  # maximum when flattening the array to 1 dimension"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can ask numpy to operate over particular axes instead of operating over the whole array.  For example, this will generate a 3D image, where each array value is the variance the 173 values at that 3D position (the variance across time):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "var_vol = np.var(data, axis=3)  # variance across the final (time) axis\n",
      "plt.imshow(var_vol[:, :, 14], cmap='gray')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We'll come back to this in later classes"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Indexing with boolean arrays"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You probably came across this idea in the reading for the day2 homework - see https://scipy-lectures.github.io/intro/numpy/array_object.html#using-boolean-masks"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's say we have an array like this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arr = np.array([[0, 1, 3, 0], [5, 2, 0, 1]])\n",
      "arr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can get a True / False (boolean) array to tell us whether these values are above some threshold:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tf_array = arr > 2\n",
      "tf_array"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can use this boolean array to *index* into the original array (or any array with a suitable shape).  This will select only the elements where the boolean array is `True`.  The returned array may well have selected only a few members from any particular row or column or (in general) higher axis, so the returned array is always one-dimensional to reflect the loss of shape:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "selected_elements = arr[tf_array]\n",
      "selected_elements"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "selected_elements.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can use this to select values in our image as well.  For example, if we wanted to select only values less than 10 in `first_vol`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tf_lt_10 = first_vol < 10\n",
      "vals_lt_10 = first_vol[tf_lt_10]\n",
      "np.unique(vals_lt_10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<hr/>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Defining functions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A function takes input arguments, transforms them, and returns the output.  For example, the following function adds two numbers:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def add(a, b):\n",
      "    return a + b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Functions can also return multiple values, for example:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def weather():\n",
      "    return 'hot', 'north-east', 78\n",
      "\n",
      "sun, wind, temp = weather()\n",
      "\n",
      "print('Sun:', sun)\n",
      "print('Window:', wind)\n",
      "print('Temperature:', temp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}