04_basic_registration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Basic Registration</h1>\n",
    "\n",
    "\n",
    "**Summary:**\n",
    "\n",
    "1. Creating an instance of the registration framework requires selection of the following components:\n",
    "   * Optimizer.\n",
    "   * Similarity metric.\n",
    "   * Interpolator.\n",
    "2. The registration framework only supports images with sitkFloat32 and sitkFloat64 pixel types (use the SimpleITK <a href=\"http://www.itk.org/SimpleITKDoxygen/html/namespaceitk_1_1simple.html#af8c9d7cc96a299a05890e9c3db911885\">Cast()</a> function if your image's pixel type is something else).\n",
    "\n",
    "3. Successful registration is highly dependent on initialization. In general you can:\n",
    "   * Use auxiliary information or user interaction to obtain an initial transformation (avoid resampling).\n",
    "   * Center the images using the <a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1CenteredTransformInitializerFilter.html\">CenteredTransformInitializer</a>.\n",
    "   * Coarsely sample the parameter space using the <a href=\"https://itk.org/Doxygen/html/classitk_1_1ExhaustiveOptimizerv4.html\">Exhaustive Optimizer</a> to obtain one or more initial transformation estimates.\n",
    "   * Manually initialize, via direct manipulation of transformation parameters and visualization or localization of corresponding points in the two images and then use the <a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1LandmarkBasedTransformInitializerFilter.html\">LandmarkBasedTransformInitializer</a>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Registration Components \n",
    "\n",
    "<img src=\"figures/ITKv4RegistrationComponentsDiagram.svg\" style=\"width:700px\"/><br><br>\n",
    "\n",
    "There are many options for creating an instance of the registration framework, all of which are configured in SimpleITK via methods of the <a href=\"http://www.itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ImageRegistrationMethod.html\">ImageRegistrationMethod</a> class. This class encapsulates many of the components available in ITK for constructing a registration instance.\n",
    "\n",
    "Currently, the available choices from the following groups of ITK components are:\n",
    "\n",
    "### Optimizers\n",
    "\n",
    "The SimpleITK registration framework supports several optimizer types via the SetOptimizerAsX() methods, these include:\n",
    "\n",
    "<ul>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ExhaustiveOptimizerv4.html\">Exhaustive</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1AmoebaOptimizerv4.html\">Nelder-Mead downhill simplex</a>, a.k.a. Amoeba.\n",
    "  </li>\n",
    "  <li>\n",
    "   <a href=\"https://itk.org/Doxygen/html/classitk_1_1PowellOptimizerv4.html\">Powell optimizer</a>.\n",
    "  </li>\n",
    "  <li>\n",
    "   <a href=\"https://itk.org/Doxygen/html/classitk_1_1OnePlusOneEvolutionaryOptimizerv4.html\">1+1 evolutionary optimizer</a>.\n",
    "  </li>\n",
    "  <li>\n",
    "  Variations on gradient descent:\n",
    "  <ul>\n",
    "    <li>\n",
    "    <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1GradientDescentOptimizerv4Template.html\">GradientDescent</a>\n",
    "    </li>\n",
    "    <li>\n",
    "    <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1GradientDescentLineSearchOptimizerv4Template.html\">GradientDescentLineSearch</a>\n",
    "    </li>\n",
    "    <li>\n",
    "    <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1RegularStepGradientDescentOptimizerv4.html\">RegularStepGradientDescent</a>\n",
    "    </li>\n",
    "  </ul>\n",
    "  </li>\n",
    "  <li>\n",
    "    <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ConjugateGradientLineSearchOptimizerv4Template.html\">ConjugateGradientLineSearch</a> \n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"https://itk.org/Insight/Doxygen/html/classitk_1_1LBFGS2Optimizerv4.html\">L-BFGS2</a> (Limited memory Broyden,  Fletcher,Goldfarb,Shannon)\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1LBFGSBOptimizerv4.html\">L-BFGS-B</a> (Limited memory Broyden,  Fletcher,Goldfarb,Shannon-Bound Constrained) - supports the use of simple constraints ($l\\leq x \\leq u$)  \n",
    "  </li>\n",
    "    \n",
    "</ul>\n",
    "\n",
    " \n",
    "### Similarity metrics\n",
    "\n",
    "The SimpleITK registration framework supports several metric types via the SetMetricAsX() methods, these include:\n",
    "\n",
    "<ul>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1MeanSquaresImageToImageMetricv4.html\">MeanSquares</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1DemonsImageToImageMetricv4.html\">Demons</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1CorrelationImageToImageMetricv4.html\">Correlation</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ANTSNeighborhoodCorrelationImageToImageMetricv4.html\">ANTSNeighborhoodCorrelation</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1JointHistogramMutualInformationImageToImageMetricv4.html\">JointHistogramMutualInformation</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1MattesMutualInformationImageToImageMetricv4.html\">MattesMutualInformation</a>\n",
    "  </li>\n",
    "</ul>\n",
    "\n",
    "\n",
    "### Interpolators\n",
    "\n",
    "The SimpleITK registration framework supports several interpolators via the SetInterpolator() method, which receives one of\n",
    "the <a href=\"http://www.itk.org/SimpleITKDoxygen/html/namespaceitk_1_1simple.html#a7cb1ef8bd02c669c02ea2f9f5aa374e5\">following enumerations</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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "from downloaddata import fetch_data as fdata\n",
    "%matplotlib notebook\n",
    "import gui\n",
    "import registration_gui as rgui\n",
    "\n",
    "import numpy as np\n",
    "import os\n",
    "OUTPUT_DIR = 'output'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read images\n",
    "\n",
    "We first read the images, specifying the pixel type that is required for registration (Float32 or Float64) and look at them. In this notebook we use a CT and MR image from the same patient. These are part of the training data from the Retrospective Image Registration Evaluation (<a href=\"http://www.insight-journal.org/rire/\">RIRE</a>) project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fixed_image =  sitk.ReadImage(fdata(\"training_001_ct.mha\"), sitk.sitkFloat32)\n",
    "moving_image = sitk.ReadImage(fdata(\"training_001_mr_T1.mha\"), sitk.sitkFloat32)\n",
    "\n",
    "ct_window_level = [835,162]\n",
    "mr_window_level = [1036,520]\n",
    "\n",
    "gui.MultiImageDisplay(image_list = [fixed_image, moving_image],                   \n",
    "                      title_list = ['fixed', 'moving'], figure_size=(8,4), window_level_list=[ct_window_level, mr_window_level]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classic Registration\n",
    "\n",
    "Estimate a 3D rigid transformation between images of different modalities. \n",
    "\n",
    "We have made the following choices with respect to initialization and registration component settings:\n",
    "\n",
    "<ul>\n",
    "\n",
    "<li>Similarity metric, mutual information (Mattes MI):\n",
    "<ul>\n",
    "  <li>Number of histogram bins, 50.</li>\n",
    "  <li>Sampling strategy, random.</li>\n",
    "  <li>Sampling percentage, 1%.</li>\n",
    "</ul>\n",
    "</li>\n",
    "<li>Interpolator, sitkLinear.</li>\n",
    "<li>Optimizer, gradient descent: \n",
    "<ul>\n",
    "  <li>Learning rate, step size along traversal direction in parameter space, 1.0 .</li>\n",
    "  <li>Number of iterations, maximal number of iterations, 100.</li>\n",
    "  <li>Convergence minimum value, value used for convergence checking in conjunction with the energy profile of the similarity metric that is estimated in the given window size, 1e-6.</li>\n",
    "  <li>Convergence window size, number of values of the similarity metric which are used to estimate the energy profile of the similarity metric, 10.</li>\n",
    "</ul>\n",
    "</li>\n",
    "</ul>\n",
    "\n",
    "We initialize registration by aligning the centers of the two volumes. To qualitatively evaluate the result we use a linked cursor approach, click on one image and the corresponding point is added to the other image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_transform = sitk.CenteredTransformInitializer(fixed_image, \n",
    "                                                      moving_image, \n",
    "                                                      sitk.Euler3DTransform(), \n",
    "                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)\n",
    "\n",
    "gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=initial_transform, fixed_window_level=ct_window_level, moving_window_level=mr_window_level);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the next cell three times:\n",
    "1. As is.\n",
    "2. Uncomment the automated optimizer scale setting so that a change in rotation (radians) has a similar effect to a change in translation (mm).\n",
    "3. Uncomment the multi-resolution settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "registration_method = sitk.ImageRegistrationMethod()\n",
    "\n",
    "# Similarity metric settings.\n",
    "registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)\n",
    "registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)\n",
    "registration_method.SetMetricSamplingPercentage(0.01)\n",
    "\n",
    "registration_method.SetInterpolator(sitk.sitkLinear)\n",
    "\n",
    "# Optimizer settings.\n",
    "registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6, convergenceWindowSize=10)\n",
    "#registration_method.SetOptimizerScalesFromPhysicalShift()\n",
    "\n",
    "# Setup for the multi-resolution framework.            \n",
    "#registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])\n",
    "#registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])\n",
    "#registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()\n",
    "\n",
    "# Don't optimize in-place, we would possibly like to run this cell multiple times.\n",
    "registration_method.SetInitialTransform(initial_transform, inPlace=False)\n",
    "\n",
    "# Connect all of the observers so that we can perform plotting during registration.\n",
    "registration_method.AddCommand(sitk.sitkStartEvent, rgui.start_plot)\n",
    "registration_method.AddCommand(sitk.sitkEndEvent, rgui.end_plot)\n",
    "registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, rgui.update_multires_iterations) \n",
    "registration_method.AddCommand(sitk.sitkIterationEvent, lambda: rgui.plot_values(registration_method))\n",
    "\n",
    "final_transform = registration_method.Execute(fixed_image, moving_image)\n",
    "\n",
    "# Always check the reason optimization terminated.\n",
    "print('Final metric value: {0}'.format(registration_method.GetMetricValue()))\n",
    "print('Optimizer\\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Qualitatively evaluate the result using a linked cursor approach (visual evaluation):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=final_transform,fixed_window_level=ct_window_level, moving_window_level=mr_window_level);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are satisfied with the results, save them to file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())\n",
    "sitk.WriteImage(moving_resampled, os.path.join(OUTPUT_DIR, 'RIRE_training_001_mr_T1_resampled.mha'))\n",
    "sitk.WriteTransform(final_transform, os.path.join(OUTPUT_DIR, 'RIRE_training_001_CT_2_mr_T1.tfm'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ITKv4 Coordinate Systems\n",
    "\n",
    "Unlike the classical registration approach where the fixed and moving images are treated differently, the ITKv4 registration framework allows you to treat both images in the same manner. This is achieved by introducing a third coordinate system, the virtual image domain.\n",
    "\n",
    "<img src=\"figures/registrationFrameworkTransformations.svg\"/><br><br>\n",
    "\n",
    "Thus, the ITK v4 registration framework deals with three transformations:\n",
    "<ul>\n",
    "<li>\n",
    "SetInitialTransform, $T_{opt}$ - composed with the moving initial transform, maps points from the virtual image domain to the moving image domain, modified during optimization. \n",
    "</li>\n",
    "<li>\n",
    "SetFixedInitialTransform $T_f$- maps points from the virtual image domain to the fixed image domain, never modified.\n",
    "</li>\n",
    "<li>\n",
    "SetMovingInitialTransform $T_m$- maps points from the virtual image domain to the moving image domain, never modified.\n",
    "</li>\n",
    "</ul>\n",
    "\n",
    "The transformation that maps points from the fixed to moving image domains is thus: $^M\\mathbf{p}  = T_{opt}(T_m(T_f^{-1}(^F\\mathbf{p})))$\n",
    "\n",
    "We now modify the previous example to use $T_{opt}$ and $T_m$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "registration_method = sitk.ImageRegistrationMethod()\n",
    "\n",
    "# Similarity metric settings.\n",
    "registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)\n",
    "registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)\n",
    "registration_method.SetMetricSamplingPercentage(0.01)\n",
    "\n",
    "registration_method.SetInterpolator(sitk.sitkLinear)\n",
    "\n",
    "# Optimizer settings.\n",
    "registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6, convergenceWindowSize=10)\n",
    "registration_method.SetOptimizerScalesFromPhysicalShift()\n",
    "\n",
    "# Setup for the multi-resolution framework.            \n",
    "registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])\n",
    "registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])\n",
    "registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()\n",
    "\n",
    "# Set the initial moving and optimized transforms.\n",
    "optimized_transform = sitk.Euler3DTransform()    \n",
    "registration_method.SetMovingInitialTransform(initial_transform)\n",
    "registration_method.SetInitialTransform(optimized_transform, inPlace=False)\n",
    "\n",
    "# Connect all of the observers so that we can perform plotting during registration.\n",
    "registration_method.AddCommand(sitk.sitkStartEvent, rgui.start_plot)\n",
    "registration_method.AddCommand(sitk.sitkEndEvent, rgui.end_plot)\n",
    "registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, rgui.update_multires_iterations) \n",
    "registration_method.AddCommand(sitk.sitkIterationEvent, lambda: rgui.plot_values(registration_method))\n",
    "\n",
    "# Need to compose the transformations after registration.\n",
    "final_transform_v4 = registration_method.Execute(fixed_image, moving_image)\n",
    "#final_transform_v4.AddTransform(initial_transform)\n",
    "\n",
    "# Always check the reason optimization terminated.\n",
    "print('Final metric value: {0}'.format(registration_method.GetMetricValue()))\n",
    "print('Optimizer\\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=final_transform_v4, fixed_window_level=ct_window_level, moving_window_level=mr_window_level);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Initialization\n",
    "\n",
    "Initialization effects both the runtime and convergence to the correct minimum. Ideally our transformation is initialized close to the correct solution ensuring convergence in a timely manner. Problem specific initialization will often yield better results than the more generic solutions we show below. As a rule of thumb, use as much prior information (external to the image content) as you can to initialize your registration task.\n",
    "\n",
    "Common initializations in the generic setting:\n",
    "1. Do nothing (a.k.a. hope/unique setting) - initialize using the identity transformation.\n",
    "2. CenteredTransformInitializer (GEOMETRY or MOMENTS) - translation based initialization, align the centers of the images or their centers of mass (intensity based).\n",
    "3. Use the exhaustive optimizer as a first step - never underestimate brute force.\n",
    "4. Manual initialization - allow an operator to control parameter settings using a GUI with visual feedback or identify multiple corresponding points in the two images. \n",
    "\n",
    "\n",
    "We start by loading our data, CT and MR scans of the CIRS (Norfolk, VA, USA) abdominal phantom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_directory = os.path.dirname(fdata(\"CIRS057A_MR_CT_DICOM/readme.txt\"))\n",
    "\n",
    "ct_window_level = [1727,-320]\n",
    "mr_window_level = [355,178]\n",
    "\n",
    "fixed_series_ID = \"1.2.840.113619.2.290.3.3233817346.783.1399004564.515\"\n",
    "moving_series_ID = \"1.3.12.2.1107.5.2.18.41548.30000014030519285935000000933\"\n",
    "\n",
    "fixed_series_filenames = sitk.ImageSeriesReader_GetGDCMSeriesFileNames(data_directory, fixed_series_ID)\n",
    "moving_series_filenames = sitk.ImageSeriesReader_GetGDCMSeriesFileNames(data_directory, moving_series_ID)\n",
    "\n",
    "fixed_image = sitk.ReadImage(fixed_series_filenames, sitk.sitkFloat32)\n",
    "moving_image = sitk.ReadImage(moving_series_filenames, sitk.sitkFloat32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identity transform initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_transform = sitk.Transform()\n",
    "gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=initial_transform, fixed_window_level=ct_window_level, moving_window_level=mr_window_level);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When working with clinical images, the DICOM tags define the orientation and position of the anatomy in the volume. The tags of interest are:\n",
    "<ul>\n",
    "  <li> (0020|0032) Image Position (Patient) : coordinates of the the first transmitted voxel. </li>\n",
    "  <li>(0020|0037) Image Orientation (Patient): directions of first row and column in 3D space. </li>\n",
    "  <li>(0018|5100) Patient Position: Patient placement on the table \n",
    "  <ul>\n",
    "  <li> Head First Prone (HFP)</li>\n",
    "  <li> Head First Supine (HFS)</li>\n",
    "  <li> Head First Decubitus Right (HFDR)</li>\n",
    "  <li> Head First Decubitus Left (HFDL)</li>\n",
    "  <li> Feet First Prone (FFP)</li>\n",
    "  <li> Feet First Supine (FFS)</li>\n",
    "  <li> Feet First Decubitus Right (FFDR)</li>\n",
    "  <li> Feet First Decubitus Left (FFDL)</li>\n",
    "  </ul>\n",
    "  </li>\n",
    "</ul>\n",
    "\n",
    "SimpleITK/ITK takes this information into account when loading DICOM images. \n",
    "\n",
    "But we are working with DICOM images, so why aren't the images oriented correctly using the identity transformation?\n",
    "\n",
    "Well, the patient position in the scanner is manually entered by the technician meaning that errors may occur, though rarely. For our data, a phantom, it is unclear which side is the \"head\" and which is the \"feet\" so the technicians entered reasonable values for each scan. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reader = sitk.ImageFileReader()\n",
    "reader.SetFileName(fixed_series_filenames[0])\n",
    "reader.ReadImageInformation()\n",
    "print('Patient name: ' + reader.GetMetaData('0010|0010') + ', Patient position:' + reader.GetMetaData('0018|5100'))\n",
    "reader.SetFileName(moving_series_filenames[0])\n",
    "reader.ReadImageInformation()\n",
    "print('Patient name: ' + reader.GetMetaData('0010|0010') + ', Patient position:' + reader.GetMetaData('0018|5100'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CenteredTransformInitializer initialization \n",
    "Compare GEOMETRY and MOMENTS based approaches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_transform = sitk.CenteredTransformInitializer(fixed_image, \n",
    "                                                      moving_image, \n",
    "                                                      sitk.Euler3DTransform(), \n",
    "                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)\n",
    "gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=initial_transform, fixed_window_level=ct_window_level, moving_window_level=mr_window_level);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exhaustive optimizer initialization\n",
    "\n",
    "The following initialization approach is a combination of using prior knowledge and the exhaustive optimizer. We know that the scans are acquired with the \"patient\" either supine (on their back) or prone (on their stomach) and that the scan direction (head-to-feet or feet-to-head) is along the images' z axis. \n",
    "We use the CenteredTransformInitializer to initialize the translation and the exhaustive optimizer to obtain an initial rigid transformation.\n",
    "\n",
    "The exhaustive optimizer evaluates the similarity metric on a grid in parameter space centered on the parameters of the initial transform. This grid is defined using three elements:\n",
    "1. numberOfSteps.\n",
    "2. stepLength.\n",
    "3. optimizer scales.\n",
    "\n",
    "The similarity metric is evaluated on the resulting parameter grid:\n",
    "initial_parameters &plusmn; numberOfSteps &times; stepLength &times; optimizerScales\n",
    "\n",
    "***Example***:\n",
    "1. numberOfSteps=[1,0,2,0,0,0]\n",
    "2. stepLength = np.pi\n",
    "3. optimizer scales = [1,1,0.5,1,1,1]\n",
    "\n",
    "Will perform 15 metric evaluations ($\\displaystyle\\prod_i (2*numberOfSteps[i] + 1)$).\n",
    "\n",
    "The parameter values for the second parameter and the last three parameters are the initial parameter values. The parameter values for the first parameter are $v_{init}-\\pi, v_{init}, v_{init}+\\pi$ and the parameter values for the third parameter are $v_{init}-\\pi, v_{init}-\\pi/2, v_{init}, v_{init}+\\pi/2, v_{init}+\\pi$.\n",
    "\n",
    "The transformation corresponding to the lowest similarity metric is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_transform = sitk.CenteredTransformInitializer(fixed_image, \n",
    "                                                      moving_image, \n",
    "                                                      sitk.Euler3DTransform(), \n",
    "                                                      sitk.CenteredTransformInitializerFilter.MOMENTS)\n",
    "registration_method = sitk.ImageRegistrationMethod()\n",
    "registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)\n",
    "registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)\n",
    "registration_method.SetMetricSamplingPercentage(0.01)\n",
    "registration_method.SetInterpolator(sitk.sitkLinear)\n",
    "# The order of parameters for the Euler3DTransform is [angle_x, angle_y, angle_z, t_x, t_y, t_z]. The parameter \n",
    "# sampling grid is centered on the initial_transform parameter values, that are all zero for the rotations. Given\n",
    "# the number of steps, their length and optimizer scales we have:\n",
    "# angle_x = 0\n",
    "# angle_y = -pi, 0, pi\n",
    "# angle_z = -pi, 0, pi\n",
    "registration_method.SetOptimizerAsExhaustive(numberOfSteps=[0,1,1,0,0,0], stepLength = np.pi)\n",
    "registration_method.SetOptimizerScales([1,1,1,1,1,1])\n",
    "\n",
    "#Perform the registration in-place so that the initial_transform is modified.\n",
    "registration_method.SetInitialTransform(initial_transform, inPlace=True)\n",
    "registration_method.Execute(fixed_image, moving_image)\n",
    "\n",
    "gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=initial_transform, fixed_window_level=ct_window_level, moving_window_level=mr_window_level);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Manual initialization\n",
    "\n",
    "When all else fails, a human in the loop will almost always be able to robustly initialize the registration.\n",
    "\n",
    "In the example below we identify corresponding points to compute an initial rigid transformation. \n",
    "\n",
    "**Note**: There is no correspondence between the fiducial markers on the phantom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point_acquisition_interface = gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), fixed_window_level=ct_window_level, moving_window_level=mr_window_level);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the manually specified points and compute the transformation.\n",
    "\n",
    "fixed_image_points, moving_image_points = point_acquisition_interface.get_points()\n",
    "\n",
    "# Previously localized points (here so that the testing passes):\n",
    "fixed_image_points = [(24.062587103074605, 14.594981536981521, -58.75), (6.178716135332678, 53.93949766601378, -58.75), (74.14383149714774, -69.04462737237648, -76.25), (109.74899278747029, -14.905272533666817, -76.25)]\n",
    "moving_image_points = [(4.358707846364581, 60.46357110706131, -71.53120422363281), (24.09010295252645, 98.21840981673873, -71.53120422363281), (-52.11888008581127, -26.57984635768439, -58.53120422363281), (-87.46150681392184, 28.73904765153219, -58.53120422363281)]\n",
    "\n",
    "fixed_image_points_flat = [c for p in fixed_image_points for c in p]        \n",
    "moving_image_points_flat = [c for p in moving_image_points for c in p]\n",
    "initial_transformation = sitk.LandmarkBasedTransformInitializer(sitk.VersorRigid3DTransform(), \n",
    "                                                                fixed_image_points_flat, \n",
    "                                                                moving_image_points_flat)\n",
    "gui.RegistrationPointDataAquisition(fixed_image, moving_image, figure_size=(8,4), known_transformation=initial_transform, fixed_window_level=ct_window_level, moving_window_level=mr_window_level);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"05_advanced_registration.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
}