2_SimpleITK_basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Images and Physical Space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Let's import matplotlib so we can view images\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# Also import NumPy for working with arrays\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The primary data type is the Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "\n",
    "image = sitk.Image()\n",
    "print(image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Use `ReadImage()` and `WriteImage()` for IO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "image = sitk.ReadImage('Data/qdna1.mha')\n",
    "\n",
    "sitk.WriteImage(image, 'qdna1.tif')\n",
    "!ls -l qdna1.tif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Spatial computations are always done in physical space\n",
    "\n",
    "This is very important to **avoid errors** when performing registration.\n",
    "\n",
    "![Physical space](Data/PhysicalSpace.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### *Origin*, *Spacing*, and *Direction* cosines define conversion from index to physical space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('Origin', image.GetOrigin())\n",
    "print('Spacing', image.GetSpacing())\n",
    "print('Direction', image.GetDirection())\n",
    "\n",
    "index = (0, 3)\n",
    "point = image.TransformIndexToPhysicalPoint(index)\n",
    "print('\\nindex', index, 'corresponds to point', point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "image.SetSpacing((3.0, 5.0))\n",
    "print('Spacing', image.GetSpacing())\n",
    "\n",
    "point = image.TransformIndexToPhysicalPoint(index)\n",
    "print('\\nindex', index, 'corresponds to point', point)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Conversion between NumPy and SimpleITK\n",
    "\n",
    "SimpleITK and NumPy indexing access is in opposite order! \n",
    "\n",
    "SimpleITK: image[x,y,z]<br>\n",
    "NumPy: image_numpy_array[z,y,x]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### From SimpleITK to NumPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "nda = sitk.GetArrayFromImage(image)\n",
    "print(image.GetSize())\n",
    "print(nda.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### From NumPy to SimpleITK\n",
    "\n",
    "Remember to to set the image's origin, spacing, and possibly direction cosine matrix. The default values may not match the physical dimensions of your image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "nda = np.zeros((10,20,3))\n",
    "\n",
    "# if this is supposed to be a 3D gray scale image [x=3, y=20, z=10]\n",
    "nda_as_image = sitk.GetImageFromArray(nda)\n",
    "print(nda_as_image.GetSize())\n",
    "\n",
    "# if this is supposed to be a 2D color image [x=20,y=10]\n",
    "nda_as_image = sitk.GetImageFromArray(nda, isVector=True)\n",
    "print(nda_as_image.GetSize())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Pixel Types\n",
    "\n",
    "The pixel type is represented as an enumerated type. The following is a table of the enumerated list.\n",
    "\n",
    "<table>\n",
    "  <tr><td>sitkUInt8</td><td>Unsigned 8 bit integer</td></tr>\n",
    "  <tr><td>sitkInt8</td><td>Signed 8 bit integer</td></tr>\n",
    "  <tr><td>sitkUInt16</td><td>Unsigned 16 bit integer</td></tr>\n",
    "  <tr><td>sitkInt16</td><td>Signed 16 bit integer</td></tr>\n",
    "  <tr><td>sitkUInt32</td><td>Unsigned 32 bit integer</td></tr>\n",
    "  <tr><td>sitkInt32</td><td>Signed 32 bit integer</td></tr>\n",
    "  <tr><td>sitkUInt64</td><td>Unsigned 64 bit integer</td></tr>\n",
    "  <tr><td>sitkInt64</td><td>Signed 64 bit integer</td></tr>\n",
    "  <tr><td>sitkFloat32</td><td>32 bit float</td></tr>\n",
    "  <tr><td>sitkFloat64</td><td>64 bit float</td></tr>\n",
    "  <tr><td>sitkComplexFloat32</td><td>complex number of 32 bit float</td></tr>\n",
    "  <tr><td>sitkComplexFloat64</td><td>complex number of 64 bit float</td></tr>\n",
    "  <tr><td>sitkVectorUInt8</td><td>Multi-component of unsigned 8 bit integer</td></tr>\n",
    "  <tr><td>sitkVectorInt8</td><td>Multi-component of signed 8 bit integer</td></tr>\n",
    "  <tr><td>sitkVectorUInt16</td><td>Multi-component of unsigned 16 bit integer</td></tr>\n",
    "  <tr><td>sitkVectorInt16</td><td>Multi-component of signed 16 bit integer</td></tr>\n",
    "  <tr><td>sitkVectorUInt32</td><td>Multi-component of unsigned 32 bit integer</td></tr>\n",
    "  <tr><td>sitkVectorInt32</td><td>Multi-component of signed 32 bit integer</td></tr>\n",
    "  <tr><td>sitkVectorUInt64</td><td>Multi-component of unsigned 64 bit integer</td></tr>\n",
    "  <tr><td>sitkVectorInt64</td><td>Multi-component of signed 64 bit integer</td></tr>\n",
    "  <tr><td>sitkVectorFloat32</td><td>Multi-component of 32 bit float</td></tr>\n",
    "  <tr><td>sitkVectorFloat64</td><td>Multi-component of 64 bit float</td></tr>\n",
    "  <tr><td>sitkLabelUInt8</td><td>RLE label of unsigned 8 bit integers</td></tr>\n",
    "  <tr><td>sitkLabelUInt16</td><td>RLE label of unsigned 16 bit integers</td></tr>\n",
    "  <tr><td>sitkLabelUInt32</td><td>RLE label of unsigned 32 bit integers</td></tr>\n",
    "  <tr><td>sitkLabelUInt64</td><td>RLE label of unsigned 64 bit integers</td></tr>\n",
    "</table>\n",
    "\n",
    "There is also `sitkUnknown`, which is used for undefined or erroneous pixel ID's. It has a value of -1.\n",
    "\n",
    "The 64-bit integer types are not available on all distributions. When not available the value is `sitkUnknown`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(image.GetPixelIDTypeAsString())\n",
    "print(image.GetPixelID())\n",
    "print(sitk.sitkUInt8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Image Construction\n",
    "\n",
    "There are a variety of ways to create an image. \n",
    "\n",
    "The following components are required for a complete definition of an image:\n",
    "<ol>\n",
    "<li>Pixel type [fixed on creation, no default]: unsigned 32 bit integer, sitkVectorUInt8, etc., see list above.</li>\n",
    "<li> Sizes [fixed on creation, no default]: number of pixels/voxels in each dimension. This quantity implicitly defines the image dimension.</li>\n",
    "<li> Origin [default is zero]: coordinates of the pixel/voxel with index (0,0,0) in physical units (i.e. mm).</li>\n",
    "<li> Spacing [default is one]: Distance between adjacent pixels/voxels in each dimension given in physical units.</li>\n",
    "<li> Direction matrix [default is identity]: mapping, rotation, between direction of the pixel/voxel axes and physical directions.</li>\n",
    "</ol>\n",
    "\n",
    "Initial pixel/voxel values are set to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "image_3D = sitk.Image(256, 128, 64, sitk.sitkInt16)\n",
    "image_2D = sitk.Image(64, 64, sitk.sitkFloat32)\n",
    "image_2D = sitk.Image([32,32], sitk.sitkUInt32)\n",
    "image_RGB = sitk.Image([128,64], sitk.sitkVectorUInt8, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Basic Image Attributes\n",
    "\n",
    "You can change the image origin, spacing and direction. Making such changes to an image already containing data should be done cautiously. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "image_3D.SetOrigin((78.0, 76.0, 77.0))\n",
    "image_3D.SetSpacing([0.5,0.5,3.0])\n",
    "\n",
    "print(image_3D.GetOrigin())\n",
    "print(image_3D.GetSize())\n",
    "print(image_3D.GetSpacing())\n",
    "print(image_3D.GetDirection())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Image dimension queries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(image_3D.GetDimension())\n",
    "print(image_3D.GetWidth())\n",
    "print(image_3D.GetHeight())\n",
    "print(image_3D.GetDepth())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "What is the depth of a 2D image?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(image_2D.GetSize())\n",
    "print(image_2D.GetDepth())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Pixel/voxel type queries: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(image_3D.GetPixelIDValue())\n",
    "print(image_3D.GetPixelIDTypeAsString())\n",
    "print(image_3D.GetNumberOfComponentsPerPixel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What is the dimension and size of a Vector image and its data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(image_RGB.GetDimension())\n",
    "print(image_RGB.GetSize())\n",
    "print(image_RGB.GetNumberOfComponentsPerPixel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Accessing Pixels and Slicing\n",
    "\n",
    "The Image class's member functions `GetPixel` and `SetPixel` provide an interface for pixel access."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "help(image_3D.GetPixel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(image_3D.GetPixel(0, 0, 0))\n",
    "image_3D.SetPixel(0, 0, 0, 1)\n",
    "print(image_3D.GetPixel(0, 0, 0))\n",
    "\n",
    "# This can also be done using pythonic notation.\n",
    "print(image_3D[0,0,1])\n",
    "image_3D[0,0,1] = 2\n",
    "print(image_3D[0,0,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Slicing of SimpleITK images returns a copy of the image data, but only a subset of pixels.\n",
    "\n",
    "The syntax is the same as NumPy arrays, which is very similar to MATLAB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the original image\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 2, 1)\n",
    "ax.imshow(sitk.GetArrayFromImage(image), cmap='gray')\n",
    "ax.set_title('Original')\n",
    "\n",
    "# Brute force sub-sample\n",
    "image_subsampled = image[::16,::16]\n",
    "ax = fig.add_subplot(1, 2, 2)\n",
    "ax.imshow(sitk.GetArrayFromImage(image_subsampled), cmap='gray')\n",
    "ax.set_title('Subsampled')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can take a 2D slice of a 3D image for visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Data from http://loci.wisc.edu/software/sample-data\n",
    "cells = sitk.ReadImage('Data/NESb_C2_TP1.tiff')\n",
    "print('Size', cells.GetSize())\n",
    "\n",
    "cells_slice = cells[:,:,4]\n",
    "print('Slice size', cells_slice.GetSize())\n",
    "\n",
    "plt.imshow(sitk.GetArrayFromImage(cells_slice), cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise 1: Image addition and physical space\n",
    "\n",
    "ITK does not allow pixel-based addition of images when they do not occupy the same location in physical space.\n",
    "\n",
    "Uncomment the addition operation below, and observe the error that occurs.\n",
    "\n",
    "Comment the `SetOrigin`, `SetSpacing`, and `SetDirection` lines.  Does the error go away? Why? Can any of these attributes be modified and addition still occur?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "image_a = sitk.Image(24, 24, sitk.sitkUInt8)\n",
    "image_a[0,0] = 2\n",
    "\n",
    "image_b = sitk.Image(image_a.GetSize(), sitk.sitkUInt8)\n",
    "image_b.SetDirection([0, 1, 0.5, 0.5])\n",
    "image_b.SetSpacing([0.5, 0.8])\n",
    "image_b.SetOrigin([10.0, 20.0])\n",
    "image_b[0,0] = 3\n",
    "\n",
    "#image_c = sitk.Add(image_a, image_b)\n",
    "#print(image_c[0,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise 2: Image physical extent\n",
    "\n",
    "In the figure below, two images are plotted at their physical location in space.\n",
    "\n",
    "By changing the *Origin* and *Spacing* of the image *b* only, make them occupy the same area in physical space.\n",
    "\n",
    "*Hint*: The length of an image with spacing, $\\Delta x$ and size $n$ is $n \\Delta x$. For two lengths to be the same, $\\Delta x_1 (n_1 - 1) = \\Delta x_2 (n_2  - 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def show_in_physical_space(image, column):\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(2, 2, column)\n",
    "    first_index = (0, 0)\n",
    "    left, bottom = image.TransformIndexToPhysicalPoint(first_index)\n",
    "    size = image.GetSize()\n",
    "    # Since indexes count from zero, we must subtract 1\n",
    "    last_index = (size[0] - 1, size[1] - 1)\n",
    "    right, top = image.TransformIndexToPhysicalPoint(last_index)\n",
    "    extent = (left, right, bottom, top)\n",
    "    print('extent', extent)\n",
    "    ax.imshow(sitk.GetArrayFromImage(image), extent=extent)\n",
    "    ax.set_xlim(0.0, 100.0)\n",
    "    ax.set_ylim(0.0, 100.0)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "image_a = sitk.Image(20, 20, sitk.sitkUInt8)\n",
    "image_a.SetOrigin([10.0, 20.0])\n",
    "image_a.SetSpacing([4.0, 2.0])\n",
    "\n",
    "show_in_physical_space(image_a, 1)\n",
    "\n",
    "image_b = sitk.Image(10, 10, sitk.sitkUInt8)\n",
    "image_b.SetOrigin([5.0, 5.0])\n",
    "image_b.SetSpacing([2.0, 1.0])\n",
    "\n",
    "show_in_physical_space(image_b, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Resampling Images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Define output image information\n",
    "\n",
    "* When resampling an image, we must first define output the image information, i.e. the *Origin*, *Spacing*, *Direction*, and *Size*.\n",
    "* This determines the grid of samples where the output will be sampled.\n",
    "* It does need not need to be the same as the input image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Resampling output information](Data/ResampleOutputInformation.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Finding a resampled pixel's value\n",
    "\n",
    "For each pixel in the output image:\n",
    "\n",
    "1. Transform its index value to its location in physical space.\n",
    "2. Apply an optional spatial transformation.\n",
    "3. Transform the physical point location to a continuous index on the image to be resampled.\n",
    "4. Interpolate from neighboring pixels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![Common reference system](Data/CommonReferenceSystem.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Different interpolators\n",
    "\n",
    "<a href=\"http://www.itk.org/SimpleITKDoxygen/html/namespaceitk_1_1simple.html#a7cb1ef8bd02c669c02ea2f9f5aa374e5\">There are many interpolators available</a>:\n",
    "<ul>\n",
    "<li> sitkNearestNeighbor </li>\n",
    "<li> sitkLinear </li>\n",
    "<li> sitkBSpline </li>\n",
    "<li> sitkGaussian </li>\n",
    "<li> sitkHammingWindowedSinc </li>\n",
    "<li> sitkCosineWindowedSinc </li>\n",
    "<li> sitkWelchWindowedSinc </li>\n",
    "<li> sitkLanczosWindowedSinc </li>\n",
    "<li> sitkBlackmanWindowedSinc </li>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Resampling and the registration framework\n",
    "\n",
    "**Registration** involves **sampling** the two images to be registered with at **spatial transform**. The sampled values are compared with a **similarity metric**.  An **optimizer** is used to improve the spatial transformation parameters iteratively until the samples coorespond.\n",
    "\n",
    "![Registration framework](Data/ITKv4RegistrationComponentsDiagram.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Tutorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Read a 3D input image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "input_image = sitk.ReadImage('Data/NESb_C2_TP1.tiff')\n",
    "print(input_image.GetSize())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "View one of its slices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "input_slice = input_image[:,:,4]\n",
    "plt.imshow(sitk.GetArrayFromImage(input_slice), cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Create a resampler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "resampler = sitk.ResampleImageFilter()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What are the options that can be set on the `resampler`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Uncomment the line below and place the cursor after the \"t\" in \"Set\"\n",
    "# Then, press the TAB key.\n",
    "#resampler.Set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can either set use a *reference image* to determine the output image pixel sampling grid, or we can set desired output image information manually.\n",
    "\n",
    "Here we will make the output image size half of the input image size in each direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "origin = input_image.GetOrigin()\n",
    "resampler.SetOutputOrigin(origin)\n",
    "\n",
    "direction = input_image.GetDirection()\n",
    "resampler.SetOutputDirection(direction)\n",
    "\n",
    "spacing = input_image.GetSpacing()\n",
    "resampler.SetOutputSpacing(spacing)\n",
    "\n",
    "size = list(input_image.GetSize())\n",
    "size[0] = int(size[0] / 2)\n",
    "size[1] = int(size[1] / 2)\n",
    "resampler.SetSize(size)\n",
    "print(size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "resampled_image = resampler.Execute(input_image)\n",
    "plt.imshow(sitk.GetArrayFromImage(resampled_image[:,:,4]), cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also apply a spatial transform in the resampling process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dimension = 3\n",
    "translation = (-100.0, -200.0, 0.0)\n",
    "transform = sitk.TranslationTransform(dimension, translation)\n",
    "resampler.SetTransform(transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "resampled_image = resampler.Execute(input_image)\n",
    "plt.imshow(sitk.GetArrayFromImage(resampled_image[:,:,4]), cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "When the input image is sampled outside of its domain, a *default pixel value* is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "resampler.SetDefaultPixelValue(150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "resampled_image = resampler.Execute(input_image)\n",
    "plt.imshow(sitk.GetArrayFromImage(resampled_image[:,:,4]), cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Finally, we can also change the interpolator used.  The default is a bi-linear interpolator for 2D images and a tri-linear interpolator for 3D images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "resampler.SetInterpolator(sitk.sitkNearestNeighbor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "resampled_image = resampler.Execute(input_image)\n",
    "plt.imshow(sitk.GetArrayFromImage(resampled_image[:,:,4]), cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise 1: Windowed Sinc\n",
    "\n",
    "Run the `resampler` with the `sitk.sitkWelchWindowedSinc` interpolator.  What is different?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise 2: Interpolator behavior\n",
    "\n",
    "The behavior of different interpolators is demonstrated on the <a href=\"http://www.cs.cornell.edu/~srm/publications/Vis94-filters-abstract.html\">Marschner-Lobb</a> image with the code below.\n",
    "\n",
    "Which interpolators before best? Is performance related to the experiences from Exercise 1?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def marschner_lobb(size=40, alpha=0.25, f_M=6.0):\n",
    "    img = sitk.PhysicalPointSource(sitk.sitkVectorFloat32, [size]*3, [-1]*3, [2.0/size]*3)\n",
    "    imgx = sitk.VectorIndexSelectionCast(img, 0)\n",
    "    imgy = sitk.VectorIndexSelectionCast(img, 1)\n",
    "    imgz = sitk.VectorIndexSelectionCast(img, 2)\n",
    "    del img\n",
    "    r = sitk.Sqrt(imgx**2 + imgy**2)\n",
    "    del imgx, imgy\n",
    "    pr = sitk.Cos((2.0*math.pi*f_M)*sitk.Cos((math.pi/2.0)*r))\n",
    "    ml = (1.0 - sitk.Sin((math.pi/2.0)*imgz) + alpha*(1.0+pr))/(2.0*(1.0+alpha))\n",
    "    return ml[:,:,ml.GetSize()[-1]//2]\n",
    "\n",
    "def myshow(image, title):\n",
    "    arr = sitk.GetArrayFromImage(image)\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.imshow(arr, cmap='gray')\n",
    "    ax.set_title(title)\n",
    "    \n",
    "def expand(image, interpolator):\n",
    "    return sitk.Expand(image, [5]*2, interpolator)\n",
    "\n",
    "interpolators = [('Nearest Neighbor', sitk.sitkNearestNeighbor),\n",
    "                 ('Linear', sitk.sitkLinear),\n",
    "                 ('B-Spline', sitk.sitkBSpline),\n",
    "                 ('Gaussian', sitk.sitkGaussian),\n",
    "                 ('Hamming Windowed Sinc', sitk.sitkHammingWindowedSinc),\n",
    "                 ('Blackman Windowed Sinc', sitk.sitkBlackmanWindowedSinc),\n",
    "                 ('Cosine Windowed Sinc', sitk.sitkCosineWindowedSinc),\n",
    "                 ('Welch Windowed Sinc', sitk.sitkWelchWindowedSinc),\n",
    "                 ('Lanczos Windowed Sinc', sitk.sitkLanczosWindowedSinc)]\n",
    "\n",
    "image = marschner_lobb()\n",
    "for interpolator in interpolators:\n",
    "    expanded = expand(image, interpolator[1])\n",
    "    myshow(expanded, interpolator[0])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Spatial Transformations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "from ipywidgets import interact, fixed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Transformation models\n",
    "\n",
    "There are many spatial transforms of varying complexity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Transform direction is from the resampled image to the input image\n",
    "\n",
    "![Transform Direction](Data/TransformDirection.png)\n",
    "\n",
    "This is because of the resampling process explained in the previous tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Transformation model parameter types\n",
    "\n",
    "Transform are defined by two sets of parameters, the **Parameters** and **FixedParameters**.  \n",
    "\n",
    "FixedParameters are not changed during the optimization process when performing registration.\n",
    "\n",
    "An example of FixedParameters is the center of rotation of a transformation modeling rigid body motion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Composite transforms\n",
    "\n",
    "It is possible to compose multiple transform together into a single transform object.\n",
    "\n",
    "With a **composite transform**, multiple resampling operations are prevented, so interpolation errors are not accumulated.\n",
    "\n",
    "Composite transforms are useful in registration because a series of increasing complex transformations can be incrementally optimized."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A number of different spatial transforms are available in SimpleITK.\n",
    "\n",
    "The simplest is the Identity Transform.  This transform simply returns input points unaltered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dimension = 2\n",
    "\n",
    "print('*Identity Transform*')\n",
    "identity = sitk.Transform(dimension, sitk.sitkIdentity)\n",
    "print('Dimension: ' + str(identity.GetDimension()))\n",
    "\n",
    "# Points are always defined in physical space\n",
    "point = (1.0, 1.0)\n",
    "def transform_point(transform, point):\n",
    "    transformed_point = transform.TransformPoint(point)\n",
    "    print('Point ' + str(point) + ' transformed is ' + str(transformed_point))\n",
    "\n",
    "transform_point(identity, point)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the TranslationTransform, the Parameters are the values of the translation Offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('*Translation Transform*')\n",
    "translation = sitk.TranslationTransform(dimension)\n",
    "\n",
    "print('Parameters: ' + str(translation.GetParameters()))\n",
    "print('Offset:     ' + str(translation.GetOffset()))\n",
    "print('FixedParameters: ' + str(translation.GetFixedParameters()))\n",
    "transform_point(translation, point)\n",
    "\n",
    "print('')\n",
    "translation.SetParameters((3.1, 4.4))\n",
    "print('Parameters: ' + str(translation.GetParameters()))\n",
    "transform_point(translation, point)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The affine transform is capable of representing translations, rotations, shearing, and scaling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print('*Affine Transform*')\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "print('Parameters: ' + str(affine.GetParameters()))\n",
    "print('FixedParameters: ' + str(affine.GetFixedParameters()))\n",
    "transform_point(affine, point)\n",
    "\n",
    "print('')\n",
    "affine.SetTranslation((3.1, 4.4))\n",
    "print('Parameters: ' + str(affine.GetParameters()))\n",
    "transform_point(affine, point)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Applying Transforms to Images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Create a function to display the images that is aware of image spacing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def myshow(img, title=None, margin=0.05, dpi=80):\n",
    "    nda = sitk.GetArrayFromImage(img)\n",
    "    spacing = img.GetSpacing()\n",
    "        \n",
    "    ysize = nda.shape[0]\n",
    "    xsize = nda.shape[1]\n",
    "      \n",
    "    figsize = (1 + margin) * ysize / dpi, (1 + margin) * xsize / dpi\n",
    "\n",
    "    fig = plt.figure(title, figsize=figsize, dpi=dpi)\n",
    "    ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])\n",
    "    \n",
    "    extent = (0, xsize*spacing[1], 0, ysize*spacing[0])\n",
    "    \n",
    "    t = ax.imshow(nda,\n",
    "            extent=extent,\n",
    "            interpolation='hamming',\n",
    "            cmap='gray',\n",
    "            origin='lower')\n",
    "    \n",
    "    if(title):\n",
    "        plt.title(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Create a grid image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "grid = sitk.GridSource(outputPixelType=sitk.sitkUInt16,\n",
    "    size=(250, 250),\n",
    "    sigma=(0.5, 0.5),\n",
    "    gridSpacing=(5.0, 5.0),\n",
    "    gridOffset=(0.0, 0.0),\n",
    "    spacing=(0.2,0.2))\n",
    "myshow(grid, 'Grid Input')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To apply the transform, a resampling operation is required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def resample(image, transform):\n",
    "    # Output image Origin, Spacing, Size, Direction are taken from the reference\n",
    "    # image in this call to Resample\n",
    "    reference_image = image\n",
    "    interpolator = sitk.sitkCosineWindowedSinc\n",
    "    default_value = 100.0\n",
    "    return sitk.Resample(image, reference_image, transform,\n",
    "                         interpolator, default_value)\n",
    "\n",
    "translation.SetOffset((3.1, 4.6))\n",
    "transform_point(translation, point)\n",
    "resampled = resample(grid, translation)\n",
    "myshow(resampled, 'Resampled Translation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What happened?  The translation is positive in both directions.  Why does the output image move down and to the left?  It important to keep in mind that a transform in a resampling operation defines *the transform from the output space to the input space*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "translation.SetOffset(-1*np.array(translation.GetParameters()))\n",
    "transform_point(translation, point)\n",
    "resampled = resample(grid, translation)\n",
    "myshow(resampled, 'Inverse Resampled')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "An affine (line preserving) transformation, can perform translation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def affine_translate(transform, x_translation=3.1, y_translation=4.6):\n",
    "    new_transform = sitk.AffineTransform(transform)\n",
    "    new_transform.SetTranslation((x_translation, y_translation))\n",
    "    resampled = resample(grid, new_transform)\n",
    "    myshow(resampled, 'Translated')\n",
    "    return new_transform\n",
    "    \n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "interact(affine_translate, transform=fixed(affine), x_translation=(-5.0, 5.0), y_translation=(-5.0, 5.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or scaling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def affine_scale(transform, x_scale=3.0, y_scale=0.7):\n",
    "    new_transform = sitk.AffineTransform(transform)\n",
    "    matrix = np.array(transform.GetMatrix()).reshape((dimension,dimension))\n",
    "    matrix[0,0] = x_scale\n",
    "    matrix[1,1] = y_scale\n",
    "    new_transform.SetMatrix(matrix.ravel())\n",
    "    resampled = resample(grid, new_transform)\n",
    "    myshow(resampled, 'Scaled')\n",
    "    print(matrix)\n",
    "    return new_transform\n",
    "\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "interact(affine_scale, transform=fixed(affine), x_scale=(0.2, 5.0), y_scale=(0.2, 5.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or rotation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def affine_rotate(transform, degrees=15.0):\n",
    "    parameters = np.array(transform.GetParameters())\n",
    "    new_transform = sitk.AffineTransform(transform)\n",
    "    matrix = np.array(transform.GetMatrix()).reshape((dimension,dimension))\n",
    "    radians = -np.pi * degrees / 180.\n",
    "    rotation = np.array([[np.cos(radians), -np.sin(radians)],[np.sin(radians), np.cos(radians)]])\n",
    "    new_matrix = np.dot(rotation, matrix)\n",
    "    new_transform.SetMatrix(new_matrix.ravel())\n",
    "    resampled = resample(grid, new_transform)\n",
    "    print(new_matrix)\n",
    "    myshow(resampled, 'Rotated')\n",
    "    return new_transform\n",
    "    \n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "interact(affine_rotate, transform=fixed(affine), degrees=(-90.0, 90.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or shearing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def affine_shear(transform, x_shear=0.3, y_shear=0.1):\n",
    "    new_transform = sitk.AffineTransform(transform)\n",
    "    matrix = np.array(transform.GetMatrix()).reshape((dimension,dimension))\n",
    "    matrix[0,1] = -x_shear\n",
    "    matrix[1,0] = -y_shear\n",
    "    new_transform.SetMatrix(matrix.ravel())\n",
    "    resampled = resample(grid, new_transform)\n",
    "    myshow(resampled, 'Sheared')\n",
    "    print(matrix)\n",
    "    return new_transform\n",
    "\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "\n",
    "interact(affine_shear, transform=fixed(affine), x_shear=(0.1, 2.0), y_shear=(0.1, 2.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Composite Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "It is possible to compose multiple transform together into a single transform object.  With a composite transform, multiple resampling operations are prevented, so interpolation errors are not accumulated.  For example, an affine transformation that consists of a translation and rotation,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "translate = (8.0, 16.0)\n",
    "rotate = 20.0\n",
    "\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "affine = affine_translate(affine, translate[0], translate[1])\n",
    "affine = affine_rotate(affine, rotate)\n",
    "\n",
    "resampled = resample(grid, affine)\n",
    "myshow(resampled, 'Single Transform')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "can also be represented with two Transform objects applied in sequence with a Composite Transform,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "translation = sitk.TranslationTransform(dimension)\n",
    "translation.SetOffset(-1*np.array(translate))\n",
    "composite = sitk.Transform(dimension, sitk.sitkComposite)\n",
    "composite.AddTransform(translation)\n",
    "affine = sitk.AffineTransform(dimension)\n",
    "affine = affine_rotate(affine, rotate)\n",
    "\n",
    "composite.AddTransform(translation)\n",
    "composite = sitk.Transform(dimension, sitk.sitkComposite)\n",
    "composite.AddTransform(affine)\n",
    "\n",
    "resampled = resample(grid, composite)\n",
    "myshow(resampled, 'Two Transforms')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "*Beware*, tranforms are noncommutative -- order matters!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "composite = sitk.Transform(dimension, sitk.sitkComposite)\n",
    "composite.AddTransform(affine)\n",
    "composite.AddTransform(translation)\n",
    "\n",
    "resampled = resample(grid, composite)\n",
    "myshow(resampled, 'Composite transform in reverse order')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise 1: Transform complexity\n",
    "\n",
    "Below is a list of transforms supported, ordered by their complexity. Choose three, and compare how many Parameters they have.\n",
    "\n",
    "\n",
    "<table width=\"100%\">\n",
    "<tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1TranslationTransform.html\">TranslationTransform</a></td><td>2D or 3D, translation</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1VersorTransform.html\">VersorTransform</a></td><td>3D, rotation represented by a versor</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1VersorRigid3DTransform.html\">VersorRigid3DTransform</a></td><td>3D, rigid transformation with rotation represented by a versor</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1Euler2DTransform.html\">Euler2DTransform</a></td><td>2D, rigid transformation with rotation represented by a Euler angle</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1Euler3DTransform.html\">Euler3DTransform</a></td><td>3D, rigid transformation with rotation represented by Euler angles</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1Similarity2DTransform.html\">Similarity2DTransform</a></td><td>2D, composition of isotropic scaling and rigid transformation with rotation represented by a Euler angle</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1Similarity3DTransform.html\">Similarity3DTransform</a></td><td>3D, composition of isotropic scaling and rigid transformation with rotation represented by a versor</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ScaleTransform.html\">ScaleTransform</a></td><td>2D or 3D, anisotropic scaling</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ScaleVersor3DTransform.html\">ScaleVersor3DTransform</a></td><td>3D, rigid transformation and anisotropic scale is <bf>added</bf> to the rotation matrix part (not composed as one would expect)</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ScaleSkewVersor3DTransform.html\">ScaleSkewVersor3DTransform</a></td><td>3D, rigid transformation with anisotropic scale and skew matrices <bf>added</bf> to the rotation matrix part (not composed as one would expect)</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1AffineTransform.html\">AffineTransform</a></td><td>2D or 3D, affine transformation.</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1BSplineTransform.html\">BSplineTransform</a></td><td>2D or 3D, deformable transformation represented by a sparse regular grid of control points. </td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/Doxygen/html/classitk_1_1DisplacementFieldTransform.html\">DisplacementFieldTransform</a></td><td>2D or 3D, deformable transformation represented as a dense regular grid of vectors.</td></tr>\n",
    "  <tr><td><a href=\"http://www.itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1Transform.html\">Transform</a></td>\n",
    "  <td>A generic transformation. Can represent any of the SimpleITK transformations, and a <b>composite transformation</b> (stack of transformations concatenated via composition, last added, first applied). </td></tr>\n",
    "  </table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}