[5d12a0]: / tutorials / MultiMetricRegistration.ipynb

Download this file

175 lines (174 with data), 34.2 kB

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ANTsPy and ANTsR inherit the ability to do multi-metric registration.\n",
    "\n",
    "Such registrations assume that the feature images are in the same physical space and at least roughly aligned.\n",
    "\n",
    "One may then use registration to optimize a multiple similarity metric objective function as in:\n",
    "\n",
    "https://doi.org/10.3389/fninf.2014.00044\n",
    "\n",
    "https://www.ncbi.nlm.nih.gov/pubmed/18995188\n",
    "\n",
    "First, import ants, read some images and create some features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/stnava/.pyenv/versions/3.8.1/lib/python3.8/site-packages/ipykernel/ipkernel.py:283: DeprecationWarning:\n",
      "\n",
      "`should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import ants\n",
    "image = ants.image_read(ants.get_ants_data('r16'))\n",
    "image2 = ants.image_read(ants.get_ants_data('r64'))\n",
    "aff = ants.registration( image, image2, \"Affine\" )\n",
    "g1 = ants.iMath_grad( image )\n",
    "g2 = ants.iMath_grad( image2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform a baseline registration with a single feature and create a couple of new metrics.  Each metric is defined by a name (\"CC\"), the input fixed (image), input moving (image2), a weight value (e.g. 2) and a sampling parameter ( for CC this defines a radius of 9x9 e.g. 4 extra pixels on all sides of the center pixel.  Five entries are needed in total."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg1 = ants.registration( image, image2, 'SyNOnly', initial_transform=aff['fwdtransforms'][0], verbose=False )\n",
    "demonsMetric = ['demons', g1, g2, 1, 1]\n",
    "ccMetric = ['CC', image, image2, 1.5, 4 ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Append the first metric to the metric list. In actuality this means that reg2 will be driven by both a demons metric and the default metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics = list( )\n",
    "metrics.append( demonsMetric )\n",
    "reg2 = ants.registration( image, image2, 'SyNOnly',\n",
    "    multivariate_extras = metrics, initial_transform=aff['fwdtransforms'][0] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add a third metric and run this new registration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics.append( ccMetric )\n",
    "reg3 = ants.registration( image, image2, 'SyNOnly',\n",
    "    multivariate_extras = metrics, initial_transform=aff['fwdtransforms'][0] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quantify the results in terms of mutual information of the registration results using the original image intensity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.17961406399219665\n",
      "-0.76977124858361\n",
      "-0.763288742703357\n",
      "-0.8067789848956768\n"
     ]
    }
   ],
   "source": [
    "print( ants.image_mutual_information( image, image2 ) )\n",
    "print( ants.image_mutual_information( image, reg1['warpedmovout'] ) )\n",
    "print( ants.image_mutual_information( image, reg2['warpedmovout'] ) )\n",
    "print( ants.image_mutual_information( image, reg3['warpedmovout'] ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABVYElEQVR4nO19WXdb2XH1xjzPAwdwlEiJGtvdbbvdXk6c5TzkLT80b3nzWnlwJ3bsniS11BLVGklJJEGCmOcZ3wO/Xax7CbU77TimgVNraUkCQODi8uxTVbt21XFMJhMYM2bs8pnzb30BxowZm24GnMaMXVIz4DRm7JKaAacxY5fUDDiNGbuk5v6+Jx0Oh6FyjRn7K9tkMnFMe9x4TmPGLqkZcBozdknNgNOYsUtqBpzGjF1SM+A0ZuySmgGnMWOX1Aw4jRm7pGbAaczYJTUDTmPGLqkZcBozdknNgNOYsUtqBpzGjF1SM+A0ZuySmgGnMWOX1Aw4jRm7pGbAaczYJTUDTmPGLqkZcBozdknNgNOYsUtqBpzGjF1SM+A0ZuySmgGnMWOX1Aw4jRm7pGbAaczYJTUDTmPGLqkZcBozdknNgNOYsUtqBpzGjF1SM+A0ZuySmgGnMWOX1Aw4jRm7pGbAaczYJbXvPTzX2OUxh8OByeT7zzJ2OKxnsP651xu73GY859/AHA7HBSDp56bZ+4D2vvfi651O8yv+ezXjOf8G9n0eTT9H0PFv/dy095hMJnC73RiPx5hMJhe8rdPpvPC48a6X1ww4/0o2LcS0e7nJZHIBjHze7XYjFAohEAjA5/NhMpmg3W6j3W6j3+/Lz2kg8o/T6ZT3GY/HAthp12YHPt+Drzfg/duZAedfybjIpy1uPq5fw9e53W6Ew2EsLi7iypUrSKVSiEajGI/HqNfrqNfrKBaLOD09RafTQbPZRKvVwmg0AnAGRg1Oh8OBQCCA4XCIXq9nAbAdsPqa6GUBA9C/lRlw/hVsWhjqcrksj9nBQBBlMhlcu3YNa2triMfjiEajiMfjmEwm6PV6AkiCslKpYH9/H4eHh6jX6xgOh/J5oVAI4XAY6XQanU4HL1++FEAGg0EAQLvdFmB/34Zi7P/eHN/3i3A4HOa39D+w78sR6c0YhvIxp9MpIezS0hJu3LiB9fV1JBIJ+P1+RKNRuFwuAdVwOES328VwOITX64XD4cDBwQHu3buHBw8eoFKpADgD389+9jNcv34dqVQK/X4fx8fHGI/HaDab8Pv9iEQiyOfz2N3dxfHxMTqdzgVvboD617fJZDKVBTSe80eY9i4aiPSEBJI9v3Q6nfD7/ZhMJnC5XAiHw8hkMohEIohEIkgmk1hZWcHi4iK8Xi/cbjfcbrd4Q4/HA5fLBb/fD4fDAb/fD7/fj1QqBa/XC6/Xi1evXqHX6yEej+Pjjz/G7du3EQ6HMZlMMB6P0e12MRgMEIlE4PV6kc/nsb29jS+++AJPnz5Fp9ORz+Nn6/DX2P+dGXD+Baa9jN1rMowdDodwuVyIRCLIZDJYXFyEy+VCMpmUx0KhkISWsVgMgUBA8j4CJRAIIBgMYjQawefzyR/gDERXr16F0+nE6uoq2u02IpEItra2kM1mJYTtdDrw+XwYDAbyHRYXFxEMBpFIJLCzs4MXL17g5cuXKJVKlk2GOao9R+bz9nDd2F9uBpw/wqbVELUHdblcskidTic8Hg8ymQw+/PBDrKysIJFIIBQKwefzIRQKwe/3o9lsol6vC+DdbjecTidcLhdisRiSySRCoRBcLhfcbjdGoxEmkwlarRaGw6GFpQ2FQsjlcojFYnJdwWAQPp8P7XYb5XIZw+EQ2WzWAvStrS1sb2/j0aNH+PLLL3F0dIRutyted5rpTUlHD/o+GftxZsD5I02zmqPRCA6HAx6PB5FIRMLIXq+Hfr8Pr9eL1dVVLCwsIJvNYmFhAV6vFx6PR4DhdDrR6/UwHA4F0Hw+kUggnU4jEolI/tlutyVEpSd88eIFnj9/jp2dHSm/DAYDuN1uDAYDAbvX68VwOITb7UYwGMRgMIDH4xFwMaT+/PPPcXh4KACd5iHtpJYGJ2AA+peYAeePMC5Ckjt+vx+bm5u4fv06lpaWEIvF4Ha7UavVkM/n0Wq1kM1mkc1mEYvFJHf0+XwIBoPwer3o9Xri+fhcLBZDJBIR1tbn82E0GgkgR6MRer0eut0ujo6OcHx8jNFohHA4DADo9Xrwer0SHvv9fql7ckMgYF0ul4TMuVwOhUIBS0tL4s3b7bbFIzJqsNdHDYn0v2cGnD/SdGnE4XDg6tWruHPnDjY2NpBMJhEMBtHr9fD8+XO8efNGgMxFTWB4vV74/X4Eg0G4XC5MJhMEg0GkUikkk0nE43FEIhF4PB4Mh0P0+30RIYxGIyGIhsMh4vG4vJ5M7ng8xmg0wng8hs/nk+ugd2Qu2e12BZzj8RihUAgLCwvodrtoNBooFAoSCTD3tG9S3Fz4ecb+MjPg/BFG7+B2u7G8vIxf/OIXuHr1KrLZLDKZDOLxOBwOB7xeL1ZWVgAA/X5fwEDFDskeerRIJILJZILFxUUkk0kkEglEIhF5XpdhmNe63W54vV6EQiGsr6/D6/Uik8kgHA4LQEejEaLRKLxeL4AzcsnpdArxRKD3ej20Wi2USiWMRiMsLi4iFothNBqhXq/j+PgYe3t7qFQqGAwGllyU/7arj0wO+uPNgPNH2ng8RjAYxN27d7GxsYFAIADgbEH6fD643W50u12pX3a7XfR6PQBnrKnH4xEvyjLIysoKHA4HUqkU4vE4gsGgxbtpyR4JIeCMEfb5fFhdXUUwGEQoFBLAud1u+P1+IaDcbreE1V6vV2qm/X4frVYL9Xod1WoVvV4PsVgMS0tL8Hg8aLfb2NrawsbGBl68eIFisYjxeIzBYCBet9frodfrWbyolhca+5+ZAeePMIZ0FA1sbm5awstOp4NAICALdjweIxwOw+/3o9frYTAYiOcDIABifhgMBhEIBCR01l6TJE+320W/35fPYN0zHA4jHA6LR/X5fAgEAvKYy+WCx+ORf49GIyG1tBTQ7XYjEomI9261Wuj3+1hcXMTW1hZarZYQU6VSCaVSCW/evEG1WpXXaq2wKbX8z82A8weYfUE5HA4kEgl88MEHWF5exsLCAhwOB4bDIUajERqNBgDIAgfOQleXy4VWq4VWqyXsLkESDAYRDAaFECIwCR7gTITAXJFhMbW0OsTk55JsIuvr8XjkfVwul4giQqGQPAcAPp8PXq8XiUQC0WhUnh8OhxgOh0gkEpbcuVar4fT0VED67NkzvHr1Cs1m84KQ3h4FGHu/GXBOsWm7vWYn3W43UqkUFhYWEAqFpJ7Y6XTQarVkUabTaSGHJpMJhsMhBoOBhJEswfD9SQjZVUd8jARSIBBAv98Xj0kgk6xhaDkYDOD3+8WL0ntrVtXn84lCye/3o9vtWkJhrWgKhUIYDoey+ZAVDofDSCaTQhitrq4inU7jwYMHKBaLlntq9Ls/3Aw4/wdGIEQiEWxvb4u6hgyrx+OBw+FAs9mUmmUoFEIkEoHD4UCj0cBoNEK/3xfBOT0ePSDNXpLgvz0ej4Cz0WgI8zoej/H69WscHR1hcXFRQElAer1eJJNJpFIpYYpHo5GQSfS07XYbwJn3ZOit//Z6vQgGg1KeIdFF0Pn9fmxsbEh4/+2336JWq8n3MF7zh5sB5xSb1vDMx3w+H1KpFNbX18VjsnYYCAQkz2u32wKeZDKJWCwGl8slXSXsQiFRM01wbidS7OWJ0WgkZIzD4UCpVMLR0RGuXbuGzc1NBINBERs4nU50Oh0A5+ExAGGPo9GopZzDcJuCCF02cbvPlg2VSfTSo9FIAJ/L5USk/+jRI8lBea3G/rwZcH6P6XCWZYJgMIhcLodUKiWqGoaRfr9fFrLT6US73UaxWMRoNBJPxoUdiUSklhkKhS6EegyDCUjWIrXelaUSAMLCssbJHJYhqcvlEuleOBwWoQTLMdTuulwu0QMzJyY7TADS+9Hbk5QikQSc5dhLS0vY2trC/v4+yuWyfDedd5oQ9/1mwPk9pheg3+9HIBDA8vIydnZ2sLCwIBI5nUt6PB4kEgkJNZvNJiqVCsLhsNQaKcNLJBIWb6W952g0Qrvdljyu0+lISxdDThJJjUYDfr8fy8vLWF5extLSEpLJpHwHAnAwGMh78RqYr5I0IoNLb8vv1+v10G630el0JGTn+wPnHphi/MlkgkAggFQqhXA4jHK5fEG8MI3NNXZuBpzvMU1ghMNhrK6u4tq1a1hdXcXOzg6i0ah0kkwmEzSbTQEcPSjJEwItGAzC7/cjHo8LEcPQkZ9FD9Vut3F6eiqha7lcRqPRQCAQEI9GL81yCeV+6XRaHvN4PKLBJek0rfbo8XgQDocRCoUs35+haafTkRCaOWi320Wr1cJgMBDPq+cUUfnEsN1+X3nv+H8DUKsZcGL6BDsW0iORCG7evImPPvoIV65cweLiInw+H/r9vkWT2mw2hRFlZwgAUdewDEFWlOSK7l4BgMFggEajgXq9Lqxos9nE27dvUSgUBECJRMLiQZnvEbAkgwKBgHhMhqv2GipfzxY2Asvv92MwGAhhRamhx+PBYDAQsUKr1ZKwnUSRDl0J0kajIZuDXThvgHnRDDhxEZyaGb169Sr+6Z/+CTs7O8hms6KqYUlCz9oZj8fo9/sYDAYIBoPiGcPhMFqtFgBYPB69pv5sekkO8up0OqhUKhLi8lopvWM+yNEl0WhU8uBGoyGihkgkIoBhndMutQuHw8Lm0kMyxGW7Gn8WgGwMuoOGzwGQctHVq1eRSCRw7949YYONx/zzZsD5/03v4iQ1gsEgNjY2pNOE3oRiARI2lOLpn+WijsViCIVCorzhQtb1TX4+yywEOBU7g8FAao66zYzgGgwGIr2j9woEApIfJhIJJBIJ+WyWUnSISc+YyWSkv7TT6WA8HiMej1uED51OR/JPkl0ULhD0zMNDoRBu376NRqOBd+/e4fXr1/J9AVg2N2NWm3tw2hepnmBAb8JcjqUPyuwY7hEQVP0wbOX7U53DMoLda2jPzce63a5cI2urBBfLHMCZd3I6nRI6VyoVmXzAVrTT01O0223E43HEYjFL6ca+SUQiEYRCISn5cKOgZJCga7Va4gU5vpPAJ4nkcDgQi8WQSqUwmUywt7cn/aE0073yfptrcNqVOPbnmLfR69Bz+Hw+S9kEOG8B015Jg48hoF20rsFBgQFFA7R2u21pwmYtFTgDYL/flw2Bnk3PHxqPxyiVSuK93wdQfg/grBTidDpF7EAv3Ol0JA/l9VAmyCkN/G5UM4VCIQSDQdy8eRNffPEF8vm8YWl/gM01OAmcabNeY7EYlpeXxWPp/I7MJ4kNyursuSu9CKcRMOSzS9kI6PF4LOUHAKjX6/JZBC51rqyxOhwOKfATIOyQoVKJZY52u41KpSIRAfNmbSSuGGazn7NarUqoS5ECmVh6anp1Ek9sGE8kEuj3+1hfX8fVq1dRKpWkQ2caa6v/P8821+CkESAEYDQaxaefforl5WUkEgl0Oh1ZbBwJQkAAkI4UejUuVns+qT/PDmQtWqeaiIDlv+mB6IG73S663a6QR4FAQDpJ+F1YJwUgJRGPxyMN29wYtIyQ4gR+NwAi1aMHZ/RAYb1mn4fDoXTWZDIZBINBFItFkT0+efIEo9FINgHeE94rY2c29+DUTKvOuzg/Nh6PYzgcSv7FkgEXE0NUkiLMSTXpw3BTNyNrHS0/W7Op8XhcZstSGUTP1O/3Ua/XpX7ZbDald5RN1rxG3bMJWPNq/d15TTonpuKIjKzf70e5XJb7oHNrrVhyOByIRqMyYVBPWVheXkY0GkW1WpVNYVr/pwGpAacYFwTJl8XFRezs7IgmlexkvV4XDa3D4RCyhI+Nx2OZE6TDZR02272FzlX1dQQCAZm0x02AHocArdfrKBQKaDabMoWBXpxel4IIj8cj0xpY69R1Vv35BKxuOWOvaLvdtuTQ/HkSQgyZ4/E4AMiEwE6ng7W1NfzLv/wL/u3f/g21Wm1qKmDE8Wc29+DUlD5wJmy/evUqUqkU/H6/vIbyObsoXpdARqORTCEAMDWU1f/WhJH2pFTy8DUEJ0HJjaLdbqNarSKfz8sMITLCBB9z3slkglAohHQ6LWWhaWqdaYJ71mr9fj9isZh8Vz1oTGtx9bQF1oTj8biM+ex2u7h79y7u3bsnzK2+V8bObO7BqcM4EicLCwvIZDISylEip2V2wHk4yJxQEzZ83u4Z9P/t3sEe4umfoZGwoVcmWO3MMdVKZHWZJ1JCSODyOR1K2gUZmnEmm0zvqsN1zVSTTKJwguAcDAZIp9PY3t7Gs2fPZMPTYa0W98+zB517cALWGazRaBRLS0sAIIDjItQDuRjS0ctRyxqPx6UOOm3B2z3TNLDqfku9UFlfpDiA183ckIIHEjXMQbno9QQE+8K3h5a8Hg0WAHIv7Cw1AIkgKIhoNpvSt6rnGvE4iFgsJsPE+N72uuc8E0UGnP/fqKP96U9/ilwuJ6QGWUd6RfZH0rTwgAcPkYx53+LXj+nFpz2PfaGORiPUajUUi0XUajU0m000Gg10u124XC7LOBG+HrhY3tA1VYbKWsxOAOkwm9enf9YOosFggFqthkKhgNPTU7lGhr2cb8QmgOFwiGQyiXK5jEqlIpuH9sbzCEhtcw9OXa9cWFjA7du3kUwmsbCwgFgsBuBcQxoMBi2LkuDSXSKapQXOQ1v9eu05NfnC0FHXXrVAXdcwJ5OJDNlirZNN3AwlOXme16fLIwSaFqLTW9vZUn3d+pq1p2cNtdvtSkN6KpUSIoh5MzfBjz/+GNeuXcPDhw/x29/+Fvl83jC1Npt7cNIikQg+/fRTbG9vS4cGJxmwTzEUCllaorSIAMCFXBGYXqqxL0C7nE+zs8yD2UuaTqdFtVOr1bC9vS1zf5hj6kHVnFELwJJnkjjSxJM9rOZ1a7MrnPhvzlXiWaK8R8PhUKbxcZQmPTbPann79q2A09i5zT04GdIlk0lks1k4HGfTDvL5vLCSPp8P8XgcyWRSeiTpTfXCtntBhqZ2dpbP07Qn0l7MTibpYxxCoRCy2eyFsZkEnA5Tp5FLOk+0g2KaKF57TXvITo+oSTNt0WhUBAw8PqLT6aDf70v7Gz9nWvg/rzb34CQQEokEgLPcqVgsyiKk9I6li2QyKX2SuuiuwakXPnAeRvLz9N/8OZ3/cVqenbjRgNDejc/pmqN+vWaB9XvSQ+sNRX8X/qx+vd4E+G9NEPG99CbkcDhEc8te0Mnk/KRuvZFNA/e82tyDkyFZJpNBOp0Wyp9SvUajgX6/Lz2VnU7HckzCNKE7+yC1FM4OUL6Wi3I8HqPT6aBWq2EymSCTyUhvJn9mGvOrPbOdIdZA1KE4cLZxtFot9Ho9kd7pvFSbfm8tueM1MFfVIz91eYTXyTlD3W5XmsmbzaaINswZK1abe3CyYJ/L5WTeKnWjPISWTc/D4RDVatVSWwTORQO6X5JejF0bbFYGzr0RP1+XTk5PT/HixQusra3hzp07SCaT4ln0NdvrgTR7qUYLC4Dz3s98Po98Pg8AWF5eRiaTseSk+rPs76PDXJZGOLJEzxcis8vP5bR7ekw2oC8sLGB1dRV7e3vvLT/No809OIGzociaiQ0Gg/IY+xc5jqPRaKDVaok2VOdnem4PvRQn4BHAWiZHj6rDQafTiXK5jJcvX6LVauHnP//5hZxMg0Pnlfbw1U7wDIdDlMtlPHz4EF9++SUqlQpyuZwoh7RH172n9JjU8vI96Q05hYH6Y+A8lyeLTQ1wu922XDvnDn3yySfY39+/8F3m2eYenHpxk03sdrsCMtYwOWrS7/ejVCrh9PT0AoHR6/VQr9ctIWwgEJC63mQykVxVe1GGeu12G16vF7du3cJ3332HJ0+eYDAY4Cc/+YkcV69DzGnjJaeVaICzw5NOTk7wzTff4I9//CNev36NcDiMlZUVS5O0HRAcvUK1D+ul/L6st1arVRlmxl5UO6nE6IPiBW46vV4PGxsbyGQyOD4+luufd5t7cNIrcOEzXKXwgH2Y9AJkSlnS0GHcYDCQKXW1Wg2j0UgUQ/qP/mwuWh4A1O12EQ6HcefOHRweHmJ/fx/VahV3797FysqKpb/UHnIC5+oeDcx2u40XL17gv//7v/H06VOUy2UpzayursqhS3oANd+XGw57Q5mj815xGBlDWm421NgyXKfH1UcWUhjfaDSwtLSEzc1NAafd7AzxPNhcg5MsIie3A+f1yNFoJBMFGMLxOZ7aRTaXM2vb7bZ4EnoXPYjLngvy85h/aZVMJBLBlStX4Pf78e7dO/znf/4ncrmcnAOaTqdlrg9DWYJUe7ZKpYLnz5/js88+w6tXrzAajZDL5bC+vo7NzU0px3Q6HQm9tTfu9/uo1WqoVqsiG6TnZz8pT01jiKrvLZvCGeayy4VzgCeTM2F+NpvFjRs38PjxYzSbzfd68nmyuQKnPfyz53oM7+hF6f302ZMUBLDD3+/3i9dhKWQymYi+lWEsw2QuOuZ0fL0WM/DzeWZnKpXCq1ev8OrVK+zu7iKXy+GDDz7A2tqaCCPskrvhcIjDw0P86U9/wv3793FycoJgMIj19XVsbW0hl8uJF2cDtj5jlBuKnoZAzw5ANi1GFbr+yrKJnq/E78tZQxRK8P5Wq1UsLCwgnU7LSFBeCzCfYe5cgROwlgXoZcjKsmuCND9DVg65YinA4/EgGo3KZLt6vY7T01PJK5lj0sPqU7zsfZ30epyFS20pa4HUzV67dg3BYBCPHj3CV199hZcvX+LnP/85PvroIynyU97XarVweHiIP/zhD/jiiy9QrVaxtraG69evY21tDel02hJeA9PFCLxGgp4SvV6vJx0uVCF5PB6EQiHRI+tcXiuodIrgcrlQqVQAnDUTJBIJrK6u4vDwUOYUzbPNDTg1OaFB4XCca1D5fKvVEs1prVaTWicAiySOahc9YJkAZGsVcyzWEe3X43KdHa1HwYLX67WwngBEpbS6uoper4dvvvkGe3t78tl3796VNi6WYj7//HM8e/YM1WoVV65cwZ07d+RwI3pm5o66lKNZWobZPIWMHTF2YQK9vA5ZKdDQLDXrvyypMFzmnKPBYGAB87zbXIDTXr7QhXsuqOPjYyn8c8HRW3o8Hjn4B4B0n2hvZx9yxfomu1V0axU9LBctz8MkIAk0fZARyahr164JCZPP53H//n2kUilsb2+j1Wrh888/x7fffovvvvsOvV4Pm5ubuHnzJjY2NmQKvRY+kLghUO3aWXrjcDiMdDoNp9OJRqMhQNLvR1aX+l478Fkvbrfbkp+Xy2URzp+enqJarV5QIc2rzTQ4+cu1a0VZX+TYy2QyKd5vPB5LvZJhmr3bRDdV65Yrkj8UJPDoBD0lnbkryyoEIqfZcRYsNwISPpTzUTBxeHiIer2Oo6MjvHr1Cqurq6jX69jd3cXe3h56vR4WFhawtbWFzc1NyQWn9Uxyg+A1MsTlfaNX5Pdi6xwByhqoFl5wMyHB1Gq1UC6XZd4ttbaVSgU+n0+ilYWFBbx7906m882zzTQ4NSNKIwgymQw++ugjLC4uIhqNIh6Py+BodncwJOVCY7hFr0hvoU/e4ut0AZ6gdjgc8noufi2i59kk9n5GuxeJRqOIRqPicUqlkuS9HC4dCAREzM9rpbe2HwPB69TkjWaBybLyj/2MTW4kHOPJcF+rpWq1Go6Pj+WMUIfDIRHJYDBAOBzG0tISarUanj9/jkqlMtdeE5hxcALndT+GWKFQCFeuXMEHH3yAjz/+GGtrazLRbjKZIBwOC/1PT0Hygywkn+OAZYZzBBxfRwZ3MBgIyPVGwZ+lh3E6nRIe68Wvf0YTSDzWjx0eetylz+dDOp22kD/sCeXoEu3x7Ycb2T1nMBhEs9m0jNXUOTSv0T4LSOev3IC46fl8PmGCyRRzk9Q5+bw2Xs8sOPXi1h0csVhMgLm5uSm5ET0dR1xygh0XqK536kl4XKSaUCHp0W635aiC8Xhskb+xhsg8j2DT83l0XycAUTCxrgjAkrdSxwoAoVAIS0tLCIfDosiht9KdISyh+Hw+2UjYc0lvT68aDAYt4nXdn6ln9trvH1vydLTB0J5sNMUf0WhUasj8zvNqMwtOvdPSE/h8PiwsLGBlZQULCwtCBlEJBEAWGOcEUd2iPaFeuASUnlBAAOpyCHWoXJS8LpZquEHQs+n6J2uLWiPb6/WETeXRCrocwhCUjzMK0OJ2fZSC7lZh6Yj/53fjAUo+n0++B6OHaDRq6SHlBsM/BKqeVMh5utxUEokEWq0WIpGIjFaZR49Jm1lw0nQ7VSaTwY0bN7CwsCA9mQQic0N6x2azKYvdPsaDcj0CgTlZt9sVxlV7Ma3b1fmafQHTKHSgNwcgYCA4+TkEm27W1qGlzld1qMhNgOM/dQiuJYnac+qRlzoE1tfD12txAiMJ6ofZesfvypCfOTtPDAesp5DNG4M7s+C01zM9Hg/W19dx5coVpFIpIX7oFQgOfcQ7FygAC1B0aEtPyTayXq+HUChkWdQ6xNYgt6uE9IIkO8uQUTdj82hAvhaACM7tNUJ9QC8Xt6696qFaWtxOI2D11Hd6a71RcboBj3/QOly+jz7Al7k9NzGGy7xeboi8ZnuDwjwAdCbBad9dGdIuLi4ilUohEolIOAfAMhGARy80Gg3poNBFdO2ldLjKaXMUDDidTslfmeMxN9QkiX0KAc3O+HLhl8tlHB4eolQqyfko/B7Utuoh1IB1EgNJGZZGCGYNHn3/9KbADYeklSa+WCKh8H84HErJhaUpr9crHprfdzgcijSQDeC6Hc/eHjft9zurNpPg1L9M/iLj8ThWV1eRSCQQCoWEsSSxoSeYc5FxUdJ7URnEME1/nvaUOmwkULWGV/9tX3jaCHBeV71ex8HBAY6OjlAsFqWLhCSSDiM5ja/T6QjBQmCSHbbLCLW0UT/O+0DSh10sZG0JWEYATBcY+jIkBs6F/gxx2XY2mUzQaDSwt7eHt2/fynvr36P99zvrNpPgBKzA9Pl82NzcxNWrVxGPxyVv08OMWV9kwT8cDl84m5LelbmlFitQ1scZrfqQXXofXZ6wX6v2qFo4wdCPoW2j0ZBZr9FoVD6fXogeqtfroVqtolarWUJR/lt/d+CcbCJAGT7yeYajujZL4BKgiURCCDXdt8qogRtHt9tFs9lEvV5HrVaTDXF/fx8PHz5ErVaTqETfs2ljO2fZZhacwPkc2Fgshs3NTSQSCfEa1HPyl03PQ89HJQzzsm63i3K5jFqtJi1d0WjUIs8jkcEByvTODBe1RwLOZ/zYCRZ7JwY3C/ZysgzDcZixWEwAyeHSw+EQ9XodpVJJcmuG2MD5SEv7Ibr6c7Un1XNn9XQ/fUCvvlZ9DL3uUuGGVygURIbIXDObzWJnZwdutxu7u7s4OTmRCMV+7+bBZhKcdnaSC2YymcjxAPQUurRAUOoQlY/rXZ5sKgvqZGtZq5wGBO0BpuVQdtIIsBJPoVAIqVQKy8vLSKVSKJfLcszewsKCDANLp9OIRCIoFouoVCool8tIp9PIZrOWei8/QxNIWs+rgalF7jrv5vVppRQjB3pO7V3p/fSxEmRqPR4PFhcXsbS0hNXVVUSjUXz22WeWiRP6muYBpDMBTh3C6l8a/9/r9XB8fIyDgwPJORkK8uf1sXlUxRCk4/FY5HD6hC+2gTEPY22UBI19KgFDRT6nQ0tdxiEQ9OLnoOu1tTVcvXpVAMvhXKlUSqYIRqNRFItFOReT3l/nu3o4mc6JtWyP16xBDVhroVo3y7Y1CjD0d9N65FAohEwmA6/XK+em6Fw9m81KuatUKsnvktczD8AEZgScgJVQYYhIa7Va+P3vf4/PP/8cgUAAN2/exO3bt6XLIhqNIpVKiZfijk+yBThbkAR1rVZDvV5HLBZDIBCQk6/pBXTDtZ77A8DiQRnSankdcD5XiM/Tm7lcLqTTaezs7AiTnEqlJGydTCZYXl62iCg4EFsL8O1dOrrBXEcRfI0GJjcYShbb7TbK5TLy+bzcy1arJaAlMxsOhwFAdMw+nw+RSEQ8LZsIms2mnIxNZZXedOfFawIzAk69mPl/LrxEIoG7d+9iZ2cHfr8fjUZDKP1CoYDBYIBEIiENxJSQMXzVXSYkevx+v6XliguWnkrL1shw6mvSLCZ1rQz/phlJFS7kzc1NDIdDHBwcoNPpoNlsyrWHw2Fks1k54HZlZQWJREK0wZrkYYhKD2gXvmtGl6ZDcIKUGwBBzwOWNKnlcDgsh0LxRDY9bYHEkcPhEBKO5J3+Xc+LzQQ4adzlueMzFLx16xYymQwcDodMm6vVapIjVSoVGY3BESThcBjdbhfpdFpCQA1OTkuIRCISJrLbRNdF7bkkwcy8dxp7yu9B70P5HgAhgba2tuDz+ZDP5/H48WPs7OwgFothMjkbSO3z+ZDNZrGxsYFYLCbEj13xxOtj5w2F7TR+B/2HEUUoFEK/35f7RSlfqVSCw+FAJBKBz+cTVpllJV2esSuHWDM9Pj6W81O4mWnybh7C25kCpwYmBQWcN7uysiIh5GQyQSwWk5Og6/U6qtWqFM/pQXg2J72ybh2jx2m1WjILRy/8SCQii5B9j1xkLNtwagJglRly0dKj0AvS83EzWFpawng8xv7+Pr799lvkcjmUy2U4HA7cuHEDS0tLQlLxuvQ1amYYOCfEmGdPK7PwPkQiEQBnGxbDel4ra8ShUAiBQADpdBqJRELYbYoN9B8tYCBpx/syTfw+68AEZgycNP7iOFmPu3coFEK9XpfhUgRfMBhELBaTljKH42wQFY8A5Hvp3I/5p11aR5t2XCAAC8GjhekaKJydy9piIBAQSRxDYZZBVlZWMBwO8ezZMxQKBfj9fiwvL8tITpZdtAcCzodC87vRGIZrFQ9w3jxAr07QM8emnI8dJjx/MxgMIpvNIh6Pi1KKh+ty6mClUhGSrVKpYDgcYmVlBb/5zW/w6NEjHB8fWzY43akzyx50psBpJwvsMrRIJCKhLCVjnKLHUJD5Dhef7v6g5pYhHWWAzKl0KxnBRo+nPSsXN8Gj82WWGrTnYIjMkgS1twR5PB4XMHJsJucSMY8jOHWYqkkWzVLrmqyuv+p7zNeyz1MrgHSUQeEDVUz1eh3lclmIHyqB7DOHYrEYtre3kU6n8bvf/Q75fF5IMHvP6KzazIFT/3swGKBcLstRAK1WS9hEHho0HA6lDEHgMBTl3FqSHXw9wzMyilTn0LNpwTtDQj1dQIOEC43kjNbE6notiRXt8QDIaxcWFhCPxxEOhy9sBsA5aPg+uoWM94vhMpvPeQ12tpbXSVDzO3i9Ximx6HBUbzicHcSDfpk+MPfmPed8pmvXriGfz6NcLr+XMJtVmylwApDFz0XUaDTw/Plz5HI5lEolrKysSO1zOByiVqsBOFPxcPfmdLjBYIBgMCiHv7bbbYv6h/VNKoFYr2NNDzifdkDg6/xND/rSeTI3B+C8pKG/H70hT0FrtVoSYhI4ug3LzmJzE9PeVD/O3FwDWntc5txkinWrHA80Go1GKBaLFmab4Xq73ZYxJhyryVosABEw1Ot1rK2tYWdnB0+ePJH3tW/Cs2ozB067mmQ4HOLx48f4yU9+gkAggEqlIjuz3+8XQoKMqD0sI/vIoVScQasXBRczwcnFrEsJOnSlblTL4bREbpqGVIOCnmlapEBxO2V+9GCaObZfiwanLh3RI04Dw3g8ltGW+uQwHpXInlgSPQQegcgGc16DnuRAb87Pj0QiMtGBm9Y82MyBE7Ae8e50ng3BevPmDVKplIAwHA4jHo9bVELM5SjL43sNh0PJNwHIYUXMz/SMIHsnhg7/dNlCS9oITA5Stsv56F2B8/Y2Gq+XkcBgMJCCPgHLz9S9n/zO3CAITs0G87UEJFlcXm+r1ZLIg2QZNynmkpT18d6w3qnruvSqvBfU0/Z6PRSLRRSLRfGsjASmMbizZjMJTnvoNhgM8Pr1a9y9exeTyUQK5CzOs8+SRXJd8nA4HHL6V7VaxWh0dnQAvR8ZSJY/CBY9REt7Qi4sO9uoc1LtJWmakOGitoNJDxHTnknXKPk5fJwe3l7/5B+9YfAzNWB5vcyHmbNqUQdzTX5XEkC6bNRsNgFAaqhOpxPNZhOnp6c4OjpCv9+/0Eg+6zaT4LQzkZxSR2/XbDZRq9UQj8eRTCZlEVEnS4/X6/VEehYIBGS8I4mhdrsts2ZbrRZ8Ph/C4TBSqZSErXZiRpctaLpMweu354l8XpMsmozR4z7srVYElA6fuUno5mrdUK0/n+DkewDWYdNaWMEh3Gw8p8iAn9XtdiWftR8t6HA4hNXtdrs4ODiQ4xqYh1KzOw82s+CkUUh9/fp1mX4wmUwkh6TnZA5Kz6kHT9FDhsNh9Ho9NBoNKZR3Oh1hGznyhGUWu0dkqKsbrXUIrsGptbj8owkkChk0wQTAUgbhe2nvyHuiPbD+eXpFfr7uq+RGwPCfGxs3kNFohFKpJN/rzZs3KBQKCAQCSKVS0pBNEo0iC4axfK96vY5KpYJCoYBut2vJR2eZALLbzIFTh3LAmYLlV7/6FW7duoVsNiseUvdYMqxlbROAhdAhQTQej0UPys4L1j3pPQgaXbOjt2HexM/U3lDnpoC111N7X/4sSy720Z5aNEBPqCVz/G72Gqfu29SbhNbfsq7KDYGbHfNP3ZBdqVTw8OFDvHnzBleuXLGkAVQE6TISAEu0QoFFpVJBqVSS958nmzlw6p3V6XQil8vhww8/xPb2NpLJJHq9nnRreL1eNBoNOV6eqhuGv5rJpGTN7XYjlUqh3++jWCyiUChgNBpJmBwMBi+cGWIPafUGYg9/v+97MbzUXpOPa5Z5Wt+mBppmiHXtlGDWf7hJvK+PkiGxlu4xT4zH47h+/TquXLki6qzJZCLRCgkkfV2UR1IRpcsnOleeZWUQbabAaf+Fud1ufPLJJ9ja2sLCwoKAJRwOI5PJSJ2Q5AVDLn1UH70qd36K1SeTiRzu8+rVK7x580baysi88lrsTcv6D0GhySBNGvF5TrVzOp2iu6WeVXfQ8HECUI+z5OfbW8YIYP2YvfyiWWRdk2Vdk7XLVquFd+/e4fj4GEtLS/joo4/k3gKQnJQ1WT0FUNdmeQoZJYEk2exTEWcZpDMFTnvOFgqFsLa2hlgsJmdJ8jVcUPx3p9ORkR70sJwRpGfwkCwiG+r3+9HtdpHP56VmmkgkAMBSptCeyF7HpFcjAHRphQuV2l2KGlwuF/r9Pt69eycnQ8fjcXlf3YupG601wPg5wMWJCDrv1T2gfIxzfzgPiLXiQqGAcrmMcDiM7e1trKysSH7v8XjkyL9utytCDACW7hzOP2q1WnLIFOumJJrsrPcs2kyB025ra2sCFFooFBJShbkOFy/VNicnJ1Kvi0ajIoDn4qQ3cjqdcjpWNpuV8FeHjoD1dGbNpuockv8mYAhKLQTQZRmHw4FyuYxvv/1W3kvngLonU7OsOmS1bxCahNLfF7CegUIijJsahewc2uX1erG8vIz19XWZdkCvyk2N4auegECW2OVySR26VCqh2WyKVHKebObASSB4PB7cuXMHq6urAIBKpYLBYCCUPKV37El0OBxIJpOoVqtyQCxwdpozCRE9GpLhMT0hzxpxOBwyvoSsL72ivj5NAvFv5r3aowEXySHgjJV9+/YtvvrqK+lg8Xq9WF1dFdKJIOc18L310C0dpurr0/fTzgiTZaV0jwKNTqeDUCiEeDyOTCYjOXooFLLoizmmhCEsH2fOyetvt9tydspoNBKvCcwHazuT4GRIqzsm2ObFfk3u3vSkZA/1VAL2eLLNibVEAot1TZYKNKPJo9d1CULXK/l/u4DcnmsyjNUi9fF4jNPTUzx+/Bi7u7uWY+ADgQCWlpbks/WZJ5p9tgvi7f+2e3O7WIHekF0+ZG+DwaD0brIpgGWmSCQiGwXLJiS2ut2ueFBuNqxHMzTXR9Frj6//P0s2U+C0y95qtRrevHmD8XiMbDYrB+gwnNLegAffarIjEAjIzs5fPhcQw1YW7bnb65DQzm5yUfMadO1RM6J21lbPtCUh9OjRI3zxxRfiTV69eoVYLCZAYDiv65/83sB5mEqPqokVO2OrWVveD3pkgpVjRIPBIMLhsBxsBJzL/uxKKZJmDJ3ZQlatVtFut3FwcIB2uw2n04l4PI5YLIajo6MLQJxFYAIzBk4aPeGbN28krJ1GatgBp0seJIBo9h5CPQBLEyh22Z1+juDUeab2aHxMey/dDcJNg+Hs27dvhVCpVqt48uSJjFm5ffu2SOS0qoffRRMq2qPq72Gv9eqNS4/QZFTCxnb2u3Ij0HVZhsQABKBs6eOxFnt7ezg8PEQ+n0c8HpeJDleuXMHR0RHK5bLc11m2mQKnZiq5o/N0Z0rzuPtzseniPUGlO0to9sK9Xeygi/maRNHg0+UTXd/UXtReRgFg8Vr5fB6vXr3C/fv30el0BDg8R+X169dYXV3F4uIi0um05Jb60CBetyaA3gdK+8ahw20taqC4g2UPPdRM/35YvyRY6/W6CDrG4zGKxSLevn2L/f19HBwcIJvNygwizkRi/XPWbabACVw8aj4YDMpQLoZXVJswl9EjIzVI6bV0KKqFBHbWU3si7XF4XVqCp8NU3XitAaBrn5PJBNVqFQcHB9jf35f8mKE6o4VyuYzT01OcnJxI6YdtZPbNhtelBQ6afLKzzXpD0/m61+uVsJYSRUoMee9YE+VUPQ7yKpVKFpF8IpHA7du3sbKygpOTE3z33XcoFApYXl5GJBJBLpfD/v6+jN+cZe85c+AEzhdcrVYT+p6UPun/fr8v4dd4PLYQFPb6mRYj/DnVjL1EQQ+iQccFrr0za6f6pC8dDrKz5unTp+j1evjNb36DZ8+e4c2bNxavBgCNRgMnJyeIx+OIRqMAzlu6tPrJbvYNRZdg9Pfhe2qCjaUOpgf0kBwXqoXu9XpdQlxejyaQksmkfOfV1VW8fftWNLacaF+pVC7c+1mzmQMnF9hoNJIQkGMkWZ/TUwY4VqPZbErTsNaaEkAM1fiHi1GHuvx84FxATkCS7JjWP8lr0L2aBGan00Gj0cDp6Sn29vZQq9WwtrYGAHj9+rWAieUdp9OJo6MjOR4iHo8jEAiId9K5o93L202rivg5PFFM13P1e/L7cKPjpkOBOz0na76RSAQulwuRSMSykfEeZTIZbG1tYX9/X2YOb29vo1AoWEpes2gzB047OcNDc+yaU310ug4j7TNwRqORLEiHwyFzY7lQp9UFGQpr6V2325U+Uv4cw1x9die9YKfTQbVaRT6fR7VaRalUwsnJiUxzbzQaEuox52PYzuvodDool8uIRCJSIpqmw9UhOu+h9kr6NTqC0KZz6MnkXGvL781RLSytAOeN4gyT9ZgVfgcOCfN6vTg8PES1WsXy8jKi0ag0v+tUZJZs5sCp26LYmsR+S+ZRBIbD4ZAFRJKC5AmVNmRYmSdR4aOPudOhrp00mkwmIjtjWUAffkRWmEDlrKJyuSzi+aWlJZmwx/y53+/jk08+sShwGEpysp3uf+S16lKKJoo0K8zrJqsLWHNs3RCg7zlzzH6/j1arhUKhIPedAOP311GIlhcCF9ve2GxAUQfzW3v5Z9YAOnPgBM5DTI/HYzmqjx5sMBgIQcSQFoDI3/Ri4W5O76YPn9XKFi0QoAch+I+OjnB8fCzT+ViwTyaTcvx6o9FArVZDuVxGp9NBLBbDxsYGEomEpZdRj+EksPRIEs6BZSlCh6V6MgHvBUNPfic7Y6sXvwYQH+N3plfm7NlisYh6vS75KM+VoUcmiaSjCM1e03vqrhWH43xYtxblzxooaTMJTv7SmMtQRgZA8k0d4jGH5GBpLiA9uS4UCskZmeFw2HJyM99Pe0+O6igWi3j06BGOjo7gcrmkOO/xeGQ2DsNrkhzLy8u4evUqVldXxavac0O7YIDPjcdjLC4uykldlUoFtVpN8kQAMsWBqijtyelduRFw5qwWTPDe8vP0temwlEx5OBy2vA8AyyY3Ho+l64Z5uS4r6e4aTprndEO9+cwaSGcSnMC5NI5HoPOUK44s0aMk+Uvlzq0nILjdbtF3RiIRy4nVACzA0Eyt2+1Gu93G8fEx3r59i7dv38Ln8yGXy0koys8ja+n3+3Hz5k189NFHljAOOBfn2wkTfe26NMRTvEqlkgxk5ulobJWrVCoWQQDBxDwvEAgIM8twdJqiifkr5/9QbUVvrT0yPaH2iDxQl72bJJEAyJEOrPNyU6Rw3i4OmSWbSXBy1+V8WXoq5kpkT7mw7EQQd3Ov14tYLCYH9egZQ/wcHdZywRJI9XodBwcH0jNKvWgoFJKWNX6e0+nEnTt38Omnn8oMIl0HJRAYhtrn/Ni7RyaTieh+o9EoCoUCms0mSqUSGo2GnEvCTcLtdqNWqwmjS6F/Op3G6uoq1tbWkEwm5RpoGqQEtfbSHPvCjY4bG4FFAX2tVpNcmZ9B4I5GI2n3Y/jLcSj20tUs2UyCk8aOezYpMxTidHYtKdPzeQDIwONwOCzA1J0d9hDTXn6ZTCYolUo4PT0FAPFAOpzmn+FwiFu3buGnP/2ptJ4BsIzatIeVmi3VAgKtKuLnxeNxuN1uuZ5arSYEEpU+rIOORiMp8NN77u/v49NPP5X8mJ/J9+f1aLEDv5cOd7nJ8Oc5NpNekxuGVlXxGnlv2TNbrVYvaJBnLbSdKXDa6X37EfMaYMxv2FnPkIuekXko2VG7okbv1rpuyEZoNiOPRiM5Ck+rb0jMdLtdrK+v4x/+4R+wtLQkHkUrdPgZNPvC1wDQ90BfH9vb2BDNBc/Ni8cX0rsDkL8bjQbC4TByuZwcX6/DeH3vGd4CEB2xZnn5OFldlliA8+Hc2jOPRmdjNTkrmKSZPhZx2u9+FmymwEmPxgXDBmDtIfU5IaytaRE32VoSGAxRpy3EafI8Gqf0UZ00mUyEDSYB0u12kUgk8K//+q/Y2tqS4/F4ndwANFi1Dpd/SN7o53ldfC1LS8vLy9ja2sK9e/ekJEQRg+7GocfifSmXyzg+PkYymbRodPW95+fpe6TrzmSZGblQ8KEbrSns0KUl4PyYh93dXXz99dcWfbT+/FmymQInYP0lUbup+wNZRqGXJBuqpW160oGWxWkljAaofk7XEdn6pPNUNnf3ej2kUin88z//M27duiVlBoKLhXzmXHqeLmANeactTK1A0qxqKBTCzs4Oms2mTFHQDCujBnaaAOei/1qtJuUg/bm6bsxr57XZB2DrUJbA4/kyum7K15PFpbD//v37KBaLwobPss0kOPXiPj09lVJFs9lEo9EQJjMajUpeCkC8j14klKJNmzcLnI/y0LkUQ9pisYhqtSohMwdUjUYjpFIpfPjhh/jlL38pINALkpMFmAuzhsgBWeyZDAQCFoLKfj3A+eR6AFLO+fjjj+FyufDtt99ajngn28sZvMxHG40G9vf3sbW1JSSPPZTWtVCn0ykeUoONzCw7akhA2Uk2dr7QuAFwY+UMInveOUs2c+AEznf0fr+P3d1dfPLJJ0in0/JLJy0PwKJesXsEPq8BTwDYcxy9SMfjMRqNBvb29lCv19Hv94VcCYfDWFhYwIcffoiPPvpI6naAVfzO8Farhg4ODlAoFESpFAqFZOoASw4sHfH/Whyhv38kEsHdu3fhdrvx8OFDnJ6eynOUAVITy9RgPB7j9u3bcowFAEsaofNLbo4kn/SGp0HH+6lNf3en0yletdvtIh6PW8aCArNHBNFmEpy08XiMk5MTPH/+HKlUSnIlnnLFvI+v1QtF504ECL2olqJpMkgTPgxnyThOJmdH3a+vr+PXv/41NjY2LMfs8W82Kw+HZ8fS69OeB4MB6vU68vm8MK8ulwu5XA7pdBrpdFr6V5eXl7G4uChMLa9LE0nRaBTXr1/HeDzGgwcPcHR0hHa7LdfQ7/dRr9dRKpVkE3vz5g3W19el9ghczLcZJjOt0CUjrR/WQLXXifmeJI8ODg5wcnKCe/fuWWYi8fNm0WYSnDr8bDab2N3dxfb2NhKJhHiATqcjuRVn9JAM0jVLdoaQwgcgtU+/339hQenwlnlhIBBAMpnE7du38ctf/hJra2sXQlk7ycQJ9DwrVB8P0el0cHJygnw+b/G2mgUl4dLtdhGNRiXHZt2U1xuLxXDlyhU0Gg2USiUcHx/LBgRA8mYebFupVGSz0aokglPXWXV4DpxHKQ6HQwTxzE81mUUA9/t9Oczo5OQER0dHcuiuHh0zqzaT4NS1RvZBHh4eyqE7wHmtk+DkwiFDS5CwKbher6NarYoH44Lyer3yuVygFDAQ3NlsVoC5srIik/r4OuCMjSRJwpyPi5SLUOeXNJaIOFgrl8shkUjIxLpCoYBKpYJsNivqJi1gcLlcSCQSuH79upyteXx8jNPTU4vqh+wpmwjovXi/AatscjweSzcPe1SZP/NUcTLY7FLhhsbNiKdgU/sLnGudabMa0gIzCk4af3H1eh3Pnj2Dz+cTMTlBwdyToKMyRtcAW62WkEkAxCMxtNPCdC1bi8Vi2NzcxPLyMj766CMsLy8LW6vLOpyVWywWZbwkh0QzV6UXIpsaiUTQbrct4z4ZNuoma4JeH0lINppeKxQKIZfL4Ve/+hVcLhe++eYbvHv3TuqcDEN1szq9Pf/wHmjxgG4iYEhdKpVEQxyJROQ6eP9arZaQSOPxWGYicUPhNfH1s2wzCU47adPv9/Ho0SOMx2NEIhGkUilpqaJwADgfOKVDLO7glOBp8bZmITVDSjZ1dXUVS0tLWF1dRSqVshx9x7CQZE+j0UCj0RBJXSgUEp0r5X70tNFoVELCRqNhmabQbDZRrValhut0nk2u43clk0xRAsFFqaOWMuoBYlrNpCMT3jddV2XIzM4TMrT1eh2np6colUpIJpMi1mAUwQ2SPbP8rv1+H4eHh3j06JHMGJ5llpY2k+AErCzgZHLW8fH8+XNcu3YNy8vLEiJyR2eh3+l0ytQEhr70nMD5MCt6KnvpQnex3Lp1C71eT/owdd1uPB5L8zU3B3pzp9Mpr9PninKgdTgclh7RSqWCSqUiXR4cvcJuEE7FY4M4Q0zWde1E1vLyMj799FO5jnw+L+QWm6N5nbxndm0rc1btFcfjMSqVivwMz0Kh8IPvE4vFkEqlJCcHgDdv3mB3d1dy7HkAJjDD4ASsjddcDCcnJ2g0GggEAhdCJBIpHJ/BXk+yjZFIRHowCQY9YUGLFHw+HxYWFizH4rFEQO+oi+iaTCLh0ev1pBxCYikej4vnHQwGqFQqqFarAM6blH0+n9RweWgQa628Rh1qavFAMBjElStX5BoIZm4QlDJO66e053/s6GHuGY/HkcvlUC6XLZ9LaSMHT+t5T6VSCd988w0ePHggXnOW80xtMw1ObQRfuVzGYDCQhcowj16NutOVlRVpT+LwKQoXOPaDeZD2HFw89AYM2whaXbvUxX8KDsh20rPycN5wOCzhKN+PulMevUcPurCwYCGsKH0jUPjemojh307nWSP4+vq6pUGgVCoBgNSDdb7J761LJJ1OR6IMTo1Ip9PI5XLCQOveWhJy4/FYjsRoNBp49eoVvvzyS5yenl4QPMw6QOcCnBo81KEyTCMjW6vVpL8xFoshkUgIrc/wUHtMexjLAr7uvNDCdD01gZ85GAxEMcRuDKp5tFrJ5XJJKEvCiGFvOp1GpVJBo9HAaDSSObD0jACkHYs/S2/IjYbMMhnWyWSCQCCAra0tYV4fP36McDgskxn0/dT13na7jVqthmq1KiwrIwwti9Tel9EKyzUUL0wmE7x8+RLv3r27EELPotDdbnMFzuFwaOloYI6mD8pNp9NYWFhAKpUSL8JF7HQ6pf+RBIYO7/hZ+jMpiePi58IiOO0eoNFoyPF36+vr8lwmk0E4HEa1WhXCieNMstmsKIai0Sji8TjC4bD0rlKEPx6PpR2MnimXyyEWi8nC183LbA7/9a9/LSetbWxsIBAICCC1uodHKRwcHODw8FCeZ0/tZDKRwdMkdfh70efLDAYD+X5HR0dS09S1zVn3msCcgBM4Z2IPDg5EMcTwlJ4wHo9LeJVKpeRxeg97SKXF3nzcLkqwLyoAki/qVjUCOBaLyfhO1h0TiYR0jvA0bobV3GAYSicSCTnVi9dOL6wnAjK0b7fblgkG1K8CENb42rVryGQyyOfzyOVyUpbRgnqG5IwgOMmBGxrP76TYgRsivW2j0UAmkxFxRygUQqFQEJac920eQEmbK3A6HA60Wi08evQI6XQa2WwWLpcLsVgMmUxG6H/gXOLXbDbR7/dFscOFrHs/7e1ROtTlomXIS8aVi5YeQ0+jY8dFoVBArVZDp9MRgb7u8+SALLa9AZCRKtw8SCQFAgFUq1WJHOLxOPx+v+R3PHiX/ZRa/ECwkyjTEwd1x47H45GyDZsN+Bp6cG5KDMtZCmK5pNfrYWlpCfV6HYVCQTayaVrmWbe5AKcGzWAwwLNnzxCNRrG4uIjt7W1ZVCx3tFotlMtlUQZROsfOEh6oywWrPwOwkiP0nAx/OS2ASiSqgXiCNoHndruRzWYRDAbleti0HYlE0Ov1JISt1+tIJpMCdCp67LXXWq2Gp0+fAgBu3rxp6dkMBAISrrPzBTivW1LcQKGDvZeS4TvDa93MTU/NEg4/lxLDyWSC1dVVmVDIPHNvbw9v3761TJ7Xv9NZt7kApz0Mbbfb+O6777C+vo5MJiPhVavVkhCsUChIVwYBRVKD2ledq2qPyc/kH11m0R6G70WSieGvLjEQ0Kenpzg+PhaPx+cTiYRlxIe+FtZQqS46OjrCyckJYrEYisWiHFdPqR+nRujTwXSzNsNU3X2jCRp+pp0QY2itxeyMKgjYVqsFr9crRNLbt2/x2WefyfmcdtHDPNhcgJNkBBU24/HZaVafffYZkskkPvjgAwyHQxGasxPE3oBNRlHncTrMtC8gna/qMJOKHOZ4rCnysB8uaOBcbQQAhUIBhUJB2sHoUVKpFICzyQ+sb1KjWq1WUalU8OrVKzSbTdy4cQPxeFymyIdCIfHEnU5Hxqfo6QT0nvrf/KPDWgAW0ouP8/uxbkovyjEmvJf0zB6PBzdv3sTXX3+Ner0O4FxQr69l1r3nzIOTC0NL+QjWg4MD/Pu//zsAWM5T4YnV1LVysdOrMK+jN9UL1k4S8W9dv5y20FmvrNVqokslUcR8d2FhQaaok3GmiHw8HiOVSsnwZpfLhU6ng3q9jmKxCIfDgQ8//BALCwvw+Xzo9/uSz/IEb06+56ngduWTLmXYQanvt/1+sLOHdV3+IVCZLrAPtdvtyqZh98zzZDMPTr2LaykZF9De3h7+4z/+A+/evcONGzeQTqeldEJ9pyZq6D21CIGmm5mBiydEc7HydQw56TE5upJHGHBj0AOsyWgyH+Z0BLLMDB+p2W21WkilUlhYWBCSi5pdaozr9bq0xPFzuanpUPZ9nspOiPFeEHz8bhTqs62Nf/PoPwI7EAgIWw1YWwDnRboHzAE4AevurkHCRffixQu0Wi1sbGwIA8raHBcMH9PezF4X1H2NOv9iSUafAULBAOca8TF9chnroQx1AVjqstTgMjwks0mR/Gg0QjKZxOLiohx+xMjB4/EIC+v3+0Wry+iAumMdYts7UXR4qx+jqIIevtPpWIBJ4okdMbynjCo6nY40V08D47S68izaXIDTXmPk37pGB8CyYMkiApAOEIawBIY9jNULmQuIAGWxnQtZE0/j8RjRaBTRaBSJROKC1I9yOD1gmQs7Go2i0WgIyFkrpeCd58Rwk6BskAILTfqwC4RRgp4TpL/r++4xyR4CkyNa6CFpzDF1fZeeutVq4dmzZ3j69KllKiLvqS7fzDIwgTkGpz0cLZfL2Nvbw/r6OpaXly3MLRcPw1sd6mkw2sMvu5fRUwLIBOsxKdTQ8nk9S8ftdov+Fjg/q4Q1Un2sA0ssdtE7/0+pH4834MYRDAaRSqWkpKPrtXYgcJNhHjocDkV1xHBbH9HHa9YdPcw1J5OJRAAnJyd48OCBlHzsJRv7ddhTllmyuQAnzf6LBs7zz1arhd/97ndIJBLSQ8lcaTAYiPhdjzVhCUSbVgXZmU6eCEaQer1eYWK1JG4ymUg9keeTMK/koGqHw2ER4wNAJBKRYxO0iodekCErWWmqjijEZ45LwAHnwNTRgZ2V5v/b7Taq1apM16NWlz/LJnE9vQ+AdKiUSiU8efIEX375pbDl+nc1zWvPIihpcwXO7/tFjsdjnJ6e4ve//z0WFxdx48YNS4dFtVqV/IjsqN/vRzKZRCgUsox2BHBhEdu9NxuRteCcUjZqdnW72GQykRYwPk9yiuHr4uKieFL92Vp0zy4RvobhsRZKUJan80ld35xm/BmGsLok5HK5EAwGLSQbr405da/Xw5s3b/DVV19hf39fNiBNsM2bzRU4tU3LpRwOB/b29vDb3/4W4/EYq6ursrgIGo4tqdfrIhIAIHmdLsYDsCxw1jwpauC/dUsZC/I0LZWjaIFdJPws3R7Gha/BxQ2DJBbfgwAhKQTgQl6pSxj2/E+/RofOBCkfZ+2VDLee5M6Q+ODgALu7u3jx4oUlOphH8QFtbsGpPZnOWzqdDu7fvw+fz4d//Md/lMZjhpY0ys90+xhBo+Vt9hqdZh91mMzpCL1eD5VKRfJLejpeJ9+HuSRDaw7DIugYWvO9CV7g/Dg//l83pev/k5DRuR7zRt1gzg2CgvVms4lisYhmsylifd4rHmLMOmetVsP+/j6++OIL/OlPfxKC7M/9zubB5hacNAJP79KDwQD3799HMplEJpNBNBq1hLa6y4MLkEclsN5on5zH96f3IjuqmV92ZxwfH6NWq6HRaGAymUj7mn4f7fUITtZMAcjx8+zq4BQH4JywolwPwIV6JuukfM/xeGzpXtHSRf4sezdZ7mm1WtKX6ff7pQ2MG9uDBw/w9ddf48mTJ+Ix9RzbaQTQPNlcg1MvdDvr12w2cf/+faysrOBnP/sZYrGYZTHrsyZHoxEajYbM2qGInSAGzg+31Q3Y9G6sC3L8ZqVSQalUQq1Wk6kMJFY4KkSznvTUlMR5PB70+30Ui0WZYACcdazoljGto9WsM0PrYrEoB9oyb9QMtS6DAJCcNhAIIJvNotls4uXLlzKoOhqNCsFVKBRw//59/O53v0M+nxfvDlwMn+cVoHMNTrsGVIdT4/HZtPivvvoKmUwGOzs7loZkPW+IrCcnGxAAulfUvhEQoLrbwu12I5FISE9no9GwtIPxRGrO2iHLySK/2+0W79jv91GpVKQNq91uy+AvbjA6Z2Y+zLov9cX03jzyYVq4zv9T3cNQPRwOY2lpSYZVO51OYZL39vbw6tUrlEolSweLJtTsGtp5A+pcgxOwllRoXLDdbhdPnz5FMplEIBDAxsYGms2mkDjtdhvA2QR4ljY4hJrlEBIl+vQxeh9NHHk8HkSjUUvHhs5DOcmAsrhWqyVTE2q1mgyvJth1WxZPrC6VSiIHpBclIPX90FPjCRy71pdA0qJ+/jsYDCKZTMo9jUajckrY8fExAoGAnPvC+U28N/r+T/udzJMZcNpKHfbyR71ex+effy7dLLlcDgBEysfHOdCKYGD+xFlEOqTVMj96HL4fAFHxTCZnh/ewoM/yBKV/uuBPoT69aSgUwsbGBuLxuIykbLVayOfzWFtbswBSX894PLY0iOtOHG42AOQ1Gtz8v9frFQ8Zj8flnpTLZVSrVRwfH+Pw8BDlcllCWLsCifdFSwTnzeYenNNMh2zj8dm81T/84Q9otVr46U9/ilwuh3g8Lo3NXLx6WkK73baUN/QRD7pDxe4d2OPIsJZ5H5liigv4OpI1fr9fBldzFlImkwEA6eU8OjoS8FCdo4HB/FnLGXX9k/eGOSsljfo0a4oWWH/1+XxoNBoSwnc6Hbx9+xbPnz+XubX6+9vlefMISpoBp82mlVgAoFwu449//COOj49x48YNXLlyBSsrKzJcq9vtSg8lPUm73UYgEJAuEHoDLnq9GLUyiJ/LsZ0cT6mF6BoQ/HnWS/kZ3BRCoRAmk/NBz/bjAe3RAkURPp/PMr6SwNTnbnJDIrDZFsbr1HJDtt4dHx/j3bt38jp9zLyxczPgfI/pXZtAarfbePHiBd69e4d0Oo0bN25gZWUFV69exY0bN6QflGepBINByfsIBoJIl2+0wJzjQfRk83q9btHVEiT8WwslarWagJfC+GAwKK1mGqTAeTsWPaEWLDAisH8u65R8jHOHOHOpXq/L4w7H+ZR5njH67NkzOS/FTvhwozFmwPm9ZhcPjMdjEas3m028e/cOoVAId+7cQa1Ww927dy3jOAgEXdukyJtgo+ehl9PqHx7HHolEJL9kuYThLdnbXC6HpaUluFwuOa5eH/NAoDKk1nVNbiAs92jJnj0f5cbDA5HoMTmFvt/vS08qBRI8V/Tw8BD37t3D8fGxXIduuTOgtJoB5/eYXUBA4w7P3Ourr74S9vTGjRsIh8PCiuqck6QMyRxOHeBkBQKGEjjWNcmuptNpASePgm80GgJWhtBOp1M8Ics59rM5df7Ls2DG47FMXddjRdghQw9qJ4u4+ZBpHo3ODlQ6Pj6W75vP5/H69Wvs7++LPFGXsoxdNAPOH2h6EWl1zmRydmzd7u6u5GDXrl2zHCCkj21gLbTZbMpJZwQic0CdNzL348wgbgg854Uza/VBuZFIxCIM0O9r18ZSSVQuly3HLNBjaqmdFqFPY5+73a7ohXmUwv7+Pg4ODlAsFtFut8XD2jc8YxfNgPMH2jQBti6Y93o9vHjxArlcDrlcbmpNUAvnO52ORZPK/Iyhrw4vufh16UKHwpzwTq8KQPSsHHPJ9wesE/Cr1Sry+TxOT0+RzWYthxUBkM2EOTQ9MN+HYOPGMRqNZOr77u4uXr16hVqtJjmqVkfxPYxNNwPOH2h2Wt/O6tKDPnz4EJlMBtlsFslkUk4U06MvCVTOq2UrGgDL1HS+jsOxKHGjjpc9nhsbG1LnZPjMRvHRaCTH//H9GfI2m03k83k8e/ZMRpSQUeUmwXyUXlGLKLhp6JJOvV7H3t4eHj16hO+++w71et3SYWK/lwac7zcDzj9j02pt75OTUfL3X//1XxiPz07LunXrlqhxOOG81Wqh3W7LtAIAkq8ROKxXjkYjEcRThDAej1Gr1XB8fIzxeCz5LScCspyh+z4paCBZ5HCcTb/f39/H4eGhbCCRSES+D+f/sOZKJRT/eDweGUlSLpfljJTPP/8cDx8+FAWVvmfzrvr5n5gB548w7QHsRFGv18Pr169xenqKR48e4fbt27h9+zauX78Oj8eDQqGAZrOJcDgsZ30yt+MUPAoKGD7qonyz2cTr16+xu7uLo6MjbG5u4pe//KXoe8fjsehhKcJ3uVxoNps4PT0VkT1rszwSkfK6g4MDKelQW1utVuFyueTAI5Z+6vU6yuUyXr9+jYcPH+LBgwci2KfumOEzrx+ApbZr7P1mwPkjTHtOew8k2dBarYZHjx7h9evXuHfvHnZ2dpBOp6U/c3t7W8JBljM4npLlCQCia2W9kEcqfPbZZygWizI9IRAIiKa3VqvB4/EgEolYQFgqlUQoz41kMpnIqJTRaITT01OZZ0RpII/kq9fr8p7Hx8fY39/H/v4+Xr58ia+//hrlctlSs7WXSeztbsa+3ww4/0KzLzrt5XgkXrPZxP7+vgzwunbtmkXeRrEC8zYSPXwfeqBWq4W9vT08ffoUhUJB8tl2u41UKiX1RwJOTzegJyWpxBIJr0N/Jt9H57jdbleOqOj1enj37h3u37+P3d1d1Ot1DAYDS28qcB5h2AFpwPnDzIDzLzCtSdX9mYC1NkrQ6WPii8Ui3rx5g5/85Cf4+OOP4XA4pGZJgDGfo2jhxYsXMpmOJ1MvLy/LhkB2lu1dAEQswP+TRCLRpL8HiSB+J0r0OHj67du3ePXqFRqNBnZ3d7G/vy/vob8/yymso+oNy9gPNwPO/yWzk0T2XkfmXuwo2d3dxevXr/HVV1/hyZMnWFtbk2MY2DDtcDjQ7XZRrVbx5MkTPH78GA8fPpSwlGyvHt5MAQTZXAoZWC/laBOyxTzyQLPDJI5IKFF1xB5MHhtI8kp/T+B8VtG0kZrGfrgZcP4v2LSZN3ZvoYkjEkD9fh/1eh1HR0dYWFhAt9tFKpXCL37xC1SrVWFemWMeHR1Jhwe9WT6fF7aUDdUEL4X0bCWjRwPOc1mK1enhBoMBSqUSSqUSer0eqtUqCoUCHj9+jJcvX6LZbAro7KEq67NGI/u/Ywacf6FNW4CazdWv0wV4rbWlDM/hcCCfz6PZbKLb7SKZTKJWq+HBgwd4+fKllEAAyCzaSqWCfD6ParWKzc1NGTPJGqXD4ZBZtFr9oyfat9ttlMtlHB0d4eDgAK9fvxbxAMevcPj197VxaTGCsb/cHN+3uzkcDrP1/Uh7HytpZyztITD/HYlEkEgkhCnVoyRpTqcTmUwGGxsbWFlZwbVr1xAOh0XgPpmcTSVYWloSRpZ5bz6fl1IJZxbt7+/j+PgY1WrV0vYFwHLOC7+H/XqM/TibTCZTC74GnP+H9kOK7tP0r/a8jo/r+UOxWAzJZBKj0UjOXQkEAojFYnKQEcPocrmMd+/eoVgsColEccI0ANLDT7ue931HA9ofbgacl9SmLfjvEzkQoLq7hLme1ryyGZvyQI/HI8fw6U4TPa6EpkksA7a/vhlw/h2ZBoTOUe2v0TVFbfqQWw3caSUN/Zjdsxtw/t+YAecMmwaffUC2JqLsYer7fvfvC78NQP869j5wGrb279ymNYLbn7O3uL3vuAPaNC9tgPl/bwacf8dmDzn/nETufd00f84MMP82ZsD5d2z2cox+bNrrjP192UWqztjfpf2YHNLY5TYDzr8j+7EgM57z79MMOP+OzIBsvsyA05ixS2oGnMaMXVIz4DRm7JKaAacxY5fUDDiNGbukZsBpzNglNQNOY8YuqRlwGjN2Sc2A05ixS2oGnMaMXVIz4DRm7JKaAacxY5fUDDiNGbukZsBpzNglNQNOY8YuqRlwGjN2Sc2A05ixS2oGnMaMXVIz4DRm7JKaAacxY5fUDDiNGbukZsBpzNglNQNOY8YuqRlwGjN2Sc2A05ixS2oGnMaMXVIz4DRm7JKaAacxY5fUvvfYeWPGjP3tzHhOY8YuqRlwGjN2Sc2A05ixS2oGnMaMXVIz4DRm7JKaAacxY5fU/h+rXn8CsZhtUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ants.plot(reg1['warpedmovout'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}