Python/03_Image_Details.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Images</h1>\n",
    "\n",
    "<table width=\"100%\">\n",
    "<tr style=\"background-color: red;\"><td><font color=\"white\">SimpleITK conventions:</font></td></tr>\n",
    "<tr><td>\n",
    "<ul>\n",
    "<li>Image access is in x,y,z order, image.GetPixel(x,y,z) or image[x,y,z], with zero based indexing.</li>\n",
    "<li>If the output of an ITK filter has non-zero starting index, then the index will be set to 0, and the origin adjusted accordingly.</li>\n",
    "</ul>\n",
    "</td></tr>\n",
    "</table>\n",
    "\n",
    "The unique feature of SimpleITK (derived from ITK) as a toolkit for image manipulation and analysis is that it views <b>images as physical objects occupying a bounded region in physical space</b>. In addition images can have different spacing between pixels along each axis, and the axes are not necessarily orthogonal. The following figure illustrates these concepts. \n",
    "\n",
    "<img src=\"ImageOriginAndSpacing.png\" style=\"width:700px\"/><br><br>\n",
    "\n",
    "\n",
    "### 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": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "import SimpleITK as sitk\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from ipywidgets import interact, fixed\n",
    "import os\n",
    "\n",
    "OUTPUT_DIR = 'Output'\n",
    "\n",
    "# Utility method that either downloads data from the Girder repository or\n",
    "# if already downloaded returns the file name for reading from disk (cached data).\n",
    "%run update_path_to_download_script\n",
    "from downloaddata import fetch_data as fdata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load your first image and display it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logo = sitk.ReadImage(fdata('SimpleITK.jpg'))\n",
    "\n",
    "plt.imshow(sitk.GetArrayViewFromImage(logo))\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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": {},
   "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": {},
   "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": {},
   "source": [
    "Image dimension queries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": {},
   "source": [
    "What is the depth of a 2D image?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(image_2D.GetSize())\n",
    "print(image_2D.GetDepth())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pixel/voxel type queries: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(image_3D.GetPixelIDValue())\n",
    "print(image_3D.GetPixelIDTypeAsString())\n",
    "print(image_3D.GetNumberOfComponentsPerPixel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the dimension and size of a Vector image and its data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(image_RGB.GetDimension())\n",
    "print(image_RGB.GetSize())\n",
    "print(image_RGB.GetNumberOfComponentsPerPixel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing Pixels and Slicing\n",
    "\n",
    "The Image class's member functions ``GetPixel`` and ``SetPixel`` provide an ITK-like interface for pixel access."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(image_3D.GetPixel)"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "Slicing of SimpleITK images returns a copy of the image data. \n",
    "\n",
    "This is similar to slicing Python lists and differs from the \"view\" returned by slicing numpy arrays. "
   ]
  },
  {
   "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": [
    "Draw a square on top of the logo image:\n",
    "After running this cell, uncomment \"Version 3\" and see its effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Version 0: get the numpy array and assign the value via broadcast - later on you will need to construct \n",
    "# a new image from the array \n",
    "logo_pixels = sitk.GetArrayFromImage(logo)\n",
    "logo_pixels[0:10,0:10] = [0,255,0]\n",
    "\n",
    "# Version 1: generates an error, the image slicing returns a new image and you cannot assign a value to an image \n",
    "#logo[0:10,0:10] = [255,0,0]\n",
    "\n",
    "# Version 2: image slicing returns a new image, so all assignments here will not have any effect on the original\n",
    "# 'logo' image\n",
    "logo_subimage = logo[0:10, 0:10]\n",
    "for x in range(0,10):\n",
    "    for y in range(0,10):\n",
    "        logo_subimage[x,y] = [255,0,0]\n",
    "\n",
    "# Version 3: modify the original image, iterate and assign a value to each pixel\n",
    "#for x in range(0,10):\n",
    "#    for y in range(0,10):\n",
    "#        logo[x,y] = [255,0,0]\n",
    "\n",
    "        \n",
    "plt.subplot(2,1,1)\n",
    "plt.imshow(sitk.GetArrayViewFromImage(logo))\n",
    "plt.axis('off')\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.imshow(logo_pixels)\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "source": [
    "### From SimpleITK to numpy\n",
    "\n",
    "We have two options for converting from SimpleITK to numpy:\n",
    "* 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",
    "* 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__."
   ]
  },
  {
   "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": {
    "simpleitk_error_expected": "assignment destination is read-only"
   },
   "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');\n",
    "\n",
    "# Trying to assign a value to the array view will throw an exception\n",
    "nda_view[0,0] = 255"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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": [
    "### There and back again\n",
    "\n",
    "The following code cell illustrates a situation where your code is a combination of SimpleITK methods and custom Python code which works with intensity values or labels outside of SimpleITK. This is a reasonable approach when you implement an algorithm in Python and don't care about the physical spacing of things (you are actually assuming the volume is isotropic). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_algorithm(image_as_numpy_array):\n",
    "    # res is the image result of your algorithm, has the same grid size as the original image\n",
    "    res = image_as_numpy_array\n",
    "    return res\n",
    "\n",
    "# Starting with SimpleITK\n",
    "img = sitk.ReadImage(fdata('training_001_mr_T1.mha'))\n",
    "\n",
    "# Custom Python code working on a numpy array.\n",
    "npa_res = my_algorithm(sitk.GetArrayFromImage(img))\n",
    "\n",
    "# Converting back to SimpleITK (assumes we didn't move the image in space as we copy the information from the original)\n",
    "res_img = sitk.GetImageFromArray(npa_res)\n",
    "res_img.CopyInformation(img)\n",
    "\n",
    "# Continuing to work with SimpleITK images\n",
    "res_img - img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Image operations\n",
    "\n",
    "SimpleITK supports basic arithmetic operations between images, <b>taking into account their physical space</b>.\n",
    "\n",
    "Repeatedly run this cell. Fix the error (comment out the SetDirection, then SetSpacing). Why doesn't the SetOrigin line cause a problem? How close do two physical attributes need to be in order to be considered equivalent?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_expected": "Inputs do not occupy the same physical space!"
   },
   "outputs": [],
   "source": [
    "img1 = sitk.Image(24,24, sitk.sitkUInt8)\n",
    "img1[0,0] = 0\n",
    "\n",
    "img2 = sitk.Image(img1.GetSize(), sitk.sitkUInt8)\n",
    "img2.SetDirection([0,1,0.5,0.5])\n",
    "img2.SetSpacing([0.5,0.8])\n",
    "img2.SetOrigin([0.000001,0.000001])\n",
    "img2[0,0] = 255\n",
    "\n",
    "img3 = img1 + img2\n",
    "print(img3[0,0])"
   ]
  },
  {
   "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).\n",
    "\n",
    "Images stored in the DICOM format have a meta-data dictionary associated with them, which is populated with the DICOM tags. When a DICOM series is read as a single image, the meta-data information is not available since DICOM tags are specific to each file. If you need the meta-data, you have three options:\n",
    "\n",
    "1. Using the object oriented interface's [ImageSeriesReader](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ImageSeriesReader.html) class, configure it to load the tags using the `MetaDataDictionaryArrayUpdateOn` method and possibly the `LoadPrivateTagsOn` method if you need the private tags. Once the series is read you can access the meta-data from the series reader using the `GetMetaDataKeys`, `HasMetaDataKey`, and `GetMetaData`.\n",
    "\n",
    "2. Using the object oriented interface's [ImageFileReader](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ImageFileReader.html), set a specific slice's file name and only read it's meta-data using the `ReadImageInformation` method which only reads the meta-data but not the bulk pixel information. Once the meta-data is read you can access it from the file reader using the `GetMetaDataKeys`, `HasMetaDataKey`, and `GetMetaData`.\n",
    "\n",
    "3. Using the object oriented interface's [ImageFileReader](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ImageFileReader.html), set a specific slice's file name and read it. Or using the procedural interface's, [ReadImage](https://itk.org/SimpleITKDoxygen/html/namespaceitk_1_1simple.html#ae3b678b5b043c5a8c93aa616d5ee574c) function, read a specific file. You can then access the meta-data directly from the [Image](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1Image.html) using the `GetMetaDataKeys`, `HasMetaDataKey`, and `GetMetaData`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following cell, we read an image in JPEG format, and write it as PNG and BMP. File formats are deduced from the file extension. Appropriate pixel type is also set - you can override this and force a pixel type of your choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = sitk.ReadImage(fdata('SimpleITK.jpg'))\n",
    "print(img.GetPixelIDTypeAsString())\n",
    "\n",
    "# write as PNG and BMP\n",
    "sitk.WriteImage(img, os.path.join(OUTPUT_DIR, 'SimpleITK.png'))\n",
    "sitk.WriteImage(img, os.path.join(OUTPUT_DIR, 'SimpleITK.bmp'))"
   ]
  },
  {
   "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": {
    "simpleitk_error_allowed": "Image RGB array must be uint8 or floating point; found uint16"
   },
   "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",
    "# (e.g. sitkVectorUInt16 or sitkUInt8).\n",
    "pixel_types = { 'sitkUInt8': sitk.sitkUInt8,\n",
    "                'sitkUInt16' : sitk.sitkUInt16,\n",
    "                'sitkFloat64' : sitk.sitkFloat64,\n",
    "                '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",
    "# Get the list of files belonging to a specific series ID.\n",
    "reader = sitk.ImageSeriesReader()\n",
    "# Use the functional interface to read the image series.\n",
    "original_image = sitk.ReadImage(reader.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)\n",
    "\n",
    "# Read it back again.\n",
    "written_image = sitk.ReadImage(output_file_name_3D)\n",
    "\n",
    "# Check that the original and written image are the same.\n",
    "statistics_image_filter = sitk.StatisticsImageFilter()\n",
    "statistics_image_filter.Execute(original_image - written_image)\n",
    "\n",
    "# Check that the original and written files are the same\n",
    "print('Max, Min differences are : {0}, {1}'.format(statistics_image_filter.GetMaximum(), statistics_image_filter.GetMinimum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write an image series as JPEG. 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(written_image), sitk.sitkUInt8), \n",
    "                [os.path.join(OUTPUT_DIR, 'slice{0:03d}.jpg'.format(i)) for i in range(written_image.GetSize()[2])]) "
   ]
  },
  {
   "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",
    "reader = sitk.ImageSeriesReader()\n",
    "series_file_names = {}\n",
    "series_IDs = reader.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] = reader.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": [
    "reader.SetFileNames(series_file_names[selected_series])\n",
    "img = reader.Execute()\n",
    "# Display the image slice from the middle of the stack, z axis\n",
    "z = int(img.GetDepth()/2)\n",
    "plt.imshow(sitk.GetArrayViewFromImage(img)[z,:,:], cmap=plt.cm.Greys_r)\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-channel images and color\n",
    "\n",
    "Generally speaking, SimpleITK represents color images as multi-channel images independent of a [color space](https://en.wikipedia.org/wiki/Color_space). It is up to you to interpret the channels correctly based on additional color space knowledge prior to using them for display or any other purpose.\n",
    "\n",
    "The following cells illustrate reading and interpretation of an interesting image in DICOM format. It is a photograph of an X-ray on a light box (yes, there are some strange things in the wild). The meta-data dictionary for this image contains the relevant information that will allow us to intelligently interpret it.\n",
    "\n",
    "We will look at the image's modality (0008|0060), which in our case is XC, which stands for \"External Camera Photography\". After reading the image we see that it has three channels and try to display it. At this point we realize that something is wrong, and we check the image's Photometric Interpretation (0028|0004), color space in DICOM speak. We then modify the pixel values so that we can display them using our image viewer's expected color space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_allowed": "Exception thrown in SimpleITK Show:"
   },
   "outputs": [],
   "source": [
    "xray_photo = sitk.ReadImage(fdata('photo.dcm'))\n",
    "print('Image Modality: {0}'.format(xray_photo.GetMetaData('0008|0060')))\n",
    "print('Number of channels: {0}'.format(xray_photo.GetNumberOfComponentsPerPixel()))\n",
    "\n",
    "# Display the image using Fiji which expects the channels to be in the RGB color space\n",
    "sitk.Show(xray_photo)\n",
    "\n",
    "# In what color space are the pixels?\n",
    "print('Photomertic Interpretation: {0}'.format(xray_photo.GetMetaData('0028|0004')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_allowed": "Exception thrown in SimpleITK Show:"
   },
   "outputs": [],
   "source": [
    "# Our pixels are in the YCbCr color space, and we want them in RGB, so we need to change the representation\n",
    "def full_ycbcr_to_rgb(image):\n",
    "    channels = [sitk.VectorIndexSelectionCast(image,i, sitk.sitkFloat32) for i in range(image.GetNumberOfComponentsPerPixel())]\n",
    "    chan0 = channels[0]\n",
    "    chan1 = channels[1] - 128\n",
    "    chan2 = channels[2] - 128\n",
    "    return sitk.Compose([sitk.Clamp(chan0 + 1.402*chan2, sitk.sitkUInt8),\n",
    "                         sitk.Clamp(chan0 - 0.114 * 1.772 / 0.587*chan1 - 0.299 * 1.402 / 0.587*chan2, sitk.sitkUInt8),\n",
    "                         sitk.Clamp(chan0 + 1.772*chan1, sitk.sitkUInt8)])\n",
    "\n",
    "sitk.Show(full_ycbcr_to_rgb(xray_photo))\n",
    "\n",
    "# x-rays are expected to be a single channel gray scale image and not a color image. To obtain a gray scale image\n",
    "# corresponding to the original three channel image we only need to take the luminance (Y channel).\n",
    "sitk.Show(sitk.VectorIndexSelectionCast(xray_photo,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finer control\n",
    "\n",
    "The ImageFileReader's interface provides finer control for reading, allowing us to require the use of a specific IO and allowing us to stream parts of an image to memory without reading the whole image (supported by a subset of the ImageIO components).\n",
    "\n",
    "#### Selecting a Specific Image IO\n",
    "SimpleITK relies on the registered ImageIOs to indicate whether they can read a file and then perform the reading. This is done automatically, going over the set of ImageIOs and inquiring whether they can read the given file. The first one that can is selected. If multiple ImageIOs can read a specific format, we do not know which one was used for the task (e.g. TIFFImageIO and LSMImageIO, which is derived from it, can both read tif files). In some cases you may want to use a specific IO, possibly one that reads the file faster, or supports a more complete feature set associated with the file format.\n",
    "\n",
    "The next cell shows how to find out which ImageIOs are registered and specify the one we want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_reader = sitk.ImageFileReader()\n",
    "\n",
    "# Get a tuple listing all registered ImageIOs\n",
    "image_ios_tuple = file_reader.GetRegisteredImageIOs()\n",
    "print(\"The supported image IOs are: \" + str(image_ios_tuple))\n",
    "\n",
    "# Optionally, just print the reader and see which ImageIOs are registered\n",
    "print('\\n',file_reader)\n",
    "\n",
    "# Specify the JPEGImageIO and read file\n",
    "file_reader.SetImageIO('JPEGImageIO')\n",
    "file_reader.SetFileName(fdata('SimpleITK.jpg'))\n",
    "logo = file_reader.Execute()\n",
    "\n",
    "# Unfortunately, now reading a non JPEG image will fail\n",
    "try:\n",
    "    file_reader.SetFileName(fdata('cthead1.png'))\n",
    "    ct_head = file_reader.Execute()\n",
    "except RuntimeError:\n",
    "    print('Got a RuntimeError exception.')\n",
    "    \n",
    "# We can reset the file reader to its default behaviour so that it automatically \n",
    "# selects the ImageIO\n",
    "file_reader.SetImageIO('')\n",
    "ct_head = file_reader.Execute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Streaming Image IO\n",
    "\n",
    "Some of the ImageIOs supported in SimpleITK allow you to stream in sub-regions of an image without the need to read the whole image into memory. This is very useful when you are memory constrained (either your images are large or your memory is limited). \n",
    "\n",
    "The ImageIOs that support streaming include HDF5ImageIO, VTKImageIO, NiftiImageIO, MetaImageIO...\n",
    "\n",
    "The next cell shows how to read in a sub/cropped image from a larger image. We read the central 1/3 portion of the image [1/3,2/3] of the original image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_allowed": "Exception thrown in SimpleITK Show:"
   },
   "outputs": [],
   "source": [
    "file_reader = sitk.ImageFileReader()\n",
    "file_reader.SetFileName(fdata('vm_head_rgb.mha'))\n",
    "\n",
    "file_reader.ReadImageInformation()\n",
    "image_size = file_reader.GetSize()\n",
    "start_index, extract_size = zip(*[(int(1.0/3.0*sz), int(1.0/3.0*sz)) for sz in file_reader.GetSize()])\n",
    "\n",
    "file_reader.SetExtractIndex(start_index)\n",
    "file_reader.SetExtractSize(extract_size)\n",
    "\n",
    "sitk.Show(file_reader.Execute())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cells show how to subtract two large images from each other with a smaller memory footprint than the direct approach, though the code is much more complex and slower than the direct approach:\n",
    "\n",
    "```\n",
    "sitk.ReadImage(image1_file_name) - sitk.ReadImage(image2_file_name)\n",
    "```\n",
    "\n",
    "**Note**: The code assume that the two images occupy the same spatial region (origin, spacing, direction cosine matrix).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def streaming_subtract(image1_file_name, image2_file_name, parts):\n",
    "    '''\n",
    "    Subtract image1 from image2 using 'parts' number of sub-regions.\n",
    "    '''\n",
    "    file_reader = sitk.ImageFileReader()\n",
    "    file_reader.SetFileName(image1_file_name)\n",
    "    file_reader.ReadImageInformation()\n",
    "    image_size = file_reader.GetSize()\n",
    "\n",
    "    # Create the result image, initially empty\n",
    "    result_img = sitk.Image(file_reader.GetSize(), file_reader.GetPixelID(), file_reader.GetNumberOfComponents())\n",
    "    result_img.SetSpacing(file_reader.GetSpacing())\n",
    "    result_img.SetOrigin(file_reader.GetOrigin())\n",
    "    result_img.SetDirection(file_reader.GetDirection())\n",
    "\n",
    "    extract_size = list(file_reader.GetSize())\n",
    "    extract_size[-1] = extract_size[-1]//parts\n",
    "    current_index = [0]*file_reader.GetDimension()\n",
    "    for i in range(parts):\n",
    "        if i == (parts-1): # last region may be smaller than the standard extract region\n",
    "            extract_size[-1] = image_size[-1] - current_index[-1]\n",
    "        file_reader.SetFileName(image1_file_name)\n",
    "        file_reader.SetExtractIndex(current_index)\n",
    "        file_reader.SetExtractSize(extract_size)\n",
    "        sub_image1 = file_reader.Execute()\n",
    "\n",
    "        file_reader.SetFileName(image2_file_name)\n",
    "        file_reader.SetExtractIndex(current_index)\n",
    "        file_reader.SetExtractSize(extract_size)\n",
    "        sub_image2 = file_reader.Execute()\n",
    "        # Paste the result of subtracting the two subregions into their location in the result_img\n",
    "        result_img = sitk.Paste(result_img, sub_image1 - sub_image2, extract_size, [0]*file_reader.GetDimension(), current_index)\n",
    "        current_index[-1] += extract_size[-1]        \n",
    "    return result_img\n",
    "\n",
    "# If you have the patience and RAM you can try this with the vm_head_rgb.mha image.\n",
    "image1_file_name = fdata('fib_sem_bacillus_subtilis.mha')\n",
    "image2_file_name = fdata('fib_sem_bacillus_subtilis.mha')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple way of seeing your system's memory usage is to open the appropriate monitoring program: (Windows) Resource Monitor; (Linux) top; (OS X) Activity Monitor. This will give you a rough idea of the memory used by the streaming vs. non streaming approaches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_img = streaming_subtract(image1_file_name, image2_file_name, parts=5)\n",
    "del result_img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_img = sitk.ReadImage(image1_file_name) - sitk.ReadImage(image2_file_name)\n",
    "del result_img"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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": 1
}