{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image registration with nipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nipy has some flexible routines for image registration.  You can choose many parameters such as cost function for matching the images, interpolation method for resampling the images during registration, and optimation method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# - make numpy print numbers to 4dp only\n",
    "np.set_printoptions(precision=4, suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "# - set gray colormap and nearest neighbor interpolation by default\n",
    "plt.rcParams['image.cmap'] = 'gray'\n",
    "plt.rcParams['image.interpolation'] = 'nearest'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we do a registration rather like the default coregistration routine with SPM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import nipy\n",
    "import nipy.core.api as nica\n",
    "from nipy.algorithms.registration import HistogramRegistration, resample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we make a new 3D image by taking the mean of all voxels across time in the 4D image.  We drop the first four volumes because of the T1 effect (see the diagnostics notebook):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img = nipy.load_image('ds114_sub009_t2r1.nii')\n",
    "mean_data = img.get_data()[..., 4:].mean(axis=-1)\n",
    "# Drop the time dimension from the coordinate map\n",
    "coordmap_3d = nica.drop_io_dim(img.coordmap, 't')\n",
    "mean_img = nica.Image(mean_data, coordmap_3d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"source\" image is the image we are going to move to match the \"target\" image.  Here we are moving the structural image to match the mean functional image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "source = nipy.load_image('ds114_sub009_highres_brain.nii')\n",
    "target = mean_img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AffineTransform(\n",
       "   function_domain=CoordinateSystem(coord_names=('i', 'j', 'k'), name='voxels', coord_dtype=float64),\n",
       "   function_range=CoordinateSystem(coord_names=('scanner-x=L->R', 'scanner-y=P->A', 'scanner-z=I->S'), name='scanner', coord_dtype=float64),\n",
       "   affine=array([[  -1.    ,    0.    ,    0.    ,  122.7012],\n",
       "                 [   0.    ,    1.3002,   -0.    ,  -86.6199],\n",
       "                 [  -0.    ,    0.    ,    1.    , -173.3327],\n",
       "                 [   0.    ,    0.    ,    0.    ,    1.    ]])\n",
       ")"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "source.coordmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These images are actually fairly well matched before registration, because they were taken from the same scanning session, and the subject didn't move much between the functional and structural scan.\n",
    "\n",
    "The affine array for each image records the positions of the voxels relative to the center of the scanner.\n",
    "\n",
    "The affines show us the correspondence of the voxels, even though the images have different fields of view and voxel sizes.   To show this we can use `resample` to resample the structural image to the same array grid as the functional image.  The `resample` command works out which voxel value to put in the empty voxel grid by using the affine and then resampling from the structural image array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Resample structural to functional without transform\n",
    "source_before = resample(source, np.eye(4), reference=target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x11219ca90>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAEdCAYAAAAcgiN1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuwXFWd9vFn5Rgn3CQmQAIJlyAEAgQICEEuQiDhooDg\nWCAzOhSjjIU4r/DqKEqVxdSMglAlyOCoIw6TYd6JOoi8REQIAQTkIkgQSAjhFu45CQEjMkYzyXr/\nSHPe/j1J9upOn8s+J99PVYrzS3fvXr2798qm93N+K+WcBQAAgPUbNtADAAAAqDNOlgAAACpwsgQA\nAFCBkyUAAIAKnCwBAABU4GQJAACgQkcnSyml41NKC1NKT6WUvthbgwKAvsb8BaBVaWP7LKWUuiQ9\nKWm6pJclPSjpjJzzE703PADofcxfANrRyTdLB0t6Oue8OOe8StIPJH2od4YFAH2K+QtAy97RwWPH\nSXqxqX5J0tTmO6SUaA8ObGJyzmmgx9AC5i8A67W+OayTk6W2J5KU4vMPGxa/2BqIpVfWrFmzzjjq\nhPFtvLqO7e3j4O3x1XXJobfHt2bNmoEeSl9oaaePGDFCq1at0vDhw/t6PButzuOr89gkxteJOo9N\n2vjxrVy5cr1/38nJ0suSdmyqd9Ta/zvboJzzOidMADAAWpq/Vq1apTVr1mjVqlUaNmyYurq6+m2A\nAPre6tWrW/ofwk7+t/shSbunlHZJKb1T0umSbqx6ACdKwOCTUgp/hoiW5q/hw4dr2LBhGj58OCdK\nwBDU1dWl4cOH9/zZkI3+Zinn/D8ppc9IukVSl6TvD8bfJKn75M/4Nl6dxyYxvoHUzvxVx0u5zeo8\nvjqPTWJ8najz2KTeH99Gtw5oaeMWkKxjZgnob34c1P1z7+MrjXeQBLyLUkp5xIgRAz0MAP1o5cqV\n653D6n1qCAAAMMA6CXgPeoPt//AxNNXtc1i38QDAQOObJQAAgAqcLAEAAFTgZAkAAKBCrTJLvf2r\nyp61KG2/3axG1fbafe7eNtDPX9LJvm1lewP5+E73denxvf3ekkkCgGp8swQAAFCBkyUAAIAKnCwB\nAABUqFVmqbf1dXakL5+7UwP9/CW9Pb6BfK974/ED+Vx1/6wAwEDjmyUAAIAKnCwBAABU4GQJAACg\nwpDOLAF1xNprADC48M0SAABABU6WAAAAKnCyBAAAUIGTJQAAgAqcLAEAAFTgZAkAAKACJ0sAAAAV\n6LMEoBJ9oQBs6vhmCQAAoAInSwAAABU4WQIAAKhAZgnoZ2R+AGBw4ZslAACACpwsAQAAVOBkCQAA\noEK/Zpa8XwsAAEDd8c0SAABAheLJUkrpX1NK3Smlx5r+blRKaU5KaVFK6daU0si+HSYAtI/5C0Bv\naOWbpWskHW9/d4GkOTnniZLmNmoAqBvmLwAdK54s5ZzvlvSG/fXJkmY2fp4p6ZReHhcAdIz5C0Bv\n2NjM0picc3fj525JY3ppPADQ15i/ALSl44B3XtuOmJbEAAYd5i8ArdjY1gHdKaWxOeclKaXtJS1t\n5UE5Z9oHABhoLc9fq1at6vl52LBh6urq6o/xAegnq1ev1po1a4r329hvlm6UdGbj5zMl3dDKgzhR\n6j0ppfAHQMtanr+GDx/e84cTJWDo6erqCsf5hqTSop4ppVmSjpS0jdZe3/+KpP8r6UeSdpK0WNJp\nOeffruexYePDhg3z21t4KVgf33cszoq+4p+t0v+F5Zxrc2B3On+NGDGi/wYLYMCtXLlyvXNY8WSp\nE5ws9R1OltBfBvPJUic4WQI2PRs6WaKDNwAAQIV+XRsOg4d/C9juN1d80wUAGCr4ZgkAAKACJ0sA\nAAAVOFkCAACoQGappjzz4xmid7wjvnWrV68Odek3Df03mvz5SpmjUqapVJd+M7Ld8Qzm3w6s+9j5\nrVUAmzq+WQIAAKjAyRIAAEAFTpYAAAAqkFnqJ6XMjvM1av7nf/6ncnvjxo0L9SuvvBJqzzg1LxAq\nrZtLKWWk2t1eKYfjGSrfnt/ebh8ov50cDjC01D37h8GNb5YAAAAqcLIEAABQgZMlAACACgOaWRpq\n15ircjCe8Sn1Gdpiiy1C7Rme3/3ud6H+wx/+EOrNN9+8cnueifLtjx49OtR/+tOfQr1y5cpQv/nm\nm5W194F617veFequrq7K2p/Pt+efHc84+fhdf2aa2t32UDtOgPXxOXH8+PGh3mqrrUL9Z3/2Z6H2\nOdaP+V//+teh9jnive99b8/PK1asCLc9+eSTlc+NoY9vlgAAACpwsgQAAFCBkyUAAIAKqS/zDyml\nsHHPoQz2LIaP953vfOcGbyu9Ns/g+PX5KVOmhHrPPfcM9c477xzqbbbZJtSbbbZZqP36fnd3d6jH\njBkT6rfeeivUI0aMCPUzzzwT6tdeey3UnrF67LHHQu0Zp+XLl1c+vtSnyvMI/lnz20tr7Xmfq04y\nTaXHtnscdJqvavfxnm9zOech0cQqpZT9c47W+Xzvc4ofw5MnTw71c889F+qXX3451NOmTQv1qaee\nGupZs2aF2ueU97znPaH+7W9/2/PzkUceGW777ne/G+qddtop1J4ZffHFF0O9dOlSYXBYuXLleucw\nvlkCAACowMkSAABABU6WAAAAKvRrnyXPRvR2L5v+zjx5r43m1+OZF6/9tc+YMSPUp512Wqh33XXX\nUHumqbRWmj+/9xHZcccdQ+0ZHc/4eB6hOa8lrZtr8e0dfPDBofbMkj/+jTfeCLWP39fC+9WvfhXq\nZcuWqR2+v0oZqU4+e/2dOSo9frBnCdE3/HPhvdp22GGHUHumaI899gj1woULQ/2BD3wg1J45Gjly\nZKjPP//8UF9//fWh3nbbbUM9YcKEUJ9zzjmhfv7553t+3m677cJtjz/+eKjPPvvsUHvm0+ezq666\nKtSLFy8OtWeevPZMJfof3ywBAABU4GQJAACgAidLAAAAFWq1NlwpK1HKTvR2/xq/v6+v5jme5lyO\nZ268L9InP/nJUB966KGh9v4uf/zjH0PtGSHvIVJa+620ltx///d/Vz5/aa04fz7POG299dahHjdu\nXKh///vfh3q33XarfD5/L/bdd99QL1myJNT33ntvqD1DUMpseYag+f7tZvPaXZeu3YxSu9snowRp\n3Zze9ttvH2rv5ebHnGeWXn311VD758wff/vtt4f6Ix/5SKifeuqpUPuc6etPeu83n0OPOOKInp+/\n8pWvhNsOP/zwUHvfO39tvg6dz1fnnXdeqH2+bu75JEk33XRTqJ9++ulQ+3xbyrCifXyzBAAAUIGT\nJQAAgAqcLAEAAFTo18ySXzf1HEhJp1mLdh/v15H98Z5rab5uvNdee4XbLrroolB7zxHvy+HXrL3v\nj+87X7vNM0KeQfLXVrrdM0t+f983vpabZ6S23HLLUHvfJN+e93Dx1+uZgLFjx4Z6/PjxoT7ggANC\n7X2aPMM0e/bsUHumqrnnlr83XnueoN2+Ru1m+UqPJ8MEad3PwcSJE0P9/ve/P9RTp04N9T777BNq\n73Xm60H6Me+91Lz32wknnBDq5oyRJP31X/91qH19Nl87ztd3a85hnnjiieE2n0+cZxh9rUxfh863\n77d7Xyefjzw/29wjSlq3J5XnuebPnx9qjvGy4jdLKaUdU0p3pJTmp5QeTyn9r8bfj0opzUkpLUop\n3ZpSGlnaFgD0J+YvAL2hlctwqySdn3PeW9Ihks5NKU2SdIGkOTnniZLmNmoAqBPmLwAdK54s5ZyX\n5Jwfafz8e0lPSBon6WRJMxt3mynplL4aJABsDOYvAL2hrcxSSmkXSVMkPSBpTM757QVxuiWN6dWR\nrUdfX1fdbLPNQu19ODwX4zmi5vXbvvSlL4Xb/Pq/Z3q89r5HPhbPS3kOxjNPnpvx7fvzeybH+yr5\nunj+3njfj9L9ffujRo0KtWeUvK+Uvxe+fb+G7z1i/Jq+Z5zGjIkf7x/84Aehbl4byvNeXvt70WkG\nqVObSl5hoOevuvGMkucofS21F154IdRHHnlkqD1jNGfOnFD7WnCTJ0+u3P655567vmH38GPaM0hn\nnXVWqH2O9Jziu9/97p6f3/e+94XbqnrqSdLo0aND7Rkj31e+b++7775Q+3zmfeI++MEPhvrGG28M\n9XHHHRdqn999e/76PD+2qcwRVVr+bbiU0paSfizpsznnkEbOa/ckexNALTF/AehES98spZSGa+1E\nc23O+YbGX3enlMbmnJeklLaXtLS0nZxzxyukA0A7Opm/mr8hHTZs2Dq/ZQpgcFu9enVLv5nfym/D\nJUnfl7Qg53xF0003Sjqz8fOZkm7wx65nW8UBAUBv6XT+Gj58eM8fTpSAoaerqysc5xvSyjdLh0n6\nmKRHU0rzGn/3JUmXSPpRSukTkhZLOq2zIbev1C/Gb/dr1qV+N57b8e01X+OWpM997nM9P/s6R57h\nKWWQfB0j70vkfY98rN7XyF+bj8drz2e5Ug7HM0qeSfLa/yHy2/0aful2f698f/vjne+P/fbbL9S+\nv++6666en73Hiu8bfy/8c+ivZaAN8v/Jqe38NdB23nnnUB944IGh9kzRggULQu3HrH+OPXPkucPd\nd9891J6zKWX5PvShD4Xae6f559aPKz9ODznkkJ6f2+0B+POf/zzUzfOBJP3FX/xFqD0TdM8994Ta\ne0j52H1tTc8g+fx9wQXxlz232267UF922WWh9kyn/3u0KWaYiidLOed7tOFvoKb37nAAoPcwfwHo\nDSx3AgAAUIGTJQAAgAr9ujbcQPNr6ltttVWofX025zmdv/3bvw31oYce2vJjS5kgz9T4/T2T5LWv\n3eavfenS+Ms/pcyRZ6p8PJ4P8IyVX6P3TFApw+TPX8pkeUaplEHwTIA/3p/P8xHN769n2Twb4fvC\n3xt/7aW+TEAr/Jj02nOW999/f6i9j9Ill1wSaj8GfS25ZcuWhfr0008PdWntUD9OPvnJT1befuut\nt4b62GOPDbXP181zUFXQV1p3LTjPEF144YWhnjJlSqgXLVoUaj/mfV96D8Dmnn7Suj2yPGPp6/T5\n/H/FFVeE+rnnngv1t7/97VB7H79NAd8sAQAAVOBkCQAAoAInSwAAABVqlVlqt29Su32VPIfivXJK\n2RDv6+F9QV5//fWen0t9i7bYYgtV8bF4n4tShshv92vwvlaQ5xd83/l4SpmgUq+gUg8V78nimabS\ne+/7w9/r0vj9+Xz/+f4ZO3Zsz8/en8b3hffI8tv9s1Pq0wS0wnMuM2bMCHXz/CVJ559/fqj32Wef\nUM+bNy/U11xzTainT4+dGa688spQe3bPeS6odBysWLEi1Pvvv3+ofU6YPXt2qJv3j+e3fCx+/B92\n2GGq4vOTZ5y851XzfCKtu7aoZ5JmzZpVOV7fN57B9HzX1KlTQ/3xj3881P5ebgpzEt8sAQAAVOBk\nCQAAoAInSwAAABVSX/ZsSSmFjZcWoizlUNq9v+dM/PbSemKeK7r88stD7de1ly9f3vOzXxP2nIrn\nUnzflPoylXIwfn/fXikD5Nfk/fEjRowItb9efz0+Pn9+7yPimaXSOn2+P328vv13vetdoS71jfLX\n5/unOQPmY3vooYdC/fTTT4f6pZdeUpV2j4u+5nk4l3MeEgGGlFL2z/lg4hmfY445JtR77rlnqD13\n47lG7+Xjvcb8+bz2HI1rd74vPd7v78/vc2hzXRq7b9v72vl8U+LzY+lzV3qtTz31VKh97bgddtgh\n1N4Da+HChaG++eabQ33HHXeEurTW5mCycuXK9c5hfLMEAABQgZMlAACACpwsAQAAVBjQPkula9Cd\nZjVKaw2Veu0cddRRoR4/fnyoPbvha4I180yNXwP32/0atmdkSmuZlTJPvm/9tfjz+e2eKSr1uCr1\nvPK16KryBFL5s+Hb99tLuZtS3s1vb85A+ba9P80222wTar/+7+sulV67G+hMEwaGH2N/8zd/E+p7\n7rkn1N437jvf+U6oPaPkfZlKOUevSzrt1VN6fOmYb84O7rLLLuG2Ug++djNKzjNKnf5bt9tuu4Xa\nM59///d/H2r/98czTwcddFCo//zP/zzU//Vf/xXq0r4ejPhmCQAAoAInSwAAABU4WQIAAKjQr5ml\n3s5SlHIyfrv3/vHr0J5JOuOMM0K97bbbhrq7uzvUzb14vC9PqQdJKdPj1+NL1+c9A+Tbf/PNN0Pt\nuRgfnyvtW1fKXPk1dVfKELW7jqDz/e2ZMu+5NW7cuA0+3tfY8jyA7/sJEyaE+v77768cW+m9wabJ\ncyrPPPNMqH39sUmTJoXae+mUcox+zP70pz8NtWfvzj777FB32jeppN05YcyYMT0/l157b+vrnKFn\nqk444YRQz58/P9RHH310qE866aRQz507N9Q+P1977bWhHgo5Sr5ZAgAAqMDJEgAAQAVOlgAAACr0\na2ap0z4a7fJsh+dqfD2bI444ItR+jX/FihWh9t5FzTklf61bbbVVqEt9gLyPka895Jkf532anO+L\nUuaotBadX9P3zI+/Pq99e54ZKPU58UzR5ptvriqlzJXX/nr8s9Ncb7nlluG29773vaH2Nbc8L+CP\n92xIaY0t1+kai/193KI1nkv0DNKJJ54Yaj9G/X31+cyPyfvuuy/Uzz77bKhvuOGGUHvms5S79M9d\n81qb0rr9yfwYvOWWW0J9++23h9r7RHlfvOZsoc8fvu98bbWtt9461J7h8fnM14McPXp05dg65T2z\ndt9991D76/P51j8Lu+66a6gff/zxyscPhTmEb5YAAAAqcLIEAABQgZMlAACACgO6NlynSte8S+ub\n+f1HjhwZas8NldZ3a17fx69Ze6bGM0i+bc+plPJX/nivPXNUWrvHX5vXfg27dE3a92XpvfPn8/fS\ncz1jx44NtWfCvPeR739//lLfKX98c37Ce5r452r//fcP9aOPPhrqBQsWhLqU/+rUUOiBsinyzJKv\nZ3b11VeH+gtf+EKoPUf5wAMPhNrXB5s5c2aom9dDlKSPfexjofbM0quvvhrqSy+9NNR+jF1++eWq\n4q/fX89tt90W6ptuuinU3guveU7xY3Ty5Mmh9szSdtttF2qfHzzT5D37vMdfu3w+L/Vma/ffzrvv\nvjvUV1xxRai971yp5+FgxDdLAAAAFThZAgAAqMDJEgAAQIVaZZba7fdSysmUMkujRo0K9QsvvBBq\nz5p4juWNN94IdfNaSKX8k9eeUfKeJ35N3F9LKbPk+9Kv13sGyq95l9aqazcT5X1M/Pl9vJ5H8D4k\nnp/wdalKfao8F+Tb99fvma3mvJr3f3nllVdC7f1iPA9x7733hto/x/7c/nyu3eNqKOQLNgWei/TP\n6D777BPqiy++ONR+TPnn6Iknngj1hRdeGOrp06eH2nM3Pgdcf/31ofY5yj/XPif6Me6OOuqoUPtx\ndccdd4Taj6vm+dzve9ddd4Xa93UpMzRnzpxQv+c971GV0jHovdp+8pOfhPqXv/xlqK+88spQ33nn\nnaH2927PPfcMtfeK8387S30ES33/BoPKb5ZSSiNSSg+klB5JKS1IKV3c+PtRKaU5KaVFKaVbU0oj\nq7YDAAOBOQxAb6g8Wco5r5Q0Lee8v6R9JU1LKR0u6QJJc3LOEyXNbdQAUCvMYQB6QzGzlHN++/rF\nOyV1SXpD0smS3v490pmSTumT0QFAh5jDAHSqmFlKKQ2T9LCk90j6ds55fkppTM65u3GXbklj+mJw\n7a4nU1qPp9RXyXtJHHDAAaH267auuXeGZ3A8I+O1v1a/Pu99hTzzVFrbzTM53gvIM1Ke8XG+r0sZ\nJM9DlPpw+Hh9e55JWrp0aaj9vff97eP3a+4TJkwIte9/z7c191Hx1+o9tvx6/kEHHRTq+++/P9S/\n/vWvQ+3bb/c4Kd1/KKzj1Gwg57C+5J9ZzzX62m1TpkwJta+F6cfMRRddFOo333wz1Jdcckmo/Rg+\n+eSTQ+2fa88seQbqtddeC7X3SvNczumnnx7qa665JtTHHntsqF9++eVQN88Rvi/8mPBj0F/LCSec\nEGqfT7zvkc+H3pNq4cKFofY8lq9j6hmmm2++OdTec8r7JJ111lmh9vn1U5/6VKh/8YtfhNpf31BQ\nPFnKOa+RtH9KaWtJt6SUptntOaVEIhRALTGHAehUy60Dcs4rJN0k6UBJ3SmlsZKUUtpe0tKqxzZt\nY2PGCAAd25g5bNWqVT1//NtJAIPf6tWrw3G+IaXfhtvm7d8SSSltJmmGpHmSbpR0ZuNuZ0q6oZVB\nDbWv9wHUW6dz2PDhw3v+lC51Axh8urq6wnG+IaXLcNtLmtm45j9M0rU557kppXmSfpRS+oSkxZJO\n641Bt/vNU+n+pevCfk3ct+drdu2xxx6h9j5LzdkUP0P1DJLnVpyPZcWKFZW3l9ZOc6U+TaVMledw\n/PbSvvf7V31IpXX7injGyp/P929zHyRJ2n777UPteQ9fd8rXcvLXM3r06A2Ozfsq+e1PPvlkqP0f\nZc9q+GeptM5dpwb5/+T06xzWnzxzc+utt4Z60qRJoZ41a1aor7vuulB//OMfD/VJJ50U6u7u7lD7\nGob+ufNc4UMPPRTqj3zkI6H2HOA//uM/hnqvvfYK9fLly0O99957h9ozWj/+8Y9D/aEPfSjUzbkc\nny+89rnf+xKdd955ofb5wucX7ynlGSPPe3m+drfddgu190HyvlG77rprqJ9//vlQP/7446FesmRJ\nqG+//fZQ+5xW9Q3NYFV5spRzfkzSAev5+9clTV/3EQBQH8xhAHoDy50AAABU4GQJAACgQq3WhnOl\nNatKa7+1+/jLLrss1N5nybMiixcvDnVz743S2mulPkA+Vr9mXrrd1+Lx8fjz+b7wzI+P3/sglR7v\n74Xfv5Rp8u35NfLS2nSeifL7e8bL+0Z5LsjH2zwef27vyeLvzbx580Lt2Q5/Ln8vXacZI35rdXDw\nnJ2vBfdP//RPlY/3z4n3bfLMkmc8S/PjU089FWrPEPnn2nOHnqvxDJRn+XxO8uPOj/mpU6eGurmf\nmWd0vMeUr6X24Q9/ONTet2jZsmWhHjt2bKg9s+T7zjNO3kPLM0m+tqXf7plTz3B+/etfD7XnyTzv\nNXfuXA11fLMEAABQgZMlAACACrW6DFe6bOZKl9X80pF/zeu/7ukt6f2rT78U419dNv86qb8W/9X3\nUusAH6tfFvJLPf6rrH6775vSZTR/rf61rY/PX6//6qi/F6XlTnz8Ph7/St7fK99fpct6fn9/fT4+\nf3zzV/5+OcMvGfqvYPtSBf4r0L6sQ+lyQ+kynL+WQd4aYJPlx8gjjzwSar+0vO+++4baf3X/m9/8\nZqj9UpMvJ+L8mPDP+de+9rVQ+2U7X2Jjxx13DLVfSvJLV95awI/h888/P9QPPvhgqP3yeDM/Zrxt\nwne+851QH3rooaH2S5Y+P/rSW5/4xCdC7ZcQr7766lB7q4ATTzwx1Ndff32o58+fH+oZM2aE2ttS\n+Jzmy5/4Z+WVV17RUMM3SwAAABU4WQIAAKjAyRIAAECFWmWWOv2V5XazF35N26/L+nVi/5Vtv+be\n/HjP6Pi2PHPjuRbPpfhz+Wst/Sq9v1av/VfrN9tss1D7vim1bfDXX1pOpZQ3K2WuShknf31+u7+3\nvn3nv7bdvD3PUjQvhbK+x06bNi3U/l42/0qztO5xUjpuyCQNTX6MePuQ6dNjg/LTTosruixatCjU\nvqTS97///VB7pumrX/1qqI855phQ+xzjx4X/+rxnlJwfJ85f/+GHHx7q5557rvLxzXOUZxhL+Vef\nT77yla+E+lvf+laovS2Cz1++fMlbb70Vat9Xp556aqgvvfTSUHt+bL/99gv1KaecEupvfOMbofY8\nl3+2/PmGIr5ZAgAAqMDJEgAAQAVOlgAAACoMaGap3SxFKatR6tNUuob+4osvhnqnnXYKtedsPNfS\nnDXxjIxnjkrLnTjv0+Svxfnze2bIM0DeR8kzS54H8GvofrvzHE7pvfPMgLfv9/F65siXJ/Br7n67\nt/P3/ec9uDzDsHz58p6ffRkGzwt4/xvf1w888ECo/XPq9/f30pWOMzJNg5Mvt+S9drxv0aRJk0Lt\nucmjjz461B/96EdD7b3M/HPpc9I222wT6i9/+cuh3mGHHdQJ/9z7kiN+TPtx532hmnOZpUyo5w69\n15nnwbwP0/HHHx/qI444ItT+/D7fnXfeeaF+/vnnQ+3zs382/L2dOHFiqK+66qpQP/3005X1prBE\nEt8sAQAAVOBkCQAAoAInSwAAABUGNLPUbn+Ydte8KuVinn322VAvW7Ys1Pvss0+ofb0cX/+mORfj\nGaHSWm2lPkqltdZ8bSG/hu7397yVZ448n1DqK+KZqlJfEr/dn8/zBL72kGeI/L3wnjGeWfLPhucv\nPDPlmQXPizR/Nvy1+Xvp/WV833l2zvlnoaTdNRYxOFStTyhJ3/ve90Ltx4SvV3bhhReG2tcsdO1+\nbjwD2il/vd47yI8zzzX6GozNuUjft74Wmq+z5/kvn298PvFMpG/fM1E+h/h86n2XZs6cGepSJtZ5\nptPXlvN/PzaFOYRvlgAAACpwsgQAAFCBkyUAAIAKtVobzrWbaSpllDwn5DkUz474dWPPongvi+bc\njD+XXzP2sZZ651T1dFrf2Pyau+8Lv92vgfv1fb/dt+fX6D2DVOor5bVnojxT5Jklfz5fj8235/uv\nND7Pt3nmoDkz5p+bMWPGqMrtt98eas/OdardfmT0XRqcfP3Ggw46KNTe++fBBx8Mtfcp8pzgnDlz\nQu1z1IwZM1ofbC/wOdF7D3kOs5RpbX49pbU1S2tHVm17fdsrZTzbVcoo+fZLfef89v/8z/8M9Rtv\nvBFqfz1DAd8sAQAAVOBkCQAAoAInSwAAABVqvTZcu/1h/Dqt53L8urHnXO66665QH3fccZXP7zmY\n5mvkvpaZ39ev8fo14dI1ce+1U8qhuNI1+dL1fc9YeV3a986fv5Tx8l5Entfwvibek8XH558F59fg\nPYe06667bvC5mteNk6Qnnngi1D/4wQ9C/cwzz4S60wxRu2vDkWEanHxO8NzeSSedFGrvFeRz0Gc+\n85lQe0bQu2f8AAAVEElEQVTz9NNPrxxPf3+O/Bj2OcT7TPnrbe5V5PNVaf5y7R5TnWaUPHP52muv\nhXrBggWh9h5Ul112Wai33Xbbyu37+IdiRsnxzRIAAEAFTpYAAAAqcLIEAABQodZ9ltrNNJWuA/s1\nfb/9nnvuCfVjjz0W6kmTJoXae/80r+9Wun7u/LV47qW0HlhpX/jjPRPl1+RL16T9+Uq3l/pG+Xvh\nfYx8e35/r33/lV6fZ5i23nrryu37Wn7Nr8e3/eqrr4b6l7/8Zajnzp1bObb+zgOQURqcHn744VD7\nemV+DHhu0nMqS5YsCfWhhx4a6mnTplWOx3MznvPrNKdT4q/Xc6Red6LdzGhvu+qqq0L90EMPhdoz\npVtuuWWoJ0+eHOo777wz1E899VSoX3/99VCTWQIAANjEtXSylFLqSinNSynNbtSjUkpzUkqLUkq3\nppRGlrYBAAOFOQxAJ1r9ZumzkhZIevu7xQskzck5T5Q0t1EDQF0xhwHYaMXMUkppvKQPSPqqpP/d\n+OuTJR3Z+HmmpDvVwmTT7hpV6xlL6SmCUk7It+drw82bNy/UU6dODbVf923OvXhfCs/geGbH1zHy\nzE27mSZ//lJmp5QJ8vH6NWqvfTy+fd+ej8dv9/fS16Lz8W6xxRahLmWg/Pmb13qT1t2f/vqac0n+\nXN4T6rrrrgu193vxvFupB1ZvG2p9lnpzDqszP2bGjRsX6oULF4b62muvDbX3ZfJjznN8nlvxTJKv\nkejba7cXUZ0+l5738nX3mns2SevmwXxtym222SbUPid4/sz37TnnnBPqm266KdS33XZbqP3fLr+/\n53G9T+CmkFFyrXyzdLmkv5PUPGOPyTl3N37ullS9UigADBzmMAAdqTxZSimdKGlpznmepPWexue1\np/v9G/0HgBYwhwHoDaXLcIdKOjml9AFJIyS9K6V0raTulNLYnPOSlNL2kpa28mQ550H/lT6AQaWj\nOaz5cuuwYcM2ycsPwFC2evXqltpYpFbzDymlIyV9Pud8UkrpUknLc85fTyldIGlkznmd6/0ppbDx\n0no7vZ3FKD1fKWfj15m/9rWvhfqYY44JdXNmwLftPUcWL14c6tJaZp6n8mvYpT5M/to9I1Vaa87z\nEN63wzNEfo279Py+fc8Ilfoq+WfHr8n7e+nP78/nfZRGjRoVat//zXkM70/zH//xH5W1Z5zaXTeq\n0yxHu4/3fbWe7dXy/4jancNSStk/J4OJf24OOeSQUB9wwAGhnjFjRqi9t4737rn44otD/YUvfCHU\np5xySqhPPvnkyvGW5gDnr6/d9ds6MWvWrFB7T76dd9451B/+8IdDff/994fa99Udd9wR6iuvvDLU\nzWtRSuvue5+/zjvvvFD7nON5s9mzZ4fa16ss/XszmK1cuXK9c1i7n663Z9VLJM1IKS2SdHSjBoC6\nYw4D0LaWO3jnnH8h6ReNn1+XNL2vBgUAvY05DMDGooM3AABAhX5dG66Uheh0LbhS5slvL23P+4j8\n8Ic/DPWBBx4Y6uZr7n693ftyeJ8Mz8R4xsm35xkmfy2+1plfz/drzn4Nu5QZ8vt735BSXynvK+IZ\no7feequyLq3z57W/t/583kPG3y9/vj/84Q+hbl7f7a677gq3+VpwpXX5Sv3BXKe/NMEvXQxNnpvs\n7u4O9bPPPhvq+fPnh/pzn/tcqL33z7//+7+H2jNOjz/+eKg9l1Nau7O0vqTXPmf5cdYJzyg678nn\n85vPx3fffXeoTz311FB7htVfq+/riy66KNSf/vSnQ/3Zz3421F/96ldD/dOf/jTUS5fG33fozX05\nWPHNEgAAQAVOlgAAACpwsgQAAFChXy9E9vXaPu2ud+a5Gr/G7zmcW265JdS77LJLqJuvE/tzeb+W\nnXbaKdQ+dn9uv2a+1VZbhdrHXsog+fP540sZJ187yK/J+zVur73vka+HVurB5RmmUh7NeR+Sd7/7\n3aH21+sZsptvvjnUM2fO7PnZsyClfdFKQzSgXX4MPP/886H23jr33XdfqP/t3/4t1Mcee2yoff1F\nz7n4MVVar9LnIO8153NWqRec8+PM+6E150h9W55pXL58eeVz+dqSPlbfFz6/7LHHHqHec889Q/2z\nn/0s1L7vvY+T58c88+Svh4zSuvhmCQAAoAInSwAAABU4WQIAAKgwqC5MttuHqbRWkPez8eu0nqNx\n//Iv/xLq0aNH9/z8qU99qnJbntl55ZVXQu1rq/n1fM8seSbKM0u+vXbzXF57psdfn68V59f8vc+R\nv7feh6qUwfK+JN63yu9fWnfKM1Fz5swJtb/3L7/8cs/PpT5KpXXt2v2cAxvjwQcfDPVxxx0X6nvv\nvTfUnlP81re+FerDDjss1H/1V38Vap/jli1bVjmeqVOnhtrnzB122CHUpTnK56Srr7461BMmTOj5\n+Ywzzgi3+THnPaueeOKJUN94442h3meffUL9/ve/P9Q+P3qmyN8bfy+8T9I3vvGNUHtfOM8o0Wut\njG+WAAAAKnCyBAAAUIGTJQAAgAq1Whuut7dfynb4NW3PMHkOyHMu3lvo61//es/P3jfj85//fKi9\nR8mKFStC7dfXfSx+jdszTL52nO8bX+dp7NixquJ9nnz8Pj7ft943yl9fKSPl76VnpPz+vn/8/t4H\nxfuUfPvb3w615xv89TVnpvw212nmqN3jqNM1FTE0+fvufeT22muvUPvacTNmzAj1SSedFOpLLrkk\n1Mcff3yoFy9eHGrvDdScA5Sk/fffP9QnnnhiqM8+++xQ/8M//EOofY7z8Tc/3vNVf/mXfxnqc889\nN9SbbbZZqD2z5PPxq6++Gupf/epXoS71vPIM0osvvlg5nr7ucbgp4JslAACACpwsAQAAVOBkCQAA\noMKA9lnqtJ9Mqa9Su9dpvf+NXxcuZT2ab//ud78bbvNr0hdddFGoPR8wfvz4UPu6Tp7BKa3DVMpj\njRs3LtR+jd37enhmyTNAnsHyPkieYfLx+DV3X8vN17XydaTcww8/HOq77ror1LNnzw71woULQ+37\n219PM9/XpXX3/P6lz2m7maPSmlrtHhcYmvxztGDBglD7Z96PWc8YeW+3yZMnh9o/d96ryNfe9F5D\nzesxStLTTz8d6ksvvTTUp556aqgnTpwY6uY+Uf7cd999d6h9XTnPTHrG84tf/GKofT7de++9Q/3Y\nY49Vbs/nd38vyCH2Pr5ZAgAAqMDJEgAAQAVOlgAAACqkvry2mVIKG/eshmc52s1etNs/ptOMk4/X\nrxM3Z09KORXv+TF9+vRQew8Q7xu02267hbp0zdxrX3vN12rza+LeF8RzNp5Z8j4l3pPK8wzel8rX\nZvPaexl5zxbvO+I9ZDwz4Ou5eUbK+0JVZaQ84+OfhXb1dl+mdrfn9y+tq5dzHhJNXFJK2Y9x/H+l\n7JtnnLxPkh/zvl7Zpz/96VD7HOpz4EMPPRRq7010zjnnhNqP8VmzZvX87HlVz4R6hrSUKcXgsXLl\nyvXOYXyzBAAAUIGTJQAAgAqcLAEAAFQY0MxS6Zp3p+vZ+P1L/WI6XS+nOZviuY4PfvCDofa1yHxd\nJO/zMWHChFAffvjhlbdvu+22lbVfU3/ppZcq7+99kebNmxfq3/zmN6H2a/ie8fG16byvk/cxKa3j\nV+J9oUo5nFLezT9LzZ+ddj9HflxUbXt92s0kdbpWHJkl9AX/HE6ZMqXydn9vHn300VD7epmjRo0K\ntecWFy1a1POzZyyx6SCzBAAAsBE4WQIAAKjAyRIAAECFAc0slfocDTbNr8f79ngfJO/j4Wuzeabo\nhRdeqHy88+vz3gfJ971nqLwHib8e76NUWput9N779j2f4Pf3DJPnfHx77WaQ2s0NVWn3sZ1m9drN\nJLWbWSrlxcgsoS90elwArdhQZqmlhXRTSosl/U7Sakmrcs4Hp5RGSfqhpJ0lLZZ0Ws75txvcCAAM\nAOYvAJ1q9TJclnRUznlKzvngxt9dIGlOznmipLmNGgDqhvkLQEfaySz511InS5rZ+HmmpFN6ZUQA\n0PuYvwBstJYySymlZyWt0Nqvsb+bc/5eSumNnPO7G7cnSa+/XTc9rq0+S0OZZ2A8w1S1zpy0bkbJ\n951ndEq9cDwT5dv32/298/GV1vnz8ZbG55mkTnsN1emz1ulabf2tlO9az/3rs7PV2fxFZgnYtHSU\nWZJ0WM751ZTStpLmpJQWNt+Yc85+YgQANcH8BaAjLV2Gyzm/2vjvMkk/kXSwpO6U0lhJSiltL2np\nhrfQs52NHykAbIRO5q9Vq1b1/Gm3azyA+lu9enU4zjekeLKUUto8pbRV4+ctJB0r6TFJN0o6s3G3\nMyXd0MK2Whg6APSOTuev4cOH9/zxS80ABr+urq5wnG9IMbOUUpqgtf83Jq29bPd/cs4XN3719keS\ndtIGfvWWzFLrWngfKm8v9THy2vsieTbDn++Pf/xjqEv/cPjr8f8rb3fdPgycwZxZ6nT+IrMEbFo2\nlFmq1UK6mzJOljhZqqvBfLLUCU6WgE0PC+kCAABsBE6WAAAAKrTaOgB9rNNLkn5ppHTZy+vSWnN+\n/3a37/jNSADAYME3SwAAABU4WQIAAKjAyRIAAECFfs0s1Xm9rqGm3cxQu+9Fu/cnowQAGKz4ZgkA\nAKACJ0sAAAAVOFkCAACoQJ8l9AkySgCAoYJvlgAAACpwsgQAAFCBkyUAAIAKnCwBAABU4GQJAACg\nAidLAAAAFThZAgAAqECfJQCV6JkFYFPHN0sAAAAVOFkCAACowMkSAABABU6WAAAAKnCyBAAAUIGT\nJQAAgAqcLAEAAFTgZAkAAKACJ0sAAAAVOFkCAACowMkSAABAhUG9NpyvWZVSGqCRAACAoaqlb5ZS\nSiNTStellJ5IKS1IKU1NKY1KKc1JKS1KKd2aUhrZ14MFgHYxfwHoVKuX4b4p6Wc550mS9pW0UNIF\nkubknCdKmtuoAaBumL8AdCT5pax17pDS1pLm5Zx3tb9fKOnInHN3SmmspDtzznvafbLVYdvDhnUW\nmeIyHND31qxZE+rSnJFzrs2B2On8NWLEiH4cLYCBtnLlyvXOYa1kliZIWpZSukbSfpJ+Lek8SWNy\nzt2N+3RLGtPpIDs9+WnhxK/tMfWW3n5tvf1a2t13dTtR7e/9U7U/+ntftPva6/be9bF+m78ADF2t\nfLXzDkkHSPrnnPMBkt6SfWWd186+1f/aAkD/Y/4C0LFWTpZekvRSzvnBRn2d1k4+SxpfXyultL2k\npaUNlb69AIBe1tH8tWrVqp4/q1ev7pcBA+g/q1evDsf5hhRPlnLOSyS9mFKa2Pir6ZLmS5ot6czG\n350p6YbStob41/0AaqbT+Wv48OE9f7q6uvp8vAD6V1dXVzjON6QY8JaklNJ+kq6W9E5Jz0g6S1KX\npB9J2knSYkmn5Zx/a49rK+C9KX/zVMoEYehoN1M00DzgXVKngLfU2fxFwBvYtGwo4N3SydLG4mSp\ndZwsbTo4WRocOFkCNj0bOlliuRMAAIAKnCwBAABU6Ne14fzyAr9dsmF9fVmOy34AALSGb5YAAAAq\ncLIEAABQgZMlAACACpws1VxfZYlyzuEPgPWre7ayzuOr89gkxteJOo9N6v3xcbIEABXa7TPV3+o8\nvjqPTWJ8najz2KTeHx8nSwAAABU4WQIAAKjQr8udABj6htJyJwM9BgD9r9/XhgMAABjsuAwHAABQ\ngZMlAACACn16spRSOj6ltDCl9FRK6Yt9+VwtjudfU0rdKaXHmv5uVEppTkppUUrp1pTSyAEa244p\npTtSSvNTSo+nlP5XzcY3IqX0QErpkZTSgpTSxXUaX2MsXSmleSml2TUc2+KU0qON8f2qhuMbmVK6\nLqX0ROP9nVqn8Q2UOs1hdZ6/GmOp7Rw2GOavxniYwzZubH0+f/XZyVJKqUvSVZKOl7SXpDNSSpP6\n6vladE1jPM0ukDQn5zxR0txGPRBWSTo/57y3pEMkndvYX7UYX855paRpOef9Je0raVpK6fC6jK/h\ns5IWSHo7iFensWVJR+Wcp+ScD278XZ3G901JP8s5T9La93dhzcbX72o4h9V5/pJqPIcNkvlLYg7b\nWH0/f3kn5976I+l9kn7eVF8g6YK+er42xrWLpMea6oWSxjR+Hitp4UCPsTGWGyRNr+P4JG0u6UFJ\ne9dlfJLGS7pN0jRJs+v23kp6TtJo+7tajE/S1pKeXc/f12J8A/ie1W4OGyzzV2M8tZzD6jh/NZ6f\nOWzjxtUv81dfXoYbJ+nFpvqlxt/VzZicc3fj525JYwZyMJKUUtpF0hRJD6hG40spDUspPdIYxx05\n5/k1Gt/lkv5OUnPb1rqMTVr7f2W3pZQeSimd3fi7uoxvgqRlKaVrUkoPp5S+l1LaokbjGyiDYQ6r\n5XtUxzms5vOXxBy2sfpl/urLk6VB15Mgrz0FHdBxp5S2lPRjSZ/NOb/ZfNtAjy/nvCav/Rp7vKT3\np5Sm2e0DMr6U0omSluac50lab4+fgd53kg7LOU+RdILWXp44ovnGAR7fOyQdIOmfc84HSHpL9pV1\nDfbfQBhUr7cu71Fd57C6zl8Sc1iH+mX+6suTpZcl7dhU76i1/2dWN90ppbGSlFLaXtLSgRpISmm4\n1k4y1+acb6jb+N6Wc14h6SZJB6oe4ztU0skppeckzZJ0dErp2pqMTZKUc3618d9lkn4i6eAaje8l\nSS/lnB9s1Ndp7eSzpCbjGyiDYQ6ry2dIjTHUfg6r4fwlMYd1ol/mr748WXpI0u4ppV1SSu+UdLqk\nG/vw+TbWjZLObPx8ptZeZ+93KaUk6fuSFuScr2i6qS7j2+bt3yZIKW0maYakeXUYX875yznnHXPO\nEyR9VNLtOeeP12FskpRS2jyltFXj5y0kHSvpsbqML+e8RNKLKaWJjb+aLmm+pNmqwfgG0GCYw2rx\nGZLqPYfVef6SmMM60W/zVx8Hr06Q9KSkpyV9qS+fq8XxzJL0iqQ/aW0W4SxJo7Q2VLdI0q2SRg7Q\n2A7X2mvVj2jtQTxPa3/zpS7jmyzp4cb4HpX0d42/r8X4msZ5pKQb6zQ2rb2m/kjjz+NvHwt1GV9j\nLPtpbej1N5Ku19rQZG3GN4D7pTZzWJ3nr8b4ajuHDZb5qzEm5rD2x9fn8xfLnQAAAFSggzcAAEAF\nTpYAAAAqcLIEAABQgZMlAACACpwsAQAAVOBkCQAAoAInSwAAABU4WQIAAKjw/wAoGkZpu2SaMwAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10328f190>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 10))\n",
    "axes[0].imshow(target.get_data()[:, :, 15])\n",
    "axes[1].imshow(source_before.get_data()[:, :, 15])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To show nipy in action, we first disturb the registration to give the algorithm something to do.\n",
    "\n",
    "This shift is just in the affine matrix, because the affine matrices record the current positions of the structural and functional image, relative to each other, and the world space of the scanner:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Move structural 10 mm to the left relative to the functional (mean)\n",
    "source.affine[0, -1] += 10\n",
    "source_before = resample(source, np.eye(4), reference=target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1128dcc10>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAEdCAYAAAAcgiN1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmQ3VWd9/HPoVsbJGxhSUJITECDsoWEIUBkyxgWZVi0\nGBBrJKNoOePAKJajeRxrSpnR8ZGawrFUfETJROdxkJJnIsgaExgWBSIkAQMhbJElpMMuYSAaOM8f\nucnc7yed3+mb28uvO+9XFUV/+/b93XO30yf39+nvSTlnAQAAoGfbDfYAAAAA6ozFEgAAQAUWSwAA\nABVYLAEAAFRgsQQAAFCBxRIAAECFthZLKaWTU0rLU0oPp5S+0FeDAoD+xvwFoLfS1vZZSil1SHpI\n0kxJT0taJOmcnPODfTc8AOh7zF8AWtHOJ0vTJD2Sc16Zc/6jpCsknd43wwKAfsX8BaDXOtu47lhJ\nTzbVT0k6ovkHUkq0Bwe2MTnnNNhj6AXmLwA96mkOa2ex1PJEklK8/e22ix9sDcbWK2+++eZm46gT\nxrf16jq2je+DjeOr65ZDG8f35ptvDvZQ+kOvHvQRI0Zo3bp16urq6u/xbLU6j6/OY5MYXzvqPDZp\n68e3du3aHr/fzmLpaUnjmupx2vCvsy3KOW+2YAKAQdCr+WvdunV64403tG7dOnV0dKizs50pE0Dd\nrF+/Xm+88Ubx59r5Z/dvJL0zpTQhpfRWSWdLurrqCiyUgKEnpRT+GyZ6NX91dXWpo6NDXV1dLJSA\nYaizs1NdXV2b/tviz23tDeSc16eUzpd0o6QOST8cin9JUvfJn/FtvTqPTWJ8g6mV+aujo2NAx9aq\nOo+vzmOTGF876jw2qe/Ht9WtA3p1cAtI1jGzBAw0fx/U/XXv4yuNd4gEvItSSnnEiBGDPQwAA2jt\n2rU9zmH1S78CAADUyDZ9En6o/Qsfw1PdXod1Gw8ADDY+WQIAAKjAYgkAAKACiyUAAIAKtcos9fWf\nKnvWonT8VrMaVcdr9bb72mDffkk7j21vjjeY12/3sS5dv6+fWzJJAFCNT5YAAAAqsFgCAACowGIJ\nAACgQq0yS32tv7Mj/Xnb7Rrs2y/p6/EN5nPdF9cfzNuq+2sFAAYbnywBAABUYLEEAABQgcUSAABA\nhWGdWQLqiL3XAGBo4ZMlAACACiyWAAAAKrBYAgAAqMBiCQAAoAKLJQAAgAoslgAAACqwWAIAAKhA\nnyUAlegLBWBbxydLAAAAFVgsAQAAVGCxBAAAUIHMEjDAyPwAwNDCJ0sAAAAVWCwBAABUYLEEAABQ\nYUAzS96vBQAAoO74ZAkAAKBCcbGUUro8pdSdUrq/6XsjU0rzU0orUko3pZR27d9hAkDrmL8A9IXe\nfLI0R9LJ9r3ZkubnnCdJWtCoAaBumL8AtK24WMo53ybpRfv2aZLmNr6eK+mMPh4XALSN+QtAX9ja\nzNKonHN34+tuSaP6aDwA0N+YvwC0pO2Ad97QjpiWxACGHOYvAL2xta0DulNKo3POq1NKYySt6c2V\ncs60DwAw2Ho9f61bt27T1x0dHersZIcoYDhZv3693njjjeLPbe0nS1dLmtX4epakeb25EgulvpNS\nCv8B6LVez19dXV2b/mOhBAw/nZ2d4X2+Jam0qWdK6T8kHSdpD204v/8Pkn4u6UpJ4yWtlHRWzvml\nHq4bDr7ddtv55b24K+iJP3Zszor+4q+tN998s/TztXljtzt/jRgxYuAGC2DQrV27tsc5rLhYageL\npf7DYgkDZSgvltrBYgnY9mxpsUQHbwAAgAqchEeP/FPAVj+54pMuAK3wOaP0ieb2228f6h122CHU\nr7zySku313x8n/8cZ0W2PXyyBAAAUIHFEgAAQAUWSwAAABXILNWUn0/3c+je88WbapXOqfv5/1Je\nwJUyTaW69JeRrY5nKP91YN3HTj4DfcFf1xMnTgz16NGjQ/3MM8+Eev369aH2jNKRRx4Z6nvvvTfU\n3d3dod59991D/frrr2/6er/99guX+V9FzpsXW3P5fXvLW95SeTmGHj5ZAgAAqMBiCQAAoAKLJQAA\ngApklgZIKbPj/Jy3n6/3440dOzbUq1atCrVnnP74xz+G2nMppYxUq8crnbP3DJUfzy9vtQ+UX04O\nB+hbPifstddeoT788MNDvWzZslCfdtppob7uuutCvWLFilB/8IMfDPX73//+UC9cuDDUp59+eqgP\nOeSQUK9cuXLT1x/+8IfDZeeee26oL7300lA/99xzoV6zJu7NvGDBglA/+OCDoS79PsDg4xkCAACo\nwGIJAACgAoslAACACoOaWap7f5lWVeVg/Hx+qc/QjjvuGGrP8Pz+978P9WuvvRbqt73tbZXH80yU\nH997kPzhD38IdXNPEmnzfZi89j5QO++8c6g7Ojoqa789P15pHykfvxvITFOrxx5u7xMMTf6eHTNm\nTKi7urpCfcopp4Ta+yr53m7Tp08P9RNPPBHqgw8+ONRnnHFGqL/+9a+H2nNHS5cuDfUnPvGJUP/5\nn//5pq/nzJkTLvM8lOer7rvvvlC//e1vD/Vjjz0W6vHjx4f6d7/7Xah9vnv66adD7fMf+h+fLAEA\nAFRgsQQAAFCBxRIAAECFAc0sDbfeNlX7AZX2NvNzzl6/9a1vDfWUKVNC/a53vSvUfo58jz32CLXv\no+QZKt83adSoUaF+9dVXQ+15g0cffTTU3nfEM1b3339/qD3j9Pzzz4faM1mlviSeWfL7W7rcnw/v\nczWQr+VW98VrVavXJy+xbdhpp51C7X2Tjj/++FD7nDFjxoxQX3/99aH296Bf33M7F154Yai/9KUv\nhfqEE04ItfcymjBhQqj33nvvUDfPAT7/ek8mz1Qedthhob7ssstC/cgjj4T6oosuCvWee+5Z+fO3\n3HJLqL0HlefFfL5cu3ZtqIfb7+KBwCdLAAAAFVgsAQAAVGCxBAAAUGFQM0t9fd50oPvP+Hni5vtT\nyiT5fffz7WeddVao991331B7nqC0V5rf/ssvvxzqcePGhdozOp4v8HP2fo7f+zb58aZNmxZqzyz5\n9V988cVQ+/h9L7y777471M8++6xa4Y9XKSPVzmtvoDNHpevT12nb5Lm9Aw88MNTe5+hjH/tYqL13\n27e+9a1Qey8hnxOWLFkS6vPPPz/UnsP03kN+ufdCOvTQQ0NdlVP8wAc+EC7zvnU+H/j84nmp4447\nrrL2TGdpvjvvvPNCffvtt4f6oIMOCrU/ti+88EKo/bkh07Q5PlkCAACowGIJAACgAoslAACACrXa\nG66UlWi330yr2Qv/ed9fzXM8zefgPXPjfZE+/vGPh9r3RfI+RuvWrQu1Z4S8L1Fp77fSXnL//d//\nXXn7pb3i/PY847TLLruEeuzYsaH2viDveMc7Km/Pnwvvi7J69epQ/+pXvwr1ypUrQ13KbHneofnn\nW83mtbovXat5glaPT0Zp2+Cv8cMPPzzU+++/f6h9jvHePyeffHKoL7jgglB/+ctfDvU555wTat8v\nzecEf11636SRI0eG2nvT+Zy5bNmyUDfvXTdixAhV8YyRZxyfeuqpUHu+y91www2h9h5Qnlnyffk8\nM+V5Wn8urr322lB7RsrvHxkmPlkCAACoxGIJAACgAoslAACACgOaWfJzzp4DKWk3a9Hq9T3349f3\nc/7NOZ0DDjggXObn6z0P4H2GXnrppVD7OXF/7HzvNs8IeQbJ71vpcs8s+c/7Y+PnvD0j5ZkAPyfv\nx/Nz+H5/PcM0evToUO+zzz6hnjp1aqi9T5NnmK655ppQe6aqOSPgz43X3qOl1b5GrWb5Stcnw7Rt\n8OfxiCOOCPXMmTND7Tm/hQsXhvqoo44K9e677x5qf897Jslzj95Lzvnr9Bvf+EaoPVPq7+FJkyZV\nHq95TvKxeUZx5513DvXPf/7zUPv86JlL7xH1zDPPhLo5PyVtPv/6Pnp+X/z3h+8dumbNmlAfeeSR\nofa9Ph977DFt64qfLKWUxqWUbk4pLUsp/Tal9LeN749MKc1PKa1IKd2UUtq1/4cLAL3H/AWgL/Tm\nNNwfJV2Ycz5Q0pGS/ial9G5JsyXNzzlPkrSgUQNAnTB/AWhbcbGUc16dc17S+HqtpAcljZV0mqS5\njR+bK+mM/hokAGwN5i8AfSG1kkdIKU2Q9F+SDpL0RM55t8b3k6QXNtZNPx8O7jmauvVu8N5G3pfD\nczF+Xru578dXvvKVcJnvq+QZJc/4lPYicp6h8QyQ52a8J4mfY/fj+X31Ph7+Ouru7g61P/eeAfDr\n+/j8572PiWe6/Hi77hrPsvjx/fH3+/vwww+H+oorrgh18/31vJc/d/5ctJtBGmie1XM553q9sRu2\nZv4q9dsZSny/MM/1+evWf37y5MmhXrBgQah9zvrSl74U6ttuuy3UPif6e7TVbJ33drvzzjtDfcwx\nx4Taews153z22muvcJnnoZzPPz7fuYceeijUN954Y6iPPfbYUH/nO98J9Wc/+9lQ33///aH2/Jnv\nJeoZJM9Med+mm266KdTep6nuc1Yr1q5d2+Mc1uu/hkspjZB0laRP55zDb/q84ZEaPo8WgGGF+QtA\nO3r113Appbdow0Tz45zzvMa3u1NKo3POq1NKYySt2fIRNsg51+7TJADDWzvzV/NfIXV0dGz2CSeA\noW39+vWbfTLYk978NVyS9ENJD+Scv9l00dWSZjW+niVpnl+3h2MVBwQAfaXd+aurq2vTfyyUgOGn\ns7MzvM+3pJhZSikdLelWSffpfz6q/l+S7pZ0paTxklZKOivn/JJdt18zS6Vz2n65Z0e89oyS97Yo\n5WAuuuiiTV8fffTR4TI/n+4rWe9p4pkZz0742Dxj5D1O/LHZc889Q+2ZJc88Oe8z4jmc5557LtR+\nf5y/Nvz4frn3KSlllrzPUimP5n2cfC+7O+64I9S33nrrpq+9R0lprzh/HfbmXzkDycfrr2VXp8xS\nu/PXcMos+V5qn//850N9/PHHh3qPPfYIdalP3pIlS0LtveR22GGH3g51q3jmysfruUTPPe622/9E\n1vy+9/c/9JcvXx5qv32fT32vUeeZoyuvvDLUnlHyPnbeF+ruu+8O9be+9a1QlzJaQ8mWMkvFfyrl\nnG/Xlj+BmrmF7wPAoGP+AtAX2O4EAACgAoslAACACi31WWr54DXLLHlA03tPeO+j0vVnz45Nf08/\n/fRNX5f6IvltPf/886H2DJL3MPFMkmdsPIPk4/GeUp4/8EyP99YpPZeeu/GeJp578dv3rIjf/gsv\nvFB5+z5+7ynj998fv1b3LWzuc3LvvfeGyzzD5I9FKThc6svU34ZyZqkdQz2z5HPI+eefH2rvneN7\nvXmvIe/dM23atFAP9QB8c27U97nzPNRQ48+N73XprwX/fXPdddeF2n+/eB5sKGu7zxIAAMC2iMUS\nAABABRZLAAAAFWp1krnVDFKrfZU8x+Ln9EvZkOZMkiSddNJJoW7O0Xiuw/v4lPYaKvUt8gyPH98v\n93Pufo7Zsxml/cz8sfTLS72CSvuheYbI+zSVnnt/PPy5Lo3fb6+031tzJsr3vPLHwntk+eWlPBfN\nXSGV+7695z3vCfWvfvWrUM+dOzfUnoP0Xmzep8n5nOPjq3vux3ONw8n3vve9UC9btizU/vvA95r7\nl3/5l1B73yV/LTXvszdc8MkSAABABRZLAAAAFVgsAQAAVBjUzFKrGaRWeS8gz6GU9hPzXJH3GfHc\nSnMWxc//e07F917zsZYyNaW9eDyz4302/HjeN8Ov73kDf+z8/pbuj/N9ozyz5OPx10opo1XKMJWO\n7+fg/blvvr7v6+RZEu/h4ntUOb+tge6zNNC3h97x9+Bhhx0W6kmTJoXa57PHH3881D6/rVq1qvL2\nli5dGupLL7001GPGjAl1896ZveFzQCnn6TyH6XOcZwObb8/7DLXbQ6r0u66/ffzjHw+1z5e+N93k\nyZND7X0Afa9Mf2y///3vh3o45Cz5ZAkAAKACiyUAAIAKLJYAAAAqDGpmqXQes90Mk/98q72Cjj/+\n+FDvs88+ofbzvrvtttsWx+Lnx/0cuF/uGSfPrXh+wO9LKbNTyvz47fnlnico9bgq9bzq6uoKtd8/\nr0uvDT++X+73x/njV+oZs/POO2/x2AcddFCoPdN08803h9rzUaX77sgYbRuaX3PS5n2C9t5771Cf\nffbZofb3iOdQ/D379a9/PdS33HJLqE844YRQ/+Vf/uXmg26yaNGiUI8aNSrU48ePr7y+z5E+vquv\nvjrUPgcecsghoT722GM3fe35KM8dev7V9xn1+cEfW89E7b///qEu7S3aLv/943OUzzme//r1r39d\nefzBzmj1Bz5ZAgAAqMBiCQAAoAKLJQAAgAoDmlnq6yxFKSfjl3vvH88oeSbpnHPOCbXvldTd3R3q\n5l5D3nfIb8vrUqbHz/mWzgF7BsiP/8orr4Taz1H7+FzpsXWlzJX3PXKlDFGr+wg6f7z9nL5nGMaO\nHbvF6zfvEShtnkfzx37ixImhvvPOOyvHVnpusG3wOcZfo1dddVWofb+vQw89NNQPPvhgqL13zsqV\nK0N95plnhtozUZ7zueKKK0J9xx13hNrfJ74fmd8/n+Nmz54d6h/96Eeh9izhbbfdFmrfO6+Zz2/+\n/j/qqKNC/eyzz4Z66tSpoT711FND3Wqmx+c336ttxYoVoZ41a1aoff71PJf3vbv99ttD/fOf/zzU\nRxxxRKhLffmGIj5ZAgAAqMBiCQAAoAKLJQAAgAoDmlka6F4Lpf3UvFfGMcccE+q3v/3tofb93Pw8\nb3OGwO+r9+Eo9QHyPkaeHyidA/YeJM4fi1LmyG+v1AfKMz9+/7z243m+wJ/L0j5+3sfElTJXXvv9\n8ddOc+09Sf7kT/4k1L///e9DvWzZslD79T3P1eo+VaWsYLt5LwwOf5145tJ7+3gfOH+PHHjggaH+\n4Q9/GOp//Md/DLW/Dr3v0vTp00PteyB6LsbH63Ok70/mc6Rn/zwj5X2f/D28/fbbb/Eyn5+eeOKJ\nytr35Zs7d26oS7lD/91Sypz6Y+OZotJz632WPvCBD4Taew767T/yyCOh9tcGmSUAAIBhjsUSAABA\nBRZLAAAAFQZ1b7h2lbIWpf3N/Oe9L4ifEy/t79Z8ztv7THimxjNIfmzPqZTyV359r/2ccWlvNL9v\nXvs5/FKuxR/L0nPnt+fPpec1fF8szzt47yN//P32S32n/PrNGQfPYvjryvvb3HfffaF+4IEHQl3K\nf7WLveSGJn8N+nvEX7NLliwJtb8ufc7yPQqvueaaUPvr1uc4fx379SdPnhzq97znPaG++OKLQ336\n6aeH+q677qq8/PDDDw/1X//1X4f6sssuC3Xze94fO+/p5I+97/f4yU9+MtSl56b5d4ckzZs3L9T+\nWH7wgx8M9bhx4yrH4/tP+u09/fTTofb5+ic/+Umon3nmmVA//vjjofbfP8MBnywBAABUYLEEAABQ\ngcUSAABAhVplllrt91LKyZQySyNHjgy198rwc/qeY3nxxRdD3XyOv5R/8tozSt5nY6+99gq135dS\nZskfS99XyjNQnvkp7VXXaibKe7z47ft4vU+V94zZeeedQ+15i1KfKs8F+fH9/ntmq6pHy6pVq0Lt\neYKDDz441L5Hlb+O/bb99lyr7ysyTEOD52h8DrnxxhtD/fDDD1ce7ytf+UqofT8xnzM8Y+Q5HO8n\nNnPmzFD7e3Tp0qWh9vfsN7/5zVD7e9Z7G733ve8N9ac//elQ+5zRvHee57s88+i/O3yvt2OPPTbU\nzz33XKh9n9HVq1eHeuHChaH+9a9/HWrfP9Lvuz/23gfJ51Pv9eZ5Lt8H9cknnwy1z0nDcf/Kyk+W\nUkrbp5TuSiktSSk9kFL658b3R6aU5qeUVqSUbkop7Vp1HAAYDMxhAPpC5WIp5/y6pBk550MlHSJp\nRkrpaEmzJc3POU+StKBRA0CtMIcB6AvFzFLOeeNnoW+V1CHpRUmnSdrYv32upDP6ZXQA0CbmMADt\nKmaWUkrbSbpX0n6SLs05L0spjco5bzyJ2S1pVH8MrtU9qbxPSKt9lW677bZQT506NdS+x5drPg/t\nGRzPyHjt99XPp3tfIc88lfZ28/P73gvIM1KeF3D+WJcySJ6r8RxOqc+SH8/zDmvWrAm1P/f+ePv4\n/Ry+7zPlj7/n25r3tfL76v1nPP/k/WA8j3DPPfeE2o/f6vuk9PPDbS+4wZzD+pO/5/w94DkVN2HC\nhFD/+7//e6gvueSSUHsvs9133z3U/p70Xjy/+MUvQu3vwbPOOivU/j7x7J/3efL+ZT6n3XDDDaG+\n9NJLQ92cQfX50fNgntG59957Q+3zl2ecPNPj8/1RRx0Vat+n1Pdma85bSdKRRx4Zav/d5n2Rzjgj\n/lth0aJFofa8ms9B/lgPR8XFUs75TUmHppR2kXRjSmmGXZ5TSiRCAdQScxiAdvX6r+Fyzi+nlK6V\ndJik7pTS6Jzz6pTSGElrClffeIxh969WAEPD1sxhzZ8gdHR0bPaJAYChbf369ZudeehJ6a/h9tj4\nVyIppR0knSBpsaSrJW38u9JZkub1fITNjtebHwOAPtHuHNbV1bXpPxZKwPDT2dkZ3udb/LnCccZI\nmts457+dpB/nnBeklBZLujKldJ6klZLOqjhGr7Xa36X0835e2HMy3vvCj+fnxPfff/9Qe5+l5myK\n9znyc9J+Pt75WF5++eXKy0t7p7lSn6ZSpspzOH556bH3n/f8gvO+IX6O3G/PH1/fC2nMmDGh9j5W\n/qbxvih+f5rzGz4276vklz/00EOh9nyW93jx11Jpn7t2DfF/5AzoHDaQ/D3uGSVf3H3ta18Ltb/O\nvve974X6u9/9bqhPOumkUHsfI3+de87Fe/Pst99+ofZeQaU5wzOk/p73fmWf//znQ+1zajMfu4/F\n5xPP+MyZMyfUn/3sZ0Pt98WPf84554Ta5z/viXX22WeH2p9bzyj5HOK/2zyT5a8tf+yvvfZaDXeV\ni6Wc8/2Spvbw/Rckzdz8GgBQH8xhAPoC250AAABUYLEEAABQodaJxdKeVaW931q9/sUXXxxq77Pk\n53lXrlwZ6ubeE6W910p9gHysfo68dLmfc/bx+O35Y+Hn/338nocoXb90jr6UT/Djee6ntDedZ6L8\n5z3/4T1sPBfk420ej9+29yTx52bx4sWh/s1vflN5W/5cunYzRuwNNzR4rx3vFfYP//APofbX+PXX\nXx/qRx99NNTee8f3K/Mcy+c+97lQe18m329x8uTJofbXufPeZ877k33xi18Mtc8hzb3RpDinlfaq\n9Pe0//wVV1wRas9n+WNb4r97/LH2+dDzaZ7PnTJlSqg//OEPh9rnGP/9Mn369FBfddVVPQ17WOGT\nJQAAgAoslgAAACrU6jRc6bSZK51W849G/c/L3/Wud4XaP+b1Px/3UzH+5+fNrQT8vviffpZaB/hY\n/bSQn+rxNgZ+uT82pdNofl+9VYCPz++vtyLw56K03YmP38fjWyv4c+WPV+m0nv+83z8fn1+/+WN5\nPx3iH5F3d3eH2rcyOPDAA0N9xx13hNqfu1a3P/H7MsRbA2yz/D3vtW/RcfTRR4f61FNPDbWf1vPX\nuL/u/D3q/D1/0UUXhbrd15231Jg9O+6FXIpl+GnC5vl8/Pjx4TJ/Tz/99NOh9vnd56cf/ehHoX74\n4YdD/Vd/9VeVt+f8d4+/py+88MLKy71VgM/Hn/jEJyqv76dkfb4cjvhkCQAAoAKLJQAAgAoslgAA\nACrUKrPU7p8st3oO3M+z+nli//Ny/3NK/3PO5uv7+Xo/lp/T9lyL5wP8tvy+lv6U3u+r1/6n9Tvs\nsEOo/bEptW3w+1/aTqWUNytlrkoZJ79/frk/t6V9wHx7lObjeVakeSuUnq47Y8aMUPtzec8994Ta\n3yel9w2ZpOHJ38OeWbrkkktC/Z3vfCfU/ufsp59+eqg/9rGPhdpbXpRyl66vX4cjR44MtbcO8MfH\nM05Vc7RnIEv5V58v/LHy2/b3+KuvvhrqUmbJ+WN7wAEHtHR9nx+9rcTll18eam8lsGZNj/tQDyt8\nsgQAAFCBxRIAAEAFFksAAAAVBjWz1Oo57FJWo9Snyc8ze7bkySefDLX32vCcjZ+nbj4P7eeAPXNU\n2u7EeR+PUl8Lv33PDHkGyPsoeWbJz1H7OXa/3Pk5+tJz532PfKsEH69njrxnjGcI/PInnngi1P74\neQ8uzzA8//zzm772fIL3VfJtJ/yxvuuuu0Ltr1P/+dI2EaX3GZmmocnfQ/6e+Pu///tQjx07NtTe\nK+iss84K9fz580PtuZUjjzwy1J/5zGcKI+5fpS2d/D3t/H1VpfSe88yT90UqHa/dvkWl+dXnuxtu\nuCHUEyZMCPVxxx0X6l/84heVtzcc8ckSAABABRZLAAAAFVgsAQAAVBjUzFKr/WFa3fOqdN72scce\nC/Wzzz4b6oMOOijUfo5/1apVoW7OxXhGqLRXW6mPUmmvtV133TXUvl+Y/7znrTxz5Ofc/Ry69yny\nTFWpL4lf7rfnPVz22WefUHuGyJ+Ll19+OdSeWfLXhueCPDPlmauVK1eGuvm14fettCeVP3aenXP+\nWihpdY9FDA3+ulq6dGmoR40aFepp06aF+txzzw219/P69re/HerDDjss1Keddlrl+HwOKu1H2S5/\nPPz2fA4q5UQHUrsZpcsuuyzUPh/6fHXeeeeF+re//W2ovc+Sz/f+89vCHMInSwAAABVYLAEAAFRg\nsQQAAFChVnvDuVYzTaWMkp8j9xyKZ0d8Dy/PonivoebzxH5bfn7cx1rqnVPV06mnsXnmxx8Lv9zP\nmXvPEr/cj+eZIM8glfpKee3nyD1T5Jklvz3fj82P549faXyeb/OeLc2ZMX/deHbELVy4MNSenWtX\nq/3I6Ls0NPlr9qmnngr1LrvsEuo5c+aE2l+nPkccfvjhofbsnfP3SCnH2N9a3cuuir9nfL702/L7\n7tdvN7918MEHh3rFihWh9r5J3jOrdP3HH3881AP93NUBnywBAABUYLEEAABQgcUSAABAhVrvDddq\nfxg/Z+99dfdcAAAUxUlEQVTnjf28sOdcbr311lCfdNJJlbfvOZjmXkW+l5n/7Isvvhhq36vMM0rO\ne+2UcijOHys/B13qWeUZK69Lj70r9UDx2/deRDvttFOovU+S93zx8flrwXnmwPMd++677xZvq3nf\nOEl68MEHQ33FFVeE+tFHHw11uxmiVveGI8M0NPnz5nsS+v5k1113XainT58eap+j/u3f/i3UU6dO\nDfXdd98d6qOOOirUpfeYvyf9feQG83Xq99X7Du23336h9p59nsF85zvfGeqHHnoo1FdffXWovefV\nMcccE2rvu3fxxReH2ntu3XHHHaF+6aWXQt1qb7fhiE+WAAAAKrBYAgAAqMBiCQAAoEKt+yy1mmny\nn/ccjJ939ctvv/32UN9///2hfve73x1qP+/cvL+bn58v7UPk98XP15fOGZceC79+aZ8mv36pT0ir\nfURK+zZ5jxY/nv+81/74le6f5yW8J40f3/fya74/fuxnnnkm1J4PWLBgQeXY/LHtb2SUhgfPLHkO\nZebMmaE+8cQTQ/2+972v8njee+yrX/1qqH3vuAsuuKByvKWMUkl/Zph8/vA815133hnq5557LtTe\no8rzX1/4whdC7Xtb+u+i++67L9Sec/zUpz4V6jPPPDPU//RP/xTqyZMnh9r7znkfvW0RnywBAABU\n6NViKaXUkVJanFK6plGPTCnNTymtSCndlFLatXQMABgszGEA2tHbT5Y+LekBSRs/55wtaX7OeZKk\nBY0aAOqKOQzAVitmllJK+0h6v6SvSvps49unSTqu8fVcSbeoF5NNq3tU9TCW0k0EpZyQH8/3hlu8\neHGojzjiiFCPGDEi1M25Fz/n6xkcz+w092iSNj9/32qmyW+/lNkpZYJ8vJ6j8drH48f34/l4/HJ/\nLv0cuo93xx13DHUpA+W3731K/PH0+9ecS/Lb8p5QP/vZz0LtPbY871bqgdXXhlufpb6cw4YSfx59\nPvPXuL+mzz333FBPmTIl1D/4wQ9CvWbNmlB7xslfxz4n+Xj9feTvwVZ/f7TDb9vvq2eMvM+Sz9cr\nV66sPL7/bvE+dosWLQq17902bty4UPt+k97XyXsMem849O6TpUsk/Z2k5lf6qJzzxndCt6TqnUIB\nYPAwhwFoS+ViKaX0Z5LW5JwXS+rxn5d5w3K+f/+pCwBbgTkMQF8onYabLum0lNL7JW0vaeeU0o8l\ndaeURuecV6eUxkhaU3mUhpzzkP9IH8CQ0tYc1ny6t6Ojo7gNEYChZf369cXIjlRYLOWcvyjpi5KU\nUjpO0udyzh9JKX1D0ixJ/7vx/3m9GVS7C6XSOenSOXGvSzmbn/zkJ6H2/X3e+973hrp5Ii3tTebn\nqL2vj/+856leffXVyrG7Ut+j0l5zXnuuxjNEr732WuXtlfoueQ+rUl8lv39+/T333DPUnhHz58P3\n8hs5cmSovWdN8/1dvXp1uMxfR8uXLw91KU/V6h6Jrb7PhltGqVm7c1hpP7M68+fR55Rf/vKXofac\nimeafD+yp59+OtTee+xDH/pQ5Xj8Pew5oL333jvUnuPx93jpdetzlN//5j5648ePD5d5HzbPKHlm\nyPd/9N8VkyZNCrU/NxMnTgy1P5bunnvuCfXXvva1UPscs2LFilD7fL0t6ezsDL9/tvS7tNU+Sxtn\n1a9LOiGltELSnzZqAKg75jAALev1Z8o55/+S9F+Nr1+QNLP6GgBQH8xhALYWHbwBAAAqDGhasXRO\nud294FrNdpSO98ILL4T6pz/9aaj9HH7zeU/PwPheQnvttVeoPRPjewv58TzD5PfFz7GXMkqlnial\nvdtef/31yss9k7THHnuE2vMInskqZbRKmSZ/bv32PDPmz5ffnp/jb97fzXuW+F5wpX35ehM2bNZu\nxmg4ZZTQe/468/3G/HXhly9ZsiTUs2bNCvWqVatC7b1+vNeQ53w893nccceFesaMGaG+5ZZbQn38\n8ceH2nNBvmfjtddeu+nrv/iLvwiXTZ8+vfJYnp9auHBhqE899dRQn3LKKaH2OcGv7/fN94Lz+cv3\n7fPfF7znW8cnSwAAABVYLAEAAFRgsQQAAFBhQDNL/d3PpdX9zjxX432XPIdz4403hnrChAmh/tSn\nPrXF2/K+Pt7Hw8fut+19fXbaaadQ+9hb3VfJr1/KOPk5es9I+Tl4r73vkfezKe0b5RmmVnsTeU+Y\n3XbbLdR+fz1Ddv3114d67ty5m772vEDpsfDXCjAYPIezbNmyUHuG84ADDgi1zwm+B+LnPve5UPuc\ndPnll4f65JNPDrXPeZ7jueCCC0J9/vnnh/rMM88M9YknnhjqefP+p9WWZ3y8B9RHPvKRUI8ePTrU\n3ufo5ptv3uJtSdK0adNC7Y/1vffeG2rPW7lSnzy0jk+WAAAAKrBYAgAAqMBiCQAAoMKQ2hWy1T5M\nnntx3mfEsySlfaG+//3vh3r33Xff9PUnP/nJymN5Zsd7kvjeZJ4p8vP3nonyzJIfr9U8l9ee6Snt\nFed9QLzPkT+33oeqlMHynizet8p/3mvnmaj58+eH2p/75n2ySn2USvvatfo6B/qDzxne28fnrEWL\nFoXa95bzPnXO96bbb7/9Qu05T9+b7qMf/Wio58yZE2qfo7wP0z777LPpa89r+V5qzRlFafP5xjNH\n3/jGN0Ltc4D3sNphhx1CPZz3bxwq+GQJAACgAoslAACACiyWAAAAKqT+zD+klMLBPXcz2OddSxkn\nzwGVci7Nl5933nnhMu8xMmbMmFD/7ne/C7XvTeb8/LtnmLxPiJ9Df+WVV0I9atSoytvzPk877rhj\nqP2x8kyT943yviSeSfLj+b5SL7/8cuXPN+cPpM0fL89TeB8Vzxj84Ac/CLXfv+bMlF/W6h6Gfa3d\nPRVd6X2Qcx4WgYqUUvZsHXrPX2elXm/77rtvqPfff//K43uu8cknnwz14sWLQ33ggQeG2ueM7u7u\nTV/fc8894TL/3eV95Up7UWLoWLt2bY9zGJ8sAQAAVGCxBAAAUIHFEgAAQIVBzSyV+iCVxlbKHLWa\n1Sj9fOn6zZf7+fkpU6aE+stf/nKofZ8lz/B4psn3MvPMkWd8PKPkfZDGjRsXas88+V5E3lfEM0Ce\nKfLHzs/x+/V97zbPTHkfEt/Xynsded+UW2+9NdTXXHNNqJcvXx5qzzz5a7mZ57NKr3vv91XKO7Sb\nOWo1T+HPFZklDATP/nnmyHsfOX+dl7KEze9LMkfbLjJLAAAAW4HFEgAAQAUWSwAAABUGNLPkWQ3P\ncpTG0kpmqKfL2804+Xi9T0dz9qSUU/E+QzNnzgz1CSecEGrfW+0d73hHqD3j5Bkfr71HiWczPJfi\nmSXP2XjmyPeN8r4kvu/Uc889F2rfm81rzx+sXLky1N5z5cYbbwz12rVrQ+0ZJ89MecbLM1LNPONT\nlW/qjXbfo+1mnPznySwBGK7ILAEAAGwFFksAAAAVWCwBAABUGNTMUjt9jHp5+6H2LEnp51vVnE3x\nXMcpp5wSat+L7M477wz1hAkTQj1x4sRQH3300ZWX77nnnpW19w166qmnKn/eewf5vktLly4NtfeF\n8oyP933yvk6rV68OtWeUPDNV4nvZlXI4pbxb1V5Qrb6O/H3R6j5TrWaS2t0rjswSgOGKzBIAAMBW\nYLEEAABQgcUSAABAhUHNLJX6HA01zffH+/Z4H6TXXnst1GPHjg21Z4qeeOKJyuu70t5t/th7hsr7\nDPn98T5KVX2Hero9r/34Vfs2SZtnmDzn48drNYPUam6oSqvXbTer12omqdXMUikvRmYJwFC1pcxS\nZ08/7FJKKyX9XtIbkv6Yc56WUhop6aeS3i5ppaSzcs4vbfEgADAImL8AtKu3p+GypONzzlNyztMa\n35staX7OeZKkBY0aAOqG+QtAW1rJLPnHUqdJmtv4eq6kM/pkRADQ95i/AGy1XmWWUkqPSXpZGz7G\n/j8558tSSi/mnHdrXJ4kvbCxbrpeS32WhjPPwHiGqWqfOWnzjJI/dp7RKfXC8UyUH98v9+fOx1fa\n58/HWxqfZ5La7TVUp9dau3u1DbRSvquHn6/Pg6325i8yS8C2pa3MkqT35JyfSSntKWl+Sml584U5\n5+wLIwCoCeYvAG3p1WIp5/xM4//PppT+U9I0Sd0ppdE559UppTGS1lQeZMP1a/UvfADDXzvz17p1\n6zZ93dHRsdknuACGtvXr1/dqR4hiZiml9LaU0k6Nr3eUdKKk+yVdLWlW48dmSZrXi2MVBwQAfaXd\n+aurq2vTfyyUgOGns7MzvM+3pJhZSilNlPSfG48r6f/mnP+58ae3V0oary386S2Zpd7rxfNQeXmp\nj5HX3hfJM1N+e83/wpY2zyg5vz++cm913z4MnqGcWWp3/iKzBGxbtpRZqtVGutsyFksslupqKC+W\n2sFiCdj2sJEuAADAVmCxBAAAUIHEYk20e0rST42UTnt5Xdprzn++1eO7uvcWAgBgIz5ZAgAAqMBi\nCQAAoAKLJQAAgAoDmlmq835dw02rmaFWn4tWf56MEgBgqOKTJQAAgAoslgAAACqwWAIAAKhAnyX0\nCzJKAIDhgk+WAAAAKrBYAgAAqMBiCQAAoAKLJQAAgAoslgAAACqwWAIAAKjAYgkAAKACfZYAVKJn\nFoBtHZ8sAQAAVGCxBAAAUIHFEgAAQAUWSwAAABVYLAEAAFRgsQQAAFCBxRIAAEAFFksAAAAVWCwB\nAABUYLEEAABQgcUSAABAhSG9N5zvWZVSGqSRAACA4apXnyyllHZNKf0spfRgSumBlNIRKaWRKaX5\nKaUVKaWbUkq79vdgAaBVzF8A2tXb03D/Kum6nPO7JR0iabmk2ZLm55wnSVrQqAGgbpi/ALQl+ams\nzX4gpV0kLc4572vfXy7puJxzd0pptKRbcs7vsp/JVodjb7dde5EpTsMB/e/NN98MdWnOyDnX5o3Y\n7vw1YsSIARwtgMG2du3aHuew3mSWJkp6NqU0R9JkSfdI+oykUTnn7sbPdEsa1e4g21389GLh1/KY\n+kpf37e+vi+tPnZ1W6gO9ONT9XgM9GPR6n2v23PXzwZs/gIwfPXmo51OSVMlfTfnPFXSq7KPrPOG\n2bf6ty0ADDzmLwBt681i6SlJT+WcFzXqn2nD5LO68fG1UkpjJK0pHaj06QUA9LG25q9169Zt+m/9\n+vUDMmAAA2f9+vXhfb4lxcVSznm1pCdTSpMa35opaZmkayTNanxvlqR5pWMN84/7AdRMu/NXV1fX\npv86O4d0pxUAPejs7Azv8y0pBrwlKaU0WdIPJL1V0qOSPiqpQ9KVksZLWinprJzzS3a9lgLe2/In\nT6VMEIaPVjNFg80D3iV1CnhL7c1fBLyBbcuWAt69WixtLRZLvcdiadvBYmloYLEEbHu2tFhiuxMA\nAIAKLJYAAAAqDGhi0U8vvPHGGwN580NKf5+W47QfAAC9wydLAAAAFVgsAQAAVGCxBAAAUIHFUs31\nV5Yo5xz+A9CzunfurvP46jw2ifG1o85jk/p+fCyWAKBC3f8Qpc7jq/PYJMbXjjqPTer78bFYAgAA\nqMBiCQAAoMKAbncCYPgbTtudDPYYAAy8Ad8bDgAAYKjjNBwAAEAFFksAAAAV+nWxlFI6OaW0PKX0\ncErpC/15W70cz+Uppe6U0v1N3xuZUpqfUlqRUroppbTrII1tXErp5pTSspTSb1NKf1uz8W2fUror\npbQkpfRASumf6zS+xlg6UkqLU0rX1HBsK1NK9zXGd3cNx7drSulnKaUHG8/vEXUa32Cp0xxW5/mr\nMZbazmFDYf5qjIc5bOvG1u/zV78tllJKHZK+LelkSQdIOiel9O7+ur1emtMYT7PZkubnnCdJWtCo\nB8MfJV2Ycz5Q0pGS/qbxeNVifDnn1yXNyDkfKukQSTNSSkfXZXwNn5b0gKSNQbw6jS1LOj7nPCXn\nPK3xvTqN718lXZdzfrc2PL/Laza+AVfDOazO85dU4zlsiMxfEnPY1ur/+cs7OffVf5KOknRDUz1b\n0uz+ur0WxjVB0v1N9XJJoxpfj5a0fLDH2BjLPEkz6zg+SW+TtEjSgXUZn6R9JP1S0gxJ19TtuZX0\nuKTd7Xu1GJ+kXSQ91sP3azG+QXzOajeHDZX5qzGeWs5hdZy/GrfPHLZ14xqQ+as/T8ONlfRkU/1U\n43t1Myrn3N34ulvSqMEcjCSllCZImiLpLtVofCml7VJKSxrjuDnnvKxG47tE0t9JerPpe3UZm7Th\nX2W/TCn9JqX0icb36jK+iZKeTSnNSSndm1K6LKW0Y43GN1iGwhxWy+eojnNYzecviTlsaw3I/NWf\ni6Uh15Mgb1iCDuq4U0ojJF0l6dM551eaLxvs8eWc38wbPsbeR9KxKaUZdvmgjC+l9GeS1uScF0vq\nscfPYD92kt6Tc54i6X3acHrimOYLB3l8nZKmSvpuznmqpFdlH1nX4PEbDEPq/tblOarrHFbX+Uti\nDmvTgMxf/blYelrSuKZ6nDb8y6xuulNKoyUppTRG0prBGkhK6S3aMMn8OOc8r27j2yjn/LKkayUd\npnqMb7qk01JKj0v6D0l/mlL6cU3GJknKOT/T+P+zkv5T0rQaje8pSU/lnBc16p9pw+SzuibjGyxD\nYQ6ry2tIjTHUfg6r4fwlMYe1Y0Dmr/5cLP1G0jtTShNSSm+VdLakq/vx9rbW1ZJmNb6epQ3n2Qdc\nSilJ+qGkB3LO32y6qC7j22PjXxOklHaQdIKkxXUYX875iznncTnniZI+JGlhzvkjdRibJKWU3pZS\n2qnx9Y6STpR0f13Gl3NeLenJlNKkxrdmSlom6RrVYHyDaCjMYbV4DUn1nsPqPH9JzGHtGLD5q5+D\nV++T9JCkRyT9r/68rV6O5z8krZL0B23IInxU0khtCNWtkHSTpF0HaWxHa8O56iXa8CZerA1/+VKX\n8R0s6d7G+O6T9HeN79difE3jPE7S1XUamzacU1/S+O+3G98LdRlfYyyTtSH0ulTS/9OG0GRtxjeI\nj0tt5rA6z1+N8dV2Dhsq81djTMxhrY+v3+cvtjsBAACoQAdvAACACiyWAAAAKrBYAgAAqMBiCQAA\noAKLJQAAgAoslgAAACqwWAIAAKjAYgkAAKDC/wclF0/VMNhCIgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1082a2950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 10))\n",
    "axes[0].imshow(target.get_data()[:, :, 15])\n",
    "axes[1].imshow(source_before.get_data()[:, :, 15])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set up the registration in nipy, choosing our image similarity metric (in this case Normalize Mutual Information), the interpolation to use when generating a new image at each step of the registration (Random perturbation of the voxel coordinates), and using the optimization method Powell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial guess...\n",
      "translation : [ 0.  0.  0.]\n",
      "rotation    : [ 0.  0.  0.]\n",
      "\n",
      "Optimizing using fmin_powell\n",
      "translation : [-13.0562   1.       6.4179]\n",
      "rotation    : [-0.0026 -0.0001 -0.0062]\n",
      "\n",
      "nmi = 0.267171289077\n",
      "\n",
      "translation : [-12.6742   2.       6.4392]\n",
      "rotation    : [-0.0021  0.0099  0.0038]\n",
      "\n",
      "nmi = 0.2681862398\n",
      "\n",
      "Optimization terminated successfully.\n",
      "         Current function value: -0.268246\n",
      "         Iterations: 2\n",
      "         Function evaluations: 83\n"
     ]
    }
   ],
   "source": [
    "registration = HistogramRegistration(source, target, similarity='nmi', interp='rand')\n",
    "transform = registration.optimize('rigid', optimizer='powell')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This results in an estimated affine that goes from world (scanner, millimeter) coordinates in the source images, to world coordinates in the target image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0.9999,  -0.0038,   0.0099, -12.6742],\n",
       "       [  0.0038,   1.    ,   0.0021,   2.    ],\n",
       "       [ -0.0099,  -0.002 ,   0.9999,   6.4392],\n",
       "       [  0.    ,   0.    ,   0.    ,   1.    ]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transform.as_affine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As ever you can inspect the registration object by using `?` to get help and doing tab completion in IPython:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# HistogramRegistration?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We make a new resampled copy of the structural image, at the same voxel grid as the structural.  It's nearer to its original position now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "source_after = resample(source, transform.inv(), reference=target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x11b2d6810>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAEdCAYAAAAcgiN1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUXFWZ9/HnSYsSAgIhkIRcIAgxIRBuA9GABiThIorI\nHXQmw6CLGXgBx1EHXOqKawGO4ojOKAo4AsOIGkGQcA+IIAICJoGQECKXQAKkwyVEZLy06f3+kUpP\nPb90n13V1ZdT3d/PWln001V1zq5z6uzenP3Usz2lZAAAAOjckP5uAAAAQJkxWAIAACjAYAkAAKAA\ngyUAAIACDJYAAAAKMFgCAAAo0NBgyd2PcPdl7v47d//XnmoUAPQ2+i8AtfLu1lly9xYze8rMZprZ\ni2b2iJmdklJ6sueaBwA9j/4LQD0aubN0gJk9nVJakVJqM7Mfm9lHeqZZANCr6L8A1OxtDbx2jJmt\nrIpXmdm06ie4O+XBgUEmpeT93YYa0H8B6FRnfVgjg6W6OxL3uP8hQ+KNrf5YeqW9vX2TdpQJ7eu+\nsrZt43WwsX1lXXJoY/va29v7uym9oaaDvtlmm9n69eutpaWlt9vTbWVuX5nbZkb7GlHmtpl1v31t\nbW2d/r6RwdKLZjauKh5nG/7vrEsppU0GTADQD2rqv9avX98xWHT3Ug6+AXRfe3t7Tf/D2siV/6iZ\n7ebuO7v7283sJDO7qegFDJSA5uPu4d8AUVP/1dLSYkOGDOn4L4CBZeP1vfFfV7p9Zyml9Fd3/39m\ndoeZtZjZfzXjN0nK3vnTvu4rc9vMaF9/qqf/KvtxKHP7ytw2M9rXiDK3zazn29ft0gE1bVwSJMuY\nswT0Nb0Oyv651/bl2tskCd5Z7p4222yz/m4GgD7U1tbWaR/GfWUAAIACjSR4N71m+z98DExl+xyW\nrT0A0N+4swQAAFCAwRIAAEABBksAAAAFSpWz1NNf9dNci9z2683VKNpevfvuaf29/5xGjm0t2+vP\n1zd6rHOv7+lzS04SABTjzhIAAEABBksAAAAFGCwBAAAUKFXOUk/r7dyR3tx3o/p7/zk9Xoq+H891\nT7y+P/dV9s8KAPQ37iwBAAAUYLAEAABQgMESAABAgQGdswSUEWuvAUBz4c4SAABAAQZLAAAABRgs\nAQAAFGCwBAAAUIDBEgAAQAEGSwAAAAUYLAEAABSgzhKAQtSFAjDYcWcJAACgAIMlAACAAgyWAAAA\nCpCzBPQxcn4AoLlwZwkAAKAAgyUAAIACDJYAAAAK9GnOktZrAQAAKDvuLAEAABTIDpbc/Qfu3uru\ni6t+N9zd57v7cne/09236d1mAkD96L8A9IRa7ixdaWZHyO/OM7P5KaWJZnZ3JQaAsqH/AtCw7GAp\npfQrM1srvz7azK6u/Hy1mR3Tw+0CgIbRfwHoCd3NWRqZUmqt/NxqZiN7qD0A0NvovwDUpeEE77Sh\nHDEliQE0HfovALXo7mCp1d1HmZm5+2gzW1PLi1jmAUAJ1Nx/rV+/vuNfe3t7nzUQQN9ob28P13lX\nujtYusnMZld+nm1mN9byIuos9Rx3D/8A1Kzm/qulpaXj35AhVFoBBpohQ4aE67wrnrvb4+4/MrMZ\nZjbCNszvf8nMfm5mc81svJmtMLMTU0pvdPLasHHtbPgj33167Lhrh96in63cHZaUUmku7Eb7r802\n26zvGgug37W1tXXah2UHS41gsNR7GCyhrzTzYKkRDJaAwaerwRL3lQEAAAr06dpwaB56F7DeO1fc\n6QIADBTcWQIAACjAYAkAAKAAgyUAAIAC5CyVlOb8aA7R294WT50W08p901C/0aT7y+Uc5XKacnHu\nm5H1tqeZvx1Y9rbzrVUAgx13lgAAAAowWAIAACjAYAkAAKAAOUt9JJezo7Ry8F//+tfC7Y0ZMybE\nL730Uog1x6mtrS3EmpeSy5Gqd3u5PBzNodLt6eP11oHSx8nDAQDUijtLAAAABRgsAQAAFGCwBAAA\nUKBfc5bKXl+mXkV5MJrjk6szNGzYsBBrDs/vf//7EP/xj38M8RZbbFG4Pc2J0u1vt912If7LX/4S\n4j/96U8hfvPNNwtjrQP1zne+M8QtLS2Fse5Pt6efHc1x0varvsxpqnfbA+06ATqTuwZzeYvbbrtt\niLVPGzp0aIhXr14d4uo+WnNEtT/C4MOdJQAAgAIMlgAAAAowWAIAACjQpzlLA622jc6xV+cB5dY2\n05wbjd/+9reHeJ999gnxpEmTQrzTTjuFeMSIESHW+XrNoWptbQ3xyJEjQ/zWW2+FePPNNw/xM888\nE+JXX301xJpjtXjx4hBrjtNrr70WYs3JytWp0vwGfb+5x/V8aA5DX36W610Xr171vl6PDdCZXA6S\n5gFpnqT2WVtvvXWItY/TvMSDDjooxP/7v/8b4jvuuKPL9owbNy48pv3Zc889F2Ltr9RA+9s3GHFn\nCQAAoACDJQAAgAIMlgAAAAp4b9Zscfewcc0L6en6MX1df0bziqrfT705LzNnzgzxiSeeGOJddtkl\nxFtttVWIc2ulaXvWrVsXYq3LpO3VHB/NN9CcJ61xotv7wx/+EGLNWdLXr127NsTafl0L7+GHHw7x\nK6+8YkW07pQer1yOVD2fPT33jV4HPZ2zlGtPDTWrBkSChrsn/Vyga3rNTJ06NcSa5zhr1qwQL126\nNMSnnnpqiG+77bYQa5/4vve9L8SaR/T1r389xPr36Oijj+74+f3vf394bPny5SHW/uQHP/hBiLW/\nW7NmTYg1B0rbgv7T1tbWaR/GnSUAAIACDJYAAAAKMFgCAAAoUKq14XK5Eo3Wm6k3F0Sfr+uraR5P\n9Ty15txoXaRPfOITIZ4+fXqIdX7/z3/+c4g1X0rn53Nrv+XWktOaJLr/3Fpxuj/NcdKaKWPGjAmx\n5jTtuuuuhfvTc6H5EroO1AMPPBDiFStWhDiXs6U5BtXPz32uVb3r0tWbo1Tv9ll7DrV417veFWKt\nezRjxozC+IUXXgix1llauXJliLUP0j5r1apVIZ4zZ06Ijz322BDrNf7GG290/PyOd7wjPKY5pbfc\nckth284444wQa3/185//PMTa/2j/q9ek9t+5v53UeWocd5YAAAAKMFgCAAAowGAJAACgQJ/mLOk8\nquaB5DSaa1Hv6zXvR1+vc97VeTq77757eEznz9/97neHWOsMVc+fm21aw0SPna7dpjlCOgeu7y33\nuM7J6/P12OhacDrHvuWWW4ZY6ybp9nbccccQ6/vVnIBRo0aFeOzYsSHed999Q6x1mjSHYN68eSHW\nnKrqHAc9NxprzaZ66yzVm8uXez05TOiMnnfNuzzppJNCPH78+BDfc889IZ4wYUKI9Tp49tlnQ6zr\nU+o1euSRR4b4vPPOK3z+xIkTQ6zvb/78+R0/a/938cUXh1jzrQ499NAQf/SjHw3xZZddFmLNgdJ8\nrSuvvDLEJ5xwQog1Z0r7V+3/n3/++RBr/4y87J0ldx/n7ve4+xJ3f8Ldz6n8fri7z3f35e5+p7tv\n0/vNBYDa0X8B6Am1TMO1mdk/p5SmmNl7zOwsd59sZueZ2fyU0kQzu7sSA0CZ0H8BaFh2sJRSWp1S\nWlT5+Q9m9qSZjTGzo83s6srTrjazY3qrkQDQHfRfAHpCXTlL7r6zme1jZr8xs5EppY0LgrWa2cgu\nXtZjejt3YujQoSHWWkaaF6N5RNVrFZ1//vnhMZ0v1zljjbWGiLZF86V0/l9znjRvRrev+9ecHK2r\npHVI9NxonaXc83X7w4cPD7HmKGldKT0Xuv1ttomzLFoTRutaaY6T5k/8+Mc/DnH12nia76CxnotG\nc5AaNVhykvq7/yo7/RxMmTIlxKecckqItT/UPM13vvOdIdb1LPfee+/C1+t1M23atBD/8Ic/DLHm\nPH3sYx8LsdZqu//++0M8bty4jp8152fPPfcM8TnnnBNirbum237yySdD/Hd/93eFbT/++ONDrDWt\nNMdJt/ezn/0sxJpjtddee4X4iSeeCLH2p6jj23DuvqWZXW9m56aUQjZy2nCVDY4eF0DTof8C0Iia\nBkvuvplt6GiuSSndWPl1q7uPqjw+2szWdPX6jQbL/8ECKI9G+q/169d3/Kv327sAyq+9vT1c512p\n5dtwbmb/ZWZLU0rfrHroJjObXfl5tpndqK/tZFu5pwBAj2m0/2ppaen4p1PdAJrfkCFDwnXelVpy\nlg40s4+b2ePuvrDyu/PN7N/MbK67n25mK8zsxMaaXL9cvRh9XDu7XL0bzdvR7W277bYh/pd/+ZeO\nn3VtMs3hyeUgvfrqqyHWukRa90jbmls7SNujseYjqFwejuYoaU6Sxvoh1cd1xJ97XM+VHm99vdLj\noXP8erzvu+++jp+feeaZ8JgeGz0X+jksW75Ak/9PTmn7r7Kpzrk0MzvssMNCfPDBB4f42muvDfHl\nl18e4rPOOivEeg1q/6l9ntZJytVd0vZqHuKiRYtCrHlE1dv/xS9+ER7Tukkvv/xyiEePHh1izUHS\n967H8umnnw7xPvvsE+KddtopxNU5kmab5sRqvpjma+nanJoT9dRTT4W4yfuAHpEdLKWU7reu70DN\n7OL3ANDv6L8A9ATuKwMAABRgsAQAAFCgT9eG62+aG6J1P3R9NqV5OmeffXaIp0+fXvNrczlBmlOj\nz9ecJI117TZ972vWxC//5HKONKdK26Nz2ppjtXbt2hBrTlAuh0n3n8vJ0hyB3DeZNE9IX6/707Xl\nqs+v5mJoDpMeCz03+t5zdZmAWujnRvNcTj755BDPmjUrxLq224UXXhhizWHSvB5d31FrpSnNA9L2\nDxs2LMRaS616rTezTWslaW2h6tpu//3f/x0e22GHHUJ86623hviiiy4K8ezZs0OsdZO0/zv88MND\nrO9dn685RppD+YlPfKKwfaeddlqIH3/88RBr3b3vfe97Idac2cGAO0sAAAAFGCwBAAAUYLAEAABQ\noFQ5S/XWTaq3rpLmoeg8by435CMf+UiIdZ759ddf7/g5V7dI59uVtkVrkORyiPRxzZnSOWnNMdJj\np+3J5QTlagXl1kPTteA0pyl37vV46LnOtV/3p8dPj8+oUaM6ftbcCD0WOt+vj+tnJ1enCaiFXiMz\nZswojHWtNu0zHn300RAvWbIkxLr2m6q3Irp+7nX9Rr2uNGdqwYIFIdbrdO7cuR0/a36V1h3SukjV\nddbMNq3LpnWNNC9R861yOaSf/OQnQ7zzzjsXPn/mzFglQ9f9Gz9+fIiXLVsW4qOOOirEuvbcYMCd\nJQAAgAIMlgAAAAowWAIAACjgvVmzxd3DxosWqas8P8S5tuWer3km+nhuPTHNK7rkkktCrOu/VdcN\n0bo8Op+ueSl6bHJ1mXJ5MPp83V4uB0hzcvT1uu6Svl99P9o+3f/QoUNDrDlLuXX69Hhqe3X7unZS\nrm6Uvj89PtX5HNo2ze3QfIdVq1ZZkXqvi96m+XAqpTQgkqrcPennqJlNmzYtxO973/tCrNeArkd2\nxBFHhLi6LlFPqPdzrp9DrfukeTg/+clPQqzX8Ic+9KGOnzVnR7el/Y3WXdt///1DrDlJ+t4WLlwY\nYl0HT/Ozcvm5OXqs582bF+I777wzxIsXLw6x1iTU/LBm1tbW1mkfxp0lAACAAgyWAAAACjBYAgAA\nKNCvOUv11oupt87S294Wy0hp7QrNR9DHDzvssBCfddZZIda8l+p5Y13/S2uU5Grr6OM6J63z9fp4\nru6Q7i+3PZWrSaWvz61/tsUWW4RY17bTHKbcudc6SrnPhr5/fX2u7lL1+9FtPffccyFevXp1iO+5\n554Qv/HGGyHWz0LuuuntnCZylpqDrh+m/dmPfvSjEH/xi18M8fLly0Osa2l++tOfLtx/Lo+w0ZzU\n3PO/+93vhljXbDzmmGNCXN0naP+QqxunOZXa3+VyVDVnaNKkSSH+m7/5m8L913vNaw7St7/97RDr\nun3bb799iPVv3zXXXBPiZq4FR84SAABANzBYAgAAKMBgCQAAoECfrg3X07kUuTwZfVzniTXvRmtZ\nnHLKKSHWedvW1tYQV9fi0bo8ui+Nc+vY6Rxwbk5Y1wbS7WudDJ2Tz63blDu2Svev+9M5f5XLd6h3\nHUGlx1tzFrTm1pgxY7p8ffUagWab5ivosZ8wYUKIH3roocK21bumFgaH3XbbLcSTJ08Osa5lqeuV\njRs3LsQnn3xy4f70Glu3bl2I58yZE+I99tgjxKeffnpd26/38X/8x38sfFxV90n1rseoOTz1OvXU\nUwv319N/O1944YUQa/+rn40HHnggxPvtt1+Id9lllxA/++yzIW7mHKaNuLMEAABQgMESAABAAQZL\nAAAABfo0Z6mv5y113lnzarR2jq6VpGsj6Zy81i6qzlPS96o1SnJ1gLSukNYdytX90LXYlB6LXM5R\nbi06zUnSnB99fxrr9rSWT66uk+YUad0mlcu50ljfj352quMtt9wyPKY1UrTm1pIlS0Ksr9d8Am1L\nTqP1bAZCvsFgoDmXmlOp9bvOPPPMEGufc8stt4T4wQcfDPHf//3fh1jXE9NrVvNgcp877VM0D3Tm\nzJkhnj59eog///nPh1ivG73OqvtY7a9XrlwZ4u222y7E2v/Um2NU7zXWaE6T5lGedtppIb777rtD\nrOsC6vvVz8ZAxJ0lAACAAgyWAAAACjBYAgAAKNCnOUs9rd45b82D0edvs802IdY5fJ3n1TnwzTff\nvONnrQukc7yag6Tb1jyVXP6Vvj631lxufS99bxprjlFuzr3etd10f3ouNd9g1KhRIdacMK19pMdf\n95+rO6Wvr85ZGjp0aHhMP1d77713iB9//PEQL126NMS5/K9G9fZacugd+hnXa17Xsly2bFmItQ94\n8cUXQ/zLX/4yxPo5njt3bohvvfXWEB988MEhnjZtWojr/dxV969mZq+++mqIf/rTnxbGl19+eYg1\nD6e6T9X+SWtQ5Wj/rfmuWudIr3Hd34IFC0J84IEHhljzwXI+8pGPhFjfr+br6rHTukvanttuu62u\n9jQD7iwBAAAUYLAEAABQgMESAABAgVLlLNVb7yWXJ5PLWRo+fHiIdR5Z5+g1j2Xt2rUhrq5jkst/\n0ljnuLWG0w477BBifS+5nCU9llqzRHOgNB8it1ZdvTlRWgdJ96/t1bon2267bYh1bSatKZOrU6U5\nA7p9ff+as1WdT6E1mF566aUQjxgxIsR77rlniDUfQD/Hum/dn6r3uiKHqTno+oS6HqTmnaxevTrE\n+rnUPMDPfe5zIV60aFGIb7/99hBfcMEFIdbPtfa3Obk1EC+++OIQf+pTnwqxXjdnnHFGiL/61a+G\nuDrHSvsbXc9Rj4XmT6nc2pfaX/7nf/5niGfNmhXiXI5So7XTdt111xDrsX7kkUdC/Nhjj4X4+eef\nD7HmYTajwjtL7r65u//G3Re5+1J3/0rl98Pdfb67L3f3O919m6LtAEB/oA8D0BMKB0sppT+Z2SEp\npb3NbKqZHeLuB5nZeWY2P6U00czursQAUCr0YQB6QjZnKaW0cf7i7WbWYmZrzexoM7u68vurzeyY\nXmkdADSIPgxAo7yGdaOGmNkCM3uXmX03pfQ5d1+bUtq28rib2esbY3lt2Hhu/bFG18fR2kY6D6x5\nKhMmTAixrtl17rnnhljX+NI59er26HvVfWvOUK5ukeYTaM6TzqlrnKvLpDlS+nzNMdIcHz32uXX4\nNJ8gV2dJt6fnVrefO/faXs15mjhxYoj1+Gt+W3UOgbZFa2xp/pPO72u+wm9/+9sQa82Wvl67LZeP\nllIq1WJy3e3D3D3p56RMdttttxDrNarx2WefHWLt//T5mnf5sY99LMSaszllypTC/el6Y7q2XK7P\n0Jws/dxrXSfN0dLt77HHHiGurgOl1+wPf/jDEGv/u8suu4R4//33D7H279r/3nTTTSG+7LLLumxb\nZ7R/+8Y3vhHiI488MsTvete7Qqx/+7797W+HWM/VjjvuGOI5c+aE+L777gux9pdl1tbW1mkflk3w\nTim1m9ne7r61md3h7ofI40kHRQBQFvRhABpVc+mAlNI6M7vFzPYzs1Z3H2Vm5u6jzWxNjdvoThsB\noGHd6cPWr1/f8S/37SwAzae9vT1c513JfRtuxMZvibj7UDObZWYLzewmM5tdedpsM7uxlkb19XQB\ngMGt0T6spaWl459OHQNofkOGDAnXeVdy03Cjzezqypz/EDO7JqV0t7svNLO57n66ma0wsxN7otH1\n3nnKPV//T1A7O62NodvTNbve/e53h1jn7KvnuXVOWnNiNG9FaVs0TyW3Dp7m2KhcnSYdYetAV+f0\n9fHcsc/lmymdU9ccK92fHl9dV2r06NEh1jpWmh+x/fbbh1jfz3bbbddl27Sukj7+1FNPhVgvWF3X\nTj9LuXXuGtXk/5PTp31YX9KcJf2Mjx07NsRXXXVViDXv5pBDwuzkJnkpWqfoe9/7XogPOOCAEOta\ndMuXLw/xvHnzQrzXXnuF+J577gnxP/3TP4VY+6yf//znhe3R607XW6uuO5Xrf/XYLVy4MMRLliwp\n3LeeK61r9J73vCfE2h7t35588snC9l577bUh1rXc7rrrrhBr/td//Md/hFhrcGn/PHLkyBA3U85S\nVwoHSymlxWa2bye/f93MZvZWowCgJ9CHAegJ3FcGAAAowGAJAACgQKnWhlO59W1ya7/V+3pd/2bf\nfePde80VWbFiRYir64Lk1l7L1QHStuocd+5xreuh7dH96bHQOXFtv9ZByr1ez4U+P5fTpNvTHIDc\n2nSaE6XP1xwvrTmjeUHa3ur25GpA6bnRfIdHH320cF96LlWjOUZ8a7U5PPHEEyHeb7/9Qrz77ruH\n+Prrrw+x5pVccsklIT7ttNNCrHl8uj7Zt771rRBrnzRu3LgQ33vvvSHWNRF1PbHp06eHWN+fXsPH\nH398iLX2kV531etJPvfcc+Gx1tbWEOfWitQ4dw3/wz/8Q4hz16Dmg2n+mOar3XLLLSF+7bXXQvy7\n3/0uxJrD+fDDD4f43//930Os+b+6dtxAwJ0lAACAAgyWAAAACpRqGi43baZy02o6daRfL580aVKI\ntfy/3nbWqRj9+nl1KQF9L/rVylzpAG2rTgvpVI+WMdDH9djkptH0vWqpgKKlXsw2/Vqvnovccifa\nfm1P9S1zs03PlR6v3LSePl/fX+6ru9VTbVtttVV4TKcM9Zb+TjvtFGJdNuLXv/51iPXc6TRfbhpO\n30uTlwYYNPS8aUmKPffcM8RPP/10iHWaTa8pfb5eox/84AdDrMvy6NffteSFLiF08sknh1jTHHTq\nSEu3KP36/G233Rbin/70pyHWa7j6/eqx0WtM+zNtu5YW0GOp23vsscdCfMwxxUsV6lIuujSNHjtd\n2kWnyfTY6hSp9reaKqDTcAMRd5YAAAAKMFgCAAAowGAJAACgQKlylhr9ynK9uReal6K5Jvr1cv26\np85TV79e57R1WzoHrHktmpei+9L3mvsqvb5XjfWr9UOHDg2xHptc2QZ9/7nlVHL5Zrmcq1yOk74/\nfVzPrW5f6Vdrq7f3hz/8ITxWvRRKZ6/VZSb0XP72t78NsV4nueuGnKSBSa8RzVm64447Qqx5Ld/5\nzndC/NnPfjbEn/70p0O89957h/j8888P8TXXXBPi3/zmNyHWvEpd/kRLH+hX/3N9xFFHHRVi7cO0\nz9O8oeq8zFypFL2mNEdJaX+n2587d26I58yZE2K9xqdNmxbiBx98MMSao3TYYYeFWEsHHHrooSGe\nP39+iPX9nXvuuSG+4oorQqxLhQ2EPog7SwAAAAUYLAEAABRgsAQAAFCgX3OW6p3HzOVq5Oo06Zy1\n5pasXLkyxOPHjw+xzjvrPHZ1ronmyGjOUW65E6V1mvS9KN2/5gzp/L/WUdL5fl26QOew9XGleTi5\nc6d1j7bddtvC9mrO0ZtvvhliXX5AH3/hhRdCrMdPa3BpzkF1DoDWl9G6Srosgx5rzfXQz6k+X8+l\nyl1nAyGfYDDQ86R5hPoZP/HEE0OsS2Jsv/32IdbcusmTJ4f40ksvDfFxxx0X4nPOOSfEN9xwQ4j1\nGps6dWqIR40aZUW0z9TlS3K13nS5lcMPPzzE1bmEuhSK5leNGTMmxFqzT/tHrWv0zDPPhFjf25e/\n/OUQf+lLXwqx9j9nnnlmiLVP0Hy03XbbLcSaz6b98yc/+cnC9urSXwMRd5YAAAAKMFgCAAAowGAJ\nAACgQL/mLNVbH6beNa9yeTHPPvtsiF955ZUQa62KF198McQvvfRSiKtzBnS+PLdWW66OUm6tNZ1z\nzq1lpPlWmnOka61pjpTWKdKcqty6fPq47k/XbRo7dmyINYdIz8W6detCrPkc+tnQvCDNmdKcK52j\nr/5s6HvTc6m5GXrsNHdO6Wchp941FtEc7r777hBrTs7SpUtDPHPmzBB/+MMfDrGuUfje9743xNqH\nXXTRRSHWz63m0eT6b+0j9JrVemTaB2kfqn3iBRdcEGKtG1WdZ6THStdi0xxJXRtNcyJ33HHHEK9a\ntSrE119/fYi1ZtUXvvCFEOvfD81j1GN34IEHhvjYY48NsR57rVml505rbGn/PRBxZwkAAKAAgyUA\nAIACDJYAAAAKlGptOFVvTlMuR0nneTUPRXNHdA0vzUXRWhrVeTO6L62jpG3N1c4pqunUWds050eP\nhT6uc9Y6J6+P6/Z0jlznsHN1pTTWfATNKdKcJd2f1ozR7enxy7VP89u0DlN1foR+bkaOHGlFfvGL\nX4RYc+caVW89MuouNQc9T5pTqbV0rrrqqhA///zzIf7ABz4Q4tmzZ4dY+4AjjzwyxLk+R+k1qNew\n9sdaa23GjBkhPuGEE0KsOUnah+r2i65T7W80J1RpHaVFixaFePTo0SHWmlHz5s0LsZ670047LcR6\nrPVYnnTSSYXt1XOreZGaj6Z/rwZDH8KdJQAAgAIMlgAAAAowWAIAACjgvVljxd3DxnP1Z1S9bcvl\n5eTyiLS2xoUXXli4v6I8GJ1f1+euXbs2xJqTo/PrSufM9b3rHHQuJ0efn6tZlTs3uWOvdP+5HC/d\nnq6TpXV4DjOcAAAWFElEQVSS9HhpzlHus6ifXc1v2GWXXbrcl57bJ598MsSXXXZZiG+++eYQ67HQ\ntvS0XP6BfpY6ef2ASFhw96SfkzKbMmVKiHU9r2nTpoVYP4faJx122GEhHj58eKNNDPSa12s616fV\nu+Zhbn3M6tpHuvab9mfan+fqDDX6XvSay/19yHnggQdCfOWVV4b40EMPDbHW4Lr22mtDrDW/3njj\njRA3Uw5TW1tbp30Yd5YAAAAKMFgCAAAowGAJAACgQKnrLNWb05SbF9Y5a338/vvvD/HixYtDPHny\n5BBrLkr12kQ6h615J0rfi+a95NYDyx0Lfb3OeefyBTRPRveXezxXN0rPheaK6Pb0+Rrr8cu9P81J\n2HrrrQu3r+tQVb8f3fbLL78c4l//+tch1vn+3LHvbc2UX4D/ozlIN9xwQ4i1ds8Xv/jFEGtdphz9\nXOb6OJXLY9RrXusyaS24ESNGhFiv2Vyez/jx4zt+ztWBq3ctNM05yvWXr7/+eoi1blyj9tlnnxA/\n+OCDIdbab/vtt1+INYdJ891y57YZDbx3BAAA0INqGiy5e4u7L3T3eZV4uLvPd/fl7n6nu2+T2wYA\n9Bf6MACNqPXO0rlmttTMNt4rPM/M5qeUJprZ3ZUYAMqKPgxAt2Vzltx9rJl90MwuNLNPV359tJlt\nXJjnajP7pdXQ2dS7RlUnbcntIsjNoev2dK2ghQsXhljrlGy55ZYhrs570TlqzcHR+fPNN988xJpz\nU29Ok+4/l7OTywnS9uqcu8baHt2+bk/bo4/rudS16LS9w4YNC3EuB0r3X73Wm9mmx1PfX3Veku5L\n11W67rrrQvzmm2+GWPMhcjWwetpAW+epJ/uwMtPPiea9fPSjHw2x1iLT/k/zWPQa1DyW3BqI9dKc\npIsuuijEWutIr9Gjjz46xNrH6vpm1bWV9BrI1anTfKkVK1aEWHOO9JrSukVr1qwJ8ec+97kQax0/\npe9N+0tdu+6II44IsX4WLr300hDrZ6fZ+4ha1HJn6RIz+6yZVV+JI1NKrZWfW82sZ68SAOg59GEA\nGlI4WHL3D5nZmpTSQjPrdOiYNgzBe/d/dQGgG+jDAPSE3DTcdDM72t0/aGabm9k73f0aM2t191Ep\npdXuPtrM1hRupSKlNChu1wEojYb6sOrpFncfkF+JBgaz9vb2mlIbCgdLKaXPm9nnzczcfYaZfSal\n9Lfu/jUzm21mX63898ZaGtXoQCn3hnJrDWmcy7PReeQ99tgjxLp+TnWeTW5tMp0D17o++nydQ37r\nrbcK265ydY/09XquNNa8Gp0T13yD3LnQHCWtYZWrq6TvT1+//fbbh1jzF/R8aD6Eroulax9Vv9/V\nq1eHx/RztGzZshDn8qlyn/tGc4wGWo5StUb7sL6ucdUIPY9PPfVUiPVzqX3K2LFjQ6xrw+XqKmnu\n3YsvvhjiSZMmddbsDi+99FKIn3vuuRDr2nf/8z//E+IZM2aEWNc70zwgzWus7s8feuih8NioUaNC\nfPjhh4d45513DvHUqVOtiPY3xx9/fIi/8Y1vhFj/HuR86UtfCrH2V2eddVaI9VyecMIJIda1Ny+5\n5JIQN3OfoX+L9G9Lx/Pq3O7Gq/HfzGyWuy83sw9UYgAoO/owAHWruYJ3SuleM7u38vPrZjaztxoF\nAD2NPgxAdzEBDwAAUKBP14bLzWs2uhZcvbkdue1pnZKf/OQnIdY6I9V5NzonrWvn7LDDDiHWnJhX\nX301xLo9zTfQ96JrneVylDRvRveXW7tN123K1ZXSuiRas0pzsnI5WrmcJj23uj/NGdPzpfvTnKzq\n9d3uu+++8JiuBZdbl6/eNbYazRdo5nwD/B89j5rncvPNN4dY1zfU2mJay0evCc2B+spXvhJirbv0\nta99LcRaC2jBggUhXrJkSYjPOOOMEGutIM3L+fjHPx7iX/3qVyHWPKTq9hxwwAHhMV0376CDDgpx\nbq1MPRd6rjRHUv8+qNzfri9/+cshvuCCC0L8ne98J8SzZs0KseZM6bkeyHmOXeHOEgAAQAEGSwAA\nAAUYLAEAABTo05yl3p7nrHe9M82r0XlmzcO54447Qqy1Nc4888wu96V1fcaPHx9ibbvuW+fjte6F\ntj2Xg6T709fncpy0jpHmSGlejsZa90jrNun+9VxqDlO9tYk0X0PzM/T9ag7ZbbfdFuKrr7664+dn\nn302PJY7Fl3V9QB6ktYx0lw6vYZvv/32EH/9618P8RVXXBFizWk67rjjQqw5SmrHHXcM8TXXXBPi\n1tbWEH/mM58JseaU3nvvvSHWOk+PPfZYiOfMmdPxsx4L7b9z/YvmKOVqAOrfFu2fdP96LjVP8sMf\n/nCIJ0+eHGI9NkrzdTVnaTDkKCnuLAEAABRgsAQAAFCAwRIAAECBPs1ZalS9dZhyi15qPRvNJdE8\nGnX55ZeHeLvttuv4WWuC6LY0Z0fnoLXuhuYUac6SzmlrzpJur958Lo01pye3VpzWMdI6R3putc5I\nLgdL1+LTulX6fI2V5kTNnz8/xHruq9fBytVRyq1rV+/nHOiOBx98MMS6llt1DqaZ2cMPPxziRx55\nJMSnnnpqiHW9s1WrVoX4sssuC7HW+tH2aM7SihUrQjx9+vQQL1y4MMQHH3xwiLX91dvXukrav2md\nIs0JGj16dIh1LTXN/7rrrrtCrDlHuv/qHEkzs9deey3E06ZNC7Eeq+XLl4f4iSeeCHGuLtRgxJ0l\nAACAAgyWAAAACjBYAgAAKFCqteF6evu53A7Nu9EcJs0D0jwXrcXx1a9+teNnrcujNUGGDRsW4nXr\n1oVYc360LZoDpDlMunacHps333wzxLpOktI6T9p+bZ8eW60bpe8vlyOl51Ln8PX5enz0+boO1po1\na0L83e9+N8Tf//73Q6zvrzpnSh9TjeYc1XsdNbqmIgYHrQ925ZVXhljzLGfOnBni97///SH+whe+\nEGLNm9S8QP1cau7fo48+GmK9xhctWhTiCy+8sPD1Z599dohffvnljp8nTpwYHtNrWmtCXXfddSHe\nfffdQ7xs2bIQa/+56667hlhzhrSu3QsvvBBiPZaaE3X//feHWGtiaf+JTXFnCQAAoACDJQAAgAIM\nlgAAAAr0a52lRuvJ5Ooq1Vu/Ruvf/PGPfyx8fdH2tYaI1iipXofIbNM57rFjx4b4+eefD7Hm4Oja\nZppflcvHGjNmTIg150nnzHXOXXOANAdL6yBpDpO2R9eR0rWSRo4cGWKd41cLFiwIsa6lNG/evBBr\njoEeb30/1fRY59bd0+fnPqf15hzp83OxYu26wUFzhPQa0DpJ+rn95je/Wfh87RNGjBgRYs0B1fUz\ntQ+aOnVqiPWa1vXMdHsTJkzocnuaU3nVVVeF+Nhjjw2x1jm6+OKLQ3zkkUcWbk9rRN16660h1v5R\n87M0h0mPlSJHqX7cWQIAACjAYAkAAKAAgyUAAIAC3ps1Vtw9bFxzNTSXo97ci3rrxzSa46Tt1Xnk\n6jn8XJ6K1jXSmiW6TpLWFNG6HDoHrXWRNNa113StNs15qq5BYrZpvoLmLOlad5qPoDVXtC6V1g3R\nWOue6NpHK1euDPEdd9wRYp3T13wNzZHSHIaiHCnN8dHPQr16ui5TvdvT5+fW1UspDYiFpNw9kdvR\nNb2mdT3G6rUyzcymTJkSYl0/TfMgb7/99hDvv//+Id5vv/1CrHWVzjrrrBAfd9xxIdY6Uu9973s7\nftY6RK+//nqIb7755hBr/tTcuXND/Morr4RY+wS9RjW/S/tbPfbUSus5bW1tnfZh3FkCAAAowGAJ\nAACgAIMlAACAAv2as5Sr91JvnaRO9h/iXL2YRteuq56H1ryOo446KsS6FtlDDz0U4p133jnEWhPk\noIMOKnxc13HSWOsGaU0Ufb7WRVq4cGGIH3vssRBrXSjN8dG16bSuk9ZIya3jl6P5ELk8nFy+m36W\nqj879X6O9Loo2nZn6s1JanStOHKWUAvNy9E8y9w11traGuLjjz++8PnXX399iE8//fTC52uttOo+\nWHMeNWdTc4pUo3mK6D/kLAEAAHQDgyUAAIACDJYAAAAK9GvOUq7OUbOpfj9at0dzH3TdOV2bTXOK\ndO0ffb3SmidaB0mPveZQaZ0hfT9aRym3Nlvu3Ov2NY9Gn685TJrno9urNwep3ryhIvW+ttFcvXpz\nkurNWcrli5GzhM7UW0dPc+O01pzmOWqekK4vOWnSpMLXP/300122BYNHVzlLNS2k6+4rzOz3Zrbe\nzNpSSge4+3Az+4mZ7WRmK8zsxJTSG11uBAD6Af0XgEbVOg2XzOzglNI+KaUDKr87z8zmp5Qmmtnd\nlRgAyob+C0BD6slZ0ttSR5vZ1ZWfrzazY3qkRQDQ8+i/AHRbTTlL7v6sma2zDbexL0spXeHua1NK\n21YedzN7fWNc9bq66iwNZJoDo7kQRevMmW2ao6THTnN0crVwNCdKt6+P67nT9uXW+cvlIyjNSWq0\n1lCZPmuNrtXW13L5XZ08vzwH2xrrv8hZGjjqzZnC4NRQzpKZHZhSetndtzez+e6+rPrBlFLSgREA\nlAT9F4CG1DQNl1J6ufLfV8zsBjM7wMxa3X2UmZm7jzazNV1voWM73W8pAHRDI/3X+vXrO/7l7qgB\naD7t7e3hOu9KdrDk7lu4+1aVn4eZ2WFmttjMbjKz2ZWnzTazG2vYVg1NB4Ce0Wj/1dLS0vFPp6IB\nNL8hQ4aE67wr2Zwld59gG/5vzGzDtN0PU0pfqXz1dq6ZjbcuvnpLzlLtGp1Pz9Ux0ljrImnOlO7v\nz3/+c4hzax/p+9ERe73r9qH/NHPOUqP9FzlLwODSVc5SqRbSHcwYLDFYKqtmHiw1gsESMPiwkC4A\nAEA3MFgCAAAoUGvpAPSyRqckdWokN+2lcW6tOX1+vdtXfDMSANAsuLMEAABQgMESAABAAQZLAAAA\nBfo0Z6nM63UNNPXmDNV7Lup9PjlKAIBmxZ0lAACAAgyWAAAACjBYAgAAKECdJfQKcpQAAAMFd5YA\nAAAKMFgCAAAowGAJAACgAIMlAACAAgyWAAAACjBYAgAAKMBgCQAAoAB1lgAUomYWgMGOO0sAAAAF\nGCwBAAAUYLAEAABQgMESAABAAQZLAAAABRgsAQAAFGCwBAAAUIDBEgAAQAEGSwAAAAUYLAEAABRg\nsAQAAFCgqdeG0zWr3L2fWgIAAAaqmu4sufs27n6duz/p7kvdfZq7D3f3+e6+3N3vdPdteruxAFAv\n+i8Ajap1Gu5bZnZrSmmymU01s2Vmdp6ZzU8pTTSzuysxAJQN/ReAhrhOZW3yBPetzWxhSmkX+f0y\nM5uRUmp191Fm9suU0iR5TpI4bHvIkMZSppiGA3pfe3t7iHN9RkqpNBdio/3XZptt1oetBdDf2tra\nOu3DaslZmmBmr7j7lWa2l5n91sw+ZWYjU0qtlee0mtnIRhvZ6OCnhoFf3W3qKT393nr6vdR77Mo2\nUO3r41N0PPr6WNT73st27npZn/VfAAauWm7tvM3M9jWzS1NK+5rZWya3rNOG3rf4ry0A9D36LwAN\nq2WwtMrMVqWUHqnE19mGzmd15fa1uftoM1uT21Du7gUA9LCG+q/169d3/NPpSADNr729PVznXckO\nllJKq81spbtPrPxqppktMbN5Zja78rvZZnZjblsD/HY/gJJptP9qaWnp+NdojiWA8hkyZEi4zruS\nTfA2M3P3vczs+2b2djN7xsxOM7MWM5trZuPNbIWZnZhSekNeV1eC92C+85TLCcLAUW9OUX+r945K\nmRK8zRrrv0jwBgaXrhK8axosdReDpdoxWBo8GCw1BwZLwODT1WCJ+8oAAAAFGCwBAAAU6NO14XR6\noSjzfLDr7Wk5pv0AAKgNd5YAAAAKMFgCAAAowGAJAACgAIOlkuutXKKUUvgHoHNlr9xd5vaVuW1m\ntK8RZW6bWc+3j8ESABQo+/9MlLl9ZW6bGe1rRJnbZtbz7WOwBAAAUIDBEgAAQIE+Xe4EwMA3kJY7\n6e82AOh7fb42HAAAQLNjGg4AAKAAgyUAAIACvTpYcvcj3H2Zu//O3f+1N/dVY3t+4O6t7r646nfD\n3X2+uy939zvdfZt+ats4d7/H3Ze4+xPufk7J2re5u//G3Re5+1J3/0qZ2ldpS4u7L3T3eSVs2wp3\nf7zSvodL2L5t3P06d3+ycn6nlal9/aVMfViZ+69KW0rbhzVD/1VpD31Y99rW6/1Xrw2W3L3FzL5t\nZkeY2e5mdoq7T+6t/dXoykp7qp1nZvNTShPN7O5K3B/azOyfU0pTzOw9ZnZW5XiVon0ppT+Z2SEp\npb3NbKqZHeLuB5WlfRXnmtlSM9uYiFemtiUzOziltE9K6YDK78rUvm+Z2a0ppcm24fwuK1n7+lwJ\n+7Ay919mJe7DmqT/MqMP667e77+0knNP/TOz95rZ7VXxeWZ2Xm/tr4527Wxmi6viZWY2svLzKDNb\n1t9trLTlRjObWcb2mdkWZvaImU0pS/vMbKyZ3WVmh5jZvLKdWzN7zsy2k9+Von1mtrWZPdvJ70vR\nvn48Z6Xrw5ql/6q0p5R9WBn7r8r+6cO6164+6b96cxpujJmtrIpXVX5XNiNTSq2Vn1vNbGR/NsbM\nzN13NrN9zOw3VqL2ufsQd19Uacc9KaUlJWrfJWb2WTOrrnFflraZbfi/srvc/VF3/2Tld2Vp3wQz\ne8Xdr3T3Be5+hbsPK1H7+ksz9GGlPEdl7MNK3n+Z0Yd1V5/0X705WGq6mgRpwxC0X9vt7lua2fVm\ndm5K6c3qx/q7fSml9rThNvZYM3u/ux8ij/dL+9z9Q2a2JqW00Mw6rfHT38fOzA5MKe1jZkfahumJ\n91U/2M/te5uZ7Wtml6aU9jWzt0xuWZfg+PWHpnq/ZTlHZe3Dytp/mdGHNahP+q/eHCy9aGbjquJx\ntuH/zMqm1d1HmZm5+2gzW9NfDXH3zWxDJ3NNSunGsrVvo5TSOjO7xcz2s3K0b7qZHe3uz5nZj8zs\nA+5+TUnaZmZmKaWXK/99xcxuMLMDStS+VWa2KqX0SCW+zjZ0PqtL0r7+0gx9WFk+Q1ZpQ+n7sBL2\nX2b0YY3ok/6rNwdLj5rZbu6+s7u/3cxOMrObenF/3XWTmc2u/DzbNsyz9zl3dzP7LzNbmlL6ZtVD\nZWnfiI3fJnD3oWY2y8wWlqF9KaXPp5TGpZQmmNnJZvaLlNLflqFtZmbuvoW7b1X5eZiZHWZmi8vS\nvpTSajNb6e4TK7+aaWZLzGyelaB9/agZ+rBSfIbMyt2Hlbn/MqMPa0Sf9V+9nHh1pJk9ZWZPm9n5\nvbmvGtvzIzN7ycz+YhtyEU4zs+G2IaluuZndaWbb9FPbDrINc9WLbMNFvNA2fPOlLO3b08wWVNr3\nuJl9tvL7UrSvqp0zzOymMrXNNsypL6r8e2LjtVCW9lXaspdtSHp9zMx+ZhuSJkvTvn48LqXpw8rc\nf1XaV9o+rFn6r0qb6MPqb1+v918sdwIAAFCACt4AAAAFGCwBAAAUYLAEAABQgMESAABAAQZLAAAA\nBRgsAQAAFGCwBAAAUIDBEgAAQIH/DwQG9wDTMW6AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1082b31d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 10))\n",
    "axes[0].imshow(target.get_data()[:, :, 15])\n",
    "axes[1].imshow(source_after.get_data()[:, :, 15])"
   ]
  }
 ],
 "metadata": {
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}