{ "metadata": { "name": "", "signature": "sha256:2a8d4d512d137fcf6818eed998bad870bb540dffb7c1385a0bb5f50af05eecdd" }, "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": "code", "collapsed": false, "input": [ "# - import numpy.linalg with a shorter name\n", "import numpy.linalg as npl" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Making floating points numbers print nicely" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default, when numpy prints an array, it looks for very small or very large numbers. If it finds either, it uses exponents to show the numbers. This can be annoying:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.pi" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 6, "text": [ "3.141592653589793" ] } ], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "np.array([np.pi, 0.000001])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "array([ 3.14159265e+00, 1.00000000e-06])" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to avoid this, you can tell numpy not to use exponential notation for small numbers:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.set_printoptions(suppress=True)\n", "np.array([np.pi, 0.000001])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 8, "text": [ "array([ 3.14159265, 0.000001 ])" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This setting stays in place until you change it:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.array([np.pi, 0.000001])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 9, "text": [ "array([ 3.14159265, 0.000001 ])" ] } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "It can also be annoying to see many digits after the decimal point, if you know that these are not important. You can set the number of digits after the decimal point for numpy printing like this:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.set_printoptions(precision=4)\n", "np.array([np.pi, 0.000001])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 10, "text": [ "array([ 3.1416, 0. ])" ] } ], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This only affects printing, not calculations:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.cos(np.pi)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 11, "text": [ "-1.0" ] } ], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "np.cos(3.1416)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 12, "text": [ "-0.99999999997301514" ] } ], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The image affine" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far we have not payed much attention to the image *header*.\n", "\n", "An image consists of:\n", "\n", "* the array data;\n", "* metadata (data about the array data).\n", "\n", "The header contains the metadata for the image.\n", "\n", "One piece of metadata, is the image affine." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import nibabel as nib\n", "img = nib.load('ds107_sub012_highres.nii')\n", "img.affine" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 13, "text": [ "array([[ 1. , 0. , 0. , -127. ],\n", " [ 0. , 1. , 0. , -83.3253],\n", " [ 0. , 0. , 1. , -90.0533],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can imagine, nibabel is getting the affine from the header:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print(img.header)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " object, endian='<'\n", "sizeof_hdr : 348\n", "data_type : \n", "db_name : \n", "extents : 0\n", "session_error : 0\n", "regular : r\n", "dim_info : 0\n", "dim : [ 3 256 208 192 1 1 1 1]\n", "intent_p1 : 0.0\n", "intent_p2 : 0.0\n", "intent_p3 : 0.0\n", "intent_code : none\n", "datatype : int16\n", "bitpix : 16\n", "slice_start : 0\n", "pixdim : [ 1. 1. 1. 1. 0. 0. 0. 0.]\n", "vox_offset : 0.0\n", "scl_slope : nan\n", "scl_inter : nan\n", "slice_end : 0\n", "slice_code : unknown\n", "xyzt_units : 10\n", "cal_max : 0.0\n", "cal_min : 0.0\n", "slice_duration : 0.0\n", "toffset : 0.0\n", "glmax : 0\n", "glmin : 0\n", "descrip : FSL4.0\n", "aux_file : \n", "qform_code : scanner\n", "sform_code : scanner\n", "quatern_b : 0.0\n", "quatern_c : 0.0\n", "quatern_d : 0.0\n", "qoffset_x : -127.0\n", "qoffset_y : -83.325302124\n", "qoffset_z : -90.0532836914\n", "srow_x : [ 1. 0. 0. -127.]\n", "srow_y : [ 0. 1. 0. -83.3253]\n", "srow_z : [ 0. 0. 1. -90.0533]\n", "intent_name : \n", "magic : n+1\n" ] } ], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "img.header.extensions" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 15, "text": [ "Nifti1Extensions()" ] } ], "prompt_number": 15 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice the `srow_x, srow_y, srow_z` fields in the header, that contain the affine for this image. It is not always this simple though - see http://nifti.nimh.nih.gov/nifti-1 for more details. In general, nibabel will take care of this for you, by extracting the affine from the header, and returning it via `img.affine`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Nifti images can also be `.img, .hdr` pairs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far, all the images we have seen have been NIfTI format images, stored in a single file with a `.nii` extension. The single file contains the header information, and the image array data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The NIfTI format also allows the image to be stored as two files, one with extension `.img` storing the image array data, and another with extension `.hdr` storing the header. These are called *NIfTI pair* images." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, there are two files in the same directory as this notebook, called `ds114_sub009_highres_moved.img` and `ds114_sub009_highres_moved.hdr`. These two files together form one NIfTI image. You can load these with nibabel in the usual way:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "pair_img = nib.load('ds114_sub009_highres_moved.img')\n", "pair_img.affine" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "array([[ 0.9416, -0.4311, -0.0586, -98.8336],\n", " [ 0.336 , 1.1887, 0.2264, -164.1377],\n", " [ -0.0215, -0.3028, 0.9723, -158.4178],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This form of the NIfTI image is getting less common, because it is inconvenient to have to keep the `.img` and `.hdr` files together, but you may still find them used. They have only one advantage, which is that, if some software wants to change only the header information, it only has to rewrite a small `.hdr` file, rather than the whole `.nii` file containing the image data and the header." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Identity matrices with `np.eye`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An identity matrix is a matrix with zeros off the diagonals, and ones on the leading diagonal. These are so common, than numpy (and MATLAB) have a function for making them. The parameter is the numbers of rows = number of columns:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.eye(3)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 17, "text": [ "array([[ 1., 0., 0.],\n", " [ 0., 1., 0.],\n", " [ 0., 0., 1.]])" ] } ], "prompt_number": 17 }, { "cell_type": "code", "collapsed": false, "input": [ "np.eye(6)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 18, "text": [ "array([[ 1., 0., 0., 0., 0., 0.],\n", " [ 0., 1., 0., 0., 0., 0.],\n", " [ 0., 0., 1., 0., 0., 0.],\n", " [ 0., 0., 0., 1., 0., 0.],\n", " [ 0., 0., 0., 0., 1., 0.],\n", " [ 0., 0., 0., 0., 0., 1.]])" ] } ], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Affines, inverses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We often have the situation where we compose an affine of several transformations. We do the composing using matrix multiplication. For example, here we compose two rotations and a translation. Remember - matrix multiplication works right to left:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# - get our functions to make rotation matrices\n", "from rotations import x_rotmat, y_rotmat, z_rotmat" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 19 }, { "cell_type": "markdown", "metadata": {}, "source": [ "A rotation matrix (3 x 3) for rotation -0.2 radians around the x axis:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "first_rotation = x_rotmat(-0.2)\n", "first_rotation" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 20, "text": [ "array([[ 1. , 0. , 0. ],\n", " [ 0. , 0.9801, 0.1987],\n", " [ 0. , -0.1987, 0.9801]])" ] } ], "prompt_number": 20 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can make this rotation matrix into an affine transformation, by putting it into the top left of a 4 x 4 identity matrix:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "first_affine = np.eye(4) # The identity affine\n", "first_affine[:3, :3] = first_rotation\n", "first_affine" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 21, "text": [ "array([[ 1. , 0. , 0. , 0. ],\n", " [ 0. , 0.9801, 0.1987, 0. ],\n", " [ 0. , -0.1987, 0.9801, 0. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we made a second affine matrix for a rotation around y of 0.4 radians:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "second_affine = np.eye(4)\n", "second_affine[:3, :3] = y_rotmat(0.4)\n", "second_affine" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 22, "text": [ "array([[ 0.9211, 0. , 0.3894, 0. ],\n", " [ 0. , 1. , 0. , 0. ],\n", " [-0.3894, 0. , 0.9211, 0. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 22 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we make a translation of 10 in x, 20 in y and 30 in z:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "third_affine = np.eye(4)\n", "third_affine[:3, 3] = [10, 20, 30]\n", "third_affine" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 23, "text": [ "array([[ 1., 0., 0., 10.],\n", " [ 0., 1., 0., 20.],\n", " [ 0., 0., 1., 30.],\n", " [ 0., 0., 0., 1.]])" ] } ], "prompt_number": 23 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We compose these three affine matrices to give an affine implementing *first* a rotation of -0.2 around the x axis, *then* a rotation of 0.4 around the y axis, and *finally* a translation [10, 20, 30] in [x, y, z]. Note the order though - matrix multiplication goes from right to left:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "combined = third_affine.dot(second_affine.dot(first_affine))\n", "combined" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 24, "text": [ "array([[ 0.9211, -0.0774, 0.3817, 10. ],\n", " [ 0. , 0.9801, 0.1987, 20. ],\n", " [ -0.3894, -0.183 , 0.9027, 30. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 24 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## nibabel.affines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In fact, nibabel has a short-cut routine to make a 4x4 affine matrix from a 3x3 matrix and an (optional) vector of translations:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Affine from a 3x3 matrix (the 'mat' in 'matvec')\n", "nib.affines.from_matvec(y_rotmat(0.4))" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 25, "text": [ "array([[ 0.9211, 0. , 0.3894, 0. ],\n", " [ 0. , 1. , 0. , 0. ],\n", " [-0.3894, 0. , 0.9211, 0. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 25 }, { "cell_type": "code", "collapsed": false, "input": [ "# Affine from a 3x3 matrix ('mat') and a translation vector ('vec')\n", "nib.affines.from_matvec(y_rotmat(0.4), [10, 20, 30])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 26, "text": [ "array([[ 0.9211, 0. , 0.3894, 10. ],\n", " [ 0. , 1. , 0. , 20. ],\n", " [ -0.3894, 0. , 0.9211, 30. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 26 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Manipulating affines with inverses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us say we have an affine, like the one we just made:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "combined" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 27, "text": [ "array([[ 0.9211, -0.0774, 0.3817, 10. ],\n", " [ 0. , 0.9801, 0.1987, 20. ],\n", " [ -0.3894, -0.183 , 0.9027, 30. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 27 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Imagine that we knew that this affine was composed of three affines, and we knew the last two, but not the first. How would we find what the first affine was?\n", "\n", "Call our combined affine $\\mathbf{D}$. We know that $\\mathbf{D} = \\mathbf{C} \\cdot \\mathbf{B} \\cdot \\mathbf{A}$. We know $\\mathbf{C}$ and $\\mathbf{B}$ but we want to find $\\mathbf{A}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Above I've written matrix multiplication with a dot - as in $\\mathbf{B} \\cdot \\mathbf{A}$, but in what follows I'll omit the dot, just writing $\\mathbf{B} \\mathbf{A}$ to mean matrix multiplication." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We find $\\mathbf{A}$ using matrix inverses. Call $\\mathbf{E} = \\mathbf{C} \\mathbf{B}$. Then $\\mathbf{D} = \\mathbf{E} \\mathbf{A}$. If we can find the inverse of $\\mathbf{E}$ (written as $\\mathbf{E^{-1}}$) then (by the definition of the inverse):\n", "\n", "$$\n", "\\mathbf{E^{-1}} \\mathbf{E} = \\mathbf{I}\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\mathbf{E^{-1}} \\mathbf{D} = \\mathbf{E^{-1}} \\mathbf{E} \\mathbf{A} \\\\\n", "\\mathbf{E^{-1}} \\mathbf{D} = \\mathbf{I} \\mathbf{A} \\\\\n", "\\mathbf{E^{-1}} \\mathbf{D} = \\mathbf{A}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For reasons we do not have time to go into, our affine matrices are almost invariably invertible.\n", "\n", "Let's see if we can reconstruct our `first_affine` from the `combined` affine, given we know the `third_affine` and `second_affine`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "E = third_affine.dot(second_affine)\n", "E_inv = npl.inv(E)\n", "E_inv.dot(combined)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 28, "text": [ "array([[ 1. , 0. , -0. , 0. ],\n", " [ 0. , 0.9801, 0.1987, 0. ],\n", " [ 0. , -0.1987, 0.9801, 0. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 28 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the same as our first affine:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "first_affine" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 29, "text": [ "array([[ 1. , 0. , 0. , 0. ],\n", " [ 0. , 0.9801, 0.1987, 0. ],\n", " [ 0. , -0.1987, 0.9801, 0. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 29 }, { "cell_type": "markdown", "metadata": {}, "source": [ "What about the situation where we know the first part of the affine, but we want to find the rest?\n", "\n", "To solve this problem, we will need the *right inverse*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inverse we have used so far is the *left inverse* - so called because we apply it multiplying on the left of the original matrix:\n", "\n", "$$\n", "\\mathbf{E^{-1}} \\mathbf{E} = \\mathbf{I}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Luckily, it turns out that, for square matrices, if there is a *left inverse* $\\mathbf{E^{-1}}$ then this is also the right inverse:\n", "\n", "$$\n", "\\mathbf{E^{-1}} \\mathbf{E} = \\mathbf{E} \\mathbf{E^{-1}} = \\mathbf{I}\n", "$$\n", "\n", "It is a little tricky to prove that there must be a right inverse, but it is easy to prove that, if there is a right inverse, it must be the same as the left inverse. Call the left inverse $\\mathbf{L}$ and the right inverse $\\mathbf{R}$:\n", "\n", "$$\n", "\\mathbf{LA} = \\mathbf{I}\\\\\n", "\\mathbf{AR} = \\mathbf{I}\\\\\n", "$$\n", "\n", "then:\n", "\n", "$$\n", "\\mathbf{LAR} = \\mathbf{LAR}\\\\\n", "\\mathbf{L(AR)} = \\mathbf{(LA)R}\\\\\n", "\\mathbf{L} = \\mathbf{R}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, in our case, where we want to find the transformations *following* the first affine, we can do this:\n", "\n", "$$\n", "\\mathbf{D} = \\mathbf{C} \\mathbf{B} \\mathbf{A} \\\\\n", "\\mathbf{D} \\mathbf{A^{-1}} = \\mathbf{C} \\mathbf{B} \\mathbf{A} \\mathbf{A^{-1}} \\\\\n", "\\mathbf{D} \\mathbf{A^{-1}} = \\mathbf{C} \\mathbf{B} \n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For our actual affines:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "third_with_second = combined.dot(npl.inv(first_affine))\n", "third_with_second" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 30, "text": [ "array([[ 0.9211, -0. , 0.3894, 10. ],\n", " [ 0. , 1. , 0. , 20. ],\n", " [ -0.3894, -0. , 0.9211, 30. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 30 }, { "cell_type": "code", "collapsed": false, "input": [ "# This is the same as\n", "third_affine.dot(second_affine)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 31, "text": [ "array([[ 0.9211, 0. , 0.3894, 10. ],\n", " [ 0. , 1. , 0. , 20. ],\n", " [ -0.3894, 0. , 0.9211, 30. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] } ], "prompt_number": 31 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }