2 - Driftcorrection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2 - Drift Correction\n",
    "This notebook implements and show cases the drift correction algorithm as described in section 4 of the paper. It uses cross correlation of all pairs of images after applying digital smoothing and edge detection filters to align Low Energy Electron Microscopy images with each other. When applied correctly, this allows for sub-pixel accurate image registration.\n",
    "This process is applied to the output of notebook 1, the detector corrected images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import numba\n",
    "import time\n",
    "import dask\n",
    "\n",
    "import dask.array as da\n",
    "from dask.delayed import delayed\n",
    "from dask.distributed import Client, LocalCluster\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interactive\n",
    "\n",
    "from scipy.optimize import least_squares\n",
    "import scipy.ndimage as ndi\n",
    "import scipy.sparse as sp\n",
    "from scipy.interpolate import interp1d\n",
    "\n",
    "from skimage import filters\n",
    "\n",
    "# Most relevant functions can be found in registration\n",
    "from registration import *\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = [12., 8.]\n",
    "SAVEFIG = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster = LocalCluster(n_workers=1, threads_per_worker=8)\n",
    "client = Client(cluster)\n",
    "client.upload_file('registration.py')\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def plot_stack(images, n, grid=False):\n",
    "    \"\"\"Plot the n-th image from a stack of n images.\n",
    "    For interactive use with ipython widgets\"\"\"\n",
    "    im = images[n, :, :].compute()\n",
    "    plt.figure(figsize=[12,10])\n",
    "    plt.imshow(im.T, cmap='gray', vmax=im.max())\n",
    "    if grid:\n",
    "        plt.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A bunch of constants\n",
    "\n",
    "folder = r'./data'\n",
    "\n",
    "# Pick the dataset to apply driftcorrection to and the range\n",
    "name = '20171120_160356_3.5um_591.4_IVhdr'\n",
    "start, stop = 42, -1 #BF\n",
    "\n",
    "# name = '20171120_215555_3.5um_583.1_IVhdr_DF2'\n",
    "# start, stop = 0, -1 #DF\n",
    "\n",
    "# A stride larger than 1 takes 1 every stride images of the total dataset. \n",
    "# This decreases computation time by a factor of stride**2, but decreases accuracy\n",
    "stride = 5\n",
    "# dE is the blocksize used by dask, the number of images computed for at once.\n",
    "dE = 10\n",
    "\n",
    "Eslice = slice(start, stop, stride)\n",
    "\n",
    "# Grab a window of 2*fftsize around the center of the picture\n",
    "fftsize=256 * 2 // 2\n",
    "z_factor = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = xr.open_dataset(os.path.join(folder, name + '_detectorcorrected.nc'), chunks={'time': dE})\n",
    "original = dataset.Intensity.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c442a250afa146a4b0a4adf5a5ddfd12",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=206, continuous_update=False, description='n', max=411), Output()), _dom…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# For interactive use we can view the original data\n",
    "interactive(lambda n: plot_stack(original, n), \n",
    "            n=widgets.IntSlider(original.shape[0]//2, 0, original.shape[0]-1, 1, continuous_update=False)\n",
    "           ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tr>\n",
       "<td>\n",
       "<table>\n",
       "  <thead>\n",
       "    <tr><td> </td><th> Array </th><th> Chunk </th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th> Bytes </th><td> 861.93 MB </td> <td> 20.97 MB </td></tr>\n",
       "    <tr><th> Shape </th><td> (411, 512, 512) </td> <td> (10, 512, 512) </td></tr>\n",
       "    <tr><th> Count </th><td> 337 Tasks </td><td> 42 Chunks </td></tr>\n",
       "    <tr><th> Type </th><td> float64 </td><td> numpy.ndarray </td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</td>\n",
       "<td>\n",
       "<svg width=\"236\" height=\"226\" style=\"stroke:rgb(0,0,0);stroke-width:1\" >\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"10\" y1=\"0\" x2=\"66\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"10\" y1=\"120\" x2=\"66\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"10\" y1=\"0\" x2=\"10\" y2=\"120\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"11\" y1=\"1\" x2=\"11\" y2=\"121\" />\n",
       "  <line x1=\"12\" y1=\"2\" x2=\"12\" y2=\"122\" />\n",
       "  <line x1=\"14\" y1=\"4\" x2=\"14\" y2=\"124\" />\n",
       "  <line x1=\"15\" y1=\"5\" x2=\"15\" y2=\"125\" />\n",
       "  <line x1=\"16\" y1=\"6\" x2=\"16\" y2=\"126\" />\n",
       "  <line x1=\"18\" y1=\"8\" x2=\"18\" y2=\"128\" />\n",
       "  <line x1=\"19\" y1=\"9\" x2=\"19\" y2=\"129\" />\n",
       "  <line x1=\"21\" y1=\"11\" x2=\"21\" y2=\"131\" />\n",
       "  <line x1=\"22\" y1=\"12\" x2=\"22\" y2=\"132\" />\n",
       "  <line x1=\"23\" y1=\"13\" x2=\"23\" y2=\"133\" />\n",
       "  <line x1=\"25\" y1=\"15\" x2=\"25\" y2=\"135\" />\n",
       "  <line x1=\"26\" y1=\"16\" x2=\"26\" y2=\"136\" />\n",
       "  <line x1=\"27\" y1=\"17\" x2=\"27\" y2=\"137\" />\n",
       "  <line x1=\"29\" y1=\"19\" x2=\"29\" y2=\"139\" />\n",
       "  <line x1=\"30\" y1=\"20\" x2=\"30\" y2=\"140\" />\n",
       "  <line x1=\"32\" y1=\"22\" x2=\"32\" y2=\"142\" />\n",
       "  <line x1=\"33\" y1=\"23\" x2=\"33\" y2=\"143\" />\n",
       "  <line x1=\"34\" y1=\"24\" x2=\"34\" y2=\"144\" />\n",
       "  <line x1=\"36\" y1=\"26\" x2=\"36\" y2=\"146\" />\n",
       "  <line x1=\"37\" y1=\"27\" x2=\"37\" y2=\"147\" />\n",
       "  <line x1=\"38\" y1=\"28\" x2=\"38\" y2=\"148\" />\n",
       "  <line x1=\"40\" y1=\"30\" x2=\"40\" y2=\"150\" />\n",
       "  <line x1=\"41\" y1=\"31\" x2=\"41\" y2=\"151\" />\n",
       "  <line x1=\"43\" y1=\"33\" x2=\"43\" y2=\"153\" />\n",
       "  <line x1=\"44\" y1=\"34\" x2=\"44\" y2=\"154\" />\n",
       "  <line x1=\"45\" y1=\"35\" x2=\"45\" y2=\"155\" />\n",
       "  <line x1=\"47\" y1=\"37\" x2=\"47\" y2=\"157\" />\n",
       "  <line x1=\"48\" y1=\"38\" x2=\"48\" y2=\"158\" />\n",
       "  <line x1=\"49\" y1=\"39\" x2=\"49\" y2=\"159\" />\n",
       "  <line x1=\"51\" y1=\"41\" x2=\"51\" y2=\"161\" />\n",
       "  <line x1=\"52\" y1=\"42\" x2=\"52\" y2=\"162\" />\n",
       "  <line x1=\"54\" y1=\"44\" x2=\"54\" y2=\"164\" />\n",
       "  <line x1=\"55\" y1=\"45\" x2=\"55\" y2=\"165\" />\n",
       "  <line x1=\"56\" y1=\"46\" x2=\"56\" y2=\"166\" />\n",
       "  <line x1=\"58\" y1=\"48\" x2=\"58\" y2=\"168\" />\n",
       "  <line x1=\"59\" y1=\"49\" x2=\"59\" y2=\"169\" />\n",
       "  <line x1=\"61\" y1=\"51\" x2=\"61\" y2=\"171\" />\n",
       "  <line x1=\"62\" y1=\"52\" x2=\"62\" y2=\"172\" />\n",
       "  <line x1=\"63\" y1=\"53\" x2=\"63\" y2=\"173\" />\n",
       "  <line x1=\"65\" y1=\"55\" x2=\"65\" y2=\"175\" />\n",
       "  <line x1=\"66\" y1=\"56\" x2=\"66\" y2=\"176\" />\n",
       "  <line x1=\"66\" y1=\"56\" x2=\"66\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"10.000000,0.000000 66.663603,56.663603 66.663603,176.663603 10.000000,120.000000\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"10\" y1=\"0\" x2=\"130\" y2=\"0\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"11\" y1=\"1\" x2=\"131\" y2=\"1\" />\n",
       "  <line x1=\"12\" y1=\"2\" x2=\"132\" y2=\"2\" />\n",
       "  <line x1=\"14\" y1=\"4\" x2=\"134\" y2=\"4\" />\n",
       "  <line x1=\"15\" y1=\"5\" x2=\"135\" y2=\"5\" />\n",
       "  <line x1=\"16\" y1=\"6\" x2=\"136\" y2=\"6\" />\n",
       "  <line x1=\"18\" y1=\"8\" x2=\"138\" y2=\"8\" />\n",
       "  <line x1=\"19\" y1=\"9\" x2=\"139\" y2=\"9\" />\n",
       "  <line x1=\"21\" y1=\"11\" x2=\"141\" y2=\"11\" />\n",
       "  <line x1=\"22\" y1=\"12\" x2=\"142\" y2=\"12\" />\n",
       "  <line x1=\"23\" y1=\"13\" x2=\"143\" y2=\"13\" />\n",
       "  <line x1=\"25\" y1=\"15\" x2=\"145\" y2=\"15\" />\n",
       "  <line x1=\"26\" y1=\"16\" x2=\"146\" y2=\"16\" />\n",
       "  <line x1=\"27\" y1=\"17\" x2=\"147\" y2=\"17\" />\n",
       "  <line x1=\"29\" y1=\"19\" x2=\"149\" y2=\"19\" />\n",
       "  <line x1=\"30\" y1=\"20\" x2=\"150\" y2=\"20\" />\n",
       "  <line x1=\"32\" y1=\"22\" x2=\"152\" y2=\"22\" />\n",
       "  <line x1=\"33\" y1=\"23\" x2=\"153\" y2=\"23\" />\n",
       "  <line x1=\"34\" y1=\"24\" x2=\"154\" y2=\"24\" />\n",
       "  <line x1=\"36\" y1=\"26\" x2=\"156\" y2=\"26\" />\n",
       "  <line x1=\"37\" y1=\"27\" x2=\"157\" y2=\"27\" />\n",
       "  <line x1=\"38\" y1=\"28\" x2=\"158\" y2=\"28\" />\n",
       "  <line x1=\"40\" y1=\"30\" x2=\"160\" y2=\"30\" />\n",
       "  <line x1=\"41\" y1=\"31\" x2=\"161\" y2=\"31\" />\n",
       "  <line x1=\"43\" y1=\"33\" x2=\"163\" y2=\"33\" />\n",
       "  <line x1=\"44\" y1=\"34\" x2=\"164\" y2=\"34\" />\n",
       "  <line x1=\"45\" y1=\"35\" x2=\"165\" y2=\"35\" />\n",
       "  <line x1=\"47\" y1=\"37\" x2=\"167\" y2=\"37\" />\n",
       "  <line x1=\"48\" y1=\"38\" x2=\"168\" y2=\"38\" />\n",
       "  <line x1=\"49\" y1=\"39\" x2=\"169\" y2=\"39\" />\n",
       "  <line x1=\"51\" y1=\"41\" x2=\"171\" y2=\"41\" />\n",
       "  <line x1=\"52\" y1=\"42\" x2=\"172\" y2=\"42\" />\n",
       "  <line x1=\"54\" y1=\"44\" x2=\"174\" y2=\"44\" />\n",
       "  <line x1=\"55\" y1=\"45\" x2=\"175\" y2=\"45\" />\n",
       "  <line x1=\"56\" y1=\"46\" x2=\"176\" y2=\"46\" />\n",
       "  <line x1=\"58\" y1=\"48\" x2=\"178\" y2=\"48\" />\n",
       "  <line x1=\"59\" y1=\"49\" x2=\"179\" y2=\"49\" />\n",
       "  <line x1=\"61\" y1=\"51\" x2=\"181\" y2=\"51\" />\n",
       "  <line x1=\"62\" y1=\"52\" x2=\"182\" y2=\"52\" />\n",
       "  <line x1=\"63\" y1=\"53\" x2=\"183\" y2=\"53\" />\n",
       "  <line x1=\"65\" y1=\"55\" x2=\"185\" y2=\"55\" />\n",
       "  <line x1=\"66\" y1=\"56\" x2=\"186\" y2=\"56\" />\n",
       "  <line x1=\"66\" y1=\"56\" x2=\"186\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"10\" y1=\"0\" x2=\"66\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"130\" y1=\"0\" x2=\"186\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"10.000000,0.000000 130.000000,0.000000 186.663603,56.663603 66.663603,56.663603\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"66\" y1=\"56\" x2=\"186\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"66\" y1=\"176\" x2=\"186\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"66\" y1=\"56\" x2=\"66\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"186\" y1=\"56\" x2=\"186\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"66.663603,56.663603 186.663603,56.663603 186.663603,176.663603 66.663603,176.663603\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Text -->\n",
       "  <text x=\"126.663603\" y=\"196.663603\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" >512</text>\n",
       "  <text x=\"206.663603\" y=\"116.663603\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" transform=\"rotate(-90,206.663603,116.663603)\">512</text>\n",
       "  <text x=\"28.331801\" y=\"168.331801\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" transform=\"rotate(45,28.331801,168.331801)\">411</text>\n",
       "</svg>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "dask.array<sub, shape=(411, 512, 512), dtype=float64, chunksize=(10, 512, 512), chunktype=numpy.ndarray>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Step 1 to 3 of the algorithm as described in section 4 of the paper.\n",
    "sobel = crop_and_filter(original[Eslice,...].rechunk({0:dE}), sigma=3, finalsize=2*fftsize)\n",
    "sobel = (sobel - sobel.mean(axis=(1,2), keepdims=True)) #.persist()  \n",
    "sobel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\dask\\array\\core.py:3855: PerformanceWarning: Increasing number of chunks by factor of 42\n",
      "  **blockwise_kwargs\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tr>\n",
       "<td>\n",
       "<table>\n",
       "  <thead>\n",
       "    <tr><td> </td><th> Array </th><th> Chunk </th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th> Bytes </th><td> 354.25 GB </td> <td> 209.72 MB </td></tr>\n",
       "    <tr><th> Shape </th><td> (411, 411, 512, 512) </td> <td> (10, 10, 512, 512) </td></tr>\n",
       "    <tr><th> Count </th><td> 23437 Tasks </td><td> 1764 Chunks </td></tr>\n",
       "    <tr><th> Type </th><td> float64 </td><td> numpy.ndarray </td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</td>\n",
       "<td>\n",
       "<svg width=\"558\" height=\"226\" style=\"stroke:rgb(0,0,0);stroke-width:1\" >\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"0\" y1=\"0\" x2=\"96\" y2=\"0\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"0\" y1=\"25\" x2=\"96\" y2=\"25\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"0\" y1=\"0\" x2=\"0\" y2=\"25\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"2\" y1=\"0\" x2=\"2\" y2=\"25\" />\n",
       "  <line x1=\"4\" y1=\"0\" x2=\"4\" y2=\"25\" />\n",
       "  <line x1=\"7\" y1=\"0\" x2=\"7\" y2=\"25\" />\n",
       "  <line x1=\"9\" y1=\"0\" x2=\"9\" y2=\"25\" />\n",
       "  <line x1=\"11\" y1=\"0\" x2=\"11\" y2=\"25\" />\n",
       "  <line x1=\"14\" y1=\"0\" x2=\"14\" y2=\"25\" />\n",
       "  <line x1=\"16\" y1=\"0\" x2=\"16\" y2=\"25\" />\n",
       "  <line x1=\"18\" y1=\"0\" x2=\"18\" y2=\"25\" />\n",
       "  <line x1=\"21\" y1=\"0\" x2=\"21\" y2=\"25\" />\n",
       "  <line x1=\"23\" y1=\"0\" x2=\"23\" y2=\"25\" />\n",
       "  <line x1=\"25\" y1=\"0\" x2=\"25\" y2=\"25\" />\n",
       "  <line x1=\"28\" y1=\"0\" x2=\"28\" y2=\"25\" />\n",
       "  <line x1=\"30\" y1=\"0\" x2=\"30\" y2=\"25\" />\n",
       "  <line x1=\"32\" y1=\"0\" x2=\"32\" y2=\"25\" />\n",
       "  <line x1=\"35\" y1=\"0\" x2=\"35\" y2=\"25\" />\n",
       "  <line x1=\"37\" y1=\"0\" x2=\"37\" y2=\"25\" />\n",
       "  <line x1=\"39\" y1=\"0\" x2=\"39\" y2=\"25\" />\n",
       "  <line x1=\"42\" y1=\"0\" x2=\"42\" y2=\"25\" />\n",
       "  <line x1=\"44\" y1=\"0\" x2=\"44\" y2=\"25\" />\n",
       "  <line x1=\"46\" y1=\"0\" x2=\"46\" y2=\"25\" />\n",
       "  <line x1=\"49\" y1=\"0\" x2=\"49\" y2=\"25\" />\n",
       "  <line x1=\"51\" y1=\"0\" x2=\"51\" y2=\"25\" />\n",
       "  <line x1=\"53\" y1=\"0\" x2=\"53\" y2=\"25\" />\n",
       "  <line x1=\"56\" y1=\"0\" x2=\"56\" y2=\"25\" />\n",
       "  <line x1=\"58\" y1=\"0\" x2=\"58\" y2=\"25\" />\n",
       "  <line x1=\"60\" y1=\"0\" x2=\"60\" y2=\"25\" />\n",
       "  <line x1=\"63\" y1=\"0\" x2=\"63\" y2=\"25\" />\n",
       "  <line x1=\"65\" y1=\"0\" x2=\"65\" y2=\"25\" />\n",
       "  <line x1=\"67\" y1=\"0\" x2=\"67\" y2=\"25\" />\n",
       "  <line x1=\"70\" y1=\"0\" x2=\"70\" y2=\"25\" />\n",
       "  <line x1=\"72\" y1=\"0\" x2=\"72\" y2=\"25\" />\n",
       "  <line x1=\"75\" y1=\"0\" x2=\"75\" y2=\"25\" />\n",
       "  <line x1=\"77\" y1=\"0\" x2=\"77\" y2=\"25\" />\n",
       "  <line x1=\"79\" y1=\"0\" x2=\"79\" y2=\"25\" />\n",
       "  <line x1=\"82\" y1=\"0\" x2=\"82\" y2=\"25\" />\n",
       "  <line x1=\"84\" y1=\"0\" x2=\"84\" y2=\"25\" />\n",
       "  <line x1=\"86\" y1=\"0\" x2=\"86\" y2=\"25\" />\n",
       "  <line x1=\"89\" y1=\"0\" x2=\"89\" y2=\"25\" />\n",
       "  <line x1=\"91\" y1=\"0\" x2=\"91\" y2=\"25\" />\n",
       "  <line x1=\"93\" y1=\"0\" x2=\"93\" y2=\"25\" />\n",
       "  <line x1=\"96\" y1=\"0\" x2=\"96\" y2=\"25\" />\n",
       "  <line x1=\"96\" y1=\"0\" x2=\"96\" y2=\"25\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"0.000000,0.000000 96.328125,0.000000 96.328125,25.412617 0.000000,25.412617\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Text -->\n",
       "  <text x=\"48.164062\" y=\"45.412617\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" >411</text>\n",
       "  <text x=\"116.328125\" y=\"12.706308\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" transform=\"rotate(0,116.328125,12.706308)\">1</text>\n",
       "\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"166\" y1=\"0\" x2=\"222\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"166\" y1=\"120\" x2=\"222\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"166\" y1=\"0\" x2=\"166\" y2=\"120\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"167\" y1=\"1\" x2=\"167\" y2=\"121\" />\n",
       "  <line x1=\"168\" y1=\"2\" x2=\"168\" y2=\"122\" />\n",
       "  <line x1=\"170\" y1=\"4\" x2=\"170\" y2=\"124\" />\n",
       "  <line x1=\"171\" y1=\"5\" x2=\"171\" y2=\"125\" />\n",
       "  <line x1=\"172\" y1=\"6\" x2=\"172\" y2=\"126\" />\n",
       "  <line x1=\"174\" y1=\"8\" x2=\"174\" y2=\"128\" />\n",
       "  <line x1=\"175\" y1=\"9\" x2=\"175\" y2=\"129\" />\n",
       "  <line x1=\"177\" y1=\"11\" x2=\"177\" y2=\"131\" />\n",
       "  <line x1=\"178\" y1=\"12\" x2=\"178\" y2=\"132\" />\n",
       "  <line x1=\"179\" y1=\"13\" x2=\"179\" y2=\"133\" />\n",
       "  <line x1=\"181\" y1=\"15\" x2=\"181\" y2=\"135\" />\n",
       "  <line x1=\"182\" y1=\"16\" x2=\"182\" y2=\"136\" />\n",
       "  <line x1=\"183\" y1=\"17\" x2=\"183\" y2=\"137\" />\n",
       "  <line x1=\"185\" y1=\"19\" x2=\"185\" y2=\"139\" />\n",
       "  <line x1=\"186\" y1=\"20\" x2=\"186\" y2=\"140\" />\n",
       "  <line x1=\"188\" y1=\"22\" x2=\"188\" y2=\"142\" />\n",
       "  <line x1=\"189\" y1=\"23\" x2=\"189\" y2=\"143\" />\n",
       "  <line x1=\"190\" y1=\"24\" x2=\"190\" y2=\"144\" />\n",
       "  <line x1=\"192\" y1=\"26\" x2=\"192\" y2=\"146\" />\n",
       "  <line x1=\"193\" y1=\"27\" x2=\"193\" y2=\"147\" />\n",
       "  <line x1=\"194\" y1=\"28\" x2=\"194\" y2=\"148\" />\n",
       "  <line x1=\"196\" y1=\"30\" x2=\"196\" y2=\"150\" />\n",
       "  <line x1=\"197\" y1=\"31\" x2=\"197\" y2=\"151\" />\n",
       "  <line x1=\"199\" y1=\"33\" x2=\"199\" y2=\"153\" />\n",
       "  <line x1=\"200\" y1=\"34\" x2=\"200\" y2=\"154\" />\n",
       "  <line x1=\"201\" y1=\"35\" x2=\"201\" y2=\"155\" />\n",
       "  <line x1=\"203\" y1=\"37\" x2=\"203\" y2=\"157\" />\n",
       "  <line x1=\"204\" y1=\"38\" x2=\"204\" y2=\"158\" />\n",
       "  <line x1=\"205\" y1=\"39\" x2=\"205\" y2=\"159\" />\n",
       "  <line x1=\"207\" y1=\"41\" x2=\"207\" y2=\"161\" />\n",
       "  <line x1=\"208\" y1=\"42\" x2=\"208\" y2=\"162\" />\n",
       "  <line x1=\"210\" y1=\"44\" x2=\"210\" y2=\"164\" />\n",
       "  <line x1=\"211\" y1=\"45\" x2=\"211\" y2=\"165\" />\n",
       "  <line x1=\"212\" y1=\"46\" x2=\"212\" y2=\"166\" />\n",
       "  <line x1=\"214\" y1=\"48\" x2=\"214\" y2=\"168\" />\n",
       "  <line x1=\"215\" y1=\"49\" x2=\"215\" y2=\"169\" />\n",
       "  <line x1=\"217\" y1=\"51\" x2=\"217\" y2=\"171\" />\n",
       "  <line x1=\"218\" y1=\"52\" x2=\"218\" y2=\"172\" />\n",
       "  <line x1=\"219\" y1=\"53\" x2=\"219\" y2=\"173\" />\n",
       "  <line x1=\"221\" y1=\"55\" x2=\"221\" y2=\"175\" />\n",
       "  <line x1=\"222\" y1=\"56\" x2=\"222\" y2=\"176\" />\n",
       "  <line x1=\"222\" y1=\"56\" x2=\"222\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"166.000000,0.000000 222.663603,56.663603 222.663603,176.663603 166.000000,120.000000\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"166\" y1=\"0\" x2=\"286\" y2=\"0\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"167\" y1=\"1\" x2=\"287\" y2=\"1\" />\n",
       "  <line x1=\"168\" y1=\"2\" x2=\"288\" y2=\"2\" />\n",
       "  <line x1=\"170\" y1=\"4\" x2=\"290\" y2=\"4\" />\n",
       "  <line x1=\"171\" y1=\"5\" x2=\"291\" y2=\"5\" />\n",
       "  <line x1=\"172\" y1=\"6\" x2=\"292\" y2=\"6\" />\n",
       "  <line x1=\"174\" y1=\"8\" x2=\"294\" y2=\"8\" />\n",
       "  <line x1=\"175\" y1=\"9\" x2=\"295\" y2=\"9\" />\n",
       "  <line x1=\"177\" y1=\"11\" x2=\"297\" y2=\"11\" />\n",
       "  <line x1=\"178\" y1=\"12\" x2=\"298\" y2=\"12\" />\n",
       "  <line x1=\"179\" y1=\"13\" x2=\"299\" y2=\"13\" />\n",
       "  <line x1=\"181\" y1=\"15\" x2=\"301\" y2=\"15\" />\n",
       "  <line x1=\"182\" y1=\"16\" x2=\"302\" y2=\"16\" />\n",
       "  <line x1=\"183\" y1=\"17\" x2=\"303\" y2=\"17\" />\n",
       "  <line x1=\"185\" y1=\"19\" x2=\"305\" y2=\"19\" />\n",
       "  <line x1=\"186\" y1=\"20\" x2=\"306\" y2=\"20\" />\n",
       "  <line x1=\"188\" y1=\"22\" x2=\"308\" y2=\"22\" />\n",
       "  <line x1=\"189\" y1=\"23\" x2=\"309\" y2=\"23\" />\n",
       "  <line x1=\"190\" y1=\"24\" x2=\"310\" y2=\"24\" />\n",
       "  <line x1=\"192\" y1=\"26\" x2=\"312\" y2=\"26\" />\n",
       "  <line x1=\"193\" y1=\"27\" x2=\"313\" y2=\"27\" />\n",
       "  <line x1=\"194\" y1=\"28\" x2=\"314\" y2=\"28\" />\n",
       "  <line x1=\"196\" y1=\"30\" x2=\"316\" y2=\"30\" />\n",
       "  <line x1=\"197\" y1=\"31\" x2=\"317\" y2=\"31\" />\n",
       "  <line x1=\"199\" y1=\"33\" x2=\"319\" y2=\"33\" />\n",
       "  <line x1=\"200\" y1=\"34\" x2=\"320\" y2=\"34\" />\n",
       "  <line x1=\"201\" y1=\"35\" x2=\"321\" y2=\"35\" />\n",
       "  <line x1=\"203\" y1=\"37\" x2=\"323\" y2=\"37\" />\n",
       "  <line x1=\"204\" y1=\"38\" x2=\"324\" y2=\"38\" />\n",
       "  <line x1=\"205\" y1=\"39\" x2=\"325\" y2=\"39\" />\n",
       "  <line x1=\"207\" y1=\"41\" x2=\"327\" y2=\"41\" />\n",
       "  <line x1=\"208\" y1=\"42\" x2=\"328\" y2=\"42\" />\n",
       "  <line x1=\"210\" y1=\"44\" x2=\"330\" y2=\"44\" />\n",
       "  <line x1=\"211\" y1=\"45\" x2=\"331\" y2=\"45\" />\n",
       "  <line x1=\"212\" y1=\"46\" x2=\"332\" y2=\"46\" />\n",
       "  <line x1=\"214\" y1=\"48\" x2=\"334\" y2=\"48\" />\n",
       "  <line x1=\"215\" y1=\"49\" x2=\"335\" y2=\"49\" />\n",
       "  <line x1=\"217\" y1=\"51\" x2=\"337\" y2=\"51\" />\n",
       "  <line x1=\"218\" y1=\"52\" x2=\"338\" y2=\"52\" />\n",
       "  <line x1=\"219\" y1=\"53\" x2=\"339\" y2=\"53\" />\n",
       "  <line x1=\"221\" y1=\"55\" x2=\"341\" y2=\"55\" />\n",
       "  <line x1=\"222\" y1=\"56\" x2=\"342\" y2=\"56\" />\n",
       "  <line x1=\"222\" y1=\"56\" x2=\"342\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"166\" y1=\"0\" x2=\"222\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"286\" y1=\"0\" x2=\"342\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"166.000000,0.000000 286.000000,0.000000 342.663603,56.663603 222.663603,56.663603\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"222\" y1=\"56\" x2=\"342\" y2=\"56\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"222\" y1=\"176\" x2=\"342\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"222\" y1=\"56\" x2=\"222\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"342\" y1=\"56\" x2=\"342\" y2=\"176\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"222.663603,56.663603 342.663603,56.663603 342.663603,176.663603 222.663603,176.663603\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Text -->\n",
       "  <text x=\"282.663603\" y=\"196.663603\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" >512</text>\n",
       "  <text x=\"362.663603\" y=\"116.663603\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" transform=\"rotate(-90,362.663603,116.663603)\">512</text>\n",
       "  <text x=\"184.331801\" y=\"168.331801\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" transform=\"rotate(45,184.331801,168.331801)\">411</text>\n",
       "</svg>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "dask.array<rechunk-merge, shape=(411, 411, 512, 512), dtype=float64, chunksize=(10, 10, 512, 512), chunktype=numpy.ndarray>"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Step 4 of the algorithm as described in paper.\n",
    "Corr = dask_cross_corr(sobel)\n",
    "Corr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cb39b27f48714a2896fb83f978903554",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=58, continuous_update=False, description='i', max=410), IntSlider(value=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot combination of original images, filtered images, crosscorrelation \n",
    "# for illustration purposes\n",
    "def plot_corr(i,j, save=SAVEFIG):\n",
    "    #fig = plt.figure(figsize=(8.2, 3.5), constrained_layout=True)\n",
    "    fig = plt.figure(figsize=(4, 7), constrained_layout=True)\n",
    "    fig.set_constrained_layout_pads(hspace=0.0, wspace=0.06)\n",
    "    #gs = mpl.gridspec.GridSpec(2, 3,\n",
    "    #                   width_ratios=[1, 1, 2.9],\n",
    "    #                   #height_ratios=[4, 1]\n",
    "    #                   )\n",
    "    \n",
    "    gs = mpl.gridspec.GridSpec(3, 2,\n",
    "                       height_ratios=[1, 1, 1.8],\n",
    "                               figure=fig,\n",
    "                       )\n",
    "    \n",
    "    ax0 = plt.subplot(gs[0, 0])\n",
    "    ax1 = plt.subplot(gs[1, 0])\n",
    "    ax2 = plt.subplot(gs[0, 1])\n",
    "    ax3 = plt.subplot(gs[1, 1])\n",
    "    ax4 = plt.subplot(gs[2, :]) #2grid((2, 4), (0, 2), rowspan=2, colspan=2)\n",
    "    ax0.imshow(original[i*stride + start,(640-fftsize):(640+fftsize),(512-fftsize):(512+fftsize)].T, \n",
    "               cmap='gray', interpolation='none')\n",
    "    ax0.set_title(f'i={i*stride + start}')\n",
    "    ax1.imshow(sobel[i,...].T, cmap='gray')\n",
    "    ax2.imshow(original[j*stride + start,(640-fftsize):(640+fftsize),(512-fftsize):(512+fftsize)].T,\n",
    "               cmap='gray', interpolation='none')\n",
    "    ax2.set_title(f'j={j*stride + start}')\n",
    "    ax3.imshow(sobel[j,...].T,\n",
    "               cmap='gray', interpolation='none')\n",
    "    im = ax4.imshow(Corr[i,j,...].compute().T, extent=[-fftsize, fftsize, -fftsize, fftsize], interpolation='none')\n",
    "    ax4.axhline(0, color='white', alpha=0.5)\n",
    "    ax4.axvline(0, color='white', alpha=0.5)\n",
    "    for ax in [ax2, ax3]:\n",
    "        ax.yaxis.set_label_position(\"right\")\n",
    "        ax.tick_params(axis='y', labelright=True, labelleft=False)\n",
    "    plt.colorbar(im, ax=ax4)\n",
    "    if save:\n",
    "        #Saving Figure for paper.\n",
    "        plt.savefig('autocorrelation.pdf', dpi=300)\n",
    "    plt.show()\n",
    "    return fig\n",
    "\n",
    "widget = interactive(plot_corr, \n",
    "            i=widgets.IntSlider(58-start,0,sobel.shape[0]-1,1, continuous_update=False), \n",
    "            j=widgets.IntSlider(100-start,0,sobel.shape[0]-1,1, continuous_update=False),\n",
    "            save=SAVEFIG\n",
    "           )\n",
    "display(widget)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 5 of the algorithm\n",
    "weights, argmax = max_and_argmax(Corr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "806.594000000041\n"
     ]
    }
   ],
   "source": [
    "# Do actual computations; get a cup of coffee. If this takes to long, consider increasing stride to reduce the workload, at the cost of accuracy\n",
    "t = time.monotonic()\n",
    "W, DX_DY = calculate_halfmatrices(weights, argmax, fftsize=fftsize)\n",
    "print(time.monotonic()-t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "# Step 6 of the algorithm\n",
    "w_diag = np.atleast_2d(np.diag(W))\n",
    "W_n = W / np.sqrt(w_diag.T*w_diag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eec7f0be4cca4fe09b53950ac962240d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.15, continuous_update=False, description='min_normed_weight', max=1.…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot W, DX and DY to pick a value for W_{min} (Step 7 of algorithm)\n",
    "def plot_masking(min_normed_weight, save=SAVEFIG):\n",
    "    extent = [start, stop, stop, start]\n",
    "    fig, axs = plt.subplots(1, 3, figsize=(8, 2.5), constrained_layout=True)\n",
    "    im = {}\n",
    "    im[0] = axs[0].imshow(DX_DY[0], cmap='seismic', extent=extent, interpolation='none')\n",
    "    im[1] = axs[1].imshow(DX_DY[1], cmap='seismic', extent=extent, interpolation='none')\n",
    "    im[2] = axs[2].imshow(W_n - np.diag(np.diag(W_n)), cmap='inferno', \n",
    "                          extent=extent, clim=(0.0, None), interpolation='none')\n",
    "    axs[0].set_ylabel('$j$')\n",
    "    fig.colorbar(im[0], ax=axs[:2], shrink=0.82, fraction=0.1)\n",
    "    axs[0].contourf(W_n, [0, min_normed_weight], \n",
    "                    colors='black', alpha=0.6, \n",
    "                    extent=extent, origin='upper')\n",
    "    axs[1].contourf(W_n, [0, min_normed_weight], \n",
    "                    colors='black', alpha=0.6, \n",
    "                    extent=extent, origin='upper')\n",
    "    CF = axs[2].contourf(W_n, [0, min_normed_weight], \n",
    "                         colors='white', alpha=0.2, \n",
    "                         extent=extent, origin='upper')\n",
    "    cbar = fig.colorbar(im[2], ax=axs[2], shrink=0.82, fraction=0.1)\n",
    "    cbar.ax.fill_between([0,1], 0, min_normed_weight, color='white', alpha=0.2)\n",
    "    for i in range(3):            \n",
    "        axs[i].set_xlabel('$i$')\n",
    "        axs[i].tick_params(labelbottom=False, labelleft=False)  \n",
    "    axs[0].set_title('$DX_{ij}$')\n",
    "    axs[1].set_title('$DY_{ij}$')\n",
    "    axs[2].set_title('$W_{ij}$')\n",
    "    if save:\n",
    "        plt.savefig('shiftsandweights.pdf', dpi=300)\n",
    "    plt.show()\n",
    "    return min_normed_weight\n",
    "\n",
    "widget = interactive(plot_masking, \n",
    "                     min_normed_weight=widgets.FloatSlider(value=0.15, min=0., max=1, \n",
    "                                                      step=0.01, continuous_update=False),\n",
    "                     save=SAVEFIG\n",
    "               )\n",
    "display(widget)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2ad8e438438>"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAAHNCAYAAAD42XiHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzde3Rc9Xn/+893z1iO3TiyMGAulm0MxikmcWKpQhSc4MRNSGEghXAxpAnJMUpXyWr4Ne1ZSRq8WP41l985sEp7fpw0OqwE2thOCiYh47Y0oUBqUoQiGQh2EsBRLEs4gFEGxcTG0sz+nj/2fGf2jEbSSBppbu/XWl7S7Ll9rRl7Hj37+T6Psdaq3E4++WTb1NSkxsbGOX/uI0eO6JRTTpnT5xweHi7L37Wenrdcf9dyvJ8kXttafl7+j6rN562n95NUXz/jenveqvs/ylpb9j8tLS325ptvtuXQ0tIy589Zrr9rPT1vPb2frOW1reXn5f+o2nzeeno/WVtfP+N6e95q+z/KK20MPn2xWKzcS5gz5fq71tPz1tP7SeK1reXnLYd6+xnX0/u4XOrpZ1xvz1sOM/m7GlsBpQ6rV6+2GzduVCwWm/MXrrW1VT09PXP6nKhdvJ9QarynUEq8n1Bq1fSeisfjitx+++3lXof+6Z/+6fZ4PK41a9aU5flbWlrK8ryoTbyfUGq8p1BKvJ9QatXynlqzZk1lZHxbW1tttfy2AAAAgOpUMTW+AAAAqD+9/Qnd/dgB9fYnpnX/66+/vujbEvgCAACgLHr7E7rxni7d+YPndeM9XUUHv88//7yuueYa3Xbbbfrd736nzZs3a3h4WNdee61GRkbGvV+0VAsHAAAApqKrb0gjSV++lUaTvrr6htSyomnS+3V2duqOO+7QGWecocsuu0x/9Vd/pQsvvFAPPPCAGhoaxr0fGV8AAACURfuqJWqIeooYaV7UU/uqJUXdz1qr+fPnKxqNKhqN6uDBg1q6dKneeOONCe9HxhcAAABl0bKiSdu3tKurb0jtq5YUle2VpJtvvlnbtm3TWWedJUl64okn9IMf/EAf/ehH9Y1vfEO/93u/V/B+dHUAAABAXaDUAQAAAHWBwBcAAAB1gcAXAAAAdYHAFwAAAHWBwBcAAADlM9At7bkz+DoNVTe5bXh4WB0dHYrH4+VeCgAAQN2a6fjgKRvolu67Qnr0S8HXIoPf8OS24eFh3XHHHZKkT3/603rllVcK3icej1dGH9/GxkZ1dnaWexkAAAB1y40PHkn6aoh62r6lvei+utN2cI+UGpFsKvh6cI/U3Dbp3fIntz377LMaGhrSyMiIli5dWvA+sVisMjK+AABg9sx5Fg9Tlv8aleM1KzQ+eNat3CBFGiQTCb6u3FDU3fInt8ViMd1www36xCc+MeH9KiLjCwAApqa3P5GZdiVp3MlXZcniISP8Oo33c89/jbZevlbbdu+f89fMjQ8eTfpTGh88I81t0se/H2R6V24oKtsrjZ3c9uEPf1h///d/rwsvvHDC+zG5DQCACuOCpaaFDUocG8kETeHjLjDyjGRlZK1V1DO6prVZV61flgmU7n7sgO78wfPyrRQx0nVty3Xm4gVjHjM/MJvq8cmuq0eFAtp9h4dlpHFfI0/S8iULdeg3xzKXL1p9sm7ddO6c/EwLvfek8X+xqgTHjh3TZz7zGV1xxRWKxWIT3pbAFwCAMsvP3t54T5dOjPqykoykiGe05eKzdO+TB9PBrpFvrfwCH+FG0vx52SyhC75Gk74inpGMyXx/+TtP1+6f/lq+tTmZxUIBW+LYSE7AnR/IrT2jsSxZykrV25/QXY+8oB8feC3zOnlGme+jEaP3rTlVpyyan/3ZjfryCzyWZzTmZzqbv2SEX/9o3ntm25Xn64YLlpf0+eYSgS8AAGWUH2RevX6ZdnYfGhPUGknukCfJ84xSvlWhT3EjafMFy/XlP3lH5jm6+oZ0+PXj2vHUoXHvc3E6s9jVN5STgXTPZYxkbWgdoUAu4gVZZ5dZ/ssPrNEtG8+Z2Q+nSu146pC2PrRv3NcnnwuCX/ntm/rp4HDmPib91SqbrZek146e0OMvHClZMJofRIcz0OH3nSRFPaPvfOrCqv2lhhpfAADKKH9DkVWQ3cvP/uUEH1FPt8eCbOsDvYNKJn0ZI8lIKT+47b/0DOi1oyd0yqL5umr9MrWvWqJdewczQWw+K+mJF1/Tfx94TS0rmhSNeEomfclIyfTt83Nl4YdJ+UGphZFVxDM6/Ppx9fYnxpRTuL9zsZnKSi+fCJcG7Ds8rNeOntB//uLVgj/j8SRTVj/82SuaFzFjXp+IZ+SnL3/nJ4eUyksJJ32r2773nPYfHs4pnyhmzeEzDOFMvav1HRkN3ldW2dfat1ZdfUMV+VoUg4wvAABlFC5FmJdXPrBoflT3PPGrnMyhkbTpvKV6V/PiMYHkg3sHC2Z0oxEjzxglU9maYD/9mEaSMRqTYY54wbP5U8xaStLjLxxRMhUEUjdduFL3PPEr+ekaZKXXUUw5RLEb8/LrUvNro/NvV6q6Vbc+V5ZSSMRIkYinVCrIzl6S/hk9+vyrSqbsmNu+7/eX6rFfvJopP3E/v8myx/klLpOtudAZhnA98fMvH9XWh/bJtzanjrzay1jI+AIAUGLF7uTftXdQRsoEu68dPaHb4/tzAsM/Wnuadu0d1AO9g5ng6UcvHNF//vyVzG3CJQW79g6OCcSSKSujbKB7XVuzzly8ICdQ3PrQvkxmVwoyx+H7mPR5dxe85gdyj79wRI/8/JWc+uM3R319fU9fJlM8krJyueuRUV93PfJCwU1b4dKMcDZ8197BMYHr0eOjOYGhOzWfXxfryg+mE4CPx2XrxwtIo+kyhDWnLRrzfnCvvytbSKWCX3z+7L1n68/ee3bm9l19Q/Lt2KA3/xcTq8I/I/ceW3tGoxLHRnJ+pidGfb169ETOGYYfH3hNPzn4G129flnmdQy/Zyo1814sAl8AAEoonAWM5AU+4UDz9u/vSweC2YzsaCiIcn1Ub9l4jlpWNOnq9csywaDL0LnbuECkZUWTtm9p1669g/qXnoFMRtE9vguurh7nlPht33tOLgmZfx+3wa1QtvTuxw7okZ+/EmSNrc0EydLY8ggnHGSFHzs/02iMkafgMb/zkwHZUAayUDbaXfZtNrj+0Pmn5wT2o+kA3AWL0z1137SwQZ4JMqE5pSgRo+vyumvkP37LiqYJN6qFb+9ajLlfNFz5iqTML0XJZFADcX/6dfdCpS+OkdI/0+AHZSU99vyrmfri514aHlNy485EjPeeqTaUOgAAUALjbSBzp7pdMOEZBQFlXtAW3kQ00Wnr/NKIiU7/u2yfC5KmkoUu9j6F1vSe1afohz97Jefv5+VtjAtzdaxeXj1p5vpxjhejUCmHJ8l42VP34cC72JZuLuB0AemWi8/Sb08kx7QqK4XJziBMdTNd/oY1I2leJJvJd+8rqbLbmE0HgS8AoK6Vou4z/zR6yiqzQSm8Mz8sv7WVFyofyO/FO9GaKyUgKbRhygWF17Q2Z1p2uWO+NKbGdTz5gVqh68Pt2VwHikKdMTwTdKlIpuyYlm7hPsju75BfX1woyCx3F4twF4ZiuJ9X/t8hv8fzeIrp8SxVZtBMqQMAoC4U+rAOb/SZ7gae3v5Ezmn0pG/1/tAGpfFOy3ue0aa809bFBgrh0+SVIn9N27e0j/n7hGtdx9uI50TSAarv20z/Yde9wvOMUikrL51pXbRgXqge9nBQ62tzf7lw7dqWn7RQO7uD5/V9q+8/ezhzm5GU1fanDuk7PxnQ+uWLM7XS45VMuMedsyln4whPXDN5A01caYTbKOmHstyuK8hkJTBh+b/kuV8UwiUqpaqhng0EvgCAmpffkN99WIdbiQWxTHrjVZF1n25IQX7rqo1rTtWfvffsMafD1y9frJ8cTASbkXyrdc2Lc7KElRIclEKh4Dz/2K69gwWDtXDWtVAGMfx9/nOE61LDHSUaop5u3XRu3vMWbu2W9K26DyYyl30Frd7CwyikIDi/vm15yUsbpsrVdk/28/mjtaeNOe5qx4v5hSv/l7yRlNWOpw7p/p6BnLMcE9VQl/tsRUWUOqxevdpu3LhRsVhs0lFzAIDak19bWuoPxPxTwa6Gduvla7Vt9/4xXRDcbvz81lhSNqBwGa5CtbquPjc8CMKdSn4wHXRNVJ9bL2bj1Hh+YDVepj9n9HOBfrUTcaUC1T7FbKrufuyA7viP58dk6fNLUVxdezLpywtt8Mz/RXCuf37xeLwyAl9qfAGgssxGVmaiusDNnU9mOhw0RD3tvHl6AeF4zzFeXWb+JCz3Qe3GA7uA2DPKOX1bqI53xZKFOvSbYzmTy9pXLRmzEU2qzNrHehV+z+T/MuNe9/zX2k24q7fXL3PmZDS35KRQ+ckvX/udHv3FqxN24SjHFDhKHQAAOYodGjDZY+Rn3cKPGR7ScOToiUzQKxXfXqrY58j0SU13DQjXjT7QO5ipQ7w9lt3Z/2BeL1zf5vWgzR884Bl1vOfszOYtV/OZfwp6vNZWKJ9w+UXLiqYxrefyW6y5kol6fA3z38/hn4sbk+02Meb/+ym0PbEcU+AIfAGgxhQ7PCH/w90Fj3c98kLO0IBi6/PGnD4OBaD/vu/Xmcd8c9TXF7/33LinlIvZKFQoOA/X646M+jk9acOuaQ0a8ef3w00cG9EtG89Rb39C9/cMFN0WKnzKttCggkrciIbxjVebXOi1rUfhn48bruGnM+JnLF6gxLGRCYd65LeSm+tNgQS+AFBDisnW5p/2dxOu8kejeiY3CJ0oaxveMR6e3DUy6hesgx0v6F23rFFbY2sn7aWaH+Rui+/X0re9RdH0aNhCm5bc7vurQ43/3e3Df8+uvqHM5h0j6exT36pfvvpGzlQwKXe8K0Fu7eO1HSvcTSL8byh/4Ea4jGi8SXZzhcAXAGpIfkAYHgnrNpB95ycDOUGhu23nnr6czV8XnZMb1BWbUXUtvDwTPE44iBw3C5QOvrfG1koKNtGE+8GGg+3EsREdPT6amZjlS3p2cFjSsKIRo+vbgtOt+ZPRrivQnzXqmTG78vM/zD950Vk5/WfDTf7r9ZQ3II0tfXD/FvKPTTaZbi6xuQ0Aakh484mvUEA5TvcCqfBkq/yuBjm1sqGM6nifIG5QQHhkav6QBtdfdO0ZjTldE8KB7tXrl2XKEdw63fSvQut2dYZnLl6QWXd+p4hwh4fxBg9M1BlAYnMaUK3I+AJADXEZmLseeSHTc/TEqK9v/PhXBVt2hSddSdma1S0Xn1W4zVfE6P2/v1SS9Ojzr+ZM3soPTP3sXjAZSde2Nk/aM/Tuxw5kssonRn298MrRTGbXKnhsxw0piBhlMs+RSO6GtUKlHuOdns3/OYbvV+gygOpD4AsANaZlRZNu3XSunuob0kgqCBgPvPpG5vrwaf/wpKtwzaora8jP6CZTVo+mJ5JFPaMPnLc0J2vrNrYVKgtwU6EmChrbVy1R1DOZdf/kYCInoA4LZ7NdZldSzoa1QjvGxzs9C6D2EfgCQA0Jn5K/prV5zEhYl3n90p+8I3MsnP0M16w2RL1MyYQT8UzOLu78yWNS7lhaaWplAS0rmsas22V2o+k2ZPljavM3wu0KDYgYb8c4G5WA+kSNLwDUiEJdF9xUqnC9b/7p/4kGS7j2ZC6j6np0zubksYnqlMOt1ya6P9lcAIUQ+AJAjfib7z6XyZSGJ4cV6tc7E3MRWI7XZxgAZmJGpQ7GmBZJd0nyJb0i6UZJP5P0UvomX7LW/tAYc56kr6ef7zZr7SMzeV4AQK78oQuRSHZyWKkDxrkoE6AUAcBsmGmN70uSPmitPWaM+bKkD0sattZekne7L0v6pILg+GFJBL4AUEL5Qxc+0rKMwBEA8ngzubO19mVr7bH0xVFJSUlvNcb8yBizwxhzUvq60621L1prfytpyBhz8kyeFwCQ1duf0EuvH1c04ilipPnzstPJAABZMwp8HWPMckmbJO2WdJG19r0KMru3u5uEbj4s6aTw/Y8cOaLW1tbMn87OzlIsCwBqntsI9u3uQ5K1ur5t+axsOAOAatbZ2anW1taZtzMzxrxN0j9L+oS1dlTSUPqq+yVtSX8f7oazWNJvwo9xyimniM1tADB14THCKd/qjMULCHoBIE9HR4c6OjpmlvE1xkQkbZe0zVr7gjGmwRgzP331eyQdSH//sjFmdTpIPsla+9pMnhcAEHBTyCJGE/atBQDMfHPbtZL+UNIiY8xtkr4m6f80xvxO0gkFG9ok6QuSvpF+vq0zfE4AqFqlbgXGFDIAKB59fAFgjuQPmCimFrdQoMyABgCYHkYWA8AceXDvoE6M+rKSRpO+du0dnDCALRQoSxoznY0BDwBQHDK+AFBi4YysJO3aO6jXjp7Qo8+/qmQq+D834kmSkbU2J/sbvu+Dewczk9g8Se9Y1ihJeu6lYfk2OOZ5Rn7eYwAACiPjCwAlFM7SRj0jX8oEu2G+lVziYSSd/d21d1AP9A4qmcre193Tl/Ts4HDm/p6RPBMEvb4NMshdfUMEvgAwAQJfABjHVGppe/sT2rV3UPtfGs60FxtNWeWHvEZSxDNK+aFrrPSdnwzI97O3HykQLIe948xGLX3bW/T4C0eUSvl0dACAIhD4AoDGBrk7njqkrQ/tK6qMoLc/oc2dT+YEq0ZBVlZGSqU7mUcjRte1NmvtGY3atnu/RkZ9GRNkdXMC4RDPaEzmOBox+vnLR/XcS8OKekbXty3XVesZUQwAkyHwBVD38ssTLllzqv7zF69mgtGR0Ea0poUN2nd4WEbS2jMate/wsPa/NKzRvAytC2ijxuj9552qUxbNzwlO15y2SF19Qzr8+nHteOpQ9n7Kljd4ki4652TduulcSUGtsBuDubP7EEMrAGCKCHwB1IWJ2oIdfv14pjxhJGX1g5+9knvnAqUIE4l4waY1F5iua16sWzaek3OblhVNmc1su/YOamTUl+cZbbn4LN375EGNJoPyhVs3nZtZb3jdu/YOZm5DiQMAFIfAF0BNGS/AzW8Btu/wcM5GMs8z8gvU1UYmKUUIW7esUeef2ZgpZSgmMC00gOKP1p42YW0xQysAYHpoZwagZhQKcBPHRnT49eOZ0gDXAiwVyt4aSWef+lYdePWNzGNFjHR923JJyrQUm0hD1NPOm9sZMgEAFYyML4Ca0dU3lClZODHq64vfe05SsDksGvGUSu8yS+Zlb62kg0O/U0PUUzIZlBxsu/J83XDB8oKlCL89kcyp8TXSmM1lrpQBAFA5KiLju3r1artx40bFYjHFYrFyLwdAlcpkfEd9+aHjRtI7lwXtv8JDJDwT9NOVggzvdW3LdebiBWOytGRvAaD6xePxygh8KXUAMBUTBaK9/Qnd9cgLeuLF13LKE/IHPhhJm85bqj0vHsnU4jL5DABqG6UOAKqGKztwm9IK9ddtWdGkWzedq58c/I1GRn3JSNamM7vWyvOMjKzmRT392XvP1p+992yyuQBQJ8j4AqgKrozhxKifyeTmlydIygSx7vumhQ05HRbchjcCXQCoP2R8AVSFB/cO5gS9zv09A0qmbFDKkO7W4DLBrneuGxZBsAsA9Y3AF0DFcTW8bkraa0dP6NHnX81mej1JMvJ9q1R6o1rKKvP9aNJXV99QzuAHAl4AAIEvgFk11Y4IO546pK0P7cvpsxtmJJ1/RqOee2l43OuZZgYAKITAF8CsyR8oMVnXhN7+hLY+tG9Mn13HSJo/z9N1f7Bcz78S1O0aI1kFI4KjntE1rc1jeuoCACAR+AIokUKZ3fBAifzyg0K6+obGHQ0cjRhdFwpqw3W77r7U8AIAJkLgC2DGxsvstq9aooaol+moMFn5QfuqJZo/zys4JW2yyWgEvACAyRD4Apix8TK7LSuatH1Le8FsbH6G2F2m3RgAYLYQ+AKYtnD3hXBmt2lhg77w3ecymVrXVix8P9eTN+IZXf7O07X7p7+Wb21RtcAAAEwHgS+AaQmXN0Q9o0vWnJq5buv39ymZbi12f++gdt7cnpPVPfz68UxP3qRv9b1nDmfuO1JELTAAANNB4AtgSsLBqytvGElZ/eBnr0gKOi+Et6eNJn3t2juYGTU8mvTlGQXdGArsY/OMoRUZAGBWVMTI4tWrV9uNGzcqFospFouVezkAxpGf5ZUxGk2OnaYWFo0YeQVu50kyXjCEwioImCOe0bYrz9cNFyyf1b8HAKD+xOPxygh8W1tbbU9PT7mXAWASf/Pd57TjqUOykiJGuq4tCFAf6B1UMunLl+QZZUofTlk0X5K0s/uQ8ruUeUa6vm25zly8QE0LG9jQBgCYdZQ6AHWuUP/d8DFJmQ1s9/cMZLK2xjOZzWtXr1+WuU1+ANvbn9CuvYNjhk00RD1dzaAJAMAcIuML1LFC/XclFSxncDW57n+MiJcNYIuZyMawCQBAuZHxBerYg3sHM90VXP9dSdmevCkrmw513e/Ingk2oPnWFj2RjWETAIBK4JV7AQDKo7c/kVO6EIkE/Xdfev24ohFPnlyQm72PkXTROSdr25XnqyHqKWJU1EQ2AAAqARlfoA4UquPt6htSMr3jzEhat6xRWx/aJ9/aoN1YupQhqHbIljXcuulctaxo0prTFlGyAACoKgS+QA0Yr4ZWUqZ/bjLl59Tjtq9akpm2FvGM9h56XSk/XNaQbTN2XVuzzly8ICfIzS9fAACg0hH4AlWuUG/dZKpwn90To77+8Ue/1CmL5stI2nr5WiWOjejw68e146lDmcf0TFD6kEoFI4jpvgAAqAUEvkCV6+obytmM5jK14e8dK+mH6QlrktQQ9bTz5qCTw669gxoZ9eWlh0hQygAAqDUEvkCVa1rYIM8Eg4K9UJ/c8PdGUqpA50LXkeGWjedo+5b2MYEuAS8AoJYQ+AJVyg2GeKB3UCk/uwktlcr9PuIZbbn4LN375EGNjAbT1ZxwRwZqdgEAtY7AF6hCrq7X9eCVJN9KfioobQh/b63VogXzMhndpoUN2nd4ODN1jWAXAFAvCHyBKhQePCEFnRfmRYLNbKlU0KXBfe+yumR0AQD1jsAXqDL5gyeiEaPrWpt11fplksRoYAAAxlERge/w8LA6OjoUi8UUi8XKvRygouUPnri2tVlf+pN3ZK5nNDAAAGPF4/HKCHwbGxvV2dlZ7mUAFa+3P5EZKRzusQsAACYWi8UqI/AFMLn8QRXXty1ncxoAAFPglXsBAMbq7U/o7scOqLc/kTkWHlSR8q3OWLyAoBcAgCkg4wtUmHBmtyHqafuWdrWsaFL7qiVqiHoaTfo5/XcBAEBxCHyBCtDbn8h0YAi3KnOT1Vxm96r1y+i/CwDANBH4AmUWHkbhGUlGmVZlkYinpoUN+sJ3n9MDvYNKpoIs8FVsaAMAYMoIfIEyyM/wvjkaDBJOWWWiXiNp3bJGbX1on1K+zQTD+VlgAABQHAJfYA719ie0a+9gJnsb9Yz8Ardzk9j2HnpdKd/mHqe+FwCAaSHwBeZIuKTBhbIjKZtzGxfwXtPaLEna8dShzHURI1qYAQAwAwS+wBxx7chsges8I0W9IOB1ga3LDo+M+vI8o21Xnq8bLlg+5+sGAKBWEPgCc6RpYYM8YyRZGaXreRU0077onJN166Zzx4wb3r6lPVMLTJYXAICZIfAF5kBvf0Lbdu9XyreKeEZbLj5L9z55MNOTNz/odVpWNBHwAgBQIgS+wCxy3RsOv348U+ZgrdWiBfPI5gIAMMcIfIFZUKh7QzTiKZXKTl0jmwsAwNwy1hbaajO3Vq9ebTdu3KhYLKZYLFbu5QAzUqh7g5G06bylelfzYjK8AACUQTwer4zAt7W11fb09JR7GUBJ3P3YAd35g+fl5/3Taoh62nlzO0EvAABl4pV7AUCtaV+1RA1RTxET9N51Uqlg4hoAACgPanyBEgu3IWta2KBtu/dnujcwcQ0AgPIh8AVmQXjj2prTFtG9AQCACkDgC0yBa0+WH8S6Lg5G0tozGpU4NpLTuYGAFwCA8iPwBYrkujWMJH01RD1tvXytEsdG1LSwQbd/f59GUtndbJ4JNrNt38JmNgAAKgWBL1Ckrr4hjSR9+VYaGfW19aF9SvlWxmhMBwffSqPJYDMbgS8AAJWBrg5AEXr7E3rp9eOKRoJuDZ5nlPKtrMYGvVKQ8WUzGwAAlYWMLzCJHU8d0taH9sm3VlHP6P2/v1SS9J+/eFWpdNRrJL1zWaPOP7NxTI0vAACoDAS+QAFuE1vTwgZtfWifkukAdzRl9egvXpVvrTwjRTwja21Q8xtbS6ALAEAFm1Hga4xpkXSXJF/SK5JulHSVpM9IelPSx621A8aY8yR9Pf18t1lrH5nRqoFZFN7E5hmTyepKStfzWvk2yPJe19asMxcvILsLAEAVmGnG9yVJH7TWHjPGfFnShyX9paQNkv5A0m2SOiR9WdInFQTHD0si8EXFCm9ik7WKeEa+b+V5RlsuPkv3PnkwM5Di6vXLCHgBAKgSMwp8rbUvhy6OSjpX0n5r7YikHxtj7khfd7q19kVJMsYMGWNOtta+NpPnBmaLGznsglvXtsxldf9o7WkMpAAAoAqVpMbXGLNc0iZJX5B0SuiqiLtJ6NiwpJMkZQLfI0eOqLW1NXODjo4OdXR0lGJpwJS5kcNuIMWa0xblBLgMpAAAoLp0dnaqs7Nz5oGvMeZtkv5Z0icUBLpvC12dSn/1Q8cWS/pN+DFOOeUU9fT0zHQpQEk9uHdQI0lfu/YOMogCAIAq5pKqM+rja4yJSNouaZu19gVJBySdZ4xpMMZcJOmn6Zu+bIxZnQ6ST6LMAZUuXOfrBlEAAIDqNtOM77WS/lDSImPMbZK+JunvJP1IQVeHj6Vv9wVJ30g/39YZPicw6/LrfBlEAQBA9TPWFhg7NcdaW1tt2UsdBrqlg3uklRuk5rbyrgUVwfXyZRMbAAC1gQEWUhD03neFlBqRIg3Sx79P8As2sQEAUGNmVONbMw7uCYJemwq+HtxT7hWhTHr7E7r7sQPq7U+UeykAAKDEyPhKQXlDpEFKnghGcy2gnrMehSe2NUQ9OjkAAFBjyGMNz3EAACAASURBVPhKQVnDpV+VPE+yvvTw54LyB9QVOjkAAFDbCHyd40OStUHgS7lDXXKdHCJGdHIAAKAGUerguHKH1IjkRaThwSDryya3uuEmttHJAQCA2kQ7s7CBbunZHdLTOyQ/SYcHAACAGkKpQ1hzm9TYHAS9NiUl3wwCYQAAAFQ9At98KzcEpQ6SJBtkf9noBgAAUPUIfPM1t0nv/qgkE1z2k2x0AwAAqAEVEfgODw+ro6ND8Xi83EsJrNssRd8imUhQ57tyQ7lXBAAAgBmIx+NsbhvXQHeQ6V25gc1tAAAANYB2ZuNpbiPgrTO9/QlamQEAUMMIfAExrhgAgHpA4Iua0tuf0K69gzKSrlq/LBO8jpfNdccPv358zLhiAl8AAGoLgS9qRm9/Qps7n9RIKqhbv793UDtvbpck3XhPl06M+op4RtuuPF83XLA8J8sb9YyiEU+plM+4YgAAahSBL2pGV9+QRlPZzZoucytJJ0Z9WUlJ3+q27z2n/YeHJSmb5U1ZbTrvVL2reTE1vgAA1KiKaGcGzFRvf0IvvX5ckYjJHItEjA6/flxNCxsU8bLHU1ba8dQh3d8zIC993Er60QtHCHoBAKhhZHwnQ1uzipdfsvCB85ZKkh5/4Yh2PHVIEc/o8neert0//bVSvpVVEOiOpqzOPvWt+uWrb8hKSqWo7QUAoJYR+E5koFu67wopNRIMsvj49wl+K9CDewczpQwp32pd82JJ0g9/9kqmvGH3T3+tbVeer32Hh/UvPQNKpoIA+ODQ7zQvSm0vAAD1gMB3Igf3BEGvTQVfD+4h8K0A4Q4Nz798VN/uPiRX2RuJZIPXiGeU9INrfGuVODaiL//JO2QUlDpYSda3+khbs85cvIAyBwAAahyB70RWbggyvS7jy+jisnNlDSdGfXkmKFlIx7Yykj7Skm1htu3K87X1oX3yrVVDKJt71fpl2rV3UKPJIMt7dajtGQAAqF0EvhNpbgvKG6jxrRhdfUPZsoa8adsRz+jq9csyl2+4YLnWnLZoTP/elhVN2r6lnSltAADUGQLfyTC6eNYVOyrYdW7wPKOUnxv1RtP9efPv37KiqeBjjnccAADULgLfqaDDQ8lNNCo4HBBLytzOM0F211qrqGd0TWtzzpQ2AACAQioi8B0eHlZHR4disZhisVi5l1MYHR5mRVff0JhRwZK0a++gHugN6nAjntH65YszJQ5G0nVsSAMAAFMQj8crI/BtbGxUZ2dnuZcxMTo8zIr2VUvUEPUyG82aFjZkNq+5Yoakb9V9MJG5TyTChjQAADA1sVisMgLfijfQLQ0PSF5U8kWHhxLK32jmMsB2nNvnd24AAAAoFoHvZMIlDl5Eavm4tG4z2d4SckFsV9+Qjh4flWeMJJtuV2bkp6eteUZqSLcfAwAAmCoC38mESxx8SY3LCHpLLNyb19XwRtJdGlw7sqaFDUocG6GmFwAATBuB72QYYjFrXNeGw68fzylvsJJsetIabccAAECpEPhOhiEWsyLcxizqGUUjnpJJX76CkoZ5oUlrAAAApUDgWwyGWJRMfpbXt1LKt5n2ZJQ0AACA2ULgizkTzvJ6RjLGyJPVvCjtyQAAwOwj8MWcCQ+rCCYOB5PXtl6+lqAXAADMOq/cC0D9cMMqTOiY28QGAAAw2wh8MaeuWr9Mf3TeUjVEPUXYxAYAAOYQpQ6Ydb39Ce3aO6gHegeVTPlqiHq6PbaWTWwAAGBOEfhO1UA3rc2K4Lo3HD0+qnue+JVS6elrkjSa9JU4NqJbNp5T1jUCAID6QuA7FeHxxZGGoL8vwe8Y+ZPYwowobwAAAOVRETW+w8PD6ujoUDweL/dSJhYeX5waCS5jDNe9IT/ojRjphguWa/uWdsobAADAnIrH45WR8W1sbFRnZ2e5lzE5xhcXxXVvGBkNJrEZSRHPaNuV5+uGC5aXe3kAAKAOxWIxGWvz83Jzr7W11fb09Mz9E0+nXpca36K4Gl8msQEAgEpRv4Ev9boAAAB1pSJqfMuCel0AAIC6Ur+Br6vXNRHqdUugtz+hux87oN7+RLmXAgAAUFBFbG4ri+a2oLzh2R1SzhBdTFW4fVnEM9py8VlatGAedb0AAKCi1G/g6zzz7aDU4ZmdU6vzZZNbRlffUKZnb9K3+sf/6pNnpIaoR+syAABQMeo78C1U51tMEMvGOEnZUcSvHT0hzzNK+dmNkr4NJrR19Q0R+AIAgIpQ34HvdPvyTjdgriG9/Qlt7nxSI6kg2I14Qa9ePz2a2DNMaAMAAJWlvgNfV+c71ZIFBlmoq29Io6lQhteXNl/QrDMXL6B3LwAAqEj1HfhKQbA71WztdAPmGhAeTDEvYjIZ33lRT1evX0agCwAAKhaB73RNJ2CucoW6N/z2RFJG0lUEvQAAoMIR+KJo+d0b7nniV/rOpy4k4AUAAFWhIgZYDA8Pq6OjQ/F4vNxLmZqBbmnPncHXOtC+aokiXrbnsW+tuvqGyrgiAACA4sTjcRlr7eS3nGWtra22p6en3MuYmjpraeZqe48eH9U9T/xKvrX06QUAAFWFUofpqqOWZq62dyTpqyHqaduV59O1AQAAVB0C3+mqo5ZmXX1DGkn6maEUiWMjumXjOeVeFgAAwJQQ+E5XHbU0a1+1RA1RT6NJn6EUAACgalHjWwoD3TUfALsaX8obAABAtSLjO1N1ssmtZUUTAS8AAKhqFdHOrKoV2uRWY3r7E7r7sQPq7U+UeykAAADTRsZ3ptwmt+QJyRhpQe3Uv/b2J7Rr76Ae6B1UMuXTvgwAAFQ1T5KMMZ4x5t3GmMuMMe8zxiwt98KqRnObdOlXJc+TrC89/LmaGGjhWpjtfOpQTkcHBlYAAIBqFTXGdEraJOlFSUckvUXSucaYY5K+Luk+a61fxjVWvuNDkrVB4FsjPX0f3DuYGU8sSUaiowMAAKhqUUnfkvQpm9fewRhzqqQbJP2ppPvKsLbqUWM9fXv7E7q/ZyAT9EYjRte1Nuuq9csocwAAAFUraq39r0JXWGtflXTXRHc2xiyS9IiktZLarbX7jDEvSnopfZMvWWt/aIw5T0H2OCrpNmvtIyX7G1SCGuvp29U3pKQfhL1G0rWtzfrSn7yjvIsCAACYoagkGWOukfSwtfaoMeaLktZL+ltr7d5J7n9c0uWS/u/QsWFr7SV5t/uypE9KekXSwwqC5drS3Fb1Aa+TP7Di6vXLyr0kAACAGXNdHW6z1t5vjLlY0gcl3SHpa5IumOjO1tqkpCPGmPDhtxpjfqQg6/tpa+1vJJ1urX1RkowxQ8aYk621r5X474ISaVnRpO1b2hlYAQAAaorr45tKf71M0testQ9JapjmY15krX2vgszu7elj4ch4WNJJ4TscOXJEra2tmT+dnZ3TfOoyG+iW9txZ1V0dXM9eSbpl4zkEvQAAoOp1dnaqtbU1GFlsjNmtIEO7SVKLghKGbmvtumIezBhzr6Q7rLX7QscWSPoPa+17jDFd1tr29PHdkm4KZ3yrfmSxVBMT3HY8dUhbH9on31p69gIAgJrjMr7XSvoPSZdaa19XkJH966k+mDGmwRgzP33xPZIOpL9/2Riz2hjzNkkn1WSZQxVPcOvtT+gL331Otz20T0nfyrfSCD17AQBAjYkaY1zZweOSlL58QlJRKVhjzL9JepekNZL+UdJfGGN+l36MT6Zv9gVJ31BQU7y1VIuvKFXa0swNqgj37JUkzxh69gIAgJoSVRDgmgLXWUmrJnsAa+0f5x0a0/PXWvszSdURCU5XlbY0yx9UIUlRz2jbledT5gAAAGpK1Fo7aXCLIrmWZm6TW4UHwAyqAAAA9cT18TWSbpR0lrX2fxpjlks6zVpbve0JyqWKNrkxqAIAANQTt7nt/5V0oYIRxZJ0VNLdZVlRtauiTW5NCxvkGSPPSPPnMagCAADUNjfA4gJr7XpjzNOSZK1NGGOm28e3vlXJJrfe/oS27d6vlG8V8Yy2Xr6W8gYAAFDTXOA7aoyJKNjQJmPMKZL8sq2qmlXJJreuviGNJINNbdZaJY6NlHtJAAAAs8oFvv8g6buSTjXGfEnSRyR9sWyrqnYu2HVlDhUW/Pb2J/TS68cVjXhKpXzNi3q0LgMAADUvKknW2u3GmF5J71ewz+nD1tqfl3Vl1ayCN7i5vr0jSV9Rz+j6tuV0cQAAAHXBdXXYZK19RNIv3BXGmI9ba8f05EURCm1wq5DA15U4+FZK+VZnLF5A0AsAAOqC6+qw1RjzNWPM7xljlhpj4pJi5VxYVXMb3Eyk4ja4ta9aooaop4gRJQ4AAKCuGGut6+P7WUmfSh/faq3dOVeLaG1ttT09RU1Irh4D3RW7wa23P6GuviG1r1pCthcAANQNF/ieJOnrkhZJWibpW5L+l7XWTnjvElm9erXduHGjYrGYYjESzQAAACiteDyeCXxfkPRVa+03jDELJP0vSa3W2j+ci4XUZMa3QpHtBQAA9cq1M9tkrT0kSdba45L+whjznvItC7Mh3NGhIepp+5Z2gl8AAFA3XOB7sjHm5Lzr3pjrxdSkCqn17e1P6K5HXsh0dBhN+urqGyLwBQAAdcMFvncWuM5Ket8crqX2VEg/X5fpPTEaTGrz6OgAAADqkBtgsbHcC6lJFdLPNzye2JN00Tkn69ZN55LtBQAAdcUNsHiLpD+XdLGCTO8eSf9orX2zjGurfq6fr8v4lqmfb9PCBnnGSLJqiHoEvQAAoC65Uod/knRU0v+TvrxZ0j9LuqYci6oZzW1BeUMZa3x7+xPatnu/Ur5VxDPaevlagl4AAFCXXOC7xlq7LnT8MWPMs+VYUM1pbivrprZwmYO1VoljI2VbCwAAQDm5kcVPG2Pa3UFjzAWSflyeJdWwgW5pz53B1znCiGIAAICAG2Dxc0lrJB1KH18u6eeSfEnWWvvO2VxEXQywKGOHB4ZWAAAAZEsdLi3rKmrdQLf0+Fek1AnJ+nPW4SEc8N6y8ZxZfS4AAIBK59qZ9Y93A2PMW+duOTXIZXqTJyT5kvHmpMMDU9oAAAByecaYO40x7zHG/J47aIxZZYz5P4wx/yGywTPjevnKl+RJqy6ZkzIHt6ktPKUNAACgnnmS/lPSpyTtN8YMG2OGJH1L0mmSPm6tfWC2FzE8PKyOjg7F4/HZfqq553r5mogUnS9d8vk5qe1lUxsAAEBWPB4PNreVW81vbhvozvbyleakr29vf0K79g7KSLpq/TLKHAAAQN2LTn4TzJjr5TtHnR3y63uvWr+s5M8BAABQbQh855Kr97WpWens4LK8+18aHlPfS8YXAADUOwLfueTqfV3Gt4SdHXr7E9rc+aRGUtnSFY/6XgAAgIyoJBlj/tla+6fhKwodwww1twXlDbNQ49vVN6TRVG699kXnnKxbN51LthcAAEDZkcVrwweNMRFJLXO/nDrQ3CZt+GzwfQnHF7evWqJ5EZO53BD1CHoBAABCosaYo5IWGGN+mz5mJI1I6izfsmpciTe5uQltt19xvvYdHqaTAwAAQAFRa+0iY8xXrbWfK/di6kYJN7kxoQ0AAKA4rtThA2VdRb0JD7WY4SY3JrQBAAAUJ2qM+aakM40x/5B/pbX2L8qwptpXwk1uTQsb5BkjydLBAQAAYAJRSf8h6X2Sesu8lvrihlrMQG9/Qtt271fKt4p4RlsvX0uZAwAAwDii1tpvG2N+bq19ttyLqTvhUcbTCIJdmYOVZK1V4thI6dcIAABQI9wAiy3GGJt/JaUOsyjc2cGLSO/+qLRu85QC4PZVS9QQ9TSa9ClzAAAAmIQLfHvKuop6lNPZISX1fFN6ZueUWpu1rGjS9i3t6uobUvuqJZQ5AAAATMBYOybRO+dWr15tN27cqFgsplgsVu7lzA2X8U2+KSn9GpiI9L6/yQ64AAAAQEnE4/Eg8DXGnCvpryStVDYLLGvt++ZiIa2trbanpw6TzgPd0rM7pKd3SH6yJMMsAAAAUJgLcu+X9I+S7pGUKt9y6ozr7LDuhpK0NgMAAMD4XOCbtNZ+rawrqWfTaG3mxhRT2wsAAFCcqDHmJElxY8yfS/qupBPuSmvtb8q2MoyLMcUAAABTF1XQ0cGkL/916DoradWcrwiTenDvoE6MBv173ZhiAl8AAICJRa21BLeVoMhhFr39Cd3fM+D6QCgSoX8vAABAMaKSZIy5RtLD1tqjxpgvSlov6X9aa58u6+rqRXiYxQSdHXr7E7rrkReU9IOw10j6SMsysr0AAABF8NJfb0sHvRdL+qCk+xR0ecBcCA+zSJ6QHv9KEAyHuLreJ158Tb6VPCPNn+fp6vXLyrRoAACA6uICX9fC7DJJX7PWPiSpoTxLqkMrNwSZXnmSfKnv8SADHAp+u/qGNJIM6no9SRedczKb2gAAAKbABb4vGWO+LulaSf9mjJkfug6zrbktKG84+xLJeJL1gwzwwT2Zm7SvWqKGqKeIkRrmebp107kEvQAAAFPg+vheK+lSSXdYa183xpyu3A4PmG3NbdIln5f6n8zW+q7ckLm6ZUWTtm9pp3cvAADANBlr7eS3mmV1O7K4kCK7OwAAAGBqopPfBHNqGlPcAAAAMDnqeAEAAFAXKiLwHR4eVkdHh+LxeLmXUhkGuqU9d0oD3ertT+juxw6otz9R7lUBAABUrXg8To1vxQkNs/C9ebph5AvqTp6jqGd0TWuzrlrPwAoAAIDpqIiML0LCwyxSo2qx++VbaSRlteOpQ7rxni6yvwAAANNA4Ftp3DALE5G8iJZ5Q1pvXpAkWUmjSV9dfUPlXSMAAEAVIvCtNG6YRcvHJEnXeo/q22/5itqiBxQx0ryop/ZVS8q8SAAAgOpDO7NK1Nyml575gZamkooaX9Yf1RfOG9KPT7+M4RUAAADTROBboZ5MnafLFJVsUqOK6sDCd+mWjeeUe1kAAABVi8C3Qp317o36RO8X1WL3q9es1V+/e2O5lwQAAFDVCHwrVMuKJv31lo/pV08/phsiP9OZ3ouSmOgGAAAwXQS+Fai3P6GuviG9/60H9ZF9fx60N3vm76V3f1Rat5mRxgAAANNA4FthevsTuvGeLo0kfZ2Y932tiZyQsb6USkk935Se2Rl0fSD4BQAAmJIZtTMzxiwyxjxljHnDGHN++th1xpj/NsY8aoxpTh87zxizxxjzpDFmUykWXqu6+oY0kvTlW+m/k29X0syTZNLX2iD7e3BPOZcIAABQlWbax/e4pMslPSBJxph5kv5S0iWSbkv/kaQvS/qkpA9K2jbD56xp7auWqCHqKWKkfZG365cf2iG13iRF5gdDLSINwZALAAAATMmMSh2stUlJR4xxGUmtlrTfWjsi6cfGmDvSx0+31r4oScaYIWPMydba12by3LXsqvXLZNJf376iSdImad0NQaZ35QbKHAAAAKah1DW+iyX9NnQ5kv5qQseGJZ0kicA3T7i+tyHq6ar1y7JXumD34B7plZ9Jx4cIggEAAKag1IFvQtLbQpdT6a9+6NhiSb8J3+nIkSNqbW3NXO7o6FBHR0eJl1b5wvW9o0lfXX1D2SltA93SfVdIyROSfMl4QfkDG90AAAAm1NnZqc7OzpIHvgcknWeMaZD0B5J+mj7+sjFmtaRXJJ2UX+ZwyimnqKenp8RLqT6uvnc06Wte1FP7qiXZKw/uCTa2ud8hrJ/d6EbgCwAAMC6XVJ1x4GuM+TdJ75K0RtLXJP2dpB9JelPSx9I3+4KkbyjIMG+d6XPWqpYVTdq+pV1dfUNqX7Ukm+2VgrKGSENexpeNbgAAAMUy1tpyr0Gtra2WjG8RBrqDDO+CJdLLz0gyDLQAAAAoEgMsqklzW/BnoFt6+HPpiW4MtAAAACjGTPv4okR6+xO6+7ED6u1PTH5jV+9rUwy0AAAAKBIZ3wqQ38Zs+5b23PrefK7eNzVCnS8AAECRCHwrwIRtzAppbgvKG1y9r8v4Uu4AAAAwLgLfCjBhG7PxuCD3viuymV9qfQEAAMZF4FsBJmxjNpFCtb4EvgAAAAUR+FaIlhVNxQe8DrW+AAAARavrwLe3PzH1LGslcbW+z+6QZMq9GgAAgIpWt4HvlDspVLJnvk1PXwAAgEnUbR/fQp0UqhI9fQEAAIpSEYHv8PCwOjo6FI/H5+w5XSeFiFHxnRQqkavzNRHqfAEAAMYRj8dlrLXlXodaW1ttT0/PnD9v1df4OgPdQaZ35QbKHAAAAMZR14FvpaiZABwAAKCC1e3mtkpR0k12ZH4BAADGReBbRr39Cd31yAtTG1c8noFuprgBAABMgMC3TFym98SoLyvJm+kmO6a4AQAATIjAt0xcOzWroLXGReecrFs3nTv9MgfX3SF5QjJGWlClXSoAAABmSUW0M6tH4XZqDfO8mQW9UpDdvfSrkudJ1pce/lxQ/gAAAABJZHzL6qr1y2TSX0vSzeH4kGRtEPgmTwTB7+nrpHWbKXsAAAB1j8C3DPI7OVy1fllpHjhc7iBfeqk3+PP0dumm3QS/AACgrlHqUAazNi65uS3o5nD2JZJM9jijjAEAAAh8y2FWxyU3t0mXfF6KzMseY5QxAAAAk9vKZdantQ10S8/ukGSCGl+J4RYAAKCuUeM7x8IB7y0bz5m9J2puywa4PfdK//bZYNNbZD7DLQAAQF0i8J1DJR1PXKyB7iDo9ZPB5dQJhlsAAIC6VBE1vsPDw+ro6FA8Hi/3UmbVrG1qm8jBPZLvZy8bj3pfAABQd+LxeGVkfBsbG9XZ2VnuZcw6t6ltNOmXflPbeFZukKLzgxZnnif98Z1kewEAQN2JxWJsbptrs76prZCBbja2AQCAulcRGd9alx/szlnA64Q3ugEAANQpAt9ZFt7QFvWMrmltLt2I4ukg+wsAAOoUge8sC29oG0lZ7XjqkHbtHZybjg75Brql+64IJrlFGmhrBgAA6kpFdHWoZW5DmxsgbDWHHR3yHdwTBL02JSXflH58l7TnziAgBgAAqHEEvrOsZUWTtm9p1+YLls/emOJirdwgeZH0BSv94l+lR/82yAIT/AIAgBpHqcMccBvarl6/bO47OoQ1t0nv/qjU800FuWcF09xSIwy1AAAANY/Adw6VpaNDvnWbpWd2Bn195QcDLSINDLUAAAA1j8C33jS3BZvaDu6RFiyRXn5GylQgAwAA1C4C33rk+voOdEsPfy4odXhmJ10eAABATWNzWz0Ld3lwdb4AAAA1isB3lvX2J3T3YwfU258o91LGWrkhqO81Eep8AQBAzaPUYRaFp7Y1RL3yDK2YiKv3fXaHqPMFAAC1riIyvsPDw+ro6FA8Hi/3UkoqPLWtbEMrivHMt6Xe+6R7L5N2/w96+gIAgJoTj8crI+Pb2Niozs7Oci+j5NzUttGkX76hFZPJqfNNBT1+2egGAABqTCwWq4zAt1a5qW1lHVoxGVfnm3xTwVALG3z/7A4CXwAAUFOMtbbca1Bra6vt6ekp9zLq10B3EOju/ZbkjwbHIvOlm3YT/AIAgJpRETW+KLPmNunyu6T1f6rMJrfUqPT4V6j3BQAANYPAdxb09if0he8+p7/57nOV2cZsPOs2S9G3KHhb+FLf49J9VxD8AgCAmkCpQ4n19ie0ufNJjaSCn2tD1NPOmyusjdlEBrqDTG/f45L1gx6/LR+TGpuDemBKHwAAQJUi41tiXX1DGk1lf5mo6DZmhTS3SZd8PqjxdW+Pvd+SHv0S2V8AAFDVCHxLrH3VEs2LZIdBVGwbs4k0t0mXflXyvCDr648y1hgAAFQ92pmVWMuKJu3suFC79g7KSLpq/bLqKXMIOz4k2XR7M0mSYawxAACoagS+s6BlRVN1Brthrr9vakTyItLqD0hvPbXcqwIAAJg2Al8U1twWTG87uEdasER6+HNBEMxUNwAAUKWo8cX4mtukDZ8Nyh4yY42p8wUAANWJwBeTc2UPJkKdLwAAqFoV0cd39erVduPGjYrFYorFYuVeDgoZ6A4yvfTyBQAAVSgej1dG4FtLAyxq3kC39OwOSSaY9EYQDAAAqgSb21C8gW7p3suCOl9Jenq7dNNugl8AAFAVqPFF8Q7ukVKj2cupE0G3h93/g4luAACg4pHxLaHe/oS6+obUvmpJ9ffxLWTlBikyL5vxlaSXeoM/e/9ZOveDQa9fSiAAAEAFosa3RHr7E7rxni6NJH01RD1t39Jem8Gvq/H99U+ll/YqO9ktJDKfEggAAFBxKHUoka6+IY0kfflWGk366uobKveSZkdzm3T5XdKlXw2yv4XQ6xcAAFQgAt8SaV+1RA1RTxEjzYt6al+1pNxLml3NbdJN/yq1fkJ6++WSFwqCvag0PEjdLwAAqCiUOpRQzdf4TsSVQLxxRHrxh5KfDIZdMN4YAABUCDa3lVDLiqb6C3id5rbgz547pef/PXe8MYEvAACoAJQ6lEhvf0J3P3ZAvf2Jci+lvBhvDAAAKhQZ3xKom44OxWhuC8obGG8MAAAqTMkDX2PMSkk/kbQ/fegaSe+T9BlJb0r6uLV2oNTPW06FOjrUbeArZcseAAAAKshsZXx/ZK39iCQZY+ZJ+ktJGyT9gaTbJHXM0vOWhevoMJr066OjQzEGusdmfQsdAwAAmCOzFfheZIzZI2mPpO2S9ltrRyT92Bhzxyw9Z9m0rGjS9i3t9dvRId9At3TfFcHmNtfZQco9dulXpeNDBMEAAGDOzEbg+2tJ50g6Jun/k3SlpN+Gro/k3+HIkSNqbW3NXO7o6FBHR/Ukheu6jVkhB/cEAa5NSckT0uNfkZpWho4dl/71L4Pb0vIMAADMss7OTnV2ds5uH19jzIckfVLSUWvtJ9PHnrTWXhi+XTX38WVjWwEu45s8IcmXZCSTbiBiU7m3NRHpfX8jbfjsXK8SAADUmZK3MzPGLApdfI+k3ZLOM8Y0GGMukvTTUj9nOdXNqOKpcJ0dzr5EkpFkg4DXpoPgMGOkBUuCYHnPnUx7ktRRugAAH5FJREFUAwCg1lTQZ/xslDpcbIz5WwWlDr9SsJntTUk/Sn/92Cw8Z9mwsW0czW3SJZ+XfvVfwRQ3x4tI1gYBr7VBMPzvfy3JMO0NAIBak7/vp8x7fBhZPAOutrdpYYMSx0ao8S2k517p3z4bBLiR+dk3/PCg1HtfuvTBZYGtJC/IFF/y+eAfBJ0gAACoXJN9Tu+5U3r0S9lSR5NOgHme9Md3Sq03zelyGWAxTdT2Fqn1JmnpeYVbmz2zM/gNMJz9lS/98rEgU3zhp6Wnvp7bHYLgFwCAylCoi1P4c3qgWxoeCPb5uMDXffX9IDG29Lw5/Wwn8J0mhlZMQaGBFq4O+Nkd0tM7JH80GwDLBmUP//0PwW2tn+0O4TLBAABg7g10B5/d7mxtfhenSz4fHM98vrtyx/SenzDrB4kxAt/KR21vCTS3BW94PynJl6yX/q3QD663frobhFFOJrgMp0YAAKh7A93SvZcFwa4UlC0E30jypb7HpYNPBJdTI8oGut7YPT6yQQnkyg1z+lcg8J0mhlaUyMoNwekRd5rkgk9JT/5vyU8p6AZhx2aCy3BqBACAundwj5QazV7O1O16QYxr/dD1Lug1UnR+7qY291hl2L9D4DsDLSuaCHhnypU8hP8BvP2y4HRJ3+PBPyJrxmaC5/jUCAAAdW/lBikyL5vxdayyGV0voiDjmy5hXPMh6aLPFC55LAO6OkwT09pmWX7BvMsEW1/yotLqD0hvPVVatzm4vas3WreZgBgAgFJz3RsWLJFefkZ644j04g+zrUgv/WpwXEaa/7bsZ3ZkfkVtTifwnQY6OsyR/BYprqB+77eCzXCS5M0LvrrLkfnSTbtphQYAwEyFg92HPze2e0P4c1bKJqyMSZ+x9StuQiulDtNAR4c5kt8NImczXJo/mnuf1EhwGyk7Njm/VyABMQAAY40XyErpckOb/Zx1n9Huc3TPndkOD9YLPntlgkB5jjewTYTAdxro6FBG49UXOZGG4DfTx78iJd9UsCEu1CtQmrjnIAAA9WS8rO67NmcD2TAvWjiQzd+sXuYJbeMh8J0GOjqUUXObdNO/BiUPv/6pdPjp7Ka3M1ukd38s+IebPKGcfoFuQ5yU/Ycc/q0VAIB6E95PEy5PSI1Ib7ySPhbuv2ukd99Q+HOz0Gb1ClQRge/w8LA6OjoUi8UUi8XKvZyi0NGhjNyplULzv5/dmc30ygT/aKXcXoGRhiAwNib4DZfSBwBAPTq4J1SeoHRfXi+47oUfBGdMTajHfqRBWnfD+I9XaGBVBYnH42xuQ5XLr0cKN9aOzJc+9H+NPdXSc29Q+uAGZIQbaVP6AACoB4U2jLuBFK6e1x1r+ZjU2FwTCaKKyPgC05ZfWO+7WqT06ZjwhrbdtyozYtHa7CkdJ/kmY5EBALVpoo1rbpxw5jMxVNrgsrw18rlI4IvakV9Y707H5I9Y9OYFxfmp0G+0UvB93+NS/5NkfgEAtcOVBrpOR+deGtq45o0dPuEng+/f/dGa649P4IvaMV5hff6IRT+ZzQQ/vSM9XUbZLDCb3gAAtWKge2yno+f/PUgA+RrbgUGq6X0vBL6oLYUK6/NboIVP26y7oXAblwrqOQgAwLS4PS1+SmM6HZ12vnT6usIZ3RoMeB02t6E+uCL+icYa090BAFArBrqlb34od+hTeEO38epyUzcZX9SHYlqsVHgbFgAAiuLKG/zQ8AkvGkwx/flDwX6WOi3tI/Cdot7+BIMrAABAZRpT3mCCjWp/fGewv2XpecEm7jot7SPwnYLe/oRuvKdLI0lfDVFP27e0E/wCAIDycyV9vf8UGjNspLM35rbprJIJa7OFwHcKuvqGNJL05VtpNOmrq2+IwBcAgFoy2X6PStkPEt67ctq6YIN2ZnJpmhcp3Ju+jkv7CHynoH3VEjVEPY0mfc2LempftaTcSwIAAKXi+t2mRsb2sXWB5tM7gg1jkYbybQzL708f3rTmuJreOg1wx0PgOwUtK5q0fUs7Nb4AANSig3uygx1SKannm9IzO6ULPiU9+b9z24Il3wwC4bkOLN3GtXB/+vAUUm+etP5Pa27wRKkQ+E5Ry4omAt5aVimnsAAAcyP8//6CJZIxkk2P8JUNpp399z/kBpdScN3eb2nCNpmlNl5f3gwjrf+odPnfzf5aqhSBL+CET3GV8xQWAGBuhEf5GhOUDPh+8NV4QbBrTG5bMJn0Vyv5o9ms8Gx+Zoy3ce3M9dKi06UXf5gtv1h3w+ysoUZUROA7PDysjo4OxWIxxWKxci8H9SrnFFf99TYEgLqSP8rXKhRURoLMqSS9cSQILFOjkudJF35aeurroY1kdnY/M8bL8nqRYNSwqz/mbOWk4vF4ZQS+jY2N6uzsLPcyUO9Wbgh+W67T3oYAUFUKBXvFBoDhTO+YkgETfAac9q7sKHsvEvTAdSUNb79s7Ea3lRtKE4CGH0NKB73J3Nvkb1yr4y4NUxGLxSoj8AUqQp33NgSAquGyoNYPxu5e+lXp5WeCQNRlZv/4/2/v/mOsKtM7gH+fc4dh2boCArsqsk5RwOI2iDNVqCEL1rKI0iY26UaWIiSGJq1/oNZEY7M1Tbdrlk7iJm42O93sql3A2qbp7sCC2lW3I/4g4ACrq7CI4qDgUhYHjTIw97z9433fe97z3nPv3Htn5p577/l+EsL9NXPPmTlz73Of87zP062HNbhBpL1sz/AhBCC6nAHQAaXt5OCeBQwBTL6sONBcsDpqKfbWdr0Azm5TLaUPfsndNbfr0gtLckDnHVy4NgoMfIlc/NRMRNTYBnbHs6DDnwHb74638wpDYPs9OpDNn3fqd/M6KL76NhPsBlHg/Nmp4qRHJWcB9z1pMsdOgJofqq30wQ22h4eA4/uB3IR4MN+1rrrvSTEMfImIiKg52Lrc2GIzJHRcMLflh8xlRPW7YQj86ikURvmueDg5mKzkLGAsc+yQoPpyuYHdwOCACdZD/T3ff00Hvm6ZBY0KA18iIiKqP7+W1ZYMuAGe/5hYXa4pUYgFvUF0mwT6cUlBMWDuUzrTW8pIZwHt2pBCxtcE04vv0vuzf6veH6C45GLStCjLbPctfy4qu7DbGA7HyyxoVBj4EpXCVbJEROPDbyMGRBnZ/s3Aum36sn1MEOjX4UInhQC4YinwB39uRvUORaUAQFT/KwF0+zEVTTdTYXR9tAuZ3aywDWQnTQN23BdNVdv7uHlwqW2YCFz5J053CdtD2Kgle0wlMfAlSsKevkREY88mFAaPRVlSv6lCfkiXM0ztiB4ThsDRl6LH5CYASx/Qr8vuArZZ1wF93U5wKUDnWmDyrNLZ1tG+tvtZ4b5ub6qaU5bhZ59VqGuU3/o5Cj+IoE0Hu25dL99/xgwDX6Ik+7dGn76Hh3RG4ZIFxaes+GJERJTMP2tWNCwCycPHAODt580QiaTHCLBwdelWXn5rygWr4/eP9+t2xxKzIO1cFV9kA2IzeW3Bar7PjBNRqtRRVz9dXV1qz549aW8GkTawG3jsluQXrWCCWRk8rOu4bNsbvjAREUXcs2b2tRLQbcgKQZ6px4XSWc45y4GPj+sFXYVo16nZdUsDRjoLl3apmp205g6/8DtLzLoOOPoyon0VoO1zPMM4zpjxrcDeo6fxypFTWDR7Gjovn5r25tB4e7eveMWwFZ5HIQWRz0ejKku1wiEiyqLYJMw8sOdHugetm8EVxMsQbFb4xzdHrcrcx1RTnpB2a0r3+f0Fev5iPffDARMp446B7wj2Hj2Nb/zwFZwbDtHeFmDznYsY/LY69zRZkNO1YqGp15Kc/t9mH6BMD8l79O25dgbBRNSaqsmiTppmMrVOpKvy0BncnL4tNzG5DGFld3w4hf+YZuMH4e5lDk2qOwa+I3jlyCmcGw4RKuD8cIhXjpxi4Nvq/N6NQMIpq0Dfbhct2P/9IJinrIioFVSz4Hdgt14XESa0FEvK8vq61hUvWGtVaWemM4iB7wgWzZ6G9rYA54dD5ALBBx99hr1HTzP4bXVJn9D7uoGDO6Dr03LApdd4tWhGIQg+qwNmvqgRUbNzF/zmz5WfShYb6pAD5t0MHHq6ugwuA0IaJ7mHHnoo7W3Apk2bHnr99dcBAPPmzUt5a+IunTIJi6+YjiAQ/Pr4xzhw7CP8dN/7WHzFdFw6ZVLam0f1FLQBB57SpQ+5dmDJvcCRF/R1mwH2ffgmMPurwOSZdd1UIqKaDOwG+jbpQHXSRfq1y2Zw7Yf6oA1o/z3g3Rejx7j818pV39X1q1M7gKX3M6Cl1PT29rKrQ6W+9/xhdD9zEKECcgLcs3we/nbZlWlvFtVbUnse2xNy5/1RJ4jCCuQccOODOkgmImo0/sIrt6NNbiJw83eAN39qPuQ73RhsZ4ZgAnDtXxUvykq7qwJRCSx1GIHt6DD18+2FkocJbQEWzZ6W9qZRGpJKIOx1W5PmBsG5dn29r3v0bwB8IyGiapR7zbDttvq36A4KuXY9Pcxt45gf0ovMwhCAqde17bis8Hzp7jZ8naIGxIxvGX5Hh2/eejVOf3qObc1oZPYN5+wZ4OVHdaYkaKu9XQ0nyRFRpZKCWvc1ozBIwo7/tZxMLhB1X1AhYiOC3XG87tcGgQ6S7bSxrnXjt49ENcp8xrdcj16/o8PpT8+xvIEqY99g3H6U+XNRZuSOn+nbRsrg2jew4wd09kWFIy8sIaLWN1LZlRvU+q8Z7kK1GCfondmpv/fLj+rruYnxEcFFnW4QvdaFoc4Uf2k+X6eo4WQ68B2pR6/b0YHlDVS1d/vMKUKXWRG965H4KucVDwMn9gEQ4OIF+nThpGnFmRUJdPbG1uMRUXMoNcSglsDQPwO04uGovMqdcmaJ6NcT+7X9P4nul0AvQvOnhy1cG7UkCwL9HEkjgvc8ZsohvKE/KuQHdGpImQ58R+rR23n5VGy+cxGntlFtOpYAbROdufSBKXnIAQd3xnv/btuI+BtVQi0dAMxeGmVdiKhxuBlXf4CNP74XUlyCUE1g7E5FGz4L9D8RXVcSvd7YoFaFOogF9EI1m5mF6FZjh58rnh7mtiRTovcpyWennMDZPDegP9DzAzo1oEwHvpVkdDsvn8qAl2pTahDG8QOm/6/L7wUcRnV19j73VCMRja9qFpMWamaHUFgEFrQBc5YDF3xRP6YwvtfJxtoSBGDkwNjVsUQ/Lp/X3+uDfeYOMyVNQd8/d4XuPa5CHSBvvzsKUiXQryk3bNT/kvbVTrAsd5bJnXTJyZXUBDId+DKjS+POn9e+70nz5uhkR5RfDuEIcsDclfrNkzPcieqj2sWksYENiGrx39qmrwcTdCCcV4gvHgv0h+C3tjs1/HkUAli/NtcNxheu0WsGoKKzRxKYmFrpAPeCGXr7bUCunHKG2UvjH6T9/fM/uJfa/0ofR9QgMh34AszoUh3F3hydFdI774/KIebdrB/71nYUVlPPvDbqA8yWZkTjzy0lSAo+929BUT1+rt1ZMGYCVys8rxeLAcAH/dGH3TAfBccxzgfjs2d0KZS7kCwIgMV36Vpcd5GazfTawRELVut/L3wbePv56HFBrrKzR5W2JGPrMmoimQ98ierGPyXorpD2V2cf/kXxKUa2NCOqzGg/INq/VfuB9OwZ3Ys7acEpRAeSV98GvPFfOpi1kxyVU6P/Qb/O+gYTdBmDSHENv/s97f27Him+Owx1t4WV3XpRrNu2LKnUYOkDwNGX9f7YVmN87aCMYh9fonqq9A056XF93cBz39JvppwIR5Ss3AfEagLiWLeCCkqTCsFuqP8+O9fq68cPRFlee/vkWcVBdDDBLGgdjp6nXBmUBMCNf69fAyrZL54tIgLAjC9RfY3m1KGfMeaKacqCatuAlSpTSGoBZssUkhZjxboVAIW62VKUiham2TKDpOe1twNRP1yIruG321/oxWtqcy3JRdvidk2o5HWF5QhEABj4EjUPdxHJpGnRanAgORhwaxG5MI6ajT99zO924J7SB6K/AfcDYpADBo9FwXOhBVjCKF7bT9t+T7fcwQ0+Ibru/pIFwMQLo8mM/te7i8ZKLf5KCkaTRqDbHt9ucMzMLVFNWOpA1Gzc1kluf2A/i+WeRs1NBNZtK+4X2spvnPXaz6z8POupqMzAsovGzHhcN8tqF32t7I4yqf1botttDa4K9eP9IQ+AzqgqFf8+/ujx3MTayyeIKHUNEfjOmTNHLVu2DKtWrcKqVavS3hwA5UcZE6Wqrxv4xT8hnoUy7KnQwsIZ5+97ZifwhUviK8NXdgNd6+qw0XVW60LAamslP/y1DtCSAiKqzcDu+KhvX2EwQxgFsbHhLzmg8w59ec9jKMrWJi08SxK0Aet3MMAlaiG9vb2NUeowefJk9PT0pL0ZBSONMiZKVccSHbQWjUNG9GaugNjwCwB4f2/8sWGog7YvzW/ON/Kk2k9brzl4LF7nuX9L6aDFnbhlx76WCpb9bLttOQfoPqyleq6O9c+3FYKwUvtQNOrbTAKzmd5LF8Y/wPm/B5XX/W1zE8wiMfdJncfZ7z1jLnDyEBKHyLi/T9bIEjW9VatWNUbg22hGGmVMlKpZ1+lMrc00ljpta4dffHzcTIpLOLvjvrn7gUi9g6taJmW5k67y5xGbmlUYGADgtZ/Ey0FszaRbo+lmEUv1bj1+IOqb6v84JdDBc193+VPjo/0Z+afwF98FfO7C9IPgUr+/pOPK3wf3zIM76tvu36s/iH7XJ17Xk8pEoqBXAsT/DpQ+4zFvBXDo6ejshz023Frh6/8m+sAjEi1o48hdopbEwDdBJaOMiVLVtS6++MW+cQPRm78dftGxBHjsluh+d2V40KaD4ie/Udwc3wYb5bKfYxUY25rOSgPF2Mp9L+hXYXHQYzPhw0PA9nuST3PbelEFHQBNmhbtp/vzSxK06RrSpLpUPxNsJQ1CGKnE4vE/iw8sCEPT59X0kk2rdKVUaYmbIXePK38f3DMPSYvBrrpFXx88Bux9XP/+lJsJlqh9mNvT1h3H63ZvAOLf3+2l7d9HRC2FgW8Jt117GcT8z2wvNST31GtSEGxbns26Dli3vbht0v4tOhPqT44KQ2DXd82VhLGpwNgO0xjYbQJGU9NZKlB0TZpmsnNefSeAQr9VONk7e3sQlB4aEOR0YGYztTvuAw4/qzPm+fPeU5g60aBNj469eEF8H/zH+pnDpGDadha442f6ur+i/5PfxgPGGKWfO43SlYHdejJYYeSuU1oyeCzqiuAfV64wr7+HHerilxXY6wO7gX1b45l+G+TaNmELVhcHrkk/j3LdFRjwErUsBr4ev773tmsvS3uTiEaWFAT7b/z+m/m7faUXEMUWCwU6K7zt7iho9gOdkQJVt47Wbc0067rimk43UPRbsgFOi6vQqeE0Wb/LF+uvUSHwm2d0YBpCB0lzluuvP/S0HiEbPWEU9J44EC93SBonm5sI3PydeIa2rzuh5trJwib97P1g2g0a9z2Z0EbL/da56Gtio3G9ANI31uUrbka3UGaSi7cgi9XZeovQ3H048oKeLlbuQ5SfDQYqO9aJiAwGvo69R0/jkf85xPpeam6VvvF3LNELgPwSiKRgygaAex+PP0aC0sM0/EVjfiDXv1kHkIMDZjucek+b3XOzov5zA4AygZZSejtmXAW892pU7uBOybKZ8CAHXHUrcMEXoxKDom1MyCTP7NS9W5N6IifVpbp1t/5CvMEBE5R7AbgIcPJgmcyueVzn2iizaeuJbYlFqQCy2nIS93foB5n2sl1EqD9dALOXAlM7onKEEMV1tv4+vPBtvc0q1Pu9f0v57WJ2lohGgYGvYTO9Q+dDvXZYwPpeam3lSiBsxq7QFs2I1caKDnSWPqCvbtsY/z62FKIw5tXLXuaHgO1368tBm65NtV+7baNeSOZmRYvqckUHm/4gA3sq3D393dcd1QSH0LXP7rhne78N4GYu1M9vA1M7nKDSTGSpgNNdWOUG4HaRXZgHjr6U/BxAVA5h98s+z1W3lA8g/XKS4bPlM8P2a2K/Q1M24u+Dzarn2qNjYd/WqPPFlct1na17XLm/m6UPAO++aH7+Ctj7BHDxNa3ZZo+IUsfA17CdHEzDHNxw5XRsvGkus73U2kpNjrLZRH8QhivIRYGOn5n98vVO9wM7bAAoCn5ta6lwGJh8WfH3KsX2ak3KviYFoCONe/bvX/Gwvr2ayXdJP8ui+mVnv9wAvK/bq0eGft6rVurAeKQFcEkBZP+WKLgsahGmgLefB9753/igBvd7uwsI3c1y9yGv9IcEmwm3X3f9X0e10jvv17+TWx8pXX+7cA2w50dm0/LN3WaPiBoaA1/D7+TAoJcyza8Z3r8F+ORk8fALwNT7eplZN2sZ5IC5X4tKC07s09lUt8WaretNqn2tdvBGqWC+VEa23P2jDbyKAk5H0BYF4P54XJvZvWFj5dtQCCB/jMJiNzvW2i0nEUTlIuGw7nIRtEUZ3IVrdBBb2KZyZRch8EG/HuRx8QKvLVhCa7hSZTgLbgdeeyL6gOD30CUiGiMNMbmtUUYWc1ob0Qj8mk93YVMptna4XJsrG8j6db3NPmo5NvACTlZX9P7e+kj8sX7brWr30++2seLheG1zYgcKty2YEbQlj/6N7YP9WqV/x7O/Chz5pSlJcUYKV9r1o5YaZCKiKjHwJaLa9HUDz30rCnSmzwFOHY4vfHOHQkgOuPHBqLa23MCDakoMGl01k+HG8vk6lgD7t0YZYPfnn1R3bEskLHdkb9I++O3E3CDbXq82gG/GDzdE1FTqFviKyCYAiwC8B2C9UqpQKMbAl6gJJfXyBeJZy3oFe82inoFduey5uy22xdwnJ4GDO5xFhAJcsax4AVypTg/NmpUnokypS+ArIgsB3KuUWiMiDwJ4Rym1xd7fCIEvyxyIalBJoMNgqP7sUAnb5SGptMI+zv3wYheluSN+WXZARC2kXovbFgN4xlzeCWA9gC2lH15f/tCKzXcuYvBLVIlKegZzoEB9JQ2VsC3QfLHRz+d07+H1O+JBcyUDSoiImkRQp+eZAuCMuTwI4CL3zpMnT6Krq6vwr6enp06bpdlWZu7QCiKipmSDWXeoRKmMre3cILn4iOulD+hg2b2diKiJ9fT0oKurq24Z39MALjSXpwD4nXvnjBkzkGapg9/KjEMriKhp+T2Jyw2pKNfGrVz7NyKiJrNhwwZs2LAhtRrfI0qprfZ+1vgSEY0h1lUTESViVwciIiIiyoS6TW5TSt1Xr+ciIiIiIvLVa3EbEREREVGqGPgSERERUSYw8CUiIiKiTGDgS0RERESZwMCXiIiIiDKBgS8RERERZQIDXyIiIiLKBAa+RERERJQJDHyJiIiIKBMaIvAdHBzEhg0b0Nvbm/amEBEREVEL6u3thSil0t4OdHV1qT179qS9GURERETUwhoi40tERERENN4yH/j29PSkvQnUQng80VjjMUVjiccTjbVmO6YaJvBNq743jV9YWvuapefN0vEE8Hfbys/L16jWfN4sHU9Atn7GWXveZnuNynzgm4as/VFk6QUvLfzdtu7zpiFrP+MsHcdpydLPOGvPm4bR7GtDLG4TkZMAPgIwmMLTTwfwf3V+zslIZ1+z9Lxp7WsaxxPA320rPy9fo1rzebN0PAHZ+hln7Xmb6jWqIQJfIiIiIqLx1jClDkRERERE4ynTga+IbBKRPhHZLCLtaW8PNQcR+YKIvCoin4jIV8xtXxeRl0TkORGZZW6bb46vl0XkpnS3mhqViHSa4+SXIvKUiEzg8USjISJfEZFd5pjaLiIX8Jii0RKR201palO/52W21EFEFgK4Vym1RkQeBPCOUmpL2ttFjU9E2gBMBbAJwL8AOAjgRQBLAPwRgDuUUhtE5L8B3AfgQwA7lVJ/nNImUwMTkYsBnFFKfSoi/wygH8DfgccT1UhEJiilzpvL/wDgCIC7wGOKaiQiAYD/APD7AK5HE7/nZTnjuxjAM+byTgAN+QuixqOUGlZKnXRumgPgDaXUOaXULgB/aG6/RCn1G6XUGQCnRGR63TeWGp5S6oRS6lNz9TyAueDxRKNgg17j8wDeA48pGp3VAP4TQIgmf8/LcuA7BcAZc3kQwEUpbgs1N/dYAoCc+V+c23iMUVki8mUAN0FnUng80aiIyJ+KSD+AZdAfqHhMUU1EJAfgLwH8u7mpqd/zshz4ngZwobk8BcDvUtwWam7usQQAefN/6NzGY4xKEpELAfwbgPUAfgseTzRKSqlnlVILobN0XwWPKardGgBPKaXs8dLU73lZDnxfAbDcXP4agF0pbgs1t8MA5otIu4jcAOCAuf2EiMwxQc1FSqk0emdSgzPZlM0A/lEpdQg8nmiURGSic3UQwCfgMUW1mw9grYjshC5z2IAmPp4yu7gN0F0dACyCrn9ar5Q6l/ImUZMQkZ8DuAbAUQDfBzAEYCOAswDWKqUGRGQ+gB8AaAPwTaXUs2ltLzUuEbkdwKMAfmVu+r75n8cT1UREboVeZBQCOAlgHYBV4DFFoyQie5RSXSLydTTp8ZTpwJeIiIiIsiPLpQ5ERERElCEMfImIiIgoExj4EhEREVEmMPAlIiIiokxg4EtEREREmcDAl4iIiIgygYEvEREREWUCA18ioiYjIh0i8qaI/KuIvCEiz4jIpLS3i4io0THwJSJqTnMAfE8pdTWAjwD8RcrbQ0TU8Bj4EhE1p3eUUvvM5b0AOlLcFiKipsDAl4ioOQ05l/MA2tLaECKiZsHAl4iIiIgygYEvEREREWWCKKXS3gYiIiIionHHjC8RERERZcL/A/mE+bnhS9pWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Part two of step 7 of the algorithm\n",
    "min_norm = widget.result\n",
    "nr = np.arange(W.shape[0])*stride + start\n",
    "coords, weightmatrix, DX, DY, row_mask = threshold_and_mask(min_norm, W, DX_DY, nr)\n",
    "#Step 8 of the algorithm: reduce the shift matrix to two vectors of absolute shifts\n",
    "dx, dy = calc_shift_vectors(DX, DY, weightmatrix)\n",
    "plt.plot(coords, dx, '.', label='dx')\n",
    "plt.plot(coords, dy, '.', label='dy')\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('shift (pixels)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<tr>\n",
       "<td>\n",
       "<table>\n",
       "  <thead>\n",
       "    <tr><td> </td><th> Array </th><th> Chunk </th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th> Bytes </th><td> 6.59 kB </td> <td> 160 B </td></tr>\n",
       "    <tr><th> Shape </th><td> (412, 2) </td> <td> (10, 2) </td></tr>\n",
       "    <tr><th> Count </th><td> 43 Tasks </td><td> 42 Chunks </td></tr>\n",
       "    <tr><th> Type </th><td> float64 </td><td> numpy.ndarray </td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</td>\n",
       "<td>\n",
       "<svg width=\"75\" height=\"170\" style=\"stroke:rgb(0,0,0);stroke-width:1\" >\n",
       "\n",
       "  <!-- Horizontal lines -->\n",
       "  <line x1=\"0\" y1=\"0\" x2=\"25\" y2=\"0\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"0\" y1=\"2\" x2=\"25\" y2=\"2\" />\n",
       "  <line x1=\"0\" y1=\"5\" x2=\"25\" y2=\"5\" />\n",
       "  <line x1=\"0\" y1=\"8\" x2=\"25\" y2=\"8\" />\n",
       "  <line x1=\"0\" y1=\"11\" x2=\"25\" y2=\"11\" />\n",
       "  <line x1=\"0\" y1=\"14\" x2=\"25\" y2=\"14\" />\n",
       "  <line x1=\"0\" y1=\"17\" x2=\"25\" y2=\"17\" />\n",
       "  <line x1=\"0\" y1=\"20\" x2=\"25\" y2=\"20\" />\n",
       "  <line x1=\"0\" y1=\"23\" x2=\"25\" y2=\"23\" />\n",
       "  <line x1=\"0\" y1=\"26\" x2=\"25\" y2=\"26\" />\n",
       "  <line x1=\"0\" y1=\"29\" x2=\"25\" y2=\"29\" />\n",
       "  <line x1=\"0\" y1=\"32\" x2=\"25\" y2=\"32\" />\n",
       "  <line x1=\"0\" y1=\"34\" x2=\"25\" y2=\"34\" />\n",
       "  <line x1=\"0\" y1=\"37\" x2=\"25\" y2=\"37\" />\n",
       "  <line x1=\"0\" y1=\"40\" x2=\"25\" y2=\"40\" />\n",
       "  <line x1=\"0\" y1=\"43\" x2=\"25\" y2=\"43\" />\n",
       "  <line x1=\"0\" y1=\"46\" x2=\"25\" y2=\"46\" />\n",
       "  <line x1=\"0\" y1=\"49\" x2=\"25\" y2=\"49\" />\n",
       "  <line x1=\"0\" y1=\"52\" x2=\"25\" y2=\"52\" />\n",
       "  <line x1=\"0\" y1=\"55\" x2=\"25\" y2=\"55\" />\n",
       "  <line x1=\"0\" y1=\"58\" x2=\"25\" y2=\"58\" />\n",
       "  <line x1=\"0\" y1=\"61\" x2=\"25\" y2=\"61\" />\n",
       "  <line x1=\"0\" y1=\"64\" x2=\"25\" y2=\"64\" />\n",
       "  <line x1=\"0\" y1=\"66\" x2=\"25\" y2=\"66\" />\n",
       "  <line x1=\"0\" y1=\"69\" x2=\"25\" y2=\"69\" />\n",
       "  <line x1=\"0\" y1=\"72\" x2=\"25\" y2=\"72\" />\n",
       "  <line x1=\"0\" y1=\"75\" x2=\"25\" y2=\"75\" />\n",
       "  <line x1=\"0\" y1=\"78\" x2=\"25\" y2=\"78\" />\n",
       "  <line x1=\"0\" y1=\"81\" x2=\"25\" y2=\"81\" />\n",
       "  <line x1=\"0\" y1=\"84\" x2=\"25\" y2=\"84\" />\n",
       "  <line x1=\"0\" y1=\"87\" x2=\"25\" y2=\"87\" />\n",
       "  <line x1=\"0\" y1=\"90\" x2=\"25\" y2=\"90\" />\n",
       "  <line x1=\"0\" y1=\"93\" x2=\"25\" y2=\"93\" />\n",
       "  <line x1=\"0\" y1=\"96\" x2=\"25\" y2=\"96\" />\n",
       "  <line x1=\"0\" y1=\"99\" x2=\"25\" y2=\"99\" />\n",
       "  <line x1=\"0\" y1=\"101\" x2=\"25\" y2=\"101\" />\n",
       "  <line x1=\"0\" y1=\"104\" x2=\"25\" y2=\"104\" />\n",
       "  <line x1=\"0\" y1=\"107\" x2=\"25\" y2=\"107\" />\n",
       "  <line x1=\"0\" y1=\"110\" x2=\"25\" y2=\"110\" />\n",
       "  <line x1=\"0\" y1=\"113\" x2=\"25\" y2=\"113\" />\n",
       "  <line x1=\"0\" y1=\"116\" x2=\"25\" y2=\"116\" />\n",
       "  <line x1=\"0\" y1=\"119\" x2=\"25\" y2=\"119\" />\n",
       "  <line x1=\"0\" y1=\"120\" x2=\"25\" y2=\"120\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Vertical lines -->\n",
       "  <line x1=\"0\" y1=\"0\" x2=\"0\" y2=\"120\" style=\"stroke-width:2\" />\n",
       "  <line x1=\"25\" y1=\"0\" x2=\"25\" y2=\"120\" style=\"stroke-width:2\" />\n",
       "\n",
       "  <!-- Colored Rectangle -->\n",
       "  <polygon points=\"0.000000,0.000000 25.412617,0.000000 25.412617,120.000000 0.000000,120.000000\" style=\"fill:#ECB172A0;stroke-width:0\"/>\n",
       "\n",
       "  <!-- Text -->\n",
       "  <text x=\"12.706308\" y=\"140.000000\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" >2</text>\n",
       "  <text x=\"45.412617\" y=\"60.000000\" font-size=\"1.0rem\" font-weight=\"100\" text-anchor=\"middle\" transform=\"rotate(-90,45.412617,60.000000)\">412</text>\n",
       "</svg>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "dask.array<array, shape=(412, 2), dtype=float64, chunksize=(10, 2), chunktype=numpy.ndarray>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Interpolate the shifts for all values not in coords\n",
    "shifts = np.stack(interp_shifts(coords, [dx, dy], n=original.shape[0]), axis=1)\n",
    "neededMargins = np.ceil(shifts.max(axis=0)).astype(int)\n",
    "shifts = da.from_array(shifts, chunks=(dE,-1))\n",
    "shifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 9, the actual shifting of the original images\n",
    "\n",
    "#Inferring output dtype is not supported in dask yet, so we need original.dtype here.\n",
    "@da.as_gufunc(signature=\"(i,j),(2)->(i,j)\", output_dtypes=original.dtype, vectorize=True)\n",
    "def shift_images(image, shift):\n",
    "    \"\"\"Shift `image` by `shift` pixels.\"\"\"\n",
    "    return ndi.shift(image, shift=shift, order=1)\n",
    "\n",
    "padded = da.pad(original.rechunk({0:dE}), \n",
    "                ((0, 0), \n",
    "                 (0, neededMargins[0]), \n",
    "                 (0, neededMargins[1])\n",
    "                ),\n",
    "                mode='constant'\n",
    "               )\n",
    "corrected = shift_images(padded.rechunk({1:-1, 2:-1}), shifts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "727e8efc6c2e4f4bb659e8a80bcdaa5a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=103, continuous_update=False, description='n', max=411), Output()), _dom…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Do an interactive viewer to inspect the results\n",
    "interactive(lambda n: plot_stack(corrected, n, grid=True), \n",
    "            n=widgets.IntSlider(corrected.shape[0]//4,0,corrected.shape[0]-1,1, continuous_update=False)\n",
    "           ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving data\n",
    "Save the resulting data in a new netCDF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "xrcorrected = dataset.reindex({'x': np.arange(0, dataset.x[1]*corrected.shape[1], dataset.x[1]), \n",
    "                               'y': np.arange(0, dataset.y[1]*corrected.shape[2], dataset.y[1])})\n",
    "xrcorrected.Intensity.data = corrected\n",
    "xrcorrected.Intensity.attrs['DriftCorrected'] = 'True'\n",
    "xrcorrected.to_netcdf(os.path.join(folder, name + '_driftcorrected.nc'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Or, save the results to zarr\n",
    "import zarr\n",
    "from numcodecs import Blosc\n",
    "compressor = Blosc(cname='zstd', clevel=1, shuffle=Blosc.SHUFFLE)\n",
    "corrected.to_zarr(os.path.join(folder, name + '_driftcorrected.zarr'), \n",
    "                  overwrite=True, compressor=compressor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Or, although parallel access to HDF5 is hard, so go single process, save to hdf5\n",
    "with dask.config.set(scheduler='threads'):  \n",
    "    corrected.to_hdf5(os.path.join(folder, name + '_driftCorrected.h5', '/Intensity')"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:light"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}