02_images_and_resampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">SimpleITK Images and Resampling</h1>\n",
    "\n",
    "\n",
    "**Summary:**    \n",
    "\n",
    "1. Images occupy a region in physical space which is defined by:\n",
    "  * Origin.\n",
    "  * Size (number of pixels per dimension).\n",
    "  * Spacing (unknown consistent units: nm, mm, m, km...).\n",
    "  * Direction cosine matrix (axis directions in physical space).\n",
    "  \n",
    "  These attributes are the image's meta-data. Computing the physical coordinates from image indexes requires all four   components.\n",
    "\n",
    "2. An image may contain a meta-data dictionary. This supplemental information often includes the image modality (e.g. CT), patient name, and information with respect to the image acquisition. \n",
    "3. Image initialization: user specified pixel type, user specified dimensionality (2,3), origin at zero, unit spacing in all dimensions and identity direction cosine matrix, intensities set to zero.\n",
    "4. Data transfer to/from numpy: GetArrayFromImage (copy), GetArrayViewFromImage (immutable), GetImageFromArray (copy) + set the meta-data yourself. \n",
    "5. A common issue with resampling resulting in an all black image is due to (a) incorrect specification of the \n",
    "desired output image's spatial domain (its meta-data); or (b) the use of the inverse of the transformation mapping from the output spatial domain to the resampled image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Images are Physical Objects\n",
    "\n",
    "<img src=\"figures/ImageOriginAndSpacing.png\" style=\"width:700px\"/><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.\n",
    "\n",
    "Some filters only work with images with a specific pixel type. The primary example is the registration framework which works with sitkFloat32 or sitkFloat64. To address this issue you can either specify the appropriate pixel type when reading or creating the image, or use the <a href=\"https://itk.org/SimpleITKDoxygen/html/namespaceitk_1_1simple.html#af8c9d7cc96a299a05890e9c3db911885\">Cast function</a>.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "\n",
    "import numpy as np\n",
    "import os\n",
    "from ipywidgets import interact, fixed\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from downloaddata import fetch_data as fdata\n",
    "\n",
    "OUTPUT_DIR = 'output'\n",
    "\n",
    "image_viewer = sitk.ImageViewer()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image Creation\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": {
    "simpleitk_error_allowed": "Exception thrown in SimpleITK ImageViewer_Execute:"
   },
   "outputs": [],
   "source": [
    "image_3D = sitk.Image(256, 128, 64, sitk.sitkInt16)\n",
    "image_2D = sitk.Image(64, 64, sitk.sitkFloat32)\n",
    "image_RGB = sitk.Image([128,64], sitk.sitkVectorUInt8, 3)\n",
    "\n",
    "image_viewer.Execute(image_3D)\n",
    "image_viewer.Execute(image_RGB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Or, creation from file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logo = sitk.ReadImage(fdata('SimpleITK.jpg'))\n",
    "\n",
    "# GetArrayViewFromImage returns an immutable numpy array view to the data.\n",
    "plt.imshow(sitk.GetArrayViewFromImage(logo))\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Image Attributes (Meta-Data)\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": {},
   "outputs": [],
   "source": [
    "selected_image = image_3D\n",
    "print('Before modification:')\n",
    "print('origin: ' + str(selected_image.GetOrigin()))\n",
    "print('size: ' + str(selected_image.GetSize()))\n",
    "print('spacing: ' + str(selected_image.GetSpacing()))\n",
    "print('direction: ' + str(selected_image.GetDirection()))\n",
    "print('pixel type: ' + str(selected_image.GetPixelIDTypeAsString()))\n",
    "print('number of pixel components: ' + str(selected_image.GetNumberOfComponentsPerPixel()))\n",
    "\n",
    "selected_image.SetOrigin((78.0, 76.0, 77.0))\n",
    "selected_image.SetSpacing([0.5,0.5,3.0])\n",
    "\n",
    "print('\\nAfter modification:')\n",
    "print('origin: ' + str(selected_image.GetOrigin()))\n",
    "print('spacing: ' + str(selected_image.GetSpacing()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing Pixels and Slicing\n",
    "\n",
    "Either use the ``GetPixel`` and ``SetPixel`` functions or the Pythonic slicing operator. The access functions and image slicing operator are in [x,y,z] order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Brute force sub-sampling \n",
    "logo_subsampled = logo[::2,::2]\n",
    "\n",
    "# Get the sub-image containing the word Simple\n",
    "simple = logo[0:115,:]\n",
    "\n",
    "# Get the sub-image containing the word Simple and flip it\n",
    "simple_flipped = logo[115:0:-1,:]\n",
    "\n",
    "n = 4\n",
    "\n",
    "plt.subplot(n,1,1)\n",
    "plt.imshow(sitk.GetArrayViewFromImage(logo))\n",
    "plt.axis('off');\n",
    "\n",
    "plt.subplot(n,1,2)\n",
    "plt.imshow(sitk.GetArrayViewFromImage(logo_subsampled))\n",
    "plt.axis('off');\n",
    "\n",
    "plt.subplot(n,1,3)\n",
    "plt.imshow(sitk.GetArrayViewFromImage(simple))\n",
    "plt.axis('off')\n",
    "\n",
    "plt.subplot(n,1,4)\n",
    "plt.imshow(sitk.GetArrayViewFromImage(simple_flipped))\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image operations\n",
    "\n",
    "SimpleITK supports basic arithmetic operations between images while taking into account their meta-data. Images must physically overlap (pixel by pixel).\n",
    "\n",
    "How close do physical attributes (meta-data values) need to be in order to be considered equivalent?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_width = 128\n",
    "img_height = 64\n",
    "img1 = sitk.Image((img_width, img_height), sitk.sitkUInt8)\n",
    "for i in range(img_width):\n",
    "    img1[i,1] = 5\n",
    "\n",
    "img2 = sitk.Image(img1.GetSize(), sitk.sitkUInt8)\n",
    "#img2.SetDirection([0,1,0.5,0.5])\n",
    "img2.SetOrigin([0.000001,0.000001])\n",
    "for i in range(img_width):\n",
    "    img2[i,1] = 120\n",
    "    img2[i,img_height//2] = 60\n",
    "\n",
    "# Add the two images (upper line has a value of 125, and lower of 60) \n",
    "img3 = img1 + img2\n",
    "\n",
    "plt.imshow(sitk.GetArrayViewFromImage(img3), cmap=plt.cm.Greys_r)\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparative operators (&gt;, &gt;=, &lt;, &lt;=, ==) are also supported, returning binary images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholded_image = img3>50\n",
    "plt.imshow(sitk.GetArrayViewFromImage(thresholded_image), cmap=plt.cm.Greys_r)\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SimpleITK2Numpy and Numpy2SimpleITK\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]\n",
    "\n",
    "### SimpleITK2Numpy\n",
    "\n",
    "1. ```GetArrayFromImage()```: returns a copy of the image data. You can then freely modify the data as it has no effect on the original SimpleITK image.\n",
    "2. ```GetArrayViewFromImage()```: returns a view on the image data which is useful for display in a memory efficient manner. You cannot modify the data and __the view will be invalid if the original SimpleITK image is deleted__.\n",
    "\n",
    "### Numpy2SimpleITK\n",
    "1. ```GetImageFromArray()```: returns a SimpleITK image with origin set to zero, spacing set to one for all dimensions, and the direction cosine matrix set to identity. Intensity data is copied from the numpy array. __In most cases you will need to set appropriate meta-data values.__ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nda = sitk.GetArrayFromImage(image_3D)\n",
    "print(image_3D.GetSize())\n",
    "print(nda.shape)\n",
    "\n",
    "nda = sitk.GetArrayFromImage(image_RGB)\n",
    "print(image_RGB.GetSize())\n",
    "print(nda.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gabor_image = sitk.GaborSource(size=[64,64], frequency=.03)\n",
    "# Getting a numpy array view on the image data doesn't copy the data\n",
    "nda_view = sitk.GetArrayViewFromImage(gabor_image)\n",
    "plt.imshow(nda_view, cmap=plt.cm.Greys_r)\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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",
    "img = sitk.GetImageFromArray(nda)\n",
    "print(img.GetSize())\n",
    "\n",
    "      #if this is supposed to be a 2D color image [x=20,y=10]\n",
    "img = sitk.GetImageFromArray(nda, isVector=True)\n",
    "print(img.GetSize())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading and Writing\n",
    "\n",
    "SimpleITK can read and write images stored in a single file, or a set of files (e.g. DICOM series). The toolkit provides both an object oriented and a procedural interface. The primary difference being that the object oriented approach provides more control while the procedural interface is more convenient.\n",
    "\n",
    "We look at DICOM images as an example illustrating this difference. Images stored in the DICOM format have a meta-data dictionary associated with them, which is populated with the DICOM tags. When a DICOM image series is read as a single image volume, the resulting image's meta-data dictionary is not populated since DICOM tags are specific to each of the files in the series. If you use the procedural method for reading the series then you do not have access to the set of meta-data dictionaries associated with each of the files. To obtain each dictionary you will have to access each of the files separately. On the other hand, if you use the object oriented interface, the set of dictionaries will be accessible from the ```ImageSeriesReader``` which you used to read the DICOM series. The meta-data dictionary for each file is available using the <a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ImageSeriesReader.html#a337b19b6bc101f5571455afb46514b6d\">HasMetaDataKey</a> and <a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ImageSeriesReader.html#a19995f33b86c60e2ae4878cb4d8c30ee\">GetMetaData</a> methods. \n",
    "\n",
    "We start with reading and writing an image using the procedural interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = sitk.ReadImage(fdata('SimpleITK.jpg'))\n",
    "sitk.WriteImage(img, os.path.join(OUTPUT_DIR, 'SimpleITK.png'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read an image in JPEG format and cast the pixel type according to user selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Several pixel types, some make sense in this case (vector types) and some are just show\n",
    "# that the user's choice will force the pixel type even when it doesn't make sense.\n",
    "pixel_types = { 'sitkVectorUInt8' : sitk.sitkVectorUInt8,\n",
    "                'sitkVectorUInt16' : sitk.sitkVectorUInt16,\n",
    "                'sitkVectorFloat64' : sitk.sitkVectorFloat64}\n",
    "\n",
    "def pixel_type_dropdown_callback(pixel_type, pixel_types_dict):\n",
    "    #specify the file location and the pixel type we want\n",
    "    img = sitk.ReadImage(fdata('SimpleITK.jpg'), pixel_types_dict[pixel_type])\n",
    "    \n",
    "    print(img.GetPixelIDTypeAsString())\n",
    "    print(img[0,0])\n",
    "    plt.imshow(sitk.GetArrayViewFromImage(img))\n",
    "    plt.axis('off')\n",
    " \n",
    "interact(pixel_type_dropdown_callback, pixel_type=list(pixel_types.keys()), pixel_types_dict=fixed(pixel_types));     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read a DICOM series and write it as a single mha file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_directory = os.path.dirname(fdata(\"CIRS057A_MR_CT_DICOM/readme.txt\"))\n",
    "series_ID = '1.2.840.113619.2.290.3.3233817346.783.1399004564.515'\n",
    "\n",
    "# Use the functional interface to read the image series.\n",
    "original_image = sitk.ReadImage(sitk.ImageSeriesReader_GetGDCMSeriesFileNames(data_directory,series_ID))\n",
    "\n",
    "# Write the image.\n",
    "output_file_name_3D = os.path.join(OUTPUT_DIR, '3DImage.mha')\n",
    "sitk.WriteImage(original_image, output_file_name_3D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a specific DICOM series from a directory and only then load user selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_directory = os.path.dirname(fdata(\"CIRS057A_MR_CT_DICOM/readme.txt\"))\n",
    "# Global variable 'selected_series' is updated by the interact function\n",
    "selected_series = ''\n",
    "file_reader = sitk.ImageFileReader()\n",
    "def DICOM_series_dropdown_callback(series_to_load, series_dictionary):\n",
    "    global selected_series\n",
    "               # Print some information about the series from the meta-data dictionary\n",
    "               # DICOM standard part 6, Data Dictionary: http://medical.nema.org/medical/dicom/current/output/pdf/part06.pdf\n",
    "    file_reader.SetFileName(series_dictionary[series_to_load][0])\n",
    "    file_reader.ReadImageInformation()\n",
    "    tags_to_print = {'0010|0010': 'Patient name: ', \n",
    "                     '0008|0060' : 'Modality: ',\n",
    "                     '0008|0021' : 'Series date: ',\n",
    "                     '0008|0080' : 'Institution name: ',\n",
    "                     '0008|1050' : 'Performing physician\\'s name: '}\n",
    "    for tag in tags_to_print:\n",
    "        try:\n",
    "            print(tags_to_print[tag] + file_reader.GetMetaData(tag))\n",
    "        except: # Ignore if the tag isn't in the dictionary\n",
    "            pass\n",
    "    selected_series = series_to_load                    \n",
    "\n",
    "# Directory contains multiple DICOM studies/series, store\n",
    "# in dictionary with key being the series ID\n",
    "series_file_names = {}\n",
    "series_IDs = sitk.ImageSeriesReader_GetGDCMSeriesIDs(data_directory)\n",
    "            # Check that we have at least one series\n",
    "if series_IDs:\n",
    "    for series in series_IDs:\n",
    "        series_file_names[series] = sitk.ImageSeriesReader_GetGDCMSeriesFileNames(data_directory, series)\n",
    "    \n",
    "    interact(DICOM_series_dropdown_callback, series_to_load=list(series_IDs), series_dictionary=fixed(series_file_names)); \n",
    "else:\n",
    "    print('Data directory does not contain any DICOM series.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = sitk.ReadImage(series_file_names[selected_series])\n",
    "# Display the image slice from the middle of the stack, z axis\n",
    "z = img.GetDepth()//2\n",
    "plt.imshow(sitk.GetArrayViewFromImage(img)[z,:,:], cmap=plt.cm.Greys_r)\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Write the volume as a series of JPEGs. The WriteImage function receives a volume and a list of images names and writes the volume according to the z axis. For a displayable result we need to rescale the image intensities (default is [0,255]) since the JPEG format requires a cast to the UInt8 pixel type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sitk.WriteImage(sitk.Cast(sitk.RescaleIntensity(img), sitk.sitkUInt8), \n",
    "                [os.path.join(OUTPUT_DIR, 'slice{0:03d}.jpg'.format(i)) for i in range(img.GetSize()[2])]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Advanced Reading\n",
    "\n",
    "SimpleITK relies on the registered ImageIO components to indicate whether they can read a file and then perform the reading. This is done automatically. You can select a specific ImageIO if multiple components can read the same format and you want to use a specific one.\n",
    "\n",
    "Additionally, some of the ImageIO components support streaming. That is, they can read a portion of the image into memory. This is useful when working with large images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_expected": "Error JPEGImageIO could not open file"
   },
   "outputs": [],
   "source": [
    "file_reader = sitk.ImageFileReader()\n",
    "\n",
    "# print the reader and see which ImageIO components are registered\n",
    "print('\\n',file_reader)\n",
    "\n",
    "# select the JPEG ImageIO and read\n",
    "file_reader.SetImageIO('JPEGImageIO')\n",
    "file_reader.SetFileName(fdata('SimpleITK.jpg'))\n",
    "logo = file_reader.Execute()\n",
    "\n",
    "# when the ImageIO doesn't match the image type we get an exception\n",
    "file_reader.SetFileName(fdata('training_001_ct.mha'))\n",
    "ct_head = file_reader.Execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_allowed": "Exception thrown in SimpleITK ImageViewer_Execute:"
   },
   "outputs": [],
   "source": [
    "# reset the reader's behavior so that it automatically selects the ImageIO    \n",
    "file_reader.SetImageIO('')\n",
    "ct_head = file_reader.Execute()\n",
    "image_viewer.Execute(ct_head)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If our region of interest is only in the central region of the image we can just read that region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_allowed": "Exception thrown in SimpleITK ImageViewer_Execute:"
   },
   "outputs": [],
   "source": [
    "file_reader = sitk.ImageFileReader()\n",
    "file_reader.SetFileName(fdata('training_001_ct.mha'))\n",
    "\n",
    "# read the image information without reading the bulk data, compute ROI start and size and read it.\n",
    "file_reader.ReadImageInformation()\n",
    "start_index, extract_size = zip(*[(int(0.25*sz), int(0.5*sz)) for sz in file_reader.GetSize()])\n",
    "file_reader.SetExtractIndex(start_index)\n",
    "file_reader.SetExtractSize(extract_size)\n",
    "\n",
    "image_viewer.Execute(file_reader.Execute())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resampling\n",
    "\n",
    "<img src=\"figures/resampling.svg\"/><br><br>\n",
    "\n",
    "Resampling as the verb implies is the action of sampling an image, which itself is a sampling of an original continuous signal.\n",
    "\n",
    "Generally speaking, resampling in SimpleITK involves four components:\n",
    "1. Image - the image we resample, given in coordinate system $m$.\n",
    "2. Resampling grid - a regular grid of points given in coordinate system $f$ which will be mapped to coordinate system $m$.\n",
    "2. Transformation $T_f^m$ - maps points from coordinate system $f$ to coordinate system $m$, $^mp = T_f^m(^fp)$.\n",
    "3. Interpolator - method for obtaining the intensity values at arbitrary points in coordinate system $m$ from the values of the points defined by the Image.\n",
    "\n",
    "\n",
    "While SimpleITK provides a large number of interpolation methods, the two most commonly used are ```sitkLinear``` and ```sitkNearestNeighbor```. The former is used for most interpolation tasks, a compromise between accuracy and computational efficiency. The later is used to interpolate labeled images representing a segmentation, it is the only interpolation approach which will not introduce new labels into the result.\n",
    "\n",
    "SimpleITK's procedural API provides three methods for performing resampling, with the difference being the way you specify the resampling grid:\n",
    "\n",
    "1. ```Resample(const Image &image1, Transform transform, InterpolatorEnum interpolator, double defaultPixelValue, PixelIDValueEnum outputPixelType)```\n",
    "2. ```Resample(const Image &image1, const Image &referenceImage, Transform transform, InterpolatorEnum interpolator, double defaultPixelValue, PixelIDValueEnum outputPixelType)```\n",
    "3. ```Resample(const Image &image1, std::vector< uint32_t > size, Transform transform, InterpolatorEnum interpolator, std::vector< double > outputOrigin, std::vector< double > outputSpacing, std::vector< double > outputDirection, double defaultPixelValue, PixelIDValueEnum outputPixelType)```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def resample_display(image, euler2d_transform, tx, ty, theta):\n",
    "    euler2d_transform.SetTranslation((tx, ty))\n",
    "    euler2d_transform.SetAngle(theta)\n",
    "    \n",
    "    resampled_image = sitk.Resample(image, euler2d_transform)\n",
    "    plt.imshow(sitk.GetArrayFromImage(resampled_image))\n",
    "    plt.axis('off')    \n",
    "    plt.show()\n",
    "\n",
    "euler2d = sitk.Euler2DTransform()\n",
    "# Why do we set the center?\n",
    "euler2d.SetCenter(logo.TransformContinuousIndexToPhysicalPoint(np.array(logo.GetSize())/2.0))\n",
    "interact(resample_display, image=fixed(logo), euler2d_transform=fixed(euler2d), tx=(-128.0, 128.0,2.5), ty=(-64.0, 64.0), theta=(-np.pi/4.0,np.pi/4.0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Common Errors\n",
    "\n",
    "It is not uncommon to end up with an empty (all black) image after resampling. This is due to:\n",
    "1. Using wrong settings for the resampling grid, not too common, but does happen.\n",
    "2. Using the inverse of the transformation $T_f^m$. This is a relatively common error, which is readily addressed by invoking the transformations ```GetInverse``` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the Resampling Grid\n",
    "\n",
    "In the example above we arbitrarily used the original image grid as the resampling grid. As a result, for many of the transformations the resulting image contained black pixels, pixels which were mapped outside the spatial domain of the original image and a partial view of the original image.\n",
    "\n",
    "If we want the resulting image to contain all of the original image no matter the transformation, we will need to define the resampling grid using our knowledge of the original image's spatial domain and the **inverse** of the given transformation. \n",
    "\n",
    "Computing the bounds of the resampling grid when dealing with an affine transformation is straightforward. An affine transformation preserves convexity with extreme points mapped to extreme points. Thus we only need to apply the **inverse** transformation to the corners of the original image to obtain the bounds of the resampling grid.\n",
    "\n",
    "Computing the bounds of the resampling grid when dealing with a BSplineTransform or DisplacementFieldTransform is more involved as we are not guaranteed that extreme points are mapped to extreme points. This requires that we apply the **inverse** transformation to all points in the original image to obtain the bounds of the resampling grid.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "euler2d = sitk.Euler2DTransform()\n",
    "# Why do we set the center?\n",
    "euler2d.SetCenter(logo.TransformContinuousIndexToPhysicalPoint(np.array(logo.GetSize())/2.0))\n",
    "\n",
    "tx = 64\n",
    "ty = 32\n",
    "euler2d.SetTranslation((tx, ty))\n",
    "\n",
    "extreme_points = [logo.TransformIndexToPhysicalPoint((0,0)), \n",
    "                  logo.TransformIndexToPhysicalPoint((logo.GetWidth(),0)),\n",
    "                  logo.TransformIndexToPhysicalPoint((logo.GetWidth(),logo.GetHeight())),\n",
    "                  logo.TransformIndexToPhysicalPoint((0,logo.GetHeight()))]\n",
    "inv_euler2d = euler2d.GetInverse()\n",
    "\n",
    "extreme_points_transformed = [inv_euler2d.TransformPoint(pnt) for pnt in extreme_points]\n",
    "min_x = min(extreme_points_transformed)[0]\n",
    "min_y = min(extreme_points_transformed, key=lambda p: p[1])[1]\n",
    "max_x = max(extreme_points_transformed)[0]\n",
    "max_y = max(extreme_points_transformed, key=lambda p: p[1])[1]\n",
    "\n",
    "# Use the original spacing (arbitrary decision).\n",
    "output_spacing = logo.GetSpacing()\n",
    "# Identity cosine matrix (arbitrary decision).   \n",
    "output_direction = [1.0, 0.0, 0.0, 1.0]\n",
    "# Minimal x,y coordinates are the new origin.\n",
    "output_origin = [min_x, min_y]\n",
    "# Compute grid size based on the physical size and spacing.\n",
    "output_size = [int((max_x-min_x)/output_spacing[0]), int((max_y-min_y)/output_spacing[1])]\n",
    "\n",
    "resampled_image = sitk.Resample(logo, output_size, euler2d, sitk.sitkLinear, output_origin, output_spacing, output_direction)\n",
    "plt.imshow(sitk.GetArrayViewFromImage(resampled_image))\n",
    "plt.axis('off')    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are you puzzled by the result? Is the output just a copy of the input? Add a rotation to the code above and see what happens (```euler2d.SetAngle(0.79)```)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"03_data_augmentation.ipynb\"><h2 align=right>Next &raquo;</h2></a>"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}