{
 "metadata": {
  "name": "",
  "signature": "sha256:3f4b7bc771f0128cbe73102ba120f3716e53d18193f73312189e00c2ca171cbb"
 },
 "name": "optimizing_rotation",
 "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 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": []
    },
    {
     "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": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - nibabel package\n",
      "import nibabel as nib"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - our rotations module, now on the PYTHONPATH\n",
      "from rotations import x_rotmat, y_rotmat, z_rotmat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Finding a secret rotation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the last notebook - [moreon_rotation](moreon_rotation.ipynb) - we undid a rotation.\n",
      "\n",
      "In that case, I told you what the rotation was.  What if I didn't tell you?  Could you work it out?\n",
      "\n",
      "This is a good job for optimization."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - the first volume of ds107_sub012_t1r2.nii\n",
      "img_4d = nib.load('ds107_sub012_t1r2.nii')\n",
      "data = img_4d.get_data()\n",
      "vol0 = data[..., 0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here is the volume with some secret rotations:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - the secretly rotated image\n",
      "rotated_img = nib.load('secret_rotated_volume.nii')\n",
      "rotated_vol0 = rotated_img.get_data()\n",
      "rotated_vol0.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Displayed side by side with the original (not-rotated) image:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - slices on z, y, and x axis from \n",
      "fig, axes = plt.subplots(3, 2, figsize=(10, 15))\n",
      "axes[0, 0].imshow(vol0[:, :, 17])\n",
      "axes[0, 1].imshow(rotated_vol0[:, :, 17])\n",
      "axes[1, 0].imshow(vol0[:, 31, :])\n",
      "axes[1, 1].imshow(rotated_vol0[:, 31, :])\n",
      "axes[2, 0].imshow(vol0[31, :, :])\n",
      "axes[2, 1].imshow(rotated_vol0[31, :, :])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You might be able to work out what the transformations are, at least roughly, but can we work them out by optimization?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first thing we need for optimization, is a mismatch metric.\n",
      "\n",
      "I suggest you use the correlation mismatch function for the metric.  You can get that from the [optimizing_space](optimizing_space.ipynb) notebook, and paste it below.  But feel free to try another mismatch metric if you like.  Remember it is a function that accepts two images and returns a scalar that is low when the images are well matched."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the correlation mismatch metric from `optimizing_space`, paste here"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we need a function that will transform a given image by a given set of rotations.  The arguments will be:\n",
      "\n",
      "* `vol_arr` : the image that we will transform;\n",
      "* `rotations` - a vector of rotations.\n",
      "\n",
      "The function then returns a copy of `vol_arr` with those rotations applied.\n",
      "\n",
      "Specifically, the `rotations` give the *resampling transform* that gives the mapping from coordinates in a new empty copy of `vol_arr` - call this `K` - to coordinates in `vol_arr`.\n",
      "\n",
      "Our vector of rotations is length 3, containing the $r_x, r_y, r_z$, which are, respectively, the rotations about the x, y and z axis.\n",
      "\n",
      "In order to do this transformation, you will need to take the three parameters, and convert them to the corresponding rotation matrix.\n",
      "\n",
      "You then need to apply this transformation matrix to the coordinates of `K` to return a new copy of `vol_arr` with $r_x, r_y, r_z$ applied.  Of course, this will be a job for `affine_transform`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Your new function will look something like this:\n",
      "\n",
      "```\n",
      "def apply_rotations(vol_arr, rotations):\n",
      "    r_x, r_y, r_z = rotations\n",
      "    rotation_matrix = ?\n",
      "    # apply rotations with affine_transform to make new image\n",
      "    # return new image\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make apply_rotations function, accepting `vol_arr` and `rotations` vector\n",
      "# Returning image "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - a quick check that 0 rotations give the same output back\n",
      "not_changed = apply_rotations(rotated_vol0, [0, 0, 0])\n",
      "assert np.allclose(not_changed, rotated_vol0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we have the function to apply the rotations, and the matching function, we can make the cost function to optimize.\n",
      "\n",
      "Use the global scope to pick up and use `rotated_vol0` and `vol0` in the function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make a cost function, that\n",
      "# * is called 'cost_function'\n",
      "# * accepts a vector of rotations as input\n",
      "# * applies the vector of rotations to `rotated_vol0` from the global scope\n",
      "# * returns the mismatch metric for the transformed copy of `rotated_vol0` and `vol0`"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - a quick check the cost function returns the current value without rotations\n",
      "current = correl_mismatch(vol0, rotated_vol0)\n",
      "redone = cost_function([0, 0, 0])\n",
      "assert np.allclose(current, redone)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we are ready to optimize.  We are going to need at least one of the cost functions from `scipy.optimize`.\n",
      "\n",
      "`fmin_powell` is a good place to start:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - get fmin_powell\n",
      "from scipy.optimize import fmin_powell"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now call `fmin_powell` with some starting guess for the rotations:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Call optimizing function and collect best estimates for rotations"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, use these parameters to:\n",
      "\n",
      "* compile the rotation matrix from the optmized parameters.  This gives the matrix mapping from coordinates in `vol0` to coordinates in `rotated_vol0`;\n",
      "* Use `affine_transform` on `rotated_vol0` to get an un-rotated version of that image.\n",
      "\n",
      "You should be able to do this with your `apply_rotations` function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use 'apply_rotations' and the estimated parameters to un-rotate the rotated image\n",
      "# Put the new un-rotated image into a variable `best_vol0`"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now you can look at the original and the un-rotated image side by side:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - slices on z, y, and x axis from \n",
      "fig, axes = plt.subplots(3, 2, figsize=(10, 15))\n",
      "axes[0, 0].imshow(vol0[:, :, 17])\n",
      "axes[0, 1].imshow(best_vol0[:, :, 17])\n",
      "axes[1, 0].imshow(vol0[:, 31, :])\n",
      "axes[1, 1].imshow(best_vol0[:, 31, :])\n",
      "axes[2, 0].imshow(vol0[31, :, :])\n",
      "axes[2, 1].imshow(best_vol0[31, :, :])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you like, you can check in the file `make_rotated.py` to see what rotation I actually applied."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## If you raced ahead"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Why not try another mismatch function?\n",
      "\n",
      "And / or try another `fmin_` function from `scipy.optimize`?\n",
      "\n",
      "Is it easy for the optimization to go wrong, or is it nearly always right whatever parameter you choose?"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}