{
 "metadata": {
  "name": "",
  "signature": "sha256:b410fb909e8f368098b206d4b4ac67bcbfb5ac14202fd619442acda6c2f47d57"
 },
 "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": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - show figures inside the notebook\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import common modules\n",
      "import numpy as np  # the Python array package\n",
      "import matplotlib.pyplot as plt  # the Python plotting package"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - set gray colormap and nearest neighbor interpolation by default\n",
      "plt.rcParams['image.cmap'] = 'gray'\n",
      "plt.rcParams['image.interpolation'] = 'nearest'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Rotations and rotation matrices"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "See [wikipedia on rotation matrices](https://en.wikipedia.org/wiki/Rotation_matrix)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In two dimensions, rotating a vector $\\theta$ around the origin can be expressed as a 2 by 2 transformation matrix:\n",
      "\n",
      "$$\n",
      "R(\\beta) = \\begin{bmatrix}\n",
      "\\cos \\beta & -\\sin \\beta \\\\\n",
      "\\sin \\beta & \\cos \\beta \\\\\n",
      "\\end{bmatrix}\n",
      "$$\n",
      "\n",
      "This rotates column vectors by matrix multiplication:\n",
      "\n",
      "$$\n",
      "\\begin{bmatrix}\n",
      "x' \\\\\n",
      "y' \\\\\n",
      "\\end{bmatrix} = \\begin{bmatrix}\n",
      "\\cos \\beta & -\\sin \\beta \\\\\n",
      "\\sin \\beta & \\cos \\beta \\\\\n",
      "\\end{bmatrix}\\begin{bmatrix}\n",
      "x \\\\\n",
      "y \\\\\n",
      "\\end{bmatrix}\n",
      "$$\n",
      "\n",
      "So the coordinates $(x',y')$ of the point $(x,y)$ after rotation are\n",
      "\n",
      "$$\n",
      "x' = x \\cos \\beta - y \\sin \\beta \\\\\n",
      "y' = x \\sin \\beta + y \\cos \\beta\n",
      "$$\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You might be able to see why this is true by looking at this diagram, and remembering:\n",
      "\n",
      "$$\n",
      "\\cos(\\alpha) = x \\\\\n",
      "\\sin(\\alpha) = y\n",
      "$$\n",
      "\n",
      "<img src=\"images/trig_addition.png\" width=\"40%\"/>\n",
      "\n",
      "(Image copyright [Stackexchange user Blue](http://math.stackexchange.com/users/409/blue) licensed with [cc-by-sa](http://creativecommons.org/licenses/by-sa/3.0))"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rotations in three dimensions simply extend from 2D.  For example, in 3D, the rotation above would be a rotation around the z axis (z stays constant, (x, y) change):"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "\\begin{alignat}{1}\n",
      "R_x(\\theta) &= \\begin{bmatrix}\n",
      "1 & 0 & 0 \\\\\n",
      "0 & \\cos \\theta &  -\\sin \\theta \\\\[3pt]\n",
      "0 & \\sin \\theta  &  \\cos \\theta \\\\[3pt]\n",
      "\\end{bmatrix} \\\\[6pt]\n",
      "R_y(\\theta) &= \\begin{bmatrix}\n",
      "\\cos \\theta & 0 & \\sin \\theta \\\\[3pt]\n",
      "0 & 1 & 0 \\\\[3pt]\n",
      "-\\sin \\theta & 0 & \\cos \\theta \\\\\n",
      "\\end{bmatrix} \\\\[6pt]\n",
      "R_z(\\theta) &= \\begin{bmatrix}\n",
      "\\cos \\theta &  -\\sin \\theta & 0 \\\\[3pt]\n",
      "\\sin \\theta & \\cos \\theta & 0\\\\[3pt]\n",
      "0 & 0 & 1\\\\\n",
      "\\end{bmatrix}\n",
      "\\end{alignat}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can combine rotations with matrix multiplication.  For example, here is an rotation around the x axis of $\\gamma$ radians:\n",
      "\n",
      "$$\n",
      "\\begin{bmatrix}\n",
      "x'\\\\\n",
      "y'\\\\\n",
      "z'\\\\\n",
      "\\end{bmatrix} =\n",
      "\\begin{bmatrix}\n",
      "1 & 0 & 0 \\\\\n",
      "0 & \\cos(\\gamma) & -\\sin(\\gamma) \\\\\n",
      "0 & \\sin(\\gamma) & \\cos(\\gamma) \\\\\n",
      "\\end{bmatrix}\n",
      "\\begin{bmatrix}\n",
      "x\\\\\n",
      "y\\\\\n",
      "z\\\\\n",
      "\\end{bmatrix}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We could then apply a rotation around the y axis of $\\phi$ radians:\n",
      "\n",
      "$$\n",
      "\\begin{bmatrix}\n",
      "x''\\\\\n",
      "y''\\\\\n",
      "z''\\\\\n",
      "\\end{bmatrix} =\n",
      "\\begin{bmatrix}\n",
      "\\cos(\\phi) & 0 & \\sin(\\phi) \\\\\n",
      "0 & 1 & 0 \\\\\n",
      "-\\sin(\\phi) & 0 & \\cos(\\phi) \\\\\n",
      "\\end{bmatrix}\n",
      "\\begin{bmatrix}\n",
      "x'\\\\\n",
      "j'\\\\\n",
      "k'\\\\\n",
      "\\end{bmatrix}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "But we could also write this as:\n",
      "\n",
      "$$\n",
      "\\begin{bmatrix}\n",
      "x''\\\\\n",
      "y''\\\\\n",
      "z''\\\\\n",
      "\\end{bmatrix} =\n",
      "\\begin{bmatrix}\n",
      "\\cos(\\phi) & 0 & \\sin(\\phi) \\\\\n",
      "0 & 1 & 0 \\\\\n",
      "-\\sin(\\phi) & 0 & \\cos(\\phi) \\\\\n",
      "\\end{bmatrix}\n",
      "\\begin{bmatrix}\n",
      "1 & 0 & 0 \\\\\n",
      "0 & \\cos(\\gamma) & -\\sin(\\gamma) \\\\\n",
      "0 & \\sin(\\gamma) & \\cos(\\gamma) \\\\\n",
      "\\end{bmatrix}\n",
      "\\begin{bmatrix}\n",
      "x\\\\\n",
      "y\\\\\n",
      "z\\\\\n",
      "\\end{bmatrix}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And, because matrix multiplication is associative:\n",
      "\n",
      "$$\n",
      "\\mathbf{Q} = \\begin{bmatrix}\n",
      "1 & 0 & 0 \\\\\n",
      "0 & \\cos(\\gamma) & -\\sin(\\gamma) \\\\\n",
      "0 & \\sin(\\gamma) & \\cos(\\gamma) \\\\\n",
      "\\end{bmatrix}\n",
      "$$\n",
      "\n",
      "$$\n",
      "\\mathbf{P} = \\begin{bmatrix}\n",
      "\\cos(\\phi) & 0 & \\sin(\\phi) \\\\\n",
      "0 & 1 & 0 \\\\\n",
      "-\\sin(\\phi) & 0 & \\cos(\\phi) \\\\\n",
      "\\end{bmatrix}\n",
      "$$\n",
      "\n",
      "$$\n",
      "\\mathbf{M} = \\mathbf{P} \\cdot \\mathbf{Q}\n",
      "$$\n",
      "\n",
      "$$\n",
      "\\begin{bmatrix}\n",
      "x''\\\\\n",
      "y''\\\\\n",
      "z''\\\\\n",
      "\\end{bmatrix} =\n",
      "\\mathbf{M}\n",
      "\\begin{bmatrix}\n",
      "x\\\\\n",
      "y\\\\\n",
      "z\\\\\n",
      "\\end{bmatrix}\n",
      "$$\n",
      "\n",
      "Here $\\mathbf{M}$ is our rotation matrix that encodes a rotation by $\\gamma$ radians around the x axis *followed by* (right to left matrix multiplication) a rotation by $\\phi$ radians around the y axis."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Have a look at the file `rotations.py` in `pna_code`.  It has routines that will make these 3 by 3 matrices for given angles of rotation around the x, y, and z axes."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Reampling, pull, push"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's say I have two images $\\mathbf{I}$ and $\\mathbf{J}$.\n",
      "\n",
      "There is some spatial transform between them, such as a translation, or rotation.\n",
      "\n",
      "We could either think of the transformation that maps $\\mathbf{I} \\to \\mathbf{J}$ or $\\mathbf{J} \\to \\mathbf{I}$.\n",
      "\n",
      "The transformation maps voxel coordinates in one image to coordinates in the other.\n",
      "\n",
      "For example, write a coordinate in $\\mathbf{I}$ as $(x_i, y_i, z_i)$, and a coordinate in $\\mathbf{J}$ as $(x_j, y_j, z_j)$.\n",
      "\n",
      "The mapping $\\mathbf{I} \\to \\mathbf{J}$ maps $(x_i, y_i, z_i) \\to (x_j, y_j, z_j)$.\n",
      "\n",
      "Now let's say that I want to move image $\\mathbf{J}$ to match image $\\mathbf{I}$.\n",
      "\n",
      "To do this moving, I need to resample $\\mathbf{J}$ onto the same voxel grid as $\\mathbf{I}$.\n",
      "\n",
      "Specifically, I am going to do the following:\n",
      "\n",
      "* make a new empty image $\\mathbf{K}$ that has the same voxel grid as $\\mathbf{I}$;\n",
      "* for each coordinate in $\\mathbf{I} : (x_i, y_i, z_i)$ I will apply the transform $\\mathbf{I} \\to \\mathbf{J}$ to get $(x_j, y_j, z_j)$;\n",
      "* I will probably need to estimate a value $v$ from $\\mathbf{J}$ at $(x_j, y_j, z_j)$ because the coordinate values $(x_j, y_j, z_j)$ will probably not be integers, and so there is no exactly matching value in $\\mathbf{J}$;\n",
      "* I place $v$ into $\\mathbf{K}$ at coordinate $(x_i, y_i, z_i)$.\n",
      "\n",
      "Notice that, in order to move $\\mathbf{J}$ to match $\\mathbf{I}$ I needed the opposite transform - that is $\\mathbf{I} \\to \\mathbf{J}$.  This is called *pull resampling*.\n",
      "\n",
      "I will call the $\\mathbf{I} \\to \\mathbf{J}$ transform the *resampling transform*."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Scipy ndimage and affine_transform"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Scipy has a function for doing reampling with transformations, called `scipy.ndimage.affine_transform`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.ndimage import affine_transform"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It does all the heavy work of resampling for us.\n",
      "\n",
      "For example, lets say I have an image volume $\\mathbf{I}$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import nibabel as nib\n",
      "img = nib.load('ds107_sub012_t1r2.nii')\n",
      "data = img.get_data()\n",
      "I = data[..., 0]  # I is the first volume"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And another image volume $\\mathbf{J}$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "J = data[..., 1]  # I is the second volume"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's say I know that the resampling transformation $\\mathbf{I} \\to \\mathbf{J}$ is given by:\n",
      "\n",
      "* A rotation by 0.2 radians about the x axis\n",
      "* A translation of [1, 2, 3] voxels."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from rotations import x_rotmat  # rotations from `pna_code`\n",
      "M = x_rotmat(0.2)  # rotation matrix for rotation of 0.2 radians around x axis\n",
      "translation = [1, 2, 3]  # Translation from I to J\n",
      "M, translation"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "(array([[ 1.        ,  0.        ,  0.        ],\n",
        "        [ 0.        ,  0.98006658, -0.19866933],\n",
        "        [ 0.        ,  0.19866933,  0.98006658]]), [1, 2, 3])"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`affine_transform` now does the work for me:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "K = affine_transform(J, M, translation, order=1)  # order=1 for linear interpolation\n",
      "K.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "(64, 64, 35)"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here `affine_transform` \n",
      "\n",
      "* makes the new empty volume `K`, assuming it will be the same shape as `J`;\n",
      "* for each coordinate $(x_i, y_i, z_i)$ implied by the volume `K`:\n",
      "    * applies the transformations implied by `M` and `translation` to $(x_i, y_i, z_i)$ to get the corresponding point in `J` : $(x_j, y_j, z_j)$;\n",
      "    * resamples `J` at $(x_j, y_j, z_j)$ to get $v$;\n",
      "    * places $v$ at $(x_i, y_i, z_i)$ in `K`"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}