4_ITK_in_Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# ITK in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Learning Objectives\n",
    "\n",
    "* Becoming familiar with the Pythonic way of using ITK\n",
    "* Understand how to use ITK object oriented interface\n",
    "* Learn how to write simple Python code with ITK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Two ways of using ITK in Python\n",
    "\n",
    "* Pythonic way\n",
    " * Functional programming\n",
    " * More concise\n",
    " * Some functions and filters are not accessible this way\n",
    "* Object oriented way\n",
    " * Full access to ITK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Let's start the Pythonic way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import itk\n",
    "\n",
    "from itkwidgets import view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "image = itk.imread(\"data/CBCT-TextureInput.png\", itk.ctype('float'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "filtered_image = itk.median_image_filter(image, radius = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a17b4011e40e4dd08a7d9b6e1812064e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Vmlld2VyKGdyYWRpZW50X29wYWNpdHk9MC4yMiwgcmVuZGVyZWRfaW1hZ2U9PGl0a0ltYWdlUHl0aG9uLml0a0ltYWdlRjI7IHByb3h5IG9mIDxTd2lnIE9iamVjdCBvZiB0eXBlICdpdGtJbWHigKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(filtered_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Pythonic exercises\n",
    "\n",
    "* In the example above, change the radius of the filter and observe the result.\n",
    "* Replace filter with `mean_image_filter`\n",
    "* Replace filter with `otsu_threshold_image_filter`\n",
    "* Visualize results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Uncomment and change the radius of the filter and observe the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# median_filtered_image = itk.median_image_filter(image, radius = XX)\n",
    "# view(median_filtered_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/4_Using_ITK_in_Python_real_world_filters_median.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Uncomment and replace filter with `mean_image_filter`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# mean_filtered_image = itk.filter(image, radius = 5)\n",
    "# view(mean_filtered_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/4_Using_ITK_in_Python_real_world_filters_mean.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Uncomment and replace filter with `otsu_threshold_image_filter`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# otsu_filtered_image = itk.filter(image)\n",
    "# view(otsu_filtered_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/4_Using_ITK_in_Python_real_world_filters_otsu.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Object oriented\n",
    "\n",
    "* Two types of C++ ITK objects\n",
    " * Smart pointers (hint: most ITK objects are)\n",
    " * \"Simple\" objects\n",
    "* Translates in two ways of creating objects in Python\n",
    " * `obj = itk.SmartPointerObjectType.New()` (use auto-completion to see if `New()` method exists)\n",
    " * `obj = itk.SimpleObjectType()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Examples of objects\n",
    "\n",
    "* With `New()` method:\n",
    " * `itk.Image`\n",
    " * `itk.MedianImageFilter`\n",
    "* Without `New()` method:\n",
    " * `itk.Index`\n",
    " * `itk.RGBPixel`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Templated types\n",
    "\n",
    "* Is my ITK type templated? (hint: yes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on itkTemplate in module itkTemplate object:\n",
      "\n",
      "itk::Image = class itkTemplate(__builtin__.object)\n",
      " |  This class manages access to available template arguments of a C++ class.\n",
      " |  \n",
      " |  This class is generic and does not give help on the methods available in\n",
      " |  the instantiated class. To get help on a specific ITK class, instantiate an\n",
      " |  object of that class.\n",
      " |  \n",
      " |  e.g.: median = itk.MedianImageFilter[ImageType, ImageType].New()\n",
      " |        help(median)\n",
      " |  \n",
      " |  There are two ways to access types:\n",
      " |  \n",
      " |  1. With a dict interface. The user can manipulate template parameters\n",
      " |  similarly to C++, with the exception that the available parameters sets are\n",
      " |  chosen at compile time. It is also possible, with the dict interface, to\n",
      " |  explore the available parameters sets.\n",
      " |  2. With object attributes. The user can easily find the available parameters\n",
      " |  sets by pressing tab in interperter like ipython\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  GetTypes(self)\n",
      " |      Helper method which prints out the available template parameters.\n",
      " |  \n",
      " |  GetTypesAsList(self)\n",
      " |      Helper method which returns the available template parameters.\n",
      " |  \n",
      " |  New(self, *args, **kwargs)\n",
      " |      Instantiate the template with a type implied from its input.\n",
      " |      \n",
      " |      Template type specification can be avoided by assuming that the type's\n",
      " |      first template argument should have the same type as its primary input.\n",
      " |      This is generally true. If it is not true, then specify the types\n",
      " |      explicitly.\n",
      " |      \n",
      " |      For example, instead of the explicit type specification::\n",
      " |      \n",
      " |        median = itk.MedianImageFilter[ImageType, ImageType].New()\n",
      " |        median.SetInput(reader.GetOutput())\n",
      " |      \n",
      " |      call::\n",
      " |      \n",
      " |        median = itk.MedianImageFilter.New(Input=reader.GetOutput())\n",
      " |      \n",
      " |      or, the shortened::\n",
      " |      \n",
      " |        median = itk.MedianImageFilter.New(reader.GetOutput())\n",
      " |      \n",
      " |      or:\n",
      " |      \n",
      " |        median = itk.MedianImageFilter.New(reader)\n",
      " |  \n",
      " |  __add__(self, paramSetString, cl)\n",
      " |      Add a new argument set and the resulting class to the template.\n",
      " |      \n",
      " |      paramSetString is the C++ string which defines the parameters set.\n",
      " |      cl is the class which corresponds to the couple template-argument set.\n",
      " |  \n",
      " |  __call__(self, *args, **kwargs)\n",
      " |      Create a process object, update with the inputs and\n",
      " |      attributes, and return the result.\n",
      " |      \n",
      " |      The syntax is the same as the one used in New().\n",
      " |      \n",
      " |      UpdateLargestPossibleRegion() is execute and the current output,\n",
      " |      or tuple of outputs if there is more than\n",
      " |      one, is returned.\n",
      " |      \n",
      " |      For example,\n",
      " |      \n",
      " |        outputImage = itk.MedianImageFilter(inputImage, Radius=(1,2))\n",
      " |  \n",
      " |  __contains__(self, key)\n",
      " |  \n",
      " |  __dir__(self)\n",
      " |      Returns the list of the attributes available in the current template.\n",
      " |      \n",
      " |      This loads all the modules that might be required by this template first,\n",
      " |      and then returns the list of attributes. It is used when dir() is called\n",
      " |      or when it tries to autocomplete attribute names.\n",
      " |  \n",
      " |  __find_param__(self, paramSetString)\n",
      " |      Find the parameters of the template.\n",
      " |      \n",
      " |      paramSetString is the C++ string which defines the parameters set.\n",
      " |      \n",
      " |      __find_param__ returns a list of itk classes, itkCType, and/or numbers\n",
      " |      which correspond to the parameters described in paramSetString.\n",
      " |      The parameters MUST have been registered before calling this method,\n",
      " |      or __find_param__ will return a string and not the wanted object, and\n",
      " |      will display a warning. Registration order is important.\n",
      " |      \n",
      " |      This method is not static only to be able to display the template name\n",
      " |      in the warning.\n",
      " |  \n",
      " |  __getattr__(self, attr)\n",
      " |      Support for lazy loading.\n",
      " |  \n",
      " |  __getitem__(self, parameters)\n",
      " |      Return the class which corresponds to the given template parameters.\n",
      " |      \n",
      " |      parameters can be:\n",
      " |          - a single parameter (Ex: itk.Index[2])\n",
      " |          - a list of elements (Ex: itk.Image[itk.UC, 2])\n",
      " |  \n",
      " |  __iter__(self)\n",
      " |      # everything after this comment is for dict interface\n",
      " |      # and is a copy/paste from DictMixin\n",
      " |      # only methods to edit dictionary are not there\n",
      " |  \n",
      " |  __len__(self)\n",
      " |  \n",
      " |  __repr__(self)\n",
      " |  \n",
      " |  get(self, key, default=None)\n",
      " |  \n",
      " |  items(self)\n",
      " |  \n",
      " |  iteritems(self)\n",
      " |      # third level takes advantage of second level definitions\n",
      " |  \n",
      " |  iterkeys(self)\n",
      " |  \n",
      " |  itervalues(self)\n",
      " |      # fourth level uses definitions from lower levels\n",
      " |  \n",
      " |  keys(self)\n",
      " |  \n",
      " |  values(self)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Static methods defined here:\n",
      " |  \n",
      " |  __new__(cls, name)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors defined here:\n",
      " |  \n",
      " |  __dict__\n",
      " |      dictionary for instance variables (if defined)\n",
      " |  \n",
      " |  __weakref__\n",
      " |      list of weak references to the object (if defined)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data and other attributes defined here:\n",
      " |  \n",
      " |  __class_to_template__ = {<class 'pyBasePython.mapDD'>: (<itkTemplate s...\n",
      " |  \n",
      " |  __doxygen_root__ = '/home/francois.budin/.virtualenvs/itk_course/local...\n",
      " |  \n",
      " |  __named_templates__ = {'itk::AbsImageFilter': <itkTemplate itk::AbsIma...\n",
      " |  \n",
      " |  __templates__ = OrderedDict([('std::cerr', <pyBasePython.ostream...hol...\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(itk.Image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* If so, you need to specify the data type\n",
    "* This is similar to `dtype` in NumPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0.])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy\n",
    "numpy.array([0,0], dtype=numpy.float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Define and create a simple object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "itkIndex3 ([0, 0, 0])\n"
     ]
    }
   ],
   "source": [
    "IndexType = itk.Index[3]\n",
    "index = IndexType()\n",
    "print(index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Define and use smart pointer object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image (0x555bfce47790)\n",
      "  RTTI typeinfo:   itk::Image<float, 2u>\n",
      "  Reference Count: 1\n",
      "  Modified Time: 1031\n",
      "  Debug: Off\n",
      "  Object Name: \n",
      "  Observers: \n",
      "    none\n",
      "  Source: (none)\n",
      "  Source output name: (none)\n",
      "  Release Data: Off\n",
      "  Data Released: False\n",
      "  Global Release Data: Off\n",
      "  PipelineMTime: 0\n",
      "  UpdateMTime: 0\n",
      "  RealTimeStamp: 0 seconds \n",
      "  LargestPossibleRegion: \n",
      "    Dimension: 2\n",
      "    Index: [0, 0]\n",
      "    Size: [0, 0]\n",
      "  BufferedRegion: \n",
      "    Dimension: 2\n",
      "    Index: [0, 0]\n",
      "    Size: [0, 0]\n",
      "  RequestedRegion: \n",
      "    Dimension: 2\n",
      "    Index: [0, 0]\n",
      "    Size: [0, 0]\n",
      "  Spacing: [1, 1]\n",
      "  Origin: [0, 0]\n",
      "  Direction: \n",
      "1 0\n",
      "0 1\n",
      "\n",
      "  IndexToPointMatrix: \n",
      "1 0\n",
      "0 1\n",
      "\n",
      "  PointToIndexMatrix: \n",
      "1 0\n",
      "0 1\n",
      "\n",
      "  Inverse Direction: \n",
      "1 0\n",
      "0 1\n",
      "\n",
      "  PixelContainer: \n",
      "    ImportImageContainer (0x555bfcf22c60)\n",
      "      RTTI typeinfo:   itk::ImportImageContainer<unsigned long, float>\n",
      "      Reference Count: 1\n",
      "      Modified Time: 1032\n",
      "      Debug: Off\n",
      "      Object Name: \n",
      "      Observers: \n",
      "        none\n",
      "      Pointer: 0\n",
      "      Container manages memory: true\n",
      "      Size: 0\n",
      "      Capacity: 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ImageType = itk.Image[itk.ctype('float'), 2]\n",
    "my_image = ImageType.New()\n",
    "print(my_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Print list of available types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<itkTemplate itk::Image>\n",
      "Options:\n",
      "  [<itkCType bool>, 2]\n",
      "  [<itkCType bool>, 3]\n",
      "  [<itkCType double>, 2]\n",
      "  [<itkCType double>, 3]\n",
      "  [<itkCType float>, 2]\n",
      "  [<itkCType float>, 3]\n",
      "  [<itkCType signed short>, 2]\n",
      "  [<itkCType signed short>, 3]\n",
      "  [<itkCType unsigned char>, 2]\n",
      "  [<itkCType unsigned char>, 3]\n",
      "  [<itkCType unsigned long long>, 2]\n",
      "  [<itkCType unsigned long long>, 3]\n",
      "  [<itkCType unsigned long>, 2]\n",
      "  [<itkCType unsigned long>, 3]\n",
      "  [<itkCType unsigned short>, 2]\n",
      "  [<itkCType unsigned short>, 3]\n",
      "  [<class 'itkCoocurrenceTextureFeaturesImageFilterPython.itkVectorF8'>, 2]\n",
      "  [<class 'itkCoocurrenceTextureFeaturesImageFilterPython.itkVectorF8'>, 3]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorD2'>, 2]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorD2'>, 3]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorD3'>, 2]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorD3'>, 3]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorD4'>, 2]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorD4'>, 3]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorF2'>, 2]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorF2'>, 3]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorF3'>, 2]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorF3'>, 3]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorF4'>, 2]\n",
      "  [<class 'itkCovariantVectorPython.itkCovariantVectorF4'>, 3]\n",
      "  [<class 'itkFixedArrayPython.itkFixedArrayF2'>, 2]\n",
      "  [<class 'itkFixedArrayPython.itkFixedArrayF3'>, 3]\n",
      "  [<class 'itkOffsetPython.itkOffset2'>, 2]\n",
      "  [<class 'itkOffsetPython.itkOffset3'>, 3]\n",
      "  [<class 'itkRGBAPixelPython.itkRGBAPixelUC'>, 2]\n",
      "  [<class 'itkRGBAPixelPython.itkRGBAPixelUC'>, 3]\n",
      "  [<class 'itkRGBPixelPython.itkRGBPixelUC'>, 2]\n",
      "  [<class 'itkRGBPixelPython.itkRGBPixelUC'>, 3]\n",
      "  [<class 'itkRunLengthTextureFeaturesImageFilterPython.itkVectorF10'>, 2]\n",
      "  [<class 'itkRunLengthTextureFeaturesImageFilterPython.itkVectorF10'>, 3]\n",
      "  [<class 'itkSparseFieldFourthOrderLevelSetImageFilterPython.itkNormalBandNodeIF2'>, 2]\n",
      "  [<class 'itkSparseFieldFourthOrderLevelSetImageFilterPython.itkNormalBandNodeIF3'>, 3]\n",
      "  [<class 'itkSymmetricSecondRankTensorPython.itkSymmetricSecondRankTensorD2'>, 2]\n",
      "  [<class 'itkSymmetricSecondRankTensorPython.itkSymmetricSecondRankTensorD3'>, 3]\n",
      "  [<class 'itkVectorPython.itkVectorD2'>, 2]\n",
      "  [<class 'itkVectorPython.itkVectorD2'>, 3]\n",
      "  [<class 'itkVectorPython.itkVectorD2'>, 4]\n",
      "  [<class 'itkVectorPython.itkVectorD3'>, 2]\n",
      "  [<class 'itkVectorPython.itkVectorD3'>, 3]\n",
      "  [<class 'itkVectorPython.itkVectorD3'>, 4]\n",
      "  [<class 'itkVectorPython.itkVectorD4'>, 2]\n",
      "  [<class 'itkVectorPython.itkVectorD4'>, 3]\n",
      "  [<class 'itkVectorPython.itkVectorD4'>, 4]\n",
      "  [<class 'itkVectorPython.itkVectorF1'>, 2]\n",
      "  [<class 'itkVectorPython.itkVectorF1'>, 3]\n",
      "  [<class 'itkVectorPython.itkVectorF2'>, 2]\n",
      "  [<class 'itkVectorPython.itkVectorF2'>, 3]\n",
      "  [<class 'itkVectorPython.itkVectorF3'>, 2]\n",
      "  [<class 'itkVectorPython.itkVectorF3'>, 3]\n",
      "  [<class 'itkVectorPython.itkVectorF4'>, 2]\n",
      "  [<class 'itkVectorPython.itkVectorF4'>, 3]\n",
      "  [<class 'stdcomplexPython.stdcomplexF'>, 2]\n",
      "  [<class 'stdcomplexPython.stdcomplexF'>, 3]\n"
     ]
    }
   ],
   "source": [
    "itk.Image.GetTypes()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Non-templated types\n",
    "\n",
    "* `MetaDataDictionary`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'new_key' value: 5\n"
     ]
    }
   ],
   "source": [
    "d = itk.MetaDataDictionary()\n",
    "d['new_key'] = 5\n",
    "print(\"'new_key' value: %d\" % d['new_key'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Note for C++ developers\n",
    "* If the type required using a smart pointer in C++\n",
    " * Use the method `New()` to create the object.\n",
    " * Objects are allocated on the heap.\n",
    " * C++ class has `Pointer` type and `New()` method defined ([see doxygen documentation](https://itk.org/Insight/Doxygen/html/index.html)) \n",
    "* If it was not a smart pointer (e.g. `itk.Index`)\n",
    " * use `()`\n",
    " * Allocated on the stack."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Filters with object oriented syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3e7cad0c9e0941eb8d8a8f947078e693",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Vmlld2VyKGdyYWRpZW50X29wYWNpdHk9MC4yMiwgcmVuZGVyZWRfaW1hZ2U9PGl0a0ltYWdlUHl0aG9uLml0a0ltYWdlRjI7IHByb3h5IG9mIDxTd2lnIE9iamVjdCBvZiB0eXBlICdpdGtJbWHigKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "PixelType = itk.ctype('float')\n",
    "image = itk.imread(\"data/CBCT-TextureInput.png\", PixelType)\n",
    "ImageType = itk.Image[PixelType, 2]\n",
    "median_filter = itk.MedianImageFilter[ImageType, ImageType].New()\n",
    "median_filter.SetInput(image)\n",
    "median_filter.SetRadius(4)\n",
    "median_filter.Update()\n",
    "view(median_filter.GetOutput())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Object exercises\n",
    "\n",
    "* In the example above, change the radius of the filter and observe the result.\n",
    "* Replace filter with `MeanImageFilter`\n",
    "* Replace filter with `OtsuThresholdImageFilter`\n",
    "* Visualize results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Uncomment and change the radius of the filter and observe the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# median_filter = itk.MedianImageFilter[ImageType, ImageType].New()\n",
    "# median_filter.SetInput(image)\n",
    "# median_filter.SetRadius(XX)\n",
    "# median_filter.Update()\n",
    "# median_filtered_image = median_filter.GetOutput()\n",
    "# view(median_filtered_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/4_Using_ITK_in_Python_real_world_filters_MedianFilter.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Uncomment and edit to use `MeanImageFilter`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# mean_filter = itk.XX[ImageType, ImageType].New()\n",
    "# mean_filter.SetInput(XX)\n",
    "# mean_filter.SetRadius(XX)\n",
    "# mean_filter.Update()\n",
    "# mean_filtered_image = mean_filter.GetOutput()\n",
    "# view(mean_filtered_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/4_Using_ITK_in_Python_real_world_filters_MeanFilter.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Uncomment and replace filter with `OtsuThresholdImageFilter`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# InputImageType = itk.Image[itk.ctype('float'), 2]\n",
    "# OutputImageType = itk.Image[itk.ctype('short'), 2]\n",
    "\n",
    "# otsu_filter = itk.OtsuThresholdImageFilter[XX]\n",
    "# XX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/4_Using_ITK_in_Python_real_world_filters_OtsuFilter.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## ITK object summary\n",
    "* Has `New()` method?\n",
    "* Is templated?\n",
    "* Call `Update()` with filters!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Automatic type deduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Often the filter type can automatically be deduced from the input of the filter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "383d5e5380f6412196dd1cf3982d9b08",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Vmlld2VyKGdyYWRpZW50X29wYWNpdHk9MC4yMiwgcmVuZGVyZWRfaW1hZ2U9PGl0a0ltYWdlUHl0aG9uLml0a0ltYWdlRjI7IHByb3h5IG9mIDxTd2lnIE9iamVjdCBvZiB0eXBlICdpdGtJbWHigKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "image = itk.imread(\"data/KitwareITK.jpg\", itk.F)\n",
    "# Filter types deduced from the image type.\n",
    "median_filter = itk.MedianImageFilter.New(image, radius = 3)\n",
    "median_filter.Update()\n",
    "view(median_filter.GetOutput())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Note the absence of image type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Type deduction exercises\n",
    "\n",
    "* Using the previous `MeanImageFilter` exercise, use the type deduction method to instantiate filter.\n",
    "* Using the previous `OtsuThresholdImageFilter` exercise, use the type deduction method to instantiate filter.\n",
    "* What is the image type of the `OtsuThresholdImageFilter` output?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Use the type deduction method to instantiate filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# mean_filter = itk. XX\n",
    "# mean_filter.Update()\n",
    "# mean_filtered_image = mean_filter.GetOutput()\n",
    "# view(mean_filtered_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/4_Using_ITK_in_Python_real_world_filters_MeanFilter-TypeDeduction.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Using the previous `OtsuThresholdImageFilter` exercise, use the type deduction method to instantiate filter.\n",
    "* What is the image type of the `OtsuThresholdImageFilter` output?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# otsu_filter = itk.XX\n",
    "# otsu_filter.XX\n",
    "# otsu_filtered_image = XX\n",
    "# print(otsu_filtered_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/4_Using_ITK_in_Python_real_world_filters_OtsuFilter-TypeDeduction.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simplify your code (and your life)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Objects accessible by index (non exhaustive)\n",
    "* MetaDataDictionary\n",
    "* Index\n",
    "* Point\n",
    "* Vector\n",
    "* Size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "index value at position 0: 1\n"
     ]
    }
   ],
   "source": [
    "# Create an `Index` of dimension 3\n",
    "index = itk.Index[3]()\n",
    "index[0]=1\n",
    "print(\"index value at position 0: %d\"%index[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Matrix is not, but you can get values with `()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Value at coordinate 0,0: 0.000000\n"
     ]
    }
   ],
   "source": [
    "# Create a 3x3 `Matrix` of type float\n",
    "m = itk.Matrix[itk.F,3,3]()\n",
    "print(\"Value at coordinate 0,0: %f\" % m(0,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* `itk.VariableLengthVector[itk.D]()` is neither.\n",
    "* Need to set and get values using the `SetElement` and `GetElement` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Value at coordinate 0,0: 1.000000\n"
     ]
    }
   ],
   "source": [
    "# Create a `VariableLengthVector of type `double` and initialize size to `1`\n",
    "v = itk.VariableLengthVector[itk.D](1)\n",
    "v.SetElement(0,1)\n",
    "print(\"Value at coordinate 0,0: %f\" % v.GetElement(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Use Python objects instead of ITK objects\n",
    "\n",
    "* Some filters expect inputs of a certain type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method itkImageBase2_SetOrigin in module _itkImagePython:\n",
      "\n",
      "itkImageBase2_SetOrigin(...) method of itkImagePython.itkImageF2 instance\n",
      "    SetOrigin(itkPointD2 _arg)\n",
      "    SetOrigin(double const * origin)\n",
      "    itkImageBase2_SetOrigin(itkImageBase2 self, float const * origin)\n",
      "    \n",
      "    Set the origin of the\n",
      "    image. The origin is the geometric coordinates of the image origin\n",
      "    (pixel [0,0]). It is stored internally as SpacePrecisionType but may\n",
      "    be set from float or double. See:   GetOrigin()\n",
      "\n"
     ]
    }
   ],
   "source": [
    "image = ImageType.New()\n",
    "help(image.SetOrigin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Use Python `list` instead (for `Index`, `Point`, `Vector`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image origin: itkPointD2 ([2, 10])\n"
     ]
    }
   ],
   "source": [
    "image.SetOrigin([2,10])\n",
    "print(\"Image origin: %s\" % str(image.GetOrigin()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Functions to know\n",
    "* `itk.imread(file_name [, pixel_type])`\n",
    "* `itk.imwrite(image, file_name [, compression])`\n",
    "* `itk.array_from_image(image)` and `itk.array_view_from_image(image)`\n",
    "* `itk.image_from_array(arr)` and `itk.image_view_from_array(arr)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Pixel types - Two options\n",
    "* itk.ctype('float'), itk.ctype('unsigned char')\n",
    "* itk.F, itk.UC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Convenience functions\n",
    "* `itk.size(image)`\n",
    "* `itk.spacing(image)`\n",
    "* `itk.origin(image)`\n",
    "* `itk.index(image)`\n",
    "* `itk.physical_size(image)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fixing errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Unsupported (image) types\n",
    "* typical error message:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "KeyError: \"itkTemplate : No template [<class 'itkImagePython.itkImageD2'>] for the itk::ImageFileReader class\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# image = itk.imread(\"data/BrainProtonDensitySlice.png\", itk.D)\n",
    "# print(image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Select different input image type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "image = itk.imread(\"data/KitwareITK.jpg\", itk.F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Cast to a different image type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "InputImageType = itk.Image[itk.F, 2]\n",
    "OutputImageType = itk.Image[itk.UC, 2]\n",
    "cast_filter = itk.CastImageFilter[InputImageType, OutputImageType].New(image)\n",
    "cast_filter.Update()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* Converting a NumPy array with too many components along one dimension:\n",
    " * Use `itk.VectorImage`\n",
    " * Save it using `itk.imwrite()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Create your own pipeline\n",
    "\n",
    "* Why?\n",
    " * Create a pipeline in a Python function.\n",
    " * Return the last filter output without updating it.\n",
    " * Error: Not be able to run the pipeline outside of the function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* This does not work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-56-60e88659c137>, line 3)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-56-60e88659c137>\"\u001b[0;36m, line \u001b[0;32m3\u001b[0m\n\u001b[0;31m    ...\u001b[0m\n\u001b[0m    ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "  def generateSource():\n",
    "    reader = itk.ImageReader[ImageType].New()\n",
    "    ...\n",
    "    itkToVtkFilter = itk.ImageToVTKImageFilter[ImageType].New()\n",
    "    itkToVtkFilter.SetInput(reader.GetOutput())\n",
    "    return itkToVtkFilter.GetImporter()\n",
    " \n",
    "  def main():\n",
    "    source = generateSource()\n",
    "    source.Update()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* Instead, you can create a pipeline, and use it where convenient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "24c75f05fcbd4024a490876c051534ae",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Vmlld2VyKGdyYWRpZW50X29wYWNpdHk9MC4yMiwgcmVuZGVyZWRfaW1hZ2U9PGl0a0ltYWdlUHl0aG9uLml0a0ltYWdlVUMyOyBwcm94eSBvZiA8U3dpZyBPYmplY3Qgb2YgdHlwZSAnaXRrSW3igKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def my_func(ImageType):\n",
    "    my_pipeline = itk.pipeline()\n",
    "    mean_filter = itk.MeanImageFilter[ImageType, ImageType].New()\n",
    "    my_pipeline.connect(mean_filter)\n",
    "    my_pipeline.expose(\"Radius\")\n",
    "    threshold_filter = itk.ThresholdImageFilter[ImageType].New()\n",
    "    my_pipeline.connect(threshold_filter)\n",
    "    my_pipeline.expose(\"Lower\")\n",
    "    return my_pipeline\n",
    "\n",
    "PixelType = itk.UC\n",
    "ImageType = itk.Image[PixelType, 2]\n",
    "image = itk.imread(\"data/CBCT-TextureInput.png\", PixelType)\n",
    "\n",
    "my_pipeline = my_func(ImageType)\n",
    "my_pipeline.SetInput(image)\n",
    "my_pipeline.SetRadius(5)\n",
    "my_pipeline.SetLower(100)\n",
    "my_pipeline.Update()\n",
    "\n",
    "view(my_pipeline.GetOutput())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Very Advanced topic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* disable lazy loading:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'absolute_value_difference_image_filter'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-58-13aabc83a193>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mitk\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mforce_load\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/home/francois.budin/.virtualenvs/itk_course/local/lib/python2.7/site-packages/itkExtras.pyc\u001b[0m in \u001b[0;36mforce_load\u001b[0;34m()\u001b[0m\n\u001b[1;32m    113\u001b[0m     \u001b[0;32mimport\u001b[0m \u001b[0mitk\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    114\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdir\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 115\u001b[0;31m         \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    116\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    117\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/francois.budin/.virtualenvs/itk_course/local/lib/python2.7/site-packages/itkLazy.pyc\u001b[0m in \u001b[0;36m__getattribute__\u001b[0;34m(self, attr)\u001b[0m\n\u001b[1;32m     45\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mnamespace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     46\u001b[0m                 \u001b[0msetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 47\u001b[0;31m             \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnamespace\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mattr\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     48\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'absolute_value_difference_image_filter'"
     ]
    }
   ],
   "source": [
    "itk.force_load()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* call_back while importing ITK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import itk\n",
    "itk.auto_progress(True)\n",
    "itk.force_load()\n",
    "itk.auto_progress(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* Cast and downcast ITK objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "ImageType = itk.Image[itk.F, 2]\n",
    "filter = itk.MedianImageFilter[ImageType, ImageType].New()\n",
    "casted_filter = itk.ProcessObject.cast(filter)\n",
    "itk.down_cast(casted_filter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# NumPy and scikit-learn interactions\n",
    "\n",
    "* ITK is great at reading and processing images\n",
    "* ITK does not show images\n",
    "* Some algorithms are not available in ITK\n",
    "* NumPy is great at processing arrays in simple ways\n",
    "* NumPy arrays can be read by many other Python packages\n",
    "    * matplotlib\n",
    "    * scikit-learn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's convert an ITK image to a NumPy array\n",
    "\n",
    "* Load an image with ITK\n",
    "* Convert the image to a NumPy array\n",
    "* Check some values in the image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import itk\n",
    "from itkwidgets import view\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[51 56 60]\n"
     ]
    }
   ],
   "source": [
    "im=itk.imread(\"data/KitwareITK.jpg\")\n",
    "arr=itk.array_from_image(im)\n",
    "print(arr[1,1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let go the other way around: NumPy array to an ITK image. First, we create an array with some values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def make_gaussian(size, fwhm = 3, center=None):\n",
    "    \"\"\" Make a square gaussian kernel.\n",
    "\n",
    "    size is the length of a side of the square\n",
    "    fwhm is full-width-half-maximum, which\n",
    "    can be thought of as an effective radius.\n",
    "    \"\"\"\n",
    "\n",
    "    x = np.arange(0, size, 1, np.float32)\n",
    "    y = x[:,np.newaxis]\n",
    "\n",
    "    if center is None:\n",
    "        x0 = y0 = size // 2\n",
    "    else:\n",
    "        x0 = center[0]\n",
    "        y0 = center[1]\n",
    "\n",
    "    return np.exp(-4*np.log(2) * ((x-x0)**2 + (y-y0)**2) / fwhm**2)\n",
    "\n",
    "\n",
    "arr = make_gaussian(11)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's look at the array. We use matplotlib to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5, 10.5, 10.5, -0.5)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD8CAYAAABXXhlaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAABFJJREFUeJzt3cFpGwsUhtGZh0BrowrUgltwCy7TPaQDt6AKhHYGgcy8bRYBh3DtxPrOWZtfg+yP2YSbddu2BWj5728/APD1hA9Bwocg4UOQ8CFI+BAkfAgSPgQJH4J2X/lh67r6Z4LwybZtWz/6GW98CBI+BAkfgoQPQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgQJH4K+9OYef25dPzyjdlf8L86fyxsfgoQPQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgQJH4KED0HChyA39z7Jbjf71e73+9G96ee73W6je9frdXRv+vm+O298CBI+BAkfgoQPQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgQ5vfWTyXNUDw8PY1vLsizH43F073A4jO6dz+fRvdPpNLp3uVzGtu7hjJc3PgQJH4KED0HChyDhQ5DwIUj4ECR8CBI+BAkfgoQPQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9B3/rm3rquo3v7/X5sa/pG3vPz8+je4+Pj6N7r6+vo3svLy+je29vb2Nb7+/vY1rIsy7Zto3u/wxsfgoQPQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgQJH4KED0HCh6BvfXNv2m4393UcDoexrWWZv5H39PQ0ujftx48fo3uTv9t74I0PQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgQJH4KED0HChyDhQ5BDZD+53W5jW+fzeWxrWZbl9fV1dG/a9PNNf3+Tv9t74I0PQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgQJH4KED0HChyDhQ9C6bdvXfdi6ft2H/YHdbu4E4cPDw9jWsizL8Xgc3TscDqN70zfyTqfT6N7lchnb+tfv923btn70M974ECR8CBI+BAkfgoQPQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgS5ufdJJu/3Lcuy7Pf70b3p55u+Q3e9Xkf3/vU7eZPc3AN+SfgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgQJH4KED0HChyDhQ5DwIUj4EOT01jexrh9eU7orX/l3eW+c3gJ+SfgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAgSPgQJH4KED0HChyDhQ5DwIUj4ECR8CNr97Qfg97hBxyRvfAgSPgQJH4KED0HChyDhQ5DwIUj4ECR8CBI+BAkfgoQPQcKHIOFDkPAhSPgQJHwIEj4ECR+CVrfcoMcbH4KED0HChyDhQ5DwIUj4ECR8CBI+BAkfgoQPQcKHIOFDkPAhSPgQJHwIEj4ECR+ChA9Bwocg4UOQ8CFI+BAkfAj6H6cyd3otC1Z1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.gray()\n",
    "plt.imshow(arr)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "And now, let's convert it to an ITK image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "im = itk.image_view_from_array(arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can now look at pixel values..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "im.GetPixel([5,5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    " or call any ITK filter on the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "filtered_image=itk.median_image_filter(im, radius=[2,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can compare some pixel values to see the effect of the filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original image pixel value: 1.000000\n",
      "Filtered image pixel value: 0.291632\n"
     ]
    }
   ],
   "source": [
    "print(\"Original image pixel value: %f\"%im.GetPixel([5,5]))\n",
    "print(\"Filtered image pixel value: %f\"%filtered_image.GetPixel([5,5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's modify a pixel in our image, and look at the original array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7efe2f117110>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD8CAYAAABXXhlaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAC2RJREFUeJzt3X+o3Xd9x/Hna7k22jo0iyAz6db7hziKMFqCVAsyGgfdFLs/yqhQcTLIP1OrCFL3z/7dHyL6xxBCrRMsLRLLVsqwilbG/ilNbwu2iW6lcW1qa7ssU7GQmvreH/e43YS0N5zz+d5z0vfzASX3nJ58zpt787zf8/NzUlVI6uV3lj2ApJ1n+FJDhi81ZPhSQ4YvNWT4UkOGLzVk+FJDhi81tLaTV5bElwlKE6uqbHcZj/hSQ4YvNWT4UkOGLzVk+FJDhi81tFD4SW5M8uMkTya5fdRQkqaVeXfgSbIL+HfgT4GTwMPAR6rq2Gv8HZ/HlyY29fP47wGerKqnqupl4B7gpgXWk7RDFgl/H/DMltMnZ+edI8mhJEeTHF3guiQNNPlLdqvqMHAYvKkvrYpFjvjPAlduOb1/dp6kFbdI+A8D70yynuQy4BbgvjFjSZrS3Df1q+pskk8ADwC7gDur6olhk0mazNxP5811Zd7Hlybn23IlXZDhSw0ZvtSQ4UsN7eiee5pfsu3jNa8rforztDziSw0ZvtSQ4UsNGb7UkOFLDRm+1JDhSw0ZvtSQ4UsNGb7UkOFLDRm+1JDhSw0ZvtSQ4UsNGb7UkOFLDRm+1JDhSw25595E1tbGfmt37949dL3R8509e3boemfOnBm63uj5LnUe8aWGDF9qyPClhgxfasjwpYbmDj/JlUkeTHIsyRNJbhs5mKTpLPKczlngs1W1keR3gUeSfLeqjg2aTdJE5j7iV9VzVbUx+/qXwHFg36jBJE1nyH38JFcB1wAPjVhP0rQWfvlWkjcD3wI+XVW/uMD/PwQcWvR6JI2TRT6OOMkbgPuBB6rqixdx+TaffexLdhfjS3bnV1Xbfqb6Io/qB/gqcPxiope0Oha5j3898FHghiSPzf7780FzSZrQ3Lf3qurfgG1vUkhaPb5yT2rI8KWGDF9qyPClhtx6a4uRz23v2bNn2FoA6+vrQ9fbu3fv0PVOnTo1dL0TJ04MXe/06dPD1no9vCbAI77UkOFLDRm+1JDhSw0ZvtSQ4UsNGb7UkOFLDRm+1JDhSw0ZvtSQ4UsNGb7UkOFLDRm+1JDhSw0ZvtSQ4UsNGb7U0CW9597mp3iNM/Lz6UbvkXfzzTcPXe/aa68dut7GxsbQ9Y4cOTJ0vZdeemnYWq+88sqwtQAW+fzKeXnElxoyfKkhw5caMnypIcOXGjJ8qaGFw0+yK8mjSe4fMZCk6Y044t8GHB+wjqQdslD4SfYDHwTuGDOOpJ2w6BH/S8DngN+82gWSHEpyNMnRBa9L0iBzh5/kQ8ALVfXIa12uqg5X1YGqOjDvdUkaa5Ej/vXAh5P8BLgHuCHJN4ZMJWlSc4dfVZ+vqv1VdRVwC/D9qrp12GSSJuPz+FJDQ96WW1U/AH4wYi1J0/OILzVk+FJDhi81ZPhSQ5f0nnujra2N+3bs3bt32Fowfo+8gwcPDl1vtAcffHDoeiN/tq8HHvGlhgxfasjwpYYMX2rI8KWGDF9qyPClhgxfasjwpYYMX2rI8KWGDF9qyPClhgxfasjwpYYMX2rI8KWGDF9qyPClhtyIbIuzZ88OW+vUqVPD1gLY2NgYut5oo+cb/f0b+bN9PfCILzVk+FJDhi81ZPhSQ4YvNbRQ+EnemuRIkh8lOZ7kvaMGkzSdRZ/O+zLw7aq6OcllwOUDZpI0sbnDT/IW4P3AXwFU1cvAy2PGkjSlRW7qrwMvAl9L8miSO5JcMWguSRNaJPw14FrgK1V1DfAr4PbzL5TkUJKjSY4ucF2SBlok/JPAyap6aHb6CJu/CM5RVYer6kBVHVjguiQNNHf4VfU88EySd83OOggcGzKVpEkt+qj+J4G7Zo/oPwV8fPGRJE1tofCr6jHAm/DSJcZX7kkNGb7UkOFLDRm+1JDhSw2lqnbuypKdu7I5rK2N24Jwz549w9YCWF9fH7re3r17h643eo+8EydODF3v9OnTw9Za9f37qirbXcYjvtSQ4UsNGb7UkOFLDRm+1JDhSw0ZvtSQ4UsNGb7UkOFLDRm+1JDhSw0ZvtSQ4UsNGb7UkOFLDRm+1JDhSw0ZvtSQe+5NZOT+fQC7d+8eut7o+UbvQ3fmzJmh6636PnkjueeepAsyfKkhw5caMnypIcOXGloo/CSfSfJEkseT3J3kjaMGkzSducNPsg/4FHCgqt4N7AJuGTWYpOkselN/DXhTkjXgcuCni48kaWpzh19VzwJfAJ4GngN+XlXfOf9ySQ4lOZrk6PxjShppkZv6e4CbgHXgHcAVSW49/3JVdbiqDlTVgfnHlDTSIjf1PwCcqKoXq+rXwL3A+8aMJWlKi4T/NHBdksuTBDgIHB8zlqQpLXIf/yHgCLAB/HC21uFBc0makO/Om4jvzluM786bn+/Ok3RBhi81ZPhSQ4YvNeSDe5eIzWdM+9jJf5evNz64J+mCDF9qyPClhgxfasjwpYYMX2rI8KWGDF9qyPClhgxfasjwpYYMX2rI8KWGDF9qyPClhgxfasjwpYYMX2rI8KWGxn6qgibjHnQaySO+1JDhSw0ZvtSQ4UsNGb7U0LbhJ7kzyQtJHt9y3u8l+W6S/5j9uWfaMSWNdDFH/H8EbjzvvNuB71XVO4HvzU5LukRsG35V/Svw3+edfRPw9dnXXwf+YvBckiY07338t1fVc7OvnwfePmgeSTtg4VfuVVW91qfgJjkEHFr0eiSNM+8R/2dJfh9g9ucLr3bBqjpcVQeq6sCc1yVpsHnDvw/42OzrjwH/PGYcSTsh2735I8ndwJ8AbwN+Bvwd8E/AN4E/AP4T+MuqOv8BwAut5TtNpIlVVba7zLbhj2T40vQuJnxfuSc1ZPhSQ4YvNWT4UkOGLzW003vu/RebT/9t522zy66iVZ4NVnu+VZ4NVnu+i53tDy9msR19Ou9iJTm6qq/0W+XZYLXnW+XZYLXnGz2bN/WlhgxfamhVwz+87AFewyrPBqs93yrPBqs939DZVvI+vqRpreoRX9KEVir8JDcm+XGSJ5Os1D5+Sa5M8mCSY0meSHLbsmc6X5JdSR5Ncv+yZzlfkrcmOZLkR0mOJ3nvsmf6rSSfmf1MH09yd5I3LnmeyTe4XZnwk+wC/gH4M+Bq4CNJrl7uVOc4C3y2qq4GrgP+ZsXmA7gNOL7sIV7Fl4FvV9UfAX/MisyZZB/wKeBAVb0b2AXcstyppt/gdmXCB94DPFlVT1XVy8A9bG7quRKq6rmq2ph9/Us2/+HuW+5U/y/JfuCDwB3LnuV8Sd4CvB/4KkBVvVxV/7Pcqc6xBrwpyRpwOfDTZQ6zExvcrlL4+4Bntpw+yQqFtVWSq4BrgIeWO8k5vgR8DvjNsge5gHXgReBrs7sidyS5YtlDAVTVs8AXgKeB54CfV9V3ljvVBQ3d4HaVwr8kJHkz8C3g01X1i2XPA5DkQ8ALVfXIsmd5FWvAtcBXquoa4FesyGcxzO4r38TmL6d3AFckuXW5U7222nwqbqGn41Yp/GeBK7ec3j87b2UkeQOb0d9VVfcue54trgc+nOQnbN5FuiHJN5Y70jlOAier6re3kI6w+YtgFXwAOFFVL1bVr4F7gfcteaYLuegNbi/GKoX/MPDOJOtJLmPzAZb7ljzT/0kSNu+jHq+qLy57nq2q6vNVtb+qrmLz+/b9qlqZo1ZVPQ88k+Rds7MOAseWONJWTwPXJbl89jM+yIo88HieoRvc7vS7815VVZ1N8gngATYfWb2zqp5Y8lhbXQ98FPhhksdm5/1tVf3LEme6lHwSuGv2S/0p4ONLngeAqnooyRFgg81nbh5lya/g27rBbZKTbG5w+/fAN5P8NbMNbhe6Dl+5J/WzSjf1Je0Qw5caMnypIcOXGjJ8qSHDlxoyfKkhw5ca+l8Y4tueAcRfjQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "im.SetPixel([5,5], 0)\n",
    "plt.imshow(arr)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now let's load an image and show it with matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7efe2f108d90>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADiCAYAAABA4gwKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzsvWtsXGea3/k7db/fbyxeRZHUzbe22y2P7Wl3p92ddiMTpzFBIzPAzvTsYHo/ZLBYIB/Smy9ZYBFgPuxmESBAgF4kSAbYZDKYmc40GsmOez0xuttuy23ZlmXRlESRoihei8ViFet+O/uBel6fokiJkkhRlN4fIJAsnqpzqiT93+c87/P8H8M0TTQajUbz+GI77AvQaDQazcGihV6j0Wgec7TQazQazWOOFnqNRqN5zNFCr9FoNI85Wug1Go3mMefAhN4wjG8bhnHZMIxpwzB+eFDn0Wg0Gs2dMQ6ijt4wDDtwBfgmcBP4NfA7pmlO7vvJNBqNRnNHDiqi/wowbZrmjGmaTeDPgDcP6FwajUajuQOOA3rdfmDe8vNN4OxuBxuGodtzNRqN5t5ZM00zebeDDkro74phGD8AfnBY59doNJrHgLm9HHRQQr8ADFp+Hrj1mMI0zR8BPwId0WseX2y2reyoYRjqT6fTAcC6P6Y9pzQHyUEJ/a+BccMwjrEl8P8I+N0DOpdG88jS7XbveoxhGA/hSjRPMgci9KZptg3D+GPgbwA78O9M07x0EOfSaB5lfD4fNpsNv99PIBDA7XbT7Xap1+tsbGxQr9ep1+uHfZmax5wDKa+854vQqRvNQ0RSKKZp9qRMDMNQqRZJrwhOp5NWq3Xba9lsNoLBIMlkEo/HQyQSYXR0lGg0ytDQEMPDwzidTsLhMD6fD8MwcLlcmKZJpVKhWq1SLBb57LPPWFtbY35+nvn5eebm5qjValSrVXUe+OIOQa5f88Rz3jTNL9/tIC30micWEXZrekXSKHa7nXa73SOmTqcTr9fL0NAQg4ODDA0NEY1GyWQyDA8P4/f7GRgYIJvN4na7cbvdt6VlZHExTVP9zrq4NJtN1tbWyOVyVKtVrl27xi9/+Uvm5ua4fPkyuVxO3QHY7Xb1vEfh/7HmUNBCr9Fsx+Px0G63lbjbbDYlktb/C1bxt9vthEIhBgcHicfjnDlzhtHRUex2O4ODgxw7doxsNovH48Hv92O329VryZ2B9Q7C6XT2RObW6FzEv1arYRgGjUaDK1eusLCwwFtvvcXFixdZWlrixo0bPXcY2xcszRODFnqN5k6IMDscDiXCIpYOhwOv10sikSAYDDI0NMTY2BiJRILR0VGOHz/O0NAQTqcT+ELQrdG5YRgq6rYKe7fb7YnsRfwBtQh5PB71OtbIfXJyktnZWf7bf/tv/PrXv2ZmZoZisUi329Vi/2SihV7zZGOz2bDZbLTbbfVYJBLh+PHjJBIJstkskUhEHdNsNqlWqzSbTVwuF263m0gkQjAYJJPJkM1mSSQSTExMEIlEHnq6xJoGmpyc5G//9m+5cOECc3NznDt3jlKppO4QtOA/MexJ6A+tYUqjOSgkFdLtdpXgud1uQqEQ/f39HDt2jHg8Tn9/P8FgEJvNRqvVolarUS6XaTabwFZOPhgMEgwGCQQC6qvDcTj/bazRfV9fH2fOnKFQKGCaJgsLC1y7do1Go6Hz9Zrb0EKveWyw2Ww9DUmwJfrpdJqRkRH6+vqYmJhgZGSEWCxGKpXC7/djmqYqd1xbW6Ner1Or1TBNE5fLhdfrJRwOk0gkSCQSBAKBQxFTu91Ot9vFMAyi0SgvvfQSAMPDwySTST788ENWV1eZnp7uqdbR0b1GC73msUIE2OFwKFEeHR0lm80SDocZGBigv7+fcDhMLBbD6/XSaDSw2WxUq9WeDlapynE4HDQaDZrNZs9G7mE1OnW7Xex2O16vl9OnTxMKhbDZbDidTorFIvF4nI8++ohqtdqTttI8uWih1zw2WKPsUChEKpXC5/Ph9/vxer14vV7cbjcOhwOHw6HuAOS5kurpdDoqcrbb7djtdlwuF3a7vack8jCietnsFaSGf21tjY2NDfx+P5ubm+RyOVZWVsjlcg/9GjWPHlroNUcSu91+W1OTCG84HGZ0dJRkMqk2VMPhMMFgsCe/3u12qVQqSuAbjQaNRoN2u02n08HhcKjKHJ/PRyQSIRAIqPTQdtE9aHZaWLxeL4ODg1SrVer1Oq1WSy1G+Xyed999l3w+ryqL5H3pSP/JQgu95khiLWe0br76fD76+/tJJpNEIhFcLhc+n09F8DabjU6nQ6vVUhuXtVqNer1OpVJRzUc+nw+fz0cymSQWizE6OkosFlNdrY+KP41cR39/P5VKhVKpRLPZpNlskkwmaTabvPfee5RKJYDbqpA0TwZa6DVHlu0bjZFIhEwmw6lTp8hkMng8HpxOJ36/n0gkgt/vJxwO4/V6Vd26iLbdbicQCOB0OjFNE5vNRl9fH0NDQ0QiEXw+X0/NvETIh4107waDQfr7+/H5fFQqFeLxOD6fjxMnTuB2u3n//ffJ5XJ6Y/YJRQu95shijartdjuJRIK+vj5SqZTynrHb7fh8PoLBIF6vl2AwiM/nw+VyqQ7VVquF3W5X+XsAv99PX18fsVgMv9//yNoMWP16vF4vrVaLQCCgOnSj0SiVSoV0Ok2n02Ftbe2Qr1hzGGih1xxZrBuj2WyW1157jcHBQcbHx4nH4wSDQXWMmInJxqqkcrrdLm63m3a7rXL4NpuNWCxGNpsFtsTU4XD0bNy2221V036YWDeIZQHL5/Pk83lsNhuhUIh4PM7Zs2dZWVnhww8/ZHFx8ZCvWvOw0UKveeSx2giYpqnqydvttupw/a3f+i1effVVBgYGGB4exuPxqBSMeMZIRY00REmjlJRMihWBYRj4/f4eqwLrJq6kbURgpUJnJwOznXL5+3lnYF18pJoolUqRy+Ww2WzU63WcTic+n49oNIphGFy4cIGFhQW1cas7aR9/tNBrHnmsTpISWRuGQbvdZmxsjDNnzvDss89y6tQpYrEYwWDwttfw+XyYpkmz2VTVJ/CF6MpmrnjUWC2Lt2M1IWs0GtRqNRVZ+/3+HnHfyTDtoDdy/X4/fr+fer2OzWZTdzHdbpdYLEYmk6FarbK6uqo+Uy30jzda6DWPPNYI3Jorf+mll/je977HU089xUsvvaRKH615a1kgxDK41WpRKpXUwI96vU6n06HT6ah0jAh2MpkkHA7jdrtxuVw911Mul8nlcuRyOTY2NnA4HNjtdoLBIOl0Go/HQzgcxuPxqEVDFpiDTvkEAgGVk9/c3FTnFvsGsVAul8u6qeoJQQu95kggKZhut0sqlWJ0dJQ/+qM/4vXXXyeVSuHxeNQiYLfbVTQuAi2iNzMzw/LyMpubmxSLRarVqhL5VqulzhMIBMhkMqTTaTUdSkoTFxYWWF5eZmNjg0KhQKFQUNfocrmIxWL4fD4GBgYYGRnB6XQqF8yH4ZNjt9vp7+9XXb+rq6uqCSwcDlOv1zl9+jTFYpGLFy/e1o+gefzQQq85Ekg0D5BIJHj22WcZHR0lkUjg8XioVqsqeramRiQVIxH4zZs3WVxcpFQqUSgUqFartFotut0uzWZTCb3NZlPdpYFAQG10mqapHi8UCmxsbFAqlVTqIxgMUigU1DW1Wi3cbjf1el1VAsXj8QP/vDweD8FgEI/HQygUUukZp9NJrVaj3W4TiURwOBzqM3oUq4o0+4MWes0jj4zxs9vtnDhxgu9///ucOnWK1157TYm/x+NRUXy73Vbdn3a7nfn5eT755BNu3LjBpUuXuHr1KsVikXK5TKlUUukdSeHAlmBLXb7P5yMcDiuxLBQK5PN5FhYWqFarVKtVtUBYSzkHBgaYnp4mFAoxNjampk+tra0Rj8fVJqnU5+8XsgkcDAY5deoU8/PzrK+vq6je4/EQjUap1WoUCgVmZmbY3Nzs6TbebXSi5miihV7zyCOC09fXx5tvvsk3v/lNnnrqqZ4otNlsYrfbVXlku90ml8tRKBR45513eO+991hcXOTTTz9lY2OjJ12xfR4rQLFY5ObNm3z++efYbDa8Xq/KrcsGrDzPOrCkWCyytLQEwOXLlzl+/DjhcFh1rLpcLgKBgLJWECuF/czbyxASwzBUmenm5qby3BcvHLmLSafTvPfee5TLZbVAapF/vNBCrzkSeL1exsfHeeGFFzh+/LgSeRFI2YBttVq4XC6q1SpTU1MsLCzwwQcfcOHCBXK5HOvr62pxkI3XnXLUUmNv3ai1LgjWyh1rqsjqI1Ov11lcXKRWq5HNZtnY2KBYLKoI3ul09tTn7xfynuT9uVwuQqGQeh+y31AoFEilUtTrdRKJBOVyWefrH1MeSOgNw7gObAIdoG2a5pcNw4gB/xkYAa4D3zNNs/Bgl6l5knE6nUxMTPA7v/M7vPbaawQCAfU7EV+v1wtsiXAul+Nv/uZv+OUvf8n169f56KOPVOrCynahtpZTbh8Yvn1BkLmw22e9WitYDMOgUChQLBbV6/t8Pmw2G4ODg8qKodFo9FT1PChi6SApHOkPEHsHwzCUv77D4SCbzSpTtOXlZfUaOmf/+LAfEf3XTdO09lX/EHjbNM0/MQzjh7d+/qf7cB7NE4TP56PVatFqtXj99dd54YUXeOONN0gkEmrz1O12K+ESEV5aWuKdd97hv/yX/8Jf/uVf7ulcO9W6b//99t/tJfKVa+t2u9y4cYNIJMKJEyeo1+uqzFGaug4Cq1hbr1/uRmR4ydLSEqurq6yvr5PP53Xa5jHkIFI3bwJfu/X9fwDeQQu95h6RCUnPPfccb775Jk8//TQDAwMqzy3CJUM4TNNkcnKSyclJ/uqv/oq33377MC9fIWLbbreZm5tjcXGRYDCoxFSi7YcdPUsJqpifvfbaa1SrVZaXl1lcXGR9ff2hXo/mYHnQUMIE3jIM47xhGD+49VjaNM2lW98vA+mdnmgYxg8Mw/jQMIwPH/AaNI8xJ06c4MyZMwwMDABf+MxsT7d0Oh2uX7/O1atXmZubo16vH9o1WxF7BIByuUy5XFYbo4eJtfHM7XaTyWTo7++nr6+PUCh0qNem2X8eNKJ/1TTNBcMwUsDPDMOYsv7SNE3TMIwdQxXTNH8E/Ahgt2M0Ty5er5ehoSHeeOMNzp49q2wQ3G43rVZLpVOkvv7ixYu8//77fPLJJ0xNTalu2kcBSeFIs1U0GmVxcZHR0dFDieaBnpJOu91ONpvl1KlTLC4uEg6He/YWNEefB4roTdNcuPV1Ffgx8BVgxTCMPoBbX1cf9CI1TxaGYTA0NMTXvvY1vv71rytBl3JAMR4zDIN6vc7Nmzd5++23+eyzz5iZmaFare57bfr9sr156/r169y8eZPLly9TqVSA/TU52yv1er2ncsjlcnHixAmee+45XnjhBfr6+h76NWkOjvsWesMw/IZhBOV74FvAZ8BPgN+/ddjvA3/9oBepebwRUbZ6woyPjzMxMYHX61WGXE6nE5fLRaPRoF6vYxgGMzMzXL16lXK5rDYUgYe6oSjXvdOmqixSNpsNt9tNtVqlUCgonxyJ9iWXbx0FeJC43W7gi9m3TqeTdDpNX18f6XRaWTfAF2WoD8O+QXMwPMjfXBr48a1/lA7gP5qm+f8ahvFr4M8Nw/hDYA743oNfpuZxxeFwKFGWJp9MJsPLL7/M888/r+wCRJDq9ToejweAqakpzp07x9LSEp988gn5fF5t4h4kVt8dSSHtFpVb9xCktn1lZYWZmRk+/PBDqtUqzz77rBJ8WdQOY1Sh0+nk2LFjeL1ePv/8cz777DM2NjaoVCo9TWGao8d9C71pmjPAszs8nge+8SAXpXlykE1Ja8v9V77yFX7zN3+T8fHxHtdHKUm8ceMGN27c4Pz580xPT7O6usr09DRLS0sPReh3ErzdfN2tPzscDnw+n6rAOXfuHPl8nnQ6TTqdVs8/zHm0qVRKjSAcGhrCNE3K5fIjMyNXc3883DH2Gs0uSEQcCAQ4efKkmvlqnekKUKlUuHr1KlevXmV2dpa5uTlmZmZYW1tTeefDuv47RbyS+pA0T6VSYXV1lbm5OaamptSCt1un7sNArBM8Hg+JRIJ0Ok00GlXXpRuoji466aY5VNxuN41Gg3a7zdDQECdPnuT1119naGhIpUikFHBjY4N3332X999/n/n5eS5dusT8/DyFQgGXy6UE0mpDsN84HA48Hk+Pp069Xt+1yicUCuF0OlXnrhiHra6uYrPZlIVwtVrl+eefJ5VKHZqgymLkcDg4efIkzz//PDabjStXrqj9A83RRAu95lBpNBrAloD+/b//93n11Vc5e/YsLperJ8qt1WqcP3+ec+fO8c4773Djxg01+1RG5gn7KfLWaVOxWIxAIEAsFsPj8dBqtahUKpTL5Z67CRlCEovFiMfjdLtdut0ulUqFXC6HYRjKC9/hcCjnyJWVFf7BP/gHRKPRQxN72Ufo7+/n9OnT2Gw2Pv74Y2ZmZlTlk+booYVec6h4PB7q9Tp/9+/+Xb7//e/zzDPP4HQ6VVWNw+GgWq2ytLTEr371K372s59x4cIFFb1L1H/QyGzacDhMOBxWG7GlUglAbaK6XC7lA59Op0kmk3Q6HWq1GktLS2pRsPrgF4tFut0uCwsLhMNhvvvd7x74+9kJq0dOIBDgq1/9KsFgkMnJSZaXl9WirDl66By95lBptVp4PB5eeuklhoaGVEmfzWZT5XytVkvZ/66uripvekA5TB40YgomNfxSaSMe79ZI1+l04vF41LASj8eD1+vF4/HgdrvpdrtK5GHrrmZtbY1cLsfMzIxaPA4LcbkMhUJkMhkGBgYIh8Pa2fIIoyN6zaHS6XT4h//wH/J7v/d7JJNJ9bi4OZqmyfz8PJ9//jmffvop8/Pz2O322/LFB5mXl+uUBieJxlutFuvr6xSLRVUvHwqFiMViOBwOJe5yfX6/X3X2ikePLBLz8/NUq1WOHTvGzMwMzz333IG9lzu9R/H0l6EtZ86c4eWXX+batWv8+Mc/1mJ/RNFCr3loSJrGSjab5Q/+4A/o7+9XjzUaDSWITqeT2dlZ5ufnuXbtWk83p5WDEnmJ3HO5HKVSCYfDoTYmdyrllDp/0zRJpVJ0Oh1arRa1Wo1KpUKpVFKDSeR9SD5+c3OTUqlErVY7kPdyN6SUFb5okgJ4+eWX+fjjj/nkk0+YmZkBdi4x1Ty6aKHXPDS2V6a43W6efvppRkdHVSpGIkn4Ii1TLBbZ2Ng4NAEElL+OXNtuFSgS+csmq0yRkqlUd3oPnU6HZrP5yIioRPgul4toNKpmz+rqm6OHFnrNQ0NSFYFAgHK5zN/5O3+H3/7t32Z8fJxarYbX61UbmjLRqVKpMDU1xaefftoz6u4wsHbC7ibGzWaTmzdvKnuDWq2GYRg0m03m5+dZWVnZ9fWltFHSPYeNpJdCoRDPPfcc77zzDktLS+qORHN00EKveWiIQJfLZXw+H9/5znf4zne+o35vHd0ng7ivXr3K1NQUs7Oz6pjDYKeu150QN816vc6NGzdU9VC322Vtbe2O0bDP58Pv9z8yQm+9sxoZGSEej5NKpZTQ6yaqo4MWes1DwyoKg4ODnD59mr6+PkzTxOPxKIGX9E2r1WJubo6lpSVlVvaopDV2Q7pLu90u1WqVtbU11SRVqVR2LQV1u90EAgHC4TB+v/8hX/Wd6Xa7BIPBnkVIWyIcLbTQaw6Fb37zm3z5y1/uEX+JIEVE8vk87733Hr/+9a975p8eBns9r9XgTPYXpOFKHhPht6ahUqkUx44dY2JigkQicTBv4j4xTROv10swGCQUCqnrl/eko/pHHy30moeGCEQ2m+WNN94gFArtKhLLy8t8+umnfPjhh0daSKymbIJ8LyKfSCQ4c+YMQ0NDvPjiiwSDwUfiPcs1WG2Wtzt16vTN0UALveahYZom8Xhc2RxYo8Ltx01PT6tyPmu6RqL9oyIucu3SROV0OtV79nq9uN1unnrqKZ5//nm+9KUvcebMmcO83B0RGwhZqK0LADz66TSNFnrNQ6Tb7TIwMMDJkycJh8Oq3X47zWaTXC7H6urqjuWIRzGK9Pv9RCIR3G63akqSzdexsTFOnDjBxMTEofrR74a1C1h+Fh6l69TsjhZ6zUPDbrfzrW99i29/+9s4HI5dI/rNzU0uX77MtWvXWFtb6/mdRJdHBa/Xi9/vZ2JigsHBQTweD6FQCL/fTzAYxOv18uUvf5kvfelLBAIBVVb6KCxk8jlLNVSj0bhN5Hf6+9M8emih1xwYcqsvQ0VisRivvvoqw8PD6vc7CdrS0hIrKytsbm7e9rv7SRNY7wCklt86wm8vrynXutM0KeviY/WT9/v9nDp1ilgsxrPPPks2myUejxMOh4lEIiSTSYLBIOl0WnXc7vaZHAYi4g6Hg3w+31NDL3+nsoGuebTRQq85MESwRASHhoYYGBggFArd8Xli/SveMveDVdytwmndFL3TCMDt3OlY6++sdxvieBmLxchmswwODpJKpXC73aRSKRKJBC6XC6fT+ciI+3akXHR9fZ1yuaxSaTu9X82jixZ6zYEh9eRiffDKK68wMTFBIBBQEfVOrK+vq1mlu7FTnt5qWbyTH45EyzI0BL5IS9xtQtROG8JOp1OZksnvTdMkGo0SDAZ55plnOH36NNFolLGxMYaHh8lkMgQCgUemKepOSF5eho+srKyQz+cB1D6C3og9Gmih1xwYViEfGRnhN3/zNwkEAsqsbDfW19dpNBp3dEoUkReHyGAwqBwv2+02lUqFWq2mhFwshsPhMG63W431a7fb6k+pVKJcLqvu3Du9L7E1EJdKt9uNYRj09fUpMR8bG2NwcFA5Wspx1o5TWXgeRdrtNk6nk7m5OT777DMmJyeVkZt1YdM8+tz1X5hhGP8O+HvAqmmaT916LAb8Z2AEuA58zzTNgrH1P/tfAd8BqsD3TdP86GAuXXNUMAyDSCSiGoF2SqlYabVatNvtO0aL4gsjnaQy2xS+ME8Tr/h2u61SJA6HQ32V+nCfz0ez2VR/2u32nux4E4kEHo+HQCBAIBDAbrfj8/mIRCJ4vV58Pl9Po1Sr1VLXJlOoHmXk76dUKlEsFikUCofeuKa5P/YSSvx74F8Df2p57IfA26Zp/olhGD+89fM/Bd4Axm/9OQv8m1tfNU8grVYLh8NBf38/X/3qV3n11VdVRyh8sXEpKZd6vc7GxoZycBQhljSBiIvH4yESiRCLxUgmk7hcLhwOB7VaTS0SItrWwSSmaVKv17Hb7arE0el0KssFn8+nNmslneN0OtWxXq8Xh8OB2+1WQ0XkexlIYrPZVGrG6XRSrVZpt9u43W5VNhqNRgkEAoRCIdxuN36/X6WO5LM5zLmxgLrrajQanDt3jvPnzyufHi3yR4+7Cr1pmj83DGNk28NvAl+79f1/AN5hS+jfBP7U3PqX8L5hGBHDMPpM09R2d08gMjz7u9/9Lr/927/dI2KSthBPG6fTydraGvPz8ywvL9NsNnuEXsQlHA6TTqeVba6IZLPZVG6RuVyO9fX1Hu/7RqOB3W5nc3MTr9eLy+XC6/X2VNNYUzYi4i6XC8MwCAaDBAIB1fgkfjROp1OlY0Toxa9+fX1dVaZUq1Vu3rypricUCpFKpYhGo4yPj5NOp9XndViRvnVyltPpZHNzk88++4x3332XDz74QEfyR5j7TQ6mLeK9DKRvfd8PzFuOu3nrsduE3jCMHwA/uM/za44A3W6XoaEhvvOd7/DMM8+ox62CIlF/vV7n2rVrzMzMMD09zfr6OtVqVVWkSMQdj8dVtYrH48EwDDWk+8aNGxQKBYrF4o7XIyK+vWxzp1pwl8ul7hokqpfrbrfb6rqlvlxeQ6ZPbd/gtQqkDBpPJBIkEgmq1Sq/8Ru/oawPDlNQG42Gsol+++23mZyc5OLFi1SrVS30R5gH3gUyTdM0DOOe//ZN0/wR8COA+3m+5mhw+vRpxsfHCQQC6jFrC32r1cIwDMrlMouLiywsLJDP5ykWi1QqFZW2sdvtSnAlwpcUTaVSoVAoUCgUdqzUsQqU3CFY2WkvoNlsUiqVsNlseDweNefVZrOpCF3y/NY7lUqlQrPZVEK/vZtUNoBlL6DdbrO2tka1WiUYDO7DJ/5gyF3U6uoqly9f5saNGywsLGiBP+Lcr9CvSErGMIw+YPXW4wvAoOW4gVuPaZ4Atg8FCQQC/O7v/i4jIyPAzq3zIsKff/45Fy9e5Pr168zOzrKwsKDKFmVBqNVq+P1+6vU6NpuNQqHAxsYGpVKJ1dXVXQeSWM+713JAsRmWKhNJ81gXnJ3Er16v9wwo2emaYrEY9XpdbdDmcjmKxSLpdPpQ69JbrZaqXHrnnXc4d+4cV65cYXl5GdgqTz2soS+aB+N+hf4nwO8Df3Lr619bHv9jwzD+jK1N2KLOzz85bBeBl19+ma9//euq8mWnkkqXy8XS0hIXLlzg/fffZ2Fhgfn5ebXxZ+00FX/3SqWi0gu1Wk1FxwfJds+dB0lj5HI5Wq0WkUgEv9+v8v7yuoeF0+kkn89z8eJF3nrrLf77f//vVKtVdRekRf7ospfyyv/E1sZrwjCMm8A/Z0vg/9wwjD8E5oDv3Tr8v7JVWjnNVnnlHxzANWseYaypkVOnTtHf33/XGaMrKyssLi6Sz+fvOBy72+1Sq9XUhqXkyh9GWsG6QbrdvdHK9g7ana5Nyj6Bnjr8OzWRPSwWFhaYnZ1laWmJjY0NgEe+DFRzd/ZSdfM7u/zqGzscawL/+EEvSnN0MU0Tl8tFKpXizTffVJGvCNl26vU6n3zyCefPn2d6elrNWLXWn1sHXOzWzLRT7n0/2emcO4ny9vTUbsd4PB6SySSZTIaBgQGCwaBKUx0WN2/e5Cc/+QmTk5P86le/Ar7wtIGj6Rqq2eLRbMnTHFlM02RiYoIXX3yRF154Adgak7cbk5OTfPDBB0xPT/dUdljLL3drXrJ60z/MVvy9mKHJYrXTtQeDQfr7+3nmmWcYHx/nK1/5itrwPcyI/he/+AXnzp1jampKVS5Zexes5aqao4X2GNXsGxKNHj9g5nGNAAAgAElEQVR+nNdee02VC1pLKQ3DoNFoYBgGq6urzM3Nsb6+zsrKCvBFmkAE5m42CIcRYe5lYZG7D6nIsdoc+Hw+kskkyWSS0dFRQqGQqix6WO9JziOLVj6fZ3Z2lunpadbW1tSCI3l5LfJHGx3Ra/YFEQav18vrr7/Ot771rduOkYoO2ZS9cuUK09PTXLp06TYPlccFEUoR78HBQTUY/fnnn2diYuJQrstms9FsNnG5XLRaLX7605/yzjvvMDMzQ7vd1imaxwwt9Jp9QSLE0dFRXn31Vfr6+nYVC7vdTqlUYmpqimvXrjE/P6+OfdyEXpCh4P39/fT39/P0009z5swZYrHYoXTD1ut1PB4Pq6urXLlyhZ///Od89tlnqiNZ83ihhV6zr7zyyiucOHFiR5Fvt9sqhTE5OcmlS5eYmpqiUqk8tht9sknscrkYGRnh2LFjTExM8MILLxCJRIDD8Y5xuVw0m01+8YtfcOHCBT7++GOVPntcF9snGb10a/YFSd0kk0m8Xu8dxavb7SrP+WKxqPLw0nX6OOJwOAiFQoRCIWXjIJHzYW3AFgoFFhcXVUeyeNyIhbPm8UFH9Jp9wTRN4vE4L7/8MkDPRqQg7f8LCwtMT0+ztLTE5ubmXatrdsP6HJvNprzmJTW03dPGajpmjVpjsZjytJcB3vF4HI/HQ61Wo1arsba2xsrKCuVymc3NzT1fq5RMjo6OcurUKV566SXOnDlDJpNRZm4Hnbqxzn6VSqBSqcRf//Vf89ZbbzE1NUUulwO4a8+D5miihV6zL7jdbk6dOsWLL74IsON4PLvdzvz8PFevXmVmZobV1VVKpRLwRbrgXurhrSIv4/oCgQBut5tCoaA6bOUcUgljt9vxeDzEYjFSqRTJZFJVwmSzWbxeL5FIhEgkojxuisUiN2/eZG1tjWvXrrG4uKjuSHYTR6fTSSwWIx6P88ILL3D27Flee+01jh07dttxB4kYlUl67MqVK5w/f56f/vSnnDt3jkKhcKDn1xw+Wug1+0IwGGR4eLhnAMhOVCoVNSaw2Wzui8e5eMJ7vV418EMsBhqNBvV6nWazqSLbQCCAz+cjm82SzWbVc+PxOMlkEr/fTzqdJpVK4XK56HQ65HI5QqEQ6+vruFwu/H4/KysrqoO00WjcNm4wEAgQjUZJJBKkUilGRkbo7+8H6PGeP+hGKal2MgyDWq2mpkXNzc2xublJu91+bPdINFtoodfsC2fPnuXNN99U1r07pSK63S7Ly8sqGpaSyr00IO2E2Af4/X411cnlcuFyuZQTZKvVIhAI0G63sdlsSnjlOaFQSI0hjMViDA4O4vP56O/vJ5PJ4PV66XQ6tNtt1tfX2dzcZGJiguvXr7OxscHMzAzLy8usrKwoEzR5/6lUiuPHj5PNZvnqV7/K17/+dVwuV88glZ1SXAeBYRgUCgX+9m//lr/4i79gcnKSyclJPRLwCUELvea+8Xq9qpHmt37rt/jqV796x+5O0zTJ5/NqYLjH4+nZ+LtXGwOn04nf78flclGtVpXRGXxRzihuk06nU6Vn+vr6VFTvcDjw+/3Kjtjv9+Pz+QiHwwQCAWU37PP58Pl8dDod1dm6vr5OJpNR/jClUolqtar2HFKpFGNjY4yNjXH27FncbrfysQdUk9R+I01p1o7ka9eucfnyZX784x9z7tw55ubm9v28mkcXLfSaB8br9ZJIJNSMVImet2MdGyiRt9PpvC+fGrvdrjZWxe5X/OJrtRqmadJoNNS4wGAwiN/vV/NdA4EAXq+3Z36sDBeRaxXBluuTY0KhkPq5XC6rwSder1dZBzidTpLJJP39/aTTafx+/4EJ+05Y76yazSYLCwtcv36dubk5ZTuseXLQQq+5b2R8n2zCSrpkt3SEzWYjHo9TLpdJJpMquq9Wq6r6Zq/IABKprJG7BcMwcDgcKn3jcrmIRqOq4iYYDCrRF4/5er2OaZpUq1Xy+bzKaXe7XTU+0PqYYRhq8Ugmk2rot7wPOX82m+WFF14glUrh8/lUtc9Bi70sULD1d/HBBx/w05/+lE8//ZRLly7RaDQO9PyaRw8t9Jr7Rqpe3njjDdLptHp8t41Fh8PB2NgYLpeLXC7HxsYGDoeDQqFAuVy+p6i+2Wze9pjkvUXsRbSr1apKvdTrdeWQKXcYspFar9dZW1vD4XCwtrbG4uIifr8ft9uNx+NRQ8hlkRHLZJvNpo4zDAOPx4PH4yGVShGLxQiFQur6HA7HgefDrYPNP//8c37+85/zi1/8gtnZWWU9rHmy0EKveSC8Xi8TExN3dKgUTNMkGAwSjUZJJpNEIhEqlcq+lhda7Y3FAbNYLLK+vo7dblc5+Ha7jcvlwmazqalQMrrQZrNRLBZZXV1VAi/DxN1uN4FAQE3Tqtfrak6sRNKSLvJ4PGo0n7z/h9GIJAPXS6US09PTXLlyhaWlJQqFgt50fULRQq+5ZxwOh4q+pRHIKiC7Ra2yqTk8PEw+n+f69etUKpV9E5+dnB87nQ6FQoGpqSny+bzKp8tkJ7vdTrFYVKItgi0TlazCHYlE8Pl8KvUj0X2r1VIpGVnM3G43brdbjUO0XuN+Y/Xsl6/r6+tMT0/z/vvvc/78eVZWVpSPjbY4ePLQQq+5Z6SyJRwOk81myWQye3qe1I7b7XZVu+7xeHrGBe43ks6x1u3n83kl9FKh0ul0qFQqlMtlOp0O1WpVRfqS35f0jNTF+/1+9ZjP51N5+2azqcYbymaxtaRyv6N6WXitG9PvvvsuV65c4f3332d2dlZVR0kpq47snyy00GvuGbEeiEajDA0N3XPqRSLvbreratThYCYYmaapolgZxA1QLpdVlU273abdbivbhJ1seg3DoFwuA5DP51lfX1fds9FolHg83lPdEw6HlbWCnEPKP/ezbt5aJSTdyOfOnePnP/858/PzTE5OUqlUeo7VIv/koYVec8+IX4rUne/VH0UETqLOTqdDvV4/8KHT1qEeIuqSboGtjd1Go6FKJXdCZtTC1qBwqa6RyVA2m03V6kejUaLRKKFQCI/H07N/cRAia/WwWVtbY3Z2lvn5eRYWFlhfX1clrXKs5slDC73mnhExDAQChMPhHStgdsI0TRU137x5k+XlZYrFonr+QUWa1rSQlFXKRimg0ixSh29Nb8hzO51Oz5hDSYVYO2r9fj/PPvssExMTBINBwuGwOq8IrETe+43dbmdqaoqPPvqIt956i1/96lcUCoUDS4lpjhZ3FXrDMP4d8PeAVdM0n7r12P8G/BGQu3XYPzNN87/e+t3/Cvwh0AH+Z9M0/+YArltziIjgZTIZgsGgym/LJmylUqHVatFut3s2JRuNBjdv3iSXy3Hp0iVWV1dZWlpSkfK9bhQahqGiZem23e04aZKKx+NkMhnsdjsul0sZepmmqapSrMPJraMN5drEW76vr4+nnnqKM2fO8OKLL+L1ehkZGVE191ZXSvnMPB7PvX7ct72X7dhsNqanp/mLv/gL5ubmuHDhAqurqyrS11G8Zi8R/b8H/jXwp9se/79M0/w/rA8YhnEa+EfAGSAL/H+GYUyYpqnDiscIydEHg0G8Xi+NRoNGo0Gz2aRer7O+vk6z2aTT6Sih73Q6FItFFhcXyefzrK6usrGxoapcABVN7zUKleYn6ZCtVCp0Oh0ajQbtdltV1YhhmTQ4pdNplc6Rxanb7eJ2u2k0GiqSt+4ZWCuNQqEQIyMjZDIZBgcHeeqppzhx4oTq+D0odtrINU2TcrnMxYsXuXTpEisrKywvL2sPG00PdxV60zR/bhjGyB5f703gz0zTbACzhmFMA18BfnXfV6h55Gi323i9XoaGhohEInQ6HcrlskrNLC8v02g0qNVqAGojNJ/Ps7CwQKFQ4Pr16ywvL/dsFO5V4D0eD4FAgNHRUbLZLE6nUzUrSYeq3CW4XC5CoRCJREJ9FfvhRqPB8vIyN27cUOWH169fVxUy1hROu91WdwSnT5/mlVdeIZVK8Y1vfINsNqsEXu4qpHt2P4VW8vByx9Htdvnkk0+YmppS6ZpcLtdjFqeFXgMPlqP/Y8Mwfg/4EPgnpmkWgH7gfcsxN289pnnM6OvrY2RkRJUaSvWMCG29XmdjY4PNzU0KhQKVSkW5PFYqFTY3N1WHajQaxel0UiqVdhV7q2hJiWN/f79ym0wmk6oBStIuYkUQiUSIxWLKNjgQCNBsNpXHfDQaZXl5WXneFItFZc8AWwKbzWYZGRlhaGiI559/nrNnzzI2NkZfX58qbZSRgfBFqmc/Symt9gn5fJ4LFy6owSHvv/++GgUoyB2AFnvN/Qr9vwH+d8C89fX/BP7He3kBwzB+APzgPs+vOWR8Pl9PDXm1WlWCJwZjjUaDarXKxsYG5XKZXC5HPp9XjUkijNJlKp2pOwmTtTRQcu5Sxx4IBIhEIsqXRvL2kkqJx+PE43HlWCmpJKfTSa1WI5VK0e12yefzrKysqOqgarWqFotMJkMqlSKTyZDNZhkcHCSRSKhNXK/XC2xF9A9jFN/S0hLnz5/nypUrat8DtszUrJVDchegebK5L6E3TVOFDoZh/N/AT2/9uAAMWg4duPXYTq/xI+BHt15DhxyPKFKWuH1Db3x8XA27FtGVKF7y9rVajZWVFebn51lfX2d5eZn19XVgS4RjsRgOh4N4PK5KEUulEmtraz35eqtXvXSiWksh7XY7wWCQRCKB1+slEAio9I7D4VA/C+LqGI/HaTabyoZBXm9zc1N5z4uVwfDwMOPj46RSKZ555hmGh4dv+4ysdxPyHu8Xa6MWfNEPsLm5ybVr1/jJT37Cz372My5evEipVFIL4PbyUC3yGrhPoTcMo880zaVbP34X+OzW9z8B/qNhGP+Src3YceCDB75KzaGx04xVn8/H888/z7Fjx0in0yqdEAgEGBwcpFQqUavVsNlsVKtV1tfXyefz5PN5Njc3VTolHo/j9/tVRB0MBimVSpRKJdWuv/0aOp0OtVqNpaUlZRhmmiYbGxu43W5M01SulB6PR224whd3A9aySsMwCAaDagNWPPYl/SR+9vF4nJGREdUkddApEeuCIe+92+3y1ltvcfHiRd5++20++OADtSegq2s0d2Iv5ZX/CfgakDAM4ybwz4GvGYbxHFupm+vA/wRgmuYlwzD+HJgE2sA/1hU3R5vtzUyGYTAxMcFv/MZvEIlEsNlsNJtNVRseCoXIZrMAyjZY0jEul0s5QcZiMTKZDPF4nHQ6rewGYMsoTSyQ4faotNvtUiwWmZycpNlskk6nqdVqlMtlgsEg5XKZaDRKrVbD4XAQCoXwer1KDGU/oVwuU6lU1KZxIBBQE6DkPLIhK4uT2+1W6ZqDnAxlvRuQTuJr167x4YcfMjk5ydTUVE85qfUOQEfxmu3spermd3Z4+N/e4fh/AfyLB7kozaODNWqV1MTQ0BADAwNq89NqaeBwOAiHwzQaDSKRCMFgUI3ykwja5/ORSCTo6+tTo/02NzfJ5/MqL74XNjc3WVxcpNVqqU3QYDBIvV5XjVkyKMTr9fZU9oifjZSBynuVBUF8cORxSZ3IAlEul4lEIvv5Ue+KnH9tbY3l5WWWl5d3tRvWG6+andCdsZo7Irlrj8ejBnQ8/fTTnDhxgna7rcRdLHpN0yQWixGLxSiVSkqIC4WCGtnn8/mYmJhQqRCfz0epVCKfz1MqlfbkySKRq9TlLy0tMT8/j8fjURbIcschs2TldUXEZTNYcvhSkeNwOFQlkET/9Xqdzc1NldOv1WqEQiHGx8fV9VjvbB4U2UeQr+12m88//5zr168zMzOzo23ETu6dGg1oodfcBbEFkKh3dHSUb3zjGwDKFGy3MsLTp08DEIlEWF1dpVwuq4HeIyMjZLNZvF4v3W6XhYUFZRsgYnmnVIT1TqLRaJDL5VhdXVWe816vV5VL+v1+NZ9W/sgCEIlEyGazqgkqFAqpWbRSnSMTpBqNBpubm1SrVex2O4lEQhma7bcrpex7SEnl4uIis7OzLC4uUiqV1OejhV2zF7TQa+6KVLoAJBIJRkZGen6/m8h5vV76+/tVGkdSKV6vl2w2SyKRUHcKhUJBDfawVprstsG4PaK15vNLpZJqxJKh3xLRi+DLWEGZENVsNolEIj1DteU4uVsRYZUU1MbGBpVKRQn9bpO1HgSrzbKkjHQOXnOvaKHX3BUR+lQqxbPPPksymVQCbHWG3Ol5AwMDZDIZWq0WtVoNu92uom4pB5ScupRpBgIB9Rr3ErFaFwgRw263q7p2t+NyuQiHwwwNDZFIJABU92w4HO65CxBPeWudv/jrWFMs+xVhi0+NDDuXcs9isairazT3jBZ6zV0R0RwbG+Ppp59WkS7ceTyeNE9JHl+GdMvvpGlJpjcFAgG8Xi/RaBSPx6MshXdC5raKsFpr7bez/Q5Bqm+azSa5XE41a7ndboaGhtTAcNjK59tsNmKxmBqqLeeJx+MPbFK2G1af+eXlZZaWllhfX++xjNAVNpq9sv/3mprHCmvaJpVKMTo6qurVRWR3i2KtfjOyAWpNx0geWmyCbTYboVCIdDpNNBpV1TcS2UqpI3wxAFvY7RpkQZGIHG73ZJeu3c3NTcrlMo1GQ404lHmxVoM0h8NBLBbD5XKRTqd7KnL2C+tCKB3GjUaj5w5Ki7xmr+iIXnNHOp0OLpeLZrPJs88+y9DQ0J6fKyWPUrsuUej2qppKpUK9XlcVMAMDAxSLRTY2NlhbW2Nzc1O9pjxXxE5cM/c6/GQn2u02jUaDUqnE+vo6Pp9Pva6cKxwOU61W1QZtKBTC7/eruntJSe0X1o3Wer1OrVa7bXHTaPaKFnrNnnnhhRcYGxvb8/HWPDOghFksDKQDVXLPMkt2YGAAh8NBtVplbW1Nibj1bqDdbrO5uUmpVKJcLrO4uHjfEbXcndjtdrUQybms0bqkqTweD8FgkFQq1RPJ72fljdy9yJCTWq3W44yp0dwLWug1d6XZbDI4OMjExAQej2fPUaXNZqPRaLC+vk6tVlMWxZVKRW0qiqCJKRigOlDFHlg2QcUoTXL36+vrFAoF5Wu/tLR0h6vZnWAwqAaSSJQuM2TFoExKK2XxcrvdPakl2Y/Ybz96a7pqr5O8NJrtaKHX7IlwOIzf77+n54hP/crKCpubm8zNzbG+vq7KBCXSd7lcqkFJBoBsNwmTyFacLxuNBvV6XQ09eRCBFRvjSCSiSjzF5kCslKvVqtoIFVEXDnrottWfR6O5H7TQa+6Kw+Hg+PHjDAwM9IzHsyJi12q1lGBfv36dGzduqNmwN27cYHFxUQl9s9nE5/MRi8VIJBIkk0k1dtDqLW/Ny1srYsShstPpEI1Ge4ZuSMmjXNtOk6KCwSCZTIbx8XHi8biq7Q8EAmxubrKwsEAul6PZbFKr1Wg0GthsNuV7Ixu8UgLp9Xr3tbxSqpW2u2Ie9DB1zeOHFnrNXel0OqpJyjr8YiecTifNZpNyuUyhUFCbqhsbG6pMsFgsUiwWlTtkPB5nbW2NQqGgrAisNeySshC/HNnYFd/7breL3+8nFotht9vZ3NzsGQUouXPJeYfDYTWhanR0VA32ttoZt1otFhcXld1BqVRSdxsyjjAajaoB4AcRzUsZqFgvb/cA0p2xmr2ihV6zJ6LR6J6P7Xa7Ko++vQRSKnBko1NKGZ1OpxryUa1WlcWwjOOTckJ5PVlQms2minx9Ph+tVmvHASY2m00ZnA0MDBCNRhkbGyObzRKLxdRsW7n+crnM+vo66+vrangKbFk0m7cGiWcymZ7X30/RldSVaZr4/X58Pl/PPgZoodfsHS30mrtiGAbpdJpOp7NrrlhE3Fq54nQ61fg+eY1qtYppmioV0mq1KJVKFItFFhYWVNesNXUjqQqxM3C73aq+XYjFYjidTur1OvF4XKVw4IvBJG63m1gsRn9/P5FIhP7+fuLxuFrExDOnXq8zNzfH3Nwcq6urVCqVHu/6TqfD0NAQyWSS8fFxTNPs6RPYD6xjABOJBJlMRg1WEXM57UGv2Sta6DV3RTzoJWWzk5hZBUdG/QWDQWAr+g6Hw6qKxu/393SjSqPSnWrhxQzN7/eTyWRwOBwkk0mCwSAej0f5zXQ6HRqNhqrokWv1+XzK8iCZTOL3+9Xma6fTUX47N27coFKpcOXKFS5duqRSTLKf4PV68Xg8apiKfB53S2ndKw6HQ+13OJ1OJfShUIhGo6Ejec09oYVec1e2bwbuFYfDoUbxwVb6JxaL0e12iUQianj4nczLBLEKlgoY8cPxeDz4fD7lOQ8oa2H4okZejM1kvqzX61WTqNrtNp1Oh0qloqySFxcXKRQK6q4DUJvA0qUqHHRdu6SwHA6HumOSzVqNZi9oodfcFRHI+3meDOMWT3vDMIhGo2okoGEYrK2t9aSFdhIw2XyVMkvZoBwcHCQWizE8PNyTa5fOVtkTAHoGihiGoYaIrK6uUiwWuXz5Mp9++imlUonp6emeBiVZjOr1uuoWFh/+gxZ6p9OpFkmfz4fT6VT9BRrNXtBCr7krbrebYDC4q+/8bkjUKVGo3W4nnU7j8XjY2NhQ5ZDWNMRu4iXHWMcQjoyMcPLkSZLJJMPDw8TjcVwuF3a7nVqtRq1WU+ZlYiHQ7XYplUq0Wi21L7C4uMjKygoff/wxV69e7YnWreeGLdENhUKkUilSqZT6PPa7zn375+Dz+ZSjptPpVOMPNZq9oIVec1ccDodqYoK9pyqsRmTSJStRtdvtVkZh8np3cp8UM6++vj7GxsbIZDKMjY1x/PhxYrEYg4ODalygVNxYbYVF6Ov1uhqCMj8/z+XLl7l27Rpra2ssLCz07BNI1Yu8puT1h4aG1GYsoOwTDhLZ3Ja6evk8dE29Zi9oodfclWg0SjAYpNlsqtr27Vj9bOQY61APQOXCraWP1oao3apWZNZrLBYjlUqpOvbBwUHC4TDhcFgNFpENWYnipaa/Uqn0TIgql8ssLy+zsrLC3NwcxWKxp/ZexN1aNun1eslkMmQyGdLpNAMDA3d18LxfRNhhawH0eDzqq3zWWuQ1e0ULveaOSE7d7/ffs0NjMBhUPjeyQVosFtV82PX1dXK53F3TEJK/T6VSTExMcPr0aTKZDIODg2pjVWa5NhoNyuWyOo9U91SrVer1Omtra+RyOTY2Nrh06RJXrlwhl8up92ptsrKmY2w2G9lslpGREcbGxkin0/f3gWo0h8Bdhd4wjEHgT4E0YAI/Mk3zXxmGEQP+MzACXAe+Z5pmwdgKz/4V8B2gCnzfNM2PDubyNQeN3W4nmUyqTcx7aQzy+/0UCgU1RWplZYX5+Xny+TxXrlzh+vXr3Lx5c09mXdFolKGhIY4fP87IyIjyphGbYOms3dzcVN7ykqOXCL9er1MsFtUQD7E4sL5Xa5RsTSVFIhF1bnHX1GiOCnv519oG/olpmh8ZhhEEzhuG8TPg+8Dbpmn+iWEYPwR+CPxT4A1g/Nafs8C/ufVVcwQRZ0erd8xekaoUq9d7LpdTNejFYrEnmt9tAZEmKknT+P1+gsGgapjaXj2zvLzc0x0rTVci9mKPLJuu1ruUndJHTqeTZDJJKpUiEomo6VcazVHhrkJvmuYSsHTr+03DMD4H+oE3ga/dOuw/AO+wJfRvAn9qbv1ved8wjIhhGH23XkdzxDh+/DgnT54E7u5zsx2ZTrW6usrCwgKXLl1iZmaGfD7P3NwchUJBHXunLk+Hw0EgEFBWwmJ+1ul0VOQ+PT3N0tIShUKB5eVl6vU6brcbwzDw+Xxqr0DSM3a7XeX1rQuY1KrLIhAMBkkmkzzzzDOcPHmSWCxGOBzWnvCaI8U93X8ahjECfAk4B6Qt4r3MVmoHthaBecvTbt56rEfoDcP4AfCDe75izUNBhPfEiRN86UtfAlD123sV+3a7TT6fZ2pqirm5OX79618zOzurxvRZ2d5ZC19E+DK2T2wMpFRTcu4yeGR2dpZ8Ps/q6irtdptAIKDm0Xq9XtX45fV6iUQiypJhbW1NnUvq06WD9vjx42SzWU6ePKn2BNrtNsvLy6yuruJ0OvH7/aprVxYX6/D0+0G87WXoeK1WIxKJEI1GCYVCPYukRnM39iz0hmEEgL8E/hfTNEvWf8CmaZqGYdxT2YFpmj8CfnTrtXU/9yNKJpNRXjDbjcLuRqPRoFAosLCwwNzcHDdu3KBUKt31ebIPIGkX2QS2/ptrtVqYpqkqaEqlEhsbG8o3R2r+pTNXDMvELwe28u59fX0qhy8llE6nk76+PgYGBjh27BjJZJKBgQG1YMjoQTFVE7EPBAKMjIwQCARUuuh+yy6tJad2ux2Hw0EoFMLn86k7FD0zVrNX9iT0hmE42RL5/8c0zb+69fCKpGQMw+gDVm89vgAMWp4+cOsxzRFCBP3UqVNkMhnlu3Ivm7GlUomrV6/y4YcfMj8/r6JQEe/t7pbC9lSK3+9Xk61EXEXoKpVKz0hB+SONUpKekV4Al8uF1+tVpYvRaJRjx47R6XSUBbHb7SaTyZDNZolGo2owigwcKRQKasM3l8vRarWIxWJEo1Gq1SovvvjiAxucyQIhdwZOp1PdXWxsbDA1NaXdKzV7Zi9VNwbwb4HPTdP8l5Zf/QT4feBPbn39a8vjf2wYxp+xtQlb1Pn5o4f4w/T19SlDMnl8r0xNTXH+/HmuXr2qInkpWWy1WrumNax5dImiZbyfdb6s1dnSKowygapSqWC329Xi4PF4iEQiBINBnE4n6XSaQCCgouS+vj6i0ajyyJcNX1lQpNHq+vXrXLt2TaVvOp0O6XSaeDxOs9kkm82SSqXUYnI/yIIq5zcMg7GxMcbGxmg0Grz77rvKe0ejuRt7iehfAf4H4KJhGJ/ceuyfsSXwf24Yxh8Cc8D3bv3uv7JVWjnNVnnlH+zrFWsOHGvTkLTcC5Iz3gsbGxsUCgXK5fJtXbUSZe/kWGlN1+A41B8AACAASURBVLjdbvW9NENJ2sZ6d2HtHJVFoFarYbPZVJVNp9NRuXTJ1YfDYdLpNOFwmOHhYYLBIHa7nWg0qoaeyKLRbDaVN/3q6iqrq6vk83k1kESGnkhZ535E3NbXENtnSS9pF0vNXtlL1c0vgd12lL6xw/Em8I8f8Lo0h4iYgiUSCYaHhwGU2EvJ5HYkVy12AEtLS3z66afMz8/3dJxKrl0MygA1nk8Mu8SP3jRNvF6visClPLJSqRAIBJRoG4ZBPB5nY2ODZrOpJlvJJuvq6iq5XE41VrVaLRKJBPF4nJGREU6dOqUsgK2dup1OB5/PR7fbZWNjg/n5eZaWlpicnGRycpKNjQ3y+bxaOGq1GseOHcNmsz1wnb31DkoW1k6nw/j4OPPz86TTaTY3Nx/oHJonB931oelBXB0B0uk0fX19eL3euxqaSTNVu92mWq2qiFfESIRaonHB5XKRzWbp7+8nlUqpBiixB5Z8vcPhoNFoqPp3yb2Lz3y1WlWdsAsLC8pD3up1I2kcsTuemJigv7+fY8eOqVmyYntgTQGtrq4yNTXFxx9/zM2bN7lw4QIzMzNqwZJ0lETc4rmz30jzWiaTIRQK7fvrax5ftNBrerCKcCQSwefzqcfvZMkrEawsCJK+kKYmSfmIcIfDYfr7+0mn0wwNDZHJZFQ1SbVaZX19nUajoTxopJxSZrjKIG4xRwsGg4TDYcrlMg6HQ02YsqY+Op0Oi4uLynbZNE1CoZC6dqvlgfQATE1NcfXqVT766CM++ugjlpaWWFpaUouh0+kkHo8zOjpKX18fzz//PH19farEcj8FX6yeZeHVaPaKFnrNrojwiTjfSeytIikOi1KvHgwGlfh7vV7cbjeDg4OMjIzQ19fH4OAgiUQCh8OBYRgUCgV1Z1CtVlUev9FoqBGEEvHL8BB5XRH+7Y1Q8IXoS97d+r7k91aKxSKzs7OqNHRhYYG1tTW1iNhsNvx+P6lUiuHhYQYHBzl27NiBiLxc3/0OgdE82Wih19yGpDskat6rsIhzZaVSUQ6TYv4lVS/ihDk+Ps7o6CiJREJVv1QqFZV6kY3UfD6v0jDNZpO1tTVisRg2m410Oq02X30+H8FgkFAopHLm25GSTLFBaLfbPdU/7XZbpXDq9ToXLlzg3LlzzM7Ocv78eRXJdzodtWdx+vRpxsfHefXVVzl+/DgTExMHNoxEegBcLpfqBdBo9oIWes2+Ifltr9dLLBbjxIkT1Go1SqUSnU6HUCikatNHRkbUHFS320273SaXy9FoNNjY2FBpHqlwaTQaqtJEZs6m02lCoRBut1t1unY6HZLJJOFwmFKp1JOKkghehFyGn8idhoh8o9Fgbm6OK1euMDMzw/z8PMViUdkiuN1uhoeHGRkZYWJigjfffJOzZ8+qUYYHicPhUAuaRrNXtNBr9g0RVWlyEhsB6UqNxWIMDQ0RjUbJZrMqrQMo7xkZ7l2pVFTH69ramvLCh62yTa/Xq/YBbDYbbrdbmZ1JVC+pJuhNy4hlspRdVqtVVZrZ7XaVSdrGxoayPhaktl/2F44fP87x48fx+/0H/vlKNY9s+Go0e0ULvWZfkC5XEVSHw8GpU6fweDxsbm6qJqS+vj7VpGQ1MjMMg3q9rjpcb9y4wc2bN7l8+TL5fB5A1anLRur6+jqdTkd1vUqZ5smTJ5mZmWFlZYVqtapm00ruXGyT8/k8y8vLyio5EokoX/uVlRUqlYqqzw+FQkpgR0ZGeOWVVzhx4gSvvfYamUzmtiElB4XYLIyOjpLJZFheXj7Q82keD7TQa/YFq69Lo9FQc2b7+vpUeiYYDBKJRHqmUVlLMqXb9Ny5c5w/f561tTUl8rC1GFQqFebm5vB6vRSLRWArlRKNRonH44TDYer1Ojdv3lQNTdPT0yq6F+sFMUK7du0ayWSSVqulyjGlozYajTIwMKA6ZyORCPF4nKeeeorXX3+dkydPqqqih7FBKnsI/f39PP3000xOTmqh1+wJLfSaHmQEIHBb5Yh12Pd2rJuD0uwEW+WHMhvW6XSq5impinE4HNTrda5fv878/DwbGxssLy+ztLSkhFwif3lN8beRDVXTNJUPj4jz6OgoJ0+eZG5ujoWFBVXPL++n1WpRq9WUe6XYI7jdbjweDz6fj1gsRiKRUJU1gUCAaDTKl7/8ZYaHh1WD2MNCUks2m43x8XHV9JXP59XvrOkqjUbQQq/pYfuEJWs6RgT6TmwXGmuk22w2ldBLBF0ul1XUffHiRTY2Nrh69SqFQkEtNNstf625dbkrkPMYhsHAwAAvvvgia2treL1elpeXmZ2dVSkccassFAosLi4yMDCA1+slGo2SSCTUa7jdbrLZrLIM7uvrIxaL0d/frywZHiamadJutzEMg7Nnz/LJJ58wOzvbk956GOkjzdFDC71mV+r1uhJbsQbYLUUhFS1Wi4N8Pq8iaREpsfctlUqqzFEGkUxNTbG5ucnKysodG4KsUbnNZlOljjJ43Ol0Mjg4yNNPP43f72dhYYFiscj6+jrNZlOVgYro1+t1VZMv5aAej4fjx4+r62i32z31+bVa7aGXOIr1g7hyfvvb32Zubg6n08l7772nPg/dTKXZjhZ6zf/f3plHx3mX9/7z0ywazaJZNNpXy7KtOI5xCAkpCTSkCemN7yGH01Lo5RAOl7AcyiE5JOWW0pay/UEP0NsLaW6AwG1oLmFp2tC0PVxKgmkaO9gJtiM7si3JI2tfZtUsmtHy3j80z8/vKNocL7LH7/ccH0mvZuZ9f+PR83ve7/N9vk8JzF7wMo7PHOBXM+qS44uLi1qOGI1GS+a2mu2EJycn9dDumZkZPcN1dnaWdDqtXxOW6B+zHUIgECAUCumBH8KTi/nZwsICXq+X7du343Q6GRkZ0eccHx/HMAzNuft8Pp3JB4NBPblKrJDNHj1yNyIS0kudOUsRWja2zs5O7rjjDvx+P9FolJ6eHmvEoYUVYQV6CyUwT3c6l6AhAXFubk6bfSUSCU2viCY+kUiQTqcZHx8nnU7rwJ5MJkmn0+RyuZJh4eZNRmgap9NJZWWl7sAVGkWuvVAoaDM0sUYQRY04PlZVVWkZZiAQ0LNo5XzLTcnM9YqFhQV9V3CpIZp/WNqUW1tbmZiYoLa29nXN9bVwdcAK9BZKIIHC6XQyMDBAb28v1113nR6xJ9nzcoh//MmTJ+nv7yeZTGrLXjEzk0KrDPCOxWKaNjEPIZEGKji74djtdlwuFy0tLWzdupW2tjba29vx+/3aCE2KxBIIfT4fDQ0NtLe3a/8cyeglSG7fvp3u7m4aGxvxeDwsLCysyHObaxNyHtlcpI4gG50Uny90xi8blHwPsH37dkKhENPT08zMzDA8PMz4+LhF4VgogRXoLayIxcVFHZC7uro0pbEWcrmc9oPJZDKkUilyuRyZTIbZ2Vmmpqa0skaoG3MwlAw9FApRXV2N0+nUTVCSyYdCIRoaGvQ/cds0q4UEFRUVuN1u2trasNvtugkKluig5uZmtm7dSmdnp27qOtfgKIVaoaVmZ2f1nYJQLRcTdrudhoYGfu/3fo9sNsuxY8d47rnnmJiYsIK9BQ0r0FtYFQsLC9pFciOFR5EsiuGYyDGlQUnUMjLEQzJrMy0ijpTCwdfW1hIIBHTjksfj0RuBuDguLCyUSDoFQq9UV1czOzur9fBOp1OraJqamkruUiRD30iAFtoI0JOv8vl8yRDyi62vl/e4rq6Oa6+9FsMw6O3t1QVtK9hbACvQW1gFDoeD+fl5IpEIyWQSv9+/ZtCSAuZyv3lxi8zlciSTSeLxuB4QYpZOSiFXKUUwGNTDQLZv305NTY1W88hGsLi4qAN4KpVicXHxNbYAYl0QDofxer36OmVEoDhdwlkOXq5jIwFairRin1woFLRR2qUqisp1VlVVceedd/KmN71Ja+rHxsaYmpq6JNdh4fKGFegtrAgJVIVCgampKdra2tZtEHI4HIRCIRKJBC6XSxuEmRusUqmULraKmyQs1QQKhQINDQ20traybds2mpqauOGGG2hvb9fdrFK4dTqdmt9Pp9MopfB6vSVSUPHHEQtjh8OB3+8HznLcEtSXq4rMnL/5sebisPDyhUIBm82Gy+XSHbWrOWheaEidwDAMPB4PHo+Hd73rXUSjUU6dOsW+fftIJpMl1JZVtL36YAV6CyvCPKxjObWyUtOUFGPFRlc4deHdzbLN9eB2u/H7/QSDQS2BFGtgsTCYnZ0ln8+XNFXF43E9i1aKt/IcoVlEa24O7ubuWlmLucAqFJB5wxKpJZwdhQiU3ClciJmxrwcdHR3s3r2bxcVFjh8/vmqHsYWrB1agt7AihFKRblbz8ZUsEJRSVFVV0dLSoqkayaplgHY+ny9R16wWcNxuN42NjbS0tOjO1IqKCrxer+b+c7mcHlcoQ7sBvdl4PB4qKytxOp3U1NToaVm1tbV6A6isrCwZgRiPx5mamtJNXDMzM7ohy+/309DQoP16JCs22zIvLCzoc27mcBCXy8Xtt99OfX09w8PDTExMlFhAXAyvfAuXN9YN9EqpVuBxoB4wgG8ZhvE3Sqm/BD4MCAn4p4Zh/GvxOZ8BPgQsAJ80DONnF+HaLVxECKXh8/lKuO+1skGbzUZ9fT2Li4taPy/OlXLMzF1vJLNcbnGQy+WIx+Mlk5/Es0aye/GtEeWONEO53W791el06rsDKRBPTk4yNDRENptlamqKTCaD2+3G4/EQDofZsWMHdXV1dHZ20tLSoq/JZrNpW2S5c1lu23CpITbQ/f39nD59mkOHDun30wr0Vx82ktHPAw8ahvGyUsoHvKSU+nnxd39tGMZXzQ9WSu0E3gtcCzQB/66U2m4YhkUKXmGoqKggEAhomeBaQUJ+73A4CAaDVFVVae262+3G5XKVNELJc1aC6O0dDoeeQCUNVz09PRw9epTh4WEikQjDw8Pk83ltcCZaf1HmOBwOHfRdLhehUEiPHxQKKJVKkUgkGBwcZHh4WKuDDMMgFArh9XppaWlhfn6eaDSq3TJl+IdZT29+79Za48WEbNJer5dbbrmFnp4e+vv79WB1q3v26sO6gd4wjDFgrPj9jFLqVaB5jafcAzxpGEYeOK2U6gNuAvZfgOu1cIkgtgNCVYhf/Gqt/xJcxJFSOk99Ph+VlZV6YIfZgdEccGQTGB8fJ5FIMD09rZ8j9geyAUxPT2s7g7GxMebm5vRj7Ha75toLhYJuvspkMjidTlKplObwFxcXSaVSxGIx3TMwPT2t5aEOh4PZ2Vkd1NPptJabirmYYCVp52ZB3lu73U5dXR3V1dU0NzfT19dnyS2vUpwTR6+U6gCuB14EbgE+oZS6FzjEUtYfZ2kTOGB62jArbAxKqY8AH3ldV23hksDpdGprAFibhjDTK0L1+P1+ampqqKurI5/PMzExgd/v18NAVsPQ0BDHjx/XIwSHhoaYn58nkUhw4sQJent7GR8fZ3R0lHg8/prnK6XIZrNaRplOp0sKpHKt8pqJRIJkMvmaa5qfn2dmZgbDMAgGgySTSdxud8lmcrlCbBoaGxvZu3cvkUgEu93O4cOHN/vSLGwCNhzolVJe4B+ABwzDSCmlHgG+yBJv/0Xga8B/3+jrGYbxLeBbxde2ZACXGUSXHggEdOBeqQgrMKtUBIZh4HK5CIfDFAoF2tvbmZqaIh6PMzIysiqFMDIyoumXZDJJIBDQWfnIyAiDg4PaHE3OLXcTonQRWmZxcRG73a79dqqqqvR5xWdnenq6RHpo1vTLMZ/PRzgcpqmpiZaWlst+Zqu5k3n37t3cfPPNVFdXMzIyQjwet+ibqwwbCvRKKQdLQf4JwzCeAjAMY8L0+28DzxR/HAFaTU9vKR6zcAVBCoti3LUeFSEctUgeAR1kPR6PHsUXCARYXFxkampKZ53LUSgUiMVi2ho4Ho9rvboYpUlzklybufFKvso1ydATmQtrrjeYzy8WxALZ7LxeL6FQiJqaGsLhsJZ7rrXxbSbMvQHSINbd3U0ul6OpqclqoroKsRHVjQIeA141DOPrpuONRf4e4F1AT/H7nwL/Vyn1dZaKsduAX1/Qq7Zw0SHZudAxYtu7GqQgmkgkmJ2d1cFXZJZiXyCDPUZHRzW3vhwLCwuMj48zNTXF+Pg4jY2NALqAOj8/r3Xx5mAvwU2knqFQSBdOpVbg9/t10TaVSlFZWcnc3Jz2qRGIkqapqYm2tjZ27tzJrl276OjooK6uTp/zcoTZ/0c6d2+88UaUUkxOTjIxMWGNILzKsJGM/hbg/cArSikh+P4U+EOl1B6WqJsI8FEAwzCOKaV+BBxnSbHzR5bi5sqDcNiS4UrxcjWkUiltZib6c5vNpjcAc6FUXnOlIG/GwsKCDsB2ux23213SABUOh8lms1rPLncgbrebmpoa7S9fX19PQ0MDXq9Xu13m83mSyST9/f0AxGIxfS6n00kgEKC5uZnOzk5aW1t5wxvewO7du6mrq1szwJutkoVeWlxc1MVbuUMyc/wulwuv16s3UhmKYg7YcG4bi9k90+FwUFNTw1vf+lby+TyDg4OcOnWKnp6edV7FQrlgI6qb54GVqnD/usZzvgx8+Tyuy8JlADEMM/vArIX5+XlSqRTZbFYHXmmeyufzWs2Sz+fJ5XIbVoAITeNwOHC5XLrRSZqjhEaRDlYZC9jQ0IDL5Srh1auqqnA4HOTzeTweD7lcjlQqhcfjIR6P60AfDAZ105ZYI4dCIW3tYL7bMcMwDPL5POl0moGBAaLRKAsLC+TzeV1DsNlsuqnKZrPpWojD4cDpdOL1evVQlQuheReKye12s2fPHjo6Opifn6evr6/kLsZC+cLqjLWwIioqKrTrpNgPrBWUq6qqqKqqIpPJkMlkND2Sz+dJJBLMzMzoBqfR0VEymcw5dY8uLi5qC2Cfz4fdbsfr9ZZkyPJPhnnX19fjdrtpbW2lublZF1QrKip0U5Tf78fr9WqJ5ezsrNbeb9myhWuvvZbGxka2bt36GpfLlRCPxzl69CiTk5McO3ZM1yJmZma0hHR+fl4PIXe5XASDQWpqavD5fIRCITo6Omhra9Ndt+cLsx9PW1sbf/AHf8CLL77I6dOnOX78uK6lWAXa8oUV6C2sCKFuNhroXS6XDlbpdJpMJqOz5bGxMRKJBH19ffT29pJOp89Jz22321lYWNDXIUqg2tpagsEgDodD6/sdDgc+n4/6+nra2trw+Xx0dHRoT3uhQrxeL1VVVZpGSiQSZDIZ7S8fDAbZsmULnZ2dWqYpgXAtKuXkyZMcOHCA0dFRenp6OHnyZElGL52zUuT1er36ziEYDJLL5QiFQprmEjXT+VoqmFU4N9xwA263m4GBAd07YAX58oYV6C2UQAKKSBFlxJ9QJnNzcyuO0FNKEQgE6Orq0vK9aDSqm5qmpqZ0Gz6c24APCUIVFRVkMhk9HrC5uZm6urrXeNh4vV4aGhpoa2srWRMsBWlz1l5bW8v09LT28JE7go6ODrZu3aq5c7EkhtIAb5aUDg8P09fXx8mTJxkfH+fEiROMjo6u+B4nk0nC4TCAnsQl7pfV1dXU1NToc10IczSzPYPdbmfHjh3cfvvtHD16lJMnTzI6Olqy+VoOl+UFK9BbeA3MQSGXy2k9OrBqZmnmgV0ul24+Eq16KpU67+sy1wwkizfPfRVvG6/XSzAYLNlMzBYFwo8LHSTrrays1JOh5LXWCrLmISWiMhL7hJmZGTKZTInDpcgdzcFUunQrKyt1lu92u8/7vVoPUqQOh8MMDw/r61kuU7VQHrACvYUSmK0J5ufnGR4eJhaL0dTUtOIs1eWQUX+yOSQSCUZGRpiYmFjzeWtB+GMxIBMDs8rKSsLhsFbISLHV5XKVqISWDzyXzUosDWR2rWj+g8Egzc3NeqrWWmuW4C0KpUKhQCaT0ddptns2Z8h1dXV0dHQQCoVoamqio6ODmpoatm/fTkdHx+t+rzaKhYUFdu3axR133EFlZSWxWIx4PK4zeSvQlxesQG+hBOYBHHNzc7z88svs2bOH6667TjcUrdYoJFxyTU0NiUQCp9NZEgRfL8yZ+ezsLPF4nJqaGq3tF/WK0EtCLa109yGKm0QiwcmTJ+nt7dU2B263W3vLmzP5tbJ6OS76feH/5U5DZtQKfD4fwWCQnTt30t3drTcVqSe0t7dfkozeMAy8Xi9vfvObWVhY4NixY8TjcYuuKVNYgd5CCST7la8DAwOMjIzoIuV62a10k0oBE9AKnNcL83OFHpGB49lslsrKStLptJ5SJY1Cbre7hFaBJZvj6elpEokE/f39jI+P6ylYItEUnb450K8G4dDNc2IDgQCzs7P4fL6SQG+326mtraWpqYlt27axY8eOklqDy+XSIxsvdkYtm3V9fT1btmyhpaWF/v5+3QNhGZ+VF6xAb6EEMmdVpID9/f08//zzfOADH9ABcCWYuXO73U59fT0ejwe326218yvB7EMjSpN8Pr/qNCpx0YzH4wwNDeHz+UgmkywuLhKLxXTTlMxuhbPadhk7GIvFyGazTE5OEovFtPZeRiWKjYN5dKL5qzmLl/cjk8no7Ly9vZ2qqiptpiYNUy0tLezYsYOWlhZuvPFGrrnmGnw+H3V1dSWbyUoqG7P3zoWA3Jk1Njbi9/s5fvw4Q0NDHDlyxAryZQgr0FtYE/l8nn379vHCCy9w1113AStz1uYg5PF4tKolGAxSX19PNpslmUzqx4laZX5+XmeXwsOvBLO1gah/YrEYkUiEqampkhqA6MZzuZwu4Mp8WfHJWVxcpFAoYBgG1dXVOhgL1bRcQmm2FBCYB6G7XC4Mw6C1tZVsNksikSAYDNLd3a2vv6mpic7OTrZs2UJ3d7ceWL5cxWMO8hu5q3g9kGHrCwsLuN1u7rnnHsbGxrQjqMXTlxesQG9hXUSjUU6ePMkdd9xRIjM0wxyIZHpTOBzG7XbjdrupqqrSlsJmt0mzGkV+ZxiG9rBXSmmVTEVFBS6Xi6qqKp21ZzIZXQSVMYPShSuKmkKhoEcimjtBRWkj5zY3XZmzdjhLdcgGZZ4ZK92uAMFgkKamJh3EZUh4VVUVjY2NtLe3l/wezm4isvFJF7DZMuFCZ/QiJ5XXbG1tZefOndTW1pLL5bQzqIXygBXoLayLubk5nnrqKfbu3cuWLVtWfIxIFCUjtdlsXHfddVpLvri4SDKZJJPJlGSpi4uLOog1NTXh8Xhobm4mEAjoAC+0kAwXkc3GMAzS6TQzMzNa5SLGahKsJFOXc1RWVmptvMfjobGxEYfDgdvt1oPIhds3Z9c2m02fR9YgryXNWDKRKhQKYRgGc3NzJXcscNar36wAko1lfn6e8fFxIpEIyWQSl8tFe3s77e3tVFZWXtAMWwK9WVt/22230dPTw8svv8yBAwdWpdssXHmwAr2FdWG329m/fz8vvPACTqeT5ubXDhiTwGHOdsPhMG984xuZmJjQ053EX0W4bp/PR1NTE11dXWzfvh2fz0dLSwu1tbV6WAgsbTbJZJJ4PK6tiuPxuG7okq+ZTIbp6WndfSvnkc0nmUwC6MJnIBCgtbWV+vp6Ojo6aGlpwe/36/PKHcLp06cZHR0lmUySzWb1APK2tjY6Ozu1FNM8QlAycvNEK4F08cJSRp9MJhkbG+PgwYNEIhFt1ub1ernhhhtoa2tjz549F+z/1OFw6I1G3p/Ozk7uuusu/H4/kUiEwcHBC3Y+C5sLK9BbWBfCXe/bt4+amhpaWlpe8xjJWpfPTr3mmmuYmJjA5/MxPz9PIBCgUCgwOztLKBSirq6O1tZWrrnmGnbs2IHX66W5ubkk2Aq1Itn6zMwMo6OjnD59mv7+fqLRKDMzMzqDF3pmeWAV/t/pdGoLhYaGBrq7u2lpaaGrq4v6+nodrOfm5hgdHSUajXLixAlOnTqlDdqkKWthYYFAIKDfE3PWbfbGMXe4mjdFmYZ16NAhBgcHOXr0KMePH2dmZob5+XlqamqIxWK0trZit9vZuXMnNpuNQqGwrgpqLZgpOPP/1w033IDL5WL//v0MDw+XKJYsXLmwAr2FdSEBamxsjLGxsfWfYILP52PLli3Mzc3phpx8Pk+hUKCuro5AIEBbWxtbt26lq6tLOzeag4/w1OJPEwqFCAaD2O12stksSimmpqZIJBLk83nN35sVPeaiqcfj0WMOw+EwjY2NNDY2EgqFSgqh+XyeWCzG5OQkg4ODRCIRPYS8urqamZkZAoHAmmMRV4J0oMoaY7GYfv0zZ84wMjKizyNzeysqKujt7WXLli14vV4d6C80vF4v9fX12nPfQnnACvQW1oVknqdPn2ZwcFBr1c00xUoQaqCrq4tAIIDX62V4eFirXurq6vB4PGzdupW2tjZNf0Dp4JPl06KkwOv3+6mvr2doaEh3xE5NTWEYhvaXl2xUaKD6+npaW1t109J1113HNddcQzgc1gVSGQw+MzOjKYxXX32Vl19+mVQqxfz8PI2NjXqgydatW9myZcuKHkCrQeoZ8XicAwcOcOzYMYaGhhgYGGBoaEjTW6lUCsMwNGff0NDA9u3bqa+vfw3/f76QzaOhoYH6+voL9roWNh9WoLewJsxeLplMhng8ri1vBau5K0rgc7lctLS0EAgEtIGYNBTZbDbtQCna/ZVee3lmLt70O3bs0A1H9fX1jI6O0tDQwOjoqKY/FhcXcblceqBIa2sr4XCY3bt3s3v3bm0JLOcSyiiVSjEwMEBvby+/+tWvGB4e1vRPLBbTxdu2tjZ27dpFIBDYkDJGCrUnTpxgbGyMV199VdtEpFIp/f7KsBKRiiaTSU6cOKE3rAstuZTNezn9ZuHKhxXoLawJ87g+sduFjZleSeCUgCVqFxkxKE1K0lQkkkIz5Lmy4Yg0U66hoqJCF3RbWlqoqKjQ+nmv16vpDzlPMBikrq6Ompoaamtrqa2t1a+Xz+epOqNFtAAAFglJREFUrKwEzvrxS6FXZtiKF0w+n9eF4VQqpaWdG8mwJXhPTU0xPT2tLYyFzjG/56IakjqJnFvWdKGVOICWo1ooH1wWgV7+OBwOB7Ozs9rHxOv1MjY2hsvl0gU24VGnp6dxu91ks1lqa2tZXFwkGo3qx8okIPMHtrq6WhcC18NyfxO/368VGzU1NeTzecLhMJFIRAequbk5HQjEHx2WiohyrfJ8kRvC2T8wv9+v/Vsymcz5vq0XBOYCpmTUiUSipCBrHltnhryHy+WF4rMCZwOZSPzgLFUEZztnodS7Zrl9bzgc5tZbbyUej9PS0sKZM2e0l768pzLNqbGxkZqaGrq7uwkEAvo1xJ9egm5fXx99fX0MDAxo900pDEtx9+jRo3R2dvLKK6+wa9cuGhoa9OdgtS7ibDbLsWPH6O/vJxaLkclksNvt+P1+ZmdndfevUgqPx6OVQDL/dXJykpmZGc3dm62SDcNYkbtfvoGm02mqqqpK7pYcDgdTU1McPXqU3/zmN5bvTRnhsgj0EhhzuRy1tbV6Sv3MzAyNjY26AOhyuUin0xQKBWpra3UDTiKR0H9U4mooRTm3243NZtOzTM2zOpfD7MEtc0UF27Zt49ChQ9TX1+suzHw+T0dHB3Nzc4yMjABLG5HT6WRkZASPx6M3FXMgjEajJedtaGhgZmZG87Bzc3N0dHQQiUQuiBf5hYBsZOeT6UmBVPTz8/Pz5HI5nbU7nU7sdnsJV7+a//1KEIvhhoYGUqmUNiuTZilAq2Vk/qzA3PzkcDgYHBykt7eXnp4eTp8+veL5FhYWiEajvPjii+zcuRNY2nCWSynlTkCGokciEXp7exkZGaFQKGgv/6qqKtxuN/X19fpzY7fbCQaD2j9IVDpnzpzRcsvldxErfV7keqRvQe525P9Egvzjjz/O888/z8GDBzf0nlu4MnBZBHpY+gB2dXXx2c9+Fo/HQ3d3N29729sYGxsjEAiQTCaprKwkmUzS2trKN77xDXbs2MHQ0BDveMc79B+tORDNzc2VDKB2u92rZik2m40dO3bw4osvcuLECV566SU+85nPEI1Geeyxx7jlllvYu3cvAwMD7Nmzh4cffhibzcYvfvELHnroIX75y1/qLPATn/gESikymYx2REyn07z//e/n8ccf15nY+973Pj7+8Y/z3HPP8clPfpKWlhYee+wxstksbrebv/qrv+KXv/zlZZFZiczuQl1LPp8nGo2SSCS0FUFVVRU2m426ujqd8Z8LPWHm2EOhELW1tZoPNxd3zYVkKLUxluPT09NMTk5qJc9KkMQgHo8zPT3N1NQUhUIBt9utX9Pc8Ts7O6vlmslksuRuSe58ZKMz33WK7bJ0AldUVDA5OUk6ndYKpbXmBcjvpRPZ6XTqzldzAfrQoUP87Gc/o7e316JuygzrBnqllAv4FVBZfPxPDMP4nFJqC/AkUAO8BLzfMIyCUqoSeBy4AYgC7zEMI7Luhdjt9PX18dWvfpVjx47R0NBAIpGgqqqKVCpFbW0tk5OTAPT19fGe97yHn/70p9x///3AWdmczOeEJc9veQ6w5of3a1/7GidOnMDn8wHw7ne/m8bGRgA+9alPkclk2LZtG1/+8peJRCLccssttLW1ccstt/DhD3+Yxx9/XN99/Md//Ie+QxF3RbvdzhNPPMHjjz9OOp3G5/ORyWS47rrr8Hg87Nmzh8OHD3PvvfcSj8fxer2aitrMQG/Wfi/fOM8VIoecmZnh1KlTjI6Okkql9F2a2B5UVlbS2NioFTmtra0ben1zhm7elIT+EXrI7Ggpfi9yJyhjD/ft28fBgwe1c+dKkMA5PT3NwMAAADt27OD666/Xj5Hi5uzsLEeOHKGnp4dkMsnk5KRWJYnfjNBR4stj3qRg6c5gdHQUr9fL7OwshUKBXbt2UVNTg9frXbWAKgZmor9PJpNUV1czNTXF4cOHiUQiPPvss+zbt0+vw0J5Qa2XLamlT4/HMIy0UsoBPA/cD3wKeMowjCeVUv8bOGIYxiNKqY8Duw3D+JhS6r3AuwzDeM865zAkoPh8PkZHR7HZbHg8HgzD4Itf/CJ//ud/rvnIt7zlLXz605/miSee4Mc//jEA3/72t/nwhz+MYRj4/X7e8Y53cOrUKQ4cOEBVVRUf/OAHeeqpp/jgBz/I9773PU2TzM7OYrPZyGaztLe3Mz4+XpJRyR+bZIJzc3O0tLQwMTGhMzq5ru7ubnp7e1dcYygUIp1Ok0gktOlXoVDg4YcfZseOHdxxxx0AbN++nZMnT/Ktb32LX/3qV/z93/89sBTEzNOJLiXMhdDf/u3f5pFHHqGzs7Nkjupq/jcSZAxjaVpVT08PY2Nj9Pb2Mjg4qDPTXC6ng5G4QPp8Prq6uti9ezehUIi2traSzPz1YDWTMKUU+Xye//zP/2R8fJx/+7d/45lnniGRSKxKn8lr2Gw2brrpJrq6uvid3/kd3vrWtxIOh3XSkE6n+fWvf83hw4cZHh7WxmqSxTudTk1bVVRU6OEl4s0js2zNG4DX68Xv97N161aqq6u1RLWyslIPT5eidywWI5FIkEgkGBgYYGpqiqGhId1wJm6hl8Odo4VzxkuGYbxpvQetm9EbS59wcThyFP8ZwO3Afyse/zvgL4FHgHuK3wP8BPimUkoZ6+wo8mG/9tpr2bJlC3feeSf/8i//wt13383XvvY14CzP+MILL/CpT32Kf/7nf+bgwYNEo1Huu+8+7r//fvx+P06nk/3797N37149bPno0aMkk0m+853v6NtWCeZyKytc8OzsrFZfwFLwDYfD/OhHP8JutzMxMYHH4wGWsiyzWmP5bTecvW0vFAq4XC6qq6t1Ifa73/0uH/vYx4AlaikUChGPx/nOd77DE088oYddZ7PZ88qmLxTMhT8oLZyu9FhRkORyOQYHBxkYGGB0dJRTp04xMDCgi5FmbX46nUYphc/nw+Px0NDQgGEYNDU1lRhxnS/MTprz8/PaBdMsc5THrfXxFZsEsWUYHx+nUChQU1OjxyiOj4/raVbZbFZvkCITNRdw5+bmdO0iFotpSigWi+kNQrp7+/r6cLlcWkkkpmzSlDU7O8vs7CzJZJJUKsXQ0BCxWKxk44CVN2oL5YMNcfRKKRtL9EwX8DDQDyQMwxBP2WFADFCagSEAwzDmlVJJluid6bXOIcVPKVD94Ac/4HOf+xyApmIAuru7+fSnP819993H9ddfTz6fZ9u2bUSjUZRSJbNJf/KTn/Doo48Si8Xwer3U1dXhcDhKnPmEYnnooYf4/Oc/z4c+9CEAnn76ae6++24AkskkkUiEG2+8kXvuuUffPQQCAb761a/yx3/8x3qA9kqKHo/HQzqd1kFKpH7XX389L730Eg8++CBdXV0ADA4O0tHRQTabJRKJcO2115JMJjf1D9Gs6JDAYQ6SqxVL5XkyIPz48eMcPnyYkZERXnnlFd18JZmkOFQODg4yPT2N3+8nk8mwuLhIU1MTjY2NNDQ0lLz2uUJUJlKUhLP1gsOHD3Pw4EGmp6fp7+/XXbfrvTeGYTAyMsLCwoK2eqiursbr9ZJMJkkmkxw5coRXX32ViYkJnZWLEkg+D3CW9xc7B6H/zG6Z5vVHIpFzWr8YxUmBWDxvrGy+vLGhQG8YxgKwRykVAP4R6D7fEyulPgJ8RH5Op9O84Q1v4KMf/SjV1dWk02n27NlDW1sbX/rSl7j33nux2Wz09vYSDod59tlnmZiY0B/+rVu3Eo1G+bM/+zP+9m//ljNnzvCWt7wFgAceeICbbrqJb37zm3zhC1/gBz/4AY899hiVlZXkcjmcTiff+973ePLJJ9m3bx8VFRU8/fTT3HzzzQD8+Mc/JpfLsXv3bmZmZujt7eXRRx9l165dmnN/73vfy1133UUoFOKee+4pWasMoHjggQcA+MIXvsDnP/953vnOd/Lggw8SiUT4yle+Qnt7O9/4xjc4duwYDz30EG9/+9u57777eOSRRzbNSdBM2xiGoR0izSqn1YKhzWYjl8sxNjbGmTNnOHLkCIcOHWJ4eJjR0VGtODHD7OI4MzOD1+ulo6MDt9utN+jz6QiVuy/JnvP5PK+88grj4+McOXKE06dPk0wmmZ6e1mtYLQiaN5upqSn9vsTjcR1AR0dHSafTnD59WstrJds+183K/Pi16KTlx80NUOY6ldn7X4q9FsoT63L0r3mCUn8B5ID/ATQUs/bfAv7SMIy7lFI/K36/XyllB8aB2rWoG6WU4fV6SafT2g9bKaW1wuaMXrLmyspK5ufnXzOuzQyRNzY3N+s/unA4zPT0tPZMF58Ss05eMlSR/Ykmf3mwDQaDAFrmucb6dPZWKBTweDwlOnlzMAkGg6TTaX1LLbWBzVJBSP1BAkh3dzePPvoot956q6ZZzCoTM5RSzM/Pc+DAAU6cOMHPf/5zXnjhBcbGxkq4fcmKzfy5UbQBvuaaa7jtttuora1l79697Nq16zXa/HOBnCOfz9PX18fp06fp6+tjcnKS3t5eent7SaVSTExMsLi4qOsQq9VGJIgKTWUO4DLaULA8mMoaVjIOW60Gsfxx8v6tdn3mGcBrncf8/lu4onBhOHqlVC0wZxhGQilVBdwJfAV4Dvh9lpQ3HwCeLj7lp8Wf9xd//+x6/DygszXR0AvMQR7O8uoSdFcL8oAOpmfOnNHHJFMz/wF6PB4d5KE005Ei4kpYL8ALqqqqyGazegNxu90lgd6cMS5/zeVNX5caErwA3XVqLoba7Xby+XyJ9l0gG4TL5dK0jHjFCFbKlperY/L5vOadzRvP8sak4eFhMpkMuVxOywilxiHqFulkzWazjI6OMjIyQiQS0dOqotGo9uKBtZVacq3y8V6+luXPXZ4xr0WXbLTovh7lst7rWM6UVwc2Qt00An9X5OkrgB8ZhvGMUuo48KRS6kvAb4DHio9/DPi+UqoPiAHvvQjXfd6Q7CoQCKy5WZzv6wM0NjbS399PMBhkfHxcbzZXGmTTMReFheddaS83Z4vmcX6wfoFT7qakKJnL5YhGo1o7Dmc91SUbP3XqFLFYjImJiZLRfzLmT8YEzszM6A7TqakpTbuIj48FC+WGjahujgLXr3B8ALhpheOzwLsvyNVdREhLfyKReM0t9oU+D6DVOus5Pl5uWH6dopCR35nnq64GcZaUZiYzXbPeuaPRKJFIRKunAoEAjY2NWik1MDDAmTNnePnllzl8+DCxWEyrXuAsbSa6erFEkPPLhrIR/tuChSsVl01n7KWEuZHK7Jp4IeH3+zUN0N/fD6DrA1ciJPDJgG0o5X1X4pSliCuZtHDd58IHZ7NZYrEYNpuN0dFRhoaGdINRLBajv7+f4eFhzpw5Q29vrx5A8nrWttrPFixc6bgqA73QBxczk8/n81RUVGi1iNfrvWKDvBmFQkF3izocDk2PrOavopQiHo9rjfrymstqEI4+l8vx0ksvEQwGicfjjI6OEgwGMYpe7adPn2Z4eJi+vj793LUK2JLdrwSrMGmhXHFVBnopeF5MPlayyuVfrzRIm74ER7vdzqFDh0gkEoTDYW07YNbVS4HQ4XBw5swZjh8/zqlTp/SdjbzOWnI+8+8qKipIJpMcPHiQQ4cOrajwMb/eWgXstYqPVmHSQrni9fWRW7iqYNZh5/N5enp69OBokVDC2UxYpKQLCwuMjIxw6tQp+vr6SpRN5wJzg5OZBjL7tFsZuAULq8MK9BbWxXIJ3+HDh/n6179OT0+PthGWoCt+NYZhsH//fr7//e/zT//0T+zfv19n2dKdulHY7XZtb7warEBvwcLquCqpGwvnjuWeNj/84Q9JpVJcf/31vO1tb6OtrU1r3o8ePcrU1BTPPPMM+/fvL7GFMHfUbrQb0yzJlBZ+eQ0ZWmLRLhYsrI5z7oy9KBeh1OZfhIUNQ2gZaYby+/1s375dN0RJu/9m2TZYsHAVYUOdsVagt3DOMFs2yHg+cZy8HD5PFixcRdhQoLc4egvnDDNnL06L5qYpp9O5oiWCBQsWNgcWR2/hdcNsygVLAV5UOJYTogULlw8ul0CfBk5s9kVcIoRZx5v/SoE5s5fC6DKUzVo3gKtlrVfLOuHKWGv7Rh50uQT6ExvhmcoBSqlD1lrLD1fLWq+WdUJ5rdXi6C1YsGChzGEFegsWLFgoc1wugf5bm30BlxDWWssTV8tar5Z1Qhmt9bLQ0VuwYMGChYuHyyWjt2DBggULFwmbHuiVUr+rlDqhlOpTSv3JZl/P+UIp9V2l1KRSqsd0LKSU+rlS6lTxa7B4XCml/ldx7UeVUm/cvCs/NyilWpVSzymljiuljiml7i8eL8e1upRSv1ZKHSmu9fPF41uUUi8W1/RDpZSzeLyy+HNf8fcdm3n9rwdKKZtS6jdKqWeKP5flWpVSEaXUK0qpw0qpQ8VjZfcZ3tRAX5xD+zDwX4CdwB8qpXZu5jVdAPwf4HeXHfsT4BeGYWwDflH8GZbWva347yPAI5foGi8E5oEHDcPYCdwM/FHx/64c15oHbjcM4w3AHuB3lVI3A18B/towjC4gDnyo+PgPAfHi8b8uPu5Kw/3Aq6afy3mtbzcMY49JSll+n2HxEN+Mf8BvAT8z/fwZ4DObeU0XaF0dQI/p5xNAY/H7Rpb6BgAeBf5wpcddaf+Ap4E7y32tgBt4GXgzS8009uJx/VkGfgb8VvF7e/FxarOv/RzW2MJSgLsdeAZQZbzWCBBedqzsPsObTd00A0Omn4eLx8oN9YZhjBW/Hwfqi9+XxfqLt+vXAy9SpmstUhmHgUng50A/kDAMQ9qBzevRay3+PgnUXNorPi/8T+DTgHg/11C+azWA/6eUekkp9ZHisbL7DF8unbFXDQzDMMrJrVMp5QX+AXjAMIyU2bO+nNZqGMYCsEcpFQD+Eeje5Eu6KFBK/Vdg0jCMl5RSt2329VwC3GoYxohSqg74uVKq1/zLcvkMb3ZGPwK0mn5uKR4rN0wopRoBil8ni8ev6PUrpRwsBfknDMN4qni4LNcqMAwjATzHEn0RUEpJsmRej15r8fd+IHqJL/X14hbgnUqpCPAkS/TN31Cea8UwjJHi10mWNvCbKMPP8GYH+oPAtmJF3wm8F/jpJl/TxcBPgQ8Uv/8AS3y2HL+3WM2/GUiabhkva6il1P0x4FXDML5u+lU5rrW2mMmjlKpiqRbxKksB//eLD1u+VnkPfh941iiSupc7DMP4jGEYLYZhdLD09/isYRjvowzXqpTyKKV88j3wDqCHMvwMb3qRALgbOMkS5/nZzb6eC7CeHwBjwBxLHN6HWOIsfwGcAv4dCBUfq1hSHfUDrwBv2uzrP4d13soSv3kUOFz8d3eZrnU38JviWnuAvyge7wR+DfQBPwYqi8ddxZ/7ir/v3Ow1vM513wY8U65rLa7pSPHfMYk/5fgZtjpjLViwYKHMsdnUjQULFixYuMiwAr0FCxYslDmsQG/BggULZQ4r0FuwYMFCmcMK9BYsWLBQ5rACvQULFiyUOaxAb8GCBQtlDivQW7BgwUKZ4/8D8w1vWzSvwU8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "im=itk.imread(\"data/CBCT-TextureInput.png\", itk.UC)\n",
    "arr = itk.array_view_from_image(im)\n",
    "plt.imshow(arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can modify a value in the array and look at the value of this pixel in the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original image pixel value:146\n",
      "Modified image pixel value:0\n"
     ]
    }
   ],
   "source": [
    "print(\"Original image pixel value:%d\"%im.GetPixel([270,270]))\n",
    "arr[270,270] = 0\n",
    "print(\"Modified image pixel value:%d\"%im.GetPixel([270,270]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Views vs copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So far we have used `itk.array_view_from_image()` and `itk.image_view_from_array()`. You can see the keyword `view` in both the names of these functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's compare the result of `itk.array_view_from_image()` and `itk.array_from_image()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Copy')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADHCAYAAAAJSqg8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAE25JREFUeJzt3X2QXXV9x/H3Z3cTkqy7ZAmg4UnCFNOxtBVmS0FFqVBEdISxTgcrSK1jZFqttp06WJ9Gqp3WaR3aqVVTtUChaEGwaH3CJ9RRU0IEBYLyFCDJQkJM2CUPkt399o97opfNbnb3d889955fPq+Zndx7zvnm973nfve7Z889+zuKCMzMrP56Op2AmZmVww3dzCwTbuhmZplwQzczy4QbuplZJtzQzcwy4YbeBSR9XNJ7O52HmdWbG3oFJH1F0uXTLD9f0mPAWyPibzuQmlnLJP2RpLWSnpI0IunLkl7c6bwORm7o1bgKuEiSpiy/GLg2IsY7kJNZyyT9JXAF8HfAs4HjgH8Dzu9kXgcrN/RqfB5YBpyxb4GkIeBVwNWSrpT0waZ1r5J0h6Qdkr4v6beK5W+U9IWm7e6TdH3T80clvaCKF2Qm6VDgcuDPIuLGiNgZEXsj4gsR8deSDpF0haTNxdcVkg4pYs+UtFHS30h6QtIGSa8v1v2OpMcl9TaN9RpJd3bmldaHG3oFImI38N/AG5oW/yFwb0Q8o0glnQx8GngLjR8CnwBuLr4RbgXOkNQj6ShgIXB6EXcC8Czgx21+OWb7nA4sAm6aYf27gdOAFwC/DZwKvKdp/XOAw4GjgUuA1ZJWRsRtwDbgnKZtLwauLjX7DLmhV+cq4LWSFhXP31Asm2oV8ImIWBMRExFxFfAL4LSIeBAYo/EN8hLgq8BmSb8OvBT4bkRMtvuFmBWWAU8c4JTh64HLI2JLRGwFPkCjMTd7b0T8IiJuBf6XxoEOFKcpASQdBrwc+K+yX0Bu+jqdwMEiIr4n6QngAkm30Thaec00mz4XuETS25qWLQSOKh7fCpwJ/FrxeAeNZn568dysKtuAwyX1zdDUjwIebnr+ML+qY4DtEbFzhvXXAOsl9dNo8t+NiJHyUs+Tj9CrdTWNI/OLgK9GxOPTbPMo8KGIWNr0tSQirivW72voZxSPb6XR0F+KG7pV6wc0fnu8YIb1m2kcoOxzXLFsn6GiYe+3PiI2Ff//a2gc1f9nSTlnzQ29WlcDZwNvZvrTLQD/Dlwq6XfV0C/plZIGivW3Ar8HLI6IjcB3gXNp/Pr7o/amb/YrEfEk8D7go5IukLRE0gJJr5D0YeA64D2SjpB0eLHtNVP+mw9IWijpDBoXCVzftO5q4J3AbwI3tv0FZcCnXCoUERskfZ/GB0Q3z7DNWklvBv4VOBHYDXwP+E6x/meSnqLRyImIUUkPAlsjYqKCl2H2SxHxT8XfUrwHuJbGZzy3Ax8C1gGD/OqD+uuBDzaFPwZsp3FUvgu4NCLubVp/E/Ax4KaI2NXO15EL+QYXZlY1SWcC10TEMbNs9wDwloj4eiWJ1ZxPuZhZV5L0B0AA3+x0LnXhUy5m1nUkfRt4PnCxL8WdO59yMTPLhE+5mJllwg3dzCwTlZ5D7+npiZ6e+f8MSYkB6O3tnX2jEsdLjdt/EsbZpZ4qm5xMOx2ZGjcxkXYlZcp4k5OTTE5Ozn9ntkhS0ptRdZ1V/f1QpTrUZytxNKZYOGK2japu6AwMDMy+4RT9/f2zbzSNlLFaGS81LuUbLbUQd+7cOftGJcaNjY1VNl7qWJ2yaNGi2TeaxuDgYFJc1XWd8oMgteFVXdejo6NJcbt2JV9O//Dsm/iUi5lZNlpq6JLOlfRTSfdLuqyspMw6zbVtdZTc0IvJ5z8KvILG9aKvk/T8shIz6xTXttVVK0fopwL3R8SDEfE08Bl82ynLg2vbaqmVhn40jale99lYLDOrO9e21VLbr3KRtIrGXXiSLs8z60bNdW3WLVpp6JuAY5ueH1Mse4aIWA2sBujr6/M8A1YHs9Z2c12nXoduVrZWTrncBpwoaYWkhcCFzDDHt1nNuLatlpKP0CNiXNJbadyouBf4dETcXVpmZh3i2ra6aukcekR8CfhSSbmYdQ3XttWR/1LUzCwTlc/lkjIvxNKlS5PGO/LIIyuNS80zZU6PPXv2JI21Y8eOpLgtW7YkxaVOCJUidV6OVvX09CS9h0NDQ0njHXHErHM0TSu1rlPzrLKut2/fnhRXh7qGuc8B4yN0M7NMuKGbmWXCDd3MLBNu6GZmmXBDNzPLhBu6mVkm3NDNzDLhhm5mlgk3dDOzTLihm5llwg3dzCwTbuhmZplwQzczy0Slsy329vYyMDAw77jUWeKOP/74pLgVK1YkxS1fvjwprspZ6UZGRpLilixZkhSXamJiYt4x27Zta0Mms+vp6WFwcHDecamzJqbWZ+r3Q2pdL168eN4xu3fvThorta5TZn9tRUpdg2dbNDM76Lihm5llwg3dzCwTyQ1d0rGSviXpHkl3S3p7mYmZdYpr2+qqlQ9Fx4G/ioh1kgaA2yXdEhH3lJSbWae4tq2Wko/QI2IkItYVj8eA9cDRZSVm1imubaurUi5blHQ8cDKwZpp1q4BVAH19lV4ladaymWq7ua57evxRlHWHlitR0rOAzwHviIjRqesjYnVEDEfEsBu61cmBaru5rt3QrVu0VImSFtAo+Gsj4sZyUjLrPNe21VErV7kI+BSwPiI+Ul5KZp3l2ra6auUI/UXAxcDLJN1RfJ1XUl5mneTatlpKPqkdEd8DVGIuZl3BtW115U9zzMwyUellJz09PUmzm6XOtpg6K91JJ52UFHfcccclxaXMQDk2NpY01rJly5LiUs11lripRkf3u2BqVp262qS3t7fSuk6dNbHquk6ZgTLlfYfq63rnzp1Jcanft3PlI3Qzs0y4oZuZZcIN3cwsE27oZmaZcEM3M8uEG7qZWSbc0M3MMuGGbmaWCTd0M7NMuKGbmWXCDd3MLBNu6GZmmXBDNzPLRC1mW1y6dGnSeMuXL0+KS51d7nnPe15SXJWzLabatm1bUtymTZuS4lLqpFOzLabW9dDQUNJ4Vdf1ypUrk+JyruvNmzcnxaXUyXz4CN3MLBNu6GZmmWi5oUvqlfQjSV8sIyGzbuHatrop4wj97cD6Ev4fs27j2rZaaamhSzoGeCXwyXLSMesOrm2ro1aP0K8A3glMlpCLWTdxbVvtJDd0Sa8CtkTE7bNst0rSWklr9+7dmzqcWWXmUtvNdT0+Pl5hdmYza+UI/UXAqyVtAD4DvEzSNVM3iojVETEcEcMLFixoYTizysxa28113ddX6Z9zmM0ouaFHxLsi4piIOB64EPhmRFxUWmZmHeLatrrydehmZpko5XfFiPg28O0y/i+zbuLatjrxEbqZWSbc0M3MMlHpx/OS6O3tnXfcokWLksZLjUuZJa6VOEnzjomIpLFSc6z6PUipk5T9WJaUmR5T983ixYuT4gYHB5PiUmsmdZbUFKmvLXVfpr537Z4R1EfoZmaZcEM3M8uEG7qZWSbc0M3MMuGGbmaWCTd0M7NMuKGbmWXCDd3MLBNu6GZmmXBDNzPLhBu6mVkm3NDNzDJR6eRcEcHExMS84/bs2ZM0Xmrc2NhYpXEpE21VnWPV70FKnaROWFaGycn530s6dd/s3r07KW50dDQpLrVmqhwr9bWl7svU9y6lTubDR+hmZplwQzczy4QbuplZJlpq6JKWSrpB0r2S1ks6vazEzDrJtW111OqHov8MfCUiXitpIbCkhJzMuoFr22onuaFLOhR4CfDHABHxNPB0OWmZdY5r2+qqlVMuK4CtwH9I+pGkT0rqLykvs05ybVsttdLQ+4BTgI9FxMnATuCyqRtJWiVpraS1e/fubWE4s8rMWtvNdT0+Pt6JHM3200pD3whsjIg1xfMbaHwTPENErI6I4YgYXrBgQQvDmVVm1tpuruu+vkr/Ps9sRskNPSIeAx6VtLJYdBZwTylZmXWQa9vqqtVDi7cB1xZXATwIvLH1lMy6gmvbaqelhh4RdwDDJeVi1jVc21ZH/ktRM7NMVPppzuTkJDt37px33I4dO5LGGxkZSYpbtmxZUlyqgYGBecekzkr3yCOPJMWl7svU9y6lTto9k92Bxk3Jd/v27Unj1aWuBwcH5x2TOmti1XWd+t6l1Ml8+AjdzCwTbuhmZplwQzczy4QbuplZJtzQzcwy4YZuZpYJN3Qzs0y4oZuZZcIN3cwsE27oZmaZcEM3M8uEG7qZWSbc0M3MMlGL2Ra3bNmSNN6SJUuS4lJt27YtKW7RokXzjtmzZ0/SWKmzyz300ENJcanvXZ1mW5yYmKi0rvv7q71fdWpdL168eN4xu3fvThorta43bNiQFFdlXc+Hj9DNzDLhhm5mlgk3dDOzTLTU0CX9haS7Jd0l6TpJ8z8ZbNaFXNtWR8kNXdLRwJ8DwxFxEtALXFhWYmad4tq2umr1lEsfsFhSH7AE2Nx6SmZdwbVttZPc0CNiE/CPwCPACPBkRHxt6naSVklaK2nt+Ph4eqZmFZlLbTfXdaculzSbqpVTLkPA+cAK4CigX9JFU7eLiNURMRwRw319lV72bpZkLrXdXNc9Pb62wLpDK5V4NvBQRGyNiL3AjcALy0nLrKNc21ZLrTT0R4DTJC2RJOAsYH05aZl1lGvbaqmVc+hrgBuAdcBPiv9rdUl5mXWMa9vqqqWT2hHxfuD9JeVi1jVc21ZH/jTHzCwTlV52MjExwdjY2Lzjent725DNzHbt2pUUt2nTpqS4Kmdb3LFjR1Jc6uxyqXEpdTIxMZE0VqsmJycZHR2dd1zVdZ0609/mzWmX4FdZ19u3b0+KS63PrVu3JsWl1Ml8+AjdzCwTbuhmZplwQzczy4QbuplZJtzQzcwy4YZuZpYJN3Qzs0y4oZuZZcIN3cwsE27oZmaZcEM3M8uEG7qZWSbc0M3MMlHpbIuTk5PJM76lSJ19L3VGtP7+/qS4lFn3Ul9b6v5PjUuZNTF1vE7drHlycjJ5hs4Uqe996nuRWtcp91pNfQ+rruvUHtHuOvERuplZJtzQzcwyMWtDl/RpSVsk3dW07DBJt0i6r/h3qL1pmpXPtW25mcsR+pXAuVOWXQZ8IyJOBL5RPDermytxbVtGZm3oEfEd4OdTFp8PXFU8vgq4oOS8zNrOtW25ST2H/uyIGCkePwY8u6R8zDrNtW211fJlixERkmKm9ZJWAauKx60OZ1aZA9V2c12bdYvUI/THJS0HKP6d8dbZEbE6IoYjYjjlulSzis2ptpvrutLszA4gtcPeDFxSPL4E+J9y0jHrONe21dZcLlu8DvgBsFLSRklvAv4e+H1J9wFnF8/NasW1bbmZ9Rx6RLxuhlVnlZyLWaVc25Ybn9Q2M8uEG7qZWSYUMeMVh+UPJm0FHp5h9eHAE5UlM7NuyQOcy3QOlMdzI+KIKpOB2tQ1OJfpdEseUEJtV9rQD0TS2m64BKxb8gDn0s15zFU35etcujcPKCcXn3IxM8uEG7qZWSa6qaGv7nQChW7JA5zLdLolj7nqpnydy/66JQ8oIZeuOYduZmat6aYjdDMza0GlDV3SuZJ+Kul+SfvdOEDSIZI+W6xfI+n4NuVxrKRvSbpH0t2S3j7NNmdKelLSHcXX+9qRSzHWBkk/KcZZO816SfqXYr/8WNIpbchhZdNrvUPSqKR3TNmmbfuklbsHSbqk2OY+SZdMt027ubanzaXjdV2Mc/DUdkRU8gX0Ag8AJwALgTuB50/Z5k+BjxePLwQ+26ZclgOnFI8HgJ9Nk8uZwBcr2jcbgMMPsP484MuAgNOANRW8V4/RuPa1kn0CvAQ4BbiradmHgcuKx5cB/zBN3GHAg8W/Q8XjoSretyn7y7W9fy5dVddN71W2tV3lEfqpwP0R8WBEPA18hsbdYZo13y3mBuAsqfxJ1CNiJCLWFY/HgPXA0WWPU6Lzgauj4YfAUhVTvLbJWcADETHTH8uULtLvHvRy4JaI+HlEbAduYf/byrWbaztN1XUNmdd2lQ39aODRpucb2b/QfrlNRIwDTwLL2plU8avvycCaaVafLulOSV+W9BttTCOAr0m6XY0bJ0w1l31XpguB62ZYV9U+gbndPajqfTMd1/b0uq2uIfPabvmORXUm6VnA54B3RMTolNXraPxa9pSk84DPAye2KZUXR8QmSUcCt0i6t/ipXjlJC4FXA++aZnWV++QZIg58Zyx7pi6p7a6pazg4arvKI/RNwLFNz48plk27jaQ+4FBgWzuSkbSARsFfGxE3Tl0fEaMR8VTx+EvAAkmHtyOXiNhU/LsFuInGr/DN5rLvyvIKYF1EPD5NnpXtk8Jc7h5U5b6ZiWt7Gl1W13AQ1HaVDf024ERJK4qflBfSuDtMs+a7xbwW+GYUnw6UqTh3+SlgfUR8ZIZtnrPvHKekU2nsq9K/ASX1SxrY9xg4B7hrymY3A28orgo4DXiy6de1sr2OGX4lrWqfNJnL3YO+Cpwjaai4UuCcYlmVXNv7j9FtdQ0HQ22341PdA3zaex6NT90fAN5dLLsceHXxeBFwPXA/8H/ACW3K48U0zu/9GLij+DoPuBS4tNjmrcDdNK5Y+CHwwjblckIxxp3FePv2S3MuAj5a7LefAMNtyqWfRhEf2rSskn1C4xttBNhL41zhm2icY/4GcB/wdeCwYtth4JNNsX9S1Mz9wBurrGnXdvfX9cFU2/5LUTOzTPgvRc3MMuGGbmaWCTd0M7NMuKGbmWXCDd3MLBNu6GZmmXBDNzPLhBu6mVkm/h8qdCdpaCqDYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gaussian=itk.gaussian_image_source(Size=11, Sigma=3, Scale=100, Mean=[5,5])\n",
    "arr_view_gaussian = itk.array_view_from_image(gaussian)\n",
    "arr_gaussian = itk.array_from_image(gaussian)\n",
    "gaussian.SetPixel([5,5],0)\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(arr_view_gaussian)\n",
    "plt.title(\"View\")\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(arr_gaussian)\n",
    "plt.title(\"Copy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The `view` functions does not copy the memory. It shares the same information. In some cases, it is better not to copy the memory (speed, memory limitation). In other cases, it is better to copy the memory (original image needs to be used later)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Copies and views work in both directions: From NumPy to ITK and from ITK to NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "41a885c7691c4cf2a4ab7642ed4d2abd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Vmlld2VyKGdyYWRpZW50X29wYWNpdHk9MC4yMiwgcmVuZGVyZWRfaW1hZ2U9PGl0a0ltYWdlUHl0aG9uLml0a0ltYWdlU1MyOyBwcm94eSBvZiA8U3dpZyBPYmplY3Qgb2YgdHlwZSAnaXRrSW3igKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "im = itk.image_from_array(arr_gaussian)\n",
    "im_view = itk.image_view_from_array(arr_gaussian)\n",
    "arr_gaussian[5,5] = 0\n",
    "view(im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "52a1e2d28ac745839cbb327af9b374b2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Vmlld2VyKGdyYWRpZW50X29wYWNpdHk9MC4yMiwgcmVuZGVyZWRfaW1hZ2U9PGl0a0ltYWdlUHl0aG9uLml0a0ltYWdlU1MyOyBwcm94eSBvZiA8U3dpZyBPYmplY3Qgb2YgdHlwZSAnaXRrSW3igKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view(im_view)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise 1: Visualize an image\n",
    "* Read an image with ITK\n",
    "* Apply a filter\n",
    "* Show both original image and filtered images with matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/4_ITK_and_NumPy_answers_Exercise1.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exercise 2: ITK image to NumPy array\n",
    "\n",
    "* Read an image with ITK\n",
    "* Convert image to NumPy array\n",
    "* Modify a pixel in the image\n",
    "* Has the array been modified?\n",
    "* Why did you make that choice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/4_ITK_and_NumPy_answers_Exercise2.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# ITK Image and NumPy array index order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's look at the size of an ITK image and the size of a NumPy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image size: itkSize2 ([570, 326])\n",
      "Array size: (326, 570)\n"
     ]
    }
   ],
   "source": [
    "im=itk.imread(\"data/CBCT-TextureInput.png\", itk.UC)\n",
    "arr = itk.array_view_from_image(im)\n",
    "print(\"Image size: \" + str(itk.size(im)))\n",
    "print(\"Array size: \" + str(arr.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The sizes appear to be inverted!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let's see what happens if we get the value of a point at a certain location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pixel value at image location [300,240]:241\n",
      "Pixel value at array location [300,240]:254\n",
      "Pixel value at array location [240,300]:241\n"
     ]
    }
   ],
   "source": [
    "print(\"Pixel value at image location [300,240]:\" + str(im.GetPixel([300,240])))\n",
    "print(\"Pixel value at array location [300,240]:\" + str(arr[300,240]))\n",
    "print(\"Pixel value at array location [240,300]:\" + str(arr[240,300]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Indeed, it is inverted. Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "ITK and NumPy use the opposite convention for their indices order. It is not important when you use filters and work on the image as a whole. The algorithms are aware of the order if they need to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Could the ITK axes order be kept when converting an image to an array?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Yes. An option is being added to the function `itk.array_from_image()` to keep the axis order.\n",
    "When *keepAxes* is *False* (default), the NumPy array will have C-order indexing.\n",
    "This is the reverse of how indices are specified in ITK, i.e. k,j,i versus i,j,k.\n",
    "However C-order indexing is expected by most algorithms in NumPy / SciPy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# scikit-learn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* scikit-learn is a machine learning package in Python.\n",
    "* It is used to illustrate solving a problem using ITK and NumPy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import sklearn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "First, we load 10 2D-images of circles with different radii and center position to which some noise has been added, and their corresponding ground truth segmentations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Image')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADHCAYAAAAJSqg8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJztnXu0FNW17r8JKIj4AhERjICoSFRQCYKAIA8FFSFGo+hNIHJCkntyYxLHMOpNMpKcxJuc5B5jRjzXYYYPTHy/EfAN+OIhG1FU8AFGxQcighFFA8q6f3T16l91drN7s/fu3rv2/I3BYO7q6lpV1atXV331zbkshCDHcRyn5dOm2jvgOI7jNA4+oDuO42QEH9Adx3Eygg/ojuM4GcEHdMdxnIzgA7rjOE5G8AG9FWNmV5nZz6q9H47jNA6tdkA3s+FmttDM/mFmG83sKTP7SrX3q1zM7Hoz+3U91p9mZk9yWQjhuyGE/2j8vXOyipm9bmZjq70fTu20q/YOVAMz21PSbEnfk3SbpF0ljZD0z2rul+M4TkNorVfoh0pSCOHmEMIXIYRPQwgPhRBWSJKZnW9mq8xsk5k9aGYH5d9oZieZ2cvJlf1/m9ljZvZvyWvTkiv9y83sQzN7zcyOT5avNbP1ZjYV22pvZn8wszfN7L1EAtkteW2Umb1lZhcm73vXzL6VvDZD0nmSLjKzj83svmT5xWa2xsw2m9lKM/tqsvxwSVdJGpqs/2GyPHWVb2bfNrPVyR3LLDM7AK8FM/uumb2aHNuVZmZN8/E4zZ2d6OunmtlyM/soef0XRdv7ppm9YWYfmNnPeCdgZm3Qtz8ws9vMrHOFD7lF0FoH9FckfWFmM81sgpntk3/BzCZJulTSGZK6SnpC0s3Ja/tKukPSJZK6SHpZ0vFF2z5O0ork9Zsk3SLpK5L6Svofkv5sZp2SdX+r3I/LwOT1HpJ+jm3tL2mvZPl0SVea2T4hhKsl3SjpP0MInUIIE5P11yh3p7GXpF9K+puZdQ8hrJL0XUmLkvX3Lj4hZjZa0v+R9HVJ3SW9kew7OS05lqOS9U7+11PrtCLq09c/kfRNSXtLOlXS98xssiSZWX9J/63cRUp3Ffp8nv8labKkkZIOkLRJ0pVNeWAtlhBCq/wn6XBJ10t6S9LnkmZJ6ibpfknTsV4bSVskHaRch1yE10zSWkn/lvw9TdKreP1ISUFSNyz7QLkB3JTr5AfjtaGS/p7EoyR9KqkdXl8vaUgSXy/p13Uc47OSJmHfnix6PW5D0jXK/UDkX+skaZukXsnfQdJwvH6bpIur/Tn6v8r+k/S6pLH16esltvNHSZcn8c8l3YzXOkraKmls8vcqSWPwevekb7ZrjGPK0r/WeoWuEMKqEMK0EEJPSUco98v/R+UG7iuS28gPJW1UbvDtkayzFtsIyv0gkPcQf5qsV7ysk3JX/x0lLUNbDyTL83wQQvgcf29J3lsryW3rs9jeEZL2reNU5DlAuavy/LF9rNwXkldK68rdF6dVUG5fl5kdZ2bzzex9M/uHcneM+b5Z/L3aolzfy3OQpLvRr1dJ+kK5CzAHtNoBnYQQXlLuavUI5TrWd0IIe+PfbiGEhZLeldQz/75EQ+5Z2zbLYINyHf7LaGevEEK5g2SqTGai8/9F0vcldQk5WeUF5X6M/mX9WnhHuS9Ofnu7K3cr/XaZ++M4O+Im5e6CDwwh7KXcM5183yz+Xu2mXN/Ls1bShKLvZIcQgvfNIlrlgG5m/ZKHjT2Tvw+UNEXSYuU62iVm9uXktb3M7KzkrXMkHWlmk82snaR/V07nrjchhO3KDcCXm9l+SVs9zKxcXfo9SX3w9+7KDdrvJ9v6lnI/UFy/p5ntWmJ7N0v6lpkNNLP2ki6TtCSE8HqZ++M4O2IPSRtDCJ+Z2WBJ5+K1OyRNTB6q7irpFyoM9lLuO/mb5KJFZtY1edblFNEqB3RJm5V7oLPEzD5RbiB/QdKFIYS7Jf1O0i1m9lGyfIIkhRA2SDpL0n8qd0vYX1KNdt7u+BNJqyUtTtp6RNJhZb73Gkn9k9vQe0IIKyX9X0mLlBu8j5T0FNafJ+lFSevMbEPxxkIIj0j6maQ7lbtiOljSOTt1VI7zr/xPSb8ys83Kaea35V8IIbyo3IPPW5Trex8r97wo/726Qrmr+4eS9y9W7vvrFGHJQwZnJzCzNspp6OeFEOZXe38cJwskzpgPJR0SQvh7tfenJdFar9B3GjM72cz2TmSJS5W7NVxc5d1ynBaNmU00s47Js5s/SHpeOUeNUw98QK8/Q5Xze2+QNFHS5BDCp9XdJcdp8UxS7sH8O5IOkXROcPmg3rjk4jiOkxEadIVuZuMtlwa/2swubqydcpxq433baYns9BW6mbVVLoV+nHIPBpdKmpK4LRynxeJ922mpNKTa4mBJq0MIr0mSmd2inA5WstObWavWd4499thq70JFWLZsWdXaDiE0RsGwevXt9u3bh913312StMcee8Tln35aeLSyefPmGO+zTywdpH/+s+B4/eKLL2Lctm3bGBdfdPH969YVknc7duwY4zZtCjffn39eSDbmdllbrdRy7hO3s+eee9a6TqljKNWWJO2yyy4x/vDDD2PcoUOHGPMcfPLJJzHu1KmQh8ft8rzuumsh9YLbZFul1uFnyO3zXBT/zXPP/WYfyPeXYtg2ef311zeEELrW+iJoyIDeQ0jXVe5K5l+8oZarDDijAe1khpqammrvQkXIQBHGOvs2+3XHjh01ZswYSdLYsYVS4S+88EKMH3nkkRifddZZMX799ddjvGnTphhz0N62bVtq5772ta/F+Pe//32MBwwYEOPddtstxu+//36Mu3QpJGByIOUPUbt2hWGBg9CGDYX0BR7nxx9/XOcxcJt77bVX6ni6d+8e43vuuSfGhx56aIy3bt0a4yVLlsR4xIgRtR7D3/9ecDv27FlI5u7fv3+tbfXu3TvGhxxySIxffPHFGHMQ5o+BJH3wQaFSQfv27WPMgX7evHkxHjJkiGrjsMNqT0OZNm3aG7W+UEST10MPucqAV0ut5wq9tT9o3tHxZ2Cwl5Tu17169QqnnHKKpPQXmwPS2rWF34czzzwzxosWLYrxli1bYswvPwczSVq4cGGMv/KVwpwsvFr/0pe+FOOjjjoqxhxw9923UObnoYceijF/DA44IFZQ1t57F4p07rfffrXu9zHHHBPjl19+OcYcJF977bXU8axatSrG48aNizH70eLFBWfwaaedFmPe9T722GMx5g/OnDlzat1v/gDMnTu31mOYMmVKjHneiwf0559/Psb8oeSV/8iRI2PMHwr+8H322Wcx/stf/qL60pCHom9LOhB/95TX/XCygfdtp0XSkAF9qaRDzKx3Un/hHOXScx2npeN922mR7LTkEkL43My+L+lBSW0lXZvUZGiVtHaZpVx4npqr/FLfvr1169aohVNzplTCh3S33357jPkA+eCDD44xb+kpH0jp23LCB2rUhCknHH744TGmzsz2+GDziCMK9d0odSxYsCDGnTsXJg96++3CjQz3k88BunVLV71ds2ZNjEs9zKR81a9fvxjPn1+ouLF69epa95X7sWLFihhTBjr11FNjfN9998WYElqfPoVaeP/4xz9SxzBjRuEx4ezZs2Pco0eh+jQlHko/PE4+7zjxxBNjXK7RoEEaeghhrqS5da7oOC0M79tOS8RT/x3HcTJCRVP/s+BycWml6WmIFNNIPvR6cfDBB4fLLrtMUtrNwlv9lSsLFvZXX301xrQa0v5HLzO90JI0fvz4GH/00Ucxpie71LbocqF1juvQz05JgzLGe+8VJibiNmltpPxEuaJYcuExvPFGwZ03aVKh5DnboAWU8su7774bY/YhykmUOtgunSb771+Y4uCBBx6I8dChQ2Nc/JlQKqK0wvNEG+vGjRtrPYann346xrRz/u1vf1sWQhikOvArdMdxnIzgA7rjOE5GaPLEoizgMktlaQlOGLJx40bdeuutktK35XRCMCln1KhRMaYzg8dKCWTChAmp9piYw0xTyjfPPvtsjJkcRCmmVEkASiWUkOg0oRuDMgmTjNguJRNuX0q7PA46KE5rm0oIYrIOHSnXXnttjM89tzCr3c033xxjJlkNHjy41mOgXMVkr4kTJ8b4ueeeizGzT6V0hifbu+GGG2LMz5HJZRdeeGGM8wlqUvqzLRe/Qnccx8kIPqA7juNkBHe5lMBlluZHOfJLNVwuffr0Cb/61a8kSXfddVdcztt13t5v3749xpRPKFcw6aX49v7ee++NMWu5UOKhrEEHBm/1KTO88sorMabLhW4MOkQoIbGGzLBhw2LM4ldMsOnbt2/qeLivlIpOOumkGDOxhsWvWHmR+0QnDeum8JjpVKGjhHVZuG90BRUX2qN81atXrxhTKmJ7lGjozmHS2fr162O8ePFid7k4juO0JnxAdxzHyQjucgEuszRvmqv7Zdu2bVFSYRILb93pfHjnnXdi/Oabb8b46KOPjjFLzz711FOp9o47rlCanbfrlGwo04wePTrGlAboeDnwwEJxySeeeKLW7bAtShd06tClQomGJWkpQ0jphJtBgwqqApN6mKxD+YZuILptKBWVmjTjy1/+cq1tUTYqVd+9uJbL8ccfH2PKJqzXztLAXIeuHdba4WfC8sE7wq/QHcdxMoIP6I7jOBmh1btcXGZp+VB+qYbLpWfPnuEHP/iBpHT52FI1RFiLhLf3dJc8/PDDMeatupSuBcM6IJwibfjw4TGeNatQyp3lcFkml+/lfKEsH1s8F2geluGlo4YyBN08bFdKT8lHOYq1Vnheu3YtTK25dOnSGNNdQnmEkg6To/iZ0HnE+jUsaUwnDJ02Uvoc02FD2YTv56xQlIp4nJyqb9q0ae5ycRzHaU34gO44jpMR3OXiOA3k888/j9IEa5/QtULHy/Tp02PMyZlZbpd1P4prejDhhtIMHSmUR+jsoHuGNVSOPPLIGHOyZSbTsAQuZYxp06bVejyUWbidYofIzJkzY0w5hcdAqYiyTv/+/WNMiefGG2+MMZO0mHxEhwylHib9MHGJyVft2qWHTkpWfA+dREyO4ixPPPcst8uEqHLxK3THcZyM4AO64zhORvAB3XEcJyO0StuiWxWzyaBBg1RTU1Nx22L37t1DXhen7kstlVmJ1L2pn7KAF7Mhi2d879KlS4xpf/vhD38Y4z/96U8x5uzxtExySjlqurQhsm1aLKn3cx3q7Dw2Tk3HYl7FUF+nVk5LJ7fF5wnPPPNMrftBeyLbpqWQ557TwHE6OtZ3L9bQqc3z2QGzfPmZMnOYGj+fl3Cbo0ePbhzboplda2brzewFLOtsZg+b2avJ//vsaBuO0xzxvu1kjXIkl+sljS9adrGkR0MIh0h6NPnbcVoa18v7tpMh6rQthhAeN7NeRYsnSRqVxDMlLZD0k0bcr0bHZRanmMbq223atImSwNatW+NyZvrRIsca5ry9p5TAW3oWpiqGdjnWH6dlcuXKlTFmTe65c+fGmFmWjCkzsDgXv088Htr07rvvvhhTVqAcIqXPGW2LhBmUtGHS2kcZY8GCBTHmOaL9kRZBSjEskMWsVko6lGKkdGE0ylE8HspGLG7Gc8bPZ9u2baovO/tQtFsIIf/prpPUbUcrO04Lwvu202JpsMsl5H6qS17+mtkMM6sxs5pS6zhOc2RHfZv9mg8/Haea7Gym6Htm1j2E8K6ZdZe0vtSKIYSrJV0tNR+Xi+PsgLL6Nvv1AQccEPKywUsvvRTXoZTAglqUNOhYYUEpZpZyKjJJGjBgQIxPP/30GDPDkZID5Q7WJee+MuOStbdLOWTYLtvi8YwfX3g8cdNNN8WY9dwl6e67745xsZSRhw4RTgtHuYKyDF1FlFBYMI0SEt9LpxIlEO5b8X5SyuGxcv+4H5RfmIHKz4dZt+Wys1fosyRNTeKpku7dwbqO05Lwvu20WMqxLd4saZGkw8zsLTObLum3ksaZ2auSxiZ/O06Lwvu2kzUynVjkzpbWRbUSi/bbb79w9tlnS0q7U9j/WNubUgzrmTOphHW7mWBTvF1OKTd27NgYc5q3UrPY9+7dO8Z0whDuH+UUFiGjG6PULPecmo7yk5SWkFhn/M4774wxE4jonuH0enQSMbmH9fJZW57HzJhFvlatWhVjJvrw2KS0a4XOGCZ+URJjkTRKLnQA0Tlz/vnnez10x3Gc1oQP6I7jOBnB66E7TgPp0KFDdJI8+uijcTlvy88444wYM6Hl/vvvjzHlGrogGEvp23K+tmbNmhjTMUPXxsiRI2NcXJc8z8knnxxjShGUijht2qhRo2L83HPPxZj11plkRelGkjZu3FjrazxOShqUWZhkRGmJUhS3s3DhwhizfjxlH067xxovXJ9SlJSWy/g50NEze/bsGFNCevPNN2NM+eqJJ55QffErdMdxnIzgA7rjOE5GcMnFcRoIa7lwlnjW9FixYkWMeUs+ZsyYGDPphRIFa7xI6aQWzjBPaeWtt96qtT06W+iioGxCSeiEE06I8bPPPhtjJiItX748xnRyUK6gdEMniyQ9/vjjtW6XsgbboHOE54wuFJbGLSUt0ZHDmjqUd7gOt0mZREp/JswcZv0X7jddLpSN+HlS+rruuutqPYZi/ArdcRwnI/iA7jiOkxEyJbl4IpFTDb744ot4W0+nymOPPRbjSZMmxZguDc4ORJfG9773vRgPHz481R7rjlDWYNucxX7w4MExpiuGbg5KNKw/QmnhqKOOijHdGNwflqrlvlFKqKlJ1+lj3RV+hymnUJp68sknY8xZkbic8ghLADNRiAk9XJ+yGeUqunaYuCSlE8foYLniiitizPPHksO33357jDnDEeWncvErdMdxnIzgA7rjOE5GyJTk4jjV4LPPPov1Sfbcc8+4/Kc//WmM6RChlECZgM4WJqowoUVKSyKUb+ikYRlWzuTD5BbKI6xfwuQo1kdhzJmCKDewlC6lDsosxRMsT5kyJcYsGUtZg/VsWCuF+8EEIu4rywTThUO5i6VteQxMxKIriBKNJO27774xZlleSm10IV111VUxpsTD5K1NmzapvvgVuuM4TkbwAd1xHCcjuOTiOI0IXSGUVoplhjxDhgyJMaWBF154IcbFU9x95zvfiTETdpgEw7YpA7EsL+udUK4hdKdQ6qArhnVjuK90clA2WrRoUaqNO+64I8aUmigbUdLgjECcDJvHSZdPqZLGnAmKEg2ThHge6dopruVCeYjyGvebSV1Tp06NMV1CdA/xXJSLX6E7juNkBB/QHcdxMoJLLo7TiHBCZkoUlAMIpYju3bvHmG4PlnaV0rf0nIGI9UhYB4Qz/wwcOLDWtrmvnBiajo/OnTvHuGvXrjFmHRhKFKxBw+0UO0QomwwaVJiUh8dDlwxdO5R79thjjxiXmvSasglnUaKjhPILZSye6+KSxtu3b4/xfvvtF2M6huh+YZlgxnQMFUtt5eBX6I7jOBnBB3THcZyM4JKL4zSQdu3aRTmCjgW6NOiCYPII3SyUElgKdujQoan2uC2u941vfCPGLNvKeidM9uFMQXR8cJuUQ1hnhMdA5wiTo+bMmRNjJkBR9pDS5XpZX+a0006L8axZs2rdPyZZUfags2XEiBEx5oxFLKtLNwrL+dK1dMQRR8R46dKlqWOgZMWEpfPPPz/GnMGKE1fPmzcvxpR46BIqF79CdxzHyQh1DuhmdqCZzTezlWb2opldkCzvbGYPm9mryf/71LUtx2lOeN92sobVVXLWzLpL6h5CeMbM9pC0TNJkSdMkbQwh/NbMLpa0TwjhJ3Vsq0nr23r53NbNoEGDVFNTY3WvmaOx+vaXvvSlcNFFF0mSrrzyyrj829/+dozpoqDThDIGb+8nT54cY0oJUvq2nM4OSjGUaW688cYYl3Kq0NnRo0cP1QbrqbBduje4HSbJ0PHC2iXFULpg4hOlGCbo8Fzk6+lI6YQrunaYNMTxgtthyVxKOiy9O2zYsNR+U7Kie4ZuGDpYeC65HzyvPBeTJ09eFkIoWIBKUOcVegjh3RDCM0m8WdIqST0kTZI0M1ltpnJfBMdpMXjfdrJGvR6KmlkvSUdLWiKpWwgh/9O8TlK3Eu+ZIWnGzu+i4zQ99e3b7Nf0QjtONSl7QDezTpLulPTDEMJHfEobQgil5JQQwtWSrk620eiaiMssTkPZmb7Nft21a9eQT/aZMGFCXIeSAeO77747xn379o0xpRVOHMy6JFI64YSukl69esWYTouxY8fGmNLCn//85xizhG2p8rSlJkmm64S1X1avXh1jlq2l7COlv8M8Tyw3y+3yB3T+/Pkx5sxMPAbu92233RZjziDEc0zZg7MpcXLr4omuuR6dTpy1auLEiTHm7FQsXczJsOmKKZeyXC5mtotyHf7GEMJdyeL3Eg0yr0WuL/V+x2mueN92skQ5LheTdI2kVSGE/8JLsyTlS4ZNlXRv4++e4zQd3redrFGO5DJM0jckPW9m+SISl0r6raTbzGy6pDckfb1pdnHHFN0eV2MXnJZLo/XtfD9kCVwm39D5wFtySgmsrcJ1eDtf3AYdM3RkcEJjJi9R1vn+978fY8o4lFnoFmHNEe43j5NSRKdOnWJMx8sjjzySOp433ngjxl/96ldVG3TklJoA+8EHH4wxxwU6TcaNGxfjJUuWxJhyFbdPdw0lELpaJOnYY4+NMWUjyil0yVB2KjU7Fc99udQ5oIcQnpRUygo2psRyx2n2eN92soZnijqO42QEr+XiOA1k8+bNsf4H62+sXbs2xrytZj0V3obztp+yx2677ZZqj7McUZqh9EHJhftx+umnx5hyAkvyUiZgEhDbYpngl19+Ocbr1q2LMcvfsp5M8aTXbIP1VSiblIqfe+65GFN+euKJJ2JMCYUTY9MhtGzZshhT4qJcxWQtJkpJabcNXTJMKKNkRQmK2507d26M6UgqF79CdxzHyQg+oDuO42QEl1wcp4F06dJF5513nqS0jMGysKz9wRl0WJ6Wt/es+8HEHUl65ZVXYsz6JXTScDnLx1LKobTA93L7lEqYrMR9LTU584wZhQRxOnCKJz+mw4by0MiRI2PMc8nzwbozdP9Q+uKE1nwvE5c4wxOPh/IQJRCWOpbS54xuuy1btsSY54aleLl/nLWKyVHl4lfojuM4GcEHdMdxnIzgkovjNJCtW7dGqYCTJLNOCyctZm0R1jvhrDl9+vSJMeugSOkEFybNUB7hOnSb0AnCmXJYVvbwww+PMWUWunOKXR556ITheymZUOqQ0lIEXTisF0MnDNvmZNqUr1gLhxLXpEmTYszJnClxcV/pUqF7hedUSteF4edI+YbnnpIYpRWuw4mxy8Wv0B3HcTKCD+iO4zgZwQd0x3GcjOAauuM0kDZt2sRsThanYjYl9VAWyGLWJ3VVTk132mmnpdqjtkxdm5ZJavPMGmUmZymNn9ZBFsWizsxZ788888wY04LYrl1heKFuTuukJI0fPz7GrHXOY2AGKbMyTz311BjTDvnZZ5/FmM8ZWJBr0KDCjG6l9nXMmEJJH061V1wIkPvNz5Hnnpm8LFzG5yXMcN2ZiVP8Ct1xHCcj+IDuOI6TEVxycZwGYmZRcuFs7gsXLowxizExy5AZl7TCUeoozqxkMSzW5T7hhBNiTGmFNklmO3K6M9YJf+aZZ2Jcqi47pQ7KLF//eqF0PK2TlDGKbZh8P2ua0/63YsWKGDN7kwXNWGCLchflJMpStD/yHLOeOYt5MaOTslnx/lEqoZxCqYzng9mrjLkf5eJX6I7jOBnBB3THcZyMkCnJhXWSJZ+SzqkMIYQotQwYMCAuZ9Yoa28PHz48xr/5zW9izNt7Fm9atGhRqj3KCXRwULpgpiSdFpR+zj///BhTuqGcQGmAmaicNo7vff/992uNOWUd64JLaUcPJQq+h3IHpRgeG7//lHvoKKGbpdQ2H3vssRgzg5ZyCF03ktSvX78YU3biNHKUhDhlHafko0w3ffp01Re/Qnccx8kIPqA7juNkhExJLo5TDbZt2xaTTugooVOF8h+LYo0dOzbGlEyYnNK7d++SbTPZhVPH0c1BNwulHEoU9913X4yZcEN5hNIIC1VRSjj00ENjTEcN5SeeF0kaNWpUjCkV8RgmTpwYYzpvOHUek6x4znj8LPJF9wvlIX5WdOSwUBcTl6R0ca45c+bEmMlbPG6+n5IV29sZybjOK3Qz62BmT5vZc2b2opn9Mlne28yWmNlqM7vVzHata1uO05zwvu1kjXIkl39KGh1CGCBpoKTxZjZE0u8kXR5C6Ctpk6T6K/iOU128bzuZok7JJeSu+/PzR+2S/AuSRks6N1k+U9IvJP2/xt9Fx2kaGqtvt23bNta9ZnIQb59ZH4XTjDEBiC4X3qq3aZO+7uJ7eOtO1wZlBibfcCq42bNnx5iSAWUc1kZnQtStt94aY0orrO1NhwzlBtYbl9IOEx4rHSaUVniemMhFWGuGcg/dKXQIUd6gnMRjY313umukdIIYXUKUX5iYxTrudLxw6rxiaaocynooamZtzexZSeslPSxpjaQPQwj5VKa3JPUo8d4ZZlZjZvWv1u44TczO9m3262I91XGqRVkDegjhixDCQEk9JQ2W1K+Ot/C9V4cQBoUQBtW9tuNUlp3t2+zXfJDnONWkXi6XEMKHZjZf0lBJe5tZu+RKpqekt3f87srDp/ieZOTsiIb27by8QKcF5QDehrNuCH8M6N6g84GSiZSuo8LX6Iw56KCDYkzJgUk2TJShjMNtUkqgc4TOFO4338vkHtY3eeWVV1LHQ5mFiTiUgXjOKJuce+65MaYsQ3mIziM6Z+jIueWWW2LcuXPnGHOaOsJzJ6UdSjwHw4YNizElOLpqCMvnFrdRDuW4XLqa2d5JvJukcZJWSZovKV8Ieaqke+vduuNUEe/bTtYo5wq9u6SZZtZWuR+A20IIs81spaRbzOzXkpZLuqYJ99NxmgLv206msEpKEWZWNd3DJZfsM2jQINXU1FjdazYu3bp1C1OmTJGUdkWw3glLoVJioOuEyzlbPGc1kqTHH388xieeeGKMKfcwCYZyB+ujcOYgygxch/KDDHXQAAAP2ElEQVQIXSusCUNpiTIB3TLcN26n+DXKNDxn+fLExcvZNuvAUH7hd59OIDpKOJsQpSs6YShd8fxKaecSXTKUspjUxNmYfvSjH8X4ySefjDGlon79+i0r5zmkp/47juNkBB/QHcdxMkKrqeXijhenqWjXrl2ULHgrvXHjxhjTscFSulu2bInxa6+9FmMmFnGGIimdcMOZctiv6WZhshNrotC1QhmHUgTbopTACaMpV1D2YKlayjjFiUV0ybCULMsEU3JhPRrWkaEMQrmH8hVr0/Acs11+VnSmPPDAAzHmZyulpSY6lzjB90svvRRjHts999wTY56zzZs3q774FbrjOE5G8AHdcRwnI7QayYW4/JItimeqqgb5W2UmqxD2s+uuuy7GdKlwJiNO1EwHipROrOFMOc8//3yMf/zjH8eYEsLixYtjfMwxx8SY0sBtt90WY0oATIbhzEyUGCit0GnC2Y4o6UhpaYFOEsod7du3jzGTjMaMGRPjmppCdZGTTz45xkwmorxDx8tZZ50V4zVr1sT42muvjTFrzhTD80R3CmUWOoN4bDx/lN1YZrlc/ArdcRwnI/iA7jiOkxFaTWJRKVxyafkUSWgV1186d+4c8jMPUXKgBDBixIgY81afiTRMTqHE0KNHutgjJQdui+4PyhpMdqL88uCDD8aYEgBdKHSL0EXCGif8DlEe4na4TnFtGko/119/fYxZopalZOfOnRtjSiulasKwbC1rvLDULytm0sHDc8/zXlxhky4cumeYFEY3C5O9KPFwuw8//HCMZ8+e7YlFjuM4rQkf0B3HcTJCq3S5EHe8tEyag7MlT8eOHePEykz0OeWUU2JMCYXSCm+x6aKgk6U4iYVSBrfF2XHo7Dj77LNjfMcdd8SYE1QvWbIkxizhevrpp8eYEsiGDRtiTJcKJQ26XFgThfVrpLQkxDov/Iwpp4wePTrGLFs7f/78GF922WUxZu2XhQsXxpgSFeUUOpVYK4YOHLqTpPRsRlyPyUis68LzxM+Qrpgmm7HIcRzHaf74gO44jpMRWr3kQlx+ad40J5mFmFmUDZhwQwcKHRVMJOGEx3RXUBqhC0RKz0bEEr0jR46MMW/vr7rqqhhTiqC7gi4N1nuhtNK/f/8YU0Jhsg6lHsokfO8ZZ5yROp4//OEPtW6XDiDWYGH9m+OOOy7GPXv2jPGll15a6/Fwm0z0YYIWJS3O2MRzRIlGSss6rOfDz4rJW3QY8XOgq4j1eMrFr9Adx3Eygg/ojuM4GaHVJxaVg8sv1aO+Mks1Eot23333kJcUWMeDCS28laabhbfbdMJQriieVJllbyll0GkxcODAGC9atCjGdKpQHmIi0vLly2Ocd+9I0qOPPhpjJkpRili6dGmMOdE1ZY+77747dTxMDlq7dm2MKa2wFgxlHcomlDdYW4X1azg7EtuiU4fuJJYVZpIRj1+SbrrpphhPnjw5xjwf/Kzp5uHny/LDLJt8zTXXeGKR4zhOa8IHdMdxnIzgLpcycPdLZWmubpZSdOjQIc5Aw5lrmBjC23s6IjgL0DvvvBPj448/PsbFZVt5u84Zew4//PAY03lCKYJJL5xpiC4cJjVR+mECzKpVq2JM+YFuFC6fNWtWjCdNmpQ6HkozTDqiA4jSCt0sLPVLBwsdRpSE2BblFy7n5Nl8Lz9blkCW0k6k+++/P8Ys7/vmm2/GmKWO2d727dtrjcul7Ct0M2trZsvNbHbyd28zW2Jmq83sVjPbta5tOE5zw/u1kyXqI7lcIGkV/v6dpMtDCH0lbZI0vTF3zHEqhPdrJzOU5XIxs56SZkr6jaQfS5oo6X1J+4cQPjezoZJ+EUI4eQebabEul3JwKab+NIW0Uh+XS2P164MOOijkE1k42fA555wTYzpEWGKWTg7KIXRy7GjGIsoPLGnLGY9YGpfbojOD26HMMGPGjBhT+qEzg9JAqaQnyjuUcaS0lMHyw0y+4XuOPvroGFN+osuHrpi77rorxpSNhg0bFmPKXfxM6DShBMTZpaR0zRfKYNwWnTc833TScB0mFl1++eWN6nL5o6SLJOU/uS6SPgwh5MXAtyT1qO2NjtOM8X7tZIo6B3QzO03S+hDCsrrWLfH+GWZWY2Y1da/tOJWhMft1cfVAx6kW5bhchkk63cxOkdRB0p6SrpC0t5m1S65mekp6u7Y3hxCulnS1lG3JxZ0w5dGMHCyN1q/79OkT8rfjlA/oXuCgT0cJJRBKJlxO6aL4Nc58w7K3THCiI2X8+PExPvXUU2vdJzpkKCHRFcIkG5abZflgJjExGYqTTUvSX//611r3j+eSPPTQQzFmGdtjjz02xjNnzowxa6gwcYdSDGdRorzDMrd77rlnjHmOpPTnyyQyHvfTTz8dY0pieYeUlJZZWA64XOq8Qg8hXBJC6BlC6CXpHEnzQgjnSZov6cxktamS7q13645TJbxfO1mkIYlFP5H0YzNbrZz2eE3j7JLjVBXv106LpV6JRSGEBZIWJPFrkgY3/i45TmXxfu1kBS/OVUFai7ZeTZ28GsW5OnXqFPI2OWqmzAKl5kqrHYtfsY4216EWLUn7779/jGkfZD116tqbNm2KMW1x1H1p4WOWJdeh1Y5ZmbQRcn1aGDn1G/eteJ/43IF6MvV02iqpj99+++0xpl2wVAEzTv9Hy2ipZwJ8lsFnCFJ6Gj3aSvld4LnhOaZVk+eCmb8XXHCBF+dyHMdpTfiA7jiOkxG8OFcFKUeKaO6yTDOyHTYbQgix4BZv42kF5G0/Z5u/8cYbY0zrHItRMctSSssdjCnNUGbhZ8aY0g2terRCMmOVdcUpgdCGyCJXU6dOjXHfvn1r3WdJ2rhxY4wp61BCYr1yTgtHKyG/Ozx+Zm7SzklZhoWzaPPkZ8W2uM9S+rxSNuF5YkYtM1ZfeumlGHMKPmaflotfoTuO42QEH9Adx3EygksuzYzGlDR4C+pSSdOx6667xlt53mJ/8sknMWYW6IYNG2JMJwhdMUOHDo0xMwyl9K07b8sp2dA9w3XowKDM0KFDhxjfcccdMT7uuONq3Q9mdFK6oGRAqYPFr3j8UtrRQ7mHNcYpY7BIFveVDhZKJcysZYYm95sSyLp162JcypnC45TSnzX3g1LbXnvtFePXX389xnTVcB32mXLxK3THcZyM4AO64zhORnDJJcO4zFIZunTpom9+85uS0u6UefPmxZh1tZmEwlt6yg28vackIaWdIMuXL48xXRhMiGnfvn2M6eCgI4cuF8pATOLhOtwnSgOsQ05XDNstPh5O0zZnzpwYM7GIkgulKTpsWGzsvffeizElJzpk+DnQkcTjoVREVxCXF7+HThruK4unUQ6lg4cuoWJ3Uzn4FbrjOE5G8AHdcRwnI7jk4jgNZPPmzVFeoXOCSTms200XBd0lhC6K4tt7/k0nyMqVK2NM+YauEMovlFDoeGFtEU6dN2HChBiXqv/NhB46P7p06RJjyjLF7dEltGbNmhhToqBr5eCDD45xTU1hDp3jjz8+xnTbsKY7zwXrtFACofQ1bty4GFN+kdKfKc89t8XjZj8ZPLhQC+7555+PMWv7lItfoTuO42QEH9Adx3EygksujtMI5OULOkcoh3A6sZNOOinGlGXoKKEbY+TIkam2brjhhhjTCcJkIsodlHV4S0/nCBOI6P6g5EC5gg6WESNGxJgyyYABA2JM1wmloeL3cDo2yjpMiKJEwbowdLPQLUIHChOLeDxDhgyJMV0+CxYsiDHPRbEDhc4dfo5099Dxwn7Cz5ryEKWycvErdMdxnIzgA7rjOE5G8BmLnExRjRmL9t9//5BPLKLL45577okxZYnDDjssxnRXUG5gwszYsWNT7dGdwpnrCZOPeEtPhwjrl7D+CPejlIzBJCbKB1zOcr6UJFivRErPFkRZglB+oVREKYbv5Tmi5ETXCWUpOl4o0VB+ohOI50tK12PJl1KWpIceeijGrC/D5COe4379+sWYtXOuueYan7HIcRynNeEDuuM4TkZwl4vjNJDt27fHW3mWS6W8QWmATgtOsEznBBNsKMtIaccMnTSc4YeJL1zOmiVHHXVUjOkWYWlXShSlSts+9dRTte4r5QPKFcVlYXme6PI44YQTYsxEK26L54nyMSfr5oxAlGL4mTDxicfAUsVr166NMSf9ltISCie6ZjIR3S+UhygJ3XzzzTEePny46otfoTuO42QEH9Adx3EyQqVdLu9L+kTShrrWbSL2rVLb1Wq3mm1Xo92DQghdK9xmvl+/odZ1rqvddms75rL6dkUHdEkys5py7DdZatuPuXXQGs+1H3PzwiUXx3GcjOADuuM4TkaoxoB+dRXarHbbfsytg9Z4rv2YmxEV19Adx3GcpsElF8dxnIxQ0QHdzMab2ctmttrMLm7itq41s/Vm9gKWdTazh83s1eT/fXa0jZ1s90Azm29mK83sRTO7oBJtm1kHM3vazJ5L2v1lsry3mS1JzvmtZrZrXdvayfbbmtlyM5tdyXabA96vm7xt79tlUrEB3czaSrpS0gRJ/SVNMbP+O35Xg7he0viiZRdLejSEcIikR5O/G5vPJV0YQugvaYikf0+Os6nb/qek0SGEAZIGShpvZkMk/U7S5SGEvpI2SZreyO3muUDSKvxdqXarivfrirTtfbtcQggV+SdpqKQH8fclki5p4jZ7SXoBf78sqXsSd5f0cgWO+15J4yrZtqSOkp6RdJxyCRDtavsMGrG9nsp9mUdLmi3JKtFuc/jn/bqybXvf3vG/SkouPSStxd9vJcsqSbcQQn6K8nWSuu1o5YZiZr0kHS1pSSXaTm4Nn5W0XtLDktZI+jCEkC/Q3FTn/I+SLpKUr2TUpULtNge8X1egbe/b5dFqH4qG3M9rk1l8zKyTpDsl/TCE8BFfa6q2QwhfhBAGKndVMVhSvzre0mDM7DRJ60MIy+pc2Wlystivk2173y6DSpbPfVvSgfi7Z7KskrxnZt1DCO+aWXflfu0bHTPbRblOf2MIIT+lTEXalqQQwodmNl+528G9zaxdckXRFOd8mKTTzewUSR0k7Snpigq021zwfl2hfi15366LSl6hL5V0SPKEeFdJ50iaVcH2lbQ3NYmnKqcDNipmZpKukbQqhPBflWrbzLqa2d5JvJty+uYqSfMlndlU7YYQLgkh9Awh9FLuM50XQjivqdttRni/bvq2vW+XSyUFe0mnSHpFOf3rfzdxWzdLelfSNuV0runK6V+PSnpV0iOSOjdBu8OVu+1cIenZ5N8pTd22pKMkLU/afUHSz5PlfSQ9LWm1pNsltW/Ccz5K0uxKt1vtf96vm7xt79tl/vNMUcdxnIzQah+KOo7jZA0f0B3HcTKCD+iO4zgZwQd0x3GcjOADuuM4TkbwAd1xHCcj+IDuOI6TEXxAdxzHyQj/H9u4qbMWS72VAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "l_label=[]\n",
    "l_im=[]\n",
    "for i in range(0,10):\n",
    "    imageName='data/sklearn/im%d.nrrd'%i\n",
    "    im=itk.imread(imageName,itk.F)\n",
    "    arr = itk.array_from_image(im)\n",
    "    l_im.append(arr)\n",
    "    labelName='data/sklearn/im%d_label.nrrd'%i\n",
    "    im=itk.imread(labelName,itk.UC)\n",
    "    arr = itk.array_from_image(im)\n",
    "    l_label.append(arr)\n",
    "size=itk.size(im)\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(l_label[0])\n",
    "plt.title(\"Segmentation\")\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(l_im[0])\n",
    "plt.title(\"Image\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The goal is to find the segmentation based on the input image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We create arrays of data:\n",
    "* X - the input samples\n",
    "* Y - the target values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "X0=l_im[0].flatten()\n",
    "X=X0\n",
    "Y=l_label[0].flatten()\n",
    "for i in range(1,10):\n",
    "    X=np.concatenate((X, l_im[i].flatten()))\n",
    "    Y=np.concatenate((Y, l_label[i].flatten()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* We use a supervised learning methods based on applying Bayes’ theorem.\n",
    "* The only information provided to the algorithm is the image intensity value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7efe141a6bd0>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD8CAYAAABXXhlaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAEBtJREFUeJzt3V+oHOd5x/Hvr7JVB0KQnR6EkEztYtGgi8bmHIyNe1GUmqpOiHVhik0ouhDoJgWHBlK5hdJAofFNnFz0RsQmugiJUycgYQJFVRRCochexU5iWyRSjENkZOuYWCS5Savk6cWOxNFqdWbO7Dsz75z394FFO3N2Z54zZx69+777/lFEYGZl+YOhAzCz/jnxzQrkxDcrkBPfrEBOfLMCOfHNCuTENyuQE9+sQAslvqR9kn4i6bykw6mCMrNuqW3PPUlbgJ8CDwMXgJeBJyLijXXec93JlpeXb3jNmTNn1j1vm/c00dVx2557o3HUHaNLs/E1iaXNezabNvdX3XV66623eO+991R3nEUS/0HgXyLir6rtpwAi4t/Wec91J5t3bmn9mNu8p4mujtv23BuNY8iu17PxNYmlzXs2mzb3V911WllZYTKZ1B54kY/6O4FfrNm+UO0zs8zd0vUJJB0CDnV9HjNrbpHEfxu4c832rmrfdSLiCHAEYGVlJSaTyboHnf0ok+pjbt1xUn2sr4t/yCpFV9p8TB/qo/28a50illQf29scp817Fvmo/zKwW9LdkrYCjwPHFziemfWkdYkfEVck/R3wn8AW4LmIeD1ZZGbWmYXq+BHxHeA7iWIxs564555ZgTpv1V9PXw0Z87RpRKw7RpPjdPHdbZ+6ahxrcq4mf7O6WLqKNUVsfXKJb1YgJ75ZgZz4ZgUatI6fouNNW131FW/TgWdM+ox/qPp6ir9ZV/d2qnxwiW9WICe+WYGc+GYFaj0ev9XJEozHTyXFYKC+vpNP9Z1wzu0NTX7HnONPJcV8FBHR6Xh8MxspJ75ZgZz4ZgVy4psVqNfGvdkZeLpqtJqnzYSWQzXm2bBynjy04b3uxj0zu5ET36xATnyzAg06SGeeNpMvtDnurFSdibqYsXWBRU+SHKc0Ker0nojDzLLjxDcrkBPfrEBOfLMCZde4N6uv0XpdLdU15PJSJY5u68tQ17bJarlNuMQ3K5AT36xATnyzAg06A888Qy0z3ea4XXX6cV3cNsIz8JhZI058swI58c0KlN1KOikGP7R5TapZdlPU112n35z6WL3J3+Ob2U058c0KVJv4kp6TdEnSa2v23SHphKRz1b+3dxummaXUpMT/KrBvZt9h4GRE7AZOVttmNhK1iR8R3wd+ObP7UeBo9fwosL/JyZaXl4mIa4+bnG/d18z+PCKQdN1jniavSRFLndk4+hqEtBmM/bpt9F7pUts6/vaIuFg9fwfYnigeM+vBwo17Mf2v66b/fUk6JGkiabK6urro6cwsgbaJ/66kHQDVv5du9sKIOBIRKxGxsrS01PJ0ZpZS28Q/Dhyonh8AjrU5yLz6bt1jbOa1A+RU15uV8/XO+bo1kdO1bfJ13teB/wH+VNIFSQeBLwAPSzoH/GW1bWYjUdtlNyKeuMmPPpY4FjPriXvumRWo10E6Z86cWbhuk2qQjs3X10pAJU42ktPv6BLfrEBOfLMCOfHNCuTENytQr417y8vLTCaTa9tdrV7T5jU5NbyMnWcZ6lebBnOX+GYFcuKbFciJb1agTbFabpMOPH3Z7B1T3DkqP55l18waceKbFciJb1ag7Ov4s9quUNtXPXSz13c3++/XpZzaf1zimxXIiW9WICe+WYGc+GYFym4GnhTLWZvlqM2As64aAF3imxXIiW9WICe+WYGy68DTVR1+s0/EMbbBMzl1ZslJX9fBJb5ZgZz4ZgVy4psVKLs6/qw2dZ6uBu3kvEJMqvP0FX+Jdfqc2jVc4psVyIlvViAnvlmBnPhmBcqucS+nGXNnlbBCzNjjb6OLRrdUHaq6ahB0iW9WICe+WYFqE1/SnZJOSXpD0uuSnqz23yHphKRz1b+3dx+umaXQpMS/Anw2IvYADwCflrQHOAycjIjdwMlqe13Ly8tExLVHE2tff7NHm+N0RdJ1D0unq2tbd2/MnrdpB7EU91xX921t4kfExYj4QfX818BZYCfwKHC0etlRYH+yqMysUxuq40u6C7gPOA1sj4iL1Y/eAbYnjczMOtM48SV9EPgW8JmI+NXan8X0M8jczyGSDkmaSJqsrq4uFKyZpdEo8SXdyjTpvxYR3652vytpR/XzHcClee+NiCMRsRIRK0tLSyliNrMFNWnVF/AscDYivrjmR8eBA9XzA8CxumNdnWV3I40z8xpWZh99Ndw10UVDkU01+Tt3cW1TNdTl9Hdv0nPvIeBvgR9LerXa94/AF4BvSjoI/Bz4m25CNLPUahM/Iv4buNl/Tx9LG46Z9cE998wKlN0gnbpBCUPXjRaVakahodoy+oylzQCVodt41pNTbC7xzQrkxDcrkBPfrEDZ1fFnNanT5zT7bZvz1L2nr1ibnKvPempd+87Y20uG5BLfrEBOfLMCOfHNCuTENytQr4mfYgaeJoYcCJPz7EB158198EkXM9qU2LAHLvHNiuTENyuQE9+sQL124Lk6EccimnTOafq+0nTViWmolWeG7IyT05LXbbjENyuQE9+sQE58swJlN0inzcQbKQbCjF1X9d2urtvY/x5ji3eWS3yzAjnxzQrkxDcrkBPfrEDZNe6lGOSR80y8XTXCdTUDbVeNWH39zn0ZW2OlS3yzAjnxzQrkxDcr0KB1/DYDbtoO0ql7T191stzrfnU24yy1Ke6FoQYHtT23S3yzAjnxzQrkxDcrUHbf49fJ+Tv63HTRjjH2+vw8Y/qdUsXqEt+sQE58swLVJr6k2yS9JOmHkl6X9Plq/92STks6L+l5SVu7D9fMUmhS4v8W2BsRHwXuBfZJegB4GngmIu4B3gcOdhemmaVUm/gx9Ztq89bqEcBe4IVq/1Fg/0ZPPm/Fm5xWORlqNZ42hlw9yNLp62/YqI4vaYukV4FLwAngZ8DliLhSveQCsLObEM0stUaJHxG/i4h7gV3A/cBHmp5A0iFJE0mTljGaWWIbatWPiMvAKeBBYJukq/0AdgFv3+Q9RyJiJSJWForUzJJp0qq/JGlb9fwDwMPAWab/ATxWvewAcCxFQDnVU3Nqb6gzbxXYnOPP6e+ck77+Zk167u0AjkrawvQ/im9GxIuS3gC+IelfgVeAZzuL0sySqk38iPgRcN+c/W8yre+b2ci4555ZgZz4ZgUa3ei8VJos1ZVbg9hmMvZrO7ZZdWe5xDcrkBPfrEBOfLMC9Zr4y8vLC3dOSNUxpUkHkjadTPp6T07GHn+JXOKbFciJb1YgJ75ZgUa3kk6TOmSq1Xbqvuufd54Uq7CkiLWJVH0XxvYd9mbT5n5xiW9WICe+WYGc+GYFcuKbFUh9NsxIyqYVaKgGNZtqM8hl7ANjujLnutTezC7xzQrkxDcrkBPfrEDZDdJpMltsikeKgSVdrV6T6ph1xxly9Z0mA61yXmVp7FzimxXIiW9WICe+WYFGN0in6XFSHLeNFN81p6q/jr0e3GaQ1GbT1X3rEt+sQE58swI58c0K5MQ3K1Cxg3SaaNO4VNcYM7YGqRIb1HLS5p7zIB0zm8uJb1YgJ75ZgbJbLTfFBBmpjptixt+x66tO79WKp/q651zimxXIiW9WoMaJL2mLpFckvVht3y3ptKTzkp6XtLW7MM0spY2U+E8CZ9dsPw08ExH3AO8DB+sOkGK13LZSTDjRZiKIuglAmjxK0GYV5M14rfqabKRR4kvaBXwc+Eq1LWAv8EL1kqPA/i4CNLP0mpb4XwI+B/y+2v4wcDkirlTbF4Cd894o6ZCkiaTJ6urqQsGaWRq1iS/pE8CliDjT5gQRcSQiViJiZWlpqc0hzCyxJt/jPwR8UtIjwG3Ah4AvA9sk3VKV+ruAt7sL08xSqi3xI+KpiNgVEXcBjwPfjYhPAaeAx6qXHQCOdRblQLqYmbeNnBqxUsXS5hh1s+7m1OGnyXVq07Cb6nde5Hv8fwD+XtJ5pnX+Zxc4lpn1aENddiPie8D3qudvAvenD8nMuuaee2YF2hQTcXS18m1fdek28edUn93s2gwg6ntVoqtWVlaYTCaeiMPMbuTENyuQE9+sQNlNxFGnq7p4TgM8mtT5U8Q7ZDvBUJN4pmo/SdGOlHJ15Y1yiW9WICe+WYGc+GYFcuKbFSj7xr2uloweqnEsxey+XerrugzV6WrIjlpdnWe2A08TLvHNCuTENyuQE9+sQL0O0llZWYnJZLLua/qoB807z5B1vzptY+tiIEmq69Skjt9Vh5cCzuNBOmZ2Iye+WYGc+GYFcuKbFSj7GXhSNYj01bCy2Y1tpqIuzjtPV43FLWdncuOemd3IiW9WICe+WYGyG6TTpqNNE3XHadMhJvdOJ110Wkr1+3Q1y00X5x37oKl5XOKbFciJb1YgJ75ZgbL/Hr+NJnXXoeqYTQzZ56DNwJ4U91Cq46YaeLToMbvS8N729/hmdiMnvlmBnPhmBXLimxVo0A48OS2H1VcHmCbvyXnG2a4aTvtqyOvyOCn01dDoEt+sQE58swI58c0K1HcHnlXg58AfAe/1duLFjClWGFe8Y4oVxhHvH0fEUt2Lek38ayeVJhHRbK2fgY0pVhhXvGOKFcYX73r8Ud+sQE58swINlfhHBjpvG2OKFcYV75hihfHFe1OD1PHNbFj+qG9WoF4TX9I+ST+RdF7S4T7P3YSk5yRdkvTamn13SDoh6Vz17+1DxniVpDslnZL0hqTXJT1Z7c813tskvSTph1W8n6/23y3pdHVPPC9p69CxXiVpi6RXJL1YbWcb60b1lviStgD/Dvw1sAd4QtKevs7f0FeBfTP7DgMnI2I3cLLazsEV4LMRsQd4APh0dT1zjfe3wN6I+ChwL7BP0gPA08AzEXEP8D5wcMAYZz0JnF2znXOsG9JniX8/cD4i3oyI/wW+ATza4/lrRcT3gV/O7H4UOFo9Pwrs7zWom4iIixHxg+r5r5neoDvJN96IiN9Um7dWjwD2Ai9U+7OJV9Iu4OPAV6ptkWmsbfSZ+DuBX6zZvlDty932iLhYPX8H2D5kMPNIugu4DzhNxvFWH51fBS4BJ4CfAZcj4kr1kpzuiS8BnwN+X21/mHxj3TA37m1ATL8CyeprEEkfBL4FfCYifrX2Z7nFGxG/i4h7gV1MPwF+ZOCQ5pL0CeBSRJwZOpau9Dke/23gzjXbu6p9uXtX0o6IuChpB9PSKguSbmWa9F+LiG9Xu7ON96qIuCzpFPAgsE3SLVVJmss98RDwSUmPALcBHwK+TJ6xttJnif8ysLtqGd0KPA4c7/H8bR0HDlTPDwDHBozlmqrO+SxwNiK+uOZHuca7JGlb9fwDwMNM2yVOAY9VL8si3oh4KiJ2RcRdTO/T70bEp8gw1tYiorcH8AjwU6Z1u3/q89wN4/s6cBH4P6Z1uINM63YngXPAfwF3DB1nFeufM/0Y/yPg1erxSMbx/hnwShXva8A/V/v/BHgJOA/8B/CHQ8c6E/dfAC+OIdaNPNxzz6xAbtwzK5AT36xATnyzAjnxzQrkxDcrkBPfrEBOfLMCOfHNCvT/C6omFjUX9rcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.naive_bayes import GaussianNB\n",
    "clf = GaussianNB()\n",
    "clf.fit(X.reshape(-1,1), Y)\n",
    "result=clf.predict(X0.reshape(-1,1)).reshape(size[0],size[1])\n",
    "plt.imshow(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To improve our segmentation, we filter the input image with a median image filter and add this information as a second sample vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "l_median=[]\n",
    "for i in range(0,10):\n",
    "    imageName='data/sklearn/im%d.nrrd'%i\n",
    "    im=itk.imread(imageName,itk.F)\n",
    "    median=itk.median_image_filter(im, radius=3)\n",
    "    arr = itk.array_from_image(median)\n",
    "    l_median.append(arr)\n",
    "\n",
    "M0=l_median[0].flatten()\n",
    "M=M0\n",
    "X0=np.concatenate((X0.reshape(-1,1),M0.reshape(-1,1)), axis=1)\n",
    "for i in range(1,10):\n",
    "    M=np.concatenate((M, l_median[i].flatten()))\n",
    "\n",
    "X = np.concatenate((X.reshape(-1,1),M.reshape(-1,1)), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7efe14112dd0>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD8CAYAAABXXhlaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAC6BJREFUeJzt3UGoXOd5xvH/U9muAyHYToMQkqldbBq0aGx8MQ7uoig1qE6ItTDFJhQtBNok4NBA6rRQCHQRb+Jk0Y2ITbQIsVMnIONNcRVBKBTZcuyktkUixRAiI1uUWCTZpFXydnGPzM31le/ce8/MndH7/8FwzzlzxufFM4/e831zZiZVhaRe/mi7C5A0ewZfasjgSw0ZfKkhgy81ZPClhgy+1JDBlxraUvCT7E/ykyRnkzw6VlGSpiubvXIvyQ7gp8B9wDngReDhqnr9fR7jZYJTctddd213Ca299NJL213Cu6oq6+1zzRb++3cDZ6vqDYAkTwEPAFcMvqbn1KlT211Ca8m6WZsrWznV3w38YsX6uWGbpDm3lY4/kSSHgcPTPo6kyW0l+G8CN69Y3zNs+wNVdQQ4Ao7xpXmxlVP9F4Hbk9ya5DrgIeDZccqSNE2b7vhVdSnJ54B/B3YAT1bVa6NVJmlqNv123qYO5qn+1PiFKttrnmb1J3k7zyv3pIamPquvrbObz7/Vz9E8nQGsxY4vNWTwpYYMvtSQY/w55Jh+8a31HM7TuN+OLzVk8KWGDL7UkGP8GXP83tckz/2s5gHs+FJDBl9qyOBLDRl8qSEn96bMyTxtxKw+7GPHlxoy+FJDBl9qyDH+yBzTa0zTGvPb8aWGDL7UkMGXGjL4UkNO7o1s9eSLk30a03qvp6WlpYn+O3Z8qSGDLzVk8KWGHOOPzDG9FoEdX2rI4EsNGXypIYMvNWTwpYYMvtTQusFP8mSSC0leXbHtpiTPJzkz/L1xumVKGtMkHf+bwP5V2x4FjlfV7cDxYV3Sglg3+FX1A+CXqzY/ABwdlo8CB0auayFU1Xtu0iLY7Bh/Z1WdH5bfAnaOVI+kGdjyJbtVVUmu2OqSHAYOb/U4ksaz2Y7/dpJdAMPfC1fasaqOVNVSVU32QWFJU7fZ4D8LHByWDwLHxilnvqw1hnc8r6vBJG/nfRv4L+DPk5xLcgj4CnBfkjPAXw/rkhbEumP8qnr4Cnd9YuRaJM2IV+5JDRl8qSGDLzVk8KWGDL7UkMGXGvJbdgdekKNO7PhSQwZfasjgSw0ZfKkhgy81ZPClhgy+1JDBlxoy+FJDBl9qyOBLDRl8qaG2H9LxQznqzI4vNWTwpYYMvtSQwZcaMvhSQwZfasjgSw0ZfKkhgy81ZPClhgy+1JDBlxoy+FJDBl9qyOBLDa0b/CQ3JzmR5PUkryV5ZNh+U5Lnk5wZ/t44/XIljWGSjn8J+EJV7QXuAT6bZC/wKHC8qm4Hjg/rc6uq/uAmdbZu8KvqfFX9cFj+NXAa2A08ABwddjsKHJhWkZLGtaExfpJbgDuBk8DOqjo/3PUWsHPUyiRNzcTfuZfkg8B3gc9X1a+SvHtfVVWSNc+fkxwGDm+1UEnjmajjJ7mW5dB/q6q+N2x+O8mu4f5dwIW1HltVR6pqqaqWxihY0tZNMqsf4AngdFV9dcVdzwIHh+WDwLHxy5M0DZOc6t8L/B3w30leGbb9I/AV4DtJDgE/B/52OiVKGtu6wa+q/wRyhbs/MW45kmbBK/ekhgy+1JDBlxoy+FJDBl9qqM2v5a680hD8tVz1ZseXGjL4UkMGX2rI4EsNGXypIYMvNWTwpYYMvtSQwZcaMvhSQwZfasjgSw0ZfKkhgy81ZPClhgy+1JDBlxpq8w08q63+Rp61+C09ulrZ8aWGDL7UkMGXGmo7xl/N8bw6seNLDRl8qSGDLzVk8KWGDL7UkMGXGlo3+EmuT/JCkh8leS3Jl4fttyY5meRskqeTXDf9ciWNYZKO/1tgX1V9DLgD2J/kHuAx4PGqug14Bzg0vTIljWnd4Ney3wyr1w63AvYBzwzbjwIHplLhjCR5z026Wk00xk+yI8krwAXgeeBnwMWqujTscg7YPZ0SJY1touBX1e+q6g5gD3A38NFJD5DkcJJTSU5tskZJI9vQrH5VXQROAB8Hbkhy+Vr/PcCbV3jMkapaqqqlLVUqaTSTzOp/JMkNw/IHgPuA0yz/A/DgsNtB4Ni0ipQ0rkk+nbcLOJpkB8v/UHynqp5L8jrwVJJ/AV4GnphinZJGtG7wq+rHwJ1rbH+D5fG+pAXjlXtSQwZfashv4Hkfqy/i8Vt6dLWw40sNGXypIYMvNWTwpYYMvtSQwZcaMvhSQ76PvwG+r6+rhR1fasjgSw0ZfKkhgy815OSetOA2843QdnypIYMvNWTwpYYc42/BWmMrL+rRNI31C092fKkhgy81ZPClhhzjj8wP8mhM0/rVZju+1JDBlxoy+FJDBl9qyMm9KXOyTxsxrcm81ez4UkMGX2rI4EsNOcafMT/Yo8tmNZ5fix1fasjgSw1NHPwkO5K8nOS5Yf3WJCeTnE3ydJLrplempDFtpOM/Apxesf4Y8HhV3Qa8Axwas7BOkrzvTYtnved0u5/XiYKfZA/wSeAbw3qAfcAzwy5HgQPTKFDS+Cbt+F8Dvgj8flj/MHCxqi4N6+eA3Ws9MMnhJKeSnNpSpZJGs27wk3wKuFBVL23mAFV1pKqWqmppM4+XNL5J3se/F/h0kvuB64EPAV8HbkhyzdD19wBvTq9MSWNat+NX1Zeqak9V3QI8BHy/qj4DnAAeHHY7CBybWpXNzdvEkN5r0Z6frbyP/w/A3yc5y/KY/4lxSpI0bZnl5aJJvDZ1JF7mO1/mqctX1brFeOWe1JAf0llQk3QYzwrea5IPSc1T954WO77UkMGXGjL4UkOO8a9i641Vx5oDmNUXim5mXmOSx3QY069mx5caMvhSQwZfasjgSw05udfYtCa1tnOyrONE3WbY8aWGDL7UkMGXGjL4UkMGX2rI4EsNGXypIYMvNWTwpYYMvtSQwZcaMvhSQwZfasjgSw0ZfKkhgy81ZPClhgy+1JDBlxoy+FJDBl9qyOBLDRl8qSGDLzVk8KWGZv1LOv8D/Bz4k2F5ESxSrbBY9S5SrbAY9f7pJDtlWr9l/r4HTU5V1dLMD7wJi1QrLFa9i1QrLF6978dTfakhgy81tF3BP7JNx92MRaoVFqveRaoVFq/eK9qWMb6k7eWpvtTQTIOfZH+SnyQ5m+TRWR57EkmeTHIhyasrtt2U5PkkZ4a/N25njZcluTnJiSSvJ3ktySPD9nmt9/okLyT50VDvl4fttyY5Obwmnk5y3XbXelmSHUleTvLcsD63tW7UzIKfZAfwr8DfAHuBh5PsndXxJ/RNYP+qbY8Cx6vqduD4sD4PLgFfqKq9wD3AZ4f/n/Na72+BfVX1MeAOYH+Se4DHgMer6jbgHeDQNta42iPA6RXr81zrhsyy498NnK2qN6rqf4GngAdmePx1VdUPgF+u2vwAcHRYPgocmGlRV1BV56vqh8Pyr1l+ge5mfuutqvrNsHrtcCtgH/DMsH1u6k2yB/gk8I1hPcxprZsxy+DvBn6xYv3csG3e7ayq88PyW8DO7SxmLUluAe4ETjLH9Q6nzq8AF4DngZ8BF6vq0rDLPL0mvgZ8Efj9sP5h5rfWDXNybwNq+S2QuXobJMkHge8Cn6+qX628b97qrarfVdUdwB6WzwA/us0lrSnJp4ALVfXSdtcyLbO8Vv9N4OYV63uGbfPu7SS7qup8kl0sd6u5kORalkP/rar63rB5buu9rKouJjkBfBy4Ick1Qyedl9fEvcCnk9wPXA98CPg681nrpsyy478I3D7MjF4HPAQ8O8Pjb9azwMFh+SBwbBtredcw5nwCOF1VX11x17zW+5EkNwzLHwDuY3le4gTw4LDbXNRbVV+qqj1VdQvLr9PvV9VnmMNaN62qZnYD7gd+yvLY7p9meewJ6/s2cB74P5bHcIdYHtsdB84A/wHctN11DrX+Jcun8T8GXhlu989xvX8BvDzU+yrwz8P2PwNeAM4C/wb88XbXuqruvwKeW4RaN3Lzyj2pISf3pIYMvtSQwZcaMvhSQwZfasjgSw0ZfKkhgy819P+1X4r2XstnnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "clf.fit(X, Y)\n",
    "result=clf.predict(X0).reshape(50,50)\n",
    "plt.imshow(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Typical processing\n",
    "\n",
    "* Image preprocessing\n",
    "    * Bias field correction (n4_bias_field_correction_image_filter)\n",
    "    * Noise reduction (smoothing_recursive_gaussian_image_filter)\n",
    "    * Feature computation (Texture, Wavelet, Edge detector)\n",
    "* Converting ITK data to NumPy and organizing the data as needed\n",
    "* Train classifier\n",
    "* Use classifier on new data\n",
    "* Convert classifier result to ITK data\n",
    "* Apply some post processing filters\n",
    "    * Fill holes (binary_fillhole_image_filter)\n",
    "    * Smoothing (median_image_filter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Enjoy ITK!"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}