{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5f4bf0ce",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Linear vs Circular, Normalized vs non-normalized correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "409bc3c0",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# we should test linear vs circular to understand the differences\n",
    "# the main difference in the new version was a bug as \n",
    "# the zero-padding for linear case does not work as long as you leave values on the \n",
    "# borders not zero. \n",
    "# see new normalize_intensity function that should take care of uneven illumination\n",
    "# both between the frames and at different regions of the image A or B\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8dfee62b",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage.registration import phase_cross_correlation\n",
    "from skimage import feature\n",
    "# from skimage.feature.phase_cross_correlation import _upsampled_dft\n",
    "from scipy.ndimage import fourier_shift\n",
    "from skimage.exposure import match_histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6a9bfb1e",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from openpiv.tools import imread\n",
    "from pylab import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e97f0258",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.fft import rfft2, rfftn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a7d672f4",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from openpiv.pyprocess import find_subpixel_peak_position\n",
    "from openpiv.pyprocess import normalize_intensity, fft_correlate_images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ba08daf9",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "a = imread('../test11/A001_1.tif')\n",
    "# b = imread('../data/PIVChallenge2001_A/A001_2.tif')\n",
    "a = a[:32,:32].copy()\n",
    "a[16:18,16:18] = 255\n",
    "# b = b[:32,:32]\n",
    "\n",
    "# should be in the order of y,x:\n",
    "# so it's about 5 pixels upwards and about 3 pixels to the right\n",
    "shift = (-12.035, -10.92)\n",
    "# The shift corresponds to the pixel offset relative to the reference image\n",
    "b = fourier_shift(np.fft.fftn(a), shift)\n",
    "b = np.fft.ifftn(b).real\n",
    "b = match_histograms(b,a).astype('uint8')\n",
    "# b = b + np.linspace(10,85,32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6cb417a6",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fbafe912a90>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAD5CAYAAABPqQIFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAArAElEQVR4nO2de3Rc1XXwf1tvS7IlWzJ+4QcYQ7AJGGweLZAQ0jQOIYW0hECzyIMsHBpYC75F0wBpG1rKatI0pHxNIHECgaSUR4AQSknAARpwv2CwjbENxsQmtrEtPyTbsiVbz9nfH/eqHVlnX81IY0nX3r+17tLM2XPuOffO1Z5zzt5nb1FVHMdx0krRcHfAcRxnMLgScxwn1bgScxwn1bgScxwn1bgScxwn1bgScxwn1ZQMprKILADuAoqBH6nqN5I+X1xdpSV144Kyona7XlF3uDxTnNC3jC3TBNU9Z8IuU/Zma7jvetDuSFGH3Ra1XaZoduUeU7Zm13hTphXGhWfErFO6z5Z1jjZFFJUaXwwgzeF7kim1z4fdjcTvM5PwFIvRRUnwLEr6zpKene5q+6TSYV+clobrSVfCDTGa6ty7m+7W1oSK/fPRD1Vp0277u81m+ar2Z1V1wWDaGywDVmIiUgx8D/gIsAV4TUSeUtW3zMbqxjHpqzcEZaM32IqgrDn8jXXU2N9V0oOYKbNlr37lblN2yiufCZZ3r6ox61Rusx9svXi33Y/5j5iy2d/7silrP/lgsDzTan/Vk16w7/328+3+V05uMWXlvxoTLD84Pumf2RRRcsCWtdXbfSxtDrdX3Gmfr2qbrTE7K+3+7z7HPmnFZvuha5sUrlfaZH9nxYZS3PSDO806udK4u5ulzx6b02dLJ22oH3SDg2QwI7GzgPWq+i6AiDwMXAKYSsxxnDSgdGvC0HeEMZg1sSnAe1nvt8RljuOkGAUyaE5Hf4jIVBF5UUTeEpE3ReSGuPw2EdkqIivj46KsOreIyHoRWSciH+2vjUGtieWCiCwEFgIUj6s93M05jlMAMhRsJNYF3KSqK0RkNLBcRBbHsu+o6j9nf1hEZgNXAHOAycCvReREVTUX6QajxLYCU7PeHxuX9UJVFwGLAMqnT/WNmo4zwlGUzgJNJ1W1AWiIX+8XkbUkz9guAR5W1Xbg9yKynmjp6rdWhcFMJ18DZonIcSJSRqQ9nxrE+RzHGQEo0I3mdAD1IrIs61honVdEZgCnA0vjoutFZJWI3CciY+OyvJepBjwSU9UuEbkeeJbIxeI+VX0zqU7pfpj4ctiqkim1Nf/J14dPu/Kn7zfrtNfa/Tg4zbYizbvtL0yZ1Ib7XpZg6dpzqn1dY58Ou2wAnLzCtkC2JfS/dHNFWDCtzazTXmM/BnUrbGtc1S8rTdnmj4Wvu6LBtoRW2N4tlBy0B/F1b9quKpv+JFw+5m3bFNp06sA8FD5x6hum7MW3zzRltavCffnRX/6LWeeyF8PPh5YVZrKTy3pXTKOqzu/vQyJSDTwO3Kiq+0TkHuB2Ip15O/Bt4OqB9HVQa2Kq+gzwzGDO4TjOyEKB7gKG6BKRUiIF9qCqPgGgqjuy5D8Eno7f5rRMlY177DuO04dMjkd/iIgA9wJrVfXOrPJJWR/7JLAmfv0UcIWIlIvIccAs4NWkNg67ddJxnHSh/7veVQjOBa4CVovIyrjsVuBKEZlLNPDbCHwJQFXfFJFHifxNu4DrkiyT4ErMcZxDUIXOAukwVV1CeEOZuQylqncAd+Tahisxx3EOQehO2sg6wnAl5jhOLxTIpMijc0iVWNco2HNy/raE/3rjfcHy+gSze1dVwubwA7aZP1Nq12uvDbc3apddp2Sffb1Jm6GTNjxPmGpHuGhuOCZYrtsM1wug3Nhg3x9Nc2w3hboV4XM2nmfvzC9fbm+S3n2h7SLy2oU/NmUnLwq7Iux7v92PY//Dfj62LLCXs1981HajSMJyH/n8PTeadcae1xQs31WSW/SJ/vCRmOM4qSVydnUl5jhOSlGgMylw2gjDlZjjOL1QhO4UuZC6EnMcpw8Z9emk4zgpxdfEEpBuKGsOy/786sVhAfBvP/1IsLy0dWDhQkpa7aFy0sZxi6Rwxm985fumbOYj15qyeWf+zpQtf22WKas5I2y1Ou2YbWadl0rnmLLqzfa9+tDlr5myl38cttT9fsGPzDq3nn6qKfv5L84zZbPX2pvlK/aGy8uabUto42mmiFGbbctl6zTbMpgUfr1lWlhhJFmoyx4dGyyX3YX4lxa6fU3McZy0EkV2dSXmOE5KURU6NCGV2AjDlZjjOH3I+JqY4zhpJVrY9+mk4zipxRf2HcdJMb6wn4Bk7Mzc9/7HH5n1Jr8Zjim/f4rd/cpdttvD6KveM2Xrd9gJjUtWVwfLD9bb6wen32Gb/8M5siNWN59kyqoNNxWAloPhuP0rdtWZdcYn3Ktp164zZU+/coYpK5kS3tT8TmerWefJx2w3io56u49lzfY/3NyrVgfL31e13azzq+2zTdmW1yebsuIDdj8yCdnNrSzxi265y6yz8B9vCJYXagDV7c6ujuOkFUXo1PSohvT01HGcIcEX9h3HSTWK+HTScZx0k6aF/fT01HGcIUEVurUop6M/RGSqiLwoIm+JyJsickNc/i0ReTvOAP5zEamNy2eIyEERWRkf9ubjGB+JOY7Ti2hhv2DbjrqAm1R1hYiMBpaLyGJgMXCLqnaJyDeBW4CvxnU2qOrcXBsYlBITkY3AfqAb6OovnXlxh1KzMewuMarRvmk754bt0xVnhyM2AOzpsi+t7afTTVn1pbvteoRdLFrPbzHrlK8I1wE4cOpBU7bwtCWm7Ee/tN1RSg6E1zIOjrfj6LfMazdlu7dMMWXTfpkURSQs+7O9XzFrtE+yz1fZYP/qt9Xb17b2u+EIHWvKTjHrdNTY60G1CfkIqq+0I4Xs+C/7Plo5Dv5m46V51ykqTIj9gi3sq2oD0BC/3i8ia4Epqvpc1sdeAS4baBuFGIl9SFUbC3Aex3FGAIrkExSxXkSWZb1fpKqLQh8UkRnA6cDSQ0RXA49kvT9ORF4H9gF/raovJ3XAp5OO4/Qhj5FYY38zMAARqQYeB25U1X1Z5V8jmnI+GBc1ANNUtUlE5gFPisic7DqHMtgxowLPichyEVk4yHM5jjMCiPJOFuV05IKIlBIpsAdV9Yms8s8DFwOfUVUFUNV2VW2KXy8HNgAnJp1/sCOx81R1q4gcAywWkbdV9aVDLmAhsBCgfFTtIJtzHOfwU7gM4CIiwL3AWlW9M6t8AfBXwAdV9UBW+Xhgt6p2i8jxwCzg3aQ2BqXEVHVr/HeniPwcOAt46ZDPLAIWAYyuPTZFeYUd5+gkStlWMOvkucBVwGoRWRmX3Qr8X6CcaPAD8IqqXgt8APh7Eekksg5dq6q2tY1BKDERqQKKYotDFfDHwN8P9HyO44wMVCXnqWL/59IlEBzWPWN8/nGiqWfODGYkNgH4eaxFS4B/V9VfJVXoGAub/iw8GBuzOv+b1nKg3JSNe7rSlHXbOSJoW2pHe+isCfe9+mXbjSKT0FbNGDsTxAMPh5OjANRtsl0RWidb0wB7ejBu3H5T1vnkeFO240z7nBW7TJFJUjSKiUttd5Smm+z7eLAh/H3un2n7IlQ02KOQzkr7mjuesd0osB9HDtaHr/vgYzPMOu0nhMu7fmO3kw9HRTwxVX0XSMgL4zhOGoniifneScdxUotHdnUcJ8VELhY+EnMcJ6UUeO/kYceVmOM4fUhTKB5XYo7j9CIKxePTySBFbUL122GfgwOTbD/Yib8NuxS07LVdG7rL7PO1JST2KN9riqh7K9yPHZ+2zf9WchGAj09da8omzbSzgdz9s4+bMivpROs02y1j1zo7OcrkhCQi5c32fdw1LyyrbLDrHPfBjfb5NtuRR9qWjjJl7cZ1162wRxrFHfY1J7lYJFHUYdcr6gx/Z11Vdp1Ru4woFuEgMXnja2KO46SWKIqFTycdx0kp0bYjV2KO46QWH4k5jpNy3GPfcZzU4tbJBFSg20jnPma9Xa9lctjxrnWKbYFMSmA85Te2CaejOmHzr2EtqnmuyqzTdIa90fixdXNN2ailtlVz4kVbTdluYxPy2FX29KCoyxTRMjkhtr29N5xoZaUvP7vxW2aNT/+THX9/3/uT4vnnH+EpKY5+lWHhBahMsNae/7e/NWVPPnaeKatdH25v78l2Pya8Ei4vVIx9n046jpNa8oyxP+y4EnMcpxcKdPlIzHGcNOPTScdx0ov6dNJxnBTjQREdx0k9PhIzKO6Eqoaw2Xj8VZvMejseCm/+7Zxou0pMeta+tKY5hp8HULUtyZQfJpNwF5NcGzKlthtFyUHbvN642I7lXtaav7m+dm3+m5MB1i68x5TN/t6Xg+UXP/iXZp2ahGset9ruY9coW2bdxzZ7z7sZ8x6gvNl+Ph56/SxTVmE/cqbrTvVm+7r2HsYY+4UMiigiU4GfEOXkUKIM4XeJyDiirN8zgI3A5aq6J07xdhdwEXAA+LyqrkhqIz2rd47jDAmK0JUpyunIgS7gJlWdDZwDXCcis4GbgedVdRbwfPwe4GNEuSZnEeWrtX8pY1yJOY7ThwyS09EfqtrQM5JS1f3AWmAKcAnwQPyxB4BL49eXAD/RiFeAWhGZlNSGr4k5jtMbPTxrYiIyAzgdWApMUNWGWLSdaLoJkYJ7L6valrisAQNXYo7j9CLPNbF6EVmW9X6Rqi469EMiUk2UFPdGVd0X56uN2lNVEcl/71iMKzHHcfqQhxJrVNX5SR8QkVIiBfagqj4RF+8QkUmq2hBPF3fG5VuBqVnVj43LTHxNzHGcXihCd6Yop6M/YmvjvcBaVb0zS/QU8Ln49eeAX2SVf1YizgGas6adQfodiYnIfcDFwE5VPSUuC5pH+zuXYrsjWG4UABOuDLtfZBLSvO/7tB2jPrO8xpTtON8OAzBqc9hOrgnm85IDtqzqgp2mbPdKO0TEsS+0m7KNl4ZvcEWDHZ2jyw7CwaiEqA0zH7nWlFVa3i+t9i9896VNpqxzSZ0pK9trz0QOXrwvWJ70DBSF00D0y4SJe03Zvg3HmLLjrn4nWL72qZPMOrXrw99LSZtZJS8K6Ox6LnAVsFpEVsZltwLfAB4VkS8Cm4DLY9kzRO4V64lcLL7QXwO5TCfvB75L5OvRQ4959BsicnP8/qs5nMtxnBGOFnBhX1WXgKkRPxz4vALX5dNGv+NBVX0J2H1IsWUedRznCEBVcjpGAgNd2LfMo47jpJ6jbAN4f+ZREVlI5HlLafXYwTbnOM4QMFJGWbkwUOvkjh4v2kPMo31Q1UWqOl9V55dUJKwgO44zIlCF7ozkdIwEBqrELPOo4zhHAIXadjQU5OJi8RBwAZFn7hbg69jm0USKuuxEC1YyEICNv54RLO9KSBRSvMo2oVcbkTQibH8JK3X8wfH2l/mJK5eYsp8tPteUTXzNdm1oq7P7OHpD+Hdp/0zbdeTs+W+ZsrXfnWPKqjfn/xvYMtOOPFKc4EbRXmt/Z0kRKaqfDz8HSdEo9k23r6usLCGaxn/abhRd57eYss3fnxUsr9tj36ui9nD/pWvAju//g5Ku6WS/SkxVrzREfcyjjuMcCRxlC/uO4xx56OAHdEOGKzHHcfpwRE0nHcc5uoisk+nZVu1KzHGcPvh00nGcVOPTyQEwqtE2eZcbASl2nmP/XLzzpz8wZZe/axtWV//ajhww87PhaAPLN04z6zx7j+1GUZrgmtGdYMrfc7I91O844WCw/N0Lf2zWuWjdRaasyLby02F7sdBRE/4+pz9l12mdYH+fRZ32/UiKtNFZabdnkRSVpLjD7mOpkaQFoGWz7ehtJWPpLhvAlK4As0Bl5OyLzIURo8Qcxxk5pGg26UrMcZxDUNARsqUoF1yJOY7TB59OOo6Tatw66ThOajni9k46jnOUoYArMQOxXQeS3AaOfTHsNrBv+iizzrzb/sKUZUoTviA7Pwcrt0wJlpevtftRlBBVIMmU33ia3ccJr9kRKbZMC0cDSUrqkanpMmWlCf2oWW9fW/MJ4XrlTXaSk/1T7MexqMMU0XyC/eysXXh3sPz8675k1umqHJgbxYHxdj/GrDdF7DX6335y+LkHmPx4OJuJSmGUj08nHcdJMVIw66SRLe0RoMchsxbYq6pz4wzha4F1sewVVbV/fWNciTmO05fCjcTu55Bsaar66Z7XIvJtINudfYOqzs2nAVdijuP0Rgu3sK+qL8UjrD7EiXUvBy4cTBvp2aruOM7QoTkeg+N8YIeq/i6r7DgReV1EfiMi5+dyEh+JOY4TIOeRWL2ILMt6v0hVF+VY90rgoaz3DcA0VW0SkXnAkyIyR1XDadxjhlSJZYqhvSZ8c6wNwwAbP14RLK/eZLdVZBvcmHDlRlP2zEnPmDLL4plkgUyi/g3bypi0Abxptp2PYPrj4ftYsW2/fb7Tx5iy5hNMUaKl7vifhdtrnVFt1hm3rs2U7Z9absqSNoefvOjLwfIJHfbO9rpVhXcvSHoeLStk8ebwcw+wd1a4j93/L69u2dj/jofSqKrz8z29iJQAfwrM6ylT1XagPX69XEQ2ACcCy4InifHppOM4venxE8vlGDh/BLytqlt6CkRkvIgUx6+PB2YB7/Z3IldijuP0QTW3oz/ibGm/BU4SkS1xhjSAK+g9lQT4ALBKRFYCjwHXquru/trwNTHHcfpSIBcLK1uaqn4+UPY48Hi+bbgScxynL77tyHGcNCO+7chxnNSiAkdSUERj79NtwDXArvhjt6qq7ZsQo8XQVRW+Od+++KdmvReaTw6Wv/zjM806bfX2l/DeL2eYsuPeucaUTWoO/zwlxV3fN8N2h9j9fruP45ebIuresl0z2saG2xu1ya6TSXgK6t+wr611sm0X2nHO6GB50nU1zbY30o/ZZNv89yW4gUz5TdiVommOvfu+bG/Sxnb7mpNcZhpPs5+DTEdYVtaRkIPB6H7BZoEpGonlYp28H1gQKP+Oqs6Nj34VmOM4KWJoPPYLQr9KTFVfAvo1czqOcwRxJCmxBK4XkVUicp+IjC1YjxzHGV6Gxtm1YAxUid0DzATmEu13+rb1QRFZKCLLRGRZ94HWATbnOM5QIprbMRIYkBJT1R2q2q2qGeCHwFkJn12kqvNVdX5xpZ1A1HGcEcSRPp0UkUlZbz8JrClMdxzHGQmkaSSWi4vFQ8AFRCE3tgBfBy4QkblEungjYAcszyJTBi3TwqbyW+//rFnPjEWfEA+/Zr1tkk+KvlD3pn3O0pZwfPit59muAaN22W2NWW//hnSX2fUOJsRyr9oWvu6WWTVmnZp37QD2L/zkXlOWhBXxY8f5djiHSS/Y59s1z15/qWi0ZVs/GH7ExyTkB6jcleTOkeBWcmbC97LVFFFyIBwvP8nVw4qKUWwH58iPEbLelQv9KjFj79PAnmzHcUY+I2iqmAvuse84Tl9ciTmOk2Yk96CIw44rMcdx+uIjMcdx0spIsjzmgisxx3H6ciRZJwtJ8UGoXWskCqm167XMDLsA/H7Bj8w6xz9he31Med5uq6PajjZgRYgoPWCfLylBRFeC72/tGjtpxpYP2S4dJa3h+ztuXdg95HBhuSlU/sp2QyjusKNAVDTaUSemLG42ZR114WQbnaPtRz8pEUtXve2OkvQ8znzETmRtudocHG8rkgkXhH02il4ukI+Fj8Qcx0kzPp10HCe9qFsnHcdJOykaiXnKNsdx+lKgDeBxqK6dIrImq+w2EdkqIivj46Is2S0isl5E1onIR3Ppqisxx3H6UMAN4PeTY2RoEZlNlI9yTlzn7p5kukm4EnMc57CRZ2ToS4CHVbVdVX8PrCchzFcPQ7omJhkoPRBW31Y5QNnesHl99oYvm3XqGuzzdVYlRDDYbpuo95xUHixvmXfQrPPouT8wZQv/8QZT1jYuHNkAkhNSbPlEWDZunVmF7lL7t2zBxz9j1xtj97HllPAPaPU2u+9F7fZq8vgVtotI0+ljTFnVjrCPy45P29/ZuKcrTVntevtfZmaz7UYxavp+U3b++e8Ey9fedIpZZ/F1vwiWn1Wx16yTF7mvidWLyLKs94tUdVEO9a4Xkc8Cy4CbVHUPMAV4JeszW+KyRHxh33Gc3uRnnWxU1fl5tnAPcHvUErcTRYa+Os9z/A8+nXQcpy+HMbJrQmTorcDUrI8eG5cl4krMcZxeCIc3smtCZOingCtEpFxEjgNmAa/2dz6fTjqO05cC+YnlExlaVd8UkUeBt4Au4DpVtRdRY1yJOY7TmwJGscg3MrSq3gHckU8bQ6rEtAg6K8ObWsd8aptZr+n1ycHy2rftO/3aP9xjyv7w/9hWpEy5PcPOGMa4RX/wE7POZc9dZ8pqS+0Nvtv/wO5H7dumiJLGcL3ifbY1rnOybY3bdIkdm3/SEttiWL8m3N6GT9kWzarNYesv2LkDACr22LKmOWHL9jsf+KFZZ9bmcH4AgLpV9jNXM3OPKSt+cpwp+++S8Lp4+US7LesZfvu9fzHr5IVvO3IcJ834BnDHcdKNKzHHcVKLZztyHCft+HTScZx040rMcZw0c0QFRRSRqcBPgAlE+nmRqt4lIuOAR4AZRA5rl8ebOE20BDpqw24Fe5+093nWHgz/LDSfYLd15l/bZvKM7TVAw4W27MqzlgTLb75joVlnTI3tRrHvBPtJqd5su1i0TDNFSFf4Xu29w3aHKHpgtH3CBPZPs90lRm8Ox6Ifuyppk4j9899dZt/HJFndm/nHnH/4U3eZspt+Y7vMdD5tu1E0nm0nWyhpDv8bHv+Enbxh/3TbLWbQpGxNLJdtR11Eu8xnA+cA18Vxf24GnlfVWcDz8XvHcVKO5HGMBPpVYqraoKor4tf7gbVE4TEuAR6IP/YAcOlh6qPjOEPNYdwAXmjyWhMTkRnA6cBSYIKqNsSi7UTTTcdxjgDSZJ3MOYqFiFQDjwM3quq+bJmqmnpZRBaKyDIRWdZ9oHVQnXUcZ4hI0UgsJyUmIqVECuxBVX0iLt7RE1Ij/rszVFdVF6nqfFWdX1yZkC3WcZyRQRwUMZdjJNCvEhMRIdp1vlZV78wSPQV8Ln79OSAcL9dxnPSRopFYLmti5wJXAatFZGVcdivwDeBREfkisAm4vL8TSTeUGDPKmo22KXzrB8PdTIrmsPd9tmzMels2dtI+U7b4X88Nlu++sM2sM+E/7MgMdavta95+tl2v5qzgoBeA7sfGB8ubO44x61TUmyKOfdGOftE60e6jFQ1k3Dr7Xu0+qcKUdSS4qpS02v9N264Ku5ac81d2JBP+vNEUtc8IR8WA5Gd4zKfs7+zFOeHf/4/85xfMOm2fCXszZZbZrhz5kKY1sX6VmKouwbamfriw3XEcZ0RwJCkxx3GOPo6okZjjOEcZigdFdBwnvfQkCkkLrsQcx+lLipSYp2xzHKcPoprT0e95RO4TkZ0isiar7Fsi8raIrBKRn4tIbVw+Q0QOisjK+Ph+Ln0dMSOxPSfapusZ/xk28xe129mcxr1lt9VZbbe1a0mdKavaE24vyY1izDt2+vok6t6yv5ombHeJt/7h7mD58U98yazz3at+ZMq+ddkVpmzfDeFIFQBjHglHxihtsV0AihK8A1rm2a4exZtt14y5x4Zzr27pnGXW6TTcVABGtdqLRfun2N9Z5d0TTdmCjZ8JlrfNsKOEdC4ZFSzXlgL8SxfWB+x+4LtEkXB6WAzcoqpdIvJN4Bbgq7Fsg6rOzacBH4k5jtOHQiXPVdWXgN2HlD2nqj0/Wa8QZfoeMK7EHMfpQx7bjup79kbHhx1cL8zVwC+z3h8nIq+LyG9E5PxcTjBippOO44wgcp9ONqpqOHFmP4jI14jiFT4YFzUA01S1SUTmAU+KyJxDA04ciisxx3F6U8AM4BYi8nngYuDDcRQcVLUdaI9fLxeRDcCJwLKkc7kScxynL4dRiYnIAuCvgA+q6oGs8vHAblXtFpHjgVnAu/2dz5WY4zi9KKSzq4g8BFxAtHa2Bfg6kTWyHFgcBcnhFVW9FvgA8Pci0km0Z+BaVd0dPHEWQ6rEMqXQZlivS+ycCLSNC5uat3zCdrGY8IJ9aZ2VA4sOXtwRNq+3jbXb2neinYSjuMN+UtrG2jaX2vW2mf+clZcFy5MSdPzNy9eYsuIZdh+Pvb3FlHXUhSM6NJ4Sdg0AyNgeBRRts90oLvvYf5uyZ+8JRx7J1NvPQOsU+5rr3zBFHBxvn3P3PPs7O/HHxcHyxGfgd+Fnv9gOEpIXkimMFlPVKwPF9xqffZwobmFe+EjMcZzejKBYYbngSsxxnD6MlKitueBKzHGcvvhIzHGcNONRLBzHSS8K5LC5e6QwpEqs5ADUvxG2qjSeFrbQAFQ2hDf/jlldbdZpOtX+Eib+1p7wd1XZ/bA2Lzcbcdz7Y/QLdvan8ma7/51VthWs+N/DAfMnLN1h1tl7hr2hfN8M+350l9mWV6uPLdPt60rKmXBgki2zch8AHJwU7scxK+14+NXbbKtg1UbbIgv2/ahdb9/H/dMrg+VJeSKqN4XPl7FjG+SFr4k5jpNaPCii4zjpRtWnk47jpBsfiTmOk25ciTmOk2Z8JOY4TnpRoDs9WqxfJSYiU4niY08gurxFqnqXiNwGXAPsij96q6o+k3wy6C4Lm7yn/6LZrLbjD2uC5aN22XbgpM24B8bbJvQ3vhKOUQ9wwTXhjdKzvmqnvd/+0YFF3k3aHF7RZLsHbD87HO+/5t/siCaZs20Xi6pt9j0ua7E34LfXhB+tulX2dVVut6+r7nV7Z3Pb5LCLAkBxR9gVYceZti9CTcIG+11njjFlo7cmJAlIoLss/DzWrbLrVOwJt5X03OTDkTYS6wJuUtUVIjIaWC4ii2PZd1T1nw9f9xzHGRaOJOukqjYQhY1FVfeLyFpgyuHumOM4w0eaRmJ5JQoRkRnA6cDSuOj6OHfcfSIyttCdcxxnGNA8jhFAzkpMRKqJApbdGAfuvweYCcwlGql926i3sCcTSmd76+B77DjOYUUA6dacjpFATkpMREqJFNiDqvoEgKruUNVuVc0APwTOCtVV1UWqOl9V55eW23sFHccZORQqA/hQ0K8SkygI9r3AWlW9M6s8e0vuJ4E1h9Z1HCeFpGw6mYt18lzgKmC1iKyMy24FrhSRuUSXshH4Un8nkoyaZvmOOjuGuny0KSx41F6G04Qr2z3PNoXP/t6X7Ypzw8Vjxtp2jiSze0d1/pE7ADLldr2JS8MRNeT0OWadupe3mrJ9Z0w2ZUnUrwovHZQ07DHraLOdXrD7fdNNWdMc211i/Irw/Ri31O7He5+caMoqGu3/3PImO5pJcZMd/aJtevg53jsr7C4DULM+nJRCugqhWY6wvZOquoRomnwoyT5hjuOklgJmO7qPKL/kTlU9JS4bBzwCzCAaAF2uqnviWd9dwEXAAeDzqrqivzbysk46jnOU0BPJor+jf+4HFhxSdjPwvKrOAp6P3wN8jCjX5CxgIZHxsF9ciTmO0xstnHVSVV8CDs0deQnwQPz6AeDSrPKfaMQrQO0ha+9BXIk5jtOXw7uwPyF2ogfYTrSlESIn+veyPreFHBzrfQO44zh9yMN9ol5ElmW9X6Sqi3KtrKoqMrgVOFdijuP0JXcl1qiq8/M8+w4RmaSqDfF0cWdcvhWYmvW5Y+OyRIZUiUm3msk2Wifa5uQyw5Vi3NKdwXKA7rIJpqzu9bB5GqB7TJkpK90WjrSx/5TxZp2qNdtN2bYv2REuGj9piih+e5Qpm/arsCk/KfpCzbu2e0tSYozWGXaiFssNJDPadnjOTLJdZtrr7Odj2iPvmTLLRaQtITnKxKW2e0t3qb0Cc2CS/b1UtdsRP9rqwi4iLe/rMOv8blq4rfZvFmCFSIHDmyjkKeBzwDfiv7/IKr9eRB4Gzgaas6adJj4ScxynF0LhvPFF5CHgAqJp5xbg60TK61ER+SKwCbg8/vgzRO4V64lcLL6QSxuuxBzH6UumMEMxVb3SEH048FkFrsu3DVdijuP05vBPJwuKKzHHcfowUjZ354IrMcdx+uJKzHGc9HKEbQAvJF0VRew+KWzOr3nXNifvnxZ2ezhwYp3d1ig7UcjOv7MjS+iz4aQkAPsuDZv5T7rbdvXoSnAbOOGuDXY/xo8zZQen2wk19k8PJ82Y8Au7rR2XzDRlRca9B6h/1Y4EYdFVZ7shdFbbj+OeE+1IFa0TbFeV3ReGE4yUr7X7UfNmgovFxIHFxJNO28WidkX4+Wmvsd2EKo0kOY12IJDcOdKyHTmOc/Tha2KO46QbV2KO46QWBTKuxBzHSS2+sO84TtpxJeY4TmpRoDs9LvtDqsQy5bB/RljWVWVHKajaFr6hrRPs7r/+tbtN2QXXXGPK2sbaX97ffOyJYPm/rv0zs87B8barx7Rn7OgXRS22mT+Jzqpwe23vnxosh+RkJkmuDZlq2/2i+O1N4X7Mep9ZZ9+MhAQov7WjaSRR/2rYdWfDn9uROxrPst1iklx3jllu91HabReihj8OJyYpb7afxYrd4fMVdRVC+SioKzHHcdKMTycdx0ktbp10HCf1+EjMcZxU40rMcZzUogrd9l7PkUa/SkxEKoCXgPL484+p6tdF5DjgYaAOWA5cpaq2CQYoruxi3NxdQdnOd+3N3BC2WhUltHb6HV82ZWPb7XTznZX2LfnODy4LllclWJGSNrY3nW7Hva9/1X6Iitrt9lonha1n496yN423l9uW4VHGRmOw4+gDdM4/IVjeXZZg3VthW2RLGuzN5lpmW1C768J5ALrG2NdV/2o4lwKAltrXvO/E0aYsc0J4Yz5A/ZrwdRclxOVvNs7XvapAWRhTNBLL5YrbgQtV9TRgLrBARM4Bvgl8R1VPAPYAXzxsvXQcZ2gpXAbww06/SizOxtvjAFMaHwpcCDwWl2dn8XUcJ9VoZJ3M5RgB5LQmJiLFRFPGE4DvARuAvara4yWZU6Zex3FSgIIWyNlVRE4CHskqOh74W6AWuAboWV+6VVWfGUgbOSkxVe0G5opILfBzwHa7PgQRWQgsBCg7xl4DchxnBFGgbUequo5oGapnMLSVSId8gWg56p8H20Zeq4Cquhd4EfgDoFZEepSgmalXVRep6nxVnV9SYy9uOo4zQlCNUrblcuTHh4ENqhrekzZA+lViIjI+HoEhIqOAjwBriZRZj7kuO4uv4zhp5/As7F8BPJT1/noRWSUi94mIvWG1H3KZTk4CHoiHgkXAo6r6tIi8BTwsIv8AvA7c29+JVKGjK2yiHptgGh63LmyCPvab6806S146xZQVddoblzOltguAtRHd2nQNUNFkihj/8g5T1jnZjvVf3Gn/AlY1hB+szQvCrgYAk5bYLieV2+229s60N1Hvfn+4H0Wd9r3aMdkUUVxiP+PFb9vXNnZtuP91K+x+JLlRdI9J2PTeYf9Tj95ku4/snBfu/6Tntpt16gyXk5KWRC+nnNHcR1n1IrIs6/0iVV106IdEpAz4E+CWuOge4HYiI+HtwLeBqwfS136VmKquAk4PlL8LnDWQRh3HGcnkNcpqVNX5OXzuY8AKVd0B0PMXQER+CDyddzdjCuQZ5zjOEUPPBvDCulhcSdZUUkQmZck+CawZaHd925HjOL1QQAu47UhEqojW0r+UVfxPIjI3bm7jIbK8cCXmOE5vtLBBEVW1lWh7YnbZVYU6vysxx3H6oCPEGz8XXIk5jtOXFIWnFh3CTZwisgvocXSrBxqHrHEb70dvvB+9SVs/pquqnbwhB0TkV3F7udCoqgsG095gGVIl1qthkWU5mma9H94P78cI6cdIxF0sHMdJNa7EHMdJNcOpxPpsTRgmvB+98X70xvsxwhm2NTHHcZxC4NNJx3FSzbAoMRFZICLrRGS9iNw8HH2I+7FRRFaLyMpDduIf7nbvE5GdIrImq2yciCwWkd/FfwccmmSQ/bhNRLbG92SliFw0BP2YKiIvishbIvKmiNwQlw/pPUnox5DeExGpEJFXReSNuB9/F5cfJyJL4/+bR+LIEI6qDulBlLpoA1GY2jLgDWD2UPcj7stGoH4Y2v0AcAawJqvsn4Cb49c3A98cpn7cBvzlEN+PScAZ8evRwDvA7KG+Jwn9GNJ7AghQHb8uBZYC5wCPAlfE5d8H/mIov6eRegzHSOwsYL2qvqtRireHgUuGoR/Dhqq+BOw+pPgSooQrMESJV4x+DDmq2qCqK+LX+4mCbk5hiO9JQj+GFI3w5Dw5MhxKbArwXtb74UwyosBzIrI8zgUwnExQ1Yb49XZgwjD2pSARNweCiMwgil+3lGG8J4f0A4b4nohIsYisBHYCi/HkPCZH+8L+eap6BlHAtutE5APD3SGIfomJFOxwcA8wkyi5QwNRxM0hQUSqgceBG1V1X7ZsKO9JoB9Dfk9UtVtV5xLlrziLPJLzHG0MhxLbCkzNem8mGTncqOrW+O9OogwswxmpdkdPoLj4787h6ISq7oj/gTLADxmieyIipUSK40FVfSIuHvJ7EurHcN2TuO295Jmc52hjOJTYa8Cs2NJSRpQ84Kmh7oSIVInI6J7XwB8ziOiSBeApooQrMIyJVwoZcTOPNoUoR8NaVb0zSzSk98Tqx1DfE0/OkyfDYU0ALiKy/GwAvjZMfTieyDL6BvDmUPaDKExvA9BJtLbxRaKgcc8DvwN+DYwbpn78FFgNrCJSIpOGoB/nEU0VVwEr4+Oiob4nCf0Y0nsCnEqUfGcVkcL826xn9lVgPfAzoHyontmRfLjHvuM4qeZoX9h3HCfluBJzHCfVuBJzHCfVuBJzHCfVuBJzHCfVuBJzHCfVuBJzHCfVuBJzHCfV/H97uViw8VqZ/AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imshow(a)\n",
    "colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5dcf3d3f",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fbafe7a80a0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAD5CAYAAABPqQIFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsBUlEQVR4nO2deZhdVZXof6vGpKqSGjIUGUlCAggoCYQIoogiitgK2D4EbQS1BfqBQjdqI3QLDc9+IoLSreIXGxrx0QwySBpRCCiTMiUQAklAEhJIZahKJZVKJTXfWu+Pc/LerdRep26lbqrqhPX7vvPVvXvdffY++567ap+91l5LVBXHcZy0UjDcHXAcxxkMrsQcx0k1rsQcx0k1rsQcx0k1rsQcx0k1rsQcx0k1RYOpLCKnADcBhcB/qOr3kz5fXFqupeU1QVmm1K7XMyrsBnLImHqzTl1XlSnrrB9lyhL7URwuL9lhu6n0FIotSxj9wu6EfhTaMsmEy5Ouq7Bz4OcD6C63ZUVthqAnoa2eBHcfexjRggShIeoebVc5dJx9X23oGmPKdjaVmbLSyg5T1t4VvrEKm+05hnUvdjVvo7t1V8KA9M8nPlKuW7clfPFZLF3e8YiqnjKY9gbLXisxESkEfgqcDNQBL4rIIlVdadUpLa/hvSdfEpRtn23/MlvnhH9liz72Y7PO5es/bcre/vHBpqz5ILsfbQeEf4HTFttfeHtVwvnG2zdpeb39a2+vtu/R0u1hRdA8y26rYoOtPEqb7X7UH2Ofc/zy8DkLO+y2itrstnqK7WvOlCT8ozDqNR5p13n0b35oyq7c9DFT9swD80zZnFPWmLJVm2uD5VW/tf9L7JwS7v+6W2806+RK47YMzz8yNafPFk9aM37QDQ6SwczEFgCrVfUtABG5CzgNMJWY4zhpQMlowpR5hDGYNbEpwPqs93VxmeM4KUaBHjSnoz9EZJqI/FFEVorIChG5JC6/WkQ2iMiy+Dg1q853RGS1iLwhIp/or41BrYnlgoicD5wPUFJWta+bcxwnD/QkLV4OjG7gMlV9SUTGAEtFZHEs+5Gq9np2F5HDgLOAw4HJwGMicrCqmms2g1FiG4BpWe+nxmW9UNWFwEKAipppvlHTcUY4itKVp8dJVd0EbIpft4jIKpKf2E4D7lLVDmCtiKwmWrp61qowmMfJF4E5IjJTREqItOeiQZzPcZwRgAIZNKcDGC8iS7KO863zisgMYB7wfFx0sYgsF5FbRaQ6LhvwMtVez8RUtVtELgYeIXKxuFVVVyTV6S6DLfPCenPSs7ZPwfopYQvfSXd+y25rrG0xLJpv6+6aFfZkMVMSrlfc3G7W6axI8IdIoHWi3cc/f/vHpuzkv/96sFwTvumSnfZ/3a5y24o39i37nC1Tw/1vm2f5XsABD5SYsiSr5tYj7DGetnhXsLxlqm35O/WfvmnKLNcGgJ7JtuzV9bZw1Iqwv0dPkX3Nso+faXJZ74ppVNX5/X1IRCqA+4BLVXWHiNwMXEukM68FbgC+sjd9HdSamKo+DDw8mHM4jjOyUCCTxxBdIlJMpMDuUNX7AVS1Pkv+C+Ch+G1Oy1TZuMe+4zh96Mnx6A8REeAWYJWq3phVPinrY2cAr8WvFwFniUipiMwE5gAvJLWxz62TjuOkC/3/61354HjgHOBVEVkWl10BnC0ic4kmfuuACwBUdYWI3EPkb9oNXJRkmQRXYo7j7IEqdOVJh6nqM4Q3f5nLUKr6PeB7ubbhSsxxnD0QMkmbVUcYrsQcx+mFAkl78UcaQ6rEinZB7Yvh5cCth9tdmfZI2P2ifHWjWadp7jhTtn2O/V+moMterpxxf1OwfNfssWadijrb/QLsaBrdo+0+HvPTS03Z+LbwWFW/bl9XQbd9x3aX2rafpOgXu94bvu6it+1rbp6ZsLG9KaGPo23Z6rPD7VWtTIg8kvCr0NO2mrLXjrrLlM3/14tNWen28HeTtOm95vXw9/xOe360j8/EHMdJLZGzqysxx3FSigJdmh7vK1dijuP0QhEyKXIhdSXmOE4fetQfJx3HSSm+JpaAFkJHZXiaOmFZl1mvdWK4m2WF9sbfpJDFVW/alrqOKrveGxdUBstr/2zX2XqCHXe9cq3dj6ZDTRGT/2Rvlt8+JzxWxQl5ALbPtm+DcSttZ+ntc+zx157wmBQlhH9Psgombbz+2089ZsoaOsMx8Z947f1mne4yu4/yezsa88Gr/6cpq91qf9fFxgb85pn2gDQdFh77rhfzoXyEjK+JOY6TVqLIrq7EHMdJKapCp+5dCKnhwJWY4zh96PE1Mcdx0kq0sO+Pk47jpBZf2HccJ8X4wn5/GJb+be+xbehj1ofN/K0z7Y3XVW+2mrKf3PkzU3ZwsR17/Zh/+rtg+eYP224IUx6zXRu2zLVvlFGN9prE5gX2omtFXbi9io22W8aOg+zbYPOxdh/HrLOvbeyy0mB5UZtdp7g1IY7+e+3x+D+/OtmUZYzbKjM1oe9rbVnZFvu7Hvu2KWL0Zvt+rPtY+D4e3WD3o/Iv4fJNSfEGBkDGnV0dx0kritCVlFlmhJGenjqOMyT4wr7jOKlGEX+cdBwn3aRpYT89PXUcZ0hQhYwW5HT0h4hME5E/ishKEVkhIpfE5deLyOtxBvAHRKQqLp8hIm0isiw+ft5fGz4TcxynF9HCft62HXUDl6nqSyIyBlgqIouBxcB3VLVbRK4DvgP8Y1xnjarOzbWBQSkxEVkHtAAZoLu/dOaF7T1UvxFOK7/hxAqzXuuEsMYft9K2Jxd02qbw9d22a8ZBRXY0jZKWcLSB6Q8FiwEoMmLeA5RsD7shAExdHI7nD9BdZcep76oIf6Vb5touLJnRdoSF5z93gyk75uFLTVnVK+H2Osfaay21Z9s+Cl13HWjKMvYwsmt6+D546jP2dX1o8aWm7MB77f53Vdgzk5aZtutOYZspMrGireRL9+RrYV9VNwGb4tctIrIKmKKqj2Z97Dngc3vbRj5mYh9RVTtjh+M4qUKRgQRFHC8iS7LeL1TVhaEPisgMYB7w/B6irwB3Z72fKSIvAzuAf1LVp5M64I+TjuP0YQAzscb+nsAARKQCuA+4VFV3ZJVfSfTIeUdctAmYrqpbReRo4Dcicnh2nT0Z7JxRgUdFZKmInD/IczmOMwKI8k4W5HTkgogUEymwO1T1/qzy84C/Ar6oqgqgqh2qujV+vRRYAxycdP7BzsQ+qKobRGQisFhEXlfVp/a4gPOB8wFGlYQjozqOM5LIXwZwERHgFmCVqt6YVX4K8G3gw6ramlU+AdimqhkRmQXMAd5KamNQSkxVN8R/G0TkAWAB8NQen1kILAQYWzElRXmFHefdSZSyLW/WyeOBc4BXRWRZXHYF8G9AKdHkB+A5Vb0QOAG4RkS6gB7gQlXdltTAXisxESkHCmKLQznwceCavT2f4zgjA1XJ+VGx/3PpMxCc1j1sfP4+okfPnBnMTKwWeCDWokXAf6nq75Mq9BQLbQeE3QPGvWa7Irzz6fAErupNe6ClxP5Pcv3nzjJlP8zYrhntx4XbK2uwXRQ0IarAAS/atvWt86pMWfkm2w2k7gvhcaz6o/1VT3zJHvtjM5eZMmrseqdf8ESw/Nnz5pl1fvr1e0zZZwu/bcqsKCcARa3h7+xTa+zz1eyyHxjaa2xZxcZOU7Z9lu0HsuPIcL1HTrrJrPP5674VLBf7VhwQ74p4Yqr6FnBkHvviOM4IIIon5nsnHcdJLR7Z1XGcFBO5WPhMzHGclJLnvZP7HFdijuP0IU2heFyJOY7TiygUjz9OBumqEDZ9IDxNrXnNrldkBHQo3hmOiAGQGW1f2vqT7YgCSe4GE5aEt2+t+bwdFaOkyTatV6617eGlzbZsx4F2RIqlH/73YPnc1m+YdUZvtR8dxh5sR9MoeLDGlP1+8QnBcknYQHLeRf9gyorH264NHX9r+0JuWxfu4+Q/2P34/U22a8NHvvv3dsUEMqNspXDVBxYFy0958utmnZoOYzzy5E7ua2KO46SWKIqFP046jpNSom1HrsQcx0ktPhNzHCfluMe+4zipxa2TCRR2QMW68OA0z7bNKiXN4Tqd1zabdToXHmDKihJimrdOsC11JTtKguXVK+3z9RTa15UpsW+UznJbVv1mhylb3RXu/4un/Nis89nf2ha3krurTVnDBxI2vheEr/uQW1qD5QA7D7Stxm0T7PEofWC8KZvYHu5H45H2+T56ZYIFMsEH9K3P2lbjz3zgBVP27zf+dbC8JmEz9/bDwteVSQzBkDv+OOk4TmoZYIz9YceVmOM4vVCg22dijuOkGX+cdBwnvag/TjqOk2I8KKLjOKnHZ2IGWgDdZWFZ+bytdr2HxwXLGx+dYjc20xZNfqLFlLXMss383aPC9vXKtXYg/ZZp9gbwXZPtdYeK9bZ9veES20fkkm+FNw1vO8veLF84274NJiyz3TmKx9ux/pceH0wCzacftDc1bzrDjlE/asVoU1b9pj3+az8THv/a5+zxrX+//QMet9yWVU63XX6evekYU9Z0UniMJy4Ou/QAzLovfA80Ng0+yH4+gyKKyDTgdqKcHEqUIfwmEakhyvo9A1gHnKmqTXGKt5uAU4FW4DxVfSmpjfSs3jmOMyQoQndPQU5HDnQDl6nqYcCxwEUichhwOfC4qs4BHo/fA3ySKNfkHKJ8tTf314ArMcdx+tCD5HT0h6pu2j2TUtUWYBUwBTgN+GX8sV8Cp8evTwNu14jngCoRmZTUhq+JOY7TG903a2IiMgOYBzwP1Krqpli0mehxEyIFtz6rWl1ctgkDV2KO4/RigGti40VkSdb7haraZ0FURCqIkuJeqqo74ny1UXuqKiJ7Hc7RlZjjOH0YgBJrVNX5SR8QkWIiBXaHqt4fF9eLyCRV3RQ/LjbE5RuAaVnVp8ZlJr4m5jhOLxQh01OQ09EfsbXxFmCVqt6YJVoEnBu/Phd4MKv8SxJxLNCc9dgZpN+ZmIjcCvwV0KCqR8RlQfNof+cq7ITKteGU8+1NYTcKgOeu+Umw/NBfX2TWuf0zPzNl1zxxrinbfKpt5h/3RNhcP265beLXglGmbOd7bPeFjirbvD670jblbxlbFSyvurfCrIPaZvmdU+1+TLjHvn3O/OcvhpuaY3djzAu2G0VBl/20seV9dr2ZD4SjZhS227kUxv4lfI8CdFfZ32fDE3bOgeJd9jmn3xN23SndkuAWs21nsFy67XYGQh6dXY8HzgFeFZFlcdkVwPeBe0Tkq8DbwJmx7GEi94rVRC4WX+6vgVweJ28DfkLk67Gb3ebR74vI5fH7f8zhXI7jjHA0jwv7qvoMmBrxpMDnFbBnJwH6nQ+q6lPAnulkLPOo4zj7AaqS0zES2NuFfcs86jhO6nmXbQDvzzwqIucTed5SUlY12OYcxxkCRsosKxf21jpZv9uLdg/zaB9UdaGqzlfV+cWlCYvLjuOMCFQh0yM5HSOBvVVilnnUcZz9gHxtOxoKcnGxuBM4kcgztw64Cts8mogWQkdlWG+WbbVNw/Oe/1KwvGqlPYhXPnKBKeuebuvuQ79ve4p0TgyXN84dY9bpGmP3cdRaO8LFqac9Z8r+cMuxpqzIyGO/9Qi7HxOW2S4W7TX2WG053q43uj48Jm3j7FtuzHr7Htg+x87QkbHzc9BVGRau/rLdj3F/tk9YkODB0JPwayreleDGMjlcsb3KjqhS2hx2K+lpGLz/upKux8l+r1hVzzZEfcyjjuPsD7zLFvYdx9n/0L3eyTj0uBJzHKcP+9XjpOM47y4i62R6tlW7EnMcpw/+OOk4Tqrxx0kDyUDpjrCpubPCnr6WLRobLJ/w5EazTvNRB5iyXbVJ5no7EoGVT7Q5ITJDzQr7X1rbBLven25cYMq6ppkixmwM+wCU19v9GF0XjogAUNhhu4+Mfdv+zjrGh6NfbD/YrMLoevt7GfeqnZRky1zbJaJ0azhSyISnbfeFbUfYYzX1cdvHom283Y/iFrv/jYbHTMUEO4pF928rg+WZ4sErH2Xk7IvMBZ+JOY7ThxQ9TboScxxnDxR0hGwpygVXYo7j9MEfJx3HSTVunXQcJ7Xsd3snHcd5l6GAK7EEjGlq5V9sc/L6k8Nm/uq1b5t1Co+wg82O3mpHFCirtxOFNM8MR52YuMQ+36htdkKKcS/uMGW7Dqo2ZUXt9ly/8YjwVzp2XUIUhUlVpqxig93/ou6EZw7D+6J6lV2nYqOdOKV1op2w5MDfbDFlmz4a9mMZ3WiPR2Gb/QPeNcn+yTx+yfWm7FPbvmnKDvn6S8HyrV88yqwjhqfH3mdv7I0/TjqOk2Ikb9ZJI1va3cAh8UeqgO2qOjfOEL4KeCOWPaeqF/bXhisxx3H6kr+Z2G3skS1NVT+/+7WI3ABk5yBco6pzB9KAKzHHcXqj+VvYV9Wn4hlWH+LEumcCHx1MG+nZqu44ztChOR6D40NAvaq+mVU2U0ReFpEnReRDuZzEZ2KO4wTIeSY2XkSWZL1fqKoLc6x7NnBn1vtNwHRV3SoiRwO/EZHDVdW2gDHESqygs4eyTe1BWU+Jvfl38p/CdeTow806FS/VmbLW904xZV3lCbHXXwmPZWGjPca607a66hQjaD/QPMvuR/UbtgV15n+uD5Zv+OwMs05Rq/0vdVRDqykraLZlmfFhi/KW99kZryrX2JbQ0mb7/ugaZ2/m7jqpOVhe/GA4qABAQZf9A25cYPfx5B9+y5Rl7K+abWeFrZAVG+22uoyACWIbXQdG7udpVNX5Az29iBQBnwWO3l2mqh1AR/x6qYisAQ4GlgRPEuOPk47j9Ga3n1gux97zMeB1Vf1/sw0RmSAihfHrWcAc4K3+TuRKzHGcPqjmdvRHnC3tWeAQEamLM6QBnEXvR0mAE4DlIrIMuBe4UFW39deGr4k5jtOXvDnNhrOlqep5gbL7gPsG2oYrMcdx+uLbjhzHSTP52r40FLgScxynNyqwPwVFNPY+XQ18Ddi98/YKVX24v3N11sK6b4RVfOHrtpl8+u/DMeC3v8eO/z7qgDJTNvrtFlPWOtM2vWdGh2OoS6XddxJkbdPs/k+5397cvnOe7SKy9cTpwfKJS21XjyT3ltZpdv/LEszw7eNHhc93oO028OY5o03ZzPvtep1V9ubwsgfD55zwhO2CU/Ahe3xnPGjfO13Vdv9L395qynYdFva/KFv2jlln3ZcPCpZnnjarDIwUzcRysU7eBpwSKP+Rqs6Nj34VmOM4KWJoPPbzQr9KTFWfAvo1czqOsx+xPymxBC4WkeUicquI2MGvHMdJF0Pj7Jo39laJ3QwcBMwl2u90g/VBETlfRJaIyJLMDntdxnGckYNobsdIYK+UmKrWq2pGVXuAXwBmpldVXaiq81V1fuHYhAVwx3FGDvv746SITMp6ewbwWn664zjOSCBNM7FcXCzuBE4kCrlRB1wFnCgic4l08TrgglwaK63LMOfy7UFZ+0F2Cvid08Km63GPrTXrbP3YTPt8k+wlvIl/tk3hdIcDm3dPsF0lusbY17V9li1rHR92lQDY8kHb3WBUXfj/UvUrdvz6nvG2O4rYTSFdRqB3oOLVTcHyCeOnmnXKttiN7TjQdqMY844d1aOwMzweO46abNbpLrXXelqus8ex/d4qU7btcweYskN+Ho600T3dDn1x4K2rg+UbG+3+DYgRst6VC/0qMWPv0y37oC+O44wERtCjYi64x77jOH1xJeY4TprJW3DFIcCVmOM4ffGZmOM4aWUkWR5zwZWY4zh92Z+sk3lFQIvCERN6iuxB6yoLm8lbj5xm1qmos83uW44MR1gA6K623Q2KVoRdOrqnV5l1th1iu1FMfqLJlGmxHVliwpN2YpK3zgu7MDQusN1KehLugpqVbaasYJctqz853I+CLrutkib7O9v8fruTHVWlpqxiw8AXd35z1fWm7Mxvf9OUjWmyL+5/Xf4rU3blq18JlndW2r+JSX8K3x+6077fBoTPxBzHSTP+OOk4TnpRt046jpN2UjQT85RtjuP0JU8bwONQXQ0i8lpW2dUiskFElsXHqVmy74jIahF5Q0Q+kUtXXYk5jtOHPG4Av40cI0OLyGFE+SgPj+v8bHcy3SRciTmOs88YYGTo04C7VLVDVdcCq0kI87WboV0TKyhAR4XN4SXNtnm94ahwBIOq1XbUg+7RtgIvr9+7Vcv2+bOD5Z2V9jBOWN5uygqa7SCRWmyfs2es7QbSVRG+tvFLw5ESAKTVjnzQcsQEU7Zr6iRTVrMq7H5R0Gq7IeyaUWHKyubav4OmDZWmrKArPI5JC9enX/MtU5bkRtFaa7s3XHFD2I0CoGxruDPVb9ht7ZgZdhPKrMzTvCT3NbHxIrIk6/1CVV2YQ72LReRLwBLgMlVtAqYAz2V9pi4uS8QX9h3H6c3ArJONqjp/gC3cDFwbtcS1RJGhbS3fD/446ThOX/ZhZNeEyNAbgGwP9qlxWSKuxBzH6YWwbyO7JkSGXgScJSKlIjITmAO80N/5/HHScZy+5MlPbCCRoVV1hYjcA6wEuoGLVNUOHxzjSsxxnN7kMYrFQCNDq+r3gO8NpI0hVWJaVECXEc9d1B61sevCq4zrT7atdJOesS1uZfW2JXT7bDsVfdN7wuVJX3i3Hcqd0jLbqpZ507bUTXzJXnWteS28aViL7JWDnnF2W0mMfaPFlNUfF762A56xrYxj/mznTGj8Qq0pq1xl38aTngy3N+ZnW8w6S14OW6EBCrptC2RHtb1hu6LOnlB0VoS/m5Lt9vlqnm8IlhftSthhPxB825HjOGnGN4A7jpNuXIk5jpNaPNuR4zhpxx8nHcdJN67EHMdJM/tVUEQRmQbcDtQS6eeFqnqTiNQAdwMziBzWzow3cZrMmbmF393xi6DsfQu/btZrn264RBTY/y46VtiX1jDfdjcYdbC9Ubryt2G3gSTTemGbHc+/dLst66gxRTTNsTe3FxqeJduOGGufMIHGj9quKqVrbReRjurwr2DnDDvWf/EO+6I76+xfVUlCsJYtRm6Bhl/Z/ahK2MeSlI/gv//+B6bsCyu/ZMp2Pn1AsLx7dDjwAcCOL4fHvv3KfiPX9E/K1sRy2XbUTbTL/DDgWOCiOO7P5cDjqjoHeDx+7zhOypEBHCOBfpWYqm5S1Zfi1y3AKqLwGKcBv4w/9kvg9H3UR8dxhpp9uAE83wxoTUxEZgDzgOeBWlXdFIs2Ez1uOo6zH5Am62TOUSxEpAK4D7hUVXslPlRVUy+LyPkiskRElmzZ2u9eTsdxRgIpmonlpMREpJhIgd2hqvfHxfW7Q2rEf4ObuVR1oarOV9X5E8blYdHRcZx9SxwUMZdjJNCvEhMRIdp1vkpVb8wSLQLOjV+fCzyY/+45jjMspGgmlsua2PHAOcCrIrIsLrsC+D5wj4h8FXgbOLO/E63cPJH5Pwi7UrQfbu++H702bGruODQcxx2go9LWzzUrTBEVj9pRLN4xEkgV2aHyKa+zv2lNmJge8JydP+DtT9vXVr42fNLmOXZbSf1YMHudKXtlzaGmbKIRyu4H195s1vnKc+eZstpFtjtKUbs9Vgd8c02wfNWDh5h1JGHVo+UYOwLK5/7Zjs3fPs625ZUaaRh2HGTfO1V/CI/H5pb8xDlN05pYv0pMVZ/BtqaelN/uOI4zItiflJjjOO8+9quZmOM47zIUD4roOE562Z0oJC24EnMcpy8pUmKess1xnD6Iak5Hv+cRuVVEGkTktayy60XkdRFZLiIPiEhVXD5DRNpEZFl8/DyXvg7tTEygx9iYXzvVDoCxffOEYHlBXZLZ3R7gkp32A//mBaWmbNzycL36D9km/qJWO7FE6xS7H09/96em7OiFl5qynQeHXVWmPWyb+Dd+0P5f9upDthtFSaspontUuL1/+Ne/M+uM+bR9DzQdYieFQe3+1y0LJ/0Yt9W+P7pH22NVtMX+PrvK7XqdVXZ7pU3heoVt9vl2TQ3Leuzu5U5+fcBuA35CFAlnN4uB76hqt4hcB3wH+MdYtkZV5w6kAZ+JOY7Th3wlz1XVp4Bte5Q9qqq7//M/R5Tpe69xJeY4Th8GsO1o/O690fFx/gCb+grwu6z3M0XkZRF5UkQ+lMsJfGHfcZy+5P442aiq8/emCRG5kihe4R1x0SZguqpuFZGjgd+IyOF7BpzYE1dijuP0Jo8ZwC1E5Dzgr4CT4ig4qGoH0BG/Xioia4CDgSVJ53Il5jhOX/ahEhORU4BvAx9W1das8gnANlXNiMgsYA7wVn/ncyXmOE4v8unsKiJ3AicSrZ3VAVcRWSNLgcVRkByeU9ULgROAa0Ski2jPwIWqui144iyGVIkVVXYx8RN1Qdm6Vyeb9QoPDoeJOO7AdWadZ0YfZsqqV9phG0766xdN2dP/eUyw/A+n/NCs89DOw03Zzf/1KVP2/n+71JSN2mmKKGkO29i3zLXrjGq0Za2TbDeQsWtsu1DrpLALQFKEiDG/rjJlhaX2r+qki581ZQ/d84FgeWmz3ZHCzoTo8e/Y15yxPX7I2J471KwMR8a4+z9uMuscd99lwXLN0y9aevKjxVT17EDxLcZn7yOKWzggfCbmOE5vRlCssFxwJeY4Th9GStTWXHAl5jhOX3wm5jhOmvEoFo7jpBcFctjcPVIYUiXW1VLMxqfD26Qqttv1OpsqguWvPHuEWWdynW19qvnGOlP230vnmrKK8nB5a0KQ+pvvtC2Q3RUJ8ffbbQvZQZ8Kx40HOLIqbP19qiG8ERqg7mXbMly8I2FT8xhTRPUb4fG/4bqfmXW+9rNw/gWAsgZ7rB743XGmbOqLHcHy5llGJAJgzAZ7Q3/JOfWmbHubbZ4seb7arjc73JdTr/6mWafKuOXq7bQTA8LXxBzHSS0eFNFxnHSj6o+TjuOkG5+JOY6TblyJOY6TZnwm5jhOelEgkx4t1q8SE5FpRPGxa4kub6Gq3iQiVwNfA7bEH71CVR9OOld19U7OOP2ZoOzXi48361lm/vbx9kC3HB2ONQ/QsLHWlE1/yHYpKGwLm+u/eGN4My5A53TbVl220d5MnHRtmxfONGUN3WFZa63d1ritCfb0L9q7w5uWhnMfABS1hft/5VtnmHXKN9v9GN1ou8yUtNjX1jDf2Hm9oNmsM6rM9lPouc2+d3r+R4spS/DCYcdx4fZKV40267RNC7uBZH6XH+Wzv83EuoHLVPUlERkDLBWRxbHsR6pqh3BwHCed7E/WSVXdRBQ2FlVtEZFVwJR93THHcYaPNM3EBpQoRERmAPOA5+Oii+PccbeKiO2S7DhOetABHCOAnJWYiFQQBSy7NA7cfzNwEDCXaKZ2g1Hv/N2ZUFqbwmtKjuOMHASQjOZ0jARyUmIiUkykwO5Q1fsBVLVeVTOq2gP8AlgQqquqC1V1vqrOL6tOCG/pOM6IIV8ZwIeCfpWYREGwbwFWqeqNWeWTsj52BvDannUdx0khKXuczMU6eTxwDvCqiCyLy64AzhaRuUSXsg64oL8TNTVV8MADHwzKalfaJvQd08P2acnY7hCV42xzd2bReFPWeERCCvv28HAVttvfZkFCvPbJf7TN/Buvst0Ndm23lx/basP1SprMKnSX2f/Lih+aaMpKwsFFANgx3RirO22bUNs0e6zqP2zfH+VrbP+FotZw+c72cC4CgLF32+E5ekrtPspzlaYsM9a+R8Y+G3alKOyw65QuCKdi3FKSj/AT+9neSVV9hugxeU8SfcIcx0kvecx2dCtRfskGVT0iLqsB7gZmEE2AzlTVpvip7ybgVKAVOE9VX+qvjQFZJx3HeZewO5JFf0f/3AacskfZ5cDjqjoHeDx+D/BJolyTc4DziYyH/eJKzHGc3mj+rJOq+hSwZ+7I04Bfxq9/CZyeVX67RjwHVO2x9h7ElZjjOH3Ztwv7tbETPcBmoi2NEDnRr8/6XB05ONb7BnDHcfowAPeJ8SKyJOv9QlVdmGtlVVWRwa3AuRJzHKcvuSuxRlWdP8Cz14vIJFXdFD8uNsTlG4BpWZ+bGpclMqRKbOq4rfzvc24Pyr597zlmvco3wwPadLg90J0rbDeK2kbbDD2mzpZt+FB4uMo32mb3Wce+bcoa1k83ZQWP2+fcOcPu47hXwvUKO+w63aPttjIJ/skl222ZGB4R3eV2W2X19vdZvNN2iWidbNerfiN83UW77AgRHbanBJkEF4vCdrve5FfsqCrrzw1HpCh9tcyss6A2/NveUNRpdyJXFNi3iUIWAecC34//PphVfrGI3AW8H2jOeuw08ZmY4zi9EPLnjS8idwInEj121gFXESmve0Tkq8DbwJnxxx8mcq9YTeRi8eVc2nAl5jhOX3ryMxVT1bMN0UmBzypw0UDbcCXmOE5v9v3jZF5xJeY4Th9GyubuXHAl5jhOX1yJOY6TXvazDeD5ZOPmcVx7fdiVouo0OyFFa3M4IcX4l+2B7k5wDeioTEoskRDRYWfYvH7PZdebdT7/g2+ZstYpCf0vtxclkiJjtNaGZWPesdsq3xw28QPM/peVpuzJP7zPlE1+KnzOuo/ZEScOeNbuY9tEu96YtaaItnHh73PX1CS3Ersfk5+xXSUyJQmJX2rsn9qE/w7LGo+0+/H0WwcFy1s68hCzb3/LduQ4zrsPXxNzHCfduBJzHCe1KNDjSsxxnNTiC/uO46QdV2KO46QWBTLpcdkfUiUmQIFhze94LOxGATBqR/i/wtaEpB4V621Z+WY76cSL5/zYlL3/1n8Ilp/xc9uNYlR3QoSFVQkuFgmRJQpsKz8dNeFyy9UAoCDBnP7YK4eZslFtCZE2poRvrcq/mFXYdqjdx+o37O+s7MKN9jnvmRosb6+13UomPmu7czQdbEfTKGm2x7GkxVYK7dXh665+3axC1zvhCBeFO/MR51RBXYk5jpNm/HHScZzU4tZJx3FSj8/EHMdJNa7EHMdJLaqQsQ0pI41+lZiIjAKeAkrjz9+rqleJyEzgLmAcsBQ4R1UTA3x3j4Kmw8Iaftb9Rr554K2/DsdDH11vW2KSYsP3FNlWtQW3hS2QACXGBnBNMAh1fXK7KWt7sdqUVdQlbIY+vdmU9bxQFSwvtC85kZpJdltt79h5DKrPrguW1z01LVgOUPuCbXbdMcO+VeXmyabswG+sDgvumG3WafqkfS8WLy83ZS0z7O9sVINt8RTDEJiUA6jkA1vDgkdsq+uASNFMLBd7bAfwUVU9EpgLnCIixwLXAT9S1dlAE/DVfdZLx3GGlvxlAN/n9KvE4my8O+O3xfGhwEeBe+Py7Cy+juOkGo2sk7kcI4Cc1sREpJDokXE28FNgDbBdVXfPXXPK1Os4TgpQ0Dw5u4rIIcDdWUWzgO8CVcDXgC1x+RWq+vDetJGTElPVDDBXRKqAB4BDc21ARM4HzgcorLbXgBzHGUHkaduRqr5BtAy1ezK0gUiHfJloOeqHg21jQHsUVHU78EfgOKBKRHYrQTNTr6ouVNX5qjq/sNxeFHUcZ4SgGqVsy+UYGCcBa1TVzii9F/SrxERkQjwDQ0RGAycDq4iU2efij2Vn8XUcJ+3sm4X9s4A7s95fLCLLReRWEdnrxzTRfjoiIu8jWrgvJFJ696jqNSIyi8jFogZ4GfgbVe1IOldF9VSd+5FLgrKeBB+AnZPD5umdM+z/BB8/YZkpW/KTeaasq9zuhxpW8h2z7X6U1+2dG0hXhf29VL5p1yvdEe7Ljhm2ib90m91WYcJm8yRXhJrfht1iyjclnDCBpPj1LdPsVZHOseHvs2usfc2SSQgekOD60jHOrpe0Ody65078wotmnd/9YX6wvO6mH9GxPiH6QQ5UFo7XY0d/KqfPPrrr9reB7AQZC1V14Z6fE5ESYCNwuKrWi0htXE+Ba4FJqvqVvelvv2tiqroc6POrV9W3gAV706jjOCOZAc2yGlU1rFF780ngJVWtB9j9F0BEfgE8NOBuxuQjbofjOPsTuzeA59fF4myyHiVFZFKW7Azgtb3trm87chynFwpoHrcdiUg50Vr6BVnFPxCRuXFz6/aQDQhXYo7j9EbzGxRRVXcRbU/MLgsnoN0LXIk5jtMHHSHe+LngSsxxnL6kKDx1vy4WeW1MZAuw29FtPL1Ns8OF96M33o/epK0fB6qqnbAiB0Tk93F7udCoqqcMpr3BMqRKrFfDIktyNM16P7wf3o8R0o+RiLtYOI6TalyJOY6TaoZTifXZmjBMeD964/3ojfdjhDNsa2KO4zj5wB8nHcdJNcOixETkFBF5Q0RWi8jlw9GHuB/rRORVEVkmIkuGsN1bRaRBRF7LKqsRkcUi8mb8d59HkDT6cbWIbIjHZJmInDoE/ZgmIn8UkZUiskJELonLh3RMEvoxpGMiIqNE5AUReSXux7/E5TNF5Pn4d3N3HBnCUdUhPYhC+qwhClNbArwCHDbU/Yj7sg4YPwztngAcBbyWVfYD4PL49eXAdcPUj6uBbw7xeEwCjopfjwH+Ahw21GOS0I8hHRNAgIr4dTHwPHAscA9wVlz+c+DvhvJ7GqnHcMzEFgCrVfUtjVK83QWcNgz9GDZU9Slg2x7FpxHFbYMhSrxi9GPIUdVNqvpS/LqFKOjmFIZ4TBL6MaRohCfnyZHhUGJTgPVZ74czyYgCj4rI0jgXwHBSq6qb4tebgdph7EteIm7uDSIygyh+3fMM45js0Q8Y4jERkUIRWQY0AIvx5Dwm7/aF/Q+q6lFEAdsuEpEThrtDEP0nJlKww8HNwEFEyR02ATcMVcMiUgHcB1yqqjuyZUM5JoF+DPmYqGpGVecS5a9YwACS87zbGA4ltgHITgNtJhnZ16jqhvhvA1EGluGMVFu/O1Bc/LdhODqhqvXxD6gH+AVDNCYiUkykOO5Q1fvj4iEfk1A/hmtM4ra3M8DkPO82hkOJvQjMiS0tJUTJAxYNdSdEpFxExux+DXycQUSXzAOLiBKuwDAmXslnxM0BtCnALcAqVb0xSzSkY2L1Y6jHxJPzDJDhsCYApxJZftYAVw5TH2YRWUZfAVYMZT+IwvRuArqI1ja+ShQ07nHgTeAxoGaY+vEr4FVgOZESmTQE/fgg0aPicmBZfJw61GOS0I8hHRPgfUTJd5YTKczvZt2zLwCrgV8DpUN1z47kwz32HcdJNe/2hX3HcVKOKzHHcVKNKzHHcVKNKzHHcVKNKzHHcVKNKzHHcVKNKzHHcVKNKzHHcVLN/wUySZLfh+SmxgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imshow(b)\n",
    "colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0d164049",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# like moving window with the 0th dimension the IW no.\n",
    "a = a[np.newaxis,:,:]\n",
    "b = b[np.newaxis,:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "968cf5e9",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %%timeit \n",
    "c1 = fft_correlate_images(a,b,'circular',normalized_correlation=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "05e66484",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %%timeit\n",
    "c2 = fft_correlate_images(a,b,'linear',normalized_correlation=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a8f527d5",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %%timeit\n",
    "c3 = fft_correlate_images(a,b,'circular',normalized_correlation=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f0086113",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %%timeit\n",
    "c4 = fft_correlate_images(a,b,'linear',normalized_correlation=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3481b581",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-12.02031891 -10.95825492] 0.05293600647141261\n",
      "[ 0.03563358 -0.00205613] 22.98857744116232\n",
      "[-12.0226627  -10.95144842] 0.043785720871595046\n",
      "[-12.01144001 -10.94500473] 0.048564728262190826\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzkAAACyCAYAAABhjQGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABSpElEQVR4nO29edgcVZn3//kSeBKWkKBBVLLAC1E22VdXBsUBZl5wfqKAyuKLMgygOOpc4jgXzuj4jo7b6CsuKMiiggqomQHFDOLgBiYgWxKWsIWHxYiQEATykHD//qjqWOlUVVd119p9f66rr+fpqlN1TnX3Xed8677PfWRmOI7jOI7jOI7jDAsb1d0Ax3Ecx3Ecx3GcInGR4ziO4ziO4zjOUOEix3Ecx3Ecx3GcocJFjuM4juM4juM4Q4WLHMdxHMdxHMdxhgoXOY7jOI7jOI7jDBUucpxakXS+pOWSbs9Q9vOSbg5fd0laUUETHacxuL04juM4Tjbk6+Q4dSLptcBTwEVmtluO494D7GVm/6e0xjlOw3B7cRzHcZxsuCfHqRUzuw54PLpN0g6SfiLpRkm/kLRTzKHHAZdU0kjHaQhuL47jOI6TjY3rboDjxHAucKqZ3S3pAODLwCGdnZLmANsDP6upfY7TJNxeHMdxHKcLFzlOo5C0BfBK4PuSOpsndxU7FrjMzNZW2TbHaRpuL47jOI4Tj4scp2lsBKwwsz1TyhwLnF5Ncxyn0bi9OI7jOE4MPifHaRRm9iRwn6S3AChgj87+cL7BVsBvamqi4zQGtxfHcRzHiWcgkSPpMEl3Sloq6ayiGuWMDpIuIRiAvVzSuKSTgbcDJ0u6BVgEHBU55FjgUsuZFlDSLEnXSlosaZGkM2PK7CTpN5JWS/pgnmMztsHtxRmIhtmLJH0x/D3fKmnv/q8stg1uL04r6PVbTUrnLmnPsM9ZFNrQMWW1wXGaQgZ7eX/Y99wq6ZpwXml0/5Zh//elnnX1m0Ja0iTgLuBQYBxYABxnZov7OqHjlIiklwAvMbObJE0FbgTeFP29SnoRMAd4E/CEmX0m67EZ6nd7cVpDRns5AngPcARwAPAFMzugoPrdXpxWkPe3Gk3nLullgIVJQ15KYGc7m9mKMtvgOHWR5bcq6S+AG8zsaUl/BxxsZsdE9n8B2Bp43MzOSKtvEE/O/sBSM7vXzCaAS1n/CaLjNAYze8TMbgr/XwUsAbbtKrPczBYAz+U9NgNuL05ryPibP4pgvR4zs+uB6aE4KgK3F6ct5P2trkvnbmZ3mdnd4f8PA8sJBm9lt8Fx6qLnb9XMrjWzp8O31wMzO/sk7QNsA/w0S2WDiJxtgQcj78fJP/BznMqRtB2wF3BDhce6vTitJOU3X+Zv2u3FaQuZf6tp6dwl7Q+MAfeU2QbHqZm8v9WTgR8DSNoI+CzwwZTy61F6djVJpwCnAGw0aWyfzTbv5yGFMwo89eRDj5lZ4g/kVQdPsRWPPx+7b/Ftzy0Cno1sOtfMzu0uF6bcvRx4XzhpOzODHJujjnX2MomN99l8463KqAabMlbKeeNYO0W9C9XIpGf7C9ktCj07kfuYJ9f8IdVWAF5z8BR7IsFeFlVgL2VTla2k0W1H/XyXo0Da/SbPZ9bv512VvWQgNp176AG9GDjRzOIbMSBp9jJIf5DlO8hzfj070bN8p0/p996d1ifV3R+UTef7SvqM9exEJnt5dcJ4rEBbQdI7gH2B14WbTgOuMrNxKdu4YhCR8xAwK/J+ZrhtPcKLOxdg6rSZts9B7xmgSmeY+Z+rz3ogbf+Kx5/nO/+1Tey+PeeMP2tm+6YdL2kTggHbt83sijxtG+TYkNz2Mm2TF9krZ7ylj6qSmdhpZu9CBbJyh+4lW5rLtHtW11r/2B3jmcv+5NEvp9oKwBOPP8/lV86I3bfT7EeKsJdMv+k+6XnuDWzlhUcXVHV2JnaetcG2sSUPxpR0op/Vep/R5vnPse74jMf+5PdfKdNe8tjBBuncJW0JXAl8JAz77If89tKjb8naV/S6b5XR56zcYfJ69+vu972OzULd/UFZdL6vpO9l7I7xTP3Lisef53tXbqiDdp39cK++JZO9SHoD8BHgdWbW+TIOAl4j6TRgC2BM0lNmlphoYxCRswCYK2n7sIHHAm8b4HyOUxoKZP95wBIz+1xVx0ao3V5c4KTTaW9dndvETjNzCZ0yyfibnwecIelSgsQDK83skYKaULi9TOw8qxIBssFAvOR6spDUlqraWlQdDRWQmX6riknnLmkM+AHB3LbLym5DHsbuGO/ZZ9R1v+oWOJ2/ve7defqkuvuDMmhI/9LztyppL+BrwGFmtryz3czeHilzErBvmsCBAUSOma2RdAZwNTAJON/MFvV7PscpmVcBxwO3Sbo53PaPwGwAM/uqpBcDC4EtgeclvQ/YBdg97lgzuypr5XXbS5UCp23ipps8TwWHmJ72AlxFkFltKfA08M6iKi/aXjqD+aqETrTOvPUlejwGbEv0XN0CqcrPZdhI+q1K+hiw0MzmhUXj0rm/FXgt8MJw0AZwkpndXEQb+r6okCSh0z1Y7i6T5C1I6hv6vd92n6+Me3dcm5vQP7RVhGW0l08TeGq+H4alLTOzI/upb6A5OeEgL/NAz3Hqwsx+CaQGcZrZo0SyeEToeWzGNtRiLy5w8lNXB9L5rup+4pbRXoyu0JuC2zAU/UseAVGm+GiDkGlDG+OI+62a2dld7/855rhvAd8qqw1FEBUscfeluP6lyD4ni0cpSpLQKbJvavuDsLTPc2KnmfBoufX3shcze0OGc1wAXNCr3ECLgTqO02xc4AxGXddUdWjhsDO25MF1r7YxsfOsxLk/3a80ovub/DkkXa9TL0kPXsbuGM/8UKYfYZB0L0y7Nyd5XzqvQWmCwCnysxxmSs+u5jhOPQybwFk1Z33HwtQHqsmCU6dXp26PjjM4/Xpx4vbnETK9wt7GljzYV0hdWjsHEU8eQtdusnpcqrqPpnlbOtvb/mDOhU5v3JMzAE+8vLo0vI6Th2EXOEnbyqSODnHUOqS2UWQoWdZyWctm8fDkFTdZhFiRuEenuUzsNDNxLk7VrJqj2P5g5Q6TU+/bZQuuTv1xL6caXOT0SUfguNBxmkZVA+MqbtZJnVfW/UVTl9BxsdNcihIRdZGljXnERtVCx8Pb6qVJ96ekvqAOodOrr8jbl5QlIIc9WsBFjuMMEVUKnLLJI16qFDs+T6d52JSx0ge6aYPpNoiZKnGPznCTlGyg7ntUvyHMeYVOUX1A1vMMuxApk9aLnCdePlaLN2WrOyfWvaJtcZy6GBaBM4hgqVLoeMjB8NMRNtFBdlahk0f41JUUIe766qDu+p3s1C1kOkSFydQHLFXgZBExdQmdLJT5mTfl+yyL1iYe6BYU0fdR4VE1T7x8rNb6ndFkGAROUQKlc54qEhNUnUrUkxGkU9Rk9X6TAHRv6z5P1sQBVXilkraNUsrqUaeoOTXd/UJV98Rp96wutE/Ke74iktJk7UOGXYyURSs9Ob08JnV4d9yL49SFC5zqzhlH1R4d7+zSGVQg5EkCkHd/Vq9JmeIg6xyXpnh5stDJFNeW9taNTRnra32buMU94+5/Vd4TewmEvAKkH8HS1AVDnRaKnG4xsWqH51m1w/M1tcZx6qXtAqfsuTRVzdUZtvA1SedLWi7p9oT90yT9p6RbJC2S9M6q21gG0QHyIN9pEQNt94Jkp25hI+kwSXdKWirprIQyb5W0OLSX73Tt21LSuKQvVdPiZJL6lF59Tfe9tglCJ01opNm3C5Ry6WUvkt4f2sqtkq6RNCey7yeSVkj6ryx1tUrkRAXOSTf+jF2fv4sttl/JFtuvZNUOz7Pfg3dz0o0/iy1fNh6i5lRNFQKnzMF7lVnRhtWrUyIXAIel7D8dWGxmewAHA5+V1Ah3dr/iIDpQPupPv2KXtUt57DXPMf4XU1i5w2T2fOIejll2XeIxadv6bV+WkLi8ZF04NO2VtZ48beqXfttXFJImAecAhwO7AMdJ2qWrzFzgw8CrzGxX4H1dp/k4cB0VoWfzjVfi+prOve6EW3/GLmuXsnr2BKtnT7Bqjtjn4aWccOvPKhc6WRb77O7Tkvq4ohYO7cWoCaos9gL8DtjXzHYHLgP+PbLv08DxWetrxZycOLGyYK+ZfOmii/juHvtzz34v4olbJ/Opn17Bh954QubzFS1MXOg4VVGVwCmDqte36a637Lk6dS0eWiRmdp2k7dKKAFMlCdgCeBxYU0XbksizEGaUblGycofJLFw7m89ffhEf3+0IfjhnHw647wH+6c5L+deXH5vrXEUQXbQzS91liJAij007Z90emT7ZH1hqZvcCSLoUOApYHCnzbuAcM3sCwMyWd3ZI2gfYBvgJsG9Vjc5Kr75m4a6BrXxst7/id3vOZtslK/m3//kOH35dMBaret4ixN97s6Z07j62+32efrHXHJ829xED0NNezOzaSPnrgXdE9l0j6eCslTVe5MQJnFU7PM/t22/Nd/fYn3d/6Doee8fmbHHxBGd++K1cO3UHpt6z/vFJGdA8ScDoIGkWcBFBZ2LAuWb2ha4yOwHfBPYGPmJmn4nsOwz4AjAJ+IaZfbKqtnfjAqeYNlQhdoa4E/sSMA94GJgKHGNmtcUN69kJ2Dz4v3ugnDT4jxtQr9xhMqtmi+vmzOHHr92NT3zgR+z1Vw/y1z+6jTOPPYHFk3Yc6DstwgsTR3f2t7aGurVU6GwLRD/wceCArjIvA5D0K4I+5J/N7CeSNgI+SzCIe0MFbV3H2B3jqX1J2r5O/7Bqjrhu9hw+tttf8X8/8UNuOmY2r7jkIf7uPe/gps3nrrvH1nkvzNuX9WprZ1/0vL3K+5yd9chiL1FOBn7cb2WNEzlZQ8zeOPtOXrnHUh67bXNe/MVVzD9lFxbsvh3cF4igqfdki8RzoTMyrAE+YGY3SZoK3ChpvplFn7Y9DrwXeFP0wIh79VACg1wgaV7XsZXQVoHTBHHTTRVip87O/Ym1m3HZk3sn7L1yhqSFkQ3nmtm5OU7/l8DNwCHADsB8Sb8wsyf7a235ZB08r54zwXYz/8ADL57BMyeMceIXrufCdxzIdQfPYcYv+j9/WZ6TKjKlVUldQqdke9kYmEsQ2jkTuE7SKwjEzVVmNh44RaslSejk7Wd+t+dsnj5hjNd8YSm/OHVHHjp0GqvHJ4CxWoVOv31ZlrZGxU6vuT1tEjRZM+w9sXYzfvDkXjF7Hh7UVtYh6R0E3s3X9XM8NFDkZF135qfLXs4OC5bz7m9dx6PvncpBF9/DfrPv59qpe2QWON31OcOLmT0CPBL+v0rSEoInClEX6XJguaS/6jo8SzhC6ZQtcIbZe5NG2WKnoeFrj5nZIKEx7wQ+aWYGLJV0H7AT8NtCWtcnRQyOJz8wxv1szZxlv2TTiya48B0H8tc/uo1rNt6LxZN2XK9sVEQkDc7rFBpt9erkEToVCbo0e3kIiDZ2Zrgtyjhwg5k9B9wn6S4C0XMQ8BpJpxGEfY5JesrMYpMXNJW9bl7GZhdN8ItTd+QVlzzEttNW8vvNX7rePbUtAid6fJb5PVn3N+z+n8jETjPh0YFO0atvyWIvSHoD8BHgdWbW94fXOJETpSNAusXO1Hs2Yte7/8AxF/2Wr3/2tXxlm9ew34vu59//7Qo+9MYtWDBr7gbnSHrvtIc/rt2Ci594ZcLe72V+ehDONdgLuCFj1Xndq4XjAqd8Vs3R0Hp1SmAZ8HrgF5K2AV4O3Ftvk9LJEq4WTDaGbbd7kMMvv52PfPQofjhjH67ZeC+++O0L+fgux3HzVjvEni9pWzd55gtloZcgaILQ6eeao9fV8DC9BcBcSdsTDNaOBd7WVeaHwHHANyXNIAhfu9fM3t4pIOkkgsnWjRc4nRCsqQ8YB9z7AGdffiX/+NE3BXNypq3knP8XzMm58aU79j5ZSyk6DG6E6GkvkvYCvgYcFp2/1g+NFjkdtrpzYgOhs9/vxjnjhBO4fZuteeq+aVw7dQ8+9MYt2HX5g+tEjguaZlBWoocuMj2ZlrQFcDnwviaH1kQpU+C0be2b1bOD39DkZeUk8hpRr84GSLqEILRmhqRx4KPAJgBm9lWCbFAXSLoNEPAhM3uspuYC+b04aeLg1dfdx9+/+QRu+NNcZiwxHrp/Fh/f5Thevuohbt5qh0IyuHXeZz3XIIP4vMf281lmPU8ez0veML06MLM1ks4AriaYb3O+mS2S9DFgoZnNC/e9UdJiYC3wD2b2x/paPTgdobPvomX8/ZtP4LoZc2Acfr/5S/nw605g18eWrRM5/d7v+vWEFNWvJd2vk84f7fPi+hAXOpnt5dMEns3vh2Gcy8zsSABJvyCIGtgi7JtONrOrk+prhciBDb06F+xzCABP3ffn+TcLZs11gdNg6p7/JGkTAoHzbTO7IsehmdyrZdBGgVOmuOl+X6bYGWWvjpkd12P/w8AbK2rOwCQNlpOEzndnvxbuh6mTbN33dPNWOxQqcPo5tgqh0087BzmmSO9LXd4cM7sKuKpr29mR/w14f/hKOscFBKnbK6OI/uWi3YOx2ORlf35AdONLdxzYi9PdP2W9Z5bRr0XFThaB03lfdnKbtpLBXhKTcJjZa/LU1RqR06HbqxM3/8YFjtNNmOr2PGCJmX0u5+FZwhEKx6aUt+xIW5ILdIubuP3u1XF60e/gN/rd1B3y5WSjCeF5o0J00J90n6zq/lbFmjxF1NHETHPD3Ae1ajHQDlvdORErZJK2O82g5u/mVQQLSB0i6ebwdYSkUyWdCiDpxaH78/3AP4UrUG9pZmuAjnt1CfA9M1tU14UMShsETmdhuaLL9kPZ84qGaAHRWum1EGSSt2HQsKms1OEB6qf+ohYbrWNhzg4TO89qTDjbsBO3aGbawpxZs3floQn30KZ7bdI+oyZ8fmXROk9OFBc07SDue6o6dM3MfkkwdyCtzKMEoWhx+zZwr7aRtgicQY4rw7NTRfgaDPcTtapIm2uTFCIVd0wVC14WVUcT1pdJCwmEetJau1enOoq+d3WHhiWdv0kD9KkPWGXrsBVN00Oo+6WVnhyn3XTCDbOuieQUQ9MFTlEembI8O6vmyL06LaHXwDZNEJTtgeicu8mD77IEXh3XXLf4aypVrLmWRFZvzqApnNPo3M+7X0XQNIGT9XMaxv7HRY7jjABtEDhFU1YImwuddpBX6PQzCO83LKrqwX6TBZUzWkzsNDOXwGqbd6GtAmdYaXW4mtNOPMywWposcMqcS9M5f1nha1BuUoK2de5tZNBwpirEQ5b1YqL7ip7D0+t8dQmoJoTotZmxO8Zr9ebkocp7YdmhyVmvpbvfTsvsFnfeURc3HVzkOM4QU/SNrunem7R62jZXx4WO0yGPkCh68B+te5B5S3kW9uxnwVDHyUN0/kxV9Ctwkrbl2Z+FYexvXOQ4zhAyyt6btHrb5tVxodMu6phcH0dW708/5+2HqBCJ857FLZKatR1pHqe6v4c2MwrphqsUOoMInKLbEVfHMH2vUXxOjuMMGS5w0utv21ydOsIOJJ0vabmk21PKHBymYl8k6X+qbF8Wsj7hL8oT0D2QH/Rcg56jzhTOvSjD+9J9vVVeu6TDJN0paamks2L2nyTpD5HlC94V2Tdb0k8lLZG0WNJ2lTU8gVEKdYp7OJU1CUE0VXacSEjbF0cVAqfzt05Rk8FeXivpJklrJB3dte/fwz5niaQvhmsgJuIix3GGCBc42XCh05MLgMOSdkqaDnwZONLMdgXeUk2zstPmyf1ZBEqbw7TKzl5XscCZBJwDHA7sAhwnaZeYot81sz3D1zci2y8CPm1mOwP7A8tLbzTp2dXa8FS/yHtikhe+l9iJmzeTV9jUSR1tzWgvy4CTgO90HftKgjUPdwd2A/YDXpdWn4scxxkSyph/M4wCp0OZqabLoEqhY2bXAY+nFHkbcIWZLQvLVzIwy0uWRSnLSpdcJh2B0xah08vD0lSPU0b2B5aa2b1mNgFcChyV5cBwcLexmc0HMLOnzOzp8poaYFN6h+0meSaaQN57YZYU0WnhxnmETr9UcX9viIeup72Y2f1mdivwfNexBkwBxoDJwCbA79Mq8zk5jjMENDXBQJEiYruZf1jv/f3jWxdy3jLm6rR1QbgcvAzYRNLPganAF8zsonqblE4Vk9TbMlivY/5Kr/O35bOLYVsg2vhx4ICYcm+W9FrgLuDvzexBAjtaIekKYHvgv4GzzGxtyW3OREfUNGRwDPQncLKSNkcnLcmML+aci6z2sgFm9htJ1wKPECzu/iUzW5J2jIscx2k5oyhwotuKEDtlZWArOvtankQET66ZwvxHd0rYe+UMSQsjG841s3NzNGVjYB/g9cCmwG8kXW9md+U4R+UMQzauQdrf9msvk5Lt5T+BS8xstaS/BS4EDiGwo9cAexGE6HyXIEznvLztL5OiBu/RvirunHHhX0n7shDXlw1yT+51bN5EMXWIx7Q2rmvPz3uf58k1U7j60Z1j9gxsK4lI2hHYGejEWs6X9Boz+0XSMR6uNiRsevtDdTeh0UiaJenacGLnIklnxpRROJFtqaRbJe0d2ZdrsltVDLvA2W7mH2IFTneZomhD+FpB3/ljZrZv5JW3ExoHrjazP5nZY8B1wB5FNKxsWuwxcOojzV4eAqLqcWa4bR1m9kcz64wsv0HwgAACO7o5DN1ZA/wQ2JsWMnbHeN1NqIQi7ucrd5jcKO9YwfTqW3raSwp/A1wfhnU+BfwYOCjtABc5Q4QLnVTWAB8ws12AA4HTYya7HQ7MDV+nAF+B/ia7tZEmCpw8ZYsSO2XM1al6PYYK+BHwakkbS9qMINwgNWzAKYboPKO8gs0FXiksAOZK2l7SGHAsMC9aQNJLIm+P5M+2sgCYLqnjjj4EWFxye0shz8KiWbw4SdvyMPUB28Dz0suLk8XLM8j9vMnipqJwu572ksIy4HVhv7MJwTgstd9xkTNkuNCJx8weMbObwv9XERjGtl3FjgIusoDrCTqfl9DHZLcqKPJm2SSBM4hgabJXp01CR9IlwG+Al0sal3SypFMlnQoQxkH/BLgV+C3wDTNLTDftNAcXOsUSemDOAK4m6Fe+Z2aLJH1M0pFhsfeGkQC3AO8lCEkjnHvzQeAaSbcRzDP4etXXUDdlD/w7YqfOOZJt8d6ULXSy2Iuk/SSNE2Tt/JqkReHhlwH3ALcBtwC3mNl/ptXnc3KGkE1vf4hndusevzsdwnUI9gJu6NoVNyFu234mu5XNMAucos7RxLk6Rc/RKQszOy5DmU8Dn66gOaUStzjlsNO9gOaoXX/RmNlVwFVd286O/P9h4MMJx84niBIYeoZpYn7SvTxuzksbxE2VZLCXBfx53k20zFrgb/PU5SJnADpekyYKiiqFTp66op6mvO17cmIKP1328qTdmSa7SdoCuBx4n5k9maXefia7lYkLnOzna2IGtrYInbbiA/fsDEMiBscpmrQsa1lwUdMceoarxa18LekFkuZLujv8u1W5zWwe0cH6prc/1LgwsSoETvS6O/8nfQ5x+wr+zHpOpA5jOC8Hvm1mV8ScI2lCXObJbmXbyzAKnCLn08SduyiKDF9rU+hamZRhL1nmqnT2V71wZBMZ9etvE8M6Hkvy8AyT56dJpI0jhk2gZZmTcwEbrnx9FnCNmc0Frgnfjzy9Bvpl0hE1z+y2bekCp9c1Rj+HpgjAMBvaecASM/tcQrF5wAlhlrUDgZVm9gj5JrtdQEn2MqwCp2xc6DSaC6ipf/HBvdNCLmBIx2PT7lm9nqhxgVM9wyZwIEO4mpldF85hiHIUcHD4/4UEWbU/VGTDepEn7Kl7kF2VlyOOMuuu47o6dTZByPTgVcDxwG2Sbg63/SMwG8DMvkoQI3oEsBR4GnhnWO4ygqw3txEkIfhJ0mS3ptpLlFESONG6PHStebTBXhynKYyCvTRF3PQKWav7vl30Qq3DKHCg/zk524RPuAEeBbbp5yR5542khUIlnSfumLhtTZy/UjZ55hT1EjHP7LZto4WOmf2SIGlAWhkDTo/ZnnuyWxcD20tRN6BRFDjddTYtIcGoC50YCulfHGdEcHspiUHn5pRBtwjsV+xEEyQMq8CBAhIPmJlJSuyhJZ1CsOYIk6dMT5yXkTbQzjJ4jhMPeQbdg0yIz0vdQifvd5D1c4wen+WYpoi9KsllL5OnA8MlcOoQN3FtaJpXx4VOPGn2ErWVKRttUWm7HKeJZLWXTt8yynT6w6T7btPETYfOeKAIgTLM4qZDv+vk/L6zuFX4d3lSQTM7tzMZfPLa5MFAEV6A7mQAg5yn15yTQek1Qb+suSxFn7OOsLwW0pe9jG2yeWENGBaBc+iL7+DQF98x8HmaOE+nqZ1qDWSyl/VsZaNNK22g4zSI/PZSYN9SBUUPxqP32lVzlPo+jiY8kGrLujt1068nZx5wIvDJ8O+PimhM08Kd4jwcRbYx6tFJEz3dbUg7X9EMcs60ELYRE0F920sRN7FhEjjR/+c/utNA53OPTmMppX9xnCFlqO2lTIGTZXs3affn6Dn8Pt4MsqSQ3mDlawJjOlTS3cAbwve56GQBK3qwW4YHpEzxlbW9vcqUcc1xSQbyfmdxZYdZ4BRpL2unDC5OhkHgJHlvivLoFOXVGSaPTlyq2oRy+0laI+noPusppX9xhpcmrusj6TBJd0paKikxu5mkN0sySfuG7zeRdKGk2yQtkRS7YGjk+ELsRc9OMHbHeJ5LdFKI8/7UfR/vtClvO6pody97kfRaSTfF9S2S1kq6OXzN61VXluxqSStfv77XscNEE7xMVc3lyZKYIU9WtahHZ5gFDjTLXoZF4GTZ3xSvzhB5dC4AvgRclFRA0iTgU8BP+62kSfbiNJ+GCpxJwDnAocA4sEDSPDNb3FVuKnAmcENk81uAyWb2CkmbAYslXWJm98fV5fbSP0X0h+6dGZyM9rIMOAn4YMwpnjGzPbPWN3DigX5pgmgYlDrSJ3cLnaLD5/KSN313nZntnPzUJXDyemmaFL5WZIrpukhIVdvNewgW192v/BY5TmPXNtofWGpm9wJIupQgrfPirnIfJ3go8A+RbQZsLmljYFNgAniy9BYDEzvNrKKaSugWMFExUoVnolNfleFqvRIn9JsZrgIh19NeOiJf0vODVtZv4oFaKCvErd+2RP+vo11NF4pZ29eUBUOHiSJu7G0ROIMeF6VJoWt1hzukIWlb4G+Ar9TdFsepmW2BqPoaD7etQ9LewCwzu7Lr2MuAPwGdRac/Y2aPl9hWAGxKex7C9Fo3J+4+2W+o1qB0BEIZQiF6Td2JEtLa00DvU0976cEUSQslXS/pTb0K1+bJyUtSuBRUP9hPEzNVeHeadL1JuGipj1EUONHjh8mjM0jY2sTExmnXMUPSwsj7c83s3Byn/w/gQ2b2vNRcMea0h2goWh0em7LsRdJGwOcIwm+62R9YC7wU2Ar4haT/7jzlbgvR5ABVLubZhAdB3ffosgRO05iY2Jhl4zPidg3at/Rijpk9JOl/AT+TdJuZ3ZNUuBUip9cgu8qFKLMO+Ju+OGZZjOI1N4lRFjjd5xlE7BS1cGjdQieFx8xs3wGO3xe4NBQ4M4AjJK0xsx8W0ThntOieazOx86ymhaal2ctDQPQCZobbOkwFdgN+HtrLi4F5ko4E3gb8xMyeA5ZL+hWBbTVe5CRlPYsuMtkm+lngudPf9godSzqum37v80X2EQV8d736ll72koqZPRT+vVfSz4G9gESR04pwtaoziznr459vO2irwClq7Zu48w5KEeFrRWVdaxJmtr2ZbWdm2xGE3JzmAsfpl25B0zCB04sFwFxJ20saA44lSOsMgJmtNLMZEXu5HjjSzBYShKgdAiBpc+BAoPibYYGUvT5L5/zRVxxFezcmLxsr5IFUWorqMkPoijpvBWvvpNpLGpK2kjQ5/H8G8Co2nPu2Hq0QORA/0K56LkeesK1hEwZtvx5JsyRdK2mxpEWSzowpI0lfDNMa3hrGUXf2zZb00zDN5+IMk7Irpc0Cp0yGRehUHa4Ql6pW0qmSTq20Ic7IMLbkwXWvNmFma4AzgKuBJcD3zGyRpI+F3po0zgG2kLSIYPD3TTO7tdwW56eX4Igr3289TaXbU5J2T+6ePzPo/XtQD1GTyGIv4dIE4wTZB78W2gfAzsBCSbcA1wKf7M5i2E0rwtU6dDKL1THgbkKyg7rJkga6wWF6a4APmNlNYSrPGyXN7zKQw4G54esAgknVB4T7LgI+YWbzJW0BDJz1oyhc4GSrZ9DwtbpD16pMK52Sqjau7EklNsUZYRoYthaLmV0FXNW17eyEsgdH/n+KYCDXSKoUHXnqqiuhQNlkvcfHZXOLniNapon0shczW0AQxtZ93K+BV+SpqzWenA4NHUCPFG38DszsETO7Kfx/FcEThG61dhRwkQVcD0yX9BJJuwAbm9n88PinzOzpKttfJsMscIqss4iFQ9vm0XGcumni2jjDThEhaXmOr0PgDCICyrwPp83X6c6WlnYNdWSXayKtEzlOM2ij0OkQhprtxfqLskFyasOXASskXSHpd5I+HS5oVTuD3sRGReA0oe4OwzhHx3Ecpx/q8uDkOVeRoWd56stCmtBpsjenKloVrtYmGhy2VRhlLkwax9rVk3jqvmlJuzOlLQxDzS4H3mdmWRdd2xh4DYEwWgZ8lyAd6HkZjy+FtgmcQQXG0VvetO7/y57cO6Vktnb0G75Wd+iaP51zRoU2hKoNI50MW0V4c5KydTV5/k1Won1oGQtA5wlfS1sQdZRxkVMiZQmdfoRFWev3NEjI9UyJK2kTAoHzbTO7IqZIUmrDjYGbIyv0/pAgC06tImcQ6koT3S9RgRN9P6jYqVPoOI7jNJkixE7eY5MG9nke7vSTDjov3X1oEcsFDEJH6Li4WZ9aw9VGYTJ/kdf4zG7bxi6KmlRHZ190/yh85nEoWKDgPGCJmX0uodg84IQwy9qBwEoze4Qg6810SZ1R7SH0SFtYNoM8zW9TiNrRW960gcDp3j8Ig3iX6p6f4zhFM7HzrHUvZ3ToJUSm3bO6NYt8Ru+rZd1jk85bRn15PgsXOBtSm8hp02B70LYWca1ZFkQt6lxDyquA44FDJN0cvo7oSot7FcEibEuBrwOnAZjZWuCDwDWSbgMU7q+FURI4RZZLwoWO4zijThaPS5rQKVMIZe3z4u6nq2dPFHqf7XWutt/T27iQaxq1hat1z+doKkW1sew5OknrCMW1v0EhZpVhZr8kECdpZQw4PWHffGD3EprWCpoqcLrL9xu+5qFrjhPMgXEvzuhRxMC2qDk2edajyUOc+MgbXpZFwNQZsjYowyZwoOZwtaYPtosWYf2eL+24XguiNv0zdvJRhxen6QKnqGPr8ui0/cmfM1y0cUFOpzqGaSDc8fJk8fbUJXCS5id5iuhseArpBMryMnXOG50vkzavJo2sAiZabhDRM0hbncGpK0ytH+oQOEWcw4VOgKTzJS2XdHvC/rdLulXSbZJ+LWmPqtvoOE1B0mGS7pS0VNJZKeXeLMkk7RvZ9uHwuDsl/WU1LW4HVQ/ikwRPkzw43Z9JG4VOL3uR9FpJN0laI+noyPY9Jf1G0qKw/zmmV10ucrrIOogfVCxk3d69reO5yVt/P8dE2+DCpl7aNA+nH6HQK8FAP+frFxc6AFwAHJay/z7gdWb2CuDjwAbp2p1y8bCyZhCumXYOcDiwC3BcuIB0d7mpwJlE1mcLyx0L7Epgb19uwhpsRSwGOih5+7yi75+d8/U67+RlY6kCp3sBzzxEj0vz3ORd8yeJKr7zjPayjGCZju90bX8aOMHMOvbyH5Kmp9XnIidC1oF8RyxUHQpWV+iZh7y1l7YInDJoo9BpCmZ2HfB4yv5fm9kT4dvrCdKtOxXiIWWNYX9gqZnda2YTwKXAUTHlPg58Cng2su0o4FIzW21m9xEkvdm/7AanER3odsROZ1tV4qeN3okkBgkrKzosLct5Kvh+e9qLmd1vZrcCz3dtv8vM7g7/fxhYDqROiK09hXQRYVtFtSUL3QP+Ij0qSckDnNGmjht+0wTOC77yFJv+ev1Y8E1/vZoXfOWpQs4fx6CLl/ZDw7w5WTkZ+HHdjXCcmtgWiCrO8XDbOiTtDcwysyvzHtsU6vbs1EGW+3Fb7tnRcUTNIrKQ37yk/YEx4J60co1MIV234MlLEXNjep3DxY6TlyoTDZTtwXlmj0146Wkr1gmdTX+9mpeetoJn9tik0Hq66VfoNCFsTRNaF0rR/QJmSFoYeZ3SVx3SXxCInA8V0mjHqQlNiMkPjMW+GMBeJG0EfA74QFltL4omCJl+BuBtERplkva5lSFqkvoXCupbUuuWXgJcDLzTzJ5PK1tLCul+1nTpNciPm7tSdHuyZDHLc22dFM9tEDhFpvyOnqcJ19Zk+r05DZPAAXjmlZN5+MvTeelpK1hx/GZMv/hpHv7ydJ55ZbZO+egtb+o7vXS/NDy19GNmtm/vYslI2h34BnC4mf2xmGY5TiNJs5eHgOgEqZnhtg5Tgd2AnwdrUvNiYJ6kIzMcWxl1C5ykvq7TlxU9sb9zviIE0urZE321r6x02VnOuWqOylo8tFffMtBvXtKWwJXAR8zs+l7lh2JOziCD76IX0UwKR8sTolYXdXjP2uKtc+rnmVdOZsXxmzHjC0+x4vjNMgucNtL0J5OSZgNXAMeb2V11t8dxamQBMFfS9pLGCBIJzOvsNLOVZjbDzLYzs+0I5rAdaWYLw3LHSposaXtgLvDbqi+gToEzyJyTpNTPvQRHlWvZJLWxJIHRBlLtJY2w/A+Ai8zssizH1LYYaJFEhULRoiXpmDZ4X/ol6RrbfE2jxuRlY30NlO8f3zq3N2f+ozvl9uZc9uTeub05m/56NdMvfprHztyC6Rc/zdMHjWUWOlV7cYCBvDh1Lygn6RLgYILQg3Hgo8AmAGb2VeBs4IUE2aAA1gzqGXKcNmJmaySdAVwNTALON7NFkj4GLDSzxAFcWO57wGJgDXC6ma2tpOEV06+QSerLsmQ9iysXvbcW9TCp7vt1mmBqmpjKYi+S9iMQM1sB/1vSv4QZ1d4KvBZ4oaSTwlOeZGY3J9VXi8jpJ/QpS1hX59zdxxVRf1x74urLe3y0TUUyaPs656hC2FRVTxuZ+oDVPUmwNPIInc4cnE6I2tMHja33vlc9/TL/0Z36PrZf6u4wAczsuB773wW8q6LmOE6jMbOrgKu6tp2dUPbgrvefAD5RWuN60IZsaXmFSTR8LCp2ihI4/dyje4W09fsZNU3EZKGXvZjZAmIydprZt4Bv5alrKMLV4kgLvaozRCqpXU0N26qqXU29/jbT72C5Hw9Ev2IgqwDZ9Jbn1hM0nTk6m97yXCHnj2MQgdPguThOgUzsPGvda9DzOE7TiKYwHiQNchEkhX1lpQ6Bk1ZPEZ9tGwVO1QylyBn1AXM/XrJh/8wkzZJ0raTF4Wq5Z8aUkaQvhqvw3hqm/Yzu31LSuKQvVdfy9tzIyhQ6j//dFht4bJ555WQe/7st+qqzTNocpjZKFC1Mks7XSwR19rnQcaomzYtTRgRBv0KlyGOq9uCURVvGBXVTm8jx8KT1KeLzyCNW6kzTXdN3vwb4gJntAhwInB6zyu7hBBM/5wKnAF/p2v9x4LqyG1okVXpzoHyPTlXnHNUwtVGi6MU0+z1f5zhf3NNpElMfsHWvooikGa6EIhO45G13P/OIsuICJzu1enLyDHbLKls1nexr3a9B6EfcZGlfmVT9HZnZI2Z2U/j/KmAJGy5AdRRB1g4LUxNOD/OxI2kfYBvgpxU2ex2D3NRGUeh4mJpTNlFR0kugDLrfcepkELEzSN81aIha5xyDUoTAKQoXOPmoPVytyYIkD3V5RcrMJle22Knru5e0HbAXcEPXrtiVeMPF3D4LfLCSBo4wRQidNgoc9+K0k7ElD+YSMEXN4Wk6o3CNTrk0JY1+3jk4vdo9yHVVsa7OsNGIFNK9sp01WQh1t7vKTGFZREvexUn7PV+/15wn092k1TD1nkRdPkPSwsj7c83s3O5CkrYALgfeZ2ZPZmzmacBVZjYepsuthUEyrVWZUhr6SytdBG0LUXOGi7iB/diSB2O3T+w8K5MHJ3psHR6fTv1Z6462N+s1OtVT5do4efuuQYTA7JmPrft/2fiMdefL+yCpqAQDecrH1Znkuel8nnGfq3t7/kztnpwOaYPkJk6KryKkq0004Dt6zMz2jbziBM4mBALn22Z2Rcw5klbiPQg4Q9L9wGeAEyR9svArGDKqDlurYy2cDu7FcdLwgX59uEdpQ9ZOqfZhXZUeh46wiSOPGMmaqrofQZb1nt9L4CThHp4/0whPThbyekiKWAsn7pxp25rmxamaJq93o8AFcx6wxMw+l1BsHoGYuRQ4AFhpZo8Ab4+c5yRgXzM7q+Qmx9Imb84g5F0sdFCB4/NwnLIZNDFBXdRZf14vktM8+umvevVVSQt9dogTOlVldYsjTtT02x9npfO5j7pXp1JPzvObbpK6v9cAua6BfVavTRkD/GgWtDwZ0cpqSxFlauJVwPHAIZJuDl9HSDpV0qlhmauAe4GlwNcJwtQcqk9CANmFS5sFTlO9OJLOl7Rc0u0J+1PTrVeNTRmr/Yl93fUPG91hb01G0mGS7gztYYMHYGE/c1vY7/yyk9lT0qGSbgz33SjpkCrau7ZBt528A/3oPbPqbG39kNa+ftqfVzQ20auTwV4mS/puuP+GcB41ksYkfTO0l1skHdyrrso9OZ1BcL+D8DzegiK8OUWJhWg7ks5ZhuepwaKjUszsl0CqtZuZAaf3KHMBcEFhDauYsp8exTHI/Jy8Hp281DkPp+Gd8wXAl4CLEvZH060fQJBu/YBKWpaCzwFpP00XNN1ImgScAxxKkKxmgaR5ZrY4Uuw7ZvbVsPyRwOeAw4DHgP9tZg9L2g24mg2zfpZCFU/6ew2wO31R1jkzSWW6t7ctaUFSuVVzNHSemIz2cjLwhJntKOlY4FPAMcC7AczsFZJeBPxY0n5m9nxSfY2Zk9Mh6+T3XoP3IrKdFZXaOS45QRXio59FQeuo18lHXTe9QbwWZXl06kw0MMxhamZ2HfB4SpHEdOuOUwYNFj/7A0vN7F4zmwAuJbCPdXQludkcsHD778zs4XD7ImBTSdVlBCiJVXM0sAdhkIdAHQ9JnQ+SGv4Qq0562kv4/sLw/8uA14dTDnYBfgZgZsuBFcC+aZX1FDlJK8VLeoGk+ZLuDv9ulf0akxkkJXIRC1xWtU4MbCh2iqyz38/Ahc5gVG0v/dDGm2+cmGmzwBnkO2jIk73YdOt5T9IGe8lKnAdpFL1KedYP6nV8N/0InQq+g0y2IOl0SfcA/w68N+Y8bwZuMrPVcZW0wVbyipu8YqQpHpoqiX6e/QrHhoWsZbGXdWXMbA2wEnghcAtwpKSNJW0P7MP6yaI2IIsnJ2ml+LOAa8xsLnBN+D4TRQ/mi/La1CUyyhI6/VK10GnCNRdI4fYSx6h5c2B9UTOqAicPkybWX7W8awXzGZIWRl6nVNKoDSnUXuoWFZ36s6yf0xTK8JAMev11fHaTJmDqMot9UYC9mNk5ZrYD8CHgn6L7JO1KEJLztymnKLxvKaofKcJzM0wUfY9v4meb1L9Qbt9yPoEoWgj8B/BrYG3aAT1FTspK8VF30oXAm7K0sIkD2qLb1I9IcKFT/zUXQdH2UhaD3ITrFjp1poqumwLFbc+U6z1ISreeiybbS78LWlY1QC9iwc3O8WWHgvXT1qTPsaawtTR7yWsLlxL5PUuaCfwAOMHM7kk6qGhbKVLgVEnRaaCLpO4QuYbQq2/JYi/rykjaGJgG/NHM1pjZ35vZnmZ2FDAduCutMbnm5Gj9leK3CdPrAjwKbNPr+I2eeS51f9UhTmnem7rDrYZl0D/KDGovvagzbKlOoTMIo+LFqYB5BOtFSdKB/Dndet8Mai9liYuiB9WdAX9R5y3yXGXTcqGTxAJgrqTtJY0BxxLYxzokzY28/Svg7nD7dOBK4Cwz+1XWCge1lUkFjf3r8jD0I15mz3xsvYVCs5Dnfl32vb2J3pw+6Wkv4fsTw/+PBn5mZiZpM0mbQ5CZEFjTlbBgAzJnV1PXSvGKrPweVh474gpdVacATJk0NfbcdQuKJLrbVbXwqDs7WpPXvWk6RdjL2OblhlbXkWmtwyAZ1wapcxDqFjhVilpJlwAHE4QejAMfBTYBCLNEXQUcQZBu/WngnQPWl9teorYyefL0QapPpUzPTBHrwHSO7edcY0sebHxGuu7ri25rAma2RtIZBJnRJgHnm9kiSR8DFppZZ/21NwDPAU/w5wHcGcCOwNmSzg63vTGcVB1Lk/qWQdZtq5rutXOSBE93uaz9ZNYMcaNORns5D7hY0lKCBDjHhoe/CLha0vME3p7je9WXSeQofqX430t6iZk9EmbViTXK0FV1LsDUaTM3ML46175JIqlNWQf9dYuTIolecz/XNYoiqSh72XzrWT1HtYN2MoMInUEXCa1S6LRd4FSNmR3XY3/PdOtZ6ddeorYybZMXpdpKvwKgaJI8EP0IjV5ejqznq0IwFCGmmiRsujGzqwiEf3Tb2ZH/z0w47l+Bf81aT5V9S1Y6D1+6+6GyBVAWUVHUg7w2C52GJKpZjwz28izwlpjj7gdenqeuLNnVklaKj7qTTgR+lKXCvItaxhHNgFZ0NrQsqamzUFR7miASuucL5VmnaNQo2l6aThvSKNcZHlcUTeyoiqAKe2nLopKDtC3u2H7PV1b4W5NFShtoet8SMwF9oPtWVlHRizjRkTdsLQ9NyQAX/R5GmSxzcmJXigc+CRwq6W7gDeH70kkbOFeV+jkrdQudMrPFNe2zbhCV28ugN7I6nzyVLUCKOP+oeXEqpjR7SRqsN1noFE3ea+0WhKP0WbWARo3FsjIsA+089/G6hE63yHQyhKv1WCn+9cU2J508HoSmhIsV1Za856lKgLjQWZ8m2UsehjFsbVgEzjB3WG21F8epmjbbSt1zd6oOW6uKYe4biiJXdrW6qWKQX4aHok1CIK6tw+C1SVpIrauMJH1R0lJJt0raO9y+p6TfhMfdKumY6q8gnbpvdoOKgaI9Ok0QOEVQ9/fqFEtZIVtJC5F6iJjTVor2gBd1viznKdt7796a7GTOrtZG0rwfvZIIFO0N6nW+OkVEd91tFzQJdBZSu0nSVOBGSfO70g8eDswNXwcAXwn/Pk2whsHdkl4aHnu1ma2o9hLSqTMJATTHo9MUgeNhak4ccdnCyjh/087lOHnp3EOT+qV+77HLxmesNy+nO6PaIJR133dR0x9DLXIGpUlhb/2Qpe2DCppOHb3OU7fAC9cReCT8f5WkzkJqUZFzFHBRmDnqeknTw6w1d0XO87Ck5cDWwIrSG14xbRc6wyRwvFPLjk0Za2yK4SQ6GceKPF9dFJEKe5C62/B9O/1TRrhZL2GTdg9Pak+RAsfv/8XQqnA1yJ9yup9U0VmPz0uvxUeLyDyXtz1VkCVjXZVismshtSjbAtHecjzcFj12f2AMSFyZuk6acGOsK3TNBY7ToS0T5n1wXgxFL7DqFE+ee1rcBPrJy8ZGxjvu9//iGHpPTpYBdq/BftSjU4QwiDtHXDujbSs6dXVRHpysdRUhYiY9a2x1Z+LTnBmSFkbenxuuC7Ae3Qup5ak/XIPgYuBEM3s+z7FtooinZlV7dFzgOB3i1mSJ8zTU6X2I0qQFOaPepbrXFGpD3U528oRSr5qj9dJQR4/t9E/93mt7hcBFy/Q6T/QcRQkwv/8XS+s8OWXQhJC0NA9P1uOj5ygzPKyfz6sCr9FjZrZv5BUncOIWUovyEBB9FDgz3IakLYErgY+Y2fXFN784irhJNuGJWVbh0hSBM0xIOkzSnWESjrNi9s8OE3n8LkzGcUQd7VzXnmcnNpho3x2+lrSvCTRpoO4JC/KTwV7eHya9uVXSNZLmdO3fUtK4pC9V1+rsrJqjda8q60yiiP6p4xmKeznlksFeJkv6brj/hjD6prNv90giqNskTUmra6hFTp7BeN6Be5UhVr0ESx7vTVkCJ2voX11JDVIWUosyDzghzLJ2ILAyXEV6DPgBwXydyypqcuspQjj0EjBNEjjD4sWRNAk4hyARxy7AcZJ26Sr2T8D3zGwv4Fjgy9W2sn+6BU7TBI/TLjLay++Afc1sd+Ay4N+79n8cuK7stvZDt9joR+j0K446x6XdF1fPnlj3qprOPb9qL04dorMoMtrLycATZrYj8HngU+GxGwPfAk41s12Bg4Hn0uobapGTlyxiolvcFJ2BLU/ZvOUHpahrrUnoxC6kJulUSaeGZa4C7gWWAl8HTgu3vxV4LXBS5Ng9K25/LprizSlT6LjAKY39gaVmdq+ZTQCXEiTliGLAluH/04CHK2yf4zSJnvZiZtea2dPh2+sJogQAkLQPsA3w04raOzBlDq7ziJXusk1Yw6as9M5JwqaFQidL/3IUcGH4/2XA68MH1W8EbjWzWwDM7I9mtjatsqGfk9M2ik413dRFQavOXNdjIbVOGQNOj9n+LYKnB62iiAXYmjA/BzacozNsAqdhxCXgOKCrzD8DP5X0HmBzgpXWG0VSyFXRWc2ckSeLvUQ5GfgxgKSNgM8C76CBNpREmQ9k2nY/LbK90blIQ0wWe1lXxszWSFoJvBB4GWCSribIcHupmXV7RdfDRU4XVa+fk1RHpy1x24eFvEkVnHpoktApiqbNwSmjY5v0rDHtntVJuzMl6kjhOOACM/uspIOAiyXt1pSEHL3mlJQ556QpSQScfJRsLwBIegewL/C6cNNpwFVmNh48qG4e3Q/L+rlXFfHALe4cZU3+T6JTV5n1JIXoxX3udYqiFHspxFYS2Bh4NbAfwfqF10i60cyuSTugMRQ9iG/7OjeDtL8MQRTXnqRMcVWH0jnxFNG5QHOEThEUKXBaHKb2mJntm7AvMQFHhJOBwwDM7Dfh5M8ZwPKiG+o4DWBQe0HSG4CPAK8zs87o8CDgNZJOA7YAxiQ9ZWYbTMauk7oG0lnqLWoNnV5iKbp/kOxukK1fziJgyvheUoR+VtJsBbLZS6fMeDgPZxrwRwKvz3Vm9hiApKuAvYFEkdOYOTl5B7pVrykTpchBea9r6My9KVsI9PNZpqXCruN7GTYmPVvMDayoG2FT5ug0pf4WC5xeLADmSto+TLpxLEFSjijLgNcDSNoZmALUr2AbgofDjRQ97UXSXsDXgCPNbN2DADN7u5nNNrPtgA8SJLdplMApiiLudWnnGOR+3M/cnn6FVTQ1di9aON8mC1n6l3nAieH/RwM/C6cTXA28QtJmofh5Hesv6L4BjRA5ZQ7gyzr3KHofoqIlz1o/Lnjqx4VO8wROUzGzNcAZBB3KEoIsaoskfUzSkWGxDwDvlnQLcAlwUtgJOc5IkdFePk3gqfl+mLime1BXKWkP0Kbds7qIp/kDk6fPatL9uOisZ1UKnSq+94z2ch7wQklLgfcDZ4XHPgF8jkAo3QzcZGZXptXXqHC1PETnc2RZzLNTNut587ShKPKGedVJ9+KoWT7bNl1fU5h2z2pW7jC5kHONcuhaE5MMNNSLA4CZXUWQbTC67ezI/4sJMhY6Xfh8nNEjg730TCpgZhcAFxTdtiR6DWiL7Hs6RO95SX1RlhCtOjwcg87/6b6upnlpqhS2GezlWeAtCcfmSgRVuydn0EFvmWmXR4len0va+jhlrkfkVHvzyUqbPDoucBzHcfJTZt/TSbXc/cp7bBF039t73esHFThNI/o9Fy1s66ZWkVOX4IjOc4m2oQkCqF8RUKZ4KPrcLnTyU1RnU+TNtiihU6bYcYHjOI4zvPQrdvo5ZvKysUL6gqq8OFlC54ZZ4ECN4WpNEBQdmtSWvHQvTFpGCF0ZxJ63vV9DJRQVPlCku7+ozDZlhK+5wHGc4SWa3MFDBMujjLC1MsjTr1V5H67rnp9H3MBwChyoQeS0WVCUTRGfTR6h0z2vJu44X8tmOBlmodPUBAMucBynOFzYON1k7dfiynT6sSYlMRiEzmcR1++MisCBikXORs88B1OrrLH5lCH6sgidvB6gtP0ugKqjiYkIoDlCxwWO4zhOsWTtdzqD57YOmosWOHXf97PU39bvKiutza5WNG3N/JW0YGjS9aQlEOj3+qPHueApn1EQOkAusVP0vJ4mC5wmJqJwnCidUDL3towO0ftS2j2qTYPqLP1jP31oHaF1cd9Jm76Lfqk9u1oTGNaBefcaNb2us4jPoalCUdIsSddKWixpkaQzY8pI0hclLZV0q6S9I/tOlHR3+Dqx+9iqKXKg27RkBB2yCpeivTcucBynGHxR1NGgSfelohYdzbNoZ1kem17n7axplOXVzSgIHHCRAzR3YF40g4SkFVlPTawBPmBmuwAHAqdL2qWrzOHA3PB1CvAVAEkvAD4KHADsD3xU0lZVNbwK2ip0mhqeBi5wHMcZHpLuP8N2X8pz346WjcvyVkYWtTThkoWVO0weGYEDLnLWUdXAvAyvUVrb49JkpyUYGFbM7BEzuyn8fxXBSrvdX8ZRwEUWcD0wXdJLgL8E5pvZ4+GKu/OBwypsfixFr0zdZKHTLWaKTj3ddIFTNZIOk3Rn6NU8K6HMWyOe0e9U3UanuYwteXDdaxToZS+SXivpJklrJB3dtW+2pJ9KWhLa03aVNTwn3QPsvP1PWYPr7lTJRaeVjtue9fxJaZyT2pi2BlB3uue8gqUp4iaDvUyW9N1w/w0dm5C0v6Sbw9ctkv6mV10+J6dCusPGihQWWc4VnXdTtqhJmis0CBs981xau2dIWhh5f66ZnRtXMDSYvYAbunZtC0R75fFwW9L2RtDkOTpAIfN0oLzFQ9sgcKp8WippEnAOcCjBb32BpHlmtjhSZi7wYeBVZvaEpBdV1kDHaRBZ7AVYBpwEfDDmFBcBnzCz+ZK2AJ4vucmFkJahqynZu4rsz6Ln6ucen5TpLOu50tazWbnD5J59RIMEThZ7ORl4wsx2lHQs8CngGOB2YF8zWxM+gL5F0n+a2Zqk+lzkVESe5ABV1V82ZQidFB4zs317FQo7kcuB95nZk+U3qxqaKnSguIQEZeACJ5b9gaVmdi+ApEsJvJzRTujdwDmhZxMzW151I5318cn+tdHTXszs/nDfegImDJne2Mzmh+WeqqjNhRI36K6KaF8VJyQG6c+6Rc2g4WhJQqcXWfqAJKHTFHETIUv/chTwz+H/lwFfkiQzezpSZgrQ88P0cLUKSBvoD2vSgyYiaRMCgfNtM7sipshDQHSm7MxwW9L2RtHU0DUoXkwUQRsETk1k8Vy+DHiZpF9Jul5S7eGbToBP9q+cQTz9LwNWSLpC0u8kfTp80u0USNbwtU5oWXeIWRlza7rJkpGuV3haAwVNHFnsZV2Z0EuzEnghgKQDJC0CbgNOTfPiwBB6csoIBRuUXh6NKjw6Tfo86kCSgPOAJWb2uYRi84AzwicLBwArzewRSVcD/zeSbOCNBKE6Q80we3TaInAGEa56doKxO8aTdmcO70xgY4IEHQcTiP7rJL3CzFb001anONyT0x96diLtsxvUXpLYGHgNQfj0MuC7BGFt5xVw7pFg0DCyoum0IWkeThK9BEqedYqqIKV/KctWADCzG4BdJe0MXCjpx2b2bFL5oRI5eRe4bBJta28LeRVwPHCbpJvDbf8IzAYws68CVwFHAEuBp4F3hvsel/RxYEF43MfM7PHqmp6dIsPWIP2G3Q9Fz9Ppt/4iaaLAyUBaeGcWz+U4cIOZPQfcJ+kuAtGzAKcWXNyUyqD2ksQ4cHMkdOeHBNk/WyVyiu538lKFuOkVChe3ryrS+ooavpteUwey2EunzLikjYFpwB+jBcxsiaSngN2AhSQwNCInac4LtMeL4UKnPMzsl0DqSN3MDDg9Yd/5wPklNK1wyripDYNXxwVOZhYAcyVtT9DZHAu8ravMD4HjgG9KmkEQdnNvlY100omGrbkAKpUs9pJ27HRJW5vZH4BDSBmwOfWSJHTqpIUpvLPYyzzgROA3wNHAz8zMwmMeDBMPzAF2Au5Pq2wo5uRUschlVbSprUXi4q5YyrjxtXmejguc7IQxzmcAVxOkWv+emS2S9DFJR4bFrgb+KGkxcC3wD2b2x/gzOlXTPS/H5+mURxZ7kbSfpHHgLcDXwjkFmNlagoxr10i6jeBB3NfruI5Bqfu+NSirZ09kevBWxfycrGT9zJv03WTsX84DXihpKfB+oJNm+tUEGdVuBn4AnGZmj6XVNxSenF5zXuocQDdZtDTFc9SENgwjbfHoQHnha2UJqWEVOB3M7CqC8M3otrMj/xtB5/P+ipvmZGBsyYPuyamQDPaygCAsJ+7Y+cDupTawAvrta6oKpyoylXQbifuc6+pvMtjLswQPBLqPuxi4OE9dQyFyoJkD5X4ETtUppcsWOtFzx30eTfzehok2CB0oJ3zNBY4zynSEjgscx+lNlv6iaQKn07e30aNTFUMjcqokOlhPGqS3QeBUTef62jZXqu2MmtApMwzOBY7TJlzgOE2nLV6cpgmcKGkLsY46QzEnZ1Ce2W3bnimeo69ex/YSEXGD+zoFTpGip9dnCcG1usCpljbM0YHBBYoLHGdUmNh5ls+1cRrBIKFqbaDoe3/3OjxF0pK1cirDPTk9KMvrkWeQX4XnI2vYWlnhbc/stm0Dl9d0elGWRwfyzdNpo7iB9nTyjuM4RdP0AXkZ9/4mJS4YBXp6ciRNkfRbSbdIWiTpX8Lt20u6QdJSSd+V1LwlzWOI88gM6lnoPrZIEVC0l6Xf/dHPLM1b089n2eTkDHlpsr2UNaAuSwRkFS4ucNpLk+2l6YwtedDD0UaIJttK0+9lUx+wnvfyTpnoy2k/WTw5q4FDzOwpSZsAv5T0Y4KsOp83s0slfRU4GfhK2ome33STgRvcD/0OorMclzSg75XxrZ/6e3lRihALec7RlOxsDaMweymDsjLZFL1oaIdeXh0XOK2n0fbiOA3CbWVAovf0pEU8y6QKL06vviWp/x/WPqmnJ8cCngrfbhK+jGDRqsvC7RcCb8pSYRED/7wD8aL3deg1wB8kKUHea6zLG1K0kGs7RdtLGZR5M6vKqzN52ZgLnCGgDfbiOE2g6bYy7Z7VG7yKPFfR981h9Nb0K3A6+5oePtgPmRIPSJoULr6zHJgP3AOsCBf1ARgHYkerkk6RtFDSwucm/gT0NyjvPibLObKKie4wrCIETrRckeFbeT+DKmhCG5pEUfYy8dyfSmtjW4VO2eIGXOBUTb/2sp6tPP9MZe11nLpoQ98yKGn3yEHEzjDfe7N8LsMoYLKQKfFAuCrvnpKmE6wyulPWCszsXOBcgKnTZq43esiShasIIZOFsgfqRV5LE0VFUlrtLOm2q0DS+cBfA8vNbLeY/VsB5wM7AM8C/8fMbg/3/T3wLoKnZrcB7wwXq4qlKHuZtsmLSn3MVOYibGUkJCibsp/qtamTlXQY8AVgEvANM/tkQrk3EzxF3s/MFvZTV7/2UqWtOE4avexF0mTgImAf4I/AMWZ2fxh29g1gb4Lx2EVm9m9J9RTVt2w5dWap9lJkv1LUuZowyI+GdRfR3+TpU7L091X1Uf3aS7jvwwQhmWuB95rZ1Wl15UohbWYrgGuBg4DpkjoiaSYl5MZq4mB+EDqD/G4R0BSPTFE09JouAA5L2f+PwM1mtjtwAoEBImlb4L3AvqE4mgQcm6XCIuxl7I7xLMUaSZtCAVzg/BlJk4BzgMOBXYDjJO0SU24qcCZwQxH1Vt2/OE4RZLSXk4EnzGxH4PPAp8LtbwEmm9krCAZ0fytpu151NtlWBhUTneOHNXwKiulviuxTyggHTGIQewnLHQvsSjCe+3J4vkSyZFfbOnxqgKRNgUOBJQQGdnRY7ETgRxmub2hIWm8m6dUhyctRNk1Ym6ZO0WNm1wGPpxTZBfhZWPYOYDtJ24T7NgY2DTuSzYCHk07SNnsp+8bWdKFTRVx2mwROyP7AUjO718wmgEuBo2LKfZyg80n0avaibfbSja9V45DNXo4imC8Dgefz9ZJEEB2wedi3bApMAE/GVdJ0WylSlNQlbjrr11QdhVCVyKhirlMGBrGXo4BLzWy1md0HLA3Pl0gWT85LgGsl3QosAOab2X8BHwLeL2kp8ELgvEyXN6TkDUWrWuA0kYZ5em4B/j8ASfsDc4CZZvYQ8BlgGfAIsNLMfppynsLtpWxvzqgKnSra1UKBA0FMfzQ38QZx/pL2BmaZ2ZUD1tXa/iUqbrIKnY4oGkVhNMTX3dNeomXC+TMrCX7XlwF/IuhblgGfMbOkh3GNtZVh8Lh0C5sqBE+3yMgqOlrar3QYxF6yHLseMqtuACJpFXBnZRUOzgzgsbobkZE2tRXi2zvHzLZOOkDST8Lj4pjC+k+Uzw1jkKPHbwf8V8KcnC0JQtT2Iph3sxPwbuAB4HLgGGAF8H3gMjP7VlI7i6Jl9jIMv78m093eVFuB/u1F0tHAYWb2rvD98cABZnZG+H4jAq/nSeG8gp8DH+x3Tk4RtMxWoP2/v6bTGHsJt90elhkP398DHAC8HDgNOAnYCvgFcLiZ3ZvxOvvC7aV02t7eQewldSw2oL38M3B9Z/wl6Tzgx2bWyS64AZkSDxTInWa2b8V19o2khW1pb5vaCv2118zS5tQMhJk9CbwTIHSL3gfcC/wlcJ+Z/SHcdwXwSqB0kUOL7GUUfn91UrG9PAREH7l3x/lPBXYDfh6YCi8G5kk6skah0xpbgdH4/dVJw+wlWmY8DE2bRjCh+m3AT8zsOWC5pF8B+xL0PWXi9lIio9Demuwly7HrkSvxgOMMK5Km688rRb8LuC4UPsuAAyVtFoqf1xPEQTvOsLIAmKtgJfUxgome8zo7zWylmc0ws+3MbDvgeqBOgeM4dZJqLyHzCObLQDB/5mcWhNEsI1jnBkmbAwcCd1TSaseph0HsZR5wrKTJkrYH5gK/Tausak+O49SCpEuAg4EZksaBjxIspoaZfRXYGbhQkgGLCLJ7YGY3SLoMuAlYA/yOMA2n4wwjZrZG0hnA1QTZBM83s0WSPgYsNLPuDslxRpaM9nIecHE4b+Zx/pyh8xzgm5IWAQK+aWa3Vn8VjlMNg9hLWO57wGKC8djpYVr1RKqek3NK9zyJJtOm9raprdC+9tZBmz6jNrUVvL3DRts+H29vubStvVXTts/H21subWtvHioVOY7jOI7jOI7jOGXjc3Icx3Ecx3EcxxkqKhE5kg6TdKekpZLOqqLOPEg6X9LyMG1dZ9sLJM2XdHf4d6s62xhF0ixJ10paLGmRpDPD7Y1ss6Qpkn4r6Zawvf8Sbt9e0g3h7+K7kYn/I43bS7G4vQw3bi/F4bYy/Li9FIfbS/MpXeRImkQwue5wglXlj5O0S9n15uQCoDsd3lnANWY2F7gmfN8U1gAfMLNdCLKxnB5+pk1t82rgEDPbA9gTOEzSgQSrpX/ezHYEniCc7D/KuL2UgtvLkOL2UjhuK0OM20vhuL00nCo8OfsDS83sXjObAC4Fjqqg3syY2XUEGRyiHAVcGP5/IfCmKtuUhpk9YmY3hf+vIkhpvC0NbbMFPBW+3SR8GUHqzM4iTo1pb824vRSM28tQ4/ZSIG4rQ4/bS4G4vTSfKkTOtsCDkffj4bams42ZPRL+/yiwTZ2NSULSdsBewA00uM2SJkm6GVgOzAfuAVaY2ZqwSFt+F2Xj9lIibi9Dh9tLSbitDCVuLyXh9tJMPPFABsJFiBqXhk7SFsDlwPvChSvX0bQ2m9laM9uTYIXa/YGd6m2RUxZN++11cHtxmkjTfnvgtuI0l6b9/sDtpclUIXIeAmZF3s8MtzWd30t6CUD4d3nN7VkPSZsQGNW3zeyKcHOj2wxgZiuAa4GDgOmSOgvStuV3UTZuLyXg9jK0uL0UjNvKUOP2UjBuL82mCpGzAJgbZm8YI1i5tA0rZs8DTgz/PxH4UY1tWQ9JIlgRdomZfS6yq5FtlrS1pOnh/5sChxLErl4LHB0Wa0x7a8btpWDcXoYat5cCcVsZetxeCsTtpQWYWekv4AjgLoLYv49UUWfO9l0CPAI8RxCPeDLwQoKsGHcD/w28oO52Rtr7agL3563AzeHriKa2Gdgd+F3Y3tuBs8Pt/wv4LbAU+D4wue62NuHl9lJ4e91ehvjl9lJoW91Whvzl9lJoW91eGv5SeIGO4ziO4ziO4zhDgScecBzHcRzHcRxnqHCR4ziO4ziO4zjOUOEix3Ecx3Ecx3GcocJFjuM4juM4juM4Q4WLHMdxHMdxHMdxhgoXOY7jOI7jOI7jDBUuchzHcRzHcRzHGSpc5DiO4ziO4ziOM1T8/3LH2Tj9XQE5AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1008x180 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots(1,4,figsize=(14,2.5))\n",
    "\n",
    "counter = 0\n",
    "for c in [c1,c2,c3,c4]:\n",
    "    s = ax[counter].contourf(c[0,:,:])\n",
    "    ax[counter].invert_yaxis()\n",
    "    colorbar(s, ax=ax[counter])\n",
    "\n",
    "    default_peak_position = np.floor(np.array(c[0,:,:].shape)/2)\n",
    "    i = np.array(find_subpixel_peak_position(c[0,:,:]))\n",
    "    ax[counter].plot(i[1],i[0],'rx')\n",
    "    print(np.array(i - default_peak_position),\n",
    "          np.sum(np.abs(np.array(i - default_peak_position)-np.array(shift))))\n",
    "    counter += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "709db5c5",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAChCAYAAABKxo+2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA02klEQVR4nO2deZgV1bX23326oRmbOSioOEXFAUecjYjiQIwSruKAnxo1osZoNA5R891o4sX74b1XjF4ThxtnjXFK1IQgRpwV0EscwSkOiIAyNzP0qe+PKvSs364+XRi7jU/W73l49O2qU7Vr166qc+rda62QJIkcx3Ecx6lO6atugOM4juN8HfAHpuM4juMUwB+YjuM4jlMAf2A6juM4TgH8gek4juM4BfAHpuM4juMU4B/2gRlC6B1CeCqE0BBC+M+vuj1fhBDCr0MI//erbofz9xFSbg4hLAghTM7+dnoIYU4IYUkIoUcL7tvHkPO1I4TwRAjhlC/42Y2y66rmy27X30uLPjBDCO+HEJZnBz87hHBLCKFTwY+fKmmupPokSX7cgs1sMZIkOS1Jkl981e1wqhNCODGE8GoIYVk2Tn8VQuhascrekoZI2iBJkl1DCG0k/ZekA5Mk6ZQkybwvuN9BIYSPqq3jY+jrRQjh2BDCi9k9b1YIYVwIYe+vul3/yGTPiQPW6iRJPsyuq8avsl15tMYvzO8kSdJJ0g6SdpR0UcHP9ZP0RvIFMiuEEGrX9TPOPychhB9L+n+SzpfURdLuSsfehBBC22y1fpLeT5JkaaZ7S2on6fVWbq7zD0wI4VxJYyWNVjpGNpJ0naTDc9b92tyj+Evv69T2L50kSVrsn6T3JR1QocdI+mOF3l3Sc5IWSnpZ0qDs77dIWi1plaQlkg5Q+nD/iaR3Jc2T9DtJ3bP1N5aUSDpZ0oeSnsr+fpKkaZIWSBovqV/FvhNJp0l6O9v/f0sKFcu/n322QdIbknbK/t5H0v2SPpX0nqSzqhz/LZIuz/5/kKSPJF0g6RNJsyQNkzRU0luS5ku6uOKzu0p6PmvbLEnXSmpbsfxASW9KWqT0onxS0ikVy5s8dv/3WR/VZ+NrBP7eKTu/J2VjaoWkxmzduyUtzcbPEkmPSwqSrsrO62JJr0raNttWnaT/yMblHEm/ltReUkdJyyWVs+0skdTHx9DX85/SL1tLJB3ZxPJLJd0n6Y5sjJyS3Useys7bO5K+j3P3YrbuHEn/lf29XbaNedl5nSKpd5V2NXUf6y/piWwbr0s6DGPuV5L+lI31A5Teyy+U9IqklZJq1cT9O9vGE2vHkqTNsutkntK3hndK6potuz27BpZn/XeBPr+f12brVOunS5U+C27LjvF1Sbu02Hlu4UH0vrIHpqQNlN5Irs5036wDhyp9GA7JdK+Kk3Z5xbbOlvRCtp06SddLujtbtraDb1N6I2qv9FvdO9nAqJX0U0nPVWwvkfSIpK5Kvwl+KungbNmRkmZKGqj0Zri50l8ZJUkvSfpXSW0lbSrpb5IOauL4PzsGpTe7Ndln2ygdyJ9KuktSZ0nbZINmk2z9nbMBWZsd3zRJP8qW9VR6IQ3Plp+t9AvG2gFa9dj932fn5+DsnNTmLLu1YnydKOmZimVrx9vaC/qgbFx0zcZLf0nrZ8uuUnqxd8/O88OSrqgYEx8100YfQ1+Df9XGUrb80qx/hym9j7SX9JTSLyrtlL6B+1TS4Gz95yX9n+z/O0naPfv/UdkY6iCpJjvH9U3ss6n7WJvs3F6s9D42WOnDZsuKMbdI0l5ZW9spvZf/VdKGWdubu38/UTGWNs+W10nqlR332Ip2vi/7w2pj2eurWj9dqvQL7dCsP66Q9EKLnecWHkTvK/3W0JB1wF/0+TeLCyXdjvXHSzqh4qRVPjCnSdq/Qq+fDcC1N4NE0qYVy8dJOrlClyQtU/YtOVt/74rlv5P0k4p2nJ1zPLtJ+hB/u0jSzU0c/2fHoPRmt1xSTaY7Z23YrWL9lyQNa2JbP5L0YPb/x0t6vmJZkDSjYoBWPXb/91m/HCdpdhPL/l3ShOz/T1T1B+Zgpb/wdpdUwnlZKmmzir/tIem9ijGxrg9MH0P/gP8kjWxqLGXLL1X25ivTGyp9a9G54m9XSLol+/+nJF0mqSe2c5LSX3UDCrSpqfvYPpJmY6zeLenSijF3Gz7zvqSTKnRz9+8nVPG2AusNkzQV2859YBbop0slPVaxbGtJy1vqPLeGhzksSZLOSi/2rZR+s5XSbzpHhhAWrv2ndHLF+k1sp5+kByvWnaa0I3tXrDMD619dsf58pTeFvhXrzK74/2VKv8lJ6Ul6t4k29EGbL0YbqjEv+dzIXp79d07F8uVr2xBC2CKE8Eg2CWWxUl9kbd/1qTzWJB0plZNHihy7k74e6tmEJ7N+trxZkiR5XOnrzv+W9EkI4YYQQr3Sb9MdJL1UcS7+nP39i+Jj6B+TeWp6LK2l8v7UR9L8JEkaKv72gT7v35MlbSFpeghhSgjh0Ozvtyt9MP02hPBxCGFMCKFNCGGfbKLRkhDCWm+9qftYH0kzkiQpN7FvtjXvb4Xv31nEw29DCDOzcXiHPh+HzdFcP0nxfbxdS/msrRZWkiTJk0q/ufxH9qcZSr+hdK341zFJkn9vYhMzJB2C9dslSTKzcjdYfxTWb58kyXMFmjtD6Xv3vL+/h212TpJkaIFtriu/kjRd0jeTJKlX+mAO2bJZSl9NS0rDHiq1/r5j/2fieaV+zPDKP2YzuQ9R+kakEEmS/DJJkp2VfsPdQukkorlKH2DbVJyHLkk6CU6y47Ul8DHUeqwdS8OqrFN5vj+W1D2E0LnibxspfYWqJEneTpLkGEnfUDop7b4QQsckSVYnSXJZkiRbS9pT0qGSjk+S5OkknVnaKUmSbbLtNXUf+1jShiGEyvv/Z/vOaWve39bl/j06++x22Tg8Tp+Pw6b2VdnWJvuptWntOMyxkoaEELZX+i3jOyGEg0IINSGEdtk0+w2a+OyvJf1bCKGfJIUQeoUQDq+yr19LuiiEsE22fpcQwpEF23mTpPNCCDtnMXibZ/udLKkhhHBhCKF91u5tQwgDC253Xeis1GNaEkLYStLpFcv+KGm7EMKw7JvUDyStV7H87zn2fxqSJFmk9LXXNSGEg7Nv6hsrfT3/kdJv880SQhgYQtgtCzdZqtRTKWff4G+UdFUI4RvZun1DCAdlH50jqUcIocuXemCf42OolcjG0r9K+u+sTztk4+mQEMKYnPVnKH21ekV27xug9FflHZIUQjguhNArG0MLs4+VQwj7hRC2y2auLlZqS5W5/Yym7mOTlP4SuyBr4yBJ35H023U45HW5f3dWas0tCiH0VfplspI5SueDRDTXT61Nqz4wkyT5VOnEnH/NOuJwpd96P1X6jeX8Km26WunkiUdDCA1KJwDtVmVfDyr9Zvbb7DXAa0p/NRRp572S/k3pZIoGSb9XOiO3Uek3uh2UzpCdq3RQtsQN7zxJx2b7v1HSPRXtm6vU0B+j9FXQ1kpn1K3Mln/hY/9nI0mSMUrH4H8ovQFNUjoW90+SZGXBzdQrPUcLlL4umifpymzZhUonWLyQnYvHJG2Z7Xu6Uu/ob9lrrT5fykF9jo+hViRJkv+UdK7SCVJr72lnKr1/5HGMUr/uY0kPSvpZkiSPZcsOlvR6CGGJ0nvf0UmSLFf6peY+pWN1mtKZzblf7Krcx1YpfUAeovQedp3SX6nT1+FY1+X+fZmknZROJPqjpAew/ApJP82ugfNyPl+tn1qVkBmlzteY7NXKR5JGJkky8atuj/P1w8eQ4zTPP2xqPKc62auQriGEOn3uTb3wFTfL+RrhY8hx1g1/YH592UPpDLi5Sl+vDMte2ThOUXwMOc464K9kHcdxHKcA/gvTcRzHcQrgD0zHcRzHKUDVbAijR48272vffvtts3zRokVGd+lioytWroxn5dfV1Rl98803G33bbbcZPXXqVKNnzrTxqt/97neNPuaYY4w+88wzozZss802Ri9dutToCRMmGL3vvvsavcEGNtTokUceMfob3/iG0W3atInawH3yMwsXLjR61apVRn/88cdGd+jQwei99trL6Pfffz9qQ9++NmnLp59+WnWfV111VVArMGTIEDPuNt98c7N8q622Mvrkk082evLkydE2r7jiCqM328zGc/fpYyM6Hn30UaPr6+uN7tXLJuuZM2eO0d26dYvawHPKNnz0ka30xetr1113Nfrll182mtfbihUrojbU1tpLnmOzbdu2Ru+www5GX3rppUa/+OKLRt91111GH354HCrN4/zDH/5gNK+v6667rsXH3YG7XlbVmyrX4VZZtquXVqyJPlOzyF7jSckexsp+3Y1uO3eZ3UCw66/q3r5aE9VmcXy/5T7L7e1xhNU2hLMErbLVSRtbojKsthW4Skty2oC+K6FfGnvaa6vc1q6f1NhjCOj72gW231b3jCtIhjX2OGoabDsbO9vn0mPP/jR3zPkvTMdxHMcpgD8wHcdxHKcAVV/JPvnkk0bz9c3559sMR9dff73R3bvbVw6StPHGGxs9YsQIo7t27Wr06tWrjeYrogcffNBovo7j/qT49STX4atlvnb605/+ZPTQoTaVLF+X8nWdFL8ue+ONN4weNmyY0a+/bmsVsx8IX9/x85L0v//7v0b/8pe/NPrOO++suo+WYqeddjL6ueds+tINN9zQ6NGjRxvN8ydJAwfa7IWvvfaa0bQKGhoajO7YsWOVFsev1H/zm99E6xxyiE2Uw9ejS5YsqbrPd9+1ebT56nKXXXYx+r777ovawNe2AwYMMPovf7Hpc7mPo446ymjeE/ia/29/+1vUhldffdVo9kO53FSmt9YjqbW/JVZ3sm2sWclXmfbVpCSt3Mje/2rw2nZNB/t6s7GPfZVYs8ruo80ie+5Wd7FjdlVXqyWpzWJrqyR4zdvYAY8AvMKtbbD331Xd7Cv7dnPs69CkXWw/ldvbvyV45VrGa95GvDYuNdp+aGxr1y+1s20Kjc2Pn9U9rYVVs2x1E2ta/Bem4ziO4xTAH5iO4ziOUwB/YDqO4zhOAap6mFtvvXXVD0+ZMsVo+nIMS5DicAr6F5yO/8knnxjNqfb0jrj9fv36RW1gaAo190HYL/Q4GYYyadKkaBtDhgwx+vTTTzf6nHPOMZqe5e9//3ujv/3tbxvNcJ08eL4uvvhio/fbb79mt9ES0IdmGNC4ceOM/vDDD43OG7dbbLGF0fTGH3roIaM5buhLk/bt7ZT/H/7wh9E67O9p06YZvWaN9bgY4sHriWPivffeM3qTTTaJ2rBgwQKjGY5Ej5Led02N9Y+OOOKIqp8fOXJk1AaG1/A4p08vXDTjS4N+YM1yey7q5tlzt7rerr+iVxzyUbfAeo5rOlovr2a59T3Ldfb3S7mN1St7tLM7gB+5ulP8+6exnT1fdfNtm2qXsA32frymM9q8onqYyequaKOkNe3tOkLARt08tAl+In3XpLPd3sredgzTt5Vi7zWqwFkwcMl/YTqO4zhOAfyB6TiO4zgF8Aem4ziO4xSgqodJr4G+2rXXXms0/cM8GGuWF6tZCb0jenOXXXaZ0UwJ9+yzz0bbZIwi06Ixtoxxm9/73veMvuaaa4w+++yzjabHKcUxjkzhN3/+fKN5LuhfNRdbyrSGUuzT8fyxjUw72FIwjpIp4JiWjqnW6IlJ0ltvvWU0432ZXo/jZp999jGa8YX0TZ9++umoDUzxR6/8m9/8ptH0D7mc8Yo8f3kpGfm3o48+2ujly211L8ZMMjaU3u+bb75pNMdxXjs5VrfddtvoMy3NqnqkDISnFdohNnCV9f4YUylJn+xs4yo7fGI/02YJ09BZ2VhnG1GDrHO1S6zXF9bEbVjZHakPF9nfSGs6VfcXE+g18A/L6KjQGGcYZLq9xna2DYy7rFlh+4k+LNP9Me3dqi7x9d8If7h2uf0M0+01hf/CdBzHcZwC+APTcRzHcQrgD0zHcRzHKUBVD5N5J+lpMf8pvUCWi5Li/LP0KP/6178a3aNHD6OPPfbYphus2OujlqRly2z+Q5ZyYrt/9KMfGT1x4kSjmduU3tQtt9wStYE+HNs5atQooxmfypg7xsTyGPP8LMafXnfddUafccYZ0WdagySxfgLLec2YMcNo5jvNy1/6zDPPGE2fdPjw4UYzBpVjZMcddzSa8cK77bZb1Aa2e9asWUbTo9xjjz2M5jl/6aWXjKb/yDEgxbljOQ+B/i/HNvuW1wqPMe8eMGjQIKPZd19FHGbdQusHruwaXy+VJCXEK7aLA/mmXmKvp22vttdTr5etj7ayq/XqIo8T0Hcrt41//7RdbP1A5sTlNpgjty1iSRs72H7hPktrYi+wDM+x7SIb48ptNNbZfqhtsHGVZeSSXdOx+jFJ8XHRmy2K/8J0HMdxnAL4A9NxHMdxCuAPTMdxHMcpgD8wHcdxHKcAVSf9MKn1vHnzjGbwPJMGMMG1JN10001GH3nkkVUbOHjwYKOZmJuTL9iG3/3ud9E2GTjNBO5MbMCkACwgTTixiUV1pXjCDdtw5ZVXGs3JGJy0wolGJ5xwgtETJkyI2sDPsG+bS0LfUrzzzjtGM2kAJ0xxkkjeuOM5r6+vN5oFox9++GGjmTidScyfeOIJozt37hy1gUH/vH6YCP3CCy80mhN0WACcx3TiiSdGbbj11luNZuFsJhE466yzjOb1e8899xjdqZMN1uekISmenESYVKI1aMTEE04cKa1Kqi9fE//22O0ntqBCtwYkLkBxZgbsM4kAJyLVrKyBjicJMSED1wmYH8MkAZyQU7PMTtgJmNCzqkv8SGHfqRbJ1HGcjR1RULqdTeheu8z2U1TMe1VczHvpBjY5fhsknVdOwoU8/Bem4ziO4xTAH5iO4ziOUwB/YDqO4zhOAap6mPQXmVSA0CdikWMp9uKY5Jr+xeOPP240PReuf8MNNxhN/1GKEyyQnj17Gn3vvfcaPXDgQKPnzJlTdXt5yROYaOCCCy4wmh4ng9qZ8IH+GPuey6U4KcQdd9xRdZ+tBT1KJiE//PDDjabnmZekgUW9mbh8l112MZrB9IRJA5gwfPbs2dFnjjvuOKMfe+wxo6dOnWp0Y6P1WXi+6JWPGDHC6KuuuipqA88xEw0w2QELYfMaZ0INXgt5Bdx32mkno+m1sph6a0Avr2YlkgLAd6OHWbs09s3q5tmAewbYUzNJQFJj92nvfDEsOC3lJCJYiDZ1wCMAniSLMbOYM4s9N7aP28DjEKYY8LjLSMbeiONisvWalSiCXVskgQPXqZ6oYi3+C9NxHMdxCuAPTMdxHMcpgD8wHcdxHKcAVT1MekeMJ6T/yGKxeYV86Q/yM9znyJEjjaaPSp9t++23N/rGG2+M2kB/kLFnTCrPBNKMT6Snye0PGTIkasMHH3xgNPuSbWIR6lWrrLfAGD7CRPqStGaNjamiz9mcj9dS7L333ka/8sorRnPM0Hdrh7gtSdp///2NpgfJcUX/njHITITOMbHRRhtFbTjppJOMZrsZ+9m/f3+j6aX/+c9/NnrTTTc1eujQoVEb6IPSL3z33XeNZvzu+PHjjT7qqKOifVSSdy7Yrttvv91oxry2Bm3n2djQNZ2tYxgQj7gGXl27T+Lri3GUkUcJr21FT9tXdVifBaOFIgVK4vttCf5e5IvOtce9upttQyMKZ7N4M/3DdrPjhP8qwYNcDb8X8wlWoh/ou7JN9IJrVsXxqG0Xoi+XI4F7u7j4dh7+C9NxHMdxCuAPTMdxHMcpgD8wHcdxHKcAVT1MepBvvfWW0euvv77RzBuZFw/HnJv0zaiZt5XxcvSe6AUOGzYsagNjFLkNens8zt69e0fbrGTs2LFG58WV0WttDuaOZRwmfVbmos0rYv32228bvc022xj93HPPrVMbvyxYpHjzzTc3mnGto0ePNpo5V6XYL2Ss55IlS4x+7bXXjC7Bh+G1QQ+ZXqAkLV682GgW6H7ooYeMprfHnLps080332x0Xrwxj5u5lukfMt/tEUccYXRNjfV+vvvd7xp9wAEHRG245JJLjGZM7HbbbRd9pqVZBe+u7SLrm5U+tXplb5szt3ZRnB83LLN9We5qPxMXVrb7YLwhPctynb19t5kf+4dJnb0HN7ZFseVV9l5Hn5TFnGuQpzWKX51rx7gklTvbe3pAfDHbWDfX9ltorJ4rll5wuU3zfmTNEsaXNhflmu270FqO4ziO80+OPzAdx3EcpwD+wHQcx3GcAlT1MOmTnXLKKUbTF6Nnklc3cuutt666j6iBtbaJ9Bu5nB7MH//4x2ibzIdJ746xoPQ0maPzzDPPNJqeWl5c2csvv2w049X4GcavMp6O/hi93Dw/mX4x8+5+Vblk6QeyPukf/vAHox955BGjN9lkk2ibTz/9tNHMccocxYyTZQwyxzpzsp5zzjlRG44++mijmVuWXi3HCGumHn/88UbTm8+rNcuYVnqthx12mNH0sffYYw+jeT0yFpTXgiQdeOCBRjOe9IEHHjD6sssui7bxZcOcqYz1Cyutt8ccquX28fVVove2yN53wkrETSJPa8L4xRU2lrDU1u4zLI191KSNPY7SYvh78I/DSrvPWhx3uYNtc1iO2NByXFeyhHaxTWxDzRJ7v20OHne0/bzP4LhKDe5hOo7jOM6Xhj8wHcdxHKcA/sB0HMdxnAKEhPkIKxgxYkTTCxXXsWOtvYMPPjj6zOTJk41mzlT6aPTe6D8yTytzrrKephTHitFH5T4Y49hc7ljG6O21115RGyZOnGg0YwsPPfRQo9lvzbXxhBNOMPrJJ5+M2sAYVcbk0ccbP348smO2DIMHDzbjjjVP77vvPqMvv/xyo0eNGhVtc9CgQUb/4Ac/MJrjhPmEt9hiC6MZX8gxwPqbknTTTTcZzdhQXovt27c3mj4142jpy9JPlKTVq+HFYQ4APU7mIKaf37lzZ6MZt8kxJMV5kJmnl8c5c+bMFh93B295oel8xlAmHe25YExkaLDXoySpzt4HEnqOy5F/tga/X+BhJrXwH+HVhcU5beA26TFiedIeOXTh9SXt4Luuht+Y9zyBN0sd9TX3QdjX9Ta+Nayw4ysX+KZJvZ3P8efXR+eOOf+F6TiO4zgF8Aem4ziO4xTAH5iO4ziOU4CqHuZ+++1nFjKnavfu3Y1mXB+9Cin2NBjPxhqMjD2jb8rYM8Yf5tXk5D4Zb0jPkXlZmY+TMK6TOXglqVu3bkbTq2U+Vfo69H7pwzKuc/DgwVEb2Hf0UelfTZkypVU8zOHDh5txx/qXrCvJOEzGRErx2HznnXeM5vnYbLPNjKZXV19fb/Sxxx5r9NVXXx21gbG13MfUqVONZk3Nbbfd1uhbb73VaHqYDQ0NURs41nl9MK6SuWOfeuopo5kHltcGc/JK8X2E/cJ2P/300y0+7g7pc6YZc0knGwsc6NWRVavjv9VU9xyTBpu/OLB2KHKu0tPk9iI/UVJ5se3LUg87zqNtcp+MDWWsZ3eb4zrk9cOa6vUvo312gF9Mv5HeLfbJ3LWSVFqC+QLN+MHuYTqO4zjO34E/MB3HcRynAP7AdBzHcZwCNJ90r4Lm4qX23HNPoy+++OJoG4899pjRzPV66qmnGs1YszFjxhhNj5OeDHPXSnFMHetVMgcuY9cYF0hNb5BtzNsm4yiZT5XQ86RfzDhNxrfmtSEv3+xXAcdVp042zoq5RidNmmR0//79o23yHPBYA7wajhF6cdwe853Sf5Rin5l+/ezZs42mh8xzSN91yy23rLo9SXr11VeN3nHHHY1esGCB0axfO2DAAKMnTJhgNHMes86oFM8RYF8yV3BrwLkczE9aXrDQ6FIve5xJV+txS4o8RW4z8iwxBhvn23NRQr7jAK+vvJ7NLy1J6ma9di1Avcq2zVzzyF+r1ajZSU+TMZeStDrH16ykHfK40tNE7GjSBevDhy0tisdc5JsusffHZGWB2E35L0zHcRzHKYQ/MB3HcRynAP7AdBzHcZwC+APTcRzHcQpQddJPcxNoOCGAkxpGjBgR7xDJnnv37m00JzY0V8yZkzeY/JlJy6U46J8TNjiJh9vkRCT2EydOMFGCFE/oYEA5i1Tvt99+RjMQn5N8GECelzyB+2TCBSZjby0YDM8ix/vvv7/RnCTCQsxSnPyAk6wY9D9lyhSjmdSBxZ05ySevv7fbbjuj33vvPaPXW289o5nQnWO/BoHxvFbmzJkTtYGTcjjpjhNyONGPCeOZhOCaa64xeuTIkVEbLrroIqNZaKCRkz5aAU6gSRDwH7pg8gwnCS2Li8RHk1UwgSbBRL1SvZ04VLOenQyVcKJKF7t+UhNPuCmh0HWUBGCevVeFjnayYXmpPd8RmCwTOuVMVmSSgOXoKyY2aCZhQ/h4rtFJHzumWWhbkpIV6AdcO9EErCbwX5iO4ziOUwB/YDqO4zhOAfyB6TiO4zgFqOphduliE+vS2/v2t79tNBOG03OR4uTpQ4cONZq+Z942KmGyaPqJ9OUk6d577zWaXg+9Jnp99AsZSL/BBhsYTc9Tkvr27Ws0j/Occ84xmkkGeNzNkdeP7Bv6xzzu1oL+Is8pE+wzycC0adOibfJYmBiC/h89y+XLbYA2PeJNNtnE6DwP8+677zb6kEMOMZrjir4rkwDweuS4zSvezKQQ3AaLvt9///1GcxzuvvvuRp922mlG5yXt4HHPmjXLaPrzrUGC5NzqZu99mmvPRXmR7cfI45RiDxNjqIyk86XuSIxOb49JAeD1lRpyfFR+pr316gKSBrAfomTtICqszaLYUuT3JkzYjmQK9HajpPJMvj7TJgBRp3j8BBbKxv0w1BbL4eO/MB3HcRynAP7AdBzHcZwC+APTcRzHcQpQ9cUtvQT6j4y7ZAwWY76kODH5G2+8YfSiRYuMZlJxxlBuv/32RjN+jt6UFPs2jHdjrCj9QsZVMuk825Dno9IvpA/HGFfGwNJjo0fHNjOWVIr7ljqv+HZrcPbZZxvNc/joo48azWPPK/DNc0yfmQW3mfCdmvGH48aNM3qrrbaK2rDrrrsazXPEccd98BhYwJ0+al4CePrtLIRNz/KFF14wmjGSjONkoe28BPAsnP3mm28afdRRR0WfaXHoaa1E3F57eHU4V6FA4YJQ29VqxDgmiHENLMDQ3X6eMZV5xZsbu9nzUVqK+FH6pPAXGWfJwsthsfXd80jo99IXZfJ19EuAn0zfNPnbh/bzOfM1Aq5fFu+O2tAE/gvTcRzHcQrgD0zHcRzHKYA/MB3HcRynAIGFUyu57LLLzELGxzVXcDgv/pCxZozr4nLmtmTRXfo+KxDbxFhSKfZe6d0x3o2eJGHsKL0qFtmVpIEDBxrN42L8GuP6eC7YT+zXGTNmRG2gL0dflLGbDzzwQE512C+f9dZbz4y7DTfc0CynB/3SSy8ZTU9Zij3KF1980Wh6c9wH+4r5ahnfuO+++0Zt+PBD67UwB25z+/zoo4+MZl7XPfbYw+i5c23OTSkeVxwXnHfA64lFq4877jij77rrLqN32mmnqA033nij0Tw39D0XLFjQ4uPukL4/NGOOBaHL7ayfX1qCYtB58YcozpyUcBjzFtptMCYSeVmTjrafSosRO8oiyVIUy5m0RxFqthsxjlHeV+TcjfaZk8eVsZ9Rnl7OlcD6jMNM6lD8HbGjSWPcD1Hf4nxyG+NmXpM75vwXpuM4juMUwB+YjuM4jlMAf2A6juM4TgGqxmEy/o0eC+PfGD83evToaJv08xjryX3Q06Svxs8zj+igQYOiNnCbeV7rukA/kbFp9BOl2LOcOHGi0c35WYztZOwhzwVrfkpx/Bv9rSFDhkSfaQ1Y55GaeWAZw8p4RUkaMGCA0YwXbK42JffBfLb00idNmhS14ZhjjjF6xx13rLoNnh/6ifvss4/RHGeMaZZif725mGRy2223GX3EEUcYTR987Nix0TZYS5TzDJ5//vmqbWgJmstfWkJeWPqTeTlUI68OOU7L/WycrZZV9xPDcsREIj+q6JEqjrMMDaipyfqVyD0b6qzXlzBGkjlYca3mtTOqPYnjjLxeepaMkSU5/RCdX2wzyonb1KYLreU4juM4/+T4A9NxHMdxCuAPTMdxHMcpQFUPk3Fe1PQr6AUy3lGKPUfWoqQvSs/lyCOPNPqMM84wmp5MXt3IYcOGGc1200tim+gn9u/f32jW91u5MvY3pk6dajRrGdLjHDlypNHMyctj4LlhnlFJOvDAA41mX1577bXRZ1oD5rBtaLBeEL0/1mFlXJ8UjzPCc0iP87XXXjOasbocM3k1HbnNV155xWh6lBx3zC3bnMfJfLlSHP/L2FD6oIxPZT5pjjP64Oeff37Uhnnz5hnNGON+/fpFn2ltkln2+gvr2+MKjA1Ercv8jcIHXWPjBZN29OqwD3qaiHkM9FkVx10ms228cIk5VFm7Ep6lGOOIWPwy4lclqbTMjsukA+Iy4ZOWFiHPK/LXRrGikdcb1wWNPEzUN1XvOP90Hv4L03Ecx3EK4A9Mx3EcxymAPzAdx3EcpwBVPUzGu9FHY45V1rZk7T1JOuigg4yuq6urqg8//HCjb7nlFqPpHTGnKtskxXGT9KN69+5t9KGHHmr0M888YzRjP+mzsvalFPu7CxYsMJr1L0855ZSq26RXy7jBe++9N2oD+47+E8/V+PHjo220BPShCY/9gAMOMJrnR4rP+aabbmo041rLyJFJP5Ae5fTp041mflsp9iwZX5qX97iSd99912j6kYz9zPPv6XMy1vq5554z+vrrrzea9TEZF81cs8x/K8VzABgvmtfulibpa8d+ab71zYXYv6QT4va6xecurIBvxhhHeJBretnapLXwMOkFslZlHoH77GlrqEYxkmwza1difeZtDXltYr5Zeo7wgxlHmSDPq7rZfiojx25NTn70QG8W3nxe7GYe/gvTcRzHcQrgD0zHcRzHKYA/MB3HcRynAFXrYZ566qlmIb2+5uLG8nK00jdjnCVzvzIfLX2ejTfe2Gj6W3l+GOP8dtllF6Nvuukmo+knfvDBB0azpmBezCOh90MYh0lfjrUOGVvap08fo9nPUhxD9+CDDxpN3278+PGtUg+zXLYBZeeee65ZznNeKtnvfQ8//HC0Tfp9jLtkbUp6cfTz6Tn36NHD6LzYW14vjEckzOVLD7N7d+tHcZ95uWRZ75LXC2Ne6bOyrxm3yTkGjz76aNSGxx9/3Gj2A2uJnnbaaS1fD3Oz82w9zAZ7rwudre9aRu7RQE9MkhgXSW8OfUtPUwttzCTbwFqXUZxnzt/KvbpiG/AgcRxhoR0PCfzDaJ85saCRhwlflZrxpqUFth+ivK+I0yy3j+s013yCa6HKc0+Sxn041uthOo7jOM4XxR+YjuM4jlMAf2A6juM4TgH8gek4juM4BaiauICTfJiomRN0mCSAk2WkONFy3759jeYkH058YNJxJsXebbfdjN5zzz2jNtx9991Gr0JiXk7geOKJJ4xmgD+TLXCyEydDSfHkJRa1ZpJrJk844YQTjN5oo42M5rnKawOTQrDwNdvQWnzve98zmsWfOfmF4ywvAQDHCQuZM0lDc2Ofk9cY0H/YYYdFbWC7mVyCE2w4MYnjioWWf/Ob3xjNiUxSPLGLE5F++tOfGn3PPfcYzWTsu+++u9FMbsJrKW+fHJucKHTaaadF2/iyiYoxd+9q5OoedsJN7Xx7vUYTdqQ4UTmTiCMxOos9C0UEokTqmIAVZuYUGMA4Lc2zE2jKPTCJBxN0EujQYJMIMJF6bhHrZpLIRwkdmNgeRa7DUiS6x/Zqcop5NzfJp9nlGf4L03Ecx3EK4A9Mx3EcxymAPzAdx3EcpwBVPUxCr4GeF71AFriV4gDwiRMnGs2AcHqWLBZMmEw6DyaFp2/KorcMpD7nnHOMvvPOO42mH5mXpICJBBhAzuB8BojTh+U+mbhg7NixURvOO+88o5mogOfmmGOOibbREjCxBP3ZmTNnGs3EEfTIpNiTZLIKJvymP8/+3XDDDY2+6667jP7xj38ctYGeJBMsfOtb3zL6vvvuM5oF3Fko+8wzzzR6+PDhURu22moro5m4ZPLkyUYzqQD3ycLcTEqflziE1ziPi+e7NUhQSJk+Wu1iFEFmsedV9vqUmi++zMToSVsE3C+0wfbJgoX28xvYguL0KyVJvPesQSKKRXZcR0kFOtr7tZgonb4sl+cQFYxusG1IFlmfVTX0fpG4AG0ufxonBCn1sl56lHSe/nET+C9Mx3EcxymAPzAdx3EcpwD+wHQcx3GcAlT1MFlImT4O/QkmAGesmxR7mIReHZNgs/AvfVOun5cAnsWU6V8xDpDe3i9+8Quj2Q/0qhinKcU+HItOjxkzxmjGs9FD4zEwVjQv+frNN99sdJG+aw2+853vGD116lSjGXe5+eabG92eHofiOMzjjz/e6GuuucboUaNGGU2f9LrrrjP6jDPOMJqeqRQnu99xxx2NZsJ9xijzfHA5vXL6qlLsD1900UVGc0ywzbx+6b3T46SW4oLr9MqvvPLK6DMtDT1I+omRR1lr+zGK45Qi7y20sduMYhiZpHx9ew2HuQurtinpnlPEehGSpzOucjn8fvqg6IekDv3CotbcnqSks/W1mUw9isPkmIHn2VwcZ6m+c9QGFtsOPC56sU3gvzAdx3EcpwD+wHQcx3GcAvgD03Ecx3EKUNXDpC+2ww47GE3Ppbn1JenZZ581mnFb9EmHDh1arYlR3kkWzeX2JenWW281mu3eZpttjGZMJH0g7oM+K3PsSs3n0GUcJds4a9Ysoz/66KNoH5Xk+ZH0ZukNNhfz2lKwwDDj9ubPn2/0008/bTT9Rkm64IILjOY47NTJ5gqlH/irX/2q6vqM/fyXf/mXqA29evWq2s4pU6YYzbFcX2/zfk6fPt1oxgNffPHFURt4zTLWk21gMfSPP/7Y6JNPPtloXr/sZym+Pph/lnl6W4Okvb2mwwrrDwbECuoTW8A9dOsab7MeBZ+R4zRAl+uRM3WVHYP03SK/MM9HRXwpP5Mw321nxKPCJw2C14f4xWR1HI+aF6NaDXqeWoxi3ozLxDyWyBvOawPijxl32xT+C9NxHMdxCuAPTMdxHMcpgD8wHcdxHKcAVc0CenWPPPKI0QMHDjT6008/NZqxgFLs4zDua++99zb6kksuMZp5XQnjMul3SXHNTObEpfdK34Y+Desv0m9kvk1JevXVV41mzCuPg31NT5IxfWxzHm0QF8Z2MwavtTjuuOOMZj5henscU+wrKfYLWSuSHuTZZ59tNH1sesZsw9VXXx21gbUkec7nzrW+GOOBJ02aVLXNrMk5e/bsqA3MDdvQYOP0dt55Z6NZT5Z+/Q033GA08w3n+ZH77LOP0cz//FXkkg3LEaeHuD3GL6psPbBkGWo0Ko67pG+WLLZ9T580wItLmCOZuWhnWX9akkJn+KhodxnzA0qM5YRfmMy199PQDeujjqgkJcjtGuD/03sN9FWRzzaKw6R3m1PbsrGn9f9rFiyN1imC/8J0HMdxnAL4A9NxHMdxCuAPTMdxHMcpQFUPs7ncscz5yTgwemRSnNOUNfzoH3Kf9OboqTDOjD6sFHtBzEXKWDTGJxJ6MPSi2I+S1K9fv6r7YP5U+ouMT+U+GbtI30+S1l/f1tNjvdP33nsv+kxrcPvttxvdnMfFmMk87/z+++83ettttzWafiD1ypU2Zo6xnz/72c+MZv9LcX/TX2ee1gcffNBoxn4yn/Att9xSdbkkDR482Oi33nrLaNamXL7cenO8xpkn+ZlnnjGa/r4U9x39evqmrQLzm9JXa2/9xDWf2PkaUU5WSaUuNqdpgJ+bwCeN6GLPd+SJkk5xzDnrnYaONs9yDY6THmWZcZiMBeUx1MX3W6EWJWM3y2hTaQnqiJL59l4WGGtKj1NSabEdx5EPSo+6CfwXpuM4juMUwB+YjuM4jlMAf2A6juM4TgHWycOkN0Q/gz5bnn/I+EP6MozlpK/G+MMTTzzRaNaVpP8lxflnGcNI744xd8zByuOkF9WzZ8+oDaw9ybhAeprMb0uflfFuPFf0iSSpe/fuVffB2MPW4txzzzWaeVzpu2233XZGv/LKK9E2Oa5Yj/Rb3/qW0TznjB388MMPjR4/frzRzJcqxb40fVHW8aQ/yLhYere8lvLif1mnkzmN33zzTaN5XPQk2UZ6v3nngvGlL7/8ctU2tQoJPCzUu6TvVtPFxvUpxL89Auuywi8s0e9jfct5C+1yeKIJfVPE+ab7RN5VeI6N69t7U81sxEzyHl6q/hsrLI7rwDZXzzLMsftMcF0kOK5SVzsGE8QSK6f+MPuu3MmuU1rhuWQdx3Ec50vDH5iO4ziOUwB/YDqO4zhOAQLjdCoZMWKEWUhvj3FejIHMywnJ2M2HH37Y6C233NLoG2+80Wh6KIzLZNxmnh/C2M3JkycbPWTIEKPpFdGLImxDXk1OerOMk2Qb6MUy/y19UvpZ77//ftSG1ahdx/g3+n4TJ060yStbCI67I444wixnvUz6jXneOdliiy2Mpl9/+eWXG3366acbzXHF/qY/KcX9SR908WIbCzhmzBijec6ffPJJow888ECj82pRMuaYsaBlxKPxGm8u7pKx1+3axbUJCecQMEZ51KhRLT7uDu5xqhlz5SV2PNR8A/MQapFjdUV8vkMzOVBZg7E8294/S/XWdyPlpTZeMeT0dUA7xVjQdvZaKXdCnOZCeJKMV0Q+2yguU4pzuzK/MPuyAfvE+lHcJX3VNbGXyxqZ5Q52G/Rux828JnfM+S9Mx3EcxymAPzAdx3EcpwD+wHQcx3GcAvgD03Ecx3EKUDVxwSabbGI0J69wwg2D5e+6665omywQzQkyxx57rNGc7MJ9MLCayaLz4HFwIgSPixN0mJjggw8+MJoB/3kJADjRiHBi0YgRI4zmRAlOrmAwOCf4SHGwPiedfFUFpDnJaty4cUZzgtSAAQOMzpvgxGLKP//5z43+yU9+YvRFF11kNIs/c1xyTOT1NxNF9Ohhk1JzHxwjHNuceDR9+nSj+/fvH7WB65x11llG33HHHUavt956RnMiUXOT11igWoqTdrzxxhtG8/yNGjUq2saXTUDi8poOSDqAiSUJJ9xggo+keILMQiR4xz5KSFIeTRJaas9/6GMLkOeyzCY3KCOhewnbrFnUTGFlJkLISTofwUk6IGlvlwdOoGLSeU7yYWHuujhJffKxTeJRw/PF890E/gvTcRzHcQrgD0zHcRzHKYA/MB3HcRynAFU9zF122cVoJj4nDNamH5m3Dr2f0aNHG01PhMkTmHi7uQTXktSnTx+jmXCawdkMOGdgPL0iLmdSeiku+MyAbyZooOZxMvk6k3/Tj5akadOmGc0E7XmfaQ1Y3JmFlIcPH270rFmzjM4Llmf/XnLJJUYzKXkJPgkTbrzwwgtGz5gxw2gmW5Ck559/3ujAoG9Aj5Le3gknnGD0X/7yF6NZKF2K/fSrr77a6AsuuMDo//mf/zGa/cBxyHGXlwCeyRGo99133+gzLU3ko62KPWizHInVo4LEkkKjPa5AHxTJ2IXizPQshSLWUZHrRTmJz5EUoLTQesrRcbN4M/zHyFddhaQDHXISVVRJjiNJYgKHntbjDkyGgGMKK1HkekWcPCH0tPMHuM883zMP/4XpOI7jOAXwB6bjOI7jFMAfmI7jOI5TgKrJ1++++26zkF4SPc2xY8caPWHChGibLP5L723ePJsElzGT9IroX/Xq1SvaJ6E3x2TPTChNmJB6++23N5oxeHkJ4OmpMa6PsWnsh+YKcTO+lYnxpbjdTz31lNE8jjlz5rRK8vVBgwaZcUc/kf1J34yJ6yXpyiuvNJr9WYP4MiZnZ6J0+qSdO9tE2Yx3lOKYRPqJHPuMR6V/yKLjeX4h4fXGsVxfX19V029kv7AIAOcgSHHff//73zeaSeYvv/zylk++vuWF9kbYtpin9Rl5Sccxppjwu9zenouaJXa+xZp6O8YCc5gvQCzosjgmMmHydVxL9GrpaTYb81iL7a2Ivd/AZOjYZ7mrvXYaOyMx+nJsE23I22cEP7McRarr7Zj886v5Y85/YTqO4zhOAfyB6TiO4zgF8Aem4ziO4xSgqofpOI7jOE6K/8J0HMdxnAL4A9NxHMdxCuAPTMdxHMcpgD8wHcdxHKcA/sB0HMdxnAL4A9NxHMdxCvD/AZdUdFxh4wwgAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Detected pixel offset (y, x): [12. 11.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAChCAYAAACWEznEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiJ0lEQVR4nO2defxVc/7HX59IQoUsZZlEFFqUkFGjkGhQ0jLSQ/ZCI8ZazG+mmSbzq/lNksmSsStLZMnka2tSCmkaQraRFlMRUchW5/fHPbjv5+d0b5i+5ny9n49HD16de8/53HM+99M9r/NeQpIkchzHyRvVfugBOI7jfBd88XIcJ5f44uU4Ti7xxctxnFzii5fjOLnEFy/HcXLJf+3iFULYPoTwZAhhVQjh/37o8XwXQgjXhBB+/UOPw/l+hAI3hhBWhBCeTf/uzBDCshDCRyGEuhvw2D6H1kHYkHFeIYS3JG0vaY2kjyQ9LGlAkiQfrcd7fy2ppaTjEg9GczYgIYSTJJ0vaTdJKyVNlDQoSZIP0u3tJI2X1DhJko9DCNXT17VJkuT573Hc9pJuS5Jkp+8z/h8rlfHL6+gkSbaQtI8Ki9Gg9XxfA0kvf5eFK4Sw8bd9j/PjJIRwvqT/lXShpDqS2qgw9x4NIWySvqyBpLeSJPk41dtL2lTSS5U8XKeYJEk22B9Jb0k6rEgPl/RQkW4jaYakDyQ9L6l9+vc3SfpC0ucq/GI7TIWF9hJJ/5L0nqS7JG2dvn4XSYmkUyUtlPRk+venSJonaYWkCkkNio6dSOov6fX0+H9R+ks03X56+t5Vkl6W1Cr9+x0k3SPpXUnzJZ1T4vPfJGlo+v/tJS2WdJGkdyQtkdRVUmdJr0l6X9LgovfuL2lmOrYlkq6StEnR9sMlvSrpQ0ljJE2VdFrR9nV+dv/z9Tmqnc6vnvj7LdLre0o6pz7VN3cP4yV9nM6fjyQ9ISlIGple15WS5kpqmu6rhqQ/pfNymaRrJNWUtLmk1ZLWpvv5SNIOPoe+xfXbwJPjLaWLl6Sd0os6KtU7qrAIdVZhYeqY6m150VI9UNLT6X5qSLpW0vh02y7pZLolnRQ1JXWR9IakPSVtLOkySTOK9pdImiRpS0k/SSfrEem2HpLelrRfOjEbqfCvbzVJsyX9j6RNJO0q6U1Jndbx+TnxvkzfW12FxfFdSeMk1ZK0dzqZG6av31eFxX3j9PPNk3Ruum0bFb4k3dLtA1VY7E9Lt5f87P7n6+tzRHpNNs7YdnPR/DpJ0vSibV/Nt41T3SmdF1um82VPSfXTbSMlPSBp6/Q6Pyjp8qI5sbjMGH0OrevcbODJ8ZYK/6KsSi/245K2TLddLOlWvL5CUl9etFTPk3Roka6fnuyvLkwiadei7ZMlnVqkq0n6ROm/Hunr2xZtv0vSJUXjGJjxeQ6QtBB/N0jSjes58VZL2ijVtdIxHFD0+tmSuq5jX+dKmpj+/4mSZhZtC5IWFU28kp/d/3x9XvpIWrqObX+U9Gj6/yep9OJ1iAq/fNpIqobr8rGk3Yr+7kBJ84vmxLddvHwOpX8qw/PqmiRJrfTEN1FhxZcKv2R6hBA++OqPpLYqLEpZNJA0sei181T4Kb990WsW4fWjil7/vgoXaMei1ywt+v9PVLhdkKSdVbg9zRrDDhjzYIyhFO8lSbIm/f/V6X+XFW1f/dUYQgh7hBAmhRCWhhBWShqmb87dDsWfNSnMrsUYZ7nP7kjLJW2zDo+0frq9LEmSPKHCLdlfJL0TQrguhFBb0raSNpM0u+haPJz+/XfF51BKpYVKJEkyVYV/Rf6U/tUiFX55bVn0Z/MkSf64jl0sknQkXr9pkiRvFx8Gr++H19dMkmTGegx3kQpPnrL+fj72WStJks7rsc9vy9WSXpG0e5IktVVYJEO6bYkKt8+SCo/yi7W+32f/MTFT0mcq3Dp9TQhhC0lHqnCnsF4kSXJlkiT7StpL0h4qPABYrsJisnfRdaiTFB5gSXa+bgiq9Byq7DivKyR1DCG0kHSbpKNDCJ1CCBuFEDYNIbQPIazrsfE1kv4QQmggSSGEbUMIXUoc6xpJg0IIe6evrxNC6LGe47xe0gUhhH3TGJ9G6XGflbQqhHBxCKFmOu6mIYT91nO/34ZaKngSH4UQmkg6s2jbQ5KahRC6pr8azpZUr2j79/nsPxqSJPlQ0hBJo0MIR4QQqocQdlHBQlgs6db12U8IYb8QwgFpCMXHKhj8a5MkWStprKSRIYTt0tfuGELolL51maS6IYQ6/9EP9g1Veg5V6uKVJMm7Kpjq/5MkySIVTMHBKpiOi1T412pdYxqlgvH5SAhhlQrm/QEljjVRhUfgd6Q/mV9U4V/T9Rnn3ZL+oIIRukrSfSo82Vwj6SgVwj7mq/Av6/UqPGL/T3OBpN7p8cdKurNofMtVeKgwXIWHHHtJek6FXxHf67P/2EiSZLgKc/BPKnzRn1FhLh6aJMln67mb2ipcoxWSFqhwTUak2y5Wwfh+Or0Wj0lqnB77FRWeXr6Z3p7t8B/5UN9QpefQBg1SdSqHEEI1FX4pnJAkyZQfejxO/sjjHPqvTQ9ySpPebm8ZQqihb7yMp3/gYTk5Iu9zyBev/HKgCk9El0s6WoWnuqtLv8VxDLmeQ37b6DhOLvFfXo7j5BJfvBzHySUlqy8MGzbM3FO+/vrrZvuHH35odJ06NmLgs8/iJ801atQw+sYbbzT6lltuMXrOnDlGv/3220Yfe+yxRh9//PFGDxgwIBrD3nvvbfTHH39s9KOPPmr0wQcfbPROO9lQtEmTJhm93XbbGV29evVoDDwm3/PBBx8Y/fnnnxv973//2+jNNtvM6IMOOsjot956KxrDjjvaYOl333235DFHjhwZVAl07NjRzLtGjRqZ7U2aNDH61FNPNfrZZ5+N9nn55ZcbvdtuNgZ5hx1slMIjjzxidO3atY3edlsbJL9s2TKjt9pqq2gMvKYcw+LFi43m92v//fc3+vnnbTUeft8+/fTTaAwbb2y/8pybm2yyidH77LOP0b/97W+Nfu6554weN26c0V26xKGY/Jz333+/0fx+jRkzJnPe+S8vx3FyiS9ejuPkkpK3jVOnTjWaPzEvvPBCo6+99lqjt95662ifu+yyi9E9e/Y0essttzT6iy++MJo/YydOnGg0bxl4PCm+heJrePvLn8Z/+9vfjO7c2aY28paOtxRS/JP+5ZdfNrpr165Gv/SSrXvH80B4i8H3S9I//vEPo6+88kqjb7/99pLH2FC0atXK6BkzbDrdzjvvbPSwYcOM5vWTpP32sxlcL774otG0M1atWmX05ptvXmLE8W3/DTfcEL3myCNtgDpv4T76yBYY5jH/9S9bK4C3V61btzZ6woQJ0Rh4a9m8eXOjH3/cpnPyGL169TKaawKtiDfffDMaw9y5c43meVi7dm30niz8l5fjOLnEFy/HcXKJL16O4+SSkp7XXnvtVfLNs2bNMpo+Dh+1S3GIAO93+Yj5nXfeMZqPj+k1cP8NGjSIxsBwC2oeg/C80BNjaMUzzzwT7aNjx45Gn3nmmUafd955RtPjuu+++4z++c9/bjRDULLg9Ro8eLDRHTp0KLuPDQF9S4a2TJ482eiFCxcanTVv99hjD6PppT7wwANGc97QxyQ1a9Y0+pe//GX0Gp7vefPmGf3ll18azbAFfp84J+bPn290w4YNozGsWLHCaIbY0NOiV7rRRhsZ3b1795LvP+GEE6IxMGSEn/OVV16J3pOF//JyHCeX+OLlOE4u8cXLcZxcUtLz4r0pfZirrrrKaPpNWTCWJSsWrBh6DfRyhgwZYjTTYp566qlon4yBYmoIY1cYF3byyScbPXr0aKMHDhxoND0xKY6hYhrT+++/bzSvBf2OcrFrTO2SYl+H149jZOrVhoJxWkyDYWoO003ooUjSa6+9ZjTjCZlixHnTrl07oxm/RJ9t2rRp0RiY5kRvdffddzeafhO3Mx6K1y8rLY1/94tf/MLo1attRRzGZDH2jF7hq6++ajTncdY4OVebNm0avScL/+XlOE4u8cXLcZxc4ouX4zi5pGQl1d69e5uNzCmjL0DviCVWpDgfkp7WP//5T6Oz7plLsf32tv8rvSIpjg1jLBnHfe655xrNOCR6F4yXuvTSS6Mx0LfhOJnbyPg3xvQw5o7wM2Ydc8yYMUafddZZRt91112VUhKnS5cuZt7RF120aJHRzNfLuubTp083mr5at26mdaPatGljNEvktGzZ0mjGT9FTk+JxMz7tZz/7mdGMqeI1p7dHv+qTTz6JxsBcxuXLbV9d7nPTTTc1ml5f1ne83Pb27dsbzblZrZr9TTVhwgQvieM4TtXBFy/HcXKJL16O4+QSX7wcx8klJYNUaSi+9957RjNQkgGiTH6VpOuvv97oHj16lBzgIYccYjTNchqnHMNdd90V7bNc/XgGsTIAlMUICR9KsECbFCeDcwwjRowwmkYqDWc+NOjbt6/RrMuf9Z5yDyIqizfeeMNoBojyYQcTebPmHa85a9Kz+OCDDz5oNJOqadD//e9/N7pWrVrRGBjgye8Pk6QvvvhioxkUzmKS/EwnnXRSNIabb77ZaBZhZMDoOeecYzS/v3feeafRW2yxhdF8yCfFDxZI1sOOLPyXl+M4ucQXL8dxcokvXo7j5JKSQaqjRo0yGxlAWu4engXzpNi7YWAr73cZFMd79MMPP9xoBs5m9SvkMfk5fvrTnxp99913G00Piw02eN+f1byBHslFF11kNIMimZRbrjgePZ6sMZQLAGZw6FVXXVUpQart27c3827NmjVmO3sB0nfJ6lfIv2NSM5tXMHCSBQVYOJBJ7kuXLo3GQA/qscceM5oJ5pzbdevWNZpBrWxmQ09Mkvr162f0kiVLSh6D3h3nFc8LvwtZxUBZrJPnhYU5x40b50GqjuNUHXzxchwnl/ji5ThOLikZ58ViaPR66C/RQ8kqCrfNNtuUfA+PyQL+9N3oBbVo0cLosWPHRmNgjBVjW9gAlsmljH9iQ1Pun802JGnBggVG81xyTCxoyGYMjBEiTF6WYt+GvlhWMndl0LZtW6NfeOEFozlnmIDMZGJJOvTQQ42ePXu20ZxXjD9kjCOTpDknfvKTn0RjOOWUU4zmuBlbtueeexpN7/Xhhx82etdddzWazZCl2NOiv8TGtowPrKioMJoNN0jWteC4br31VqOzPMss/JeX4zi5xBcvx3FyiS9ejuPkkpJxXgMGDDAb6YGwYSzzmLIaADC+iT4L30PfjPE4zAmkd8SiflIcf8bYFXpBvG9nwUN+BuY20leQpMaNG0d/VwrGqzHOix4YcyNvuummaJ9sysH4thkzZhhdUVFRKXFeLIJZr149s52eSKdOnYxmDqAUN2BlET56WIxPYoE8zhF+F+hpSrFfxHjAI444wmj6mMzx5JjYnIbxjJK02267Gc0cTZ5bNuulJ8aCia1atTL6sMMOi8bw5z//2WjG3LHJ8JAhQzzOy3GcqoMvXo7j5BJfvBzHySUlPa9u3bqZjaeddprZTh+FtbSy6l7RSyiXo1cufoleA+/ZGVcmxQ0aeI/NWDPCGmFsGEsPLKvhK6GvVu5zMKds5cqVRjO2LCvHk74Ma4Txc1RWbuNpp51m5h19k/vvv99oekchxMNkE1jOQ84JxsUxxpF+EptrMD5Rir8vffr0MZpNQthst3///kafeOKJRjPXkdulOGbu6aefNvqYY44xmr7ngQceaDTj49jE9o477ojGwJxNxqtxXs6dO9c9L8dxqg6+eDmOk0t88XIcJ5eU9Lx69uy57o2K6/DcdtttRjNuRYrrpDOHj94M474Yk0VvhzmATzzxRDSGY4891mj6bjwGcxvL5TIyNo1eoCRNmTLFaHouRx11lNE8b+XGyBr2U6dOjcbAGDh6efQeKivO65BDDjHzjjXeJkyYYPTQoUONZs0qKW50evbZZxvNecL8VsY70YfjHGDNKimu/87cRn4Xa9asaTR9TXqp9PGY6yjF9f3pGdMTo/fH+DbW+2JNv6x4N8YkMm+Un/Ptt992z8txnKqDL16O4+QSX7wcx8klJT2vDh06mI3169c321njnHFBvLeV4ntgxsswf5K5i/TZmKfI3MasmmI8JmuC0aNiniD7BhLGjb322mvRa7baaiuj6e2xdj99AHqF9O0YN8b+l1J87srFOs2aNatSPC/GF7J+F+tiTZo0yWjGXEnx3GRvSF4P5gDS22GPxN69exs9atSoaAyM3eMx5syZYzRrgjVt2tRo9mCk58VelFI81/n9YNxW9+7djX7yySeNZl4ivxsvvvhiNAauIzwvHPe0adPc83Icp+rgi5fjOLnEFy/HcXJJyRr2pFw8ButBDR48ONoHe9U99NBDRp9xxhlGM5Zl+PDhRtMT4z38XnvtFY2BMTust8WcTMbGMO6Iml4Sx5i1T8ZplcvppEdGf5FxYFm5jRxDVv21HwLOK+Yu3nvvvUY/88wzRrP2uxRfA35W5kNyjtC74f6GDBlidFYNe/qS9HfZ65GeI68hfTrWiMvqHTl37lyjW7ZsaTTre7H+XvPmzY1mjftyObhS7CnzXDJ3dV34Ly/HcXKJL16O4+QSX7wcx8klvng5jpNLShr25cxvmnk0JHv27BkfEImgbGZBU5KFAZkYSuOViaFMaJbiAE+arTTguU8+ROB5oumZ1QiBZiyDB0eOHGl0hw4djGbQJQ16BgtmBcrymAyuZaJ2ZcHARzZ9YANZGryNGjWK9slAVz4gYYDnrFmzjGYALwsF0qDPOt/NmjUzmk0/2GiEyd6c+2x+we/KsmXLojHQUOcDM5rpfEjHZHIGnI4ePdrorKKMgwYNMppFCNasWRO9Jwv/5eU4Ti7xxctxnFzii5fjOLnkWxUjpCfCIEcmEzPhUoqbeHTu3JnHNJr3+YT+Ev0n+jpSnHxKb4DeBPfB9/Nz77TTTkbTL5Hihg4MSuV2ejD05cp5XlnnsVyTDp7b6667rlISs5s1a2bmHZurUtOHoU8jxdeQ3isDRnnNVq9ebTTPL72iLM+LyeBHHnmk0byG9LAY8Mk5wAYcTKKW4oRynrvzzz/f6Hvuucdo+sWcpxwzCypK0pIlS4zmuaef+9hjj3lituM4VQdfvBzHySW+eDmOk0tKxnnx3pOFABnXxRiPrIavTFpmo4MPP/zQaCYc8567RYsWRtMbymogS6+A8TT0Q+hxMW6LCekcA71CKfaXmAjMGDrG2NGTKefLMVZNis8tdVYhx8pg4MCBRvMaPvLII0bzs2cVi+Q1pi/J4o1MBqdmfNPkyZONbtKkSTSG/fff32heI847HoOfgcVAGzZsaHRWcjiT2umB0eNiU1rGYDFOjEUbs5LDWYTx1VdfNbpXr17Re7LwX16O4+QSX7wcx8klvng5jpNLSsZ5DRkyxGxkPli54nXMAZTiWBbGinE742cYE0KfgLFlderUicZAr45eD+Np6GERxqbR22DBNiluXFuu8QjjhngteJ54XrOaUtDHoY/G2LB77723UuK86tWrZ+bdzjvvbLbTs5w9e7bR9CCl2NNic2J6OTwGzxXzJ1lAMavp7MKFC41mTma5Yy5evNhoxmgdeOCBRi9fvjwaA+cV5wV9an6fWACxT58+Ro8bN87oVq1aRWMYO3as0bw29MlWrFjhcV6O41QdfPFyHCeX+OLlOE4uKRnnxfga3pMzvobxOcOGDYv2Sf+HsWTlcvbow/D98+bNM7p9+/bRGLjPLG/u20D/ibEv9J+k2OOaMmWK0eX8D8aOMbaJ14I1y6Q4voZ+CHMfKwvWqaJmHihj5BgPJcWNIxiPVK62Fo/B5sb0XtkURJKOP/54o9n8gvvg9aH/1K5dO6M5zxgzKcV+bLmYR3LLLbcYzaa09E2vuOKKaB+shUZfeubMmSXH8BX+y8txnFzii5fjOLnEFy/HcXJJyTivo446at0bFXsP9EwYTyXFMVWspUUfjff1PXr0MPqss84ymv5TFl27djX6pptuMpreA+/j6T+xySnrFWU13iSMZaInxrxS5ojS72CsE/PepPhc81zOmTOHb6mUOK+GDRuaecd5xrpV06dPN5pxQ5JUt25do+nt8D30xOgflcuHrVGjRjQG1nNnTBQ9LV4f7pNj4PuZvynFTZ0Ze8bmuox/4/eXsWdcA+hJS/F3lPXXWCNs0qRJHuflOE7VwRcvx3FyiS9ejuPkkpJxXvQaeM/OnD/W5mLtIEnq1KmT0byPp+7SpYvR5fwp5vhxTFJ8z01Pir0k2VeOHgtjy+jLsXaXFPuBrBfP+l2s/c99sjYU45LuvvvuaAw8d/QreK0qKiqifWwIGHtE+NkPO+wwo3l9pPia77rrrkYzbo69I+kf0ct55ZVXjGa+pSS98MILRjN+LSsPtxj2r2R+LGPLOCek2BdjLOeMGTOMvvbaa41mfS/GXTL3kZ6YFHvG9BOzxp2F//JyHCeX+OLlOE4u8cXLcZxcUjLO64wzzjAb6Q2Vi0vJyhmkz8JcK+YiMj+SvgBr2tMPyfJPWJu9devWRl9//fVG039asGCB0ayJlBVTRcr1o2ScF30c1mpi7Brr7PM8S3H82sSJE42mz1NRUVEpcV5r16418+5Xv/qV2c5rzj6ODz74YLRP+kOMzWNtLXo39H/pUTKO7LPPPovGwO9LVn/JYphbSs+LNex5zKzcRtbr4vdl1apVRtOX47lmfCE9afYbkKQnnnjCaJ4H1kLr37+/x3k5jlN18MXLcZxc4ouX4zi5xBcvx3FySckgVRr0TOKkuc6AUBrdUtxIgkmYNOhpWjIhmYmkBxxwgNFZzRjGjx9vNJNNab6yASyDORlYywcVfJAhxQ8eWCCRidgMlO3bt6/RbDDKa5U1BgYAs4gix1BZnHzyyUazkCCNa86zrGBPzhMWxWRAbrm5zwdPDN485phjojFw3AwkpjnOhwqcVyzad8MNNxjNhxBS/FCGDxEuu+wyo++8806jmajdpk0boxnIzu9S1jE5N2ny9+/fP9qH5L+8HMfJKb54OY6TS3zxchwnl5T0vAjvTemR0DtiAwEpDvZj4wkG/9HjevPNN0uOkYmmWTBhnD7bhRdeaDSD5s477zyjb7/9dqPLNQ2R4qBRBgsyEJPBgPTteEwGqWY1QrjggguMZlAqrw0bSGwoGERMP49F+BgkTE9Fij0sBiYzGZh+Ls8vi0eysOD5558fjYEeFoNpWWRxwoQJRrPJLAsoDhgwwOhu3bpFY2jSpInRDFJ/9tlnjWYAKY/JYqLlGupI8Xecn4vXe134Ly/HcXKJL16O4+QSX7wcx8klJT0vFuXjfT/vZ5kczFgaKfa8CL0dJsiyiBx9Nr4+KzmchfnodzDOiF7Q73//e6N5HuhtZDVjoG/DAobDhw83mvEy9Fz4GRiLlpWYfeONNxq9PueuMjj66KONZiMQxnU1atTI6Jo1a0b7ZJzXiSeeaPTo0aON7tevn9H01caMGWM0m5dkNV1hIjybzjIZnzGQvB7cTm+VPpwU+4mDBg0ymnOCY+b3l14tPbGsZigs3klvdcSIEdF7svBfXo7j5BJfvBzHySW+eDmOk0tKel70UfbZZx+jeY9e7vWS9NRTTxnNuBD6ap07dy41xCgPigXYsppe3nzzzUZz3GzGwJgr+gY8Bn25rCYg5XI6GafFMbKxbVajg2Ky/Ct6efSSysXUbShYrI5xQe+//77R06ZNM5r+lCRddNFFRnMesikE/aOrr7665OsZW3bcccdFY9h2221LjpMNbTiXa9eubTSbfjDecPDgwdEY+J1lLBnHwMKabBxz6qmnGs3vL8+zFH8/mA/JvNF14b+8HMfJJb54OY6TS3zxchwnl5S8uaS3M2nSJKPZ1ICNPRlrJMX3/Ywradu2rdGXXnqp0cwzJIz7oj8ixTW/mKNJr473+byvZ/0o+lPM/5KkuXPnGs2YOn4Onmt6WIwZ4pizqF69utEcN2N8Kos+ffoYzfxWekGcUzxXUuwvsdYVPauBAwcaTd+THiPHMGrUqGgMrIXFa86mKow3ZFNZjpk1xZYuXRqNgbmKbLix7777Gs16ePR3r7vuOqOZ/5rlX7Vr185o5iN7bqPjOFUaX7wcx8klvng5jpNLSjadbdu2rdnIppZsBso4k6w4L+YBMjeKcDt9GN6D04ugbyfF9+Hc58qVK41m/BPzLxmHRO+CHpkkNWjQwOiDDjrIaObzMa6L8W88JmupZTUgrV+/vtH0zRij89JLL1VK09mBAweaeVfOE2HuaIcOHaJ93nPPPUazYSv9I2oekzXhfvOb3xj95JNPRmPg+aYfy3nFecjYsnJNoOmpSXFNPvaUYG0tzkt+bp77evXqGZ3VT4B5vfR3WYesd+/e3nTWcZyqgy9ejuPkEl+8HMfJJSU9r5YtW5qNjNsqFyeUVUOK8U28j2esGO/RmSfI3oJ8fVZNe+ZDMpdx/vz5RtMTY04gt5ergSSVr51Fn43eHXPU6OPxvGblPjLHjD4avbrBgwdXiue1YMECM++YV0ifplmzZkYznkqSZsyYYTR9FtaP5/nq2bOn0fR3ef55bqXYx6SHzDpk9JP4XWX+ZTlPTIo9YubYMpaM85Sfk2MMwU4RXispjl97/vnnje7Vq5fRBx98sHtejuNUHXzxchwnl/ji5ThOLinpefXs2dNspBfEOur0YbJylFifnL3rGjdubPTYsWONZtwIc/gYR8Z7eim+72evuo4dOxpNb4HeBeEYsmqK0ZtjHBbHQO+O+ZjbbLON0fQ/suLpGFfEPDb6j1OmTKkUz4vzrnv37mY7632tT2wf2WOPPYymvzd06FCjzzzzTKM5r3i+6WdJ8flcuHCh0YzrYh8DXvOpU6caffjhhxudVUuL9f4Za7Z27Vqj+R2nxzV9+nSj6YNn1bAnjC9kDGS/fv3c83Icp+rgi5fjOLnEFy/HcXKJL16O4+SSksUIWXyfxjPNcgaxZjW9ZLFBmtu9e/c2mkY1j1FRUWH06tWro2MSfg6amPxcNNcZdLpgwQKjWbiOWoofEhA+FGCQJE1OGqMM/KM5L8WBmTSMf6hihHxAMnnyZKP5cKN58+ZGZz2cYGG+3/3ud0ZfcsklRrMZKwNfOS85J7LON5PB2UiYx+Ac4dzmQwM25GDz46zXnHPOOUbfdtttRjPRmg8Byj14YrFDKQ58ZeA5rx8bAH+F//JyHCeX+OLlOE4u8cXLcZxcUtLzat26tdFZjUuLYWAe/aus19ArYPNV3kMzUJZJueWSX6W4oSsL9zEQj8GFDIKkt8DtTFiX4uKBDO5jMC41P2e5xGD6l5I0b948o5msnPWeyqBp06ZGT5w40WgWq2OhxqzASJ5fNnZhEnO1avbfdQZXs9kFG3wwsFaSZs6caTSTmAk9LXpBffv2Nfrxxx83mknWUuy/slEIm/P+9a9/NZrngfOQ8y6r+QwDYanLNdn5Cv/l5ThOLvHFy3GcXOKLl+M4uaRkYvb48ePNRnoP9MCuuOIKox999NFonyxORq+GTTEZk0VvgX4Hm39mQS+HiaBMNiVMVm3RooXRjPHJSg6nB8O4Ica+8DyUK+rI+DkmzUvxuNk0gp9j2bJllZKY3b59ezPv6D/xfNJnYVK7JI0YMcJons+NNtrIaCZuM4mavlqtWrWMZjyVFMc80X/i3Ge8G/0mNs7N8pcIv2+cy7Vr1y6p6U/xvLBAQFbjF577008/3WgmoA8dOtQTsx3HqTr44uU4Ti7xxctxnFxS0vNyHMf5b8V/eTmOk0t88XIcJ5f44uU4Ti7xxctxnFzii5fjOLnEFy/HcXLJ/wOKWh4h4u0UKgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Detected subpixel offset (y, x): [-12.029 -10.932]\n"
     ]
    }
   ],
   "source": [
    "image = a[0,:,:]\n",
    "offset_image = b[0,:,:]\n",
    "# pixel precision first\n",
    "shift, error, diffphase = phase_cross_correlation(image, offset_image)\n",
    "\n",
    "fig = plt.figure(figsize=(8, 3))\n",
    "ax1 = plt.subplot(1, 3, 1)\n",
    "ax2 = plt.subplot(1, 3, 2, sharex=ax1, sharey=ax1)\n",
    "ax3 = plt.subplot(1, 3, 3)\n",
    "\n",
    "ax1.imshow(image, cmap='gray')\n",
    "ax1.set_axis_off()\n",
    "ax1.set_title('Reference image')\n",
    "\n",
    "ax2.imshow(offset_image.real, cmap='gray')\n",
    "ax2.set_axis_off()\n",
    "ax2.set_title('Offset image')\n",
    "\n",
    "# Show the output of a cross-correlation to show what the algorithm is\n",
    "# doing behind the scenes\n",
    "image_product = np.fft.fft2(image).conj() * np.fft.fft2(offset_image)\n",
    "cc_image = np.fft.fftshift(np.fft.ifft2(image_product))\n",
    "ax3.imshow(cc_image.real)\n",
    "ax3.set_axis_off()\n",
    "ax3.set_title(\"Cross-correlation\")\n",
    "\n",
    "plt.show()\n",
    "\n",
    "print(\"Detected pixel offset (y, x): {}\".format(shift))\n",
    "\n",
    "# subpixel precision\n",
    "shift, error, diffphase = phase_cross_correlation(offset_image, image, upsample_factor=1000)\n",
    "\n",
    "fig = plt.figure(figsize=(8, 3))\n",
    "ax1 = plt.subplot(1, 3, 1)\n",
    "ax2 = plt.subplot(1, 3, 2, sharex=ax1, sharey=ax1)\n",
    "# ax3 = plt.subplot(1, 3, 3)\n",
    "\n",
    "ax1.imshow(image, cmap='gray')\n",
    "ax1.set_axis_off()\n",
    "ax1.set_title('Reference image')\n",
    "\n",
    "ax2.imshow(offset_image.real, cmap='gray')\n",
    "ax2.set_axis_off()\n",
    "ax2.set_title('Offset image')\n",
    "\n",
    "# Calculate the upsampled DFT, again to show what the algorithm is doing\n",
    "# behind the scenes.  Constants correspond to calculated values in routine.\n",
    "# See source code for details.\n",
    "# cc_image = _upsampled_dft(image_product, 150, 100, (shift*100)+75).conj()\n",
    "# ax3.imshow(cc_image.real)\n",
    "# ax3.set_axis_off()\n",
    "# ax3.set_title(\"Supersampled XC sub-area\")\n",
    "\n",
    "\n",
    "plt.show()\n",
    "\n",
    "print(\"Detected subpixel offset (y, x): {}\".format(shift))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cfbb6b2",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb6f99f7",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0daa72f",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.12 ('echopiv')",
   "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.12"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 0.048653,
   "end_time": "2022-10-20T19:16:01.045215",
   "environment_variables": {},
   "exception": null,
   "input_path": "./notebooks/studying_correlations.ipynb",
   "output_path": "./notebooks/studying_correlations.ipynb",
   "parameters": {},
   "start_time": "2022-10-20T19:16:00.996562",
   "version": "2.4.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "f83b0c3a4910470a1212112b1707d582432916ed4ba8aec962241a050aa18fae"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}