32_registration_3d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>SimpleITK Tutorial</center></h1>\n",
    "\n",
    "# Registration of 3D Images\n",
    "\n",
    "Author: José D. Tascón-Vidarte \\\n",
    "Email: jota@di.ku.dk \\\n",
    "University of Copenhagen\n",
    "\n",
    "Objectives:\n",
    "- Learn how to register 3D images\n",
    "\n",
    "The documentation of SimpleITK is available at https://itk.org/SimpleITKDoxygen100/html/index.html\n",
    "\n",
    "This tutorial is based on the Notebooks provided by the Insight Software Consortium. http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Cloning files\n",
    "\n",
    "Open the files to be used from github."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !rm -rf /content/registration_tutorial/\n",
    "# !git clone https://github.com/josetascon/registration_tutorial.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import libraries\n",
    "Import the libraries to be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install simpleitk           # install simpleitk in colabs\n",
    "\n",
    "import SimpleITK as sitk         # simpleitk library\n",
    "import sys                       # sys library, used to include local py files\n",
    "import numpy as np               # array and matrix library\n",
    "import matplotlib.pyplot as plt  # plot library\n",
    "\n",
    "from ipywidgets import interact, fixed # jupyter widgets libraries\n",
    "from IPython.display import clear_output # jupyter clean out\n",
    "\n",
    "# sys.path.insert(0, '/content/registration_tutorial/') # local library path\n",
    "from local_utils import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Registration Problem\n",
    "\n",
    "Having two images, a reference (fixed) image $I_{1}(x)$ and an input (moving) image $I_{2}(x)$, find a transformation $\\phi(x)$ that maps $I_{2}$ into image $I_{1}$ in order to them to be similar.\n",
    "\n",
    "Therefore, registration consist of finding a transformation between images. In this notebook only rigid or affine transformation are used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image Registration Framework in SimpleITK\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",
    "![Registration](./images/docs/itk_registration_framework.png \"Registration Framework\")\n",
    "\n",
    "In ITK v4 the registration method accepts three transformations (if you look at the diagram above you will only see two transformations, Moving transform represents $T_{opt} \\circ T_m$):\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",
    "* SetFixedInitialTransform $T_f$- maps points from the virtual image domain to the fixed image domain, never modified.\n",
    "* SetMovingInitialTransform $T_m$- maps points from the virtual image domain to the moving image domain, never modified.\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})))$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiresolution Registration.\n",
    "\n",
    "In order to achieve a improvement in the optimization procedure it is common to register images in a multiresolution style. This algorithm start with the same images with lower resolution. Since there is a significant reduction in the amount of data this lower resolution images converge faster. The transformation result of this registration is then used as the input of the next resolution. The algorithm reach the original resolution and find the proper transformation. This scheme is depicted in the figure, where the lower resolution images are shown as a pyramid.\n",
    "\n",
    "![Registration](./images/docs/pyramid_registration.png \"Multiresolution Framework\")\n",
    "\n",
    "SimpleITK incorporate this kind of multiresolution algorithm. The features or the multiresolution framework are:\n",
    "* Multi scale - specify both scale, and how much to smooth with respect to original image.</li>\n",
    "* Initial transformation modified in place, so in the end we have the same type of transformation in hand.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Images\n",
    "\n",
    "In the following example we use images form the Retrospective Image Registration Evaluation (<a href=\"http://www.insight-journal.org/rire/\">RIRE</a>) project.\n",
    "\n",
    "In particular, the training_001 images are used. The CT Scan and the MRI T1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "===== Image Information =====\n",
      "Pixel type: \t\t64-bit float\n",
      "Pixel channels: \t1\n",
      "Dimensions: \t\t3\n",
      "Size: \t\t\t(512, 512, 29)\n",
      "Length (mm): \t\t(0.653595, 0.653595, 4.0)\n",
      "Total Elements: \t7602176\n",
      "\n",
      "\n",
      "===== Image Information =====\n",
      "Pixel type: \t\t64-bit float\n",
      "Pixel channels: \t1\n",
      "Dimensions: \t\t3\n",
      "Size: \t\t\t(256, 256, 26)\n",
      "Length (mm): \t\t(1.25, 1.25, 4.0)\n",
      "Total Elements: \t1703936\n",
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "39e708125c5540188554f18b9573c89f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=14, description='range_z', max=28), Output()), _dom_classes=('widget-int…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e8c874bab78d4719ae72847bc3381261",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=12, description='range_z', max=25), Output()), _dom_classes=('widget-int…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Read the image\n",
    "dim = 3 # dimension 2 (images in 2D)\n",
    "\n",
    "prefix = './images/'\n",
    "# prefix = '/content/registration_tutorial/images/'\n",
    "\n",
    "# Load the fixed image\n",
    "fixed_image = sitk.ReadImage(prefix + 'training_001_ct.mha', sitk.sitkFloat64) # load as double\n",
    "\n",
    "# Load the moving image\n",
    "moving_image = sitk.ReadImage(prefix + 'training_001_mr_T1.mha', sitk.sitkFloat64) # load as double\n",
    "\n",
    "# Print the information\n",
    "print(image_info(fixed_image))  # Fixed Image. CT Scan\n",
    "print(image_info(moving_image)) # Moving image. MRI T1\n",
    "\n",
    "# Plot the images with the axial point of view\n",
    "interact(imshow_axial, range_z =(0,fixed_image.GetSize()[2]-1), image_itk = fixed(fixed_image));\n",
    "interact(imshow_axial, range_z =(0,moving_image.GetSize()[2]-1), image_itk = fixed(moving_image));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions\n",
    "\n",
    "The next functions are used with the registration algorithm to plot the metric values during the optimization procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Callback invoked when the StartEvent happens, sets up our new data.\n",
    "def start_register_plot():\n",
    "    global metric_values, multires_iterations\n",
    "    \n",
    "    metric_values = []\n",
    "    multires_iterations = []\n",
    "\n",
    "# Callback invoked when the EndEvent happens, do cleanup of data and figure.\n",
    "def end_register_plot():\n",
    "    global metric_values, multires_iterations\n",
    "    \n",
    "    del metric_values\n",
    "    del multires_iterations\n",
    "    # Close figure, we don't want to get a duplicate of the plot latter on.\n",
    "    plt.close()\n",
    "\n",
    "# Callback invoked when the IterationEvent happens, update our data and display new figure.    \n",
    "def plot_register_values(registration_method):\n",
    "    global metric_values, multires_iterations\n",
    "    \n",
    "    metric_values.append(registration_method.GetMetricValue())                                       \n",
    "    # Clear the output area (wait=True, to reduce flickering), and plot current data\n",
    "    clear_output(wait=True)\n",
    "    # Plot the similarity metric values\n",
    "    plt.plot(metric_values, 'r')\n",
    "    plt.plot(multires_iterations, [metric_values[index] for index in multires_iterations], 'b*')\n",
    "    plt.xlabel('Iteration Number',fontsize=12)\n",
    "    plt.ylabel('Metric Value',fontsize=12)\n",
    "    plt.show()\n",
    "    \n",
    "# Callback invoked when the sitkMultiResolutionIterationEvent happens, update the index into the \n",
    "# metric_values list. \n",
    "def update_multires_iterations():\n",
    "    global metric_values, multires_iterations\n",
    "    multires_iterations.append(len(metric_values))\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Axial View of Both Images\n",
    "\n",
    "We can plot the images with a variable alpha acting as opacity level. When the alpha is zero the fixed image is visualized, and when the alpha is one the moving image is shown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Callback invoked by the IPython interact method for scrolling and modifying the alpha blending\n",
    "# of an image stack of two images that occupy the same physical space. \n",
    "def display_images_with_alpha(image_z, alpha, fixed, moving):\n",
    "    img = (1.0 - alpha)*fixed[:,:,image_z] + alpha*moving[:,:,image_z] \n",
    "    plt.imshow(sitk.GetArrayViewFromImage(img),cmap=plt.cm.Greys_r);\n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Registration Algorithm\n",
    "\n",
    "Again, the registration algorithm will use the two images to find the transformation. Since the images to be used have different information (histograms), the registration algorithm will use mutual information as the similarity metric.\n",
    "\n",
    "Additionally, we will add the multiresolution framework to the registration method. Another feature that is allowed in SimpleITK is to reduce the number of voxels used in the metric. This is perform by sampling the voxels (typically random sampling) and reducing the amount of data involved in the computation. This feature will be tested also\n",
    "\n",
    "The next function register two images with the mutual information as metric and the gradient descent as the optimizer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def registration_mi(fixed_image, moving_image, transform, \n",
    "                    interpolator = sitk.sitkLinear,\n",
    "                    bins = 50, sampling = 0.02,\n",
    "                    num_iterations = 50, learning_rate = 1.5, \n",
    "                    multiresolution= True, verbose = True, plot = True):\n",
    "    '''\n",
    "    Image regristration with metric mutual information (mi)\n",
    "    Input:\n",
    "        fixed_image: sitk.Image\n",
    "        moving_image: sitk.Image\n",
    "    Output:\n",
    "        sitk.Transform\n",
    "    '''\n",
    "    # Define the registration object class\n",
    "    registration_method = sitk.ImageRegistrationMethod()\n",
    "\n",
    "    # Set transform, intepolation and metric\n",
    "    registration_method.SetInitialTransform(transform)\n",
    "    registration_method.SetInterpolator(interpolator)\n",
    "    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=bins)\n",
    "    \n",
    "    # Set the sampling method\n",
    "    registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)\n",
    "    registration_method.SetMetricSamplingPercentage(sampling)\n",
    "    \n",
    "    # Set optimizer as gradient descent\n",
    "    registration_method.SetOptimizerAsGradientDescent(learningRate=learning_rate, \n",
    "        numberOfIterations=num_iterations, convergenceMinimumValue=1e-6, convergenceWindowSize=10)\n",
    "    registration_method.SetOptimizerScalesFromPhysicalShift() # Set appropiate scales\n",
    "    \n",
    "    # Setup for the multi-resolution framework.\n",
    "    if multiresolution:\n",
    "        registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])\n",
    "        registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])\n",
    "        registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()\n",
    "\n",
    "    # Enable plotting\n",
    "    if plot:\n",
    "        registration_method.AddCommand(sitk.sitkStartEvent, start_register_plot)\n",
    "        registration_method.AddCommand(sitk.sitkEndEvent, end_register_plot)\n",
    "        if multiresolution: registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, \n",
    "                                                           update_multires_iterations)\n",
    "        registration_method.AddCommand(sitk.sitkIterationEvent, \n",
    "                                   lambda: plot_register_values(registration_method))\n",
    "    \n",
    "    transform_estimated = registration_method.Execute(fixed_image, moving_image)\n",
    "    \n",
    "    # Print the information about the optimization in the registration \n",
    "    if verbose:\n",
    "        info = '\\n===== Registration Info ====='\n",
    "        info += '\\nFinal metric value: \\t' + str(registration_method.GetMetricValue())\n",
    "        info += '\\nLast iteration: \\t' + str(registration_method.GetOptimizerIteration())\n",
    "        info += '\\nStop condition: \\n'\n",
    "        info += (str(registration_method.GetOptimizerStopConditionDescription()))\n",
    "        info += '\\n'\n",
    "        print(info)\n",
    "    \n",
    "    return transform_estimated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The registration algorithm is ready to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial Transform\n",
    "\n",
    "We proceed to initilize the transformation. We will use a rigid transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "===== Transform Info =====\n",
      "Transform type: \tEuler3DTransform\n",
      "Dimensions: \t\t3\n",
      "Parameters: \t\t(0.0, 0.0, 0.0, -7.618522500000012, -7.618522500000012, -6.0)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "initial_transform = sitk.CenteredTransformInitializer(fixed_image, \n",
    "                                                      moving_image, \n",
    "                                                      sitk.Euler3DTransform(), \n",
    "                                                      sitk.CenteredTransformInitializerFilter.GEOMETRY)\n",
    "\n",
    "print(transform_info(initial_transform))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial estimation is visualized with the utility function with alpha defined before. To be able to plot both images in the same space is necesarry to put the moving image in the space of the fixed. This resampling is only made for visualization purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4abc9d75a2e5447da4d25519c8823dca",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=14, description='image_z', max=29), FloatSlider(value=0.5, description='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Resample the moving image to the fixed image space\n",
    "moving_resampled = sitk.Resample(moving_image, fixed_image, initial_transform, \n",
    "                                 sitk.sitkLinear, 0.0, moving_image.GetPixelID())\n",
    "\n",
    "# Visualize the image in the fixed image space\n",
    "interact(display_images_with_alpha, image_z=(0,fixed_image.GetSize()[2]), \n",
    "         alpha=(0.0,1.0,0.05), fixed = fixed(fixed_image), moving=fixed(moving_resampled));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we correctly initialize the transformation, we proceed to run the registration algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAENCAYAAAAhRzNRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYXGWZ9/Hv3d3ZSchKCNkTtmAggFFgxATCYkAcmBERAWFmcHhxf9FRcXAZib7iMqPDgCgDCgIuozCyiGwhgGxKgEjCGgyB7HtCQjpr3+8f9ymq0qnurnRX1anl97muuk7VqdNVd3Wn+5fnec55HnN3RERE9lRD2gWIiEh1UoCIiEinKEBERKRTFCAiItIpChAREekUBYiIiHSKAkRERDpFASIiIp2iABERkU5pSruAUho8eLCPGTMm7TJERKrK008/vdrdh3R0XE0HyJgxY5g9e3baZYiIVBUze72Q49SFJSIinaIAERGRTlGAiIhIpyhARESkUxQgIiLSKQqQNixbBlOnwvLlaVciIlKZFCBtmDEDHn0ULr887UpERCqTAqSVXr3ADK65BlpaYmsW+0VEJEsB0sqCBXDOmVvpzVsA9O4N554Lr72WcmEiIhVGAdLKsGHQb3APttCLnmxhyxanXz/Yd9+0KxMRqSwKkDxWrICLT17AkxzFxact0UC6iEgeNT0XVmfddhvQPByGvsbVg74GP/1p2iWJiFQctUDa0qsXfPCDcOutsGVL2tWIiFQcBUh7zjkH3nwTfv/7tCsREak4CpD2TJsWo+e33JJ2JSIiFUcB0p7GRvjIR+Cuu2D16rSrERGpKAqQjvzjP8L27XDzzWlXIiJSURQgHTn0UHj3u+H668E97WpERCqGAqQQF14I8+bBU0+lXYmISMVQgBTi7LNjTpPrr0+7EhGRiqEAKUS/fvChD8EvfwlvvZV2NSIiFUEBUqgLL4SNG5PL1EVERAFSqGOPjZkW77wz7UpERCqCAqRQZnDqqXDvvXFar4hInVOA7In3vz+mNnn88bQrERFJnQJkT5x4InTrprmxRERQgOyZvn1hyhQFiIgIFRAgZjbQzO43s/nJdkA7x/YzsyVmdlU5a9zF+98PL7wACxemVoKISCVIPUCAS4GZ7n4AMDN53JYZwMNlqaotp54a27vvTrUMEZG0VUKAnA7cmNy/ETgj30Fm9k5gKHBfmerK78ADYfx4dWOJSN2rhAAZ6u7LAJLtPq0PMLMG4N+BL5S5tt2ZwRlnwH33weuvp12NiEhqyhIgZvaAmc3Lczu9wJf4BHC3uy8q4L0uMrPZZjZ71apVXSu8LZ/9bATJt79dmtcXEakC5ilPUW5mLwPHufsyMxsGPOTuB7U65hbgvUALsBfQHfiRu7c3XsLkyZN99uzZpSn84x+PyRVffRVGjSrNe4iIpMDMnnb3yR0dVwldWHcAFyT3LwBub32Au5/r7qPcfQzwL8DPOwqPkvvyl2N7xRWpliEikpZKCJArgJPMbD5wUvIYM5tsZtelWll7Ro2K1Qqvvx4WddizJiJSc1IPEHdf4+4nuPsByXZtsn+2u38sz/E3uPunyl9pHpdeCtu2xTTvIiJ1JvUAqWpjx8LBB8Mjj6RdiYhI2SlAumrKFHj0Udi5M+1KRETKSgHSVVOmwIYNMHdu2pWIiJSVAqSrpkyJrbqxRKTOKEC6auTIGAtRgIhInVGAFMOUKREgKV+UKSJSTgqQYpgyBVatgpdfTrsSEZGyUYAUg8ZBRKQOKUCKYfx4GDYMHk53qRIRkXJSgBSDGUydCn/4AyxdmnY1IiJloQAplq9+NaY1Oess2L497WpEREpOAVIshxwC110Hjz0GX/xi2tWIiJScAqSYzj4bPvMZ+OEPteStiNQ8BUixfe97cNBB0QrR/FgiUsMUIMXWvTvMmAEvvAC33JJ2NSIiJaMAKYUPfhCOPBK+/vUYWBcRqUEKkFJoaIBvfxsWLoRrr027GhGRklCAlMpJJ8W1Id//ftqViIiUhAKkVMzg7/8eXn8dFi9OuxoRkaJTgJTSMcfE9okn0q1DRKQEFCClNGkS9OypABGRmqQAKaXu3WHyZAWIiNQkBUipHX00PPMMbN2adiUiIkWlACm1Y46Ja0GeeSbtSkREikoBUmoaSBeRGqUAKbVhw2D0aHjyybQrEREpKgVIORxzjFogIlJzFCDlcMwxcTGhLigUkRqiACmHo4+O7WOPpVuHiEgRKUDK4cgjYe+94f77065ERKRoFCDl0NQEJ54I99wD7mlXIyJSFAqQcpk+HZYsgeefT7sSEZGiUICUy/Tpsf3DH9KtQ0SkSBQg5TJiBEycGN1YIiI1IPUAMbOBZna/mc1PtgPaOG6nmc1JbneUu86imD4d/vhH2LQp7UpERLqs4AAxswlm9lUzuzp5fLCZHVaEGi4FZrr7AcDM5HE+ze5+eHL72yK8b/mdcgps3w6zZqVdiYhIlxUUIGb2IeBhYDjw0WT3XsB/FKGG04Ebk/s3AmcU4TUr03veA336aBxERGpCoS2Qy4GT3f1iYGey7y/ApCLUMNTdlwEk233aOK6nmc02syfNrDpDpkcPmDYN7r037UpERLqsqcDj9iECA8BztgVd1GBmDwD75nnqsgLfH2CUuy81s3HAg2Y2193/mue9LgIuAhg1atQevHyZTJsGd94Z05qMGJF2NSIinVZoC+Rpsl1XGWcDfy7ki939RHefmOd2O7DCzIYBJNuVbbzG0mS7AHgIOKKN465198nuPnnIkCGFlFdeU6bE9pFH0q1DRKSLCg2QzwDfNLOHgT5mdi8wA7ikCDXcAVyQ3L8AuL31AWY2wMx6JPcHA+8BXijCe5ffpEnQr58CRESqXkFdWO7+kpkdDJwG3AUsAu5y92Kcj3oF8D9mdiHwBvAhADObDFzs7h8DJgA/MbMWIvSucPfqDJDGRjj2WAWIiFS9QsdAcPfNwP8UuwB3XwOckGf/bOBjyf3HgUOL/d6pmTIF7r4bVq6Efdo6Z0BEpLIVFCBm9kfaGDB39ylFrageZMZB/vhH+OAH061FRKSTCm2BXNfq8b7AhcDNxS2nTrzzndC7d3RjKUBEpEoVOgZyY+t9ZnYr8DPiGhHZE927xyqFGgcRkbasWgVbt1b06f5dmQtrCVCMqUzq05Qp8Je/wLp1aVciIpXos5+FM89Mu4p2FToG8k+tdvUG/h54sugV1YspU2Jxqcceg9NOS7saEak0ixfDG2+kXUW7Ch0DaX0R4VvA48APiltOHTksabzNn59uHSJSmdavh9Wr4z+aZrFv0ybYuTOWyK4AhY6BHF/qQurOgAHQq1f8L0NEpLV162L27jffzAbGxz8eK5s++GC6tSXaDJBkzqkOJVOLyJ4yi8ExBYiI5LN+fWxXr84GyMsvV1S3VnstkFeJaz+snWMcaCxqRfVEASIi+Wzfnl14btUqGD8+7q9cCStWwLZtcTZnyto8C8vdG9y9Mdm2dVN4dIUCRETy2bAhe3/16ti6R3gALFuWfX7WLLjrrvLVliP1JW3r2ogRsHRpDIqJiGTknt6fCZBNm2DLlri/ZEn2+a99DS4pxry2e67Q03ibgE8AU4HB5HRraSqTLhgxAnbsiGbpsGFpVyMilSIz/gHZAMm0PmDXnouFC6NFsn07dOsW+zZuhL32yp69VSKFtkB+APwf4BHgncCtxCJTlXEqQLXKXGGqbiwRyZXbAlm1KrYrc5ZKyrRAtm6N+zt3woKc85nOOw+OOqrkZRYaIH8PnOLu/wnsSLZnADq9tysUICKST0ctkEyALFoUYyMAr7ySff6552Ds2NLWSOEB0ptYAwSg2cx6u/tLtLEqoBRIASIi+WRaIIMHZwMk0wLp0ycbIAsXZr8mEyAbN8b+Q0u/AkahV6K/CLyLWMJ2NvBvZvYmMR+WdNbgwXEqngJERHJlWiAHHJDtwsq0QCZNyv7NyARIY2M2QObNi20ZAqTQFshngR3J/c8BRwIfAC4qRVF1o6EBhg9XgIjIrtati/9cjhy5awtkwIDomsptgTQ2xhIRmQCZOze2h5V+rtt2WyBm1t/d17v7U5l97j4fOLHkldULXQsiIq2tXw/9+8OQIbuOgQwdGv/pXLIkxj5eey1C5pBD4L774ri5c+MMrNGjS15mRy2Q5Wb2GzN7v5nposFSGDlSASIiu1q/PlobQ4Zk58TKLIE9fHhcib56dbRAxoyBAw+Ma8o2bYoAmTgxejhKrKN3eDewELgWWGpmPzCzSSWvqp5kWiCed8VgEalH69ZFC2Tw4Hi8du2uLRCIVkgmQA46KPa98kqcgVWG8Q/oIEDc/Tl3/wIwEjgfGAI8ZmbPmdnnzWzfchRZ00aMyP5vQkQEsi2QTICsXr1rCwTiuo+lS2NM5MADY99DD0X4VEKAZLh7i7vf6+7nAUOB7wOfAipnWshqpVN5RaS11i2QJUti39Ch2b8Zjz8e2zFjYrJFM7j11thXSQGSYWb7AZ8EvgQMAm4pRVF1RQEiIq3ljoEAvPhibIcOhX33jfGNRx+NfWPGxNpCo0ZlQ6VSAsTMepnZeWZ2P/AacDJwBbCvu/9jqQuseQoQEcnlvnsL5PnnY7vPPtDUFEHyzDOxb8yY2Ga6sfbbDwYNKkup7QaImd0IrAAuA2YB+7v7ie5+k7tvLkeBNS/zD0IBIiIAmzfHJKsDBmSD4IUXYjt0aGyHD48zs5qasmMimYH0MrU+oOMr0ZuBk939yXIUU5caG+N/DAoQEYHsNCb9+0OPHtCv364tEIiei9mzo9uqMbnCItMCqZQAcfeLy1VIXdPFhCKSkZnGpH//2A4enJ1pN7cFAtnuK0glQLSgVCUYOXLXqZhFpH5lWiADBsQ2Mw7Sq1dMpAj5A+S44+CrX4UzzihHlYACpDIcdVRcEKRWiIi0boFkzsQaOjS7QFS+AOnRAy6/PLq8ykQBUgmmTYvtg1qfS6TutdUCyYx/QPbszTKs+dGeggLEzM43s8Na7ZtkZh8tTVl15tBD4x+JAkRE8o2BQHb8A+Doo+HjH4fp08tbWyuFtkBmkF1QKmMR8M3illOnGhrg+OMjQDQnlkh9y7RA9t47tvlaIL17w49+lH0uJYUGSD/gzVb7NgD9i1tOHZs2LZanfPXVtCsRkTStXw99+8Y1HrDrGEiFKTRAXgA+2Grf3xErFUoxnHBCbNWNJVLfMtOYZOTrwqoQhQbIl4DrzOxWM/uumd0GXA98vqsFmNlAM7vfzOYn2wFtHDfKzO4zsxfN7AUzG9PV964o++8fA2MKEJH6lpnGJCNfF1aFKHQ23keBicBTQB9ibfSJ7v5YEWq4FJjp7gcAM5PH+fwc+J67TyDWKVlZhPeuHGbRjfXgg9DSknY1IpKW1i2Qd78bvvIVOOWU9GpqQ8Gn8br7G+5+hbt/Mtm2HlTvrNOBG5P7NwK7XQVjZocATe5+f1LLppqci+uEE2Le/3nz0q5ERNLSugXSrRvMmJEdVK8gbU5lYmbXuvtFyf2bgLynB7n7+V2sYai7L0tea5mZ5WunHQisT7rOxgIPAJe6+848dV8EXAQwatSoLpZWZlOnxvaxx+Cww9o/VkRqU+sWSAVrby6s13Lud+nUIDN7AMi3euFlBb5EE/Be4AhiEatfA/9AjMPswt2vJZbgZfLkydV1TuyoUTBwIDz7bNqViEhaWrdAKlibAeLu3wYws0bimo9fuPuWzryJu5/Y1nNmtsLMhiWtj2HkH9tYDDzr7guSr/kdcDR5AqSqmcERRyhAROrVjh2wcWPVBEiHYyBJN9F/dDY8CnAHcEFy/wLg9jzHPAUMMLPkhGimEacW154jjoC5c2OufxGpLxs2xLZKurAKHUS/08w+UKIargBOMrP5wEnJY8xsspldB2+H2L8AM81sLmDAf5eonnQdcQRs3QovvZR2JSJSbq2nMalwHS0oldET+K2ZPUF0Z709ttDVQXR3XwOckGf/bOBjOY/vB2p/ZPnww2M7Z05Z5/UXkQqQmZG7Ai8azKfQAJmX3KTUDjoo5v1/9ln4qOaqFKkrc+bEdtKkdOsoUKEB8hN3X956p5nlO7NKuqKxMU7h1UC6SP2ZMydaH/tWx5/WQsdAXmljf20OZKftiCPiH5Jm5hWpL3PmZLuxq0ChAWK77TDrB2jOjVI4/PAYTHv99bQrEZFy2bYNnn++qgKk3S4sM8sMmPcyszdaPT0I+GWpCqtrRxwR22ef3XXJShGpXS+9FKfv10qAAOcRrY+7gdwRXQdWuPvLpSqsrh16aIyFPPss/N3fpV2NiJRDZgC9VgLE3R8GMLPBNTl5YaXq1QsOPlgD6SL1ZM6c+N0/4IC0KylYoWMgO83sW2a2wMw2AJjZyWb2qRLWVt8mTYK//CXtKkSkXObMiTMwGxvTrqRghQbID4n1QM4lexHh88DHS1GUAIccEkvcbtqUdiUiUmruVXcGFhQeIGcA57j7EyRnXrn7EmB4qQqrexMmxFZTmojUvkWLYhbeGg2QbbQaL0kmNlxT9IokZALkRS07L1Lrlj34IlN5iOUj35V2KXuk0AD5DXCjmY0FSKZdvwr4VakKq3v77w9NTQqQNixbFutvLd9tfgSR6jPj6oE8yrFc/rvqmu6v0AD5V2AhMBfoD8wHlgLfKE1ZQrducTaGAiSvGTPg0Ufh8svTrkSk83r1imWArpn9Llpo5JrrumEW+6tBQQHi7tvc/f+6+17AUKCvu1/i7ttKW16dmzBBAdLK279w10BLS2yr6RdOJNeCBXDOOdC7oRmA3r3h3HPhtdc6+MIK0W6AmNmo1jegFzAy57GUyoQJ8OqrMcWBADm/cN13ALGtpl84kVzDhkG/frClpTs9G7exZUs8rpK5FDu8En0h2dN2d5sPK3muek5arjYTJsDOnREihxySdjUV4e1fuO2N9KSZLdt7VNUvnEhrK5a3cDE/4aJ/cq7t9kmWLUu7osJ11IX1HDHe8RVgNNCt1a17Saurd5kzsV7QpMe5VqyAi6e8wJMczcVj7tVAulS12/5rKVfzSSa9s4mrr4bbbku7osJ1NJXJ4WY2kVir/FHgJeDnwG3u3lyG+urbQQfFVuMgu7jtNuA/H4CHn+Pqxs/AbfPTLkmk85Ysie2IEenW0QkdDqK7+zx3/wIwFvgP4DRgmZkdWeri6l6fPjB6tAIkn83J1Gx//auu1pfqllnGdnj1XZdd6Gm8AAcAU4FjgGeBdSWpSHalM7Hya04awO6xhoJItcoESK21QMxsoJl90sz+DPwO2ARMcffj3V3nvZTDhAnw8stxzqpkNef0oD73XHp1iHTVkiXQowcMGpR2JXuso7OwlgKvATcBTyb79jez/TMHuPuDJapNIAKkuTlWJxw7Nu1qKkdzMwwYADt2KECkui1eHN1Xlu9E18rWUYAsB3oC/5zcWnNgXLGLkhyZ03fnzVOA5GpujjGiUaMUIFLdFi+uyu4r6KALy93HuPvYdm4Kj1I74gjYay/43/9Nu5LKsnlzXH5+2GERIO4df41IJVqypDYDRCpA795w1lnwm9/AW2+lXU3laG6OAJk0Cdavzw5EilQT92wXVhVSgFSDCy6IU1Wr6QqjUmtujnA9LJm9VKs3SjVavTqmKlILRErm2GNh3Di44Ya0K6kcmRbIxInxWOMgUo2q+CJCUIBUh4YGOP98mDUrzsaS7BhIv35xcsGVV8KnPw133QWrVqVdnUhhqvgiQlCAVI/zz4/+0ptuSruSypBpgQBcdRVMngzXXw8f+ADss0+cnfW+98GFF8I3vgG33AJPPRUrUW3Zkm7tIhlV3gLp6DReqRRjx8YSfLfcAl/5StrVpC8zBgJw6qlx27IFnngCnn46bn/9K8ydG8sWtj5Lq08fGD8+Fu0aPx7GjIluwsMPh6FDy/5xpE4tXhw9DFX6b04BUk3OOAMuuSQWv6j3a0JyWyAZPXvC8cfHLdeWLbGQyCuvRJisWwcrV8Y0+XPnwh13wPbt2eOHD4+JLEeMiPvDhsV88fvtByNHxrZJvzpSBIsXx7+vKv33VJ1V16vp0yNA7r0XLr447WrSlRkDKUTPnnFBZltrqrS0RLDMnw/PPBO3BQtizGnp0liTJZcZDB4coTJsWATKfvvFvkGDYO+9oW/fuH6nX7943NQU72MW+3r2zF+HWVVekSydVMXXgIACpLocdFDMznvPPQqQfC2QzmpoyIbA1Km7PtfSAmvWxNjJkiXxP8bFi2NRkuXLY/8LL8S2ddC0p3v3qL8xWY9t8+ZoKZlF11zv3tCtWwSPWXTBNTREQA0aFAFkFs/vtVfceveO/b16RYD17Rv1b94cp4r26hW3/v1h4MAIMvc4ZtCg+PzduhXneyqFWbw4u+5PFUo9QMxsIPBrYAyxAuJZ7r6u1THHAz/I2XUwcLa7/65MZVYGMzjlFLj55viD0L1O1/PauTO6nMqxEHpDAwwZErfMNSf5tLTAm2/GGWBvvhnX7WzcGPc3bIiaGxqyx23YAFu3xlxe7jEm07t3HLd5c9y2b892rTU0xLFr18a1A1u3xtft2BEXmG7cGAHU3Nz5q/IbGiJYMvbeO/rm9947/r1t3RrPjxoV85CtXx+3Xr1if58+cVwmrPr2jYBcuTK+L01Nsa9796hzy5b4mkGD4usHDIhbr14xuWC3btnWWLdu8XWZ/d27x2s1VPl5QIsXw0knpV1Fp6UeIMClwEx3v8LMLk0efyn3AHefBRwObwfOq8B95S60IkyfDj/+MTz+OBx3XNrVpCMzE29mEL0SNDTE/+z790+3Dvf4Q58Jr8bGbGtmy5YIpg0bolW1cWP2D/Dq1bBo0a6nQK9fHy2tlSujZdO9e4y/PfJIPLf33vF5t2yJ19sR69TTvXuESEZjY3Tv7dwZ77l1a7a1tGnTrsfuicbGCPZBgyKI+vSJUOnXL4Jm69a4NTXF4549o6XWp0/U2NQU369MYEN8P7p1ywZ65vuTae01NcX7NjVlg6xbt92Pa2iI9850tR5wQATvokXRYm1uji7QjRvVhdVFpwPHJfdvBB6iVYC0cibwB3ffXNqyKtS0afEP9J576jdAMr/s5WiBVBuz+EPZs2f8cS2VlpZd//efaQ1lutx27oxw2LEjWhWtj820LNyjBbV2bZzcsG5dBNLWrdnWV+a1M4GwfXuEzpo1EW5r18ZrbNoU10m9+WY836NH/IHPfG1zcxzXnGcx1UyXYEtLvH45l08YNap871VklRAgQ919GYC7LzOzfTo4/mxiZcT61LdvXJl+zz1wxRVpV5OOzB8ABUh6Wncdme06ftLYGC2UfHJPEjDLjuGU6w9pS0uESqbF1LPn7gGXaT1kugNbWrJdpzt3xtdmuhi3bdv9uJ07s+NRmzbFCRoLF8ZZfIccEi2chQuj5Xf66eX53CVQlgAxsweAffM8ddkevs4w4FDg3naOuQi4CGBUFSd7u6ZPh0svjTOE9tsv7WrKTwEiXdHQEC2TtsYQc1txxXL00bvve8c7ivf6KSnLCJS7n+juE/PcbgdWJMGQCYiV7bzUWcD/uvv2tg5w92vdfbK7Tx5SyiZ8mqZPj+199TkMVJFjICJ1qBJOYbgDuCC5fwFwezvHfgT4ZckrqnSHHRYDcPe22RCrbWqBiFSESgiQK4CTzGw+cFLyGDObbGbXZQ4yszHASODhFGqsLGZw8snRAtmTaw9qhQbRRSpC6gHi7mvc/QR3PyDZrk32z3b3j+Uct9Ddh7t7GU+PqGDTp8fZJ08/nXYl5acWiEhFSD1ApJNOOilaIvXYjaUxEJGKoACpVoMHxxTm99yTdiXlpxaISEVQgFSz970PnnwyLr6qJxoDEakICpBqNn16XLg0c2balZSXWiAiFUEBUs2OOiqu9q23cRCNgYhUBAVINWtqivmwHnkk7UrKq7k5psrQ1OMiqVKAVLuDD44ZUuvpepA9WUxKREpGAVLtxo2LCd2WLEm7kvIp5mJSItJpCpBqN358bP/613TrKCcFiEhFUIBUu3HjYrtgQbp1lFNzswbQRSqAAqTajRwZg+n1FiBqgYikTgFS7ZqaYPTo+urC0iC6SEVQgNSCcePUAhGRslOA1IJ6DBCNgYikTgFSC8aPhzVrYMOGtCspD7VARCqCAqQW1NuZWBoDEakICpBaUG8BohaISEVQgNSCTIDUy5lYGgMRqQgKkFqw994waFB9tEDc1QIRqRAKkFoxfnx9BMj27TFxpAJEJHUKkFoxblx9dGFpMSmRiqEAqRXjxsHrr8OOHe0ft349nHkmvPFGeeoqNi0mJVIxFCC1Yvz46NpZtKj9437/e7j1Vvj1r8tTV7GpBSJSMRQgtSJzJtarr7Z/3MMPx/ahh0paTsls3hxbBYhI6hQgtWLixNj++c/tH5cJkEcf7bi7qxKpBSJSMRQgtWLwYJg0CR58sO1jli2DV16Bd74T3nwT5swpX33FojEQkYqhAKklJ5wAjz2W/SPbWqb18bWv7fq4UG+9BZdeGgPxaVELRKRiKEBqybRpsHUrPPFE/ucffhj69oVTT4UDD9zzcZDbb4fvfAduvrnLpXaaAkSkYihAasmUKdDYCDNn5n/+4Yfh2GNjEaqpU+GPf4wztwo1a1Zs77qr67V2lgbRRSqGAqSW9O0L7353/nGQFSvgxRcjOACOOy6mf//LXwp//czrzpoFGzd2udxOUQtEpGIoQGrNCSfAU0/FIHmuRx6J7XHHxTYTJIV2Y73xRkyVcsYZsG0bPPBAMardcxpEF6kYCpBaM21adEtlAmPrVrj/fvjRj6BPHzjyyNg/fDjsvz/ceWdh3ViZ7quvfhX690+vG0stEJGKoQCpNcccAz17wtVXw4c/HLP0nnxyDKx/4QvQrVv22IsvjhbIhz8MW7bEvgUL4nTf1mbNitc6/HCYPj2uaG9pKctH2oXGQEQqRlPaBUiR9ewZg+n33AP77APnngsf+EC0TFp3+3z+89DQAJ/7XMyjtXUrzJ0bV7W/9FI2bNwjQI47Lo4/7TT41a/gT3+KCxKvuioeH3NM6T9fc3PU1dhY+vcSkXal3gIxs4Fmdr+ZzU+2A9o47rtm9ryZvWhmV5qZlbttxr5gAAAMM0lEQVTWqvHTn8Yf9qVL4Sc/iT/4bY0ZXHIJ3HILLFwYXVOf/GS0Qm66KXvMa6/FGMjxx8fjU06JIDnlFPjiF2H5crjoophqvdS0mJRIxUg9QIBLgZnufgAwM3m8CzP7G+A9wGHAROBdwNRyFllVhg+H97yn8P+ln3MOrFoV4yb/9V8weTJ885vZQMicfZUJkIED4cQTwSyC5te/hnnz4Mori/9ZWtNiUiIVoxK6sE4Hjkvu3wg8BHyp1TEO9AS6AwZ0A1aUp7w6Ywb/9m/Rarnpptj+7GcwdChMmJA97re/ja6tfv1ie9pp8PWvw1lnRTfaypUxZjJkSHG7mzZvVoCIVIhKaIEMdfdlAMl2n9YHuPsTwCxgWXK7191fzPdiZnaRmc02s9mrVq0qYdk17NRToxXyr/8KBx0UpwV/4xsRLhl9+0Z4QOy/8soYVB89OsZeJk6EYcOge3c46qjoJtu2bdf3cd/z60nUAhGpGGVpgZjZA8C+eZ66rMCv3x+YAIxIdt1vZlPc/ZHWx7r7tcC1AJMnT/bOVVznzGDGjBjjmDoVrrlm19ZHPmPHwi9+EWMvo0ZFiKxdC0uWxPoj550Xg/WHHx7HrlkTxy5fHgP93/lOdL11RGMgIhXD3NP9G2tmLwPHufsyMxsGPOTuB7U65gtAT3efkTz+GrDF3b/b3mtPnjzZZ8+eXarSa9/rr0cYdPV8hZYWuPfemENr/vxYenevveC9743xlGuvjW6uz30uBvH3bfV/jZ07Y4xl7Ni4kHHHjux1LiJSdGb2tLtP7ui4SujCugO4ILl/AXB7nmPeAKaaWZOZdSMG0PN2YUkRjR7d9fCA7Blbt9wS65WsWRPhdPPN0fX1wgvx/De/GYF13nkxceNbb8Xpw5MnR8ulf/+Yv0tdWCIVoRIG0a8A/sfMLiSC4kMAZjYZuNjdPwb8FpgGzCUG1O9x9ztTqleKbdy4GJSfPz+uKbnhhgibbt3iTLBRo+LCyLVr4emnoxUiIqlLvQurlNSFVaW2bYuWxt13w377wSc+oVaHSBkV2oVVCS0QkV117x6TQp5wQtqViEg7KmEMREREqpACREREOkUBIiIinaIAERGRTlGAiIhIpyhARESkUxQgIiLSKQoQERHplJq+Et3MVgGvd+ElBgOri1ROmmrhc9TCZwB9jkqjz5HfaHcf0tFBNR0gXWVmswu5nL/S1cLnqIXPAPoclUafo2vUhSUiIp2iABERkU5RgLTv2rQLKJJa+By18BlAn6PS6HN0gcZARESkU9QCERGRTlGA5GFm083sZTN71cwuTbueQpnZSDObZWYvmtnzZvbZZP9AM7vfzOYn2wFp11oIM2s0s2fN7K7k8Vgz+1PyOX5tZt3TrrEjZtbfzH5rZi8lP5djqvHnYWaXJP+m5pnZL82sZzX8PMzsp2a20szm5ezL+/23cGXye/+cmR2ZXuW7auNzfC/5d/Wcmf2vmfXPee7Lyed42czeV6q6FCCtmFkjcDVwCnAI8BEzOyTdqgq2A/i8u08AjgY+mdR+KTDT3Q8AZiaPq8FngRdzHn8H+EHyOdYBF6ZS1Z75T2IJ5oOBScTnqaqfh5kNBz4DTHb3iUAjcDbV8fO4AZjeal9b3/9TgAOS20XANWWqsRA3sPvnuB+Y6O6HAa8AXwZIfufPBt6RfM2Pkr9rRacA2d27gVfdfYG7bwN+BZyeck0Fcfdl7v5Mcn8j8cdqOFH/jclhNwIVv6i4mY0A3g9clzw2YBrw2+SQiv8cZtYPmAJcD+Du29x9PVX48yBWL+1lZk1Ab2AZVfDzcPdHgLWtdrf1/T8d+LmHJ4H+ZjasPJW2L9/ncPf73H1H8vBJYERy/3TgV+6+1d1fA14l/q4VnQJkd8OBRTmPFyf7qoqZjQGOAP4EDHX3ZRAhA+yTXmUF+yHwRaAleTwIWJ/zC1MNP5dxwCrgZ0lX3HVm1ocq+3m4+xLg+8AbRHBsAJ6m+n4eGW19/6v5d/+fgD8k98v2ORQgu7M8+6rqVDUz2wu4Ffi/7v5m2vXsKTM7DVjp7k/n7s5zaKX/XJqAI4Fr3P0I4C0qvLsqn2SM4HRgLLAf0Ifo7mmt0n8eHanGf2OY2WVE9/UtmV15DivJ51CA7G4xMDLn8QhgaUq17DEz60aExy3ufluye0WmKZ5sV6ZVX4HeA/ytmS0kuhCnES2S/kkXClTHz2UxsNjd/5Q8/i0RKNX28zgReM3dV7n7duA24G+ovp9HRlvf/6r73TezC4DTgHM9e01G2T6HAmR3TwEHJGeYdCcGo+5IuaaCJOME1wMvuvt/5Dx1B3BBcv8C4PZy17Yn3P3L7j7C3ccQ3/8H3f1cYBZwZnJYNXyO5cAiMzso2XUC8AJV9vMguq6ONrPeyb+xzOeoqp9Hjra+/3cA5ydnYx0NbMh0dVUiM5sOfAn4W3ffnPPUHcDZZtbDzMYSJwX8uSRFuLturW7AqcRZDX8FLku7nj2o+1iiqfocMCe5nUqMH8wE5ifbgWnXugef6TjgruT+uOQX4VXgN0CPtOsroP7DgdnJz+R3wIBq/HkA3wBeAuYBNwE9quHnAfySGLfZTvzP/MK2vv9E18/Vye/9XOKss9Q/Qzuf41VirCPzu/7jnOMvSz7Hy8AppapLV6KLiEinqAtLREQ6RQEiIiKdogAREZFOUYCIiEinKEBERKRTFCAie8jMNpnZuLTr6Cozu8HMvpl2HVK9FCBSVcxsoZmdmNz/BzN7tMTv95CZfSx3n7vv5e4LSvBeC81sRTJfVmbfx8zsoWK/l0gxKECkbuVMw1FJmohp7KtKqaYLl8qmAJGqZGYTgB8DxyRdSuuT/T3M7Ptm9kbyv/kfm1mv5LnjzGyxmX3JzJYTs+QOMLO7zGyVma1L7o9Ijv8W8F7gquQ9rkr2u5ntn9zf28x+nnz962b2FTNrSJ77BzN7NKlnnZm9Zmb5JiHM9T3gX3IXB8r5zGOS927K2fd2Cyl5v8fM7Admtt7MFpjZ3yT7FyULEl3Q6mUHWyyqtNHMHjaz0TmvfXDy3NpkYaKzcp67wcyuMbO7zewt4PgCfmxSYxQgUpXc/UXgYuCJpEsp8wf3O8CBxBQi+xPTWH8t50v3BQYCo4lFgxqAnyWPRwHNwFXJe1wG/BH4VPIen8pTyn8BexPTekwFzgf+Mef5o4jpJAYD3wWuT+aTasts4CHgXzr8JuR3FDFtyiDgF8RklO8ivhfnEWG4V87x5wIzkvrmkMzomnSj3Z+8xj7AR4iFid6R87XnAN8C+gIl7UqUyqQAkZqR/GH+Z+ASd1/rsajW/yMmZMxoAb7usdhOs7uvcfdb3X1zcvy3iCAo5P0agQ8DX3b3je6+EPh34KM5h73u7v/t7juJxYuGAUM7eOmvAZ82syGF1NHKa+7+s+T9fk3Mynp58nnvA7YRYZLxe3d/xN23EvMnHWNmI4kZXhcmr7XDY6GyW8lOnghwu7s/5u4t7r6lE7VKlavEPmCRzhpCrJb3dM5/8o1YgjVjVe4fOzPrDfyAWPozszZ5XzNrTP4It2cw0B14PWff6+y6eM/yzB1335zUldsC2I27z7NYB/5Sdl3StxArcu43J6/Xel/u+7+98JC7bzKztcSaH6OBozJdg4kmYiLF3b5W6pMCRKpZ65lAVxN/IN/hsYpeIV/zeeAg4Ch3X25mhwPPkl2Up73ZRlcTs6OOJqY3h+gGa+u998TXgWeIFk3GW8m2N5BZKGzfLr7P2+tGJF1bA4m1IxYBD7v7Se18rWZirXPqwpJqtgIYkazbgru3AP8N/MDM9gEws+Fm9r52XqMvETrrzWwg8Ye79XvkveYjaaH8D/AtM+ubDEB/Dri5C58p89qvEl1Qn8nZt4oIp/PMrNHM/gkY38W3OtXMjk2+hzOAP7n7IuAu4EAz+6iZdUtu70pOXhABFCBS3R4EngeWm9nqZN+XiHUSnjSzN4EHiBZGW34I9CJaE08C97R6/j+BM5OzqK7M8/WfJloGC4iB5F8AP+3cx9nN5cTysbn+GfgCsAZ4B/B4F9/jF0RorgXeSQyqk4wHnUyMHy0luuK+Q6wDIgKg9UBERKRz1AIREZFOUYCIiEinKEBERKRTFCAiItIpChAREekUBYiIiHSKAkRERDpFASIiIp2iABERkU75/8HTO7g1I7jJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "===== Registration Info =====\n",
      "Final metric value: \t-0.5655408250240676\n",
      "Last iteration: \t9\n",
      "Stop condition: \n",
      "GradientDescentOptimizerv4Template: Convergence checker passed at iteration 9.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "final_transform = registration_mi(fixed_image, moving_image, initial_transform,\n",
    "                               bins = 50, num_iterations = 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Form the plot we can observe how the multiresolution framework works. The registration algorithm goes from images scale down 4 times, then 2 and finally at full resolution.\n",
    "\n",
    "Always remember to see why the optimizer terminated. This will help you understand whether termination is too early, either due to thresholds being too tight, early termination due to small number of iterations - numberOfIterations, or too loose, early termination due to large value for minimal change in similarity measure - convergenceMinimumValue)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The computed transform is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "===== Transform Info =====\n",
      "Transform type: \tEuler3DTransform\n",
      "Dimensions: \t\t3\n",
      "Parameters: \t\t(0.021248824107106472, -0.005986501697035448, -0.08073511199774476, 21.033862419096508, -34.196077569876664, -23.926569742540234)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(transform_info(initial_transform))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the Registration in 2D\n",
    "\n",
    "The result in visualized in here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "854781c011664aa39c6c73b5b27dcc37",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=14, description='image_z', max=29), FloatSlider(value=0.5, description='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Resample the moving image to the fixed image space\n",
    "moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, \n",
    "                                 sitk.sitkLinear, 0.0, moving_image.GetPixelID())\n",
    "\n",
    "# Visualize the image in the fixed image space\n",
    "interact(display_images_with_alpha, image_z=(0,fixed_image.GetSize()[2]), \n",
    "         alpha=(0.0,1.0,0.05), fixed = fixed(fixed_image), moving=fixed(moving_resampled));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize Registration in 3D\n",
    "\n",
    "The registration process can be visualized in a single image in 3D. We use the same principle as the axial image with alpha.\n",
    "\n",
    "Here, a function to visualize the registration is defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def display_volumes_with_alpha(alpha, fixed, moving):\n",
    "    img = (1.0 - alpha)*fixed[:,:,:] + alpha*moving[:,:,:]\n",
    "    imshow_3d(img, controls = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the visualization function is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "32dfb82537b14dc0b33cbef8a3858c0b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.5, description='alpha', max=1.0), Output()), _dom_classes=('widget-i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Resample the moving image to the fixed image space\n",
    "moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, \n",
    "                                 sitk.sitkLinear, 0.0, moving_image.GetPixelID())\n",
    "\n",
    "# Visualize the image in the fixed image space\n",
    "interact(display_volumes_with_alpha, image_z=(0,fixed_image.GetSize()[2]), \n",
    "         alpha=(0.0,1.0,0.1), fixed = fixed(fixed_image), moving=fixed(moving_resampled));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write the Results\n",
    "\n",
    "If we are satisfied with the results, we can store the solution in files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "prefix = './images/'\n",
    "# prefix = '/content/registration_tutorial/images/'\n",
    "\n",
    "sitk.WriteImage(moving_resampled, prefix + 'training_001_mr_T1_registered.mha')\n",
    "sitk.WriteTransform(final_transform, prefix + 'RIRE_training_001_CT_2_mr_T1.tfm')"
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}